-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathsensitivityAnalysis.py
executable file
·164 lines (144 loc) · 6.37 KB
/
sensitivityAnalysis.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
############################################################################
#
# Imperial College London, United Kingdom
# Multifunctional Nanomaterials Laboratory
#
# Project: ERASE
# Year: 2020
# Python: Python 3.7
# Authors: Ashwin Kumar Rajagopalan (AK)
#
# Purpose:
# Script to perform a sensitivity analysis on the sensor response and
# concentration estimate
#
# Last modified:]
# - 2020-11-26, AK: Parallel processing fix
# - 2020-11-24, AK: More fix for 3 gas mole fraction
# - 2020-11-23, AK: Fix for 3 gas mole fraction
# - 2020-11-19, AK: Modify for three gas system
# - 2020-11-12, AK: Save arrayConcentration in output
# - 2020-11-12, AK: Bug fix for multipler error
# - 2020-11-11, AK: Add multipler nosie
# - 2020-11-10, AK: Improvements to run in HPC
# - 2020-11-10, AK: Add measurement nosie
# - 2020-11-06, AK: Initial creation
#
# Input arguments:
#
#
# Output arguments:
#
#
############################################################################
import numpy as np
from numpy import savez
import multiprocessing # For parallel processing
from joblib import Parallel, delayed # For parallel processing
import auxiliaryFunctions
import os
import sys
from tqdm import tqdm # To track progress of the loop
from estimateConcentration import estimateConcentration
import argparse
# For atgument parser if run through terminal. Sensor configuration provided
# as an input using --s and sorbent ids separated by a space
# e.g. python sensitivityAnalysis.py --s 6 2
parser = argparse.ArgumentParser()
parser.add_argument('--s', nargs='+', type=int)
# Get the commit ID of the current repository
gitCommitID = auxiliaryFunctions.getCommitID()
# Get the current date and time for saving purposes
simulationDT = auxiliaryFunctions.getCurrentDateTime()
# Find out the total number of cores available for parallel processing
num_cores = multiprocessing.cpu_count()
# Number of adsorbents
numberOfAdsorbents = 30
# Number of gases
numberOfGases = 2
# Sensor combination
# Check if argument provided (from terminal)
if len(sys.argv)>1:
print("Sensor configuration provided!")
for _, value in parser.parse_args()._get_kwargs():
sensorID = value
# Use default values
else:
print("\nSensor configuration not not provided. Default used!")
sensorID = [6, 2,]
# Measurement noise (Guassian noise)
meanError = 0. # [g/kg]
stdError = 0.1 # [g/kg]
# Multipler error for the sensor measurement
multiplierError = [1., 1.,]
# Custom input mole fraction for gas 1
meanMoleFracG1 = np.array([0.001, 0.01, 0.1, 0.25, 0.50, 0.75, 0.90, 0.99])
diffMoleFracG1 = 0.00 # This plus/minus the mean is the bound for uniform dist.
numberOfMoleFrac = len(meanMoleFracG1)
# For three gases generate the input concentration from a drichlet distribution
if numberOfGases == 3:
inputMoleFracALL = np.array([[0.00, 0.20, 0.80],
[0.15, 0.20, 0.65],
[0.30, 0.20, 0.50],
[0.45, 0.20, 0.35],
[0.60, 0.20, 0.20],
[0.80, 0.20, 0.00]])
numberOfMoleFrac = inputMoleFracALL.shape[0]
# Number of iterations for the estimator
numberOfIterations = 1000
# Initialize mean and standard deviation of concentration estimates
meanConcEstimate = np.zeros([numberOfMoleFrac,numberOfGases])
stdConcEstimate = np.zeros([numberOfMoleFrac,numberOfGases])
# Initialize the arrayConcentration matrix
arrayConcentration = np.zeros([numberOfMoleFrac,numberOfIterations,
numberOfGases+len(sensorID)])
# Loop through all mole fractions
for ii in range(numberOfMoleFrac):
# Generate a uniform distribution of mole fractions
if numberOfGases == 2:
inputMoleFrac = np.zeros([numberOfIterations,2])
inputMoleFrac[:,0] = np.random.uniform(meanMoleFracG1[ii]-diffMoleFracG1,
meanMoleFracG1[ii]+diffMoleFracG1,
numberOfIterations)
inputMoleFrac[:,1] = 1. - inputMoleFrac[:,0]
elif numberOfGases == 3:
inputMoleFrac = np.zeros([numberOfIterations,3])
inputMoleFrac[:,0] = inputMoleFracALL[ii,0]
inputMoleFrac[:,1] = inputMoleFracALL[ii,1]
inputMoleFrac[:,2] = inputMoleFracALL[ii,2]
# Loop over all the sorbents for a single material sensor
# Using parallel processing to loop through all the materials
arrayConcentrationTemp = Parallel(n_jobs=num_cores)(delayed(estimateConcentration)
(numberOfAdsorbents,numberOfGases,None,sensorID,
moleFraction = inputMoleFrac[ii],
multiplierError = multiplierError,
addMeasurementNoise = [meanError,stdError])
for ii in tqdm(range(inputMoleFrac.shape[0])))
# Convert the output list to a matrix
arrayConcentration[ii,:,:] = np.array(arrayConcentrationTemp)
# Check if simulationResults directory exists or not. If not, create the folder
if not os.path.exists('simulationResults'):
os.mkdir('simulationResults')
# Save the array concentration into a native numpy file
# The .npz file is saved in a folder called simulationResults (hardcoded)
filePrefix = "sensitivityAnalysis"
sensorText = str(sensorID).replace('[','').replace(']','').replace(' ','-').replace(',','')
saveFileName = filePrefix + "_" + sensorText + "_" + simulationDT + "_" + gitCommitID;
savePath = os.path.join('simulationResults',saveFileName)
# Save the results as an array
if numberOfGases == 2:
savez (savePath, numberOfGases = numberOfGases,
numberOfIterations = numberOfIterations,
trueMoleFrac = meanMoleFracG1,
multiplierError = multiplierError,
meanError = meanError,
stdError = stdError,
arrayConcentration = arrayConcentration)
if numberOfGases == 3:
savez (savePath, numberOfGases = numberOfGases,
numberOfIterations = numberOfIterations,
trueMoleFrac = inputMoleFracALL,
multiplierError = multiplierError,
meanError = meanError,
stdError = stdError,
arrayConcentration = arrayConcentration)