{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "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": [] }, "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", "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 ] }, "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 ] }, "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 ] }, "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": [] }, "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 ] }, "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 ], "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Time for phase 1 = 0.026\n", "PCE order = 1\n", "Number of samples = 16\n", "Time for phase 2 = 0.070\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "100%|████████████████████████████████████████████████| 16/16 [00:00<00:00, 2804.15it/s]" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Time for phase 3 = 0.021\n", "Time for phase 4 = 0.032\n", "Time for phase 5 = 0.054\n", "Time for phase 6 = 0.003\n", "Time for phase 1 = 0.011\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "PCE order = 2\n", "Number of samples = 81\n", "Time for phase 2 = 0.124\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "100%|████████████████████████████████████████████████| 81/81 [00:00<00:00, 3330.87it/s]" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Time for phase 3 = 0.033\n", "Time for phase 4 = 0.004\n", "Time for phase 5 = 0.133\n", "Time for phase 6 = 0.001\n", "Time for phase 1 = 0.008\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "PCE order = 3\n", "Number of samples = 256\n", "Time for phase 2 = 0.207\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "100%|██████████████████████████████████████████████| 256/256 [00:00<00:00, 4990.55it/s]" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Time for phase 3 = 0.062\n", "Time for phase 4 = 0.005\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Time for phase 5 = 0.269\n", "Time for phase 6 = 0.001\n", "Time for phase 1 = 0.008\n", "PCE order = 4\n", "Number of samples = 625\n", "Time for phase 2 = 0.421\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "100%|██████████████████████████████████████████████| 625/625 [00:00<00:00, 5099.29it/s]" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Time for phase 3 = 0.143\n", "Time for phase 4 = 0.012\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Time for phase 5 = 0.656\n", "Time for phase 6 = 0.001\n", "Time for phase 1 = 0.007\n", "PCE order = 5\n", "Number of samples = 1296\n", "Time for phase 2 = 0.774\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "100%|████████████████████████████████████████████| 1296/1296 [00:00<00:00, 2887.71it/s]\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Time for phase 3 = 0.498\n", "Time for phase 4 = 0.109\n", "Time for phase 5 = 2.258\n", "Time for phase 6 = 0.002\n", "Time for phase 1 = 0.018\n", "PCE order = 6\n", "Number of samples = 2401\n", "Time for phase 2 = 1.281\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "100%|████████████████████████████████████████████| 2401/2401 [00:00<00:00, 4844.01it/s]\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Time for phase 3 = 0.560\n", "Time for phase 4 = 0.092\n", "Time for phase 5 = 5.140\n", "Time for phase 6 = 0.002\n", "Time for phase 1 = 0.015\n", "PCE order = 7\n", "Number of samples = 4096\n", "Time for phase 2 = 1.867\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "100%|████████████████████████████████████████████| 4096/4096 [00:00<00:00, 5006.41it/s]\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Time for phase 3 = 0.943\n", "Time for phase 4 = 0.107\n", "Time for phase 5 = 13.187\n", "Time for phase 6 = 0.014\n", "Time for phase 1 = 0.024\n", "PCE order = 8\n", "Number of samples = 6561\n", "Time for phase 2 = 3.065\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "100%|████████████████████████████████████████████| 6561/6561 [00:01<00:00, 5320.48it/s]\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Time for phase 3 = 1.469\n", "Time for phase 4 = 0.089\n", "Time for phase 5 = 34.744\n", "Time for phase 6 = 0.005\n", "Time for phase 1 = 0.019\n", "PCE order = 9\n", "Number of samples = 10000\n", "Time for phase 2 = 4.579\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "100%|██████████████████████████████████████████| 10000/10000 [00:01<00:00, 5055.38it/s]\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Time for phase 3 = 2.305\n", "Time for phase 4 = 0.256\n", "Time for phase 5 = 85.647\n", "Time for phase 6 = 0.008\n", "Time for phase 1 = 0.016\n", "PCE order = 10\n", "Number of samples = 14641\n", "Time for phase 2 = 7.401\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "100%|██████████████████████████████████████████| 14641/14641 [00:02<00:00, 5216.68it/s]\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Time for phase 3 = 3.314\n", "Time for phase 4 = 0.266\n", "Time for phase 5 = 169.084\n", "Time for phase 6 = 0.009\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 ] }, "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 ] }, "outputs": [ { "data": { "text/html": [ "
\n", " | Total | \n", "Phase 1 | \n", "Phase 2 | \n", "Phase 3 | \n", "Phase 4 | \n", "Phase 5 | \n", "Phase 6 | \n", "
---|---|---|---|---|---|---|---|
1 | \n", "0.206571 | \n", "0.026080 | \n", "0.070131 | \n", "0.020785 | \n", "0.032049 | \n", "0.054390 | \n", "0.002557 | \n", "
2 | \n", "0.305832 | \n", "0.010971 | \n", "0.124241 | \n", "0.033282 | \n", "0.003821 | \n", "0.132683 | \n", "0.000605 | \n", "
3 | \n", "0.552754 | \n", "0.007502 | \n", "0.207028 | \n", "0.062329 | \n", "0.005470 | \n", "0.269407 | \n", "0.000714 | \n", "
4 | \n", "1.240322 | \n", "0.008090 | \n", "0.421225 | \n", "0.142600 | \n", "0.011607 | \n", "0.655824 | \n", "0.000745 | \n", "
5 | \n", "3.649108 | \n", "0.007265 | \n", "0.773599 | \n", "0.498296 | \n", "0.108962 | \n", "2.258257 | \n", "0.001887 | \n", "
6 | \n", "7.094146 | \n", "0.018417 | \n", "1.280857 | \n", "0.559989 | \n", "0.092320 | \n", "5.140118 | \n", "0.001978 | \n", "
7 | \n", "16.133547 | \n", "0.015276 | \n", "1.867243 | \n", "0.943253 | \n", "0.106865 | \n", "13.186906 | \n", "0.013542 | \n", "
8 | \n", "39.397349 | \n", "0.024417 | \n", "3.064570 | \n", "1.469353 | \n", "0.088808 | \n", "34.744386 | \n", "0.005351 | \n", "
9 | \n", "92.814330 | \n", "0.019228 | \n", "4.578781 | \n", "2.305010 | \n", "0.255850 | \n", "85.646646 | \n", "0.008496 | \n", "
10 | \n", "180.090970 | \n", "0.015762 | \n", "7.400623 | \n", "3.313908 | \n", "0.265976 | \n", "169.084253 | \n", "0.009452 | \n", "