{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Run the fusion EasyVVUQ campaign using PCE\n",
"\n",
"Run an EasyVVUQ campaign to analyze the sensitivity of the temperature\n",
"profile predicted by a simplified model of heat conduction in a\n",
"tokamak plasma.\n",
"\n",
"This is done with PCE."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"ExecuteTime": {
"end_time": "2021-10-28T07:02:52.718158Z",
"start_time": "2021-10-28T07:02:51.190929Z"
},
"code_folding": [
0
],
"execution": {
"iopub.execute_input": "2025-07-19T11:39:23.493326Z",
"iopub.status.busy": "2025-07-19T11:39:23.493178Z",
"iopub.status.idle": "2025-07-19T11:39:25.228754Z",
"shell.execute_reply": "2025-07-19T11:39:25.228468Z",
"shell.execute_reply.started": "2025-07-19T11:39:23.493311Z"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/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.\n",
" import pkg_resources\n"
]
}
],
"source": [
"# import packages that we will use\n",
"%matplotlib inline\n",
"import os\n",
"import easyvvuq as uq\n",
"import chaospy as cp\n",
"import pickle\n",
"import time\n",
"import numpy as np\n",
"import pandas as pd\n",
"import matplotlib\n",
"if not os.getenv(\"DISPLAY\"): matplotlib.use('Agg')\n",
"import matplotlib.pylab as plt\n",
"from IPython.display import display\n",
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"ExecuteTime": {
"end_time": "2021-10-28T07:02:52.792236Z",
"start_time": "2021-10-28T07:02:52.719309Z"
},
"code_folding": [
0,
2,
4
],
"execution": {
"iopub.execute_input": "2025-07-19T11:39:25.229275Z",
"iopub.status.busy": "2025-07-19T11:39:25.229089Z",
"iopub.status.idle": "2025-07-19T11:39:26.044974Z",
"shell.execute_reply": "2025-07-19T11:39:26.044684Z",
"shell.execute_reply.started": "2025-07-19T11:39:25.229267Z"
}
},
"outputs": [],
"source": [
"# we need fipy -- install if not already available\n",
"\n",
"try:\n",
" import fipy\n",
"except ModuleNotFoundError:\n",
" ! pip install future\n",
" ! pip install fipy\n",
" import fipy"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"ExecuteTime": {
"end_time": "2021-10-28T07:02:52.795743Z",
"start_time": "2021-10-28T07:02:52.793135Z"
},
"code_folding": [
0,
2
],
"execution": {
"iopub.execute_input": "2025-07-19T11:39:26.045466Z",
"iopub.status.busy": "2025-07-19T11:39:26.045398Z",
"iopub.status.idle": "2025-07-19T11:39:26.047417Z",
"shell.execute_reply": "2025-07-19T11:39:26.047196Z",
"shell.execute_reply.started": "2025-07-19T11:39:26.045458Z"
}
},
"outputs": [],
"source": [
"# routine to write out (if needed) the fusion .template file\n",
"\n",
"def write_template(params):\n",
" str = \"\"\n",
" first = True\n",
" for k in params.keys():\n",
" if first:\n",
" str += '{\"%s\": \"$%s\"' % (k,k) ; first = False\n",
" else:\n",
" str += ', \"%s\": \"$%s\"' % (k,k)\n",
" str += '}'\n",
" print(str, file=open('fusion.template','w'))"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"ExecuteTime": {
"end_time": "2021-10-28T07:02:52.800893Z",
"start_time": "2021-10-28T07:02:52.796831Z"
},
"code_folding": [
0
],
"execution": {
"iopub.execute_input": "2025-07-19T11:39:26.047783Z",
"iopub.status.busy": "2025-07-19T11:39:26.047728Z",
"iopub.status.idle": "2025-07-19T11:39:26.049998Z",
"shell.execute_reply": "2025-07-19T11:39:26.049781Z",
"shell.execute_reply.started": "2025-07-19T11:39:26.047778Z"
}
},
"outputs": [],
"source": [
"# define parameters of the fusion model\n",
"def define_params():\n",
" return {\n",
" \"Qe_tot\": {\"type\": \"float\", \"min\": 1.0e6, \"max\": 50.0e6, \"default\": 2e6},\n",
" \"H0\": {\"type\": \"float\", \"min\": 0.00, \"max\": 1.0, \"default\": 0},\n",
" \"Hw\": {\"type\": \"float\", \"min\": 0.01, \"max\": 100.0, \"default\": 0.1},\n",
" \"Te_bc\": {\"type\": \"float\", \"min\": 10.0, \"max\": 1000.0, \"default\": 100},\n",
" \"chi\": {\"type\": \"float\", \"min\": 0.01, \"max\": 100.0, \"default\": 1},\n",
" \"a0\": {\"type\": \"float\", \"min\": 0.2, \"max\": 10.0, \"default\": 1},\n",
" \"R0\": {\"type\": \"float\", \"min\": 0.5, \"max\": 20.0, \"default\": 3},\n",
" \"E0\": {\"type\": \"float\", \"min\": 1.0, \"max\": 10.0, \"default\": 1.5},\n",
" \"b_pos\": {\"type\": \"float\", \"min\": 0.95, \"max\": 0.99, \"default\": 0.98},\n",
" \"b_height\": {\"type\": \"float\", \"min\": 3e19, \"max\": 10e19, \"default\": 6e19},\n",
" \"b_sol\": {\"type\": \"float\", \"min\": 2e18, \"max\": 3e19, \"default\": 2e19},\n",
" \"b_width\": {\"type\": \"float\", \"min\": 0.005, \"max\": 0.025, \"default\": 0.01},\n",
" \"b_slope\": {\"type\": \"float\", \"min\": 0.0, \"max\": 0.05, \"default\": 0.01},\n",
" \"nr\": {\"type\": \"integer\", \"min\": 10, \"max\": 1000, \"default\": 100},\n",
" \"dt\": {\"type\": \"float\", \"min\": 1e-3, \"max\": 1e3, \"default\": 100},\n",
" \"out_file\": {\"type\": \"string\", \"default\": \"output.csv\"}\n",
" }"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"ExecuteTime": {
"end_time": "2021-10-28T07:02:52.806065Z",
"start_time": "2021-10-28T07:02:52.801728Z"
},
"code_folding": [
0,
2,
17,
21,
28
],
"execution": {
"iopub.execute_input": "2025-07-19T11:39:26.050439Z",
"iopub.status.busy": "2025-07-19T11:39:26.050353Z",
"iopub.status.idle": "2025-07-19T11:39:26.052767Z",
"shell.execute_reply": "2025-07-19T11:39:26.052538Z",
"shell.execute_reply.started": "2025-07-19T11:39:26.050432Z"
}
},
"outputs": [],
"source": [
"# define varying quantities\n",
"def define_vary():\n",
" vary_all = {\n",
" \"Qe_tot\": cp.Uniform(1.8e6, 2.2e6),\n",
" \"H0\": cp.Uniform(0.0, 0.2),\n",
" \"Hw\": cp.Uniform(0.1, 0.5),\n",
" \"chi\": cp.Uniform(0.8, 1.2),\n",
" \"Te_bc\": cp.Uniform(80.0, 120.0),\n",
" \"a0\": cp.Uniform(0.9, 1.1),\n",
" \"R0\": cp.Uniform(2.7, 3.3),\n",
" \"E0\": cp.Uniform(1.4, 1.6),\n",
" \"b_pos\": cp.Uniform(0.95, 0.99),\n",
" \"b_height\": cp.Uniform(5e19, 7e19),\n",
" \"b_sol\": cp.Uniform(1e19, 3e19),\n",
" \"b_width\": cp.Uniform(0.015, 0.025),\n",
" \"b_slope\": cp.Uniform(0.005, 0.020)\n",
" }\n",
" vary_2 = {\n",
" \"Qe_tot\": cp.Uniform(1.8e6, 2.2e6),\n",
" \"Te_bc\": cp.Uniform(80.0, 120.0)\n",
" }\n",
" vary_5 = {\n",
" \"Qe_tot\": cp.Uniform(1.8e6, 2.2e6),\n",
" \"H0\": cp.Uniform(0.0, 0.2),\n",
" \"Hw\": cp.Uniform(0.1, 0.5),\n",
" \"chi\": cp.Uniform(0.8, 1.2),\n",
" \"Te_bc\": cp.Uniform(80.0, 120.0)\n",
" }\n",
" vary_10 = {\n",
" \"Qe_tot\": cp.Uniform(1.8e6, 2.2e6),\n",
" \"H0\": cp.Uniform(0.0, 0.2),\n",
" \"Hw\": cp.Uniform(0.1, 0.5),\n",
" \"chi\": cp.Uniform(0.8, 1.2),\n",
" \"Te_bc\": cp.Uniform(80.0, 120.0),\n",
" \"b_pos\": cp.Uniform(0.95, 0.99),\n",
" \"b_height\": cp.Uniform(5e19, 7e19),\n",
" \"b_sol\": cp.Uniform(1e19, 3e19),\n",
" \"b_width\": cp.Uniform(0.015, 0.025),\n",
" \"b_slope\": cp.Uniform(0.005, 0.020)\n",
" }\n",
" return vary_5"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"ExecuteTime": {
"end_time": "2021-10-28T07:02:52.810207Z",
"start_time": "2021-10-28T07:02:52.806896Z"
},
"code_folding": [
0
],
"execution": {
"iopub.execute_input": "2025-07-19T11:39:26.055757Z",
"iopub.status.busy": "2025-07-19T11:39:26.055570Z",
"iopub.status.idle": "2025-07-19T11:39:26.058208Z",
"shell.execute_reply": "2025-07-19T11:39:26.057767Z",
"shell.execute_reply.started": "2025-07-19T11:39:26.055746Z"
}
},
"outputs": [],
"source": [
"# define a model to run the fusion code directly from python, expecting a dictionary and returning a dictionary\n",
"def run_fusion_model(input):\n",
" import json\n",
" import fusion\n",
" qois = [\"te\", \"ne\", \"rho\", \"rho_norm\"]\n",
" del input['out_file']\n",
" return {q: v for q,v in zip(qois, [t.tolist() for t in fusion.solve_Te(**input, plots=False, output=False)])}"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"ExecuteTime": {
"end_time": "2021-10-28T07:02:52.819459Z",
"start_time": "2021-10-28T07:02:52.811176Z"
},
"code_folding": [
0
],
"execution": {
"iopub.execute_input": "2025-07-19T11:39:26.059039Z",
"iopub.status.busy": "2025-07-19T11:39:26.058774Z",
"iopub.status.idle": "2025-07-19T11:39:26.063378Z",
"shell.execute_reply": "2025-07-19T11:39:26.063098Z",
"shell.execute_reply.started": "2025-07-19T11:39:26.059032Z"
}
},
"outputs": [],
"source": [
"# routine to run a PCE campaign\n",
"\n",
"def run_pce_case(pce_order=2, local=True, dask=True, batch_size=os.cpu_count(), use_files=True):\n",
" \"\"\"\n",
" Inputs:\n",
" pce_order: order of the PCE expansion\n",
" local: if using Dask, whether to use the local option (True)\n",
" dask: whether to use dask (True)\n",
" batch_size: for the non Dask option, number of cases to run in parallel (16)\n",
" Outputs:\n",
" results_df: Pandas dataFrame containing inputs to and output from the model\n",
" results: Results of the PCE analysis\n",
" times: Information about the elapsed time for the various phases of the calculation\n",
" pce_order: pce_order \n",
" count: number of PCE samples\n",
" \"\"\"\n",
" \n",
" if dask:\n",
" if local:\n",
" print('Running locally')\n",
" import multiprocessing.popen_spawn_posix\n",
" from dask.distributed import Client, LocalCluster\n",
" cluster = LocalCluster(threads_per_worker=1)\n",
" client = Client(cluster) # processes=True, threads_per_worker=1)\n",
" else:\n",
" print('Running using SLURM')\n",
" from dask.distributed import Client\n",
" from dask_jobqueue import SLURMCluster\n",
" cluster = SLURMCluster(\n",
" job_extra=['--qos=p.tok.openmp.2h', '--mail-type=end', '--mail-user=dpc@rzg.mpg.de', '-t 2:00:00'], \n",
" queue='p.tok.openmp', \n",
" cores=8, \n",
" memory='8 GB',\n",
" processes=8)\n",
" cluster.scale(32)\n",
" print(cluster)\n",
" print(cluster.job_script())\n",
" client = Client(cluster)\n",
" print(client)\n",
"\n",
" else:\n",
" import concurrent.futures\n",
"# client = concurrent.futures.ProcessPoolExecutor(max_workers=batch_size)\n",
" client = concurrent.futures.ThreadPoolExecutor(max_workers=batch_size)\n",
"# client = None\n",
" \n",
" times = np.zeros(7)\n",
"\n",
" time_start = time.time()\n",
" time_start_whole = time_start\n",
" # Set up a fresh campaign called \"fusion_pce.\"\n",
" my_campaign = uq.Campaign(name='fusion_pce.') \n",
"\n",
" # Define parameter space\n",
" params = define_params()\n",
"\n",
" # Create an encoder and decoder for PCE test app\n",
" if use_files:\n",
" encoder = uq.encoders.GenericEncoder(template_fname='fusion.template',\n",
" delimiter='$',\n",
" target_filename='fusion_in.json')\n",
"\n",
"\n",
" decoder = uq.decoders.SimpleCSV(target_filename=\"output.csv\",\n",
" output_columns=[\"te\", \"ne\", \"rho\", \"rho_norm\"])\n",
"\n",
" execute = uq.actions.ExecuteLocal('python3 %s/fusion_model.py fusion_in.json' % (os.getcwd()))\n",
"\n",
" actions = uq.actions.Actions(uq.actions.CreateRunDirectory('/tmp'), \n",
" uq.actions.Encode(encoder), execute, uq.actions.Decode(decoder))\n",
" else:\n",
" actions = uq.actions.Actions(uq.actions.ExecutePython(run_fusion_model))\n",
"\n",
"\n",
" # Add the app (automatically set as current app)\n",
" my_campaign.add_app(name=\"fusion\", params=params, actions=actions)\n",
"\n",
" time_end = time.time()\n",
" times[1] = time_end-time_start\n",
" print('Time for phase 1 = %.3f' % (times[1]))\n",
"\n",
" time_start = time.time()\n",
" # Associate a sampler with the campaign\n",
" my_campaign.set_sampler(uq.sampling.PCESampler(vary=define_vary(), polynomial_order=pce_order))\n",
" my_campaign.draw_samples()\n",
" print('Number of samples = %s' % my_campaign.get_active_sampler().count)\n",
"\n",
" time_end = time.time()\n",
" times[2] = time_end-time_start\n",
" print('Time for phase 2 = %.3f' % (times[2]))\n",
"\n",
" time_start = time.time()\n",
" # Perform the actions\n",
" my_campaign.execute(pool=client).collate(progress_bar=True)\n",
"\n",
" if dask:\n",
" client.close()\n",
" client.shutdown()\n",
"\n",
" time_end = time.time()\n",
" times[3] = time_end-time_start\n",
" print('Time for phase 3 = %.3f' % (times[3]))\n",
"\n",
" time_start = time.time()\n",
" # Collate the results\n",
" results_df = my_campaign.get_collation_result()\n",
"\n",
" time_end = time.time()\n",
" times[4] = time_end-time_start\n",
" print('Time for phase 4 = %.3f' % (times[4]))\n",
"\n",
" time_start = time.time()\n",
" # Post-processing analysis\n",
" results = my_campaign.analyse(qoi_cols=[\"te\", \"ne\", \"rho\", \"rho_norm\"])\n",
"\n",
" time_end = time.time()\n",
" times[5] = time_end-time_start\n",
" print('Time for phase 5 = %.3f' % (times[5]))\n",
"\n",
" time_start = time.time()\n",
" # Save the results\n",
" pickle.dump(results, open('fusion_results.pickle','bw'))\n",
" time_end = time.time()\n",
" times[6] = time_end-time_start\n",
" print('Time for phase 6 = %.3f' % (times[6]))\n",
"\n",
" times[0] = time_end - time_start_whole\n",
"\n",
" return results_df, results, times, pce_order, my_campaign.get_active_sampler().count"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"ExecuteTime": {
"end_time": "2021-10-28T07:02:52.833748Z",
"start_time": "2021-10-28T07:02:52.820385Z"
},
"code_folding": [
0,
2,
21,
40,
52,
66,
78
],
"execution": {
"iopub.execute_input": "2025-07-19T11:39:26.064030Z",
"iopub.status.busy": "2025-07-19T11:39:26.063945Z",
"iopub.status.idle": "2025-07-19T11:39:26.070102Z",
"shell.execute_reply": "2025-07-19T11:39:26.069861Z",
"shell.execute_reply.started": "2025-07-19T11:39:26.064024Z"
}
},
"outputs": [],
"source": [
"# routines for plotting the results\n",
"\n",
"def plot_Te(results, title=None):\n",
" # plot the calculated Te: mean, with std deviation, 1, 10, 90 and 99%\n",
" plt.figure()\n",
" rho = results.describe('rho', 'mean')\n",
" plt.plot(rho, results.describe('te', 'mean'), 'b-', label='Mean')\n",
" plt.plot(rho, results.describe('te', 'mean')-results.describe('te', 'std'), 'b--', label='+1 std deviation')\n",
" plt.plot(rho, results.describe('te', 'mean')+results.describe('te', 'std'), 'b--')\n",
" plt.fill_between(rho, results.describe('te', 'mean')-results.describe('te', 'std'), results.describe('te', 'mean')+results.describe('te', 'std'), color='b', alpha=0.2)\n",
" plt.plot(rho, results.describe('te', '10%'), 'b:', label='10 and 90 percentiles')\n",
" plt.plot(rho, results.describe('te', '90%'), 'b:')\n",
" plt.fill_between(rho, results.describe('te', '10%'), results.describe('te', '90%'), color='b', alpha=0.1)\n",
" plt.fill_between(rho, results.describe('te', '1%'), results.describe('te', '99%'), color='b', alpha=0.05)\n",
" plt.legend(loc=0)\n",
" plt.xlabel('rho [$m$]')\n",
" plt.ylabel('Te [$eV$]')\n",
" if not title is None: plt.title(title)\n",
" plt.savefig('Te.png')\n",
" plt.savefig('Te.pdf')\n",
"\n",
"def plot_ne(results, title=None):\n",
" # plot the calculated ne: mean, with std deviation, 1, 10, 90 and 99%\n",
" plt.figure()\n",
" rho = results.describe('rho', 'mean')\n",
" plt.plot(rho, results.describe('ne', 'mean'), 'b-', label='Mean')\n",
" plt.plot(rho, results.describe('ne', 'mean')-results.describe('ne', 'std'), 'b--', label='+1 std deviation')\n",
" plt.plot(rho, results.describe('ne', 'mean')+results.describe('ne', 'std'), 'b--')\n",
" plt.fill_between(rho, results.describe('ne', 'mean')-results.describe('ne', 'std'), results.describe('ne', 'mean')+results.describe('ne', 'std'), color='b', alpha=0.2)\n",
" plt.plot(rho, results.describe('ne', '10%'), 'b:', label='10 and 90 percentiles')\n",
" plt.plot(rho, results.describe('ne', '90%'), 'b:')\n",
" plt.fill_between(rho, results.describe('ne', '10%'), results.describe('ne', '90%'), color='b', alpha=0.1)\n",
" plt.fill_between(rho, results.describe('ne', '1%'), results.describe('ne', '99%'), color='b', alpha=0.05)\n",
" plt.legend(loc=0)\n",
" plt.xlabel('rho [$m$]')\n",
" plt.ylabel('ne [$m^{-3}$]')\n",
" if not title is None: plt.title(title)\n",
" plt.savefig('ne.png')\n",
" plt.savefig('ne.pdf')\n",
"\n",
"def plot_sobols_first(results, title=None, field='te'):\n",
" # plot the first Sobol results\n",
" plt.figure()\n",
" rho = results.describe('rho', 'mean')\n",
" for k in results.sobols_first()[field].keys(): plt.plot(rho, results.sobols_first()[field][k], label=k)\n",
" plt.legend(loc=0)\n",
" plt.xlabel('rho [$m$]')\n",
" plt.ylabel('sobols_first')\n",
" if not title is None: plt.title(field + ': ' + title)\n",
" plt.savefig('sobols_first_%s.png' % (field))\n",
" plt.savefig('sobols_first_%s.pdf' % (field))\n",
"\n",
"def plot_sobols_second(results, title=None, field='te'):\n",
" # plot the second Sobol results\n",
" plt.figure()\n",
" rho = results.describe('rho', 'mean')\n",
" for k1 in results.sobols_second()[field].keys():\n",
" for k2 in results.sobols_second()[field][k1].keys():\n",
" plt.plot(rho, results.sobols_second()[field][k1][k2], label=k1+'/'+k2)\n",
" plt.legend(loc=0, ncol=2)\n",
" plt.xlabel('rho [$m$]')\n",
" plt.ylabel('sobols_second')\n",
" if not title is None: plt.title(field + ': ' + title)\n",
" plt.savefig('sobols_second_%s.png' % (field))\n",
" plt.savefig('sobols_second_%s.pdf' % (field))\n",
"\n",
"def plot_sobols_total(results, title=None, field='te'):\n",
" # plot the total Sobol results\n",
" plt.figure()\n",
" rho = results.describe('rho', 'mean')\n",
" for k in results.sobols_total()[field].keys(): plt.plot(rho, results.sobols_total()[field][k], label=k)\n",
" plt.legend(loc=0)\n",
" plt.xlabel('rho [$m$]')\n",
" plt.ylabel('sobols_total')\n",
" if not title is None: plt.title(field + ': ' + title)\n",
" plt.savefig('sobols_total_%s.png' % (field))\n",
" plt.savefig('sobols_total_%s.pdf' % (field))\n",
"\n",
"def plot_distribution(results, results_df, title=None):\n",
" te_dist = results.raw_data['output_distributions']['te']\n",
" rho_norm = results.describe('rho_norm', 'mean')\n",
" plt.subplots(3,3,figsize=(12,12))\n",
" ip=0\n",
" for i in [np.maximum(0, int(i-1)) \n",
" for i in np.linspace(0,1,9) * rho_norm.shape]:\n",
" ip += 1\n",
" plt.subplot(3,3,ip)\n",
" pdf_raw_samples = cp.GaussianKDE(results_df.te[i])\n",
" pdf_kde_samples = cp.GaussianKDE(te_dist.samples[i])\n",
" plt.hist(results_df.te[i], density=True, bins=50, label='histogram of raw samples', alpha=0.25)\n",
" if hasattr(te_dist, 'samples'):\n",
" plt.hist(te_dist.samples[i], density=True, bins=50, label='histogram of kde samples', alpha=0.25)\n",
"\n",
" plt.plot(np.linspace(pdf_raw_samples.lower, pdf_raw_samples.upper), pdf_raw_samples.pdf(np.linspace(pdf_raw_samples.lower, pdf_raw_samples.upper)), label='PDF (raw samples)')\n",
" plt.plot(np.linspace(pdf_kde_samples.lower, pdf_kde_samples.upper), pdf_kde_samples.pdf(np.linspace(pdf_kde_samples.lower, pdf_kde_samples.upper)), label='PDF (kde samples)')\n",
"\n",
" plt.legend(loc=0)\n",
" plt.xlabel('Te [$eV$]')\n",
" plt.title('Distributions for rho_norm = %0.4f' % (rho_norm[i]))\n",
" plt.subplots_adjust(left=0.05, right=0.95, bottom=0.05, top=0.925, wspace=0.4, hspace=0.3)\n",
" if not title is None: plt.suptitle(title)\n",
" plt.savefig('distribution_function.png')\n",
" plt.savefig('distribution_function.pdf')"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"ExecuteTime": {
"end_time": "2021-10-28T07:31:39.974149Z",
"start_time": "2021-10-28T07:02:52.834544Z"
},
"code_folding": [
0
],
"execution": {
"iopub.execute_input": "2025-07-19T11:39:26.070560Z",
"iopub.status.busy": "2025-07-19T11:39:26.070478Z",
"iopub.status.idle": "2025-07-19T12:07:43.800749Z",
"shell.execute_reply": "2025-07-19T12:07:43.800447Z",
"shell.execute_reply.started": "2025-07-19T11:39:26.070553Z"
},
"scrolled": true
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/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'}\n",
" warn(\n",
"/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'}\n",
" warn(\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Time for phase 1 = 0.018\n",
"Number of samples = 32\n",
"Time for phase 2 = 0.049\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|███████████████████████████████████████████| 32/32 [00:00<00:00, 72.46it/s]\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Time for phase 3 = 0.457\n",
"Time for phase 4 = 0.006\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/distributions/kernel/baseclass.py:84: RuntimeWarning: divide by zero encountered in matmul\n",
" self._pcovariance = numpy.matmul(numpy.matmul(\n",
"/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/distributions/kernel/baseclass.py:84: RuntimeWarning: overflow encountered in matmul\n",
" self._pcovariance = numpy.matmul(numpy.matmul(\n",
"/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/distributions/kernel/baseclass.py:84: RuntimeWarning: invalid value encountered in matmul\n",
" self._pcovariance = numpy.matmul(numpy.matmul(\n",
"Traceback (most recent call last):\n",
" File \"/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/easyvvuq/analysis/pce_analysis.py\", line 498, in analyse\n",
" dY_hat = build_surrogate_der(fit, verbose=False)\n",
" ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n",
" File \"/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/easyvvuq/analysis/pce_analysis.py\", line 338, in build_surrogate_der\n",
" assert(n1 == n2)\n",
" ^^^^^^^^\n",
"AssertionError\n",
"Traceback (most recent call last):\n",
" File \"/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/easyvvuq/analysis/pce_analysis.py\", line 498, in analyse\n",
" dY_hat = build_surrogate_der(fit, verbose=False)\n",
" ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n",
" File \"/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/easyvvuq/analysis/pce_analysis.py\", line 338, in build_surrogate_der\n",
" assert(n1 == n2)\n",
" ^^^^^^^^\n",
"AssertionError\n",
"Traceback (most recent call last):\n",
" File \"/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/easyvvuq/analysis/pce_analysis.py\", line 498, in analyse\n",
" dY_hat = build_surrogate_der(fit, verbose=False)\n",
" ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n",
" File \"/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/easyvvuq/analysis/pce_analysis.py\", line 338, in build_surrogate_der\n",
" assert(n1 == n2)\n",
" ^^^^^^^^\n",
"AssertionError\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Time for phase 5 = 0.726\n",
"Time for phase 6 = 0.014\n",
"Time for phase 1 = 0.005\n",
"Number of samples = 243\n",
"Time for phase 2 = 0.139\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|█████████████████████████████████████████| 243/243 [00:03<00:00, 78.23it/s]\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Time for phase 3 = 3.185\n",
"Time for phase 4 = 0.025\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/distributions/kernel/baseclass.py:84: RuntimeWarning: divide by zero encountered in matmul\n",
" self._pcovariance = numpy.matmul(numpy.matmul(\n",
"/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/distributions/kernel/baseclass.py:84: RuntimeWarning: overflow encountered in matmul\n",
" self._pcovariance = numpy.matmul(numpy.matmul(\n",
"/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/distributions/kernel/baseclass.py:84: RuntimeWarning: invalid value encountered in matmul\n",
" self._pcovariance = numpy.matmul(numpy.matmul(\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Time for phase 5 = 3.443\n",
"Time for phase 6 = 0.016\n",
"Time for phase 1 = 0.005\n",
"Number of samples = 1024\n",
"Time for phase 2 = 0.515\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|███████████████████████████████████████| 1024/1024 [00:14<00:00, 72.49it/s]\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Time for phase 3 = 14.287\n",
"Time for phase 4 = 0.114\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/distributions/kernel/baseclass.py:84: RuntimeWarning: divide by zero encountered in matmul\n",
" self._pcovariance = numpy.matmul(numpy.matmul(\n",
"/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/distributions/kernel/baseclass.py:84: RuntimeWarning: overflow encountered in matmul\n",
" self._pcovariance = numpy.matmul(numpy.matmul(\n",
"/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/distributions/kernel/baseclass.py:84: RuntimeWarning: invalid value encountered in matmul\n",
" self._pcovariance = numpy.matmul(numpy.matmul(\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Time for phase 5 = 13.609\n",
"Time for phase 6 = 0.019\n",
"Time for phase 1 = 0.009\n",
"Number of samples = 3125\n",
"Time for phase 2 = 1.334\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|███████████████████████████████████████| 3125/3125 [00:43<00:00, 71.03it/s]\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Time for phase 3 = 44.234\n",
"Time for phase 4 = 0.332\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/distributions/kernel/baseclass.py:84: RuntimeWarning: divide by zero encountered in matmul\n",
" self._pcovariance = numpy.matmul(numpy.matmul(\n",
"/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/distributions/kernel/baseclass.py:84: RuntimeWarning: overflow encountered in matmul\n",
" self._pcovariance = numpy.matmul(numpy.matmul(\n",
"/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/distributions/kernel/baseclass.py:84: RuntimeWarning: invalid value encountered in matmul\n",
" self._pcovariance = numpy.matmul(numpy.matmul(\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Time for phase 5 = 60.253\n",
"Time for phase 6 = 0.026\n",
"Time for phase 1 = 0.014\n",
"Number of samples = 7776\n",
"Time for phase 2 = 2.977\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|███████████████████████████████████████| 7776/7776 [01:43<00:00, 75.37it/s]\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Time for phase 3 = 103.546\n",
"Time for phase 4 = 0.885\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/distributions/kernel/baseclass.py:84: RuntimeWarning: divide by zero encountered in matmul\n",
" self._pcovariance = numpy.matmul(numpy.matmul(\n",
"/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/distributions/kernel/baseclass.py:84: RuntimeWarning: overflow encountered in matmul\n",
" self._pcovariance = numpy.matmul(numpy.matmul(\n",
"/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/distributions/kernel/baseclass.py:84: RuntimeWarning: invalid value encountered in matmul\n",
" self._pcovariance = numpy.matmul(numpy.matmul(\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Time for phase 5 = 255.867\n",
"Time for phase 6 = 0.033\n",
"Time for phase 1 = 0.010\n",
"Number of samples = 16807\n",
"Time for phase 2 = 6.281\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|█████████████████████████████████████| 16807/16807 [03:40<00:00, 76.19it/s]\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Time for phase 3 = 221.146\n",
"Time for phase 4 = 2.180\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/distributions/kernel/baseclass.py:84: RuntimeWarning: divide by zero encountered in matmul\n",
" self._pcovariance = numpy.matmul(numpy.matmul(\n",
"/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/distributions/kernel/baseclass.py:84: RuntimeWarning: overflow encountered in matmul\n",
" self._pcovariance = numpy.matmul(numpy.matmul(\n",
"/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/distributions/kernel/baseclass.py:84: RuntimeWarning: invalid value encountered in matmul\n",
" self._pcovariance = numpy.matmul(numpy.matmul(\n",
"/Volumes/UserData/dpc/GIT/EasyVVUQ/env/lib/python3.12/site-packages/chaospy/descriptives/correlation/pearson.py:45: RuntimeWarning: invalid value encountered in sqrt\n",
" vvar = numpy.sqrt(numpy.outer(var, var))\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Time for phase 5 = 961.891\n",
"Time for phase 6 = 0.065\n"
]
}
],
"source": [
"# Calculate the polynomial chaos expansion for a range of orders\n",
"\n",
"if __name__ == '__main__':\n",
" local = False # if True, use local cores; if False, use SLURM\n",
" dask = False # if True, use DASK; if False, use a fall-back non-DASK option\n",
"\n",
" R = {}\n",
" for pce_order in range(1, 7):\n",
" R[pce_order] = {}\n",
" (R[pce_order]['results_df'], \n",
" R[pce_order]['results'], \n",
" R[pce_order]['times'], \n",
" R[pce_order]['order'], \n",
" R[pce_order]['number_of_samples']) = run_pce_case(pce_order=pce_order, \n",
" local=local, dask=dask, \n",
" batch_size=7, use_files=False)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"ExecuteTime": {
"end_time": "2021-10-28T07:31:40.668315Z",
"start_time": "2021-10-28T07:31:39.989788Z"
},
"code_folding": [
0
],
"execution": {
"iopub.execute_input": "2025-07-19T12:07:43.802034Z",
"iopub.status.busy": "2025-07-19T12:07:43.801895Z",
"iopub.status.idle": "2025-07-19T12:07:44.110898Z",
"shell.execute_reply": "2025-07-19T12:07:44.110437Z",
"shell.execute_reply.started": "2025-07-19T12:07:43.802023Z"
}
},
"outputs": [],
"source": [
"# save the results\n",
"\n",
"if __name__ == '__main__':\n",
"\n",
" pickle.dump(R, open('collected_results.pickle','bw'))"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"ExecuteTime": {
"end_time": "2021-10-28T07:31:40.693513Z",
"start_time": "2021-10-28T07:31:40.669193Z"
},
"code_folding": [
0
],
"execution": {
"iopub.execute_input": "2025-07-19T12:07:44.111451Z",
"iopub.status.busy": "2025-07-19T12:07:44.111368Z",
"iopub.status.idle": "2025-07-19T12:07:44.127516Z",
"shell.execute_reply": "2025-07-19T12:07:44.127246Z",
"shell.execute_reply.started": "2025-07-19T12:07:44.111445Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" Total | \n",
" Phase 1 | \n",
" Phase 2 | \n",
" Phase 3 | \n",
" Phase 4 | \n",
" Phase 5 | \n",
" Phase 6 | \n",
"
\n",
" \n",
" \n",
" \n",
" | 1 | \n",
" 1.270956 | \n",
" 0.018249 | \n",
" 0.049477 | \n",
" 0.457251 | \n",
" 0.006136 | \n",
" 0.725576 | \n",
" 0.014101 | \n",
"
\n",
" \n",
" | 2 | \n",
" 6.813447 | \n",
" 0.004716 | \n",
" 0.139352 | \n",
" 3.184508 | \n",
" 0.024941 | \n",
" 3.443412 | \n",
" 0.016367 | \n",
"
\n",
" \n",
" | 3 | \n",
" 28.555106 | \n",
" 0.005205 | \n",
" 0.515448 | \n",
" 14.287420 | \n",
" 0.113977 | \n",
" 13.609176 | \n",
" 0.019068 | \n",
"
\n",
" \n",
" | 4 | \n",
" 106.189113 | \n",
" 0.009365 | \n",
" 1.334293 | \n",
" 44.234357 | \n",
" 0.332101 | \n",
" 60.252722 | \n",
" 0.025862 | \n",
"
\n",
" \n",
" | 5 | \n",
" 363.321228 | \n",
" 0.013717 | \n",
" 2.976511 | \n",
" 103.546055 | \n",
" 0.884662 | \n",
" 255.866825 | \n",
" 0.032947 | \n",
"
\n",
" \n",
" | 6 | \n",
" 1191.574097 | \n",
" 0.010267 | \n",
" 6.280999 | \n",
" 221.146430 | \n",
" 2.179955 | \n",
" 961.891130 | \n",
" 0.064794 | \n",
"
\n",
" \n",
"
\n",
"