Skip to content

Commit

Permalink
Doc: Add reliability problems
Browse files Browse the repository at this point in the history
  • Loading branch information
jschueller committed Sep 20, 2024
1 parent 3bb8c8e commit 2c25938
Show file tree
Hide file tree
Showing 71 changed files with 3,005 additions and 17,785 deletions.
1 change: 1 addition & 0 deletions doc/conf.py
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,7 @@
'examples_dirs': ['examples'],
'gallery_dirs': ['auto_examples'],
'show_signature': False,
'ignore_pattern': 'rp8.py|rp33'
}

if Version(sphinx_gallery.__version__) >= Version("0.17.0"):
Expand Down
2 changes: 2 additions & 0 deletions doc/examples/reliability_problems/README.txt
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
Reliability examples
====================
122 changes: 122 additions & 0 deletions doc/examples/reliability_problems/plot_case_rs.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,122 @@
"""
R-S analysis and 2D graphics
============================
"""

# %%
# The objective of this example is to present the R-S problem.
# We also present graphic elements for the visualization of the limit state surface in 2 dimensions.

# %%
import openturns as ot
import openturns.viewer as otv
import otbenchmark as otb

# %%
problem = otb.RminusSReliability()


# %%
event = problem.getEvent()
g = event.getFunction()


# %%
problem.getProbability()

# %%
# Create the Monte-Carlo algorithm
algoProb = ot.ProbabilitySimulationAlgorithm(event)
algoProb.setMaximumOuterSampling(1000)
algoProb.setMaximumCoefficientOfVariation(0.01)
algoProb.run()


# %%
# Get the results
resultAlgo = algoProb.getResult()
neval = g.getEvaluationCallsNumber()
print("Number of function calls = %d" % (neval))
pf = resultAlgo.getProbabilityEstimate()
print("Failure Probability = %.4f" % (pf))
level = 0.95
c95 = resultAlgo.getConfidenceLength(level)
pmin = pf - 0.5 * c95
pmax = pf + 0.5 * c95
print("%.1f %% confidence interval :[%.4f,%.4f] " % (level * 100, pmin, pmax))

# %%
# ## Plot the contours of the function

# %%
inputVector = event.getAntecedent()
distribution = inputVector.getDistribution()

# %%
R = distribution.getMarginal(0)
S = distribution.getMarginal(1)

# %%
alphaMin = 0.001
alphaMax = 1 - alphaMin
lowerBound = ot.Point([R.computeQuantile(alphaMin)[0], S.computeQuantile(alphaMin)[0]])
upperBound = ot.Point([R.computeQuantile(alphaMax)[0], S.computeQuantile(alphaMax)[0]])

# %%
nbPoints = [100, 100]
_ = otv.View(g.draw(lowerBound, upperBound, nbPoints))

# %%
Y = R - S
Y

# %%
_ = otv.View(Y.drawPDF())

# %%
# Print the iso-values of the distribution
# ----------------------------------------

# %%
_ = otv.View(distribution.drawPDF())

# %%
# ## Visualise the safe and unsafe regions on a sample

# %%
sampleSize = 500

# %%
drawEvent = otb.DrawEvent(event)

# %%
cloud = drawEvent.drawSampleCrossCut(sampleSize)
cloud

# %%
# ## Draw the limit state surface

# %%
bounds = ot.Interval(lowerBound, upperBound)
bounds


# %%
graph = drawEvent.drawLimitStateCrossCut(bounds)
graph.add(cloud)
graph

# %%
# ## Fill the event domain with a color

# %%
domain = drawEvent.fillEventCrossCut(bounds)
_ = otv.View(domain)

# %%
domain.setLegends(["", ""])
domain.add(cloud)
_ = otv.View(domain)

# %%
otv.View.ShowAll()
Loading

0 comments on commit 2c25938

Please sign in to comment.