Exploration of the effect of aleatoric uncertainties on Sobol Coefficients
In this notebook we will explore the role od aleatoric uncertainties on the calculation od statistical quantities such as Sobol Coefficients.
To do this we take the results of an existing run using the fusion tutorial model and then perturb the “measured” values at the various positions determined by the PCE campaign, and recalculate the Sobol indices. We do this for a number of cases and explore the impact of the chosen noise level on the calculated Sobol indices.
[1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import pickle
import easyvvuq as uq
import os
import ast
import time
import pandas as pd
/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/__init__.py:9: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
import pkg_resources
[2]:
# This tutorial builds on the fusion example
import fusion
Te, ne, rho, rho_norm = fusion.solve_Te()
[3]:
# initialize the random number generator
rng = np.random.default_rng()
[4]:
# function that will be used to add correlated noise
def randomize(x, ng=11, lcorr=0.2):
y = 0.0
for x0, r in zip(np.linspace(0,1,ng), rng.standard_normal(ng)):
y += np.exp(-((x-x0)/lcorr)**2)*r
return y
[5]:
# Demonstrate how the randomize function works
Randomize = np.array([randomize(rho_norm) for i in range(1000)])
plt.clf()
plt.plot(rho_norm, Randomize.mean(axis=0), lw=5)
plt.fill_between(rho_norm, Randomize.mean(axis=0)-Randomize.std(axis=0), Randomize.mean(axis=0)+Randomize.std(axis=0), alpha=0.5)
plt.plot(rho_norm, Randomize.T, alpha=0.05)
plt.xlabel('rho_tor_norm')
plt.ylabel('Random number')
plt.title('Based on 1000 samples')
[5]:
Text(0.5, 1.0, 'Based on 1000 samples')
[6]:
# Show what 10% noise applied to the Te profile looks like
Te_scan = np.array([Te+Te/10*randomize(rho_norm) for i in range(1000)])
Te_scan = Te*(1+0.1*Randomize)
plt.clf()
plt.plot(rho_norm, Te_scan.T, alpha=0.01)
plt.plot(rho_norm, Te_scan.mean(axis=0), lw=5)
plt.fill_between(rho_norm, Te_scan.mean(axis=0)-Te_scan.std(axis=0), Te_scan.mean(axis=0)+Te_scan.std(axis=0), alpha=0.5)
plt.plot(rho_norm, Te, '--')
[6]:
[<matplotlib.lines.Line2D at 0x15515b2c0>]
[7]:
# Now run an EasyVVUQ campaign
%run ./easyvvuq_fusion_dask_tutorial.py -l
Running locally
/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/__init__.py:9: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
import pkg_resources
/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/__init__.py:9: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
import pkg_resources
/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/__init__.py:9: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
import pkg_resources
/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/__init__.py:9: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
import pkg_resources
/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/__init__.py:9: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
import pkg_resources
/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/__init__.py:9: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
import pkg_resources
/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/__init__.py:9: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
import pkg_resources
/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/__init__.py:9: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
import pkg_resources
/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/__init__.py:9: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
import pkg_resources
/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/__init__.py:9: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
import pkg_resources
/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/__init__.py:9: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
import pkg_resources
/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/__init__.py:9: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
import pkg_resources
/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/__init__.py:9: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
import pkg_resources
/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/__init__.py:9: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
import pkg_resources
<Client: 'tcp://127.0.0.1:57155' processes=14 threads=14, memory=64.00 GiB>
Time for phase 1 = 3.422
Number of samples = 243
Time for phase 2 = 0.170
/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/cerberus/validator.py:618: UserWarning: These types are defined both with a method and in the'types_mapping' property of this validator: {'integer'}
warn(
/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/cerberus/validator.py:618: UserWarning: These types are defined both with a method and in the'types_mapping' property of this validator: {'integer'}
warn(
Time for phase 3 = 9.439
Running locally
Running locally
Running locally
Running locally
Running locally
Running locally
Running locally
Running locally
Running locally
Running locally
Running locally
Running locally
Running locally
Running locally
C
Time for phase 4 = 0.550
Time for phase 5 = 0.030
/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/distributions/kernel/baseclass.py:84: RuntimeWarning: divide by zero encountered in matmul
self._pcovariance = numpy.matmul(numpy.matmul(
/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/distributions/kernel/baseclass.py:84: RuntimeWarning: overflow encountered in matmul
self._pcovariance = numpy.matmul(numpy.matmul(
/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/distributions/kernel/baseclass.py:84: RuntimeWarning: invalid value encountered in matmul
self._pcovariance = numpy.matmul(numpy.matmul(
Time for phase 6 = 3.352
Time for phase 7 = 0.002
Time for phase 8 = 0.015
<Figure size 640x480 with 0 Axes>
[8]:
# We grab the saved results
results = pickle.load(open('easyvvuq_fusion_dask_tutorial/fusion_results.pickle', 'rb'))
[9]:
# And plot some results
rho_norm = results.describe("rho_norm", "mean")
te_mean = results.describe("te", "mean")
te_std = results.describe("te", "std")
te_10_pct = results.describe("te", "10%")
te_90_pct = results.describe("te", "90%")
te_min = results.describe("te", "min")
te_max = results.describe("te", "max")
plt.figure()
plt.plot(rho, te_mean, "b-", label="Mean")
plt.plot(rho, te_mean - te_std, "b--", label="+1 std deviation")
plt.plot(rho, te_mean + te_std, "b--")
plt.fill_between(rho, te_mean - te_std, te_mean + te_std,
color="b", alpha=0.2)
plt.plot(rho, te_10_pct, "b:", label="10 and 90 percentiles")
plt.plot(rho, te_90_pct, "b:")
plt.fill_between(rho, te_10_pct, te_90_pct, color="b", alpha=0.1)
plt.fill_between(rho, te_min, te_max, color="b", alpha=0.05)
plt.legend(loc=0)
plt.xlabel("rho [m]")
plt.ylabel("Te [eV]");
[10]:
plt.clf()
plt.contourf(rho_norm, rho_norm, results.raw_data['correlation_matrices']['te'])
plt.colorbar()
plt.xlabel('rho_tor_norm')
plt.ylabel('rho_tor_norm')
plt.title('Te correlation matrix')
plt.contour(rho_norm, rho_norm, results.raw_data['correlation_matrices']['te'], levels=[0.9, 0.99, 0.999, 0.9999, 0.99999, 0.999999], colors='k');
[11]:
# We recover the old campaign
DIR = 'easyvvuq_fusion_dask_tutorial'
old_campaign = uq.Campaign(name="fusion_pce.", db_location= f'sqlite:///{os.path.abspath(os.curdir)}/{DIR}/campaign.db')
old_runs = old_campaign.list_runs()
results_df = old_campaign.get_collation_result()
[12]:
# Define a function for perturbing an old result
def get_case(old_runs, noise=0.1):
new_runs = []
for r in old_runs:
d = r[1].copy()
d['run_id'] = r[0]
res = ast.literal_eval(d['result'])
res['te'] = list(np.array(res['te']) * (1+noise*randomize(np.array(res['rho_norm']))))
d['result'] = res
new_runs.append(d)
df_runs=[]
for d in new_runs:
D = {**d['params'], **d['result']}
D['run_name'] = d['run_name']
D['run_dir'] = d['run_dir']
D['run_id'] = d['run_id']
pd_result={}
for key in D.keys():
if not isinstance(D[key], list):
try:
pd_result[(key, 0)].append(D[key])
except KeyError:
pd_result[(key, 0)] = D[key]
else:
for i, elt in enumerate(D[key]):
try:
pd_result[(key, i)].append(D[key][i])[0]
except KeyError:
pd_result[(key, i)] = [D[key][i]][0]
df_runs.append(pd_result)
df = pd.DataFrame(df_runs)
df.columns = pd.MultiIndex.from_tuples(df.columns)
return df
[13]:
# See if we recover the old results with zero noise
df = get_case(old_runs, noise=0.0)
# Here we create an analysis instance based on the sampler from the old campaign
analysis = uq.analysis.PCEAnalysis(sampler=old_campaign.get_active_sampler(), qoi_cols=results.qois)
# And use this to perform the PCE analysis based on the just created dataframe
R = analysis.analyse(df)
/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/distributions/kernel/baseclass.py:84: RuntimeWarning: divide by zero encountered in matmul
self._pcovariance = numpy.matmul(numpy.matmul(
/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/distributions/kernel/baseclass.py:84: RuntimeWarning: overflow encountered in matmul
self._pcovariance = numpy.matmul(numpy.matmul(
/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/distributions/kernel/baseclass.py:84: RuntimeWarning: invalid value encountered in matmul
self._pcovariance = numpy.matmul(numpy.matmul(
[14]:
# Compare the mean of te
results.describe('te', 'mean') - R.describe('te', 'mean')
[14]:
array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])
[15]:
# Compare the te ifirst sobol for Qe_tot
results.raw_data['sobols_first']['te']['Qe_tot'] - R.raw_data['sobols_first']['te']['Qe_tot']
[15]:
array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])
[16]:
# Compare (graphically) the Sobol first results: stabndard and with zero noise
rho_norm = results.describe('rho_norm', 'mean')
plt.clf()
for s in results.sobols_first('te').keys():
plt.plot(rho_norm, results.sobols_first('te')[s], '-', label=f'{s} REF')
for s in R.sobols_first('te').keys():
plt.plot(rho_norm, R.sobols_first('te')[s], '--', label=s)
plt.xlabel('rho_norm')
plt.ylabel('sobols first')
plt.legend(loc=0, ncol=2)
[16]:
<matplotlib.legend.Legend at 0x169778950>
[17]:
# Now perform a campaign over four noise levels with 100 samples each
# We will also restrict the analysis to just 'te' and switch off the calculation/saving of CorrelationMatrices and CorrelationMatrices
analysis = uq.analysis.PCEAnalysis(sampler=old_campaign.get_active_sampler(), qoi_cols=['te'], CorrelationMatrices=False, OutputDistributions=False)
N = 4 * 100
time_start = time.time()
icnt = 0
collect_results={}
for noise in [0.01, 0.02, 0.05, 0.10]:
print(f'{noise = }')
collect_results[noise] = []
for i in range(100):
print(f'{i = }')
df = get_case(old_runs, noise=noise)
R = analysis.analyse(df)
collect_results[noise].append({'df': df, 'R': R})
icnt += 1
print(f'Expected remaining time = {(time.time() - time_start)/icnt*(N-icnt):0.3f}')
time_end = time.time()
print(f'Elapsed time was {time_end-time_start:0.3f} seconds')
noise = 0.01
i = 0
/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/easyvvuq/analysis/pce_analysis.py:548: RuntimeWarning: Skipping computation of cp.Corr
if self.CorrelationMatrices_messages == 0: warnings.warn(f"Skipping computation of cp.Corr", RuntimeWarning) ## print out the warning only once
/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/easyvvuq/analysis/pce_analysis.py:566: RuntimeWarning: Skipping computation of cp.QoI_Dist
if self.OutputDistributions_messages == 0: warnings.warn(f"Skipping computation of cp.QoI_Dist", RuntimeWarning) ## print out the warning only once
Expected remaining time = 354.138
i = 1
Expected remaining time = 374.669
i = 2
Expected remaining time = 364.389
i = 3
Expected remaining time = 359.279
i = 4
Expected remaining time = 356.719
i = 5
Expected remaining time = 353.436
i = 6
Expected remaining time = 351.097
i = 7
Expected remaining time = 348.859
i = 8
Expected remaining time = 347.952
i = 9
Expected remaining time = 346.262
i = 10
Expected remaining time = 344.291
i = 11
Expected remaining time = 342.746
i = 12
Expected remaining time = 341.386
i = 13
Expected remaining time = 342.657
i = 14
Expected remaining time = 341.138
i = 15
Expected remaining time = 340.292
i = 16
Expected remaining time = 339.531
i = 17
Expected remaining time = 338.422
i = 18
Expected remaining time = 337.373
i = 19
Expected remaining time = 336.161
i = 20
Expected remaining time = 334.963
i = 21
Expected remaining time = 334.032
i = 22
Expected remaining time = 332.608
i = 23
Expected remaining time = 331.373
i = 24
Expected remaining time = 330.527
i = 25
Expected remaining time = 330.595
i = 26
Expected remaining time = 329.372
i = 27
Expected remaining time = 328.450
i = 28
Expected remaining time = 327.343
i = 29
Expected remaining time = 326.562
i = 30
Expected remaining time = 325.618
i = 31
Expected remaining time = 324.544
i = 32
Expected remaining time = 323.665
i = 33
Expected remaining time = 322.627
i = 34
Expected remaining time = 321.704
i = 35
Expected remaining time = 320.608
i = 36
Expected remaining time = 319.502
i = 37
Expected remaining time = 319.622
i = 38
Expected remaining time = 318.548
i = 39
Expected remaining time = 317.649
i = 40
Expected remaining time = 316.640
i = 41
Expected remaining time = 315.711
i = 42
Expected remaining time = 314.619
i = 43
Expected remaining time = 313.495
i = 44
Expected remaining time = 312.475
i = 45
Expected remaining time = 311.661
i = 46
Expected remaining time = 310.770
i = 47
Expected remaining time = 310.437
i = 48
Expected remaining time = 309.637
i = 49
Expected remaining time = 308.807
i = 50
Expected remaining time = 308.529
i = 51
Expected remaining time = 307.526
i = 52
Expected remaining time = 306.602
i = 53
Expected remaining time = 305.603
i = 54
Expected remaining time = 304.701
i = 55
Expected remaining time = 303.640
i = 56
Expected remaining time = 302.701
i = 57
Expected remaining time = 302.130
i = 58
Expected remaining time = 301.408
i = 59
Expected remaining time = 300.418
i = 60
Expected remaining time = 299.661
i = 61
Expected remaining time = 298.696
i = 62
Expected remaining time = 298.256
i = 63
Expected remaining time = 297.231
i = 64
Expected remaining time = 296.214
i = 65
Expected remaining time = 295.205
i = 66
Expected remaining time = 294.245
i = 67
Expected remaining time = 293.192
i = 68
Expected remaining time = 292.213
i = 69
Expected remaining time = 291.169
i = 70
Expected remaining time = 290.277
i = 71
Expected remaining time = 289.488
i = 72
Expected remaining time = 288.621
i = 73
Expected remaining time = 287.736
i = 74
Expected remaining time = 286.903
i = 75
Expected remaining time = 286.511
i = 76
Expected remaining time = 285.702
i = 77
Expected remaining time = 284.848
i = 78
Expected remaining time = 283.937
i = 79
Expected remaining time = 283.065
i = 80
Expected remaining time = 282.150
i = 81
Expected remaining time = 281.283
i = 82
Expected remaining time = 280.578
i = 83
Expected remaining time = 279.763
i = 84
Expected remaining time = 279.011
i = 85
Expected remaining time = 278.316
i = 86
Expected remaining time = 277.626
i = 87
Expected remaining time = 277.181
i = 88
Expected remaining time = 276.279
i = 89
Expected remaining time = 275.419
i = 90
Expected remaining time = 274.519
i = 91
Expected remaining time = 273.616
i = 92
Expected remaining time = 272.711
i = 93
Expected remaining time = 271.801
i = 94
Expected remaining time = 270.862
i = 95
Expected remaining time = 270.008
i = 96
Expected remaining time = 269.084
i = 97
Expected remaining time = 268.203
i = 98
Expected remaining time = 267.307
i = 99
Expected remaining time = 266.399
noise = 0.02
i = 0
Expected remaining time = 265.780
i = 1
Expected remaining time = 264.873
i = 2
Expected remaining time = 264.022
i = 3
Expected remaining time = 263.472
i = 4
Expected remaining time = 262.672
i = 5
Expected remaining time = 261.785
i = 6
Expected remaining time = 260.865
i = 7
Expected remaining time = 259.927
i = 8
Expected remaining time = 259.034
i = 9
Expected remaining time = 258.334
i = 10
Expected remaining time = 257.564
i = 11
Expected remaining time = 256.646
i = 12
Expected remaining time = 255.714
i = 13
Expected remaining time = 255.162
i = 14
Expected remaining time = 254.461
i = 15
Expected remaining time = 253.585
i = 16
Expected remaining time = 252.745
i = 17
Expected remaining time = 251.856
i = 18
Expected remaining time = 251.000
i = 19
Expected remaining time = 250.120
i = 20
Expected remaining time = 249.238
i = 21
Expected remaining time = 248.435
i = 22
Expected remaining time = 247.569
i = 23
Expected remaining time = 246.645
i = 24
Expected remaining time = 245.752
i = 25
Expected remaining time = 244.852
i = 26
Expected remaining time = 244.182
i = 27
Expected remaining time = 243.248
i = 28
Expected remaining time = 242.342
i = 29
Expected remaining time = 241.450
i = 30
Expected remaining time = 240.588
i = 31
Expected remaining time = 239.676
i = 32
Expected remaining time = 238.753
i = 33
Expected remaining time = 237.811
i = 34
Expected remaining time = 236.863
i = 35
Expected remaining time = 235.911
i = 36
Expected remaining time = 235.002
i = 37
Expected remaining time = 234.104
i = 38
Expected remaining time = 233.370
i = 39
Expected remaining time = 232.692
i = 40
Expected remaining time = 231.762
i = 41
Expected remaining time = 230.885
i = 42
Expected remaining time = 230.008
i = 43
Expected remaining time = 229.143
i = 44
Expected remaining time = 228.263
i = 45
Expected remaining time = 227.362
i = 46
Expected remaining time = 226.502
i = 47
Expected remaining time = 225.599
i = 48
Expected remaining time = 224.718
i = 49
Expected remaining time = 223.877
i = 50
Expected remaining time = 223.010
i = 51
Expected remaining time = 222.138
i = 52
Expected remaining time = 221.248
i = 53
Expected remaining time = 220.511
i = 54
Expected remaining time = 219.584
i = 55
Expected remaining time = 218.670
i = 56
Expected remaining time = 217.766
i = 57
Expected remaining time = 216.842
i = 58
Expected remaining time = 215.946
i = 59
Expected remaining time = 215.045
i = 60
Expected remaining time = 214.108
i = 61
Expected remaining time = 213.192
i = 62
Expected remaining time = 212.270
i = 63
Expected remaining time = 211.386
i = 64
Expected remaining time = 210.520
i = 65
Expected remaining time = 209.611
i = 66
Expected remaining time = 208.833
i = 67
Expected remaining time = 207.930
i = 68
Expected remaining time = 207.013
i = 69
Expected remaining time = 206.142
i = 70
Expected remaining time = 205.264
i = 71
Expected remaining time = 204.333
i = 72
Expected remaining time = 203.448
i = 73
Expected remaining time = 202.556
i = 74
Expected remaining time = 201.643
i = 75
Expected remaining time = 200.762
i = 76
Expected remaining time = 199.897
i = 77
Expected remaining time = 199.013
i = 78
Expected remaining time = 198.093
i = 79
Expected remaining time = 197.181
i = 80
Expected remaining time = 196.423
i = 81
Expected remaining time = 195.664
i = 82
Expected remaining time = 194.780
i = 83
Expected remaining time = 193.880
i = 84
Expected remaining time = 192.955
i = 85
Expected remaining time = 192.028
i = 86
Expected remaining time = 191.106
i = 87
Expected remaining time = 190.187
i = 88
Expected remaining time = 189.296
i = 89
Expected remaining time = 188.393
i = 90
Expected remaining time = 187.490
i = 91
Expected remaining time = 186.590
i = 92
Expected remaining time = 185.669
i = 93
Expected remaining time = 184.872
i = 94
Expected remaining time = 183.973
i = 95
Expected remaining time = 183.036
i = 96
Expected remaining time = 182.112
i = 97
Expected remaining time = 181.184
i = 98
Expected remaining time = 180.267
i = 99
Expected remaining time = 179.365
noise = 0.05
i = 0
Expected remaining time = 178.448
i = 1
Expected remaining time = 177.539
i = 2
Expected remaining time = 176.642
i = 3
Expected remaining time = 175.739
i = 4
Expected remaining time = 174.856
i = 5
Expected remaining time = 173.972
i = 6
Expected remaining time = 173.035
i = 7
Expected remaining time = 172.213
i = 8
Expected remaining time = 171.315
i = 9
Expected remaining time = 170.443
i = 10
Expected remaining time = 169.542
i = 11
Expected remaining time = 168.638
i = 12
Expected remaining time = 167.727
i = 13
Expected remaining time = 166.806
i = 14
Expected remaining time = 165.899
i = 15
Expected remaining time = 164.987
i = 16
Expected remaining time = 164.081
i = 17
Expected remaining time = 163.173
i = 18
Expected remaining time = 162.255
i = 19
Expected remaining time = 161.354
i = 20
Expected remaining time = 160.450
i = 21
Expected remaining time = 159.635
i = 22
Expected remaining time = 158.738
i = 23
Expected remaining time = 157.830
i = 24
Expected remaining time = 156.912
i = 25
Expected remaining time = 155.997
i = 26
Expected remaining time = 155.104
i = 27
Expected remaining time = 154.195
i = 28
Expected remaining time = 153.285
i = 29
Expected remaining time = 152.373
i = 30
Expected remaining time = 151.460
i = 31
Expected remaining time = 150.548
i = 32
Expected remaining time = 149.642
i = 33
Expected remaining time = 148.734
i = 34
Expected remaining time = 147.831
i = 35
Expected remaining time = 147.015
i = 36
Expected remaining time = 146.105
i = 37
Expected remaining time = 145.185
i = 38
Expected remaining time = 144.277
i = 39
Expected remaining time = 143.376
i = 40
Expected remaining time = 142.473
i = 41
Expected remaining time = 141.566
i = 42
Expected remaining time = 140.661
i = 43
Expected remaining time = 139.768
i = 44
Expected remaining time = 138.861
i = 45
Expected remaining time = 137.951
i = 46
Expected remaining time = 137.048
i = 47
Expected remaining time = 136.189
i = 48
Expected remaining time = 135.287
i = 49
Expected remaining time = 134.437
i = 50
Expected remaining time = 133.533
i = 51
Expected remaining time = 132.623
i = 52
Expected remaining time = 131.718
i = 53
Expected remaining time = 130.817
i = 54
Expected remaining time = 129.911
i = 55
Expected remaining time = 129.007
i = 56
Expected remaining time = 128.101
i = 57
Expected remaining time = 127.204
i = 58
Expected remaining time = 126.308
i = 59
Expected remaining time = 125.405
i = 60
Expected remaining time = 124.504
i = 61
Expected remaining time = 123.600
i = 62
Expected remaining time = 122.697
i = 63
Expected remaining time = 121.794
i = 64
Expected remaining time = 120.951
i = 65
Expected remaining time = 120.042
i = 66
Expected remaining time = 119.140
i = 67
Expected remaining time = 118.238
i = 68
Expected remaining time = 117.338
i = 69
Expected remaining time = 116.445
i = 70
Expected remaining time = 115.538
i = 71
Expected remaining time = 114.634
i = 72
Expected remaining time = 113.732
i = 73
Expected remaining time = 112.830
i = 74
Expected remaining time = 111.924
i = 75
Expected remaining time = 111.035
i = 76
Expected remaining time = 110.138
i = 77
Expected remaining time = 109.247
i = 78
Expected remaining time = 108.393
i = 79
Expected remaining time = 107.490
i = 80
Expected remaining time = 106.579
i = 81
Expected remaining time = 105.676
i = 82
Expected remaining time = 104.772
i = 83
Expected remaining time = 103.862
i = 84
Expected remaining time = 102.956
i = 85
Expected remaining time = 102.052
i = 86
Expected remaining time = 101.149
i = 87
Expected remaining time = 100.249
i = 88
Expected remaining time = 99.344
i = 89
Expected remaining time = 98.440
i = 90
Expected remaining time = 97.536
i = 91
Expected remaining time = 96.632
i = 92
Expected remaining time = 95.735
i = 93
Expected remaining time = 94.881
i = 94
Expected remaining time = 93.977
i = 95
Expected remaining time = 93.076
i = 96
Expected remaining time = 92.174
i = 97
Expected remaining time = 91.275
i = 98
Expected remaining time = 90.376
i = 99
Expected remaining time = 89.473
noise = 0.1
i = 0
Expected remaining time = 88.569
i = 1
Expected remaining time = 87.668
i = 2
Expected remaining time = 86.776
i = 3
Expected remaining time = 85.876
i = 4
Expected remaining time = 84.977
i = 5
Expected remaining time = 84.070
i = 6
Expected remaining time = 83.173
i = 7
Expected remaining time = 82.304
i = 8
Expected remaining time = 81.406
i = 9
Expected remaining time = 80.503
i = 10
Expected remaining time = 79.606
i = 11
Expected remaining time = 78.710
i = 12
Expected remaining time = 77.809
i = 13
Expected remaining time = 76.910
i = 14
Expected remaining time = 76.025
i = 15
Expected remaining time = 75.131
i = 16
Expected remaining time = 74.234
i = 17
Expected remaining time = 73.333
i = 18
Expected remaining time = 72.435
i = 19
Expected remaining time = 71.537
i = 20
Expected remaining time = 70.636
i = 21
Expected remaining time = 69.738
i = 22
Expected remaining time = 68.862
i = 23
Expected remaining time = 67.962
i = 24
Expected remaining time = 67.066
i = 25
Expected remaining time = 66.169
i = 26
Expected remaining time = 65.276
i = 27
Expected remaining time = 64.374
i = 28
Expected remaining time = 63.477
i = 29
Expected remaining time = 62.579
i = 30
Expected remaining time = 61.683
i = 31
Expected remaining time = 60.785
i = 32
Expected remaining time = 59.887
i = 33
Expected remaining time = 58.989
i = 34
Expected remaining time = 58.096
i = 35
Expected remaining time = 57.198
i = 36
Expected remaining time = 56.300
i = 37
Expected remaining time = 55.424
i = 38
Expected remaining time = 54.525
i = 39
Expected remaining time = 53.626
i = 40
Expected remaining time = 52.725
i = 41
Expected remaining time = 51.826
i = 42
Expected remaining time = 50.930
i = 43
Expected remaining time = 50.032
i = 44
Expected remaining time = 49.137
i = 45
Expected remaining time = 48.241
i = 46
Expected remaining time = 47.346
i = 47
Expected remaining time = 46.452
i = 48
Expected remaining time = 45.554
i = 49
Expected remaining time = 44.656
i = 50
Expected remaining time = 43.760
i = 51
Expected remaining time = 42.870
i = 52
Expected remaining time = 41.993
i = 53
Expected remaining time = 41.099
i = 54
Expected remaining time = 40.203
i = 55
Expected remaining time = 39.307
i = 56
Expected remaining time = 38.411
i = 57
Expected remaining time = 37.515
i = 58
Expected remaining time = 36.619
i = 59
Expected remaining time = 35.724
i = 60
Expected remaining time = 34.828
i = 61
Expected remaining time = 33.933
i = 62
Expected remaining time = 33.038
i = 63
Expected remaining time = 32.145
i = 64
Expected remaining time = 31.252
i = 65
Expected remaining time = 30.357
i = 66
Expected remaining time = 29.463
i = 67
Expected remaining time = 28.569
i = 68
Expected remaining time = 27.684
i = 69
Expected remaining time = 26.789
i = 70
Expected remaining time = 25.896
i = 71
Expected remaining time = 25.000
i = 72
Expected remaining time = 24.106
i = 73
Expected remaining time = 23.211
i = 74
Expected remaining time = 22.317
i = 75
Expected remaining time = 21.422
i = 76
Expected remaining time = 20.529
i = 77
Expected remaining time = 19.635
i = 78
Expected remaining time = 18.742
i = 79
Expected remaining time = 17.848
i = 80
Expected remaining time = 16.955
i = 81
Expected remaining time = 16.061
i = 82
Expected remaining time = 15.168
i = 83
Expected remaining time = 14.282
i = 84
Expected remaining time = 13.390
i = 85
Expected remaining time = 12.496
i = 86
Expected remaining time = 11.602
i = 87
Expected remaining time = 10.709
i = 88
Expected remaining time = 9.816
i = 89
Expected remaining time = 8.923
i = 90
Expected remaining time = 8.030
i = 91
Expected remaining time = 7.137
i = 92
Expected remaining time = 6.245
i = 93
Expected remaining time = 5.352
i = 94
Expected remaining time = 4.460
i = 95
Expected remaining time = 3.567
i = 96
Expected remaining time = 2.675
i = 97
Expected remaining time = 1.783
i = 98
Expected remaining time = 0.892
i = 99
Expected remaining time = 0.000
Elapsed time was 356.761 seconds
[18]:
# Save the results
with open('collect_results_100.pickle', "bw") as f_pickle:
pickle.dump(collect_results, f_pickle)
[19]:
# Plot all of the samples, together with the reference values
for noise in collect_results.keys():
plt.figure()
for s in results.sobols_first('te').keys():
plt.plot(rho_norm, results.sobols_first('te')[s], '-', lw=3, label=f'{s} REF')
for c in collect_results[noise]:
for s in results.sobols_first('te').keys():
plt.plot(rho_norm, c['R'].sobols_first('te')[s], alpha=0.5)
plt.xlabel('rho_norm')
plt.ylabel('sobols first')
plt.title(f'{noise = }')
plt.legend(loc=0)
[20]:
# Plot the mean and standard deviation of the perturbed results, together with the original value
for noise in collect_results.keys():
plt.figure()
for s in results.sobols_first('te').keys():
plt.plot(rho_norm, results.sobols_first('te')[s], '--', lw=3, label=f'{s} REF')
for c in collect_results[noise]:
for s in results.sobols_first('te').keys():
V = np.array([c['R'].sobols_first('te')[s] for c in collect_results[noise]])
plt.plot(rho_norm, V.mean(axis=0), '-')
plt.fill_between(rho_norm, V.mean(axis=0) - V.std(axis=0), V.mean(axis=0) + V.std(axis=0), alpha=0.1)
plt.xlabel('rho_norm')
plt.ylabel('sobols first')
plt.title(f'{noise = }')
plt.legend(loc=0)
At higher levels of noise we see that the Sobol first values for H0 drop, particularly in the range 0.5 – 0.9