You said a few posts up that
if you want to use a recent sphinx you have to use a recent matplotlib too:
sphinx-doc/sphinx#7659
So the chain goes like this: compile with default Buster libs --> Bibliography error with [dimitriadis2016] --> install numpydoc>=0.9 with pip --> pip pulls sphinx 3.1.2 --> default maplotlib raises an Exception --> upgrade matplotlib to 3.3.0--> maplotib 3.3.0 raises an other (non-critical) Exception!
WeibullMinFactory
class is order to fit the distribution parameters to a given sample. I thought it was Cobyla, but that does not seem to be the case: I tried putting additional warnings in my copy of the code in Cobyla::ComputeObjectiveAndConstraint
but they are never reached.
f = ot.SymbolicFunction(['x'], ['x * sin(x)', 'x * cos(x)'])
(this is my only modification). There is no problem with python2.7 but I get a Segmentation Fault with python3.8. Can someone help me ? Have some idea about the origin of this issue :-)
Hi,
I try to compute the Sobol indices, so I consider the documentation example:
import openturns as ot
ot.RandomGenerator.SetSeed(0)
formula = ['sin(pi_*X1)+7*sin(pi_*X2)^2+0.1*(pi_*X3)^4*sin(pi_*X1)']
model = ot.SymbolicFunction(['X1', 'X2', 'X3'], formula)
distribution = ot.ComposedDistribution([ot.Uniform(-1.0, 1.0)] * 3)
size = 100
inputDesign = ot.SobolIndicesExperiment(distribution, size).generate()
outputDesign = model(inputDesign)
sensitivityAnalysis = ot.SaltelliSensitivityAlgorithm(inputDesign, outputDesign, size)
print(sensitivityAnalysis.getFirstOrderIndices())
and I've tried the four algorithms: Saltelli, Martinez, MauntzKucherenko and Jansen.
I'm wondering why the indices results are all different? And when I change the size argument in the SensitivityAlgorithm method, the results are also very differents.
What's happen?
Hi @christelle-git, these algorithms provide statistical estimators of the Sobol indices, not the true values of the indices. They become more accurate as the sample size increases. Here is a script (based on yours) to illustrate this.
import openturns as ot
import numpy as np
ot.RandomGenerator.SetSeed(1)
formula = ['sin(pi_*X1)+7*sin(pi_*X2)^2+0.1*(pi_*X3)^4*sin(pi_*X1)']
model = ot.SymbolicFunction(['X1', 'X2', 'X3'], formula)
distribution = ot.ComposedDistribution([ot.Uniform(-1.0, 1.0)] * 3)
def sobol_estimators(size):
"""
Estimate Sobol first order indices with Saltelli, Martinez, Mauntz-Kucherenko, Jansen.
Return the standard deviation of these estimates for each input random variable.
"""
inputDesign = ot.SobolIndicesExperiment(distribution, size).generate()
outputDesign = model(inputDesign)
sensitivityAnalysis = ot.SaltelliSensitivityAlgorithm(inputDesign, outputDesign, size)
saltelli = np.array(sensitivityAnalysis.getFirstOrderIndices())
print("Saltelli : {}".format(saltelli))
sensitivityAnalysis = ot.MartinezSensitivityAlgorithm(inputDesign, outputDesign, size)
martinez = np.array(sensitivityAnalysis.getFirstOrderIndices())
print("Martinez : {}".format(martinez))
sensitivityAnalysis = ot.MauntzKucherenkoSensitivityAlgorithm(inputDesign, outputDesign, size)
mauntz_kucherenko = np.array(sensitivityAnalysis.getFirstOrderIndices())
print("Mauntz-Kucherenko: {}".format(mauntz_kucherenko))
sensitivityAnalysis = ot.JansenSensitivityAlgorithm(inputDesign, outputDesign, size)
jansen = np.array(sensitivityAnalysis.getFirstOrderIndices())
print("Jansen: {}".format(jansen))
return np.vstack([saltelli, martinez, mauntz_kucherenko, jansen]).std(axis=0)
for size in [100, 1000, 10000]:
print("Size =", size)
print("Estimators std deviation:{}".format(sobol_estimators(size)))
print("############################################################")
The output is:
Size = 100
Saltelli : [ 0.18075179 0.42489663 -0.31641387]
Martinez : [ 0.15745024 0.39894551 -0.23866791]
Mauntz-Kucherenko: [0.51365383 0.75779868 0.01648818]
Jansen: [ 0.01032805 0.3205321 -0.6364488 ]
Estimators std deviation:[0.1840981 0.16742825 0.23304392]
############################################################
Size = 1000
Saltelli : [0.31550553 0.47364577 0.03620589]
Martinez : [0.30481089 0.47386356 0.03576605]
Mauntz-Kucherenko: [0.30604697 0.46418721 0.02674733]
Jansen: [0.27867508 0.47437777 0.02030996]
Estimators std deviation:[0.01367971 0.0042411 0.00663328]
############################################################
Size = 10000
Saltelli : [ 0.32205724 0.45480035 -0.01725293]
Martinez : [ 0.31796198 0.44660686 -0.01709273]
Mauntz-Kucherenko: [ 0.32709197 0.45983508 -0.0122182 ]
Jansen: [ 0.3097761 0.43638247 -0.03198465]
Estimators std deviation:[0.0063401 0.00887852 0.00741046]
############################################################
For 3 different sample sizes, the script computes the 4 different estimates of the Sobol first order indices and the empirical standard deviation of these 4 estimates. You can see that, as the sample size increases, the standard deviation vanishes, because all 4 estimators converge to the true values of the indices.
draw()
method which all SensitivityAlgorithm
classes possess to draw confidence intervals. Here is the same example with the SaltelliSensitivityAlgorithm
:import openturns as ot
from openturns.viewer import View
ot.RandomGenerator.SetSeed(0)
formula = ['sin(pi_*X1)+7*sin(pi_*X2)^2+0.1*(pi_*X3)^4*sin(pi_*X1)']
model = ot.SymbolicFunction(['X1', 'X2', 'X3'], formula)
distribution = ot.ComposedDistribution([ot.Uniform(-1.0, 1.0)] * 3)
for size in [100, 1000]:
inputDesign = ot.SobolIndicesExperiment(distribution, size).generate()
outputDesign = model(inputDesign)
sensitivityAnalysis = ot.SaltelliSensitivityAlgorithm(inputDesign, outputDesign, size)
graph = sensitivityAnalysis.draw()
graph.setTitle(graph.getTitle() + " \n Size " + str(size)
+ " - confidence intervals with level " + str(sensitivityAnalysis.getConfidenceLevel()))
View(graph)
FunctionalChaosAlgorithm
serves to compute the surrogate model and the class FunctionalChaosSobolIndices
then takes it and uses it to estimate Sobol indices. The example Sobol’ sensitivity indices from chaos shows how to do this.
Thank you, several notes:
1) when I execute your code from yesterday 22:00, the output is the same but full of warning, how can I avoid it ?
2) when I execute your code today 9:26, I have no plotted graph...
3) I've already tried the Chaos example, but also I have no graph too, instead I have:
Out[9]: class=Graph name=Sobol' indices implementation=class=GraphImplementation name=Sobol' indices title=Sobol' indices xTitle=inputs yTitle=index value axes=ON grid=ON legendposition=topright legendFontSize=1 drawables=[class=Drawable name=First order implementation=class=Cloud name=First order derived from class=DrawableImplementation name=First order legend=First order data=class=Sample name=Unnamed implementation=class=SampleImplementation name=Unnamed size=5 dimension=2 data=[[1,0],[2,0],[3,0],[4,0.000403291],[5,0]] color=red fillStyle=solid lineStyle=solid pointStyle=circle lineWidth=1,class=Drawable name=Total order implementation=class=Cloud name=Total order derived from class=DrawableImplementation name=Total order legend=Total order data=class=Sample name=Unnamed implementation=class=SampleImplementation name=Unnamed size=5 dimension=2 data=[[1.125,0],[2.125,0],[3.125,0],[4.125,0.000106662],[5.125,0.999893]] color=blue fillStyle=solid lineStyle=solid pointStyle=square lineWidth=1,class=Drawable name=Unnamed implementation=class=Text name=Unnamed derived from class=DrawableImplementation name=Unnamed legend= data=class=Sample name=Unnamed implementation=class=SampleImplementation name=Unnamed size=5 dimension=2 data=[[1.25,0],[2.25,0],[3.25,0],[4.25,0.000254976],[5.25,0.499947]] color=black fillStyle=solid lineStyle=solid pointStyle=plus lineWidth=1]
4) my input and output are csv files with discrete points (values of parameters) from geology simulations and the output is the corresponding numerical response.
1) I typically run code within Spyder. Spyder suppresses warnings. There might be a way to disable warnings directly from within OpenTURNS, but I don't know about it. Maybe @jschueller has an idea?
2) Did you run the script from an IDE like Spyder or from the command line? When I run the script from Spyder, the pictures are displayed. From the command line, they are not. However, they can be saved to disk. Here is a modified version that does this:
import openturns as ot
from openturns.viewer import View
ot.RandomGenerator.SetSeed(0)
formula = ['sin(pi_*X1)+7*sin(pi_*X2)^2+0.1*(pi_*X3)^4*sin(pi_*X1)']
model = ot.SymbolicFunction(['X1', 'X2', 'X3'], formula)
distribution = ot.ComposedDistribution([ot.Uniform(-1.0, 1.0)] * 3)
for size in [100, 1000]:
inputDesign = ot.SobolIndicesExperiment(distribution, size).generate()
outputDesign = model(inputDesign)
sensitivityAnalysis = ot.SaltelliSensitivityAlgorithm(inputDesign, outputDesign, size)
graph = sensitivityAnalysis.draw()
graph.setTitle(graph.getTitle() + " \n Size " + str(size)
+ " - confidence intervals with level " + str(sensitivityAnalysis.getConfidenceLevel()))
view = View(graph)
view.save('Sobol_size_{}.pdf'.format(size))
draw()
method inside a Jupyter notebook: there is no need to call View
.
I really think you would be better off using FunctionalChaosSobolIndices
, which works regardless of sample size (though of course, the bigger the sample, the better).
Assuming your input sample is an OpenTURNS Sample
called X and your output an OpenTURNS Sample
called Y:
algo = ot.FunctionalChaosAlgorithm(X, Y)
algo.run()
result = algo.getResult()
sensitivityAnalysis = ot.FunctionalChaosSobolIndices(result)
# draw Sobol' indices
first_order = [sensitivityAnalysis.getSobolIndex(i) for i in range(X.getDimension())]
total_order = [sensitivityAnalysis.getSobolTotalIndex(i) for i in range(X.getDimension())]
graph = ot.SobolIndicesAlgorithm.DrawSobolIndices(X.getDescription(), first_order, total_order)
# save graph
from openturns.viewer import View
view = View(graph)
view.save('Sobol_indices.pdf')
Hi All,
I just compiled OpenTURNS (Debian 10 'buster') and I noticed this warning:
/home/devel/GIT/openturns/build-hohenbichler/python/src/CMakeFiles/common.dir/common_modulePYTHON_wrap.cxx: In function ‘void SWIG_Python_addvarlink(PyObject*, char*, PyObject* (*)(), int (*)(PyObject*))’:
/home/devel/GIT/openturns/build-hohenbichler/python/src/CMakeFiles/common.dir/common_modulePYTHON_wrap.cxx:35736:16: warning: ‘char* strncpy(char*, const char*, size_t)’ specified bound depends on the length of the source argument [-Wstringop-overflow=]
strncpy(gv->name,name,size);
~~~~~~~^~~~~~~~~~~~~~~~~~~~
/home/devel/GIT/openturns/build-hohenbichler/python/src/CMakeFiles/common.dir/common_modulePYTHON_wrap.cxx:35733:27: note: length computed here
size_t size = strlen(name)+1;
However, it's not a big deal since it does not stop the compiling process.
Thanks for your help
Hi! I ran into this error during compilation of the library:
[ 25%] Building CXX object lib/src/CMakeFiles/OT.dir/Base/Stat/Last.cxx.o
/home/devel/GIT/openturns/lib/src/Base/Stat/HMatrixFactory.cxx: In member function ‘OT::HMatrix OT::HMatrixFactory::build(const OT::Sample&, OT::UnsignedInteger, OT::Bool, const OT::HMatrixParameters&)’:
/home/devel/GIT/openturns/lib/src/Base/Stat/HMatrixFactory.cxx:83:12: error: ‘struct hmat_settings_t’ has no member named ‘assemblyEpsilon’; did you mean ‘acaEpsilon’?
settings.assemblyEpsilon = parameters.getAssemblyEpsilon();
^~~~~~~~~~~~~~~
acaEpsilon
/home/devel/GIT/openturns/lib/src/Base/Stat/HMatrixFactory.cxx:84:12: error: ‘struct hmat_settings_t’ has no member named ‘recompressionEpsilon’; did you mean ‘compressionMethod’?
settings.recompressionEpsilon = parameters.getRecompressionEpsilon();
^~~~~~~~~~~~~~~~~~~~
compressionMethod
make[2]: *** [lib/src/CMakeFiles/OT.dir/build.make:4158: lib/src/CMakeFiles/OT.dir/Base/Stat/HMatrixFactory.cxx.o] Error 1
make[2]: *** Attente des tâches non terminées....
/home/devel/GIT/openturns/lib/src/Base/Stat/HMatrixImplementation.cxx: In member function ‘OT::Bool OT::HMatrixImplementation::setKey(const String&, const String&)’:
/home/devel/GIT/openturns/lib/src/Base/Stat/HMatrixImplementation.cxx:594:21: error: ‘struct hmat_settings_t’ has no member named ‘assemblyEpsilon’; did you mean ‘acaEpsilon’?
iss >> settings.assemblyEpsilon;
^~~~~~~~~~~~~~~
acaEpsilon
/home/devel/GIT/openturns/lib/src/Base/Stat/HMatrixImplementation.cxx:599:21: error: ‘struct hmat_settings_t’ has no member named ‘recompressionEpsilon’; did you mean ‘compressionMethod’?
iss >> settings.recompressionEpsilon;
^~~~~~~~~~~~~~~~~~~~
compressionMethod
make[2]: *** [lib/src/CMakeFiles/OT.dir/build.make:4184: lib/src/CMakeFiles/OT.dir/Base/Stat/HMatrixImplementation.cxx.o] Error 1
make[1]: *** [CMakeFiles/Makefile2:823: lib/src/CMakeFiles/OT.dir/all] Error 2
make: *** [Makefile:163: all] Error 2
This error happened after I compiled hmat-oss from the Github sources (master branch). Do you know what causes this error?
Hi,
I am working on the estimation of a failure probability using Monte Carlo. I have quite a large amount of input variables and some are correlated. Therefore, I would like to create the joint distribution using a copula. The copula function takes the correlation matrix as input. As far as I know, the correlation matrix has to be created using ot.Correlation Matrix and you have to manually add the correlations like this: R[0,1] = 0.8 for example, which is a lot of work for a large amount of input variables and it makes me lose the overview .
Is there a way to define the correlation matrix in e.g. a csv file instead and use this in ot.NormalCopula? If I try to load a csv file into an array, ot.NormalCopula does not recognize my correlation matrix.
Thanks!
Hi,
Another question. I am trying to create the event whose probability I want to estimate, using ot.CompositeRandomVector(model,vector). I saw some examples using ot.SymbolicFunction for the model. However, my function is quite elaborate and also contains some variables without distributions, so I am looking for an easier way to define the model. I am trying to achieve this by defining the model as a regular python function:
def Z_GEKB(Hm0,h,Tmm10,crest_height,alpha):
Rc = crest_height - h #freeboard
br = math.tan(alpha) / math.sqrt( (Hm0*2*np.pi) / (g*(Tmm10)**2) )
q = math.sqrt(g * (Hm0)**3) * ( (0.023/math.sqrt(math.tan(alpha))) * br * math.exp(-(2.7 * (Rc / (br*Hm0)))**(1.3)) )
return math.log10(m_c * q_c) - math.log10(m_a * q)
vect = ot.RandomVector(jointdist)
out = ot.CompositeRandomVector(Z_GEKB,vect)
But then I get the error as sent above.
OpenTURNS now has a forum!
This forum will allow us to discuss OpenTURNS in a more structured way, with threads dedicated to specific use cases, methods, or issues. More importantly, it will allow us to preserve discussions that would be lost on Gitter.
As a bonus feature, the forum supports LaTeX. This should make mathematical discussions easier.