Skip to content

Commit

Permalink
Migrate remaining notebooks
Browse files Browse the repository at this point in the history
  • Loading branch information
jschueller committed Sep 20, 2024
1 parent 2c25938 commit 8da05ef
Show file tree
Hide file tree
Showing 36 changed files with 797 additions and 9,643 deletions.
2 changes: 1 addition & 1 deletion doc/conf.py
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@
'examples_dirs': ['examples'],
'gallery_dirs': ['auto_examples'],
'show_signature': False,
'ignore_pattern': 'rp8.py|rp33'
'ignore_pattern': 'rp8.py|rp33|drawevent|draw_cross_cuts|plot_bbrc'
}

if Version(sphinx_gallery.__version__) >= Version("0.17.0"):
Expand Down
117 changes: 117 additions & 0 deletions doc/examples/plot_check_reliability_reference_probabilities.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,117 @@
"""
Check reliability problems reference probabilities
==================================================
"""

# %%
# In this example, we check that the reference probabilities in the reliability problems
# are consistent with confidence bounds from Monte-Carlo simulations.
# These 95% confidence bounds are stored in 'reliability_compute_reference_proba.csv'
# and required approximately than :math`10^9` function evaluations for each problem.
#
# We consider two different metrics:
#
# * we check if the reference probability is within the 95% confidence bounds,
# * we compute the number of significant digits by comparing the Monte-Carlo estimator and the reference value.
#
# The number of significant digits may be as high as 17 when all decimal digits are correct.
# However, the reference probabilities are only known up to 3 digits for most problems.
# In order to keep some safeguard, we will be happy with 2 correct digits.
#
# These metrics may fail.
#
# * On average, a fraction equal to 5% of the estimates are outside the confidence bounds.
# * The Monte-Carlo estimator may not be accurate enough, e.g. if the probability is very close to zero.

# %%
import otbenchmark as otb
import pandas as pd


# %%
df = pd.read_csv("reliability_compute_reference_proba.csv")
df


# %%
data = df.values


# %%
pf_reference = data[:, 1]
pmin = data[:, 3]
pmax = data[:, 4]


# %%
benchmarkProblemList = otb.ReliabilityBenchmarkProblemList()
numberOfProblems = len(benchmarkProblemList)
numberOfProblems


# %%
digitsMinimum = 2


# %%
categoryA = []
categoryB = []
categoryC = []
categoryD = []


# %%
for i in range(numberOfProblems):
problem = benchmarkProblemList[i]
name = problem.getName()
pf = problem.getProbability()
event = problem.getEvent()
antecedent = event.getAntecedent()
distribution = antecedent.getDistribution()
dimension = distribution.getDimension()
if pf > pmin[i] and pf < pmax[i]:
tagBounds = "In"
else:
tagBounds = "Out"
digits = otb.ComputeLogRelativeError(pf_reference[i], pf)
if tagBounds == "In" and digits >= digitsMinimum:
categoryA.append(name)
elif tagBounds == "Out" and digits >= digitsMinimum:
categoryB.append(name)
elif tagBounds == "In" and digits < digitsMinimum:
categoryC.append(name)
else:
categoryD.append(name)
print(
"#%d, %-10s, pf=%.2e, ref=%.2e, C.I.=[%.2e,%.2e], digits=%d : %s"
% (i, name[0:10], pf, pf_reference[i], pmin[i], pmax[i], digits, tagBounds)
)

# %%
# There are four different cases.
#
# * Category A: all good. For some problems, both metrics are correct in the sense
# that the reference probability is within the bounds and the number of significant digits is larger than 2.
# The RP24, RP55, RP110, RP63, R-S, Axial stressed beam problems fall in that category.
# * Category B: correct digits, not in bounds.
# We see that the RP8 problem has a reference probability outside of the 95% confidence bounds,
# but has 2 significant digits.
# * Category C: insufficient digits, in bounds. The difficult RP28 problem fall in that category.
# * Category D: insufficient digits, not in bounds. These are suspicious problems.

# %%
print(categoryA)

# %%
print(categoryB)

# %%
print(categoryC)

# %%
print(categoryD)

# %%
# The number of suspicious problems seems very large.
# However, we notice that all these cases are so that the reference probability is close,
# in absolute value, to the Monte-Carlo estimator.
File renamed without changes.
108 changes: 108 additions & 0 deletions doc/examples/plot_draw_cross_cuts.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,108 @@
"""
Draw multidimensional functions, distributions and events
=========================================================
"""

# %%
# This example shows how to represent multidimensional functions, distributions and events.
# When 2D plots are to draw, contours are used.
# We use 2D cross-sections to represent multidimensional objects when required,
# which leads to cross-cuts representations.

# %%
import otbenchmark as otb
import openturns.viewer as otv
import matplotlib.pyplot as plt

# %%
problem = otb.ReliabilityProblem33()

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


# ## Compute the bounds of the domain

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

# %%
inputDimension = distribution.getDimension()
inputDimension

# %%
alpha = 1 - 1.0e-5

bounds, marginalProb = distribution.computeMinimumVolumeIntervalWithMarginalProbability(
alpha
)

# %%
referencePoint = distribution.getMean()
referencePoint

# %%
crossCut = otb.CrossCutFunction(g, referencePoint)
fig = crossCut.draw(bounds)
# Remove the legend labels because there
# are too many for such a small figure
for ax in fig.axes:
ax.legend("")
# Increase space between sub-figures so that
# there are no overlap
plt.subplots_adjust(hspace=0.4, wspace=0.4)

# %%
# Plot cross-cuts of the distribution
crossCut = otb.CrossCutDistribution(distribution)

# %%
fig = crossCut.drawMarginalPDF()
# Remove the legend labels because there
# are too many for such a small figure
for ax in fig.axes:
ax.legend("")
# Increase space between sub-figures so that
# there are no overlap
plt.subplots_adjust(hspace=0.4, wspace=0.4)

# %%
# The correct way to represent cross-cuts of a distribution is to draw the contours
# of the PDF of the conditional distribution.
fig = crossCut.drawConditionalPDF(referencePoint)
# Remove the legend labels because there
# are too many for such a small figure
for ax in fig.axes:
ax.legend("")
# Increase space between sub-figures so that
# there are no overlap
plt.subplots_adjust(hspace=0.4, wspace=0.4)

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

# %%
drawEvent = otb.DrawEvent(event)

# %%
_ = drawEvent.drawLimitState(bounds)

# %%
# In the following figure, we present the cross-cuts of samples with size equal to 500.
# These are three different samples, each of which is plotted with the `drawSampleCrossCut` method.
# For each cross-cut plot (i,j), the current implementation uses the marginal bivariate distribution,
# then generates a sample from this distribution.
# A more rigorous method would draw the conditional distribution, but this might reduce the performance in general.
# See https://github.com/mbaudin47/otbenchmark/issues/47 for details.
sampleSize = 500
_ = drawEvent.drawSample(sampleSize)

# %%
_ = drawEvent.fillEvent(bounds)

# %%
otv.View.ShowAll()
117 changes: 117 additions & 0 deletions doc/examples/plot_drawevent.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,117 @@
"""
Demonstration of the DrawEvent class
====================================
"""

import otbenchmark as otb
import openturns.viewer as otv

# %%
# 3D problem
# ----------

# %%
problem = otb.ReliabilityProblem33()

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

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

# %%
inputDimension = distribution.getDimension()
inputDimension

# %%
alpha = 1 - 0.00001

# %%
bounds, marginalProb = distribution.computeMinimumVolumeIntervalWithMarginalProbability(
alpha
)

# %%
referencePoint = distribution.getMean()
referencePoint

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

# %%
drawEvent = otb.DrawEvent(event)

# %%
# The highest level method is the `draw` method which flags allow to gather various graphics into a single one.
_ = drawEvent.draw(bounds)

# %%
_ = drawEvent.draw(bounds, fillEvent=True)

# %%
# The `drawLimitState` method only draws the limit state.
_ = drawEvent.drawLimitState(bounds)

# %%
# The `drawSample` method plots a sample with a color code which specifies which points are inside or outside the event.
sampleSize = 500
_ = drawEvent.drawSample(sampleSize)

# %%
_ = drawEvent.fillEvent(bounds)

# %%
# 2D problem
# ----------
#
# When the problem has 2 dimensions, single cross-cuts are sufficient.
# This is why we use the `*CrossCut` methods.

# %%
problem = otb.ReliabilityProblem22()

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

# %%
inputVector = event.getAntecedent()
distribution = inputVector.getDistribution()
bounds, marginalProb = distribution.computeMinimumVolumeIntervalWithMarginalProbability(
1.0 - 1.0e-6
)

# %%
sampleSize = 10000
drawEvent = otb.DrawEvent(event)

# %%
cloud = drawEvent.drawSampleCrossCut(sampleSize)
_ = otv.View(cloud)

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

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

# %%
domain.add(cloud)
_ = otv.View(domain)

# %%
# For a 2D sample, it is sometimes handy to re-use a precomputed sample.
# In this case, we can use the `drawInputOutputSample` method.
inputSample = distribution.getSample(sampleSize)
outputSample = g(inputSample)
drawEvent.drawInputOutputSample(inputSample, outputSample)

# %%
otv.View.ShowAll()
Loading

0 comments on commit 8da05ef

Please sign in to comment.