{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Sensityivity analysis for the Ishigama function with noise using PCE\n",
"\n",
"Run an EasyVVUQ campaign to analyze the sensitivity for the Ishigami function with noise\n",
"\n",
"This is done with PCE providing a normal distributed noise value"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"ExecuteTime": {
"end_time": "2021-12-10T08:41:19.729525Z",
"start_time": "2021-12-10T08:41:14.617602Z"
},
"code_folding": [],
"execution": {
"iopub.execute_input": "2025-07-18T12:08:01.062513Z",
"iopub.status.busy": "2025-07-18T12:08:01.062304Z",
"iopub.status.idle": "2025-07-18T12:08:19.139957Z",
"shell.execute_reply": "2025-07-18T12:08:19.139397Z",
"shell.execute_reply.started": "2025-07-18T12:08:01.062502Z"
}
},
"outputs": [],
"source": [
"# Run an EasyVVUQ campaign to analyze the sensitivity for the Ishigami function with noise\n",
"# This is done with PCE providing a normal distributed noise value.\n",
"%matplotlib inline\n",
"import os\n",
"import easyvvuq as uq\n",
"import chaospy as cp\n",
"import pickle\n",
"import numpy as np\n",
"import matplotlib.pylab as plt\n",
"import time\n",
"import pandas as pd"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"ExecuteTime": {
"end_time": "2021-12-10T08:41:19.734097Z",
"start_time": "2021-12-10T08:41:19.730725Z"
},
"code_folding": [
0,
1
],
"execution": {
"iopub.execute_input": "2025-07-18T12:08:19.141306Z",
"iopub.status.busy": "2025-07-18T12:08:19.140478Z",
"iopub.status.idle": "2025-07-18T12:08:19.148294Z",
"shell.execute_reply": "2025-07-18T12:08:19.147965Z",
"shell.execute_reply.started": "2025-07-18T12:08:19.141242Z"
}
},
"outputs": [],
"source": [
"# Define the Ishigami function\n",
"def ishigamiSA(a,b):\n",
" '''Exact sensitivity indices of the Ishigami function for given a and b.\n",
" From https://openturns.github.io/openturns/master/examples/meta_modeling/chaos_ishigami.html\n",
" '''\n",
" var = 1.0/2 + a**2/8 + b*np.pi**4/5 + b**2*np.pi**8/18\n",
" S1 = (1.0/2 + b*np.pi**4/5+b**2*np.pi**8/50)/var\n",
" S2 = (a**2/8)/var\n",
" S3 = 0\n",
" S13 = b**2*np.pi**8/2*(1.0/9-1.0/25)/var\n",
" exact = {\n",
" 'expectation' : a/2,\n",
" 'variance' : var,\n",
" 'S1' : (1.0/2 + b*np.pi**4/5+b**2*np.pi**8.0/50)/var,\n",
" 'S2' : (a**2/8)/var,\n",
" 'S3' : 0,\n",
" 'S12' : 0,\n",
" 'S23' : 0,\n",
" 'S13' : S13,\n",
" 'S123' : 0,\n",
" 'ST1' : S1 + S13,\n",
" 'ST2' : S2,\n",
" 'ST3' : S3 + S13\n",
" }\n",
" return exact\n",
"\n",
"Ishigami_a = 7.0\n",
"Ishigami_b = 0.1\n",
"exact = ishigamiSA(Ishigami_a, Ishigami_b)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"ExecuteTime": {
"end_time": "2021-12-10T08:41:19.736546Z",
"start_time": "2021-12-10T08:41:19.734786Z"
},
"code_folding": [
0
],
"execution": {
"iopub.execute_input": "2025-07-18T12:08:19.148956Z",
"iopub.status.busy": "2025-07-18T12:08:19.148692Z",
"iopub.status.idle": "2025-07-18T12:08:19.151590Z",
"shell.execute_reply": "2025-07-18T12:08:19.151321Z",
"shell.execute_reply.started": "2025-07-18T12:08:19.148927Z"
}
},
"outputs": [],
"source": [
"# define a model to run the Ishigami code directly from python, expecting a dictionary and returning a dictionary\n",
"def run_ishigami_model(input):\n",
" import Ishigami\n",
" qois = [\"Ishigami\"]\n",
" del input['out_file']\n",
" N = input['N']\n",
" del input['N']\n",
" return {qois[0]: Ishigami.evaluate(**input)+N}"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"ExecuteTime": {
"end_time": "2021-12-10T08:41:19.740153Z",
"start_time": "2021-12-10T08:41:19.737318Z"
},
"code_folding": [
0
],
"execution": {
"iopub.execute_input": "2025-07-18T12:08:19.152118Z",
"iopub.status.busy": "2025-07-18T12:08:19.151997Z",
"iopub.status.idle": "2025-07-18T12:08:19.154201Z",
"shell.execute_reply": "2025-07-18T12:08:19.153972Z",
"shell.execute_reply.started": "2025-07-18T12:08:19.152110Z"
}
},
"outputs": [],
"source": [
"# Define parameter space\n",
"def define_params():\n",
" return {\n",
" \"x1\": {\"type\": \"float\", \"min\": -np.pi, \"max\": np.pi, \"default\": 0.0},\n",
" \"x2\": {\"type\": \"float\", \"min\": -np.pi, \"max\": np.pi, \"default\": 0.0},\n",
" \"x3\": {\"type\": \"float\", \"min\": -np.pi, \"max\": np.pi, \"default\": 0.0},\n",
" \"a\": {\"type\": \"float\", \"min\": Ishigami_a, \"max\": Ishigami_a, \"default\": Ishigami_a},\n",
" \"b\": {\"type\": \"float\", \"min\": Ishigami_b, \"max\": Ishigami_b, \"default\": Ishigami_b},\n",
" \"N\": {\"type\": \"float\", \"min\": -100.0, \"max\": 100.0, \"default\": 0.0},\n",
" \"out_file\": {\"type\": \"string\", \"default\": \"output.csv\"}\n",
" }"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"ExecuteTime": {
"end_time": "2021-12-10T08:41:19.744486Z",
"start_time": "2021-12-10T08:41:19.741899Z"
},
"code_folding": [],
"execution": {
"iopub.execute_input": "2025-07-18T12:08:19.154580Z",
"iopub.status.busy": "2025-07-18T12:08:19.154506Z",
"iopub.status.idle": "2025-07-18T12:08:19.156381Z",
"shell.execute_reply": "2025-07-18T12:08:19.156164Z",
"shell.execute_reply.started": "2025-07-18T12:08:19.154572Z"
}
},
"outputs": [],
"source": [
"# Define varying space\n",
"def define_vary():\n",
" return {\n",
" \"x1\": cp.Uniform(-np.pi, np.pi),\n",
" \"x2\": cp.Uniform(-np.pi, np.pi),\n",
" \"x3\": cp.Uniform(-np.pi, np.pi),\n",
" \"N\": cp.Normal(0, 10.0)\n",
" }"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"ExecuteTime": {
"end_time": "2021-12-10T08:41:19.752409Z",
"start_time": "2021-12-10T08:41:19.745560Z"
},
"code_folding": [
0
],
"execution": {
"iopub.execute_input": "2025-07-18T12:08:19.156718Z",
"iopub.status.busy": "2025-07-18T12:08:19.156649Z",
"iopub.status.idle": "2025-07-18T12:08:19.160808Z",
"shell.execute_reply": "2025-07-18T12:08:19.160591Z",
"shell.execute_reply.started": "2025-07-18T12:08:19.156711Z"
}
},
"outputs": [],
"source": [
"# Set up and run a campaign\n",
"def run_campaign(pce_order=2, use_files=False):\n",
"\n",
" times = np.zeros(7)\n",
"\n",
" time_start = time.time()\n",
" time_start_whole = time_start\n",
"\n",
" # Set up a fresh campaign called \"Ishigami_pce.\"\n",
" my_campaign = uq.Campaign(name='Ishigami_pce.')\n",
"\n",
" # Create an encoder and decoder for PCE test app\n",
" if use_files:\n",
" encoder = uq.encoders.GenericEncoder(template_fname='Ishigami.template',\n",
" delimiter='$',\n",
" target_filename='Ishigami_in.json')\n",
"\n",
" decoder = uq.decoders.SimpleCSV(target_filename=\"output.csv\",\n",
" output_columns=[\"Ishigami\"])\n",
"\n",
" execute = uq.actions.ExecuteLocal('python3 %s/Ishigami.py Ishigami_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_ishigami_model))\n",
"\n",
" # Add the app (automatically set as current app)\n",
" my_campaign.add_app(name=\"Ishigami\", params=define_params(), actions=actions)\n",
"\n",
" # Create the sampler\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",
" Sampler_PCE = uq.sampling.PCESampler(vary=define_vary(), polynomial_order=pce_order)\n",
" my_campaign.set_sampler(Sampler_PCE)\n",
"\n",
" # Will draw all (of the finite set of samples)\n",
" my_campaign.draw_samples()\n",
" print('PCE order = %s' % pce_order)\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",
" # Run the cases\n",
" my_campaign.execute(sequential=True).collate(progress_bar=True)\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",
" # Get 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=[\"Ishigami\"])\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('Ishigami_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": 7,
"metadata": {
"ExecuteTime": {
"end_time": "2021-12-10T08:47:01.246304Z",
"start_time": "2021-12-10T08:41:19.753371Z"
},
"code_folding": [
0
],
"execution": {
"iopub.execute_input": "2025-07-18T12:08:19.161301Z",
"iopub.status.busy": "2025-07-18T12:08:19.161217Z",
"iopub.status.idle": "2025-07-18T12:14:05.152296Z",
"shell.execute_reply": "2025-07-18T12:14:05.151146Z",
"shell.execute_reply.started": "2025-07-18T12:08:19.161293Z"
},
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Time for phase 1 = 0.034\n",
"PCE order = 1\n",
"Number of samples = 16\n",
"Time for phase 2 = 0.056\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|█████████████████████████████████████████| 16/16 [00:00<00:00, 3128.04it/s]"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Time for phase 3 = 0.030\n",
"Time for phase 4 = 0.004\n",
"Time for phase 5 = 0.041\n",
"Time for phase 6 = 0.002\n",
"Time for phase 1 = 0.007\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"\n",
"Traceback (most recent call last):\n",
" File \"/Volumes/UserData/dpc/GIT/EasyVVUQ/env_3.12/lib/python3.12/site-packages/easyvvuq/analysis/pce_analysis.py\", line 495, in analyse\n",
" dY_hat = build_surrogate_der(fit, verbose=False)\n",
" ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n",
" File \"/Volumes/UserData/dpc/GIT/EasyVVUQ/env_3.12/lib/python3.12/site-packages/easyvvuq/analysis/pce_analysis.py\", line 335, in build_surrogate_der\n",
" assert(n1 == n2)\n",
" ^^^^^^^^\n",
"AssertionError\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"PCE order = 2\n",
"Number of samples = 81\n",
"Time for phase 2 = 0.088\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|█████████████████████████████████████████| 81/81 [00:00<00:00, 6073.81it/s]\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Time for phase 3 = 0.019\n",
"Time for phase 4 = 0.002\n",
"Time for phase 5 = 0.133\n",
"Time for phase 6 = 0.001\n",
"Time for phase 1 = 0.007\n",
"PCE order = 3\n",
"Number of samples = 256\n",
"Time for phase 2 = 0.167\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|███████████████████████████████████████| 256/256 [00:00<00:00, 6319.55it/s]"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Time for phase 3 = 0.049\n",
"Time for phase 4 = 0.004\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Time for phase 5 = 0.265\n",
"Time for phase 6 = 0.001\n",
"Time for phase 1 = 0.006\n",
"PCE order = 4\n",
"Number of samples = 625\n",
"Time for phase 2 = 0.370\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|███████████████████████████████████████| 625/625 [00:00<00:00, 6280.84it/s]"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Time for phase 3 = 0.115\n",
"Time for phase 4 = 0.008\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Time for phase 5 = 0.642\n",
"Time for phase 6 = 0.011\n",
"Time for phase 1 = 0.007\n",
"PCE order = 5\n",
"Number of samples = 1296\n",
"Time for phase 2 = 0.520\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|█████████████████████████████████████| 1296/1296 [00:00<00:00, 6776.62it/s]\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Time for phase 3 = 0.220\n",
"Time for phase 4 = 0.016\n",
"Time for phase 5 = 1.636\n",
"Time for phase 6 = 0.001\n",
"Time for phase 1 = 0.006\n",
"PCE order = 6\n",
"Number of samples = 2401\n",
"Time for phase 2 = 0.989\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|█████████████████████████████████████| 2401/2401 [00:00<00:00, 6610.41it/s]\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Time for phase 3 = 0.419\n",
"Time for phase 4 = 0.026\n",
"Time for phase 5 = 5.114\n",
"Time for phase 6 = 0.002\n",
"Time for phase 1 = 0.018\n",
"PCE order = 7\n",
"Number of samples = 4096\n",
"Time for phase 2 = 1.599\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|█████████████████████████████████████| 4096/4096 [00:00<00:00, 6753.46it/s]\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Time for phase 3 = 0.753\n",
"Time for phase 4 = 0.040\n",
"Time for phase 5 = 15.794\n",
"Time for phase 6 = 0.004\n",
"Time for phase 1 = 0.014\n",
"PCE order = 8\n",
"Number of samples = 6561\n",
"Time for phase 2 = 3.132\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|█████████████████████████████████████| 6561/6561 [00:01<00:00, 5595.64it/s]\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Time for phase 3 = 1.328\n",
"Time for phase 4 = 0.134\n",
"Time for phase 5 = 34.943\n",
"Time for phase 6 = 0.006\n",
"Time for phase 1 = 0.022\n",
"PCE order = 9\n",
"Number of samples = 10000\n",
"Time for phase 2 = 3.888\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|███████████████████████████████████| 10000/10000 [00:01<00:00, 5809.02it/s]\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Time for phase 3 = 1.941\n",
"Time for phase 4 = 0.200\n",
"Time for phase 5 = 79.951\n",
"Time for phase 6 = 0.009\n",
"Time for phase 1 = 0.015\n",
"PCE order = 10\n",
"Number of samples = 14641\n",
"Time for phase 2 = 5.835\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|███████████████████████████████████| 14641/14641 [00:02<00:00, 5799.49it/s]\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Time for phase 3 = 2.988\n",
"Time for phase 4 = 0.259\n",
"Time for phase 5 = 182.045\n",
"Time for phase 6 = 0.016\n"
]
}
],
"source": [
"# Calculate the polynomial chaos expansion for a range of orders\n",
"\n",
"R = {}\n",
"for pce_order in range(1, 11):\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_campaign(pce_order=pce_order, use_files=False)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"ExecuteTime": {
"end_time": "2021-12-10T08:47:01.296719Z",
"start_time": "2021-12-10T08:47:01.247807Z"
},
"code_folding": [
0
],
"execution": {
"iopub.execute_input": "2025-07-18T12:14:05.157076Z",
"iopub.status.busy": "2025-07-18T12:14:05.156832Z",
"iopub.status.idle": "2025-07-18T12:14:05.217663Z",
"shell.execute_reply": "2025-07-18T12:14:05.216015Z",
"shell.execute_reply.started": "2025-07-18T12:14:05.157054Z"
}
},
"outputs": [],
"source": [
"# save the results\n",
"\n",
"pickle.dump(R, open('collected_results.pickle','bw'))"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"ExecuteTime": {
"end_time": "2021-12-10T08:47:01.372146Z",
"start_time": "2021-12-10T08:47:01.297577Z"
},
"code_folding": [
0
],
"execution": {
"iopub.execute_input": "2025-07-18T12:14:05.219493Z",
"iopub.status.busy": "2025-07-18T12:14:05.219364Z",
"iopub.status.idle": "2025-07-18T12:14:05.292215Z",
"shell.execute_reply": "2025-07-18T12:14:05.291906Z",
"shell.execute_reply.started": "2025-07-18T12:14:05.219485Z"
}
},
"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",
" 0.166434 | \n",
" 0.034119 | \n",
" 0.056325 | \n",
" 0.029776 | \n",
" 0.003714 | \n",
" 0.040676 | \n",
" 0.001628 | \n",
"
\n",
" \n",
" | 2 | \n",
" 0.250982 | \n",
" 0.007045 | \n",
" 0.088102 | \n",
" 0.019497 | \n",
" 0.002389 | \n",
" 0.133038 | \n",
" 0.000772 | \n",
"
\n",
" \n",
" | 3 | \n",
" 0.494265 | \n",
" 0.007188 | \n",
" 0.167473 | \n",
" 0.049170 | \n",
" 0.004417 | \n",
" 0.264831 | \n",
" 0.001052 | \n",
"
\n",
" \n",
" | 4 | \n",
" 1.152005 | \n",
" 0.006146 | \n",
" 0.369653 | \n",
" 0.115196 | \n",
" 0.008033 | \n",
" 0.642317 | \n",
" 0.010516 | \n",
"
\n",
" \n",
" | 5 | \n",
" 2.399874 | \n",
" 0.006595 | \n",
" 0.520325 | \n",
" 0.220133 | \n",
" 0.015788 | \n",
" 1.635806 | \n",
" 0.001044 | \n",
"
\n",
" \n",
" | 6 | \n",
" 6.555950 | \n",
" 0.006280 | \n",
" 0.988502 | \n",
" 0.419108 | \n",
" 0.025505 | \n",
" 5.114128 | \n",
" 0.002055 | \n",
"
\n",
" \n",
" | 7 | \n",
" 18.209050 | \n",
" 0.017612 | \n",
" 1.598959 | \n",
" 0.753175 | \n",
" 0.040242 | \n",
" 15.793751 | \n",
" 0.004482 | \n",
"
\n",
" \n",
" | 8 | \n",
" 39.557551 | \n",
" 0.014352 | \n",
" 3.131594 | \n",
" 1.327839 | \n",
" 0.134381 | \n",
" 34.942869 | \n",
" 0.005903 | \n",
"
\n",
" \n",
" | 9 | \n",
" 86.011995 | \n",
" 0.022007 | \n",
" 3.888345 | \n",
" 1.941350 | \n",
" 0.199872 | \n",
" 79.950801 | \n",
" 0.009002 | \n",
"
\n",
" \n",
" | 10 | \n",
" 191.162857 | \n",
" 0.014980 | \n",
" 5.835231 | \n",
" 2.987610 | \n",
" 0.258682 | \n",
" 182.044722 | \n",
" 0.016074 | \n",
"
\n",
" \n",
"
\n",
"