diff --git a/.gitignore b/.gitignore new file mode 100755 index 0000000..d1d4543 --- /dev/null +++ b/.gitignore @@ -0,0 +1,157 @@ +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[cod] +*$py.class + +# Graphic files +*.eps +*.pdf +*.png +*.tif +*.jpg + +# Logs and databases +*.log +*.sql +*.sqlite +#*.bib +*.xlsx +*.txt + +# Python related outputs +*.npy +*.npz + +# C extensions +*.so + +# Distribution / packaging +.Python +build/ +develop-eggs/ +dist/ +downloads/ +eggs/ +.eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +wheels/ +share/python-wheels/ +*.egg-info/ +.installed.cfg +*.egg +MANIFEST + +# PyInstaller +# Usually these files are written by a python script from a template +# before PyInstaller builds the exe, so as to inject date/other infos into it. +*.manifest +*.spec + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt + +# Unit test / coverage reports +htmlcov/ +.tox/ +.nox/ +.coverage +.coverage.* +.cache +nosetests.xml +coverage.xml +*.cover +*.py,cover +.hypothesis/ +.pytest_cache/ +cover/ + +# Translations +*.mo +*.pot + +# Django stuff: +*.log +local_settings.py +db.sqlite3 +db.sqlite3-journal + +# Flask stuff: +instance/ +.webassets-cache + +# Scrapy stuff: +.scrapy + +# Sphinx documentation +docs/_build/ + +# PyBuilder +.pybuilder/ +target/ + +# Jupyter Notebook +.ipynb_checkpoints + +# IPython +profile_default/ +ipython_config.py + +# pyenv +# For a library or package, you might want to ignore these files since the code is +# intended to run in multiple environments; otherwise, check them in: +# .python-version + +# pipenv +# According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control. +# However, in case of collaboration, if having platform-specific dependencies or dependencies +# having no cross-platform support, pipenv may install dependencies that don't work, or not +# install all needed dependencies. +#Pipfile.lock + +# PEP 582; used by e.g. github.com/David-OConnor/pyflow +__pypackages__/ + +# Celery stuff +celerybeat-schedule +celerybeat.pid + +# SageMath parsed files +*.sage.py + +# Environments +.env +.venv +env/ +venv/ +ENV/ +env.bak/ +venv.bak/ + +# Spyder project settings +.spyderproject +.spyproject + +# Rope project settings +.ropeproject + +# mkdocs documentation +/site + +# mypy +.mypy_cache/ +.dmypy.json +dmypy.json + +# Pyre type checker +.pyre/ + +# pytype static type analyzer +.pytype/ + +# Cython debug symbols +cython_debug/ diff --git a/auxiliaryFunctions/__init__.py b/auxiliaryFunctions/__init__.py new file mode 100755 index 0000000..4717a6e --- /dev/null +++ b/auxiliaryFunctions/__init__.py @@ -0,0 +1 @@ +# __init__.py from .getCommitID import getCommitID from .getCurrentDateTime import getCurrentDateTime \ No newline at end of file diff --git a/auxiliaryFunctions/addInertGasToMaterials.py b/auxiliaryFunctions/addInertGasToMaterials.py new file mode 100644 index 0000000..0480bd4 --- /dev/null +++ b/auxiliaryFunctions/addInertGasToMaterials.py @@ -0,0 +1,77 @@ +############################################################################ +# +# Imperial College London, United Kingdom +# Multifunctional Nanomaterials Laboratory +# +# Project: ERASE +# Year: 2020 +# Python: Python 3.7 +# Authors: Ashwin Kumar Rajagopalan (AK) +# +# Purpose: +# Adds an inert gas to the exisitng material property matrix for an n gas +# system. This is done to simulate the sensor array with a full model. The +# inert is used to clean the sensor array +# +# Last modified: +# - 2021-01-20, AK: Initial creation +# +# Input arguments: +# +# +# Output arguments: +# +# +############################################################################ + +def addInertGasToMaterials(numberOfGases): + import numpy as np + from numpy import load + from numpy import savez + import os + import auxiliaryFunctions + + # Get the commit ID of the current repository + gitCommitID = auxiliaryFunctions.getCommitID() + + # Get the current date and time for saving purposes + simulationDT = auxiliaryFunctions.getCurrentDateTime() + + # For now load a given adsorbent isotherm material file + if numberOfGases == 2: + loadFileName = "isothermParameters_20201020_1756_5f263af.npz" # Two gases + elif numberOfGases == 3: + loadFileName = "isothermParameters_20201022_1056_782efa3.npz" # Three gases + hypoAdsorbentFile = os.path.join('../inputResources',loadFileName); + + # Check if the file with the adsorbent properties exist + if os.path.exists(hypoAdsorbentFile): + loadedFileContent = load(hypoAdsorbentFile) + adsorbentIsothermTemp = loadedFileContent['adsIsotherm'] + adsorbentDensity = loadedFileContent['adsDensity'] + molecularWeightTemp = loadedFileContent['molWeight'] + else: + errorString = "Adsorbent property file " + hypoAdsorbentFile + " does not exist." + raise Exception(errorString) + + # Create adsorent isotherm matrix with the addition of an intert gas + adsorbentIsotherm = np.zeros([numberOfGases+1,3,adsorbentIsothermTemp.shape[2]]) + adsorbentIsotherm[0:numberOfGases,:,:] = adsorbentIsothermTemp + + # Add the moleuclar weight of the inert (assumed to be helium) + molecularWeight = np.concatenate((molecularWeightTemp,np.array([4.00]))) + + # Save the adsorbent isotherm parameters into a native numpy file + # The .npz file is saved in a folder called inputResources (hardcoded) + filePrefix = "isothermParameters" + saveFileName = filePrefix + "_" + simulationDT + "_" + gitCommitID + ".npz"; + savePath = os.path.join('../inputResources',saveFileName) + + # Check if inputResources directory exists or not. If not, create the folder + if not os.path.exists('../inputResources'): + os.mkdir('../inputResources') + + # Save the adsorbent material array + savez (savePath, adsIsotherm = adsorbentIsotherm, + adsDensity = adsorbentDensity, + molWeight = molecularWeight) \ No newline at end of file diff --git a/auxiliaryFunctions/getCommitID.py b/auxiliaryFunctions/getCommitID.py new file mode 100755 index 0000000..e92bb21 --- /dev/null +++ b/auxiliaryFunctions/getCommitID.py @@ -0,0 +1,36 @@ +############################################################################ +# +# Imperial College London, United Kingdom +# Multifunctional Nanomaterials Laboratory +# +# Project: ERASE +# Year: 2020 +# Python: Python 3.7 +# Authors: Ashwin Kumar Rajagopalan (AK) +# +# Purpose: +# Generates a short SHA key of the git commit id in the current branch and +# repository +# +# Last modified: +# - 2020-10-19, AK: Initial creation +# +# Input arguments: +# - N/A +# +# Output arguments: +# - short_sha: Short git commit ID +# +############################################################################ + +def getCommitID(): + # Use gitpython to get the git information of the current repository + import git + repo = git.Repo(search_parent_directories=True) + # Get the simple hashing algorithm tag (SHA) + sha = repo.head.commit.hexsha + # Parse the first six characters of the sha + short_sha = repo.git.rev_parse(sha, short=7) + + # Return the git commit id + return short_sha \ No newline at end of file diff --git a/auxiliaryFunctions/getCurrentDateTime.py b/auxiliaryFunctions/getCurrentDateTime.py new file mode 100755 index 0000000..8ca8ec7 --- /dev/null +++ b/auxiliaryFunctions/getCurrentDateTime.py @@ -0,0 +1,33 @@ +############################################################################ +# +# Imperial College London, United Kingdom +# Multifunctional Nanomaterials Laboratory +# +# Project: ERASE +# Year: 2020 +# Python: Python 3.7 +# Authors: Ashwin Kumar Rajagopalan (AK) +# +# Purpose: +# Obtain the current date and time to be used either for saving in file +# name or to enhance traceability of the simulation +# +# Last modified: +# - 2020-10-19, AK: Initial creation +# +# Input arguments: +# - N/A +# +# Output arguments: +# - simulationDT: Current date and time in YYYYmmdd_HHMM format +# +############################################################################ + +def getCurrentDateTime(): + # Get the current date and time for saving purposes + from datetime import datetime + now = datetime.now() + simulationDT = now.strftime("%Y%m%d_%H%M") + + # Return the current date and time + return simulationDT \ No newline at end of file diff --git a/estimateConcentration.py b/estimateConcentration.py new file mode 100755 index 0000000..412b04e --- /dev/null +++ b/estimateConcentration.py @@ -0,0 +1,169 @@ +############################################################################ +# +# Imperial College London, United Kingdom +# Multifunctional Nanomaterials Laboratory +# +# Project: ERASE +# Year: 2020 +# Python: Python 3.7 +# Authors: Ashwin Kumar Rajagopalan (AK) +# +# Purpose: +# Estimates the concentration of the gas mixture using the sensor response +# by minimization of the sum of square error between the "true" and the +# "estimated" differences in the change of mass of the sensor array +# +# Last modified: +# - 2021-01-25, AK: Add zero response to avoid division by zero (IMP!) +# - 2021-01-21, AK: Add full model functionality +# - 2020-11-12, AK: Bug fix for multipler error +# - 2020-11-11, AK: Add multiplier error to true sensor response +# - 2020-11-10, AK: Add measurement noise to true sensor response +# - 2020-11-09, AK: Changes to initial condition and optimizer bounds +# - 2020-11-05, AK: Introduce keyword argument for custom mole fraction +# - 2020-10-30, AK: Fix to find number of gases +# - 2020-10-22, AK: Change error to relative from absolute, add opt bounds, +# input arguments, and initial guess +# - 2020-10-21, AK: Initial creation +# +# Input arguments: +# +# +# Output arguments: +# +# +############################################################################ + +def estimateConcentration(numberOfAdsorbents, numberOfGases, moleFracID, sensorID, **kwargs): + import numpy as np + from generateTrueSensorResponse import generateTrueSensorResponse + from scipy.optimize import basinhopping + + # Total pressure of the gas [Pa] + if 'pressureTotal' in kwargs: + pressureTotal = np.array(kwargs["pressureTotal"]); + else: + pressureTotal = np.array([1.e5]); + + # Temperature of the gas [K] + # Can be a vector of temperatures + if 'temperature' in kwargs: + temperature = np.array(kwargs["temperature"]); + else: + temperature = np.array([298.15]); + + # Sensor combinations used in the array. This is a [gx1] vector that maps to + # the sorbent/material ID generated using the + # generateHypotheticalAdsorbents.py function + sensorID = np.array(sensorID) + + # Get the individual sensor reponse for all the given "experimental/test" concentrations + if 'fullModel' in kwargs: + if kwargs["fullModel"]: + fullModelFlag = True + else: + fullModelFlag = False + else: + fullModelFlag = False + + # Add measurement noise for the true measurement if the user wants it + measurementNoise = np.zeros(sensorID.shape[0]) + if 'addMeasurementNoise' in kwargs: + # The mean and the standard deviation of the Gaussian error is an + # input from the user + measurementNoise = np.random.normal(kwargs["addMeasurementNoise"][0], + kwargs["addMeasurementNoise"][1], + sensorID.shape[0]) + + # Check if it is for full model or not + # Full model condition + if fullModelFlag: + # Parse out the true sensor response from the input (time resolved) + # and add measurement noise if asked for. There is no multipler error for + # full model simualtions + # Note that using the full model here is only for comparison purposes + # When kinetics are available the other estimateConcentration function + # should be used + if 'fullModelResponse' in kwargs: + fullModelResponse = kwargs["fullModelResponse"] + multiplierError = np.ones(sensorID.shape[0]) # Always set to 1. + arrayTrueResponse = np.zeros(sensorID.shape[0]) + for ii in range(sensorID.shape[0]): + arrayTrueResponse[ii] = (multiplierError[ii]*fullModelResponse[ii] + + measurementNoise[ii]) + else: + errorString = "Sensor response from full model not available. You should not be here!" + raise Exception(errorString) + # Equilibrium condition + else: + # Get the individual sensor reponse for all the given "experimental/test" concentrations + if 'moleFraction' in kwargs: + sensorTrueResponse = generateTrueSensorResponse(numberOfAdsorbents,numberOfGases, + pressureTotal,temperature, moleFraction = kwargs["moleFraction"]) + moleFracID = 0 # Index becomes a scalar quantity + else: + sensorTrueResponse = generateTrueSensorResponse(numberOfAdsorbents,numberOfGases, + pressureTotal,temperature) + # True mole fraction index (provide the index corresponding to the true + # experimental mole fraction (0-4, check generateTrueSensorResponse.py) + moleFracID = moleFracID + + # Add a multiplier error for the true measurement if the user wants it + multiplierError = np.ones(sensorID.shape[0]) + if 'multiplierError' in kwargs: + # The mean and the standard deviation of the Gaussian error is an + # input from the user + multiplierErrorTemp = kwargs["multiplierError"] + multiplierError[0:len(multiplierErrorTemp)] = multiplierErrorTemp + # Parse out the true sensor response for a sensor array with n number of + # sensors given by sensorID + arrayTrueResponse = np.zeros(sensorID.shape[0]) + for ii in range(sensorID.shape[0]): + arrayTrueResponse[ii] = (multiplierError[ii]*sensorTrueResponse[sensorID[ii],moleFracID] + + measurementNoise[ii]) + + # Replace all negative values to eps (for physical consistency). Set to + # eps to avoid division by zero + # Print if any of the responses are negative + if any(ii<=0. for ii in arrayTrueResponse): + print("Number of negative response: " + str(sum(arrayTrueResponse<=0))) + arrayTrueResponse[arrayTrueResponse<=0.] = np.finfo(float).eps + + # Pack the input parameters/arguments useful to compute the objective + # function to estimate the mole fraction as a tuple + inputParameters = (arrayTrueResponse, pressureTotal, temperature, + sensorID, multiplierError) + + # Minimize an objective function to compute the mole fraction of the feed + # gas to the sensor + initialCondition = np.random.uniform(0,1,numberOfGases) # Initial guess + optBounds = np.tile([0.,1.], (numberOfGases,1)) # BOunding the mole fractions + optCons = {'type':'eq','fun': lambda x: sum(x) - 1} # Cinstrain the sum to 1. + # Use the basin hopping minimizer to escape local minima when evaluating + # the function. The number of iterations is hard-coded and fixed at 50 + estMoleFraction = basinhopping(concObjectiveFunction, initialCondition, + minimizer_kwargs = {"args": inputParameters, + "bounds": optBounds, + "constraints": optCons}, + niter = 50) + return np.concatenate((sensorID,estMoleFraction.x), axis=0) + +# func: concObjectiveFunction, computes the sum of square error for the +# finger print for varying gas concentrations, using the minimize function +def concObjectiveFunction(x, *inputParameters): + import numpy as np + from simulateSensorArray import simulateSensorArray + + # Unpack the tuple that contains the true response, pressure, temperature, + # and the sensor identifiers + arrayTrueResponse, pressureTotal, temperature, sensorID, multiplierError = inputParameters + + # Reshape the mole fraction to a row vector for compatibility + moleFraction = np.array([x]) # This is needed to keep the structure as a row instead of column + # moleFraction = np.array([x,1.-x]).T # This is needed to keep the structure as a row instead of column + + # Compute the sensor reponse for a given mole fraction input + arraySimResponse = simulateSensorArray(sensorID, pressureTotal, temperature, moleFraction) * multiplierError + + # Compute the sum of the error for the sensor array + return sum(np.power((arrayTrueResponse - arraySimResponse)/arrayTrueResponse,2)) \ No newline at end of file diff --git a/estimateConcentrationFullModel.py b/estimateConcentrationFullModel.py new file mode 100644 index 0000000..b36f44a --- /dev/null +++ b/estimateConcentrationFullModel.py @@ -0,0 +1,147 @@ +############################################################################ +# +# Imperial College London, United Kingdom +# Multifunctional Nanomaterials Laboratory +# +# Project: ERASE +# Year: 2021 +# Python: Python 3.7 +# Authors: Ashwin Kumar Rajagopalan (AK) +# +# Purpose: +# Estimates the concentration of the gas mixture using the sensor response +# by minimization of the sum of square error between the "true" and the +# "estimated" differences in the change of mass of the sensor array as a +# function of time. Here the full model instead of equilibrium characteristics +# is used +# +# Last modified: +# - 2021-01-25, AK: Modify the objective function from eqbm method +# - 2021-01-22, AK: Initial creation +# +# Input arguments: +# +# +# Output arguments: +# +# +############################################################################ + +def estimateConcentrationFullModel(numberOfGases, sensorID, **kwargs): + import numpy as np + from generateTrueSensorResponse import generateTrueSensorResponse + from scipy.optimize import basinhopping + + # Total pressure of the gas [Pa] + if 'pressureTotal' in kwargs: + pressureTotal = np.array(kwargs["pressureTotal"]); + else: + pressureTotal = np.array([1.e5]); + + # Temperature of the gas [K] + # Can be a vector of temperatures + if 'temperature' in kwargs: + temperature = np.array(kwargs["temperature"]); + else: + temperature = np.array([298.15]); + + # Sensor combinations used in the array. This is a [gx1] vector that maps to + # the sorbent/material ID generated using the + # generateHypotheticalAdsorbents.py function + sensorID = np.array(sensorID) + + # Parse out the true sensor response from the input (time resolved) + if 'fullModelResponse' in kwargs: + arrayTrueResponse = kwargs["fullModelResponse"] + else: + errorString = "Sensor response from full model not available. You should not be here!" + raise Exception(errorString) + + # Parse out the rate constant for the materials from the input + if 'rateConstant' in kwargs: + rateConstant = kwargs["rateConstant"] + else: + errorString = "Rate constant for materials not available. You should not be here!" + raise Exception(errorString) + + # Parse out the time of integration from the input + if 'timeInt' in kwargs: + timeInt = kwargs["timeInt"] + else: + errorString = "Integration time for the model not available. You should not be here!" + raise Exception(errorString) + + # Parse out the feed flow rate from the input + if 'flowIn' in kwargs: + flowIn = kwargs["flowIn"] + else: + errorString = "Feed flow rate not available. You should not be here!" + raise Exception(errorString) + + # Replace all negative values to eps (for physical consistency). Set to + # eps to avoid division by zero + # Print if any of the responses are negative + if (arrayTrueResponse<=0.).any(): + print("Number of zero/negative response: " + str(np.sum(arrayTrueResponse<=0.))) + arrayTrueResponse[arrayTrueResponse<=0.] = np.finfo(float).eps + + # Pack the input parameters/arguments useful to compute the objective + # function to estimate the mole fraction as a tuple + inputParameters = (arrayTrueResponse, flowIn, sensorID, rateConstant, timeInt) + + # Minimize an objective function to compute the mole fraction of the feed + # gas to the sensor + initialCondition = np.random.uniform(0,1,numberOfGases+1) # Initial guess + optBounds = np.tile([0.,1.], (numberOfGases+1,1)) # BOunding the mole fractions + optCons = {'type':'eq','fun': lambda x: sum(x) - 1} # Cinstrain the sum to 1. + # Use the basin hopping minimizer to escape local minima when evaluating + # the function. The number of iterations is hard-coded and fixed at 50 + estMoleFraction = basinhopping(concObjectiveFunction, initialCondition, + minimizer_kwargs = {"args": inputParameters, + "bounds": optBounds, + "constraints": optCons}, + niter = 50) + return np.concatenate((sensorID,estMoleFraction.x), axis=0) + +# func: concObjectiveFunction, computes the sum of square error for the +# finger print for varying gas concentrations, using the minimize function +def concObjectiveFunction(x, *inputParameters): + import numpy as np + from simulateSensorArray import simulateSensorArray + from simulateFullModel import simulateFullModel + + # Unpack the tuple that contains the true response, sensor identifiers, + # rate constant and time of integration + # IMPORTANT: Pressure and Temperature not considered here!!! + arrayTrueResponse, flowIn, sensorID, rateConstant, timeInt = inputParameters + + # Reshape the mole fraction to a row vector for compatibility + moleFraction = np.array(x) # This is needed to keep the structure as a row instead of column + + # Compute the sensor reponse for a given mole fraction input + # Calls the full sensor model for each material + outputStruct = {} + for ii in range(len(sensorID)): + timeSim , _ , sensorFingerPrint, _ = simulateFullModel(sensorID = sensorID[ii], + rateConstant = rateConstant[ii], + feedMoleFrac = moleFraction, + timeInt = timeInt, + flowIn = flowIn) + outputStruct[ii] = {'timeSim':timeSim, + 'sensorFingerPrint':sensorFingerPrint} # Check simulateFullModel.py for entries + + # Prepare the simulated response at a given mole fraction to compute the + # objective function + timeSim = [] + timeSim = outputStruct[0]["timeSim"] # Need for array initialization + arraySimResponse = np.zeros([len(timeSim),len(sensorID)]) + for ii in range(len(sensorID)): + arraySimResponse[:,ii] = outputStruct[ii]["sensorFingerPrint"] + + # Prepare the objective function for the optimizer + relError = np.divide((arrayTrueResponse - arraySimResponse), arrayTrueResponse) + relErrorPow2 = np.power(relError, 2) # Square of relative error + objFunction = np.sum(relErrorPow2) # Sum over all the sensors and times + + # Return the sum of squares of relative errors + return objFunction \ No newline at end of file diff --git a/experimental/simulateVolumetricSystem.py b/experimental/simulateVolumetricSystem.py new file mode 100644 index 0000000..5161941 --- /dev/null +++ b/experimental/simulateVolumetricSystem.py @@ -0,0 +1,267 @@ +############################################################################ +# +# Imperial College London, United Kingdom +# Multifunctional Nanomaterials Laboratory +# +# Project: ERASE +# Year: 2021 +# Python: Python 3.7 +# Authors: Ashwin Kumar Rajagopalan (AK) +# +# Purpose: +# Simulates a volumetric system that can be used to estimate pure component +# isotherm and kinetics +# +# Last modified: +# - 2021-02-19, AK: Fix for ode solver and add mass balance (pressure) +# - 2021-02-18, AK: Initial creation +# +# Input arguments: +# +# +# Output arguments: +# +# +############################################################################ + +def simulateVolumetricSystem(**kwargs): + import numpy as np + from scipy.integrate import solve_ivp + import auxiliaryFunctions + from simulateSensorArray import simulateSensorArray + + # Plot flag + plotFlag = False + + # Get the commit ID of the current repository + gitCommitID = auxiliaryFunctions.getCommitID() + + # Get the current date and time for saving purposes + currentDT = auxiliaryFunctions.getCurrentDateTime() + + # Sensor ID + if 'sensorID' in kwargs: + sensorID = np.array([kwargs["sensorID"]]) + else: + sensorID = np.array([6]) + + # Time span for integration [tuple with t0 and tf] + if 'timeInt' in kwargs: + timeInt = kwargs["timeInt"] + else: + timeInt = (0.0,500) + + # Kinetic rate constant [/s] - Only one constant (pure gas) + if 'rateConstant' in kwargs: + rateConstant = np.array(kwargs["rateConstant"]) + else: + rateConstant = np.array([0.1]) + + # Valve rate constant [mol/s/Pa] - Only one constant (pure gas) + if 'valveConstant' in kwargs: + valveConstant = np.array(kwargs["valveConstant"]) + else: + valveConstant = np.array([0.04e5]) # From 10.1021/la026451+ (different units) + + # Temperature of the gas [K] + # Can be a vector of temperatures + if 'temperature' in kwargs: + temperature = np.array(kwargs["temperature"]); + else: + temperature = np.array([298.15]); + + # Dosing volume pressure [Pa] + if 'pressDose' in kwargs: + pressDose = np.array(kwargs["pressDose"]); + else: + pressDose = np.array([1e5]); + + # Uptake volume pressure [Pa] + if 'pressUptake' in kwargs: + pressUptake = np.array(kwargs["pressUptake"]); + else: + pressUptake = np.array([0.e5]); + + # Dosing volume of the setup [m3] + if 'volDose' in kwargs: + volDose = np.array(kwargs["volDose"]); + else: + volDose = np.array([385e-6]); # From 10.1021/la026451+ (different units) + + # Uptake volume of the setup [m3] + if 'volUptake' in kwargs: + volUptake = np.array(kwargs["volUptake"]); + else: + volUptake = np.array([366e-6]); # From 10.1021/la026451+ (different units) + + # Uptake voidFrac of the uptake volume [-] + if 'voidFrac' in kwargs: + voidFrac = np.array(kwargs["voidFrac"]); + else: + voidFrac = np.array([0.98]); # From 10.1021/la026451+ (calculated) + + # Initial Gas Mole Fraction [-] + # Only pure gas system for volumetric system (DO NOT CHANGE) + initMoleFrac = np.array([1.,0.]) + + # Number of gases - Note this is only for loading purposes + # Only pure gas is simulated + # Used to chose the isotherm properties from the inputResources folder + numberOfGases = len(initMoleFrac) + + # Prepare tuple of input parameters for the ode solver + inputParameters = (sensorID, rateConstant, valveConstant, + numberOfGases, initMoleFrac, temperature, voidFrac, + volDose, volUptake) + + # Obtain the initial solid loading in the material for initialization + materialLoadingPerGasVol = simulateSensorArray(sensorID, pressUptake, + temperature, np.array([initMoleFrac])) + + # Initial conditions for the solver + initialConditions = np.zeros([4]) + initialConditions[0] = pressDose # Dosing pressure [Pa] + initialConditions[1] = pressUptake # Uptake pressure [Pa] + initialConditions[2] = materialLoadingPerGasVol[0] # Initial solid loading [mol/m3] + # 2 is the moles through the valve + + # Solve the system of ordinary differential equations + # Stiff solver used for the problem: BDF or Radau + outputSol = solve_ivp(solveVolumetricSystemEquation, timeInt, initialConditions, + method='Radau', t_eval = np.arange(timeInt[0],timeInt[1],0.01), + rtol = 1e-6, args = inputParameters) + + # Parse out the time and the pressure and loading from the output + timeSim = outputSol.t + resultMat = outputSol.y + + # Compute the mass adsrobed based on mass balance (with pressures) + massAdsorbed = estimateMassAdsorbed(resultMat[0,:], resultMat[1,:], + materialLoadingPerGasVol[0]*(1-voidFrac)*volUptake, + volDose, volUptake, voidFrac, + temperature) + + # Call the plotting function + if plotFlag: + plotInputs = (timeSim, resultMat, massAdsorbed, voidFrac, volUptake, + sensorID, gitCommitID, currentDT) + plotFullModelResult(plotInputs) + + # Return time, the output matrix, and the parameters used in the simulation + return timeSim, resultMat, massAdsorbed, inputParameters + +# func: solveVolumetricSystemEquation +# Solves the system of ODEs to simulate the volumetric system +def solveVolumetricSystemEquation(t, f, *inputParameters): + import numpy as np + from simulateSensorArray import simulateSensorArray + + # Gas constant + Rg = 8.314; # [J/mol K] + + # Unpack the tuple of input parameters used to solve equations + sensorID, rateConstant, valveConstant, numberOfGases, initMoleFrac, temperature, voidFrac, volDose, volUptake = inputParameters + + # Initialize the derivatives to zero + df = np.zeros([4]) + + materialLoadingPerGasVol = simulateSensorArray(sensorID, f[1], + temperature, np.array([initMoleFrac])) + + # Linear driving force model (derivative of solid phase loadings) + df[2] = rateConstant*(materialLoadingPerGasVol[0]-f[2]) + + # Valve flow model (derivative of gas flow through the valve) + df[3] = valveConstant*(f[0]-f[1]) + + # Rate of change of pressure in the dosing volume + df[0] = -(Rg*temperature/volDose)*df[3] + + # Rate of change of pressure in the uptake volume + df[1] = (Rg*temperature/(volUptake*voidFrac))*(df[3] - (1-voidFrac)*volUptake*df[2]) + + # Return the derivatives for the solver + return df + +# func: estimateMassAdsorbed +# Given the pressure from the two volumes, the volumes and temperature +# compute the mass adsorbed +def estimateMassAdsorbed(pressDoseALL, pressUptakeALL, initMass, + volDose, volUptake, voidFrac, temperature): + import numpy as np + + # Gas constant + Rg = 8.314; # [J/mol K] + + # Compute the pressure difference over time on dosing side + delPressDosing = pressDoseALL[0] - pressDoseALL # Pd(0) - Pd(t) + + # Compute the pressure difference over time on uptake side + delPressUptake = pressUptakeALL - pressUptakeALL[0] # Pu(t) - Pu(0) + + # Calculate mass adsorbed + # massAdsorbed = Initial mass adsorbed + uptake + massAdsorbed = initMass + (np.multiply(delPressDosing,(volDose/(Rg*temperature))) + - np.multiply(delPressUptake,((voidFrac*volUptake)/(Rg*temperature)))) + + # Return the mass adsorbed + return massAdsorbed + +# func: plotFullModelResult +# Plots the model output for the conditions simulated locally +def plotFullModelResult(plotInputs): + import numpy as np + import os + import matplotlib.pyplot as plt + + # Unpack the inputs + timeSim, resultMat, massAdsorbed, voidFrac, volUptake, sensorID, gitCommitID, currentDT = plotInputs + + # Save settings + saveFlag = False + saveFileExtension = ".png" + + os.chdir("plotFunctions") + # Plot the solid phase compositions + plt.style.use('doubleColumn.mplstyle') # Custom matplotlib style file + fig = plt.figure + ax = plt.subplot(1,2,1) + ax.plot(timeSim, resultMat[0,:]/resultMat[0,0], + linewidth=1.5,color='r', label = 'Dosing') + ax.plot(timeSim, resultMat[1,:]/resultMat[0,0], + linewidth=1.5,color='b', label = 'Uptake') + ax.set(xlabel='$t$ [s]', + ylabel='$P/P_0$ [-]', + xlim = [timeSim[0], timeSim[-1]], ylim = [0., 1.]) + ax.locator_params(axis="x", nbins=4) + ax.locator_params(axis="y", nbins=4) + ax.legend() + + ax = plt.subplot(1,2,2) + volAds = (1-voidFrac)*volUptake # Volume of adsorbent [m3] + ax.plot(timeSim, resultMat[2,:]*volAds, + linewidth=1.5,color='k', + label = 'From model') # Uptake estimated from the model + ax.plot(timeSim, massAdsorbed,'--', + linewidth=1.5,color='r', + label = 'From pressure') # Uptake estimated from mass balance with pressures + ax.set(xlabel='$t$ [s]', + ylabel='$q$ [mol]', + xlim = [timeSim[0], timeSim[-1]], + ylim = [0.9*np.min(resultMat[2,:])*volAds, + 1.1*np.max(resultMat[2,:])*volAds]) + ax.locator_params(axis="x", nbins=4) + ax.locator_params(axis="y", nbins=4) + ax.legend() + + # Save the figure + if saveFlag: + # FileName: solidLoadingFM___ + saveFileName = "volumetricSysFM_" + str(sensorID) + "_" + currentDT + "_" + gitCommitID + saveFileExtension + savePath = os.path.join('..','simulationFigures',saveFileName.replace('[','').replace(']','')) + # Check if inputResources directory exists or not. If not, create the folder + if not os.path.exists(os.path.join('..','simulationFigures')): + os.mkdir(os.path.join('..','simulationFigures')) + plt.savefig (savePath) + plt.show() + os.chdir("..") \ No newline at end of file diff --git a/fullModelConcentrationEstimatorWrapper.py b/fullModelConcentrationEstimatorWrapper.py new file mode 100644 index 0000000..245fcb6 --- /dev/null +++ b/fullModelConcentrationEstimatorWrapper.py @@ -0,0 +1,178 @@ +############################################################################ +# +# Imperial College London, United Kingdom +# Multifunctional Nanomaterials Laboratory +# +# Project: ERASE +# Year: 2021 +# Python: Python 3.7 +# Authors: Ashwin Kumar Rajagopalan (AK) +# +# Purpose: +# Wrapper function to estimate the concentration given a time-resolved +# measurement from the full model +# +# Last modified: +# - 2021-02-03, AK: Change inputParameters (add adsorbent density) +# - 2021-01-26, AK: Add noise to true measurement +# - 2021-01-25, AK: Integrate full model concentration estimator +# - 2021-01-21, AK: Initial creation +# +# Input arguments: +# +# +# Output arguments: +# +# +############################################################################ + +import auxiliaryFunctions +from simulateFullModel import simulateFullModel +from estimateConcentration import estimateConcentration +from estimateConcentrationFullModel import estimateConcentrationFullModel +from tqdm import tqdm # To track progress of the loop +import numpy as np +from numpy import savez +from joblib import Parallel, delayed # For parallel processing +import multiprocessing +import os +import time +import socket + +# Find out the total number of cores available for parallel processing +num_cores = multiprocessing.cpu_count() + +# Get the commit ID of the current repository +gitCommitID = auxiliaryFunctions.getCommitID() + +# Get the current date and time for saving purposes +currentDT = auxiliaryFunctions.getCurrentDateTime() + +# Flag to determine whether concentration estimated accounted for kinetics +# (False) or not (True) +flagIndTime = True + +# Flag to determine whether constant pressure (False) or constant flow rate +# model to be used (True) +modelConstF = False + +# Sensor ID +sensorID = [6,2] # [-] + +# Number of Gases +numberOfGases = 2 # [-] + +# Rate Constant +rateConstant = ([.01,.01,10000.0], + [.01,.01,10000.0]) # [1/s] + +# Feed mole fraction +feedMoleFrac = [0.1,0.9,0.0] # [-] + +# Time span for integration [tuple with t0 and tf] [s] +timeInt = (0,1000) # [s] + +# Volumetric flow rate [m3/s] +flowIn = 5e-7 # [s] + +# Measurement noise characteristics +meanNoise = 0.0 # [g/kg] +stdNoise = 0.0 # [g/kg] + +# Loop over all rate constants +outputStruct = {} +for ii in tqdm(range(len(sensorID))): + # Call the full model with a given rate constant + timeSim, _ , sensorFingerPrint, inputParameters = simulateFullModel(sensorID = sensorID[ii], + rateConstant = rateConstant[ii], + feedMoleFrac = feedMoleFrac, + timeInt = timeInt, + flowIn = flowIn) + + # Generate the noise data for all time instants and all materials + measurementNoise = np.random.normal(meanNoise, stdNoise, len(timeSim)) + sensorFingerPrintRaw = sensorFingerPrint # Raw sensor finger print + # Add measurement noise to the sensor finger print + sensorFingerPrint = sensorFingerPrint + measurementNoise + + outputStruct[ii] = {'timeSim':timeSim, + 'sensorFingerPrint':sensorFingerPrint, + 'sensorFingerPrintRaw':sensorFingerPrintRaw, # Without noise + 'measurementNoise':measurementNoise, + 'inputParameters':inputParameters} # Check simulateFullModel.py for entries + +# Prepare time-resolved sendor finger print +timeSim = [] +timeSim = outputStruct[0]["timeSim"] +sensorFingerPrint = np.zeros([len(timeSim),len(sensorID)]) +for ii in range(len(sensorID)): + sensorFingerPrint[:,ii] = outputStruct[ii]["sensorFingerPrint"] + +# flagIndTime - If true, each time instant evaluated without knowledge of +# actual kinetics in the estimate of the concentration (accounted only in the +# generation of the true sensor response above) - This would lead to a +# scenario of concentrations evaluated at each time instant +if flagIndTime: + # Start time for time elapsed + startTime = time.time() + # Initialize output matrix + arrayConcentration = np.zeros([len(timeSim),numberOfGases+len(sensorID)]) + # Loop over all time instants and estimate the gas composition + arrayConcentrationTemp = Parallel(n_jobs=num_cores)(delayed(estimateConcentration) + (None,numberOfGases,None,sensorID, + fullModel = True, + fullModelResponse = sensorFingerPrint[ii,:]) + for ii in tqdm(range(len(timeSim)))) + # Convert the list to array + arrayConcentration = np.array(arrayConcentrationTemp) + # Stop time for time elapsed + stopTime = time.time() + # Time elapsed [s] + timeElapsed = stopTime - startTime +# flagIndTime - If false, only one concentration estimate obtained. The full +# model is simulated in concentration estimator accounting for the kinetics. +else: + # Start time for time elapsed + startTime = time.time() + # Initialize output matrix + arrayConcentration = np.zeros([len(timeSim)-1,numberOfGases+1+len(sensorID)]) + # Loop over all time instants and estimate the gas composition + # The concentration is estimated for all the time instants in timeSim. + # This simulates using the full model as and when a new measurement is + # available. This assumes that the estimator has knowledge of the past + # measurements + # The first time instant is not simulated. This is fixed by adding a dummy + # row after the simulations are over + arrayConcentrationTemp = Parallel(n_jobs=num_cores)(delayed(estimateConcentrationFullModel) + (numberOfGases,sensorID, + fullModelResponse = sensorFingerPrint[0:ii+1,:], + rateConstant = rateConstant, + timeInt = (0,timeSim[ii+1]),flowIn = flowIn) + for ii in tqdm(range(len(timeSim)-1))) + + # Convert the list to array + arrayConcentration = np.array(arrayConcentrationTemp) + # Add dummy row to be consistent (intialized to init condition) + firstRow = np.concatenate((np.array(sensorID), inputParameters[6])) + arrayConcentration = np.vstack([firstRow,arrayConcentration]) + # Stop time for time elapsed + stopTime = time.time() + # Time elapsed [s] + timeElapsed = stopTime - startTime + +# 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 = "fullModelConcentrationEstimate" +sensorText = str(sensorID).replace('[','').replace(']','').replace(' ','-').replace(',','') +saveFileName = filePrefix + "_" + sensorText + "_" + currentDT + "_" + gitCommitID; +savePath = os.path.join('simulationResults',saveFileName) +savez (savePath, outputStruct = outputStruct, # True response + arrayConcentration = arrayConcentration, # Estimated response + eqbmModelFlag = flagIndTime, # Flag to tell whether eqbm. or full model used + noiseCharacteristics = [meanNoise, stdNoise], # Noise characteristics + timeElapsed = timeElapsed, # Time elapsed for the simulation + hostName = socket.gethostname()) # Hostname of the computer \ No newline at end of file diff --git a/generateHypotheticalAdsorbents.py b/generateHypotheticalAdsorbents.py new file mode 100755 index 0000000..86fdd1f --- /dev/null +++ b/generateHypotheticalAdsorbents.py @@ -0,0 +1,86 @@ +############################################################################ +# +# Imperial College London, United Kingdom +# Multifunctional Nanomaterials Laboratory +# +# Project: ERASE +# Year: 2020 +# Python: Python 3.7 +# Authors: Ashwin Kumar Rajagopalan (AK) +# +# Purpose: +# Generates hypothetical sorbents using latin hypercube sampling. The +# sorbents are assumed to exhibit Langmuirian behavior. +# +# Last modified: +# - 2020-10-28, AK: Add auxiliary functions as a module +# - 2020-10-19, AK: Add adsorbent density and molecular weight +# - 2020-10-19, AK: Integrate git commit and save material properties +# - 2020-10-16, AK: Initial creation +# +# Input arguments: +# +# +# Output arguments: +# +# +############################################################################ + +def generateHypotheticalAdsorbents(numberOfGases, numberOfAdsorbents): + import numpy as np + from numpy import savez + import auxiliaryFunctions + from smt.sampling_methods import LHS + import os # For OS related stuff (make directory, file separator, etc.) + + # Get the commit ID of the current repository + gitCommitID = auxiliaryFunctions.getCommitID() + + # Get the current date and time for saving purposes + simulationDT = auxiliaryFunctions.getCurrentDateTime() + + # Define a range for single site Langmuir isotherm for one sensor material + # [qsat [mol/kg] b0 [m3/mol] delH [J/mol]] + singleIsothermRange = np.array([[0.0, 10.0], [0.0, 3e-6],[0, -40000]]) + + # Define a range for adsorbent densities for all the materials + densityRange = np.array([[500.0, 1500.0]]) + + # Define the molecular weight for the gases of interest + # [CO2 N2 O2 SO2 NO2 H2O] - HARD CODED + molecularWeight = np.array([44.01, 28.01, 15.99, 64.07, 46.01, 18.02]) + + # Initialize the output matrix with zeros + adsorbentIsotherm = np.zeros((numberOfGases,3,numberOfAdsorbents)) + + # Generate latin hypercube sampled hypothethical adsorbent materials. + # Isotherm parameters obtained as a matrix with dimensions [ + # numberOfAdsorbents*numberOfGases x 3] + samplingLHS = LHS(xlimits=singleIsothermRange) + allIsothermParameter = samplingLHS(numberOfAdsorbents*numberOfGases); + + # Also generate adsorbent densities that is latin hypercube sampled + samplingLHS = LHS(xlimits=densityRange) + adsorbentDensity = samplingLHS(numberOfAdsorbents) + + # Get the isotherms for each material for the predefined number of gases + isothermCounter = 0; + for ii in range(0,numberOfAdsorbents): + for jj in range(0,numberOfGases): + adsorbentIsotherm[jj,:,ii] = allIsothermParameter[isothermCounter,:] + isothermCounter += 1 + + # Save the adsorbent isotherm parameters into a native numpy file + # The .npz file is saved in a folder called inputResources (hardcoded) + filePrefix = "isothermParameters" + saveFileName = filePrefix + "_" + simulationDT + "_" + gitCommitID + ".npz"; + savePath = os.path.join('inputResources',saveFileName) + + # Check if inputResources directory exists or not. If not, create the folder + if not os.path.exists('inputResources'): + os.mkdir('inputResources') + + # Save the adsorbent material array + savez (savePath, adsIsotherm = adsorbentIsotherm, + adsDensity = adsorbentDensity, + molWeight = molecularWeight[0:numberOfGases]) \ No newline at end of file diff --git a/generateTrueSensorResponse.py b/generateTrueSensorResponse.py new file mode 100755 index 0000000..26c37d2 --- /dev/null +++ b/generateTrueSensorResponse.py @@ -0,0 +1,76 @@ +############################################################################ +# +# Imperial College London, United Kingdom +# Multifunctional Nanomaterials Laboratory +# +# Project: ERASE +# Year: 2020 +# Python: Python 3.7 +# Authors: Ashwin Kumar Rajagopalan (AK) +# +# Purpose: +# Generates a "true" sensor response for known mole fractions of gas for the +# sorbents that are hypothetically generated. The mole fraction is hard coded +# and should be modified. +# +# Last modified: +# - 2020-11-09, AK: Introduce custom mole fraction input +# - 2020-11-05, AK: Introduce new case for mole fraction sweep +# - 2020-10-30, AK: Fix to find number of gases +# - 2020-10-22, AK: Add two/three gases +# - 2020-10-21, AK: Initial creation +# +# Input arguments: +# +# +# Output arguments: +# +# +############################################################################ + +def generateTrueSensorResponse(numberOfAdsorbents, numberOfGases, pressureTotal, temperature, **kwargs): + import numpy as np + import pdb + from simulateSensorArray import simulateSensorArray + + # Mole fraction of the gas [-] + # Can be [jxg], where j is the number of mole fractions for g gases + if numberOfGases == 2: + moleFraction = np.array([[0.05, 0.95], + [0.15, 0.85], + [0.40, 0.60], + [0.75, 0.25], + [0.90, 0.10]]) + elif numberOfGases == 3: + moleFraction = np.array([[0.05, 0.15, 0.80], + [0.15, 0.25, 0.60], + [0.40, 0.35, 0.25], + [0.75, 0.10, 0.15], + [0.90, 0.05, 0.05]]) + # To sweep through the entire mole fraction range for 2 gases + elif numberOfGases == 20000: + moleFraction = np.array([np.linspace(0,1,1001), 1 - np.linspace(0,1,1001)]).T + # To sweep through the entire mole fraction range for 3 gases + elif numberOfGases == 30000: + moleFractionTemp = np.zeros([1001,3]) + num1 = np.random.uniform(0.0,1.0,1001) + num2 = np.random.uniform(0.0,1.0,1001) + num3 = np.random.uniform(0.0,1.0,1001) + sumNum = num1 + num2 + num3 + moleFractionTemp[:,0] = num1/sumNum + moleFractionTemp[:,1] = num2/sumNum + moleFractionTemp[:,2] = num3/sumNum + moleFraction = moleFractionTemp[moleFractionTemp[:,0].argsort()] + + # Check if a custom mole fraction is provided intead of the predefined one + if 'moleFraction' in kwargs: + moleFraction = np.array([kwargs["moleFraction"]]) + + # Get the individual sensor reponse for all the five "test" concentrations + sensorTrueResponse = np.zeros((numberOfAdsorbents,moleFraction.shape[0])) + for ii in range(numberOfAdsorbents): + for jj in range(moleFraction.shape[0]): + moleFractionTemp = np.array([moleFraction[jj,:]]) # This is needed to keep the structure as a row instead of column + sensorTrueResponse[ii,jj] = simulateSensorArray(np.array([ii]), + pressureTotal, temperature, moleFractionTemp) + return sensorTrueResponse \ No newline at end of file diff --git a/inputResources/isothermParameters_20201020_1635_5f263af.npz b/inputResources/isothermParameters_20201020_1635_5f263af.npz new file mode 100644 index 0000000..b1c8357 Binary files /dev/null and b/inputResources/isothermParameters_20201020_1635_5f263af.npz differ diff --git a/inputResources/isothermParameters_20201020_1756_5f263af.npz b/inputResources/isothermParameters_20201020_1756_5f263af.npz new file mode 100644 index 0000000..bf47fe4 Binary files /dev/null and b/inputResources/isothermParameters_20201020_1756_5f263af.npz differ diff --git a/inputResources/isothermParameters_20201022_1056_782efa3.npz b/inputResources/isothermParameters_20201022_1056_782efa3.npz new file mode 100644 index 0000000..47cf9d4 Binary files /dev/null and b/inputResources/isothermParameters_20201022_1056_782efa3.npz differ diff --git a/inputResources/isothermParameters_20201028_1756_3a07f95.npz b/inputResources/isothermParameters_20201028_1756_3a07f95.npz new file mode 100644 index 0000000..1e42d9a Binary files /dev/null and b/inputResources/isothermParameters_20201028_1756_3a07f95.npz differ diff --git a/inputResources/isothermParameters_20210120_1722_fb57143.npz b/inputResources/isothermParameters_20210120_1722_fb57143.npz new file mode 100644 index 0000000..4e9335b Binary files /dev/null and b/inputResources/isothermParameters_20210120_1722_fb57143.npz differ diff --git a/inputResources/isothermParameters_20210120_1724_fb57143.npz b/inputResources/isothermParameters_20210120_1724_fb57143.npz new file mode 100644 index 0000000..067256d Binary files /dev/null and b/inputResources/isothermParameters_20210120_1724_fb57143.npz differ diff --git a/plotFunctions/analyzeSensorResponse.py b/plotFunctions/analyzeSensorResponse.py new file mode 100755 index 0000000..60e4a82 --- /dev/null +++ b/plotFunctions/analyzeSensorResponse.py @@ -0,0 +1,644 @@ +############################################################################ +# +# Imperial College London, United Kingdom +# Multifunctional Nanomaterials Laboratory +# +# Project: ERASE +# Year: 2020 +# Python: Python 3.7 +# Authors: Ashwin Kumar Rajagopalan (AK) +# +# Purpose: +# +# +# Last modified: +# - 2020-11-05, AK: Initial creation +# +# Input arguments: +# +# +# Output arguments: +# +# +############################################################################ +import pdb +import numpy as np +from numpy import load +from generateTrueSensorResponse import generateTrueSensorResponse +from simulateSensorArray import simulateSensorArray +import os +import pandas as pd +import ternary +import matplotlib as mpl +import matplotlib.pyplot as plt +plt.style.use('doubleColumn.mplstyle') # Custom matplotlib style file +import auxiliaryFunctions + +os.chdir("..") + +# Save flag for figure +saveFlag = False + +# Save file extension (png or pdf) +saveFileExtension = ".png" + +# Plotting colors +colorsForPlot = ["eac435","345995","03cea4","fb4d3d","ca1551"] + +# Number of molefractions +numMolFrac= 101 + +# Total pressure of the gas [Pa] +pressureTotal = np.array([1.e5]); + +# Temperature of the gas [K] +# Can be a vector of temperatures +temperature = np.array([298.15]); + +# Number of Adsorbents +numberOfAdsorbents = 20 + +# Number of Gases +numberOfGases = 2 + +# Multiplier Error +multiplierError = [1., 1.] + +# Sensor ID +sensorID = np.array([6,2]) + +# Get the commit ID of the current repository +gitCommitID = auxiliaryFunctions.getCommitID() + +# Get the current date and time for saving purposes +currentDT = auxiliaryFunctions.getCurrentDateTime() + +# Simulate the sensor response for all possible concentrations +moleFractionRange = {} +moleFractionRange_y1 = {} +moleFractionRange_y2 = {} +# Two gases +if numberOfGases == 2: + moleFractionRange['y3_0'] = np.array([np.linspace(0,1,numMolFrac), 1 - np.linspace(0,1,numMolFrac)]).T +# Three gases - done by fixing one gas and going through combinations of +# y1 and y2 +elif numberOfGases == 3: + thirdGasMoleFrac = np.linspace(0,1.,numMolFrac) + for ii in range(len(thirdGasMoleFrac)): + tempDictKey = 'y3_{}'.format(ii) + remainingMoleFrac = 1. - thirdGasMoleFrac[ii] + moleFractionRange[tempDictKey] = np.array([np.linspace(0,remainingMoleFrac,numMolFrac), + remainingMoleFrac - np.linspace(0,remainingMoleFrac,numMolFrac), + np.tile(thirdGasMoleFrac[ii],numMolFrac)]).T + + secondGasMoleFrac = np.linspace(0,1.,numMolFrac) + for ii in range(len(secondGasMoleFrac)): + tempDictKey = 'y2_{}'.format(ii) + remainingMoleFrac = 1. - secondGasMoleFrac[ii] + moleFractionRange_y2[tempDictKey] = np.array([np.linspace(0,remainingMoleFrac,numMolFrac), + np.tile(secondGasMoleFrac[ii],numMolFrac), + remainingMoleFrac - np.linspace(0,remainingMoleFrac,numMolFrac)]).T + + firstGasMoleFrac = np.linspace(0,1.,numMolFrac) + for ii in range(len(firstGasMoleFrac)): + tempDictKey = 'y1_{}'.format(ii) + remainingMoleFrac = 1. - firstGasMoleFrac[ii] + moleFractionRange_y1[tempDictKey] = np.array([np.tile(firstGasMoleFrac[ii],numMolFrac), + np.linspace(0,remainingMoleFrac,numMolFrac), + remainingMoleFrac - np.linspace(0,remainingMoleFrac,numMolFrac)]).T + +# Get the simulated sensor response, condition number and the derivative of +# the mole fraction wrt the sensor response +arraySimResponse = {} +firstDerivativeSimResponse_y1 = {} +firstDerivativeSimResponse_y2 = {} +conditionNumber_y1 = {} +conditionNumber_y2 = {} +# Loop through all variables +for ii in range(len(moleFractionRange)): + tempDictKey = 'y3_{}'.format(ii) + # Temporary variable + moleFractionRangeTemp = moleFractionRange[tempDictKey] + arraySimResponseTemp = np.zeros([moleFractionRangeTemp.shape[0], + sensorID.shape[0]]) + firstDerivativeSimResponse_y1Temp = np.zeros([moleFractionRangeTemp.shape[0], + sensorID.shape[0]]) + firstDerivativeSimResponse_y2Temp = np.zeros([moleFractionRangeTemp.shape[0], + sensorID.shape[0]]) + conditionNumber_y1Temp = np.zeros([moleFractionRangeTemp.shape[0], + sensorID.shape[0]]) + conditionNumber_y2Temp = np.zeros([moleFractionRangeTemp.shape[0], + sensorID.shape[0]]) + # Loop through all mole fractions (y1 and y2) to get sensor response + for kk in range(moleFractionRangeTemp.shape[0]): + arraySimResponseTemp[kk,:] = simulateSensorArray(sensorID, pressureTotal, + temperature, np.array([moleFractionRangeTemp[kk,:]])) * multiplierError + + # Loop through all sensor responses to get the derivative and the condition + # number + for jj in range(arraySimResponseTemp.shape[1]): + firstDerivativeSimResponse_y1Temp[:,jj] = np.gradient(moleFractionRangeTemp[:,0], + arraySimResponseTemp[:,jj]) + # Condition number + conditionNumber_y1Temp[:,jj] = np.abs(np.multiply(np.divide(firstDerivativeSimResponse_y1Temp[:,jj], + moleFractionRangeTemp[:,0]), + arraySimResponseTemp[:,jj])) + firstDerivativeSimResponse_y2Temp[:,jj] = np.gradient(moleFractionRangeTemp[:,1], + arraySimResponseTemp[:,jj]) + # Condition number + conditionNumber_y2Temp[:,jj] = np.abs(np.multiply(np.divide(firstDerivativeSimResponse_y2Temp[:,jj], + moleFractionRangeTemp[:,1]), + arraySimResponseTemp[:,jj])) + # Save the sensor responseand the derivative to a dictionary + arraySimResponse[tempDictKey] = arraySimResponseTemp + firstDerivativeSimResponse_y1[tempDictKey] = firstDerivativeSimResponse_y1Temp + firstDerivativeSimResponse_y2[tempDictKey] = firstDerivativeSimResponse_y2Temp + conditionNumber_y1[tempDictKey] = conditionNumber_y1Temp + conditionNumber_y2[tempDictKey] = conditionNumber_y2Temp + arraySimResponseTemp = [] + firstDerivativeSimResponse_y1Temp = [] + firstDerivativeSimResponse_y2Temp = [] + conditionNumber_y1Temp = [] + conditionNumber_y2Temp = [] + +# For three gases +if numberOfGases == 3: + arraySimResponse_y2 = {} + firstDerivativeSimResponse_y31 = {} + firstDerivativeSimResponse_y3 = {} + for ii in range(len(moleFractionRange_y2)): + tempDictKey = 'y2_{}'.format(ii) + # Temporary variable + moleFractionRangeTemp = moleFractionRange_y2[tempDictKey] + arraySimResponseTemp = np.zeros([moleFractionRangeTemp.shape[0], + sensorID.shape[0]]) + firstDerivativeSimResponse_y31Temp = np.zeros([moleFractionRangeTemp.shape[0], + sensorID.shape[0]]) + firstDerivativeSimResponse_y3Temp = np.zeros([moleFractionRangeTemp.shape[0], + sensorID.shape[0]]) + + # Loop through all mole fractions (y1 and y2) to get sensor response + for kk in range(moleFractionRangeTemp.shape[0]): + arraySimResponseTemp[kk,:] = simulateSensorArray(sensorID, pressureTotal, + temperature, np.array([moleFractionRangeTemp[kk,:]])) * multiplierError + + # Loop through all sensor responses to get the derivative + for jj in range(arraySimResponseTemp.shape[1]): + firstDerivativeSimResponse_y31Temp[:,jj] = np.gradient(moleFractionRangeTemp[:,0],arraySimResponseTemp[:,jj]) + + firstDerivativeSimResponse_y3Temp[:,jj] = np.gradient(moleFractionRangeTemp[:,2], + arraySimResponseTemp[:,jj]) + + # Save the sensor responseand the derivative to a dictionary + arraySimResponse_y2[tempDictKey] = arraySimResponseTemp + firstDerivativeSimResponse_y31[tempDictKey] = firstDerivativeSimResponse_y31Temp + firstDerivativeSimResponse_y3[tempDictKey] = firstDerivativeSimResponse_y3Temp + arraySimResponseTemp = [] + firstDerivativeSimResponse_y31Temp = [] + firstDerivativeSimResponse_y3Temp = [] + + + # arraySimResponse_y1 = {} + # firstDerivativeSimResponse_y12 = {} + # firstDerivativeSimResponse_y13 = {} + # for ii in range(len(moleFractionRange_y1)): + # tempDictKey = 'y1_{}'.format(ii) + # # Temporary variable + # moleFractionRangeTemp = moleFractionRange_y1[tempDictKey] + # arraySimResponseTemp = np.zeros([moleFractionRangeTemp.shape[0], + # sensorID.shape[0]]) + # firstDerivativeSimResponse_y12Temp = np.zeros([moleFractionRangeTemp.shape[0], + # sensorID.shape[0]]) + # firstDerivativeSimResponse_y13Temp = np.zeros([moleFractionRangeTemp.shape[0], + # sensorID.shape[0]]) + + # # Loop through all mole fractions (y1 and y2) to get sensor response + # for kk in range(moleFractionRangeTemp.shape[0]): + # arraySimResponseTemp[kk,:] = simulateSensorArray(sensorID, pressureTotal, + # temperature, np.array([moleFractionRangeTemp[kk,:]])) * multiplierError + + # # Loop through all sensor responses to get the derivative + # for jj in range(arraySimResponseTemp.shape[1]): + # firstDerivativeSimResponse_y12Temp[:,jj] = np.gradient(moleFractionRangeTemp[:,1], + # arraySimResponseTemp[:,jj]) + + # firstDerivativeSimResponse_y13Temp[:,jj] = np.gradient(moleFractionRangeTemp[:,2], + # arraySimResponseTemp[:,jj]) + + # # Save the sensor responseand the derivative to a dictionary + # arraySimResponse_y1[tempDictKey] = arraySimResponseTemp + # firstDerivativeSimResponse_y12[tempDictKey] = firstDerivativeSimResponse_y12Temp + # firstDerivativeSimResponse_y13[tempDictKey] = firstDerivativeSimResponse_y13Temp + # arraySimResponseTemp = [] + # firstDerivativeSimResponse_y12Temp = [] + # firstDerivativeSimResponse_y13Temp = [] + +# Plot the sensor response and the derivative +# Two gases +if numberOfGases == 2: + fig = plt.figure + # Parse out mole fraction and sensor response + arraySimResponseTemp = arraySimResponse['y3_0'] + moleFractionRangeTemp = moleFractionRange['y3_0'] + firstDerivativeSimResponse_y1Temp = firstDerivativeSimResponse_y1['y3_0'] + conditionNumber_y1Temp = conditionNumber_y1['y3_0'] + + ax1 = plt.subplot(1,3,1) + # Loop through all sensors and plot the sensor response + for kk in range(arraySimResponseTemp.shape[1]): + ax1.plot(arraySimResponseTemp[:,kk],moleFractionRangeTemp[:,0], + color='#'+colorsForPlot[kk], label = '$s_'+str(kk+1)+'$') # Simulated Response + ax1.set(ylabel='$y_1$ [-]', + xlabel='$m_i$ [g kg$^{-1}$]', + ylim = [0,1], xlim = [0, 1.1*np.max(arraySimResponseTemp)]) + ax1.locator_params(axis="x", nbins=4) + ax1.locator_params(axis="y", nbins=4) + ax1.legend() + + ax2 = plt.subplot(1,3,2) + # Loop through all sensors and plot the sensor derivative + for kk in range(arraySimResponseTemp.shape[1]): + ax2.plot(arraySimResponseTemp[:,kk],conditionNumber_y1Temp[:,kk], + color='#'+colorsForPlot[kk], label = '$s_'+str(kk+1)+'$') # Simulated Response + ax2.set(xlabel='$m_i$ [g kg$^{-1}$]', + ylabel='$\chi$ [-]', + xlim = [0,1.1*np.max(arraySimResponseTemp)]) + ax2.locator_params(axis="x", nbins=4) + ax2.locator_params(axis="y", nbins=4) + ax2.legend() + + ax3 = plt.subplot(1,3,3) + # Loop through all sensors and plot the sensor derivative + for kk in range(arraySimResponseTemp.shape[1]): + ax3.plot(moleFractionRangeTemp[:,0],conditionNumber_y1Temp[:,kk], + color='#'+colorsForPlot[kk], label = '$s_'+str(kk+1)+'$') # Simulated Response + ax3.set(xlabel='$y_1$ [-]', + ylabel='$\chi$ [-]', + xlim = [0,1]) + ax3.locator_params(axis="x", nbins=4) + ax3.locator_params(axis="y", nbins=4) + ax3.legend() + + plt.show() + +# Three gases +elif numberOfGases == 3: + # Sensor response - y3 + fig, tax = ternary.figure(scale=1) + fig.set_size_inches(4,3.3) + tax.boundary(linewidth=1.0) + tax.gridlines(multiple=.2, color="gray") + currentMin = np.Inf + currentMax = -np.Inf + for ii in range(len(moleFractionRange)-1): + tempDictKey = 'y3_{}'.format(ii) + arraySimResponseTemp = arraySimResponse[tempDictKey] + responseMin = min(arraySimResponseTemp) + responseMax = max(arraySimResponseTemp) + currentMin = min([currentMin,responseMin]) + currentMax = max([currentMax,responseMax]) + + for ii in range(len(moleFractionRange)-1): + tempDictKey = 'y3_{}'.format(ii) + moleFractionRangeTemp = moleFractionRange[tempDictKey] + arraySimResponseTemp = arraySimResponse[tempDictKey] + if ii == len(moleFractionRange)-2: + tax.scatter(moleFractionRangeTemp, marker='o', s=2, c=arraySimResponseTemp, + vmin=currentMin,vmax=currentMax,cmap=plt.cm.PuOr, + colorbar = True, colormap=plt.cm.PuOr) + else: + tax.scatter(moleFractionRangeTemp, marker='o', s=2, c=arraySimResponseTemp, + vmin=currentMin,vmax=currentMax,cmap=plt.cm.PuOr) + + moleFractionRangeTemp = [] + arraySimResponseTemp = [] + + tax.left_axis_label("$y_2$ [-]",offset=0.20,fontsize=10) + tax.right_axis_label("$y_1$ [-]",offset=0.20,fontsize=10) + tax.bottom_axis_label("$y_3$ [-]",offset=0.20,fontsize=10) + tax.ticks(axis='lbr', linewidth=1, multiple=0.2, tick_formats="%.1f", + offset=0.035,clockwise=True,fontsize=10) + tax.clear_matplotlib_ticks() + tax._redraw_labels() + plt.axis('off') + tax.show() + + # Sensor response - y2 + fig, tax = ternary.figure(scale=1) + fig.set_size_inches(4,3.3) + tax.boundary(linewidth=1.0) + tax.gridlines(multiple=.2, color="gray") + currentMin = np.Inf + currentMax = -np.Inf + for ii in range(len(moleFractionRange_y2)-1): + tempDictKey = 'y2_{}'.format(ii) + arraySimResponseTemp = arraySimResponse_y2[tempDictKey] + responseMin = min(arraySimResponseTemp) + responseMax = max(arraySimResponseTemp) + currentMin = min([currentMin,responseMin]) + currentMax = max([currentMax,responseMax]) + + for ii in range(len(moleFractionRange)-1): + tempDictKey = 'y2_{}'.format(ii) + moleFractionRangeTemp = moleFractionRange_y2[tempDictKey] + arraySimResponseTemp = arraySimResponse_y2[tempDictKey] + if ii == len(moleFractionRange)-2: + tax.scatter(moleFractionRangeTemp, marker='o', s=2, c=arraySimResponseTemp, + vmin=currentMin,vmax=currentMax,cmap=plt.cm.PuOr, + colorbar = True, colormap=plt.cm.PuOr) + else: + tax.scatter(moleFractionRangeTemp, marker='o', s=2, c=arraySimResponseTemp, + vmin=currentMin,vmax=currentMax,cmap=plt.cm.PuOr) + + moleFractionRangeTemp = [] + arraySimResponseTemp = [] + + tax.left_axis_label("$y_2$ [-]",offset=0.20,fontsize=10) + tax.right_axis_label("$y_1$ [-]",offset=0.20,fontsize=10) + tax.bottom_axis_label("$y_3$ [-]",offset=0.20,fontsize=10) + tax.ticks(axis='lbr', linewidth=1, multiple=0.2, tick_formats="%.1f", + offset=0.035,clockwise=True,fontsize=10) + tax.clear_matplotlib_ticks() + tax._redraw_labels() + plt.axis('off') + tax.show() + + # # Sensor response - y2 + # fig, tax = ternary.figure(scale=1) + # fig.set_size_inches(4,3.3) + # tax.boundary(linewidth=1.0) + # tax.gridlines(multiple=.2, color="gray") + # currentMin = np.Inf + # currentMax = -np.Inf + # for ii in range(len(moleFractionRange_y1)-1): + # tempDictKey = 'y1_{}'.format(ii) + # arraySimResponseTemp = arraySimResponse_y1[tempDictKey] + # responseMin = min(arraySimResponseTemp) + # responseMax = max(arraySimResponseTemp) + # currentMin = min([currentMin,responseMin]) + # currentMax = max([currentMax,responseMax]) + + # for ii in range(len(moleFractionRange)-1): + # tempDictKey = 'y1_{}'.format(ii) + # moleFractionRangeTemp = moleFractionRange_y1[tempDictKey] + # arraySimResponseTemp = arraySimResponse_y1[tempDictKey] + # if ii == len(moleFractionRange)-2: + # tax.scatter(moleFractionRangeTemp, marker='o', s=2, c=arraySimResponseTemp, + # vmin=currentMin,vmax=currentMax,cmap=plt.cm.PuOr, + # colorbar = True, colormap=plt.cm.PuOr) + # else: + # tax.scatter(moleFractionRangeTemp, marker='o', s=2, c=arraySimResponseTemp, + # vmin=currentMin,vmax=currentMax,cmap=plt.cm.PuOr) + + # moleFractionRangeTemp = [] + # arraySimResponseTemp = [] + + # tax.left_axis_label("$y_2$ [-]",offset=0.20,fontsize=10) + # tax.right_axis_label("$y_1$ [-]",offset=0.20,fontsize=10) + # tax.bottom_axis_label("$y_3$ [-]",offset=0.20,fontsize=10) + # tax.ticks(axis='lbr', linewidth=1, multiple=0.2, tick_formats="%.1f", + # offset=0.035,clockwise=True,fontsize=10) + # tax.clear_matplotlib_ticks() + # tax._redraw_labels() + # plt.axis('off') + # tax.show() + + # Derivative - dy1/dm + fig, tax = ternary.figure(scale=1) + fig.set_size_inches(4,3.3) + tax.boundary(linewidth=1.0) + tax.gridlines(multiple=.2, color="gray") + currentMin = np.Inf + currentMax = -np.Inf + for ii in range(len(moleFractionRange)-1): + tempDictKey = 'y3_{}'.format(ii) + firstDerivativeSimResponse_y1Temp = firstDerivativeSimResponse_y1[tempDictKey] + firstDerivativeMin = min(firstDerivativeSimResponse_y1Temp[:,0]) + firstDerivativeMax = max(firstDerivativeSimResponse_y1Temp[:,0]) + currentMin = min([currentMin,firstDerivativeMin]) + currentMax = max([currentMax,firstDerivativeMax]) + + for ii in range(len(moleFractionRange)-1): + tempDictKey = 'y3_{}'.format(ii) + moleFractionRangeTemp = moleFractionRange[tempDictKey] + firstDerivativeSimResponse_y1Temp = firstDerivativeSimResponse_y1[tempDictKey] + if ii == len(moleFractionRange)-2: + tax.scatter(moleFractionRangeTemp, marker='o', s=2, c=firstDerivativeSimResponse_y1Temp[:,0], + vmin=currentMin,vmax=currentMax,cmap=plt.cm.PuOr, + colorbar = True, colormap=plt.cm.PuOr) + else: + tax.scatter(moleFractionRangeTemp, marker='o', s=2, c=firstDerivativeSimResponse_y1Temp[:,0], + vmin=currentMin,vmax=currentMax,cmap=plt.cm.PuOr) + + moleFractionRangeTemp = [] + firstDerivativeSimResponse_y1Temp = [] + + tax.left_axis_label("$y_2$ [-]",offset=0.20,fontsize=10) + tax.right_axis_label("$y_1$ [-]",offset=0.20,fontsize=10) + tax.bottom_axis_label("$y_3$ [-]",offset=0.20,fontsize=10) + tax.ticks(axis='lbr', linewidth=1, multiple=0.2, tick_formats="%.1f", + offset=0.035,clockwise=True,fontsize=10) + tax.clear_matplotlib_ticks() + tax._redraw_labels() + plt.axis('off') + tax.show() + + # Derivative - dy2/dm + fig, tax = ternary.figure(scale=1) + fig.set_size_inches(4,3.3) + tax.boundary(linewidth=1.0) + tax.gridlines(multiple=.2, color="gray") + currentMin = np.Inf + currentMax = -np.Inf + for ii in range(len(moleFractionRange)-1): + tempDictKey = 'y3_{}'.format(ii) + firstDerivativeSimResponse_y2Temp = firstDerivativeSimResponse_y2[tempDictKey] + firstDerivativeMin = min(firstDerivativeSimResponse_y2Temp[:,0]) + firstDerivativeMax = max(firstDerivativeSimResponse_y2Temp[:,0]) + currentMin = min([currentMin,firstDerivativeMin]) + currentMax = max([currentMax,firstDerivativeMax]) + + for ii in range(len(moleFractionRange)-1): + tempDictKey = 'y3_{}'.format(ii) + moleFractionRangeTemp = moleFractionRange[tempDictKey] + firstDerivativeSimResponse_y2Temp = firstDerivativeSimResponse_y2[tempDictKey] + if ii == len(moleFractionRange)-2: + tax.scatter(moleFractionRangeTemp, marker='o', s=2, c=firstDerivativeSimResponse_y2Temp[:,0], + vmin=currentMin,vmax=currentMax,cmap=plt.cm.PuOr, + colorbar = True, colormap=plt.cm.PuOr) + else: + tax.scatter(moleFractionRangeTemp, marker='o', s=2, c=firstDerivativeSimResponse_y2Temp[:,0], + vmin=currentMin,vmax=currentMax,cmap=plt.cm.PuOr) + + moleFractionRangeTemp = [] + firstDerivativeSimResponse_y2Temp = [] + + tax.left_axis_label("$y_2$ [-]",offset=0.20,fontsize=10) + tax.right_axis_label("$y_1$ [-]",offset=0.20,fontsize=10) + tax.bottom_axis_label("$y_3$ [-]",offset=0.20,fontsize=10) + tax.ticks(axis='lbr', linewidth=1, multiple=0.2, tick_formats="%.1f", + offset=0.035,clockwise=True,fontsize=10) + tax.clear_matplotlib_ticks() + tax._redraw_labels() + plt.axis('off') + tax.show() + + # Derivative - dy3/dm + fig, tax = ternary.figure(scale=1) + fig.set_size_inches(4,3.3) + tax.boundary(linewidth=1.0) + tax.gridlines(multiple=.2, color="gray") + currentMin = np.Inf + currentMax = -np.Inf + for ii in range(len(moleFractionRange_y2)-1): + tempDictKey = 'y2_{}'.format(ii) + firstDerivativeSimResponse_y3Temp = firstDerivativeSimResponse_y3[tempDictKey] + firstDerivativeMin = min(firstDerivativeSimResponse_y3Temp[:,0]) + firstDerivativeMax = max(firstDerivativeSimResponse_y3Temp[:,0]) + currentMin = min([currentMin,firstDerivativeMin]) + currentMax = max([currentMax,firstDerivativeMax]) + + for ii in range(len(moleFractionRange_y2)-1): + tempDictKey = 'y2_{}'.format(ii) + moleFractionRangeTemp = moleFractionRange_y2[tempDictKey] + firstDerivativeSimResponse_y3Temp = firstDerivativeSimResponse_y3[tempDictKey] + if ii == len(moleFractionRange)-2: + tax.scatter(moleFractionRangeTemp, marker='o', s=2, c=firstDerivativeSimResponse_y3Temp[:,0], + vmin=currentMin,vmax=currentMax,cmap=plt.cm.PuOr, + colorbar = True, colormap=plt.cm.PuOr) + else: + tax.scatter(moleFractionRangeTemp, marker='o', s=2, c=firstDerivativeSimResponse_y3Temp[:,0], + vmin=currentMin,vmax=currentMax,cmap=plt.cm.PuOr) + + moleFractionRangeTemp = [] + firstDerivativeSimResponse_y3Temp = [] + + tax.left_axis_label("$y_2$ [-]",offset=0.20,fontsize=10) + tax.right_axis_label("$y_1$ [-]",offset=0.20,fontsize=10) + tax.bottom_axis_label("$y_3$ [-]",offset=0.20,fontsize=10) + tax.ticks(axis='lbr', linewidth=1, multiple=0.2, tick_formats="%.1f", + offset=0.035,clockwise=True,fontsize=10) + tax.clear_matplotlib_ticks() + tax._redraw_labels() + plt.axis('off') + tax.show() + + # Derivative - dy3/dm + fig, tax = ternary.figure(scale=1) + fig.set_size_inches(4,3.3) + tax.boundary(linewidth=1.0) + tax.gridlines(multiple=.2, color="gray") + currentMin = np.Inf + currentMax = -np.Inf + for ii in range(len(moleFractionRange_y2)-1): + tempDictKey = 'y2_{}'.format(ii) + firstDerivativeSimResponse_y31Temp = firstDerivativeSimResponse_y31[tempDictKey] + firstDerivativeMin = min(firstDerivativeSimResponse_y31Temp[:,0]) + firstDerivativeMax = max(firstDerivativeSimResponse_y31Temp[:,0]) + currentMin = min([currentMin,firstDerivativeMin]) + currentMax = max([currentMax,firstDerivativeMax]) + + for ii in range(len(moleFractionRange_y2)-1): + tempDictKey = 'y2_{}'.format(ii) + moleFractionRangeTemp = moleFractionRange_y2[tempDictKey] + firstDerivativeSimResponse_y31Temp = firstDerivativeSimResponse_y31[tempDictKey] + if ii == len(moleFractionRange)-2: + tax.scatter(moleFractionRangeTemp, marker='o', s=2, c=firstDerivativeSimResponse_y31Temp[:,0], + vmin=currentMin,vmax=currentMax,cmap=plt.cm.PuOr, + colorbar = True, colormap=plt.cm.PuOr) + else: + tax.scatter(moleFractionRangeTemp, marker='o', s=2, c=firstDerivativeSimResponse_y31Temp[:,0], + vmin=currentMin,vmax=currentMax,cmap=plt.cm.PuOr) + + moleFractionRangeTemp = [] + firstDerivativeSimResponse_y31Temp = [] + + tax.left_axis_label("$y_2$ [-]",offset=0.20,fontsize=10) + tax.right_axis_label("$y_1$ [-]",offset=0.20,fontsize=10) + tax.bottom_axis_label("$y_3$ [-]",offset=0.20,fontsize=10) + tax.ticks(axis='lbr', linewidth=1, multiple=0.2, tick_formats="%.1f", + offset=0.035,clockwise=True,fontsize=10) + tax.clear_matplotlib_ticks() + tax._redraw_labels() + plt.axis('off') + tax.show() + + # # Derivative - dy3/dm + # fig, tax = ternary.figure(scale=1) + # fig.set_size_inches(4,3.3) + # tax.boundary(linewidth=1.0) + # tax.gridlines(multiple=.2, color="gray") + # currentMin = np.Inf + # currentMax = -np.Inf + # for ii in range(len(moleFractionRange_y1)-1): + # tempDictKey = 'y1_{}'.format(ii) + # firstDerivativeSimResponse_y12Temp = firstDerivativeSimResponse_y12[tempDictKey] + # firstDerivativeMin = min(firstDerivativeSimResponse_y12Temp[:,0]) + # firstDerivativeMax = max(firstDerivativeSimResponse_y12Temp[:,0]) + # currentMin = min([currentMin,firstDerivativeMin]) + # currentMax = max([currentMax,firstDerivativeMax]) + + # for ii in range(len(moleFractionRange_y1)-1): + # tempDictKey = 'y1_{}'.format(ii) + # moleFractionRangeTemp = moleFractionRange_y1[tempDictKey] + # firstDerivativeSimResponse_y12Temp = firstDerivativeSimResponse_y12[tempDictKey] + # if ii == len(moleFractionRange)-2: + # tax.scatter(moleFractionRangeTemp, marker='o', s=2, c=firstDerivativeSimResponse_y12Temp[:,0], + # vmin=currentMin,vmax=currentMax,cmap=plt.cm.PuOr, + # colorbar = True, colormap=plt.cm.PuOr) + # else: + # tax.scatter(moleFractionRangeTemp, marker='o', s=2, c=firstDerivativeSimResponse_y12Temp[:,0], + # vmin=currentMin,vmax=currentMax,cmap=plt.cm.PuOr) + + # moleFractionRangeTemp = [] + # firstDerivativeSimResponse_y12Temp = [] + + # tax.left_axis_label("$y_2$ [-]",offset=0.20,fontsize=10) + # tax.right_axis_label("$y_1$ [-]",offset=0.20,fontsize=10) + # tax.bottom_axis_label("$y_3$ [-]",offset=0.20,fontsize=10) + # tax.ticks(axis='lbr', linewidth=1, multiple=0.2, tick_formats="%.1f", + # offset=0.035,clockwise=True,fontsize=10) + # tax.clear_matplotlib_ticks() + # tax._redraw_labels() + # plt.axis('off') + # tax.show() + + # # Derivative - dy3/dm + # fig, tax = ternary.figure(scale=1) + # fig.set_size_inches(4,3.3) + # tax.boundary(linewidth=1.0) + # tax.gridlines(multiple=.2, color="gray") + # currentMin = np.Inf + # currentMax = -np.Inf + # for ii in range(len(moleFractionRange_y1)-1): + # tempDictKey = 'y1_{}'.format(ii) + # firstDerivativeSimResponse_y13Temp = firstDerivativeSimResponse_y13[tempDictKey] + # firstDerivativeMin = min(firstDerivativeSimResponse_y13Temp[:,0]) + # firstDerivativeMax = max(firstDerivativeSimResponse_y13Temp[:,0]) + # currentMin = min([currentMin,firstDerivativeMin]) + # currentMax = max([currentMax,firstDerivativeMax]) + + # for ii in range(len(moleFractionRange_y1)-1): + # tempDictKey = 'y1_{}'.format(ii) + # moleFractionRangeTemp = moleFractionRange_y1[tempDictKey] + # firstDerivativeSimResponse_y13Temp = firstDerivativeSimResponse_y13[tempDictKey] + # if ii == len(moleFractionRange)-2: + # tax.scatter(moleFractionRangeTemp, marker='o', s=2, c=firstDerivativeSimResponse_y13Temp[:,0], + # vmin=currentMin,vmax=currentMax,cmap=plt.cm.PuOr, + # colorbar = True, colormap=plt.cm.PuOr) + # else: + # tax.scatter(moleFractionRangeTemp, marker='o', s=2, c=firstDerivativeSimResponse_y13Temp[:,0], + # vmin=currentMin,vmax=currentMax,cmap=plt.cm.PuOr) + + # moleFractionRangeTemp = [] + # firstDerivativeSimResponse_y13Temp = [] + + # tax.left_axis_label("$y_2$ [-]",offset=0.20,fontsize=10) + # tax.right_axis_label("$y_1$ [-]",offset=0.20,fontsize=10) + # tax.bottom_axis_label("$y_3$ [-]",offset=0.20,fontsize=10) + # tax.ticks(axis='lbr', linewidth=1, multiple=0.2, tick_formats="%.1f", + # offset=0.035,clockwise=True,fontsize=10) + # tax.clear_matplotlib_ticks() + # tax._redraw_labels() + # plt.axis('off') + # tax.show() \ No newline at end of file diff --git a/plotFunctions/doubleColumn.mplstyle b/plotFunctions/doubleColumn.mplstyle new file mode 100755 index 0000000..8a5f6bc --- /dev/null +++ b/plotFunctions/doubleColumn.mplstyle @@ -0,0 +1,44 @@ +# matplotlib configuration file for the ERASE project +# https://matplotlib.org/3.3.2/tutorials/introductory/customizing.html + +## Figure property +figure.figsize : 7, 3 # width, height in inches +figure.dpi : 600 # dpi +figure.autolayout : true # for labels not being cut out + +## Axes +axes.titlesize : 10 +axes.labelsize : 10 +axes.formatter.limits : -5, 3 + +## Grid +axes.grid : true +grid.color : cccccc +grid.linewidth : 0.5 + +## Lines & Scatter +lines.linewidth : 1.5 +lines.markersize : 4 +scatter.marker: o + +## Ticks +xtick.top : true +ytick.right : true +xtick.labelsize : 8 +ytick.labelsize : 8 + +## Fonts +font.family : arial +font.weight : normal +font.size : 10 + +## Legends +legend.frameon : true +legend.fontsize : 10 +legend.edgecolor : 1 +legend.framealpha : 0.6 + +## Save figure +savefig.dpi : figure +savefig.format : png +savefig.transparent : false \ No newline at end of file diff --git a/plotFunctions/doubleColumn2Row.mplstyle b/plotFunctions/doubleColumn2Row.mplstyle new file mode 100644 index 0000000..d71b19b --- /dev/null +++ b/plotFunctions/doubleColumn2Row.mplstyle @@ -0,0 +1,44 @@ +# matplotlib configuration file for the ERASE project +# https://matplotlib.org/3.3.2/tutorials/introductory/customizing.html + +## Figure property +figure.figsize : 7, 6 # width, height in inches +figure.dpi : 600 # dpi +figure.autolayout : true # for labels not being cut out + +## Axes +axes.titlesize : 10 +axes.labelsize : 10 +axes.formatter.limits : -5, 4 + +## Grid +axes.grid : true +grid.color : cccccc +grid.linewidth : 0.5 + +## Lines & Scatter +lines.linewidth : 1.5 +lines.markersize : 4 +scatter.marker: o + +## Ticks +xtick.top : true +ytick.right : true +xtick.labelsize : 8 +ytick.labelsize : 8 + +## Fonts +font.family : arial +font.weight : normal +font.size : 10 + +## Legends +legend.frameon : true +legend.fontsize : 10 +legend.edgecolor : 1 +legend.framealpha : 0.6 + +## Save figure +savefig.dpi : figure +savefig.format : png +savefig.transparent : false \ No newline at end of file diff --git a/plotFunctions/mapTrueSensorResponse.py b/plotFunctions/mapTrueSensorResponse.py new file mode 100755 index 0000000..f0fef09 --- /dev/null +++ b/plotFunctions/mapTrueSensorResponse.py @@ -0,0 +1,164 @@ +############################################################################ +# +# Imperial College London, United Kingdom +# Multifunctional Nanomaterials Laboratory +# +# Project: ERASE +# Year: 2020 +# Python: Python 3.7 +# Authors: Ashwin Kumar Rajagopalan (AK) +# +# Purpose: +# Maps the sensor response as a function of mole fraction and the sorbent +# used +# +# Last modified: +# - 2020-11-06, AK: Cosmetic changes +# - 2020-10-23, AK: Initial creation +# +# Input arguments: +# +# +# Output arguments: +# +# +############################################################################ + +import numpy as np +from generateTrueSensorResponse import generateTrueSensorResponse +import matplotlib.pyplot as plt +import auxiliaryFunctions +from scipy import stats +plt.style.use('doubleColumn.mplstyle') # Custom matplotlib style file +import os +os.chdir("..") + +# Histogram properties +nBins = 100 +rangeX = ([-10,10]) +histTypeX = 'step' +alphaX=0.5 +densityX = False + + +# Flag for saving figures +saveFlag = False + +# Save file extension (png or pdf) +saveFileExtension = ".png" + +# Sensors to plot +sensorID = range(0,20) + +# Sensor ID for plotting histogram +histSensorID = 16 + +# Total pressure of the gas [Pa] +numberOfGases = 20000 # 20000 does a complete concentraiton sweep for 2/3 gases + +# Mole fraction to be plotted +moleFraction = np.linspace(0,1,1001) + +# Total pressure of the gas [Pa] +pressureTotal = np.array([1.e5]) + +# Temperature of the gas [K] +# Can be a vector of temperatures +temperature = np.array([298.15]) + +# Total number of sensor elements/gases simulated and generated using +# generateHypotheticalAdsorbents.py function +numberOfAdsorbents = 20 + +# Generate the true sensor response +sensorTrueResponse = generateTrueSensorResponse(numberOfAdsorbents, numberOfGases, + pressureTotal, temperature) + +# Parse out sensors that need to be plotted +sensorPlotResponse = np.zeros([len(sensorID),sensorTrueResponse.shape[1]]) +for jj in range(len(sensorID)): + for ii in range(sensorTrueResponse.shape[1]): + sensorPlotResponse[jj,ii] = (sensorTrueResponse[sensorID[jj],ii]) + +# Perform a standard normal variate of the sensor output +meanResponse = np.mean(sensorPlotResponse,axis=1) +stdResponse = np.std(sensorPlotResponse,axis=1) + +sensorSNVResponse = np.zeros([len(sensorID),sensorTrueResponse.shape[1]]) +for jj in range(len(sensorID)): + sensorSNVResponse[jj,:] = (sensorPlotResponse[jj,:] - meanResponse[jj])/stdResponse[jj] + +# Perform a Box-Cpx transformation +lambdaParam = np.zeros(len(sensorID)) +for jj in range(len(sensorID)): + _ , lambdaParam[jj] = stats.boxcox(sensorPlotResponse[jj,:]) + +# Calculate the skewness of the isotherm +skewnessSNVResponse = np.zeros(len(sensorID)) +for jj in range(len(sensorID)): + skewnessNum = np.sum(np.power(sensorSNVResponse[jj,:] - np.mean(sensorSNVResponse[jj,:]),3)) + skewnessDen = (len(sensorSNVResponse[jj,:]))*np.power(np.std(sensorSNVResponse[jj,:]),3) + skewnessSNVResponse[jj] = skewnessNum/skewnessDen + +# Calculate the kurtosis of the isotherm +kurtosisSNVResponse = np.zeros(len(sensorID)) +for jj in range(len(sensorID)): + kurtosisNum = np.sum(np.power(sensorSNVResponse[jj,:] - np.mean(sensorSNVResponse[jj,:]),4)) + kurtosisDen = (len(sensorSNVResponse[jj,:]))*np.power(np.std(sensorSNVResponse[jj,:]),4) + kurtosisSNVResponse[jj] = kurtosisNum/kurtosisDen + +# Get the commit ID of the current repository +gitCommitID = auxiliaryFunctions.getCommitID() + +# Get the current date and time for saving purposes +currentDT = auxiliaryFunctions.getCurrentDateTime() + +# Plot the sensor finger print for different concentrations of gases +fig = plt.figure +ax = plt.subplot(1,2,1) +for ii in range(len(sensorID)): + plotX = np.ones(sensorTrueResponse.shape[1])*sensorID[ii] + plotY = sensorPlotResponse[ii,:] + s1 = ax.scatter(plotX,plotY,c = moleFraction, cmap='RdYlBu') + +ax.set(xlabel='Adsorbent ID [-]', + ylabel='$m_i$ [g kg$^{-1}$]', + xlim = [0, 20], ylim = [0,300]) +ax.locator_params(axis="x", nbins=4) +ax.locator_params(axis="y", nbins=4) +plt.colorbar(s1,ax=ax,label="$y_1$ [-]") + +# Save the figure +if saveFlag: + # FileName: SensorMapRaw__ + saveFileName = "SensorMapRaw_" + currentDT + "_" + gitCommitID + saveFileExtension + savePath = os.path.join('simulationFigures',saveFileName) + # Check if simulationFigures directory exists or not. If not, create the folder + if not os.path.exists(os.path.join('..','simulationFigures')): + os.mkdir(os.path.join('simulationFigures')) + plt.savefig (savePath) + +# Plot the sensor finger print for different concentrations of gases, but with +# mean centering and normalizing with standard deviation +ax = plt.subplot(1,2,2) +for ii in range(len(sensorID)): + plotX = np.ones(sensorTrueResponse.shape[1])*sensorID[ii] + plotY = sensorSNVResponse[ii,:] + s2 = ax.scatter(plotX,plotY,c = moleFraction, cmap='RdYlBu') +ax.set(xlabel='Adsorbent ID [-]', + ylabel='$\hat{m}_i$ [-]', + xlim = [0, 20], ylim = [-10,5]) +ax.locator_params(axis="x", nbins=4) +ax.locator_params(axis="y", nbins=4) +plt.colorbar(s2,ax=ax,label="$y_1$ [-]") +# Save the figure +if saveFlag: + # FileName: SensorMapSNV__ + saveFileName = "SensorMapSNV_" + currentDT + "_" + gitCommitID + saveFileExtension + savePath = os.path.join('simulationFigures',saveFileName) + # Check if simulationFigures directory exists or not. If not, create the folder + if not os.path.exists(os.path.join('..','simulationFigures')): + os.mkdir(os.path.join('simulationFigures')) + plt.savefig (savePath) + +plt.show() \ No newline at end of file diff --git a/plotFunctions/plotAdsorbentCharacteristics.py b/plotFunctions/plotAdsorbentCharacteristics.py new file mode 100755 index 0000000..fd2c91d --- /dev/null +++ b/plotFunctions/plotAdsorbentCharacteristics.py @@ -0,0 +1,109 @@ +############################################################################ +# +# Imperial College London, United Kingdom +# Multifunctional Nanomaterials Laboratory +# +# Project: ERASE +# Year: 2020 +# Python: Python 3.7 +# Authors: Ashwin Kumar Rajagopalan (AK) +# +# Purpose: +# Plots the adsorbent properties for the hypothetical materials +# +# Last modified: +# - 2020-10-28, AK: Minor fix for save file name +# - 2020-10-28, AK: Add auxiliary functions as a module +# - 2020-10-27, AK: Initial creation +# +# Input arguments: +# +# +# Output arguments: +# +# +############################################################################ + +import os +from numpy import load +import auxiliaryFunctions +import matplotlib.pyplot as plt +plt.style.use('doubleColumn.mplstyle') # Custom matplotlib style file + +# Flag for saving figure +saveFlag = False + +# Save file extension (png or pdf) +saveFileExtension = ".pdf" + +# Select the id of gas that is to be plotted +gasID = 0 + +# For now load a given adsorbent isotherm material file +# loadFileName = "isothermParameters_20201020_1756_5f263af.npz" # Two gases +loadFileName = "isothermParameters_20201022_1056_782efa3.npz" # Three gases +hypoAdsorbentFile = os.path.join('..','inputResources',loadFileName); + +# Check if the file with the adsorbent properties exist +if os.path.exists(hypoAdsorbentFile): + loadedFileContent = load(hypoAdsorbentFile) + adsorbentIsothermTemp = loadedFileContent['adsIsotherm'] + adsorbentIsotherm = adsorbentIsothermTemp[gasID,:,:] + adsorbentDensity = loadedFileContent['adsDensity'] + molecularWeight = loadedFileContent['molWeight'] +else: + errorString = "Adsorbent property file " + hypoAdsorbentFile + " does not exist." + raise Exception(errorString) + +# Get the commit ID of the current repository +gitCommitID = auxiliaryFunctions.getCommitID() + +# Get the current date and time for saving purposes +currentDT = auxiliaryFunctions.getCurrentDateTime() + +# Git commit id of the loaded isotherm file +gitCommmitID_loadedFile = hypoAdsorbentFile[-11:-4] + +# Scatter plot for qsat vs b0, delH, and rho +colorVar = range(1,101) +fig = plt.figure + +# qsat vs b0 +ax1 = plt.subplot(1,3,1) +s1 = ax1.scatter(adsorbentIsotherm[0,:], adsorbentIsotherm[1,:], c = colorVar, cmap='RdYlBu') +ax1.set(xlabel='$q_\mathregular{sat}$ [mol kg$^{\mathregular{-1}}$]', + ylabel='$b_\mathregular{0}$ [m$^{\mathregular{3}}$ mol$^{\mathregular{-1}}$]', + xlim = [0, 10], ylim = [0, 3e-6]) +ax1.locator_params(axis="x", nbins=4) +ax1.locator_params(axis="y", nbins=4) + +# qsat vs delH +ax2 = plt.subplot(1,3,2) +s2 = ax2.scatter(adsorbentIsotherm[0,:], -adsorbentIsotherm[2,:], c = colorVar, cmap='RdYlBu') +ax2.set(xlabel='$q_\mathregular{sat}$ [mol kg$^{\mathregular{-1}}$]', + ylabel='-$\Delta H$ [J mol$^{\mathregular{-1}}$]', + xlim = [0, 10], ylim = [0, 4e4]) +ax2.locator_params(axis="x", nbins=4) +ax2.locator_params(axis="y", nbins=4) + +# qsat vs rho +ax3 = plt.subplot(1,3,3) +s3 = ax3.scatter(adsorbentIsotherm[0,:], adsorbentDensity, c = colorVar, cmap='RdYlBu') +ax3.set(xlabel='$q_\mathregular{sat}$ [mol kg$^{\mathregular{-1}}$]', + ylabel='$\\rho$ [kg m$^{\mathregular{-3}}$]', + xlim = [0, 10], ylim = [500, 1500]) +ax3.locator_params(axis="x", nbins=4) +ax3.locator_params(axis="y", nbins=4) + +# Save the figure +if saveFlag: + # FileName: AdsCharac____ + saveFileName = "AdsCharac_" + str(gasID) + "_" + currentDT + "_" + gitCommitID + "_" + gitCommmitID_loadedFile + saveFileExtension + savePath = os.path.join('..','simulationFigures',saveFileName) + # Check if inputResources directory exists or not. If not, create the folder + if not os.path.exists(os.path.join('..','simulationFigures')): + os.mkdir(os.path.join('..','simulationFigures')) + plt.savefig (savePath) + +# For the figure to be saved show should appear after the save +plt.show() \ No newline at end of file diff --git a/plotFunctions/plotConcentrationEstimate.py b/plotFunctions/plotConcentrationEstimate.py new file mode 100755 index 0000000..f8215ef --- /dev/null +++ b/plotFunctions/plotConcentrationEstimate.py @@ -0,0 +1,264 @@ +############################################################################ +# +# Imperial College London, United Kingdom +# Multifunctional Nanomaterials Laboratory +# +# Project: ERASE +# Year: 2020 +# Python: Python 3.7 +# Authors: Ashwin Kumar Rajagopalan (AK) +# +# Purpose: +# Plots to visualize estimated concentration +# +# Last modified: +# - 2020-11-17, AK: Cosmetic changes +# - 2020-11-09, AK: Add functionality for .npz file +# - 2020-11-09, AK: Cosmetic changes +# - 2020-11-04, AK: Improve plotting capability for three gases/sensors +# - 2020-10-30, AK: Add zoomed in version +# - 2020-10-29, AK: Improvements to the plots +# - 2020-10-23, AK: Initial creation +# +# Input arguments: +# +# +# Output arguments: +# +# +############################################################################ + +import numpy as np +from numpy import load +import os +import auxiliaryFunctions +import matplotlib.pyplot as plt +plt.style.use('singleColumn.mplstyle') # Custom matplotlib style file + +# Number of gases +numberOfGases = 2 + +# Flag for saving figure +saveFlag = False + +# Save file extension (png or pdf) +saveFileExtension = ".png" + +# Plot zoomed plot +plotZoomDist = False + +# Mole frac ID (use this if only .npz is used - sensitivity studies) +moleFracID = 0 + +# Gas concentration +molFracG1 = 0.90 +molFracG2 = 0.10 +molFracG3 = 1 - molFracG1 - molFracG2 + +# Xlimits and Ylimits +xLimits = [0,1] +yLimits = [0,500] +xLimitsSum = [0,2] +yLimitsSum = [0,300] + +xLimitsZ = [0,1] # Limits for zoom distribution +yLimitsZ = [0,15] # Limits for zoom distribution +xLimitsSumZ = [0,2] # Limits for zoom gas sum +yLimitsSumZ = [0,15] # Limits for zoom gas sum + +# Histogram properties +nBins = 100 +rangeX = (xLimits) +rangeXS = (xLimitsSum) +histTypeX = 'stepfilled' +alphaX=0.5 +densityX = False + +# For now load a given adsorbent isotherm material file +# loadFileName = "arrayConcentration_20201030_1109_5c77a62.npy" # 3 gases, 1 sensor +# loadFileName = "arrayConcentration_20201030_0913_5c77a62.npy" # 3 gases, 2 sensors +# loadFileName = "arrayConcentration_20201029_2328_5c77a62.npy" # 3 gases, 3 sensors +# loadFileName = "arrayConcentration_20201102_1423_da1707b.npy" # 2 gases, 1 sensor +# loadFileName = "arrayConcentration_20201104_1732_cc08dc4.npy" # 2 gases, 2 sensor [0.15, 0.85] +# loadFileName = "arrayConcentration_20201030_1731_da1707b.npy" # 2 gases, 2 sensor [0.4, 0.6] +# loadFileName = "arrayConcentration_20201104_1842_cc08dc4.npy" # 2 gases, 2 sensor [0.75, 0.25] +# loadFileName = "arrayConcentration_20201104_2227_cc08dc4.npy" # 2 gases, 2 sensor [0.75, 0.25] +loadFileName = "sensitivityAnalysis_17-16_20201114_2032_c9b2a41.npz" +simResultsFile = os.path.join('..','simulationResults',loadFileName); + +# Get the commit ID of the current repository +gitCommitID = auxiliaryFunctions.getCommitID() + +# Get the current date and time for saving purposes +currentDT = auxiliaryFunctions.getCurrentDateTime() + +# Git commit id of the loaded isotherm file +simID_loadedFile = loadFileName[-21:-4] + +# Check if the file is npz and parse out the necessary information +loadFileExt = loadFileName[-3:] +if loadFileExt == 'npz': + fileNPZ = True +else: + fileNPZ = False + +# Check if the file with the adsorbent properties exist +if os.path.exists(simResultsFile): + if fileNPZ: + resultOutput = load(simResultsFile)["arrayConcentration"][moleFracID,:,:] + else: + resultOutput = load(simResultsFile) + if numberOfGases == 2: + if resultOutput.shape[1] == 3: + resultOutput = np.delete(resultOutput,[0],1) + elif resultOutput.shape[1] == 4: + resultOutput = np.delete(resultOutput,[0,1],1) + elif numberOfGases == 3: + if resultOutput.shape[1] == 4: + resultOutput = np.delete(resultOutput,0,1) + elif resultOutput.shape[1] == 5: + resultOutput = np.delete(resultOutput,[0,1],1) + elif resultOutput.shape[1] == 6: + resultOutput = np.delete(resultOutput,[0,1,2],1) +else: + errorString = "Simulation result file " + simResultsFile + " does not exist." + raise Exception(errorString) + +# Plot the pure single component isotherm for the n gases +fig = plt.figure +# Histogram for gas 1 +ax1 = plt.subplot(1,1,1) +ax1.hist(resultOutput[:,0], bins = nBins, range = rangeX, density = densityX, + linewidth=1.5, histtype = histTypeX, color='r', alpha = alphaX, label = '$g_1$') +ax1.axvline(x=molFracG1, linewidth=1, linestyle='dotted', color = 'r', alpha = 0.6) + +# Histogram for gas 2 +# ax2 = plt.subplot(1,4,2) +ax1.hist(resultOutput[:,1], bins = nBins, range = rangeX, density = densityX, + linewidth=1.5, histtype = histTypeX, color='b', alpha = alphaX, label = '$g_2$') +ax1.axvline(x=molFracG2, linewidth=1, linestyle='dotted', color = 'b', alpha = 0.6) + +if numberOfGases == 3: + # Histogram for gas 3 + # ax3 = plt.subplot(1,4,3) + ax1.hist(resultOutput[:,2], bins = nBins, range = rangeX, density = densityX, + linewidth=1.5, histtype = histTypeX, color='g', alpha = alphaX, label = '$g_3$') + ax1.axvline(x=molFracG3, linewidth=1, linestyle='dotted', color = 'g', alpha = 0.6) + +ax1.locator_params(axis="x", nbins=4) +ax1.locator_params(axis="y", nbins=4) +ax1.set(xlabel='$y_i$ [-]', + ylabel='$f$ [-]', + xlim = xLimits, ylim = yLimits) +ax1.legend() + +# Save the figure +if saveFlag: + # FileName: ConcEstimate___ + saveFileName = "ConcEstimate_" + currentDT + "_" + gitCommitID + "_" + simID_loadedFile + saveFileExtension + savePath = os.path.join('..','simulationFigures',saveFileName) + # Check if inputResources directory exists or not. If not, create the folder + if not os.path.exists(os.path.join('..','simulationFigures')): + os.mkdir(os.path.join('..','simulationFigures')) + plt.savefig (savePath) + + +# For the figure to be saved show should appear after the save +plt.show() + +# Histogram for the sum of mole fraction +fig = plt.figure +ax2 = plt.subplot(1,1,1) +if numberOfGases == 2: + ax2.hist(np.sum(resultOutput,1), bins = nBins, range = rangeXS, density = densityX, + linewidth=1.5, histtype = histTypeX, color='k', alpha = alphaX) + ax2.axvline(x=1., linewidth=1, linestyle='dotted', color = 'k') +elif numberOfGases == 3: + ax2.hist(np.sum(resultOutput,1), bins = nBins, range = rangeXS, density = densityX, + linewidth=1.5, histtype = histTypeX, color='k', alpha = alphaX) + ax2.axvline(x=1., linewidth=1, linestyle='dotted', color = 'k') + +ax2.set(xlabel='$\Sigma y_i$ [-]', + xlim = xLimitsSum, ylim = yLimitsSum) +ax2.locator_params(axis="x", nbins=4) +ax2.locator_params(axis="y", nbins=4) + +# Save the figure +if saveFlag: + # FileName: ConcEstimate___ + saveFileName = "ConcEstimateSum_" + currentDT + "_" + gitCommitID + "_" + simID_loadedFile + saveFileExtension + savePath = os.path.join('..','simulationFigures',saveFileName) + # Check if inputResources directory exists or not. If not, create the folder + if not os.path.exists(os.path.join('..','simulationFigures')): + os.mkdir(os.path.join('..','simulationFigures')) + plt.savefig (savePath) + +# For the figure to be saved show should appear after the save +plt.show() + +if plotZoomDist: + # Plot the pure single component isotherm for the n gases + fig = plt.figure + # Histogram for gas 1 + ax1 = plt.subplot(1,1,1) + ax1.hist(resultOutput[:,0], bins = nBins, range = rangeX, density = densityX, + linewidth=1.5, histtype = histTypeX, color='r', alpha = alphaX, label = '$g_1$') + ax1.axvline(x=molFracG1, linewidth=1, linestyle='dotted', color = 'r', alpha = 0.6) + + # Histogram for gas 2 + # ax2 = plt.subplot(1,4,2) + ax1.hist(resultOutput[:,1], bins = nBins, range = rangeX, density = densityX, + linewidth=1.5, histtype = histTypeX, color='b', alpha = alphaX, label = '$g_2$') + ax1.axvline(x=molFracG2, linewidth=1, linestyle='dotted', color = 'b', alpha = 0.6) + + if numberOfGases == 3: + # Histogram for gas 3 + # ax3 = plt.subplot(1,4,3) + ax1.hist(resultOutput[:,2], bins = nBins, range = rangeX, density = densityX, + linewidth=1.5, histtype = histTypeX, color='g', alpha = alphaX, label = '$g_3$') + ax1.axvline(x=molFracG3, linewidth=1, linestyle='dotted', color = 'g', alpha = 0.6) + + ax1.locator_params(axis="x", nbins=4) + ax1.locator_params(axis="y", nbins=4) + ax1.set(xlabel='$y_i$ [-]', + ylabel='$f$ [-]', + xlim = xLimitsZ, ylim = yLimitsZ) + ax1.legend() + + # Save the figure + if saveFlag: + # FileName: ConcEstimate___ + saveFileName = "ConcEstimateZoom_" + currentDT + "_" + gitCommitID + "_" + simID_loadedFile + saveFileExtension + savePath = os.path.join('..','simulationFigures',saveFileName) + # Check if inputResources directory exists or not. If not, create the folder + if not os.path.exists(os.path.join('..','simulationFigures')): + os.mkdir(os.path.join('..','simulationFigures')) + plt.savefig (savePath) + + + # For the figure to be saved show should appear after the save + plt.show() + + # Histogram for the sum of mole fraction + fig = plt.figure + ax2 = plt.subplot(1,1,1) + ax2.hist(np.sum(resultOutput,1), bins = nBins, range = rangeXS, density = densityX, + linewidth=1.5, histtype = histTypeX, color='k', alpha = alphaX) + ax2.axvline(x=1., linewidth=1, linestyle='dotted', color = 'k') + ax2.set(xlabel='$\Sigma y_i$ [-]', + xlim = xLimitsSumZ, ylim = yLimitsSumZ) + ax2.locator_params(axis="x", nbins=4) + ax2.locator_params(axis="y", nbins=4) + + # Save the figure + if saveFlag: + # FileName: ConcEstimate___ + saveFileName = "ConcEstimateSumZoom_" + currentDT + "_" + gitCommitID + "_" + simID_loadedFile + saveFileExtension + savePath = os.path.join('..','simulationFigures',saveFileName) + # Check if inputResources directory exists or not. If not, create the folder + if not os.path.exists(os.path.join('..','simulationFigures')): + os.mkdir(os.path.join('..','simulationFigures')) + plt.savefig (savePath) + + # For the figure to be saved show should appear after the save + plt.show() \ No newline at end of file diff --git a/plotFunctions/plotConcentrationFullModel.py b/plotFunctions/plotConcentrationFullModel.py new file mode 100644 index 0000000..b7830f7 --- /dev/null +++ b/plotFunctions/plotConcentrationFullModel.py @@ -0,0 +1,123 @@ +############################################################################ +# +# Imperial College London, United Kingdom +# Multifunctional Nanomaterials Laboratory +# +# Project: ERASE +# Year: 2021 +# Python: Python 3.7 +# Authors: Ashwin Kumar Rajagopalan (AK) +# +# Purpose: +# Plotting function to compare the concentration estimates obtained from using +# the full model and one from assuming instantaneous equilibrium (non kinetics) +# +# Last modified: +# - 2021-01-26, AK: Initial creation +# +# Input arguments: +# +# +# Output arguments: +# +# +############################################################################ + +import numpy as np +from numpy import load +import os +import auxiliaryFunctions +import matplotlib.pyplot as plt +plt.style.use('doubleColumn.mplstyle') # Custom matplotlib style file + +# Get the commit ID of the current repository +gitCommitID = auxiliaryFunctions.getCommitID() + +# Get the current date and time for saving purposes +currentDT = auxiliaryFunctions.getCurrentDateTime() + +# Save file extension (png or pdf) +saveFileExtension = ".png" + +# Save flag +saveFlag = True + +# Color combo +colorForMat = ["ff006e","3a86ff"] +colorForConc = ["fe7f2d","619b8a","233d4d"] + +# Y limits for the plot +X_LIMITS = [0.,1000.] +scaleLog = True + +# Legend First +legendText = ["Without Noise", "With Noise"] +legendFlag = False + +# File name for equilibrium model and full model estimates +loadFileName_E = "fullModelConcentrationEstimate_6-2_20210126_1810_848451b.npz" # Eqbm +loadFileName_F = "fullModelConcentrationEstimate_6-2_20210126_1434_848451b.npz" # Full model + +# Parse out equilbirum results file +simResultsFile = os.path.join('..','simulationResults',loadFileName_E); +loadedFile_E = load(simResultsFile, allow_pickle=True) +concentrationEstimate_E = loadedFile_E["arrayConcentration"] + +# Parse out full model results file +simResultsFile = os.path.join('..','simulationResults',loadFileName_F); +loadedFile_F = load(simResultsFile, allow_pickle=True) +concentrationEstimate_F = loadedFile_F["arrayConcentration"] + +# Parse out true responses (this should be the same for both eqbm and full +# model (here loaded from eqbm) +trueResponseStruct = loadedFile_F["outputStruct"].item() +# Parse out time +timeSim = [] +timeSim = trueResponseStruct[0]["timeSim"] +# Parse out feed mole fraction +feedMoleFrac = trueResponseStruct[0]["inputParameters"][4] +# Parse out true sensor finger print +sensorFingerPrint = np.zeros([len(timeSim),len(trueResponseStruct)]) +for ii in range(len(trueResponseStruct)): + sensorFingerPrint[:,ii] = trueResponseStruct[ii]["sensorFingerPrint"] + +# Points that will be taken for sampling (for plots) +lenSampling = 6 +fig = plt.figure +# Plot the true sensor response (using the full model) +ax = plt.subplot(1,2,1) +ax.plot(timeSim[0:len(timeSim):lenSampling], + sensorFingerPrint[0:len(timeSim):lenSampling,0], + marker = 'o', markersize = 2, linestyle = 'dotted', linewidth = 0.5, + color='#'+colorForMat[0], label = str(trueResponseStruct[0]["inputParameters"][0]).replace('[','').replace(']','')) +ax.plot(timeSim[0:len(timeSim):lenSampling], + sensorFingerPrint[0:len(timeSim):lenSampling,1], + marker = 'o', markersize = 2, linestyle = 'dotted', linewidth = 0.5, + color='#'+colorForMat[1], label = str(trueResponseStruct[1]["inputParameters"][0]).replace('[','').replace(']','')) +ax.locator_params(axis="x", nbins=4) +ax.locator_params(axis="y", nbins=4) +ax.set(xlabel='$t$ [s]', + ylabel='$m_i$ [g kg$^{\mathregular{-1}}$]', + xlim = X_LIMITS, ylim = [0, 30]) +ax.legend() + +# Plot the evolution of the gas composition with respect to time +ax = plt.subplot(1,2,2) +ax.plot([timeSim[0],timeSim[-1]], [feedMoleFrac[0],feedMoleFrac[0]], + linestyle = 'dashed', linewidth = 1., + color='#'+colorForConc[2]) +ax.plot(timeSim[0:len(timeSim):lenSampling], + concentrationEstimate_E[0:len(timeSim):lenSampling,2], + marker = 'v', markersize = 2, linestyle = 'dotted', linewidth = 0.5, + color='#'+colorForConc[0], label = 'Equilibrium') +ax.plot(timeSim[0:len(timeSim):lenSampling], + concentrationEstimate_F[0:len(timeSim):lenSampling,2], + marker = 'o', markersize = 2, linestyle = 'dotted', linewidth = 0.5, + color='#'+colorForConc[1], label = 'Full Model') +ax.locator_params(axis="x", nbins=4) +ax.locator_params(axis="y", nbins=4) +ax.set(xlabel='$t$ [s]', + ylabel='$y_1$ [-]', + xlim = X_LIMITS, ylim = [0, 0.2]) +ax.legend() +plt.show() \ No newline at end of file diff --git a/plotFunctions/plotConcentrationViolinPlots.py b/plotFunctions/plotConcentrationViolinPlots.py new file mode 100755 index 0000000..65d496c --- /dev/null +++ b/plotFunctions/plotConcentrationViolinPlots.py @@ -0,0 +1,326 @@ +############################################################################ +# +# Imperial College London, United Kingdom +# Multifunctional Nanomaterials Laboratory +# +# Project: ERASE +# Year: 2020 +# Python: Python 3.7 +# Authors: Ashwin Kumar Rajagopalan (AK) +# +# Purpose: +# Plots to visualize different sensor responses +# +# Last modified: +# - 2021-01-11, AK: Cosmetic changes +# - 2020-11-27, AK: Major modification for 3 gas system +# - 2020-11-23, AK: Add standard deviation/CV plotting +# - 2020-11-20, AK: Introduce 3 gas capability +# - 2020-11-18, AK: Changes to data reconciliation and new plots +# - 2020-11-13, AK: Initial creation +# +# Input arguments: +# +# +# Output arguments: +# +# +############################################################################ +import pdb +import numpy as np +from numpy import load +import os +import pandas as pd +import seaborn as sns +import auxiliaryFunctions +import matplotlib.pyplot as plt +plt.style.use('doubleColumn.mplstyle') # Custom matplotlib style file + +# Get the commit ID of the current repository +gitCommitID = auxiliaryFunctions.getCommitID() + +# Get the current date and time for saving purposes +currentDT = auxiliaryFunctions.getCurrentDateTime() + +# Save file extension (png or pdf) +saveFileExtension = ".png" + +# Save flag +saveFlag = False + +# Flag for comparison +flagComparison = False + +# Color combo +colorTemp = ["eac435","345995","03cea4","fb4d3d","ca1551"] +colorForPlot = ["#" + counter for counter in colorTemp] + +# Mole fraction ID +moleFracID = 3 +meanMolFrac = [0.00,0.2,0.35,0.5,0.75,0.9] + +# Y limits for the plot +Y_LIMITS = [None,None] +scaleLog = True + +# Legend First +legendText = ["Without Noise", "With Noise"] +legendFlag = False + +# Sensor ID +sensorText = ["y1", "y2", "y3", "0/1/2"] + +# Initialize x, y, and type for the plotting +concatenatedX = [] +concatenatedX2 = [] +concatenatedX3 = [] +concatenatedY1 = [] +concatenatedY2 = [] +concatenatedY3 = [] +concatenatedType = [] + +# File to be loaded for the left of violin plot +loadFileName = ["sensitivityAnalysis_6-12-55_20210112_1859_50a3ed7.npz", + "sensitivityAnalysis_6-12-55_20210112_1924_50a3ed7.npz", + "sensitivityAnalysis_6-12-55_20210112_1939_50a3ed7.npz"] + +saveFileSensorText = [3,4,1] + +if flagComparison and len(loadFileName) != 2: + errorString = "When flagComparison is True, only two files can be loaded for comparison." + raise Exception(errorString) + +# Loop through the different files to generate the violin plot +for kk in range(len(loadFileName)): + # Initialize x, y, and type for the local loop + xVar = [] + x2Var = [] + x3Var = [] + y1Var = [] + y2Var = [] + y3Var = [] + typeVar = [] + + simResultsFile = os.path.join('..','simulationResults',loadFileName[kk]); + resultOutput = load(simResultsFile)["arrayConcentration"] + numberOfGases = load(simResultsFile)["numberOfGases"] + if numberOfGases == 2: + moleFrac = load(simResultsFile)["trueMoleFrac"] + elif numberOfGases == 3: + moleFracTemp = load(simResultsFile)["trueMoleFrac"] + moleFrac = moleFracTemp[:,0] + moleFrac2 = moleFracTemp[:,1] + moleFrac3 = moleFracTemp[:,2] + + # Loop through all the molefractions + for ii in range(resultOutput.shape[0]): + # For the cases where there are two gases + if numberOfGases == 2: + if resultOutput.shape[2] == 4: + counterInd = 0 + elif resultOutput.shape[2] == 5: + counterInd = 1 + # For the cases where there are two gases + elif numberOfGases == 3: + if resultOutput.shape[2] == 6: + counterInd = 1 + + y1Var = np.concatenate((y1Var,resultOutput[ii,:,counterInd+2])) # y1 + y2Var = np.concatenate((y2Var,resultOutput[ii,:,counterInd+3])) # y2 + if numberOfGases == 3: + y3Var = np.concatenate((y3Var,resultOutput[ii,:,counterInd+4])) # y3 + xVar = xVar + ([str(moleFrac[ii])] * len(resultOutput[ii,:,counterInd+2])) # x (true mole fraction) + if numberOfGases == 3: + x2Var = x2Var + ([str(moleFrac2[ii])] * len(resultOutput[ii,:,counterInd+2])) # x (true mole fraction) + x3Var = x3Var + ([str(moleFrac3[ii])] * len(resultOutput[ii,:,counterInd+2])) # x (true mole fraction) + if not flagComparison: + typeVar = typeVar+[sensorText[kk]] * len(resultOutput[ii,:,counterInd+2]) + # Generate type for comparison + if flagComparison: + typeVar = [legendText[kk]] * len(y1Var) # Type - string + + # Concatenate all the data to form a data frame with x, y, and type + concatenatedX = concatenatedX + xVar + concatenatedY1 = np.concatenate((concatenatedY1,y1Var)) + concatenatedY2 = np.concatenate((concatenatedY2,y2Var)) + if numberOfGases == 3: + concatenatedY3 = np.concatenate((concatenatedY3,y3Var)) + concatenatedX2 = concatenatedX2 + x2Var + concatenatedX3 = concatenatedX3 + x3Var + + concatenatedType = concatenatedType + typeVar + + # Reinitialize all the loaded values to empty variable + simResultsFile = [] + resultOutput = [] + moleFrac = [] + +# Generate panda data frame +# x = molefraction (true) +# y = molefraction (estimated) +# dataType = either sensor id/comparison type +if numberOfGases == 2: + df = pd.DataFrame({'x':concatenatedX, + 'y1':concatenatedY1, + 'y2':concatenatedY2, + 'dataType':concatenatedType}) + # Compute the mean and standard deviation + meanData = df.groupby(['dataType','x'], as_index=False, sort=False).mean() + stdData = df.groupby(['dataType','x'], as_index=False, sort=False).std() + # Coefficient of variation + cvData = stdData.copy() + cvData['y1'] = stdData['y1']/meanData['y1'] + +elif numberOfGases == 3: + df = pd.DataFrame({'x':concatenatedX, + 'x2':concatenatedX2, + 'x3':concatenatedX3, + 'y1':concatenatedY1, + 'y2':concatenatedY2, + 'y3':concatenatedY3, + 'dataType':concatenatedType}) + # Compute the mean and standard deviation + meanData = df.groupby(['dataType','x','x2','x3'], as_index=False, sort=False).mean() + stdData = df.groupby(['dataType','x','x2','x3'], as_index=False, sort=False).std() + # Coefficient of variation + cvData = stdData.copy() + cvData['y1'] = stdData['y1']/meanData['y1'] + cvData['y2'] = stdData['y2']/meanData['y2'] + cvData['y3'] = stdData['y3']/meanData['y3'] + +# Plot quantities of interest +# Two gases +if numberOfGases == 2: + # Plot the figure + sns.set(style="ticks", palette="pastel", color_codes=True) + fig = plt.figure + ax1 = plt.subplot(1,1,1) + # Draw a nested violinplot for easier comparison + if flagComparison: + if scaleLog: + ax1.set_yscale('log') + sns.violinplot(data=df, x="x", y="y1", hue="dataType", inner = "box", + split=True, linewidth=1, palette={legendText[0]: colorForPlot[0], + legendText[1]: colorForPlot[1]}, + scale='width') + ax1.set(xlabel='$y_1$ [-]', ylabel='${\hat{y}_1}$ [-]', ylim = Y_LIMITS) + plt.legend(loc='upper left') + if not legendFlag: + plt.legend([],[], frameon=False) + # Draw violin plot for compaison of different sensors + else: + sns.violinplot(data=df[df.x == str(meanMolFrac[moleFracID])], + x="dataType", y="y1", inner = "box", linewidth=1, + scale='width', palette = colorForPlot[0:len(loadFileName)]) + ax1.set(xlabel='Sensor ID [-]', ylabel='${\hat{y}_1}$ [-]', ylim = Y_LIMITS) + if flagComparison: + for kk in range(len(meanMolFrac)): + ax1.axhline(meanMolFrac[kk], linestyle=':', linewidth=1, color = '#c0c0c0') + else: + ax1.axhline(meanMolFrac[moleFracID], linestyle=':', linewidth=1, color = '#c0c0c0') + ax1.locator_params(axis="y", nbins=4) + # Save the figure + if saveFlag: + # FileName: SensorViolinPlot____> + saveFileSensorText = str(saveFileSensorText).replace('[','').replace(']','').replace(' ','-').replace(',','') + saveFileName = "SensorViolinPlot_" + saveFileSensorText + "_" + currentDT + "_" + gitCommitID + saveFileExtension + savePath = os.path.join('..','simulationFigures',saveFileName) + # Check if inputResources directory exists or not. If not, create the folder + if not os.path.exists(os.path.join('..','simulationFigures')): + os.mkdir(os.path.join('..','simulationFigures')) + plt.savefig (savePath) + plt.show() + + plt.style.use('doubleColumn.mplstyle') # Custom matplotlib style file + fig = plt.figure + # Standard deviation + ax1 = plt.subplot(1,2,1) + stdData["x"] = pd.to_numeric(stdData["x"], downcast="float") + sns.lineplot(data=stdData, x='x', y='y1', hue='dataType', style='dataType', + dashes = False, markers = ['o']*len(loadFileName), + palette = colorForPlot[0:len(loadFileName)]) + ax1.set(xlabel='$y_1$ [-]', + ylabel='$\sigma (\hat{y}_1)$ [-]', + xlim = [0.,1.], ylim = [1e-10,1.]) + ax1.set_yscale('log') + ax1.locator_params(axis="x", nbins=4) + if len(loadFileName) > 1: + plt.legend(loc='best') + else: + plt.legend([],[], frameon=False) + + # CV + ax2 = plt.subplot(1,2,2) + cvData["x"] = pd.to_numeric(cvData["x"], downcast="float") + sns.lineplot(data=cvData, x='x', y='y1', hue='dataType', style='dataType', + dashes = False, markers = ['o']*len(loadFileName), + palette = colorForPlot[0:len(loadFileName)]) + ax2.set(xlabel='$y_1$ [-]', + ylabel='$CV (\hat{y}_1)$ [-]', + xlim = [0.,1.], ylim = [1e-10,None]) + ax2.locator_params(axis="x", nbins=4) + ax2.set_yscale('log') + if not legendFlag: + plt.legend([],[], frameon=False) + if saveFlag: + # FileName: SensorViolinPlot____> + saveFileSensorText = str(saveFileSensorText).replace('[','').replace(']','').replace(' ','-').replace(',','') + saveFileName = "SensorStdCV_" + saveFileSensorText + "_" + currentDT + "_" + gitCommitID + saveFileExtension + savePath = os.path.join('..','simulationFigures',saveFileName) + # Check if inputResources directory exists or not. If not, create the folder + if not os.path.exists(os.path.join('..','simulationFigures')): + os.mkdir(os.path.join('..','simulationFigures')) + plt.savefig (savePath) + plt.show() + +# Three gases +if numberOfGases == 3: + plt.style.use('doubleColumn.mplstyle') # Custom matplotlib style file + fig = plt.figure + ax1 = plt.subplot(1,3,1) + cvData["x"] = pd.to_numeric(cvData["x"], downcast="float") + sns.lineplot(data=cvData, x='x', y='y1', hue='dataType', style='dataType', + dashes = False, markers = ['o']*len(loadFileName), + palette = colorForPlot[0:len(loadFileName)]) + ax1.set(xlabel='$y_1$ [-]', + ylabel='$CV (\hat{y}_1)$ [-]', + xlim = [0.,1.], ylim = [1e-8,100.]) + ax1.locator_params(axis="x", nbins=4) + ax1.set_yscale('log') + if not legendFlag: + plt.legend([],[], frameon=False) + ax2 = plt.subplot(1,3,2) + cvData["x2"] = pd.to_numeric(cvData["x2"], downcast="float") + sns.lineplot(data=cvData, x='x2', y='y2', hue='dataType', style='dataType', + dashes = False, markers = ['o']*len(loadFileName), + palette = colorForPlot[0:len(loadFileName)]) + ax2.set(xlabel='$y_2$ [-]', + ylabel='$CV (\hat{y}_2)$ [-]', + xlim = [0.,1.], ylim = [1e-8,100.]) + ax2.locator_params(axis="x", nbins=4) + ax2.set_yscale('log') + if not legendFlag: + plt.legend([],[], frameon=False) + ax3 = plt.subplot(1,3,3) + cvData["x3"] = pd.to_numeric(cvData["x3"], downcast="float") + sns.lineplot(data=cvData, x='x3', y='y3', hue='dataType', style='dataType', + dashes = False, markers = ['o']*len(loadFileName), + palette = colorForPlot[0:len(loadFileName)]) + ax3.set(xlabel='$y_3$ [-]', + ylabel='$CV (\hat{y}_3)$ [-]', + xlim = [-0.1,1], ylim = [1e-8,100.]) + ax3.locator_params(axis="x", nbins=4) + ax3.set_yscale('log') + if not legendFlag: + plt.legend([],[], frameon=False) + + if saveFlag: + # FileName: SensorViolinPlot____> + saveFileSensorText = str(saveFileSensorText).replace('[','').replace(']','').replace(' ','-').replace(',','') + saveFileName = "SensorStdCV_" + saveFileSensorText + "_" + currentDT + "_" + gitCommitID + saveFileExtension + savePath = os.path.join('..','simulationFigures',saveFileName) + # Check if inputResources directory exists or not. If not, create the folder + if not os.path.exists(os.path.join('..','simulationFigures')): + os.mkdir(os.path.join('..','simulationFigures')) + plt.savefig (savePath) + plt.show() \ No newline at end of file diff --git a/plotFunctions/plotConstrainedObjective.py b/plotFunctions/plotConstrainedObjective.py new file mode 100644 index 0000000..2a7a284 --- /dev/null +++ b/plotFunctions/plotConstrainedObjective.py @@ -0,0 +1,245 @@ +############################################################################ +# +# Imperial College London, United Kingdom +# Multifunctional Nanomaterials Laboratory +# +# Project: ERASE +# Year: 2020 +# Python: Python 3.7 +# Authors: Ashwin Kumar Rajagopalan (AK) +# +# Purpose: +# Plots the objective function used for concentration estimation for a mole +# fraction grid. The constraint of mass conservation will show up as a line +# in the plot. The concept\visualization is similar to the one of Lagrange +# multipliers. +# +# Last modified: +# - 2021-01-08, AK: Code improvements +# - 2020-12-18, AK: Initial creation +# +# Input arguments: +# +# +# Output arguments: +# +# +############################################################################ +def plotConstrainedObjective(): + import numpy as np + import os + import matplotlib.pyplot as plt + from mpl_toolkits.mplot3d import Axes3D + plt.style.use('doubleColumn.mplstyle') # Custom matplotlib style file + import auxiliaryFunctions + + os.chdir("..") + + # Save flag for figure + saveFlag = False + + # Number of molefractions + numMolFrac= 1001 + numMesh = 21 + + # Total pressure of the gas [Pa] + pressureTotal = np.array([1.e5]); + + # Temperature of the gas [K] + # Can be a vector of temperatures + temperature = np.array([298.15]); + + # Number of Adsorbents + numberOfAdsorbents = 20 + + # Number of Gases + numberOfGases = 3 + + # Mole Fraction of interest + moleFrac = np.array([[0.2,0.8,0.0], [0.6,0.4,0.0]]) + + + # Multiplier Error + multiplierError = [1,1,1] + + # Sensor ID + sensorID = np.array([3,4,1]) + + # Get the commit ID of the current repository + gitCommitID = auxiliaryFunctions.getCommitID() + + # Get the current date and time for saving purposes + currentDT = auxiliaryFunctions.getCurrentDateTime() + + # For two gases + if numberOfGases == 2: + # Obtain the objective functino for a mole fraction grid (constraint + # shown as a line int the plot) - similar to Lagrange multiplier + objFunction = np.zeros([moleFrac.shape[0],numMesh,numMesh]) + for ii in range(moleFrac.shape[0]): + moleFracTemp = moleFrac[ii,:] + x1m,x2m,objFunctionTemp = functionForMoleFrac(numberOfAdsorbents,numberOfGases,pressureTotal, + temperature,sensorID,moleFracTemp,numMolFrac, + numMesh,multiplierError) + objFunction[ii,:,:] = objFunctionTemp[:,:] + if ii == 0: + minObj = np.max([-np.inf,np.min(objFunction[ii,:,:])]) + maxObj = np.min([np.inf,np.max(objFunction[ii,:,:])]) + else: + minObj = np.max([minObj,np.min(objFunction[ii,:,:])]) + maxObj = np.min([maxObj,np.max(objFunction[ii,:,:])]) + + + fig = plt.figure + ax = plt.subplot(1,2,1) + # Plot obbjective function for first mole fraction + plt.contourf(x1m,x2m,objFunction[0,:,:],levels = np.linspace(minObj,maxObj,200),cmap='RdYlBu') + plt.colorbar() + plt.plot([0.0,1.0],[1.0,0.0],linestyle = ':', linewidth = 1, color = 'k') + ax.set(xlabel='$y_1$ [-]', + ylabel='$y_2$ [-]', + xlim = [0,1], ylim = [0, 1.]) + ax.locator_params(axis="x", nbins=4) + ax.locator_params(axis="y", nbins=4) + + ax = plt.subplot(1,2,2) + # Plot obbjective function for second mole fraction + plt.contourf(x1m,x2m,objFunction[1,:,:],levels = np.linspace(minObj,maxObj,200),cmap='RdYlBu') + plt.colorbar() + plt.plot([0.0,1.0],[1.0,0.0],linestyle = ':', linewidth = 1, color = 'k') + ax.set(xlabel='$y_1$ [-]', + ylabel='$y_2$ [-]', + xlim = [0,1], ylim = [0, 1.]) + ax.locator_params(axis="x", nbins=4) + ax.locator_params(axis="y", nbins=4) + + # For three gases + elif numberOfGases == 3: + # Obtain the objective functino for a mole fraction grid (constraint + # shown as a line int the plot) - similar to Lagrange multiplier + objFunction = np.zeros([moleFrac.shape[0],numMesh,numMesh]) + for ii in range(moleFrac.shape[0]): + moleFracTemp = moleFrac[ii,:] + x1m,x2m,x3m,objFunctionTemp = functionForMoleFrac(numberOfAdsorbents,numberOfGases,pressureTotal, + temperature,sensorID,moleFracTemp,numMolFrac, + numMesh,multiplierError) + plotInd = np.where(x1m[:,0]==moleFracTemp[2])[0][0] + objFunction[ii,:,:] = objFunctionTemp[:,:,plotInd] + if ii == 0: + minObj = np.max([-np.inf,np.min(objFunction[ii,:,:])]) + maxObj = np.min([np.inf,np.max(objFunction[ii,:,:])]) + else: + minObj = np.max([minObj,np.min(objFunction[ii,:,:])]) + maxObj = np.min([maxObj,np.max(objFunction[ii,:,:])]) + + # Plot obbjective function stacked on top of each other + fig = plt.figure() + ax = fig.gca(projection='3d') + for ii in range(moleFrac.shape[0]): + ax.contourf(x1m[:,:,plotInd],x2m[:,:,plotInd], + objFunction[ii,:,:],zdir='z',offset = ii, + levels = np.linspace(minObj,maxObj,30), + cmap='RdYlBu') + + # Set 3D-axis-limits: + ax.set_xlim3d(0,1) + ax.set_ylim3d(0,1) + ax.set_zlim3d(0,moleFrac.shape[0]) + plt.show() + + fig = plt.figure + ax = plt.subplot(1,2,1) + # Plot obbjective function for first mole fraction + plt.contourf(x1m[:,:,0],x2m[:,:,0],objFunction[0,:,:],levels = np.linspace(minObj,maxObj,200),cmap='RdYlBu') + plt.colorbar() + plt.plot([0.0,1-moleFrac[0,2]],[1-moleFrac[0,2],0.0],linestyle = ':', linewidth = 1, color = 'k') + ax.set(xlabel='$y_1$ [-]', + ylabel='$y_2$ [-]', + xlim = [0,1-moleFrac[0,2]], ylim = [0., 1.-moleFrac[0,2]]) + ax.locator_params(axis="x", nbins=4) + ax.locator_params(axis="y", nbins=4) + + ax = plt.subplot(1,2,2) + # Plot obbjective function for second mole fraction + plt.contourf(x1m[:,:,0],x2m[:,:,0],objFunction[1,:,:],levels = np.linspace(minObj,maxObj,200),cmap='RdYlBu') + plt.colorbar() + plt.plot([0.0,1-moleFrac[1,2]],[1-moleFrac[1,2],0.0],linestyle = ':', linewidth = 1, color = 'k') + ax.set(xlabel='$y_1$ [-]', + ylabel='$y_2$ [-]', + xlim = [0,1-moleFrac[1,2]], ylim = [0, 1-moleFrac[1,2]]) + ax.locator_params(axis="x", nbins=4) + ax.locator_params(axis="y", nbins=4) + +# Function to evaluate the objective function for a given mole fraction +def functionForMoleFrac(numberOfAdsorbents,numberOfGases,pressureTotal, + temperature,sensorID,moleFrac,numMolFrac, + numMesh,multiplierError): + import numpy as np + from generateTrueSensorResponse import generateTrueSensorResponse + from simulateSensorArray import simulateSensorArray + # Simulate the sensor response for all possible concentrations + if numberOfGases == 2: + x1 = np.linspace(0,1,numMesh) + x2 = np.linspace(0,1,numMesh) + x1m, x2m = np.meshgrid(x1, x2, sparse=False, indexing='ij') + arraySimResponse = np.zeros([len(sensorID),numMesh,numMesh]) + elif numberOfGases == 3: + x1 = np.linspace(0,1,numMesh) + x2 = np.linspace(0,1,numMesh) + x3 = np.linspace(0,1,numMesh) + x1m, x2m, x3m = np.meshgrid(x1, x2, x3, sparse=False, indexing='ij') + arraySimResponse = np.zeros([len(sensorID),numMesh,numMesh,numMesh]) + + # Get simulated sensor response + if numberOfGases == 2: + for ii in range(numMesh): + for jj in range(numMesh): + arraySimResponseTemp = simulateSensorArray(sensorID, pressureTotal, + temperature, np.array([[x1m[ii,jj],x2m[ii,jj]]])) * multiplierError + for kk in range(len(sensorID)): + arraySimResponse[kk,ii,jj] = arraySimResponseTemp[kk] + elif numberOfGases == 3: + for ii in range(numMesh): + for jj in range(numMesh): + for kk in range(numMesh): + arraySimResponseTemp = simulateSensorArray(sensorID, pressureTotal, + temperature, np.array([[x1m[ii,jj,kk],x2m[ii,jj,kk],x3m[ii,jj,kk]]])) * multiplierError + for ll in range(len(sensorID)): + arraySimResponse[ll,ii,jj,kk] = arraySimResponseTemp[ll] + + + # Get the individual sensor reponse for all the given "experimental/test" concentrations + sensorTrueResponse = generateTrueSensorResponse(numberOfAdsorbents,numberOfGases, + pressureTotal,temperature,moleFraction = moleFrac) + # Parse out the true sensor response for the desired sensors in the array + arrayTrueResponseTemp = np.zeros(len(sensorID)) + for ii in range(len(sensorID)): + arrayTrueResponseTemp[ii] = sensorTrueResponse[sensorID[ii]]*multiplierError[ii] + if numberOfGases == 2: + arrayTrueResponse = np.zeros([len(sensorID),numMesh,numMesh]) + for kk in range(len(sensorID)): + arrayTrueResponse[kk,:,:] = np.tile(arrayTrueResponseTemp[kk],(numMesh,numMesh)) + elif numberOfGases == 3: + arrayTrueResponse = np.zeros([len(sensorID),numMesh,numMesh,numMesh]) + for ll in range(len(sensorID)): + arrayTrueResponse[ll,:,:,:] = np.tile(arrayTrueResponseTemp[ll],(numMesh,numMesh,numMesh)) + + # Compute the objective function over all the mole fractions + if numberOfGases == 2: + objFunction = np.zeros([numMesh,numMesh]) + for kk in range(len(sensorID)): + objFunctionTemp = np.power(np.divide(arrayTrueResponse[kk,:,:] - arraySimResponse[kk,:,:],arrayTrueResponse[kk,:,:]),2) + objFunction += objFunctionTemp + objFunctionTemp = [] + elif numberOfGases == 3: + objFunction = np.zeros([numMesh,numMesh,numMesh]) + for kk in range(len(sensorID)): + objFunctionTemp = np.power(np.divide(arrayTrueResponse[kk,:,:,:] - arraySimResponse[kk,:,:,:],arrayTrueResponse[kk,:,:,:]),2) + objFunction += objFunctionTemp + objFunctionTemp = [] + + # Function return + if numberOfGases == 2: + return x1m,x2m,objFunction + elif numberOfGases == 3: + return x1m,x2m,x3m,objFunction \ No newline at end of file diff --git a/plotFunctions/plotObjectiveFunction.py b/plotFunctions/plotObjectiveFunction.py new file mode 100755 index 0000000..4145a2d --- /dev/null +++ b/plotFunctions/plotObjectiveFunction.py @@ -0,0 +1,387 @@ +############################################################################ +# +# Imperial College London, United Kingdom +# Multifunctional Nanomaterials Laboratory +# +# Project: ERASE +# Year: 2020 +# Python: Python 3.7 +# Authors: Ashwin Kumar Rajagopalan (AK) +# +# Purpose: +# Plots the objective function used for concentration estimation +# +# Last modified: +# - 2021-01-11, AK: Cosmetic changes +# - 2020-11-27, AK: More plotting fix for 3 gas system +# - 2020-11-24, AK: Fix for 3 gas system +# - 2020-11-23, AK: Change ternary plots +# - 2020-11-20, AK: Introduce ternary plots +# - 2020-11-19, AK: Add 3 gas knee calculator +# - 2020-11-19, AK: Multigas plotting capability +# - 2020-11-17, AK: Multisensor plotting capability +# - 2020-11-11, AK: Cosmetic changes and add standard deviation plot +# - 2020-11-05, AK: Initial creation +# +# Input arguments: +# +# +# Output arguments: +# +# +############################################################################ +import pdb +import numpy as np +from numpy import load +from kneed import KneeLocator # To compute the knee/elbow of a curve +from generateTrueSensorResponse import generateTrueSensorResponse +from simulateSensorArray import simulateSensorArray +import os +from sklearn.cluster import KMeans +import pandas as pd +import ternary +import scipy +import matplotlib as mpl +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D +plt.style.use('doubleColumn.mplstyle') # Custom matplotlib style file +import auxiliaryFunctions + +os.chdir("..") + +# Save flag for figure +saveFlag = False + +# Save file extension (png or pdf) +saveFileExtension = ".png" + +# Plotting colors +colorsForPlot = ["ff499e","d264b6","a480cf","779be7","49b6ff"] +colorsForPlot = ["eac435","345995","03cea4","fb4d3d","ca1551"] +colorGroup = ["#f94144","#43aa8b"] +colorIntersection = ["ff595e","ffca3a","8ac926","1982c4","6a4c93"] + +# Number of molefractions +numMolFrac= 10001 + +# Total pressure of the gas [Pa] +pressureTotal = np.array([1.e5]); + +# Temperature of the gas [K] +# Can be a vector of temperatures +temperature = np.array([298.15]); + +# Number of Adsorbents +numberOfAdsorbents = 60 + +# Number of Gases +numberOfGases = 3 + +# Experimental mole fraction for 3 gases +if numberOfGases == 3: + inputMoleFracALL = np.array([[0.01, 0.99, 0.00], + [0.10, 0.90, 0.00], + [0.25, 0.75, 0.00], + [0.50, 0.50, 0.00], + [0.75, 0.25, 0.00], + [0.90, 0.10, 0.00], + [0.99, 0.01, 0.00]]) + + # Fix one gas + fixOneGas = False + # Third gas mole fraction + thirdGasMoleFrac = 0. + +# Mole Fraction of interest +moleFrac = [0.3, 0.7] + +# Multiplier Error +multiplierError = [1, 1.,1.] + +# Sensor ID +sensorID = np.array([6,12,55]) +saveFileSensorText = [6,12,55] + +# Acceptable SNR +signalToNoise = 25*0.1 + +# Get the commit ID of the current repository +gitCommitID = auxiliaryFunctions.getCommitID() + +# Get the current date and time for saving purposes +currentDT = auxiliaryFunctions.getCurrentDateTime() + +# Simulate the sensor response for all possible concentrations +if numberOfGases == 2: + moleFractionRange = np.array([np.linspace(0,1,numMolFrac), 1 - np.linspace(0,1,numMolFrac)]).T +elif numberOfGases == 3: + moleFractionRangeTemp = np.random.dirichlet((1,1,1),numMolFrac) + moleFractionRange = moleFractionRangeTemp[moleFractionRangeTemp[:,0].argsort()] + if fixOneGas: + remainingMoleFrac = 1. - thirdGasMoleFrac + moleFractionRange = np.array([np.linspace(0,remainingMoleFrac,numMolFrac), + remainingMoleFrac - np.linspace(0,remainingMoleFrac,numMolFrac), + np.tile(thirdGasMoleFrac,numMolFrac)]).T + +arraySimResponse = np.zeros([moleFractionRange.shape[0],sensorID.shape[0]]) +for ii in range(moleFractionRange.shape[0]): + arraySimResponse[ii,:] = simulateSensorArray(sensorID, pressureTotal, + temperature, np.array([moleFractionRange[ii,:]])) * multiplierError + +# Get the individual sensor reponse for all the given "experimental/test" concentrations +sensorTrueResponse = generateTrueSensorResponse(numberOfAdsorbents,numberOfGases, + pressureTotal,temperature,moleFraction = moleFrac) +# Parse out the true sensor response for the desired sensors in the array +arrayTrueResponse = np.zeros(sensorID.shape[0]) +for ii in range(sensorID.shape[0]): + arrayTrueResponse[ii] = sensorTrueResponse[sensorID[ii]]*multiplierError[ii] +arrayTrueResponse = np.tile(arrayTrueResponse,(moleFractionRange.shape[0],1)) + +# Compute the objective function over all the mole fractions +objFunction = np.sum(np.power((arrayTrueResponse - arraySimResponse)/arrayTrueResponse,2),1) + +# Compute the first derivative, elbow point, and the fill regions for all +# sensors for 2 gases +if numberOfGases == 2 or (numberOfGases == 3 and fixOneGas == True): + xFill = np.zeros([arraySimResponse.shape[1],2]) + # Loop through all sensors + firstDerivative = np.zeros([arraySimResponse.shape[0],arraySimResponse.shape[1]]) + firstDerivativeSimResponse_y1 = np.zeros([moleFractionRange.shape[0],arraySimResponse.shape[1]]) + firstDerivativeSimResponse_y2 = np.zeros([moleFractionRange.shape[0],arraySimResponse.shape[1]]) + secondDerivative = np.zeros([firstDerivative.shape[0],firstDerivative.shape[1]]) + for kk in range(arraySimResponse.shape[1]): + firstDerivative[:,kk] = np.gradient(arraySimResponse[:,kk],moleFractionRange[1,0]-moleFractionRange[0,0]) + secondDerivative[:,kk] = np.gradient(firstDerivative[:,kk],moleFractionRange[1,0]-moleFractionRange[0,0]) + firstDerivativeSimResponse_y1[:,kk] = np.gradient(moleFractionRange[:,0],arraySimResponse[:,kk]) + firstDerivativeSimResponse_y2[:,kk] = np.gradient(moleFractionRange[:,1],arraySimResponse[:,kk]) + # Get the sign of the first derivative for increasing/decreasing + if all(i >= 0. for i in firstDerivative[:,kk]): + slopeDir = "increasing" + elif all(i < 0. for i in firstDerivative[:,kk]): + slopeDir = "decreasing" + else: + print("Dangerous! I should not be here!!!") + # Get the sign of the second derivative for concavity/convexity + if all(i >= 0. for i in secondDerivative[:,kk]): + secondDerDir = "convex" + elif all(i < 0. for i in secondDerivative[:,kk]): + secondDerDir = "concave" + else: + print("Dangerous! I should not be here!!!") + + + kneedle = KneeLocator(moleFractionRange[:,0], arraySimResponse[:,kk], + curve=secondDerDir, direction=slopeDir) + elbowPoint = list(kneedle.all_elbows) + + # Obtain coordinates to fill working region + if secondDerDir == "concave": + if slopeDir == "increasing": + xFill[kk,:] = [0,elbowPoint[0]] + else: + if numberOfGases == 2: + xFill[kk,:] = [elbowPoint[0], 1.0] + elif numberOfGases == 3: + if fixOneGas: + xFill[kk,:] = [elbowPoint[0], 1.-thirdGasMoleFrac] + elif secondDerDir == "convex": + if slopeDir == "increasing": + if numberOfGases == 3: + if fixOneGas: + xFill[kk,:] = [elbowPoint[0],1.-thirdGasMoleFrac] + else: + xFill[kk,:] = [0,elbowPoint[0]] + else: + print("Dangerous! I should not be here!!!") + +# Start plotting +if numberOfGases == 2 or (numberOfGases == 3 and fixOneGas == True): + if numberOfGases == 2: + fig = plt.figure + ax = plt.gca() + elif numberOfGases == 3: + fig = plt.figure + ax = plt.subplot(1,2,1) + ax2 = plt.subplot(1,2,2) + # Loop through all sensors + for kk in range(arraySimResponse.shape[1]): + ax.plot(moleFractionRange[:,0],arraySimResponse[:,kk],color='#'+colorsForPlot[kk], label = '$s_'+str(kk+1)+'$') # Simulated Response + ax.fill_between(xFill[kk,:],1.1*np.max(arraySimResponse), facecolor='#'+colorsForPlot[kk], alpha=0.25) + if numberOfGases == 2: + ax.fill_between([0.,1.],signalToNoise, facecolor='#4a5759', alpha=0.25) + elif numberOfGases == 3: + mpl.rcParams['hatch.linewidth'] = 0.1 + ax.fill_between([0.,1.-thirdGasMoleFrac],signalToNoise, facecolor='#4a5759', alpha=0.25) + ax.fill_between([1.-thirdGasMoleFrac,1.],1.1*np.max(arraySimResponse), facecolor='#555b6e', alpha=0.25, hatch = 'x') + ax2.plot(moleFractionRange[:,1],arraySimResponse[:,kk],color='#'+colorsForPlot[kk], label = '$s_'+str(kk+1)+'$') # Simulated Response + + ax.set(xlabel='$y_1$ [-]', + ylabel='$m_i$ [g kg$^{-1}$]', + xlim = [0,1], ylim = [0, 1.1*np.max(arraySimResponse)]) + ax.locator_params(axis="x", nbins=4) + ax.locator_params(axis="y", nbins=4) + ax.legend() + if numberOfGases == 3: + ax2.set(xlabel='$y_2$ [-]', + ylabel='$m_i$ [g kg$^{-1}$]', + xlim = [0,1], ylim = [0, 1.1*np.max(arraySimResponse)]) + ax2.locator_params(axis="x", nbins=4) + ax2.locator_params(axis="y", nbins=4) + ax2.legend() + + # Save the figure + if saveFlag: + # FileName: SensorResponse___ + sensorText = str(sensorID).replace('[','').replace(']','').replace(' ','-') + saveFileName = "SensorResponse_" + sensorText + "_" + currentDT + "_" + gitCommitID + saveFileExtension + savePath = os.path.join('simulationFigures',saveFileName) + # Check if inputResources directory exists or not. If not, create the folder + if not os.path.exists(os.path.join('..','simulationFigures')): + os.mkdir(os.path.join('..','simulationFigures')) + plt.savefig (savePath) + plt.show() + +# Make ternanry/ternary equivalent plots +if numberOfGases == 3 and fixOneGas == False: + # Loop through all the materials in the array + sensitiveGroup = {} + sensitiveResponse = {} + skewnessResponse = np.zeros(len(sensorID)) + for ii in range(len(sensorID)): + # Key for dictionary entry + tempDictKey = 's_{}'.format(ii) + # Reshape the response for k-means clustering + reshapedArraySimResponse = np.reshape(arraySimResponse[:,ii],[-1,1]) + # Get the skewness of the sensor resposne + skewnessResponse = scipy.stats.skew(reshapedArraySimResponse) + # Perform k-means clustering to separate out the data + kMeansStruct = KMeans(n_clusters=2,random_state=None).fit(reshapedArraySimResponse) + # Obtain the group of the sensor (sensitive/non sensitive) + predictionGroup = kMeansStruct.predict(reshapedArraySimResponse) + if kMeansStruct.cluster_centers_[0] <= kMeansStruct.cluster_centers_[1]: + if skewnessResponse < 0.0: + sensitiveGroup[tempDictKey] = moleFractionRange[predictionGroup==0,:] + sensitiveResponse[tempDictKey] = reshapedArraySimResponse[predictionGroup==0,:] + colorGroup = ["#43aa8b","#f94144"] + else: + sensitiveGroup[tempDictKey] = moleFractionRange[predictionGroup==1,:] + sensitiveResponse[tempDictKey] = reshapedArraySimResponse[predictionGroup==1,:] + colorGroup = ["#f94144","#43aa8b"] + + else: + if skewnessResponse < 0.0: + sensitiveGroup[tempDictKey] = moleFractionRange[predictionGroup==1,:] + sensitiveResponse[tempDictKey] = reshapedArraySimResponse[predictionGroup==1,:] + colorGroup = ["#f94144","#43aa8b"] + else: + sensitiveGroup[tempDictKey] = moleFractionRange[predictionGroup==0,:] + sensitiveResponse[tempDictKey] = reshapedArraySimResponse[predictionGroup==0,:] + colorGroup = ["#43aa8b","#f94144"] + + # Plot raw response in a ternary plot + fig, tax = ternary.figure(scale=1) + fig.set_size_inches(4,3.3) + tax.boundary(linewidth=1.0) + tax.gridlines(multiple=.2, color="gray") + tax.scatter(moleFractionRange, marker='o', s=2, c=arraySimResponse[:,ii], + vmin=min(arraySimResponse[:,ii]),vmax=max(arraySimResponse[:,ii]), + colorbar=True,colormap=plt.cm.PuOr, cmap=plt.cm.PuOr, + cbarlabel = '$m_i$ [g kg$^{-1}$]') + tax.left_axis_label("$y_2$ [-]",offset=0.20,fontsize=10) + tax.right_axis_label("$y_1$ [-]",offset=0.20,fontsize=10) + tax.bottom_axis_label("$y_3$ [-]",offset=0.20,fontsize=10) + tax.ticks(axis='lbr', linewidth=1, multiple=0.2, tick_formats="%.1f", + offset=0.035,clockwise=True,fontsize=10) + tax.clear_matplotlib_ticks() + tax._redraw_labels() + plt.axis('off') + if saveFlag: + # FileName: SensorResponse___ + sensorText = str(sensorID[ii]).replace('[','').replace(']','').replace(' ','-') + saveFileName = "SensorResponse_" + sensorText + "_" + currentDT + "_" + gitCommitID + saveFileExtension + savePath = os.path.join('simulationFigures',saveFileName) + # Check if inputResources directory exists or not. If not, create the folder + if not os.path.exists(os.path.join('..','simulationFigures')): + os.mkdir(os.path.join('..','simulationFigures')) + plt.savefig (savePath) + tax.show() + + # Plot the region of sensitivity in a ternary plot overlayed for all the + # sensors + fig, tax = ternary.figure(scale=1) + fig.set_size_inches(4,3.3) + tax.boundary(linewidth=1.0) + tax.gridlines(multiple=.2, color="gray") + for ii in range(len(sensitiveGroup)): + tempDictKey = 's_{}'.format(ii) + tax.scatter(sensitiveGroup[tempDictKey], marker='o', s=2, + color = '#'+colorIntersection[ii], alpha = 0.15) + tax.scatter(inputMoleFracALL, marker='o', s=20, + color = '#'+colorsForPlot[0]) + inputMoleFracALL[:,[1,2]]=inputMoleFracALL[:,[2,1]] + tax.scatter(inputMoleFracALL, marker='o', s=20, + color = '#'+colorsForPlot[1]) + inputMoleFracALL[:,[0,1]]=inputMoleFracALL[:,[1,0]] + tax.scatter(inputMoleFracALL, marker='o', s=20, + color = '#'+colorsForPlot[2]) + tax.left_axis_label("$y_2$ [-]",offset=0.20,fontsize=10) + tax.right_axis_label("$y_1$ [-]",offset=0.20,fontsize=10) + tax.bottom_axis_label("$y_3$ [-]",offset=0.20,fontsize=10) + tax.ticks(axis='lbr', linewidth=1, multiple=0.2, tick_formats="%.1f", + offset=0.035,clockwise=True,fontsize=10) + tax.clear_matplotlib_ticks() + tax._redraw_labels() + plt.axis('off') + if saveFlag: + # FileName: SensorRegion___ + sensorText = str(saveFileSensorText).replace('[','').replace(']','').replace(' ','-') + saveFileName = "SensorRegion_" + sensorText + "_" + currentDT + "_" + gitCommitID + saveFileExtension + savePath = os.path.join('simulationFigures',saveFileName) + # Check if inputResources directory exists or not. If not, create the folder + if not os.path.exists(os.path.join('..','simulationFigures')): + os.mkdir(os.path.join('..','simulationFigures')) + plt.savefig (savePath) + tax.show() + +# Plot the objective function used to evaluate the concentration for individual +# sensors and the total (sum) +if numberOfGases == 3 and fixOneGas == True: + fig = plt.figure + ax = plt.subplot(1,3,1) + for kk in range(arraySimResponse.shape[1]): + ax.plot(moleFractionRange[:,0],np.power((arrayTrueResponse[:,kk] + -arraySimResponse[:,kk])/arrayTrueResponse[:,kk],2), + color='#'+colorsForPlot[kk], label = '$J_'+str(kk+1)+'$') + ax.plot(moleFractionRange[:,0],objFunction,color='#'+colorsForPlot[-1], label = '$\Sigma J_i$') # Error all sensors + ax.locator_params(axis="x", nbins=4) + ax.locator_params(axis="y", nbins=4) + ax.set(xlabel='$y_1$ [-]', + ylabel='$J$ [-]', + xlim = [0,1.], ylim = [0, 5]) + ax.legend() + + ax = plt.subplot(1,3,2) + for kk in range(arraySimResponse.shape[1]): + ax.plot(moleFractionRange[:,1],np.power((arrayTrueResponse[:,kk] + -arraySimResponse[:,kk])/arrayTrueResponse[:,kk],2), + color='#'+colorsForPlot[kk], label = '$J_'+str(kk+1)+'$') + ax.plot(moleFractionRange[:,1],objFunction,color='#'+colorsForPlot[-1], label = '$\Sigma J_i$') # Error all sensors + ax.locator_params(axis="x", nbins=4) + ax.locator_params(axis="y", nbins=4) + ax.set(xlabel='$y_2$ [-]', + ylabel='$J$ [-]', + xlim = [0,1.], ylim = [0, 5]) + + ax = plt.subplot(1,3,3) + for kk in range(arraySimResponse.shape[1]): + ax.plot(moleFractionRange[:,2],np.power((arrayTrueResponse[:,kk] + -arraySimResponse[:,kk])/arrayTrueResponse[:,kk],2), + color='#'+colorsForPlot[kk], label = '$J_'+str(kk+1)+'$') + ax.plot(moleFractionRange[:,2],objFunction,color='#'+colorsForPlot[-1], label = '$\Sigma J_i$') # Error all sensors + ax.locator_params(axis="x", nbins=4) + ax.locator_params(axis="y", nbins=4) + ax.set(xlabel='$y_3$ [-]', + ylabel='$J$ [-]', + xlim = [0,1.], ylim = [0, None]) + + plt.show() \ No newline at end of file diff --git a/plotFunctions/plotSSLIsotherm.py b/plotFunctions/plotSSLIsotherm.py new file mode 100755 index 0000000..9039743 --- /dev/null +++ b/plotFunctions/plotSSLIsotherm.py @@ -0,0 +1,126 @@ +############################################################################ +# +# Imperial College London, United Kingdom +# Multifunctional Nanomaterials Laboratory +# +# Project: ERASE +# Year: 2020 +# Python: Python 3.7 +# Authors: Ashwin Kumar Rajagopalan (AK) +# +# Purpose: +# Plots the single site Langmuir isotherms +# +# Last modified: +# - 2020-11-04, AK: Plot 3rd gas +# - 2020-10-28, AK: Minor fix for save file name +# - 2020-10-28, AK: Add auxiliary functions as a module +# - 2020-10-27, AK: Further improvements and cosmetic changes +# - 2020-10-26, AK: Initial creation +# +# Input arguments: +# +# +# Output arguments: +# +# +############################################################################ + +import numpy as np +import os +import auxiliaryFunctions +from numpy import load +from simulateSSL import simulateSSL +import matplotlib.pyplot as plt +plt.style.use('singleColumn.mplstyle') # Custom matplotlib style file + +# Flag for saving figure +saveFlag = False + +# Save file extension (png or pdf) +saveFileExtension = ".png" + +# Sensor ID to be plotted +sensorID = 17 + +# Number of gases to determine the hypothetical sorbents roster +numGases = 2 + +# Total pressure of the gas [Pa] +pressureTotal = np.array([1.e5]) + +# Temperature of the gas [K] +# Can be a vector of temperatures +temperature = np.array([298.15]) + +# Molefraction +moleFraction = np.array([np.linspace(0,1,101)]) + +# For now load a given adsorbent isotherm material file +if numGases == 2: + loadFileName = "isothermParameters_20201020_1756_5f263af.npz" # Two gases +elif numGases == 3: + loadFileName = "isothermParameters_20201022_1056_782efa3.npz" # Three gases + +hypoAdsorbentFile = os.path.join('..','inputResources',loadFileName); + +# Check if the file with the adsorbent properties exist +if os.path.exists(hypoAdsorbentFile): + loadedFileContent = load(hypoAdsorbentFile) + adsorbentIsothermTemp = loadedFileContent['adsIsotherm'] + adsorbentDensityTemp = loadedFileContent['adsDensity'] + molecularWeight = loadedFileContent['molWeight'] +else: + errorString = "Adsorbent property file " + hypoAdsorbentFile + " does not exist." + raise Exception(errorString) + +###### TO DO: SERIOUS ISSUE WITH THE ISOTHERM PLOTTING +# Evaluate the isotherms +adsorbentID = np.array([sensorID]) # Do this for consistency +adsorbentIsotherm = adsorbentIsothermTemp[:,:,adsorbentID] +adsorbentDensity = adsorbentDensityTemp[adsorbentID] +equilibriumLoadings = np.zeros([moleFraction.shape[1],adsorbentIsotherm.shape[0]]) +# Loop through all the gases so that the single component isotherm is +# generated. If not multicomponent genretaed. Additionally, several +# transpose operations are performed to be self-consistent with other codes +for ii in range(adsorbentIsotherm.shape[0]): + equilibriumLoadings[:,ii] = np.squeeze(simulateSSL(adsorbentIsotherm[ii,:,:].T,adsorbentDensity, + pressureTotal,temperature,moleFraction.T))/adsorbentDensity # [mol/m3] + +# Get the commit ID of the current repository +gitCommitID = auxiliaryFunctions.getCommitID() + +# Get the current date and time for saving purposes +currentDT = auxiliaryFunctions.getCurrentDateTime() + +# Git commit id of the loaded isotherm file +gitCommmitID_loadedFile = hypoAdsorbentFile[-11:-4] + +# Plot the pure single component isotherm for the n gases +fig = plt.figure +ax = plt.gca() +# HARD CODED for 3 gases +ax.plot(pressureTotal*moleFraction.T/1.e5, equilibriumLoadings[:,0], + linewidth=1.5,color='r', label = '$g_1$') +ax.plot(pressureTotal*moleFraction.T/1.e5, equilibriumLoadings[:,1], + linewidth=1.5,color='b', label = '$g_2$') +if numGases == 3: + ax.plot(pressureTotal*moleFraction.T/1.e5, equilibriumLoadings[:,2], + linewidth=1.5,color='g', label = '$g_3$') +ax.set(xlabel='$P$ [bar]', + ylabel='$q^*$ [mol kg$^{\mathregular{-1}}$]', + xlim = [0, 1], ylim = [0, 0.01]) +ax.legend() + +# Save the figure +if saveFlag: + # FileName: PureIsotherm____ + saveFileName = "PureIsotherm_" + str(sensorID) + "_" + currentDT + "_" + gitCommitID + "_" + gitCommmitID_loadedFile + saveFileExtension + savePath = os.path.join('..','simulationFigures',saveFileName) + # Check if inputResources directory exists or not. If not, create the folder + if not os.path.exists(os.path.join('..','simulationFigures')): + os.mkdir(os.path.join('..','simulationFigures')) + plt.savefig (savePath) + +# For the figure to be saved show should appear after the save +plt.show() \ No newline at end of file diff --git a/plotFunctions/plotsForArticle_Simulation.py b/plotFunctions/plotsForArticle_Simulation.py new file mode 100644 index 0000000..9a262d0 --- /dev/null +++ b/plotFunctions/plotsForArticle_Simulation.py @@ -0,0 +1,1286 @@ +############################################################################ +# +# Imperial College London, United Kingdom +# Multifunctional Nanomaterials Laboratory +# +# Project: ERASE +# Year: 2020 +# Python: Python 3.7 +# Authors: Ashwin Kumar Rajagopalan (AK) +# +# Purpose: +# Plots for the simulation manuscript +# +# Last modified: +# - 2021-08-13, AK: Add plot for absolute sensor response +# - 2021-05-03, AK: Cosmetic changes to all plots +# - 2021-04-07, AK: Add plot for design variables +# - 2021-03-05, AK: Add plot for full model +# - 2021-03-05, AK: Add plot for three materials +# - 2021-03-05, AK: Add plot for comparing sensor array with graphical tool +# - 2021-02-24, AK: Add function to generate sensitive region for each material +# - 2021-02-23, AK: Add mean error to sensor shape plot +# - 2021-02-11, AK: Initial creation +# +# Input arguments: +# +# +# Output arguments: +# +# +############################################################################ + +def plotsForArticle_Simulation(**kwargs): + import auxiliaryFunctions + + # Get the commit ID of the current repository + gitCommitID = auxiliaryFunctions.getCommitID() + + # Get the current date and time for saving purposes + currentDT = auxiliaryFunctions.getCurrentDateTime() + + # Flag for saving figure + if 'saveFlag' in kwargs: + if kwargs["saveFlag"]: + saveFlag = kwargs["saveFlag"] + else: + saveFlag = False + + # Save file extension (png or pdf) + if 'saveFileExtension' in kwargs: + if kwargs["saveFileExtension"]: + saveFileExtension = kwargs["saveFileExtension"] + else: + saveFileExtension = ".png" + + # If sensor array plot needs to be plotted + if 'sensorArray' in kwargs: + if kwargs["sensorArray"]: + plotForArticle_SensorArray(gitCommitID, currentDT, + saveFlag, saveFileExtension) + + # If sensor response curve needs to be plotted + if 'responseShape' in kwargs: + if kwargs["responseShape"]: + plotForArticle_ResponseShape(gitCommitID, currentDT, + saveFlag, saveFileExtension) + + # If graphical tool needs to be plotted + if 'graphicalTool' in kwargs: + if kwargs["graphicalTool"]: + plotForArticle_GraphicalTool(gitCommitID, currentDT, + saveFlag, saveFileExtension) + + # If absolute sensor response curve needs to be plotted + if 'absoluteResponse' in kwargs: + if kwargs["absoluteResponse"]: + plotForArticle_AbsoluteResponse(gitCommitID, currentDT, + saveFlag, saveFileExtension) + + # If three materials needs to be plotted + if 'threeMaterials' in kwargs: + if kwargs["threeMaterials"]: + plotForArticle_ThreeMaterials(gitCommitID, currentDT, + saveFlag, saveFileExtension) + + # If kinetic importance needs to be plotted + if 'kineticsImportance' in kwargs: + if kwargs["kineticsImportance"]: + plotForArticle_KineticsImportance(gitCommitID, currentDT, + saveFlag, saveFileExtension) + + # If kinetic importance needs to be plotted + if 'designVariables' in kwargs: + if kwargs["designVariables"]: + plotForArticle_DesignVariables(gitCommitID, currentDT, + saveFlag, saveFileExtension) + +# fun: plotForArticle_SensorArray +# Plots the histogram of gas compositions for a one and two material +# sensor array +def plotForArticle_SensorArray(gitCommitID, currentDT, + saveFlag, saveFileExtension): + import numpy as np + from numpy import load + import os + import matplotlib.pyplot as plt + plt.style.use('doubleColumn2Row.mplstyle') # Custom matplotlib style file + # For now load the results file + loadFileName = ("arrayConcentration_20210211_1818_b02f8c3.npz", # 1 material w/o constraint + "arrayConcentration_20210212_1055_b02f8c3.npz", # 2 material w/o constraint + "arrayConcentration_20210212_1050_b02f8c3.npz", # 1 material with constraint + "arrayConcentration_20210211_1822_b02f8c3.npz") # 2 material with constraint + + # Git commit id of the loaded file + simID_loadedFile = loadFileName[0][-11:-4] + + # Loop through the two files to get the histogram + for ii in range(len(loadFileName)): + # Create the file name with the path to be loaded + simResultsFile = os.path.join('..','simulationResults',loadFileName[ii]); + + # Check if the file with the adsorbent properties exist + if os.path.exists(simResultsFile): + resultOutput = load(simResultsFile)["arrayConcentration"][:,:] + if resultOutput.shape[1] == 3: + resultOutput = np.delete(resultOutput,[0],1) + elif resultOutput.shape[1] == 4: + resultOutput = np.delete(resultOutput,[0,1],1) + else: + errorString = "Simulation result file " + simResultsFile + " does not exist." + raise Exception(errorString) + + # Gas concentration + molFracG1 = 0.05 + molFracG2 = 0.95 + + # Xlimits and Ylimits + xLimits = [0,1] + yLimits = [0,60] + + # Histogram properties + nBins = 50 + rangeX = (xLimits) + histTypeX = 'stepfilled' + alphaX=0.75 + densityX = True + + # Plot the histogram of the gas compositions + ax = plt.subplot(2,2,ii+1) + # Histogram for 1 material array + ax.axvline(x=molFracG1, linewidth=1, linestyle='dotted', color = '#e5383b', alpha = 0.6) + ax.hist(resultOutput[:,0], bins = nBins, range = rangeX, density = densityX, + linewidth=1.5, histtype = histTypeX, color='#e5383b', alpha = alphaX, label = '$g_1$') + + # Histogram for 2 material array + ax.axvline(x=molFracG2, linewidth=1, linestyle='dotted', color = '#343a40', alpha = 0.6) + ax.hist(resultOutput[:,1], bins = nBins, range = rangeX, density = densityX, + linewidth=1.5, histtype = histTypeX, color='#343a40', alpha = alphaX, label = '$g_2$') + + ax.set(xlim = xLimits, ylim = yLimits) + ax.locator_params(axis="x", nbins=4) + ax.locator_params(axis="y", nbins=4) + if ii == 0 or ii == 2: + ax.set(ylabel='$f$ [-]') + ax.text(0.82, 55, "$n$ = 1", fontsize=10, + backgroundcolor = 'w', color = '#0077b6') + if ii == 0: + ax.text(0.075, 55, "(a)", fontsize=10, + backgroundcolor = 'w') + ax.text(0.53, 51, "Without Constraint", fontsize=10, + backgroundcolor = 'w', color = '#0077b6') + else: + ax.set(xlabel='$y$ [-]') + ax.text(0.075, 55, "(c)", fontsize=10, + backgroundcolor = 'w') + ax.text(0.595, 51, "With Constraint", fontsize=10, + backgroundcolor = 'w', color = '#0077b6') + ax.text(0.085, 25, "$y_1$ = 0.05", fontsize=10, + backgroundcolor = 'w', color = '#e5383b') + ax.text(0.705, 25, "$y_2$ = 0.95", fontsize=10, + backgroundcolor = 'w', color = '#343a40') + elif ii == 1 or ii == 3: + ax.text(0.81, 55, "$n$ = 2", fontsize=10, + backgroundcolor = 'w', color = '#0077b6') + if ii == 1: + ax.text(0.075, 55, "(b)", fontsize=10, + backgroundcolor = 'w') + ax.text(0.53, 51, "Without Constraint", fontsize=10, + backgroundcolor = 'w', color = '#0077b6') + else: + ax.set(xlabel='$y$ [-]') + ax.text(0.075, 55, "(d)", fontsize=10, + backgroundcolor = 'w') + ax.text(0.595, 51, "With Constraint", fontsize=10, + backgroundcolor = 'w', color = '#0077b6') + ax.text(0.085, 25, "$y_1$ = 0.05", fontsize=10, + backgroundcolor = 'w', color = '#e5383b') + ax.text(0.705, 25, "$y_2$ = 0.95", fontsize=10, + backgroundcolor = 'w', color = '#343a40') + # Remove tick labels + if ii == 0 or ii == 1: + ax.axes.xaxis.set_ticklabels([]) + if ii == 1 or ii == 3: + ax.axes.yaxis.set_ticklabels([]) + + # Save the figure + if saveFlag: + # FileName: sensorArray___ + saveFileName = "sensorArray_" + currentDT + "_" + gitCommitID + "_" + simID_loadedFile + saveFileExtension + savePath = os.path.join('..','simulationFigures','simulationManuscript',saveFileName) + # Check if inputResources directory exists or not. If not, create the folder + if not os.path.exists(os.path.join('..','simulationFigures','simulationManuscript')): + os.mkdir(os.path.join('..','simulationFigures','simulationManuscript')) + plt.savefig (savePath) + + # For the figure to be saved show should appear after the save + plt.show() + +# fun: plotForArticle_SensorArray +# Plots the sensor response for a given sensor array +def plotForArticle_ResponseShape(gitCommitID, currentDT, + saveFlag, saveFileExtension): + import numpy as np + import os + import pandas as pd + import seaborn as sns + import matplotlib.pyplot as plt + plt.style.use('doubleColumn.mplstyle') # Custom matplotlib style file + + # Materials to be plotted + sensorID = np.array([17,16,6]) + sensorText = ["A", "B", "C"] + + # File to be loaded for the simulation results + # No Noise + # loadFileName = ["sensitivityAnalysis_17_20210212_1259_b02f8c3.npz", # No Noise + # "sensitivityAnalysis_16_20210212_1300_b02f8c3.npz", # No Noise + # "sensitivityAnalysis_6_20210212_1259_b02f8c3.npz"] # No Noise + # Noise (0.1 g/kg) + loadFileName = ["sensitivityAnalysis_17_20210706_2258_ecbbb3e.npz", # Noise + "sensitivityAnalysis_16_20210707_0842_ecbbb3e.npz", # Noise + "sensitivityAnalysis_6_20210707_1125_ecbbb3e.npz"] # Noise + + # Colors for plot + colorsForPlot = ("#5fad56","#f78154","#b4436c") + + # Get the sensor response and the sensor sensitive region + os.chdir("..") + moleFractionRange, arraySimResponse, _ = getSensorSensitiveRegion(sensorID) + os.chdir("plotFunctions") + + # Plot the figure + fig = plt.figure + ax1 = plt.subplot(1,3,1) + # Loop through all sensors + for kk in range(arraySimResponse.shape[1]): + ax1.plot(moleFractionRange[:,0],arraySimResponse[:,kk], + color=colorsForPlot[kk]) # Simulated Response + + ax1.set(xlabel='$y_1$ [-]', + ylabel='$m$ [g kg$^{-1}$]', + xlim = [0,1], ylim = [0, 300]) + ax1.locator_params(axis="x", nbins=4) + ax1.locator_params(axis="y", nbins=4) + ax1.text(0.05, 270, "(a)", fontsize=10, + backgroundcolor = 'w') + + # Label for the materials + ax1.text(0.85, 225, sensorText[0], fontsize=10, + backgroundcolor = 'w', color = colorsForPlot[0]) + ax1.text(0.1, 150, sensorText[1], fontsize=10, + backgroundcolor = 'w', color = colorsForPlot[1]) + ax1.text(0.8, 75, sensorText[2], fontsize=10, + backgroundcolor = 'w', color = colorsForPlot[2]) + + # Call the concatenateConcEstimate function + meanErr, cvData = concatenateConcEstimate(loadFileName[0:3],sensorText) + + # Mean Error - No noise + ax2 = plt.subplot(1,3,2) + meanErr["x"] = pd.to_numeric(meanErr["x"], downcast="float") + sns.lineplot(data=meanErr, x='x', y='y1', hue='dataType', style='dataType', + dashes = [(1,1),(1,1),(1,1)], markers = ['o','s','D'], + palette = colorsForPlot[0:len(loadFileName)], linewidth = 0.5, + markersize = 5) + + ax2.set(xlabel='$y_1$ [-]', + ylabel='$\psi$ [-]', + xlim = [0.,1.], ylim = [1e-8,100]) + ax2.locator_params(axis="x", nbins=4) + ax2.set_yscale('log') + plt.legend([],[], frameon=False) + ax2.text(0.05, 8, "(b)", fontsize=10, + backgroundcolor = 'w') + + # Label for the materials + ax2.text(0.85, 1e-2, sensorText[0], fontsize=10, + backgroundcolor = 'w', color = colorsForPlot[0]) + ax2.text(0.3, 4e-4, sensorText[1], fontsize=10, + backgroundcolor = 'w', color = colorsForPlot[1]) + ax2.text(0.6, 3e-6, sensorText[2], fontsize=10, + backgroundcolor = 'w', color = colorsForPlot[2]) + + # Label for the formula + ax2.text(0.38, 7, "$\psi = |\mu - \hat{\mu}|/\mu$", fontsize=10, + backgroundcolor = 'w', color = '#0077b6') + + # CV - No noise + ax3 = plt.subplot(1,3,3) + cvData["x"] = pd.to_numeric(cvData["x"], downcast="float") + sns.lineplot(data=cvData, x='x', y='y1', hue='dataType', style='dataType', + dashes = [(1,1),(1,1),(1,1)], markers = ['o','s','D'], + palette = colorsForPlot[0:len(loadFileName)], linewidth = 0.5, + markersize = 5) + + ax3.set(xlabel='$y_1$ [-]', + ylabel='$\chi$ [-]', + xlim = [0.,1.], ylim = [1e-8,100]) + ax3.locator_params(axis="x", nbins=4) + ax3.set_yscale('log') + plt.legend([],[], frameon=False) + ax3.text(0.05, 8, "(c)", fontsize=10, + backgroundcolor = 'w') + + # Label for the materials + ax3.text(0.85, 3e-1, sensorText[0], fontsize=10, + backgroundcolor = 'w', color = colorsForPlot[0]) + ax3.text(0.81, 4e-5, sensorText[1], fontsize=10, + backgroundcolor = 'w', color = colorsForPlot[1]) + ax3.text(0.6, 3e-4, sensorText[2], fontsize=10, + backgroundcolor = 'w', color = colorsForPlot[2]) + + # Label for the formula + ax3.text(0.62, 7, "$\chi = \hat{\sigma}/\hat{\mu}$", fontsize=10, + backgroundcolor = 'w', color = '#0077b6') + + # Save the figure + if saveFlag: + # FileName: responseShape___ + sensorText = str(sensorID).replace('[','').replace(']','').replace(' ','-').replace(' ','-') + saveFileName = "responseShape_" + sensorText + "_" + currentDT + "_" + gitCommitID + saveFileExtension + savePath = os.path.join('..','simulationFigures','simulationManuscript',saveFileName) + # Check if inputResources directory exists or not. If not, create the folder + if not os.path.exists(os.path.join('..','simulationFigures','simulationManuscript')): + os.mkdir(os.path.join('..','simulationFigures','simulationManuscript')) + plt.savefig (savePath) + plt.show() + +# fun: plotForArticle_AbsoluteResponse +# Plots the sensor response for a given sensor array +def plotForArticle_AbsoluteResponse(gitCommitID, currentDT, + saveFlag, saveFileExtension): + import numpy as np + import os + import pandas as pd + import seaborn as sns + import matplotlib.pyplot as plt + plt.style.use('doubleColumn.mplstyle') # Custom matplotlib style file + + # Materials to be plotted + sensorID = np.array([6,2,2,2]) + multiplierError = np.array([1,1,3,10]) + sensorText = ["$\\alpha$", "$\\beta$", "$\\beta_3$", "$\\beta_{10}$"] + arrayText = ["D", "D$_3$", "D$_{10}$"] + + # File to be loaded for the simulation results + # Noise (0.1 g/kg) + loadFileName = ["sensitivityAnalysis_6-2_20210719_1117_ecbbb3e.npz", # 1 + "sensitivityAnalysis_6-2_20210719_2145_ecbbb3e.npz", # 3 + "sensitivityAnalysis_6-2_20210719_1458_ecbbb3e.npz"] # 10 + + # Colors for plot + colorsForPlot = ("#5fad56","#FF9E00","#D66612","#AD2E24") + + # Get the sensor response and the sensor sensitive region + os.chdir("..") + moleFractionRange, arraySimResponse, _ = getSensorSensitiveRegion(sensorID) + os.chdir("plotFunctions") + + # Plot the figure + fig = plt.figure + ax1 = plt.subplot(1,3,1) + # Loop through all sensors + for kk in range(arraySimResponse.shape[1]): + ax1.plot(moleFractionRange[:,0],arraySimResponse[:,kk]*multiplierError[kk], + color=colorsForPlot[kk]) # Simulated Response + + ax1.set(xlabel='$y_1$ [-]', + ylabel='$m$ [g kg$^{-1}$]', + xlim = [0,1], ylim = [0, 150]) + ax1.locator_params(axis="x", nbins=4) + ax1.locator_params(axis="y", nbins=4) + ax1.text(0.05, 135, "(a)", fontsize=10, + backgroundcolor = 'w') + + # Label for the materials + ax1.text(0.30, 75, sensorText[0], fontsize=10, + color = colorsForPlot[0]) + ax1.text(0.8, 17, sensorText[1], fontsize=10, + color = colorsForPlot[1]) + ax1.text(0.8, 40, sensorText[2], fontsize=10, + color = colorsForPlot[2]) + ax1.text(0.8, 130, sensorText[3], fontsize=10, + color = colorsForPlot[3]) + + # Call the concatenateConcEstimate function + meanErr, cvData = concatenateConcEstimate(loadFileName[0:3],sensorText) + + # Mean Error - No noise + ax2 = plt.subplot(1,3,2) + meanErr["x"] = pd.to_numeric(meanErr["x"], downcast="float") + sns.lineplot(data=meanErr, x='x', y='y1', hue='dataType', style='dataType', + dashes = [(1,1),(1,1),(1,1)], markers = ['o','s','D'], + palette = colorsForPlot[1:len(loadFileName)+1], linewidth = 0.5, + markersize = 5) + + ax2.set(xlabel='$y_1$ [-]', + ylabel='$\psi$ [-]', + xlim = [0.,1.], ylim = [1e-6,1]) + ax2.locator_params(axis="x", nbins=4) + ax2.set_yscale('log') + plt.legend([],[], frameon=False) + ax2.text(0.05, 0.25, "(b)", fontsize=10,) + plt.minorticks_off() + + # Label for the materials + ax2.text(0.85, 8e-4, arrayText[0], fontsize=10, + backgroundcolor = 'w', color = colorsForPlot[1]) + ax2.text(0.3, 1.2e-4, arrayText[1], fontsize=10, + color = colorsForPlot[2]) + ax2.text(0.63, 3e-6, arrayText[2], fontsize=10, + backgroundcolor = 'w', color = colorsForPlot[3]) + + # Label for the formula + ax2.text(0.38, 0.25, "$\psi = |\mu - \hat{\mu}|/\mu$", fontsize=10, + backgroundcolor = 'w', color = '#0077b6') + + # CV - No noise + ax3 = plt.subplot(1,3,3) + cvData["x"] = pd.to_numeric(cvData["x"], downcast="float") + sns.lineplot(data=cvData, x='x', y='y1', hue='dataType', style='dataType', + dashes = [(1,1),(1,1),(1,1)], markers = ['o','s','D'], + palette = colorsForPlot[1:len(loadFileName)+1], linewidth = 0.5, + markersize = 5) + + ax3.set(xlabel='$y_1$ [-]', + ylabel='$\chi$ [-]', + xlim = [0.,1.], ylim = [1e-6,1]) + ax3.locator_params(axis="x", nbins=4) + ax3.set_yscale('log') + plt.legend([],[], frameon=False) + ax3.text(0.05, 0.25, "(c)", fontsize=10,) + plt.minorticks_off() + + # Label for the materials + ax3.text(0.8, 1.3e-2, arrayText[0], fontsize=10, + color = colorsForPlot[1]) + ax3.text(0.8, 3e-4, arrayText[2], fontsize=10, + color = colorsForPlot[3]) + + # Label for the formula + ax3.text(0.62, 0.25, "$\chi = \hat{\sigma}/\hat{\mu}$", fontsize=10, + backgroundcolor = 'w', color = '#0077b6') + + # Save the figure + if saveFlag: + # FileName: responseShape___ + sensorText = str(sensorID).replace('[','').replace(']','').replace(' ','-').replace(' ','-') + saveFileName = "absoluteResponse_" + currentDT + "_" + gitCommitID + saveFileExtension + savePath = os.path.join('..','simulationFigures','simulationManuscript',saveFileName) + # Check if inputResources directory exists or not. If not, create the folder + if not os.path.exists(os.path.join('..','simulationFigures','simulationManuscript')): + os.mkdir(os.path.join('..','simulationFigures','simulationManuscript')) + plt.savefig (savePath) + plt.show() + +# fun: plotForArticle_GraphicalTool +# Plots the graphical tool to screen for materials +def plotForArticle_GraphicalTool(gitCommitID, currentDT, + saveFlag, saveFileExtension): + import numpy as np + import os + import pandas as pd + import seaborn as sns + import matplotlib.pyplot as plt + plt.style.use('doubleColumn2Row.mplstyle') # Custom matplotlib style file + + # File to be loaded for the simulation results + # No noise + # loadFileName = ["sensitivityAnalysis_6-2_20210305_1109_b02f8c3.npz", # 6,2 + # "sensitivityAnalysis_17-16_20210305_1050_b02f8c3.npz"] # 17,16 + # Noise + loadFileName = ["sensitivityAnalysis_6-2_20210719_1117_ecbbb3e.npz", # 6,2 + "sensitivityAnalysis_17-16_20210706_2120_ecbbb3e.npz"] # 17,16 + + # Materials to be plotted + sensorID = np.array([[6,2],[17,16]]) + arrayText = ["D", "E"] + materialText = ["$\\alpha$", "$\\beta$", "$\gamma$", "$\delta$"] + + # Colors for plot + colorsForPlot = ("#5fad56","#ff9e00") + colorLeft = ("#e5383b","#6c757d") + colorRight = ("#6c757d","#e5383b") + + # Plot the figure + fig = plt.figure + ax1 = plt.subplot(2,2,1) + # Get the sensor response and the sensor sensitive region + os.chdir("..") + moleFractionRange, arraySimResponse, sensitiveRegion = getSensorSensitiveRegion(sensorID[0,:]) + os.chdir("plotFunctions") + # Loop through all sensors + for kk in range(arraySimResponse.shape[1]): + ax1.plot(moleFractionRange[:,0],arraySimResponse[:,kk], + color=colorsForPlot[kk]) # Simulated Response + ax1.fill_between(sensitiveRegion[kk,:],1.5*np.max(arraySimResponse), + facecolor=colorsForPlot[kk], alpha=0.25) + + ax1.set(ylabel='$m$ [g kg$^{-1}$]', + xlim = [0,1], ylim = [0, 150]) + ax1.axes.xaxis.set_ticklabels([]) + ax1.locator_params(axis="x", nbins=4) + ax1.locator_params(axis="y", nbins=4) + ax1.text(0.025, 138, "(a)", fontsize=10) + ax1.text(0.78, 138, "Array D", fontsize=10, + color = '#0077b6') + + # Label for the materials + ax1.text(0.9, 120, materialText[0], fontsize=10, + backgroundcolor = 'w', color = colorsForPlot[0]) + ax1.text(0.9, 23, materialText[1], fontsize=10, + backgroundcolor = 'w', color = colorsForPlot[1]) + + ax2 = plt.subplot(2,2,2) + # Call the concatenateConcEstimate function + meanErr, cvData = concatenateConcEstimate([loadFileName[0]],arrayText[0]) + + # Mean Error - No noise + meanErr["x"] = pd.to_numeric(meanErr["x"], downcast="float") + sns.lineplot(data=meanErr, x='x', y='y1', hue='dataType', style='dataType', + dashes = [(1,1)], markers = ['o'], + palette = colorLeft, linewidth = 0.5,markersize = 6) + ax2.set(ylabel='$\psi$ [-]', + xlim = [0.,1.], ylim = [1e-8,1]) + ax2.locator_params(axis="x", nbins=4) + ax2.set(xlabel=None) + ax2.axes.xaxis.set_ticklabels([]) + ax2.set_yscale('log') + ax2.yaxis.label.set_color(colorLeft[0]) + ax2.tick_params(axis='y', colors=colorLeft[0]) + plt.legend([],[], frameon=False) + # CV - No noise + ax2r = plt.twinx() + cvData["x"] = pd.to_numeric(cvData["x"], downcast="float") + sns.lineplot(data=cvData, x='x', y='y1', hue='dataType', style='dataType', + dashes = [(1,1)], markers = ['D'], + palette = colorRight, linewidth = 0.5,markersize = 6, + ax = ax2r) + # Plot sensitive region + for kk in range(arraySimResponse.shape[1]): + ax2r.fill_between(sensitiveRegion[kk,:],1.5*np.max(arraySimResponse), + facecolor=colorsForPlot[kk], alpha=0.25) + + ax2r.set(ylabel='$\chi$ [-]',ylim = [1e-8,1]) + ax2r.locator_params(axis="x", nbins=4) + ax2r.axes.xaxis.set_ticklabels([]) + ax2r.set_yscale('log') + ax2r.yaxis.label.set_color(colorLeft[1]) + ax2r.tick_params(axis='y', colors=colorLeft[1]) + plt.legend([],[], frameon=False) + ax2r.annotate("", xy=(0.55, 1e-4), xytext=(0.65, 1e-4), + arrowprops=dict(arrowstyle="-|>", color = colorLeft[0])) + ax2r.annotate("", xy=(0.95, 3e-2), xytext=(0.85, 3e-2), + arrowprops=dict(arrowstyle="-|>", color = colorLeft[1])) + ax2r.text(0.025, 0.2, "(b)", fontsize=10) + ax2r.spines["left"].set_color(colorLeft[0]) + ax2r.spines["right"].set_color(colorLeft[1]) + + ax2r.text(0.78, 0.2, "Array D", fontsize=10, + color = '#0077b6') + + ax3 = plt.subplot(2,2,3) + # Get the sensor response and the sensor sensitive region + os.chdir("..") + moleFractionRange, arraySimResponse, sensitiveRegion = getSensorSensitiveRegion(sensorID[1,:]) + os.chdir("plotFunctions") + # Loop through all sensors + for kk in range(arraySimResponse.shape[1]): + ax3.plot(moleFractionRange[:,0],arraySimResponse[:,kk], + color=colorsForPlot[kk]) # Simulated Response + ax3.fill_between(sensitiveRegion[kk,:],1.5*np.max(arraySimResponse), + facecolor=colorsForPlot[kk], alpha=0.25) + + ax3.set(xlabel='$y_1$ [-]', + ylabel='$m$ [g kg$^{-1}$]', + xlim = [0,1], ylim = [0, 300]) + ax3.locator_params(axis="x", nbins=4) + ax3.locator_params(axis="y", nbins=4) + ax3.text(0.025, 275, "(c)", fontsize=10) + ax3.text(0.78, 275, "Array E", fontsize=10, + color = '#0077b6') + + # Label for the materials + ax3.text(0.78, 225, materialText[2], fontsize=10, + backgroundcolor = 'w', color = colorsForPlot[0]) + ax3.text(0.1, 150, materialText[3], fontsize=10, + backgroundcolor = 'w', color = colorsForPlot[1]) + + ax4 = plt.subplot(2,2,4) + # Call the concatenateConcEstimate function + meanErr, cvData = concatenateConcEstimate([loadFileName[1]],arrayText[1]) + + # Mean Error - No noise + meanErr["x"] = pd.to_numeric(meanErr["x"], downcast="float") + sns.lineplot(data=meanErr, x='x', y='y1', hue='dataType', style='dataType', + dashes = [(1,1)], markers = ['o'], + palette = colorLeft, linewidth = 0.5,markersize = 6) + ax4.set(xlabel='$y_1$ [-]', + ylabel='$\psi$ [-]', + xlim = [0.,1.], ylim = [1e-8,1]) + ax4.locator_params(axis="x", nbins=4) + ax4.set_yscale('log') + ax4.yaxis.label.set_color(colorLeft[0]) + ax4.tick_params(axis='y', colors=colorLeft[0]) + plt.legend([],[], frameon=False) + # CV - No noise + ax4r = plt.twinx() + cvData["x"] = pd.to_numeric(cvData["x"], downcast="float") + sns.lineplot(data=cvData, x='x', y='y1', hue='dataType', style='dataType', + dashes = [(1,1)], markers = ['D'], + palette = colorRight, linewidth = 0.5,markersize = 6, + ax = ax4r) + # Plot sensitive region + for kk in range(arraySimResponse.shape[1]): + ax4r.fill_between(sensitiveRegion[kk,:],1.5*np.max(arraySimResponse), + facecolor=colorsForPlot[kk], alpha=0.25) + + ax4r.set(ylabel='$\chi$ [-]',ylim = [1e-8,1]) + ax4r.locator_params(axis="x", nbins=4) + ax4r.set_yscale('log') + ax4r.yaxis.label.set_color(colorLeft[1]) + ax4r.tick_params(axis='y', colors=colorLeft[1]) + plt.legend([],[], frameon=False) + ax4r.annotate("", xy=(0.2, 5e-4), xytext=(0.3, 5e-4), + arrowprops=dict(arrowstyle="-|>", color = colorLeft[0])) + ax4r.annotate("", xy=(0.7, 3e-2), xytext=(0.6, 3e-2), + arrowprops=dict(arrowstyle="-|>", color = colorLeft[1])) + ax4r.text(0.025, 0.2, "(d)", fontsize=10) + ax4r.spines["left"].set_color(colorLeft[0]) + ax4r.spines["right"].set_color(colorLeft[1]) + + ax4r.text(0.78, 0.2, "Array E", fontsize=10, + color = '#0077b6') + + # Save the figure + if saveFlag: + # FileName: graphicalTool__ + saveFileName = "graphicalTool_" + currentDT + "_" + gitCommitID + saveFileExtension + savePath = os.path.join('..','simulationFigures','simulationManuscript',saveFileName) + # Check if inputResources directory exists or not. If not, create the folder + if not os.path.exists(os.path.join('..','simulationFigures','simulationManuscript')): + os.mkdir(os.path.join('..','simulationFigures','simulationManuscript')) + plt.savefig (savePath) + plt.show() + +# fun: plotForArticle_GraphicalTool +# Plots the analaysis for three material arrays +def plotForArticle_ThreeMaterials(gitCommitID, currentDT, + saveFlag, saveFileExtension): + import numpy as np + import os + import pandas as pd + import seaborn as sns + import matplotlib.pyplot as plt + plt.style.use('doubleColumn2Row.mplstyle') # Custom matplotlib style file + + # File to be loaded for the simulation results + # No Noise + # loadFileName = ["sensitivityAnalysis_17-15-6_20210306_1515_b02f8c3.npz", # 17,15,6 + # "sensitivityAnalysis_17-15-16_20210306_1515_b02f8c3.npz", # 17,15,16 + # "sensitivityAnalysis_17-15_20210308_1002_b02f8c3.npz"] # 17,15 + # Noise + loadFileName = ["sensitivityAnalysis_17-15-6_20210707_2036_ecbbb3e.npz", # 17,15,6 + "sensitivityAnalysis_17-15-16_20210708_0934_ecbbb3e.npz", # 17,15,16 + "sensitivityAnalysis_17-15_20210709_1042_ecbbb3e.npz"] # 17,15 + + # Materials to be plotted + sensorID = np.array([[17,15,6],[17,15,16]]) + arrayText = ["F", "G", "Ref"] + materialText = ["$\\alpha$", "$\\beta$", "$\gamma$", "$\delta$", "$\zeta$"] + + # Colors for plot + colorsForPlot = ("#5fad56","#98c1d9","#ff9e00") + colorLeft = ("#e5383b","#6c757d") + colorRight = ("#6c757d","#e5383b") + + # Plot the figure + fig = plt.figure + ax1 = plt.subplot(2,2,1) + # Get the sensor response and the sensor sensitive region + os.chdir("..") + moleFractionRange, arraySimResponse, sensitiveRegion = getSensorSensitiveRegion(sensorID[0,:]) + os.chdir("plotFunctions") + # Loop through all sensors + for kk in range(arraySimResponse.shape[1]): + ax1.plot(moleFractionRange[:,0],arraySimResponse[:,kk], + color=colorsForPlot[kk]) # Simulated Response + ax1.fill_between(sensitiveRegion[kk,:],1.5*np.max(arraySimResponse), + facecolor=colorsForPlot[kk], alpha=0.25) + + ax1.set(ylabel='$m$ [g kg$^{-1}$]', + xlim = [0,1], ylim = [0, 300]) + ax1.axes.xaxis.set_ticklabels([]) + ax1.locator_params(axis="x", nbins=4) + ax1.locator_params(axis="y", nbins=4) + ax1.text(0.025, 275, "(a)", fontsize=10) + ax1.text(0.78, 275, "Array F", fontsize=10, + color = '#0077b6') + + # Label for the materials + ax1.text(0.78, 225, materialText[2], fontsize=10, + backgroundcolor = 'w', color = colorsForPlot[0]) + ax1.text(0.78, 25, materialText[4], fontsize=10, + backgroundcolor = 'w', color = colorsForPlot[1]) + ax1.text(0.78, 80, materialText[0], fontsize=10, + backgroundcolor = 'w', color = colorsForPlot[2]) + + ax2 = plt.subplot(2,2,2) + # Call the concatenateConcEstimate function + meanErrRef, cvDataRef = concatenateConcEstimate([loadFileName[2]],arrayText[2]) + meanErr, cvData = concatenateConcEstimate([loadFileName[0]],arrayText[0]) + + # Mean Error - No noise + meanErr["x"] = pd.to_numeric(meanErr["x"], downcast="float") + sns.lineplot(data=meanErr, x='x', y='y1', hue='dataType', style='dataType', + dashes = [(1,1)], markers = ['o'], + palette = colorLeft, linewidth = 0.5,markersize = 6) + meanErrRef["x"] = pd.to_numeric(meanErrRef["x"], downcast="float") + sns.lineplot(data=meanErrRef, x='x', y='y1', hue='dataType', style='dataType', + dashes = [(5,5)], markers = ['o'], + palette = colorLeft, linewidth = 0.5, alpha = 0.35,markersize = 6) + ax2.set(ylabel='$\psi$ [-]', + xlim = [0.,1.], ylim = [1e-8,1]) + ax2.locator_params(axis="x", nbins=4) + ax2.set(xlabel=None) + ax2.axes.xaxis.set_ticklabels([]) + ax2.set_yscale('log') + ax2.yaxis.label.set_color(colorLeft[0]) + ax2.tick_params(axis='y', colors=colorLeft[0]) + plt.legend([],[], frameon=False) + # CV - No noise + ax2r = plt.twinx() + cvData["x"] = pd.to_numeric(cvData["x"], downcast="float") + sns.lineplot(data=cvData, x='x', y='y1', hue='dataType', style='dataType', + dashes = [(1,1)], markers = ['D'], + palette = colorRight, linewidth = 0.5,markersize = 6, + ax = ax2r) + cvDataRef["x"] = pd.to_numeric(cvDataRef["x"], downcast="float") + sns.lineplot(data=cvDataRef, x='x', y='y1', hue='dataType', style='dataType', + dashes = [(5,5)], markers = ['D'], + palette = colorRight, linewidth = 0.5, alpha = 0.35,markersize = 6, + ax = ax2r) + # Plot sensitive region + for kk in range(arraySimResponse.shape[1]): + ax2r.fill_between(sensitiveRegion[kk,:],1.5*np.max(arraySimResponse), + facecolor=colorsForPlot[kk], alpha=0.25) + + ax2r.set(ylabel='$\chi$ [-]',ylim = [1e-8,1]) + ax2r.locator_params(axis="x", nbins=4) + ax2r.axes.xaxis.set_ticklabels([]) + ax2r.set_yscale('log') + ax2r.yaxis.label.set_color(colorLeft[1]) + ax2r.tick_params(axis='y', colors=colorLeft[1]) + plt.legend([],[], frameon=False) + ax2r.annotate("", xy=(0.4, 3e-5), xytext=(0.5, 3e-5), + arrowprops=dict(arrowstyle="-|>", color = colorLeft[0])) + ax2r.annotate("", xy=(0.95, 5e-3), xytext=(0.85, 5e-3), + arrowprops=dict(arrowstyle="-|>", color = colorLeft[1])) + ax2r.text(0.025, 0.2, "(b)", fontsize=10) + ax2r.spines["left"].set_color(colorLeft[0]) + ax2r.spines["right"].set_color(colorLeft[1]) + + ax2r.text(0.78, 1e-6, "Array F", fontsize=10, + color = '#0077b6') + ax2r.text(0.4, 0.3, "Reference ($\gamma \zeta$)", fontsize=10, + color = '#0077b6', alpha = 0.35) + + ax3 = plt.subplot(2,2,3) + # Get the sensor response and the sensor sensitive region + os.chdir("..") + moleFractionRange, arraySimResponse, sensitiveRegion = getSensorSensitiveRegion(sensorID[1,:]) + os.chdir("plotFunctions") + # Loop through all sensors + for kk in range(arraySimResponse.shape[1]): + ax3.plot(moleFractionRange[:,0],arraySimResponse[:,kk], + color=colorsForPlot[kk]) # Simulated Response + ax3.fill_between(sensitiveRegion[kk,:],1.5*np.max(arraySimResponse), + facecolor=colorsForPlot[kk], alpha=0.25) + + ax3.set(xlabel='$y_1$ [-]', + ylabel='$m$ [g kg$^{-1}$]', + xlim = [0,1], ylim = [0, 300]) + ax3.locator_params(axis="x", nbins=4) + ax3.locator_params(axis="y", nbins=4) + ax3.text(0.025, 275, "(c)", fontsize=10) + ax3.text(0.78, 275, "Array G", fontsize=10, + color = '#0077b6') + + # Label for the materials + ax3.text(0.78, 225, materialText[2], fontsize=10, + backgroundcolor = 'w', color = colorsForPlot[0]) + ax3.text(0.78, 25, materialText[4], fontsize=10, + backgroundcolor = 'w', color = colorsForPlot[1]) + ax3.text(0.1, 150, materialText[3], fontsize=10, + backgroundcolor = 'w', color = colorsForPlot[2]) + + ax4 = plt.subplot(2,2,4) + # Call the concatenateConcEstimate function + meanErr, cvData = concatenateConcEstimate([loadFileName[1]],arrayText[1]) + + # Mean Error - No noise + meanErr["x"] = pd.to_numeric(meanErr["x"], downcast="float") + sns.lineplot(data=meanErr, x='x', y='y1', hue='dataType', style='dataType', + dashes = [(1,1)], markers = ['o'], + palette = colorLeft, linewidth = 0.5,markersize = 6) + sns.lineplot(data=meanErrRef, x='x', y='y1', hue='dataType', style='dataType', + dashes = [(5,5)], markers = ['o'], + palette = colorLeft, linewidth = 0.5, alpha = 0.35,markersize = 6) + ax4.set(xlabel='$y_1$ [-]', + ylabel='$\psi$ [-]', + xlim = [0.,1.], ylim = [1e-8,1]) + ax4.locator_params(axis="x", nbins=4) + ax4.set_yscale('log') + ax4.yaxis.label.set_color(colorLeft[0]) + ax4.tick_params(axis='y', colors=colorLeft[0]) + plt.legend([],[], frameon=False) + # CV - No noise + ax4r = plt.twinx() + cvData["x"] = pd.to_numeric(cvData["x"], downcast="float") + sns.lineplot(data=cvData, x='x', y='y1', hue='dataType', style='dataType', + dashes = [(1,1)], markers = ['D'], + palette = colorRight, linewidth = 0.5,markersize = 6, + ax = ax4r) + sns.lineplot(data=cvDataRef, x='x', y='y1', hue='dataType', style='dataType', + dashes = [(5,5)], markers = ['o'], + palette = colorRight, linewidth = 0.5, alpha = 0.35,markersize = 6, + ax = ax4r) + # Plot sensitive region + for kk in range(arraySimResponse.shape[1]): + ax4r.fill_between(sensitiveRegion[kk,:],1.5*np.max(arraySimResponse), + facecolor=colorsForPlot[kk], alpha=0.25) + + ax4r.set(ylabel='$\chi$ [-]',ylim = [1e-8,1]) + ax4r.locator_params(axis="x", nbins=4) + ax4r.set_yscale('log') + ax4r.yaxis.label.set_color(colorLeft[1]) + ax4r.tick_params(axis='y', colors=colorLeft[1]) + plt.legend([],[], frameon=False) + ax4r.annotate("", xy=(0.08, 5e-4), xytext=(0.18, 5e-4), + arrowprops=dict(arrowstyle="-|>", color = colorLeft[0])) + ax4r.annotate("", xy=(0.72, 1e-2), xytext=(0.62, 1e-2), + arrowprops=dict(arrowstyle="-|>", color = colorLeft[1])) + ax4r.text(0.025, 0.2, "(d)", fontsize=10) + ax4r.spines["left"].set_color(colorLeft[0]) + ax4r.spines["right"].set_color(colorLeft[1]) + + ax4r.text(0.6, 1e-5, "Array G", fontsize=10, + color = '#0077b6') + ax4r.text(0.3, 0.3, "Reference ($\gamma \zeta$)", fontsize=10, + color = '#0077b6', alpha = 0.35) + # Save the figure + if saveFlag: + # FileName: threeMaterials__ + saveFileName = "threeMaterials_" + currentDT + "_" + gitCommitID + saveFileExtension + savePath = os.path.join('..','simulationFigures','simulationManuscript',saveFileName) + # Check if inputResources directory exists or not. If not, create the folder + if not os.path.exists(os.path.join('..','simulationFigures','simulationManuscript')): + os.mkdir(os.path.join('..','simulationFigures','simulationManuscript')) + plt.savefig (savePath) + plt.show() + +# fun: plotForArticle_KineticsImportance +# Plots to highlight the importance of incorporating kinetics +def plotForArticle_KineticsImportance(gitCommitID, currentDT, + saveFlag, saveFileExtension): + import numpy as np + from numpy import load + import os + import matplotlib.pyplot as plt + plt.style.use('doubleColumn.mplstyle') # Custom matplotlib style file + + # Colors for plot + colorsForPlot = ("#5fad56","#ff9e00","#e5383b","#6c757d") + + # Labels for materials + materialText = ["$\\alpha$", "$\\beta$"] + + # File name for equilibrium model and full model estimates + loadFileName_E = "fullModelConcentrationEstimate_6-2_20210320_1336_4b80775.npz" # Eqbm + loadFileName_F = "fullModelConcentrationEstimate_6-2_20210320_1338_4b80775.npz" # Full model + + # Parse out equilbirum results file + simResultsFile = os.path.join('..','simulationResults',loadFileName_E); + loadedFile_E = load(simResultsFile, allow_pickle=True) + concentrationEstimate_E = loadedFile_E["arrayConcentration"] + + # Parse out full model results file + simResultsFile = os.path.join('..','simulationResults',loadFileName_F); + loadedFile_F = load(simResultsFile, allow_pickle=True) + concentrationEstimate_F = loadedFile_F["arrayConcentration"] + + # Parse out true responses (this should be the same for both eqbm and full + # model (here loaded from full model) + trueResponseStruct = loadedFile_F["outputStruct"].item() + # Parse out time + timeSim = [] + timeSim = trueResponseStruct[0]["timeSim"] + # Parse out feed mole fraction + feedMoleFrac = trueResponseStruct[0]["inputParameters"][5] + # Parse out true sensor finger print + sensorFingerPrint = np.zeros([len(timeSim),len(trueResponseStruct)]) + for ii in range(len(trueResponseStruct)): + sensorFingerPrint[:,ii] = trueResponseStruct[ii]["sensorFingerPrint"] + + # Points that will be taken for sampling (for plots) + lenSampling = 6 + fig = plt.figure + # Plot the true sensor response (using the full model) + ax1 = plt.subplot(1,2,1) + ax1.plot(timeSim[0:len(timeSim):lenSampling], + sensorFingerPrint[0:len(timeSim):lenSampling,0], + marker = 'o', markersize = 2, linestyle = 'dotted', linewidth = 0.5, + color=colorsForPlot[0]) + ax1.plot(timeSim[0:len(timeSim):lenSampling], + sensorFingerPrint[0:len(timeSim):lenSampling,1], + marker = 'D', markersize = 2, linestyle = 'dotted', linewidth = 0.5, + color=colorsForPlot[1]) + ax1.locator_params(axis="x", nbins=4) + ax1.locator_params(axis="y", nbins=4) + ax1.set(xlabel='$t$ [s]', + ylabel='$m$ [g kg$^{\mathregular{-1}}$]', + xlim = [0, 1000.], ylim = [0, 40]) + + ax1.text(20, 37, "(a)", fontsize=10) + ax1.text(800, 37, "Array D", fontsize=10, + color = '#0077b6') + ax1.text(720, 33.5, "$y^{\mathregular{in}}_{\mathregular{1}} (t)$ = 0.1", fontsize=10, + color = '#0077b6') + + # Label for the materials + ax1.text(900, 26, materialText[0], fontsize=10, + backgroundcolor = 'w', color = colorsForPlot[0]) + ax1.text(900, 4, materialText[1], fontsize=10, + backgroundcolor = 'w', color = colorsForPlot[1]) + + # Plot the evolution of the gas composition with respect to time + ax2 = plt.subplot(1,2,2) + ax2.plot(timeSim[0:len(timeSim):lenSampling], + concentrationEstimate_E[0:len(timeSim):lenSampling,2], + marker = 'v', markersize = 2, linestyle = 'dotted', linewidth = 0.5, + color=colorsForPlot[2]) + ax2.plot(timeSim[0:len(timeSim):lenSampling], + concentrationEstimate_F[0:len(timeSim):lenSampling,2], + marker = '^', markersize = 2, linestyle = 'dotted', linewidth = 0.5, + color=colorsForPlot[3]) + ax2.locator_params(axis="x", nbins=4) + ax2.locator_params(axis="y", nbins=4) + ax2.set(xlabel='$t$ [s]', + ylabel='$y_1$ [-]', + xlim = [0, 1000.], ylim = [0, 0.2]) + + ax2.text(20, 0.185, "(b)", fontsize=10) + ax2.text(800, 0.185, "Array D", fontsize=10, + color = '#0077b6') + ax2.text(720, 0.11, "$y^{\mathregular{in}}_{\mathregular{1}} (t)$ = 0.1", fontsize=10, + color = '#0077b6') + + # Label for the materials + ax2.text(280, 0.06, "Equilibrium Model", fontsize=10, + backgroundcolor = 'w', color = colorsForPlot[2]) + ax2.text(50, 0.11, "Full Model", fontsize=10, + backgroundcolor = 'w', color = colorsForPlot[3]) + + # Save the figure + if saveFlag: + # FileName: kineticsImportance__ + saveFileName = "kineticsImportance_" + currentDT + "_" + gitCommitID + saveFileExtension + savePath = os.path.join('..','simulationFigures','simulationManuscript',saveFileName) + # Check if inputResources directory exists or not. If not, create the folder + if not os.path.exists(os.path.join('..','simulationFigures','simulationManuscript')): + os.mkdir(os.path.join('..','simulationFigures','simulationManuscript')) + plt.savefig (savePath) + plt.show() + +# fun: plotForArticle_DesignVariables +# Plots to highlight the effect of different design variables +def plotForArticle_DesignVariables(gitCommitID, currentDT, + saveFlag, saveFileExtension): + import numpy as np + from numpy import load + import os + import matplotlib.pyplot as plt + import matplotlib.patches as patches + from matplotlib.transforms import Bbox + + plt.style.use('doubleColumn2Row.mplstyle') # Custom matplotlib style file + + # Curved arrow properties + style = "Simple, tail_width=0.5, head_width=4, head_length=8" + kw = dict(arrowstyle=style, color="#0077b6") + + colorsForPlot = ["#E5383B","#CD4448","#B55055", + "#9C5D63","#846970","#6C757D"] + plotIndex = ["(a)", "(b)", "(c)", "(d)"] + + # For now load the results file + loadFileName = ("fullModelSensitivity_rateConstant_20210324_0957_c9e8179.npz", # Varying rate constant + "fullModelSensitivity_flowIn_20210407_1656_c9e8179.npz", # Varying flow in + "fullModelSensitivity_volTotal_20210324_1013_c9e8179.npz", # Varying total volume + "fullModelSensitivity_voidFrac_20210324_1006_c9e8179.npz") # Varying void fraction + + # Loop through the two files to get the histogram + for ii in range(len(loadFileName)): + ax = plt.subplot(2,2,ii+1) + # Create the file name with the path to be loaded + simResultsFile = os.path.join('..','simulationResults',loadFileName[ii]); + + # Check if the file with the simultaion results exist + if os.path.exists(simResultsFile): + # Load the file + loadedFileTemp = load(simResultsFile,allow_pickle=True)["outputStruct"] + # Unpack to get the dictionary + loadedFile = loadedFileTemp[()] + + # Prepare the data for plotting + timeSim = loadedFile[0]["timeSim"]; # Time elapsed [s] + adsorbentDensity = loadedFile[0]["inputParameters"][1]; # [kg/m3] + sensorFingerPrint = np.zeros((len(timeSim),len(loadedFile))); + adsorbentVolume = np.zeros((len(loadedFile),1)); + # Loop over to get the sensor finger print and adsorbent volume + for jj in range(len(loadedFile)): + fingerPrintTemp = (loadedFile[jj]["sensorFingerPrint"] + *loadedFile[jj]["inputParameters"][9] + *adsorbentDensity) # Compute true response [g] + ax.plot(timeSim, fingerPrintTemp, + linewidth=1.5,color=colorsForPlot[jj]) + ax.set(xlim = [timeSim[0], 2000], ylim = [0, 0.12]) + ax.locator_params(axis="x", nbins=4) + ax.locator_params(axis="y", nbins=5) + if ii == 2 or ii == 3: + ax.set(xlabel = '$t$ [s]') + if ii == 0 or ii == 2: + ax.set(ylabel = '$m$ [g]') + # Remove ticks + if ii == 0 or ii == 1: + ax.axes.xaxis.set_ticklabels([]) + if ii == 1 or ii == 3: + ax.axes.yaxis.set_ticklabels([]) + # Create labels for the plots + ax.text(50, 0.11, plotIndex[ii], fontsize=10) + ax.text(1600, 0.11, "Array C", fontsize=10, + color = '#0077b6') + if ii == 0: + curvArr = patches.FancyArrowPatch((800, 0.02), (300, 0.06), + connectionstyle="arc3,rad=0.35", **kw) + ax.add_patch(curvArr) + ax.text(300, 0.065, "$k$", fontsize=10, + color = '#0077b6') + ax.text(1240, 0.101, "Varying Kinetics", fontsize=10, + color = '#0077b6') + ax.text(1700, 0.0025, "0.0001", fontsize=8, + color = colorsForPlot[0]) + ax.text(1700, 0.023, "0.0005", fontsize=8, + color = colorsForPlot[1], rotation=10) + ax.text(1700, 0.037, "0.001", fontsize=8, + color = colorsForPlot[2], rotation=8) + ax.text(200, 0.027, "0.005", fontsize=8, + color = colorsForPlot[3], rotation=45) + ax.text(60, 0.042, "0.01", fontsize=8, + color = colorsForPlot[4], rotation=45) + ax.text(1700, 0.056, "10000", fontsize=8, + color = colorsForPlot[5], rotation=0) + + if ii == 1: + curvArr = patches.FancyArrowPatch((800, 0.02), (300, 0.06), + connectionstyle="arc3,rad=0.35", **kw) + ax.add_patch(curvArr) + ax.text(300, 0.065, "$F^\mathregular{in}$", fontsize=10, + color = '#0077b6') + ax.text(1140, 0.101, "Varying Flow Rate", fontsize=10, + color = '#0077b6') + ax.text(1700, 0.005, "0.001", fontsize=8, + color = colorsForPlot[0]) + ax.text(1700, 0.017, "0.005", fontsize=8, + color = colorsForPlot[1], rotation=10) + ax.text(1700, 0.033, "0.01", fontsize=8, + color = colorsForPlot[2], rotation=12) + ax.text(330, 0.023, "0.05", fontsize=8, + color = colorsForPlot[3], rotation=45) + ax.text(230, 0.032, "0.1", fontsize=8, + color = colorsForPlot[4], rotation=60) + ax.text(1800, 0.056, "1", fontsize=8, + color = colorsForPlot[5], rotation=0) + if ii == 2: + curvArr = patches.FancyArrowPatch((800, 0.01), (300, 0.06), + connectionstyle="arc3,rad=0.35", **kw) + ax.add_patch(curvArr) + ax.text(300, 0.065, "$V_\mathregular{T}$", fontsize=10, + color = '#0077b6') + ax.text(1020, 0.101, "Varying Total Volume", fontsize=10, + color = '#0077b6') + ax.text(1820, 0.008, "0.1", fontsize=8, + color = colorsForPlot[2], rotation=0) + ax.text(1820, 0.018, "0.3", fontsize=8, + color = colorsForPlot[3], rotation=0) + ax.text(1820, 0.034, "0.6", fontsize=8, + color = colorsForPlot[4], rotation=0) + ax.text(1820, 0.056, "1.0", fontsize=8, + color = colorsForPlot[5], rotation=0) + if ii == 3: + curvArr = patches.FancyArrowPatch((30, 0.08), (800, 0.015), + connectionstyle="arc3,rad=-0.35", **kw) + ax.add_patch(curvArr) + ax.text(800, 0.035, "$\epsilon$", fontsize=10, + color = '#0077b6') + ax.text(960, 0.101, "Varying Dead Voidage", fontsize=10, + color = '#0077b6') + ax.text(1800, 0.090, "0.10", fontsize=8, + color = colorsForPlot[0]) + ax.text(1800, 0.074, "0.25", fontsize=8, + color = colorsForPlot[1]) + ax.text(1800, 0.056, "0.50", fontsize=8, + color = colorsForPlot[2]) + ax.text(1820, 0.029, "0.75", fontsize=8, + color = colorsForPlot[3]) + ax.text(1820, 0.013, "0.90", fontsize=8, + color = colorsForPlot[4]) + ax.text(1820, 0.003, "0.99", fontsize=8, + color = colorsForPlot[5]) + + # Save the figure + if saveFlag: + # FileName: designVariables__ + saveFileName = "designVariables_" + currentDT + "_" + gitCommitID + saveFileExtension + savePath = os.path.join('..','simulationFigures','simulationManuscript',saveFileName) + # Check if inputResources directory exists or not. If not, create the folder + if not os.path.exists(os.path.join('..','simulationFigures','simulationManuscript')): + os.mkdir(os.path.join('..','simulationFigures','simulationManuscript')) + plt.savefig (savePath) + plt.show() + + +# fun: getSensorSensitiveRegion +# Simulate the sensor array and obtain the region of sensitivity +def getSensorSensitiveRegion(sensorID): + import numpy as np + from simulateSensorArray import simulateSensorArray + from kneed import KneeLocator # To compute the knee/elbow of a curve + + # Total pressure of the gas [Pa] + pressureTotal = np.array([1.e5]); + + # Temperature of the gas [K] + # Can be a vector of temperatures + temperature = np.array([298.15]); + + # Number of molefractions + numMolFrac= 101 + moleFractionRange = np.array([np.linspace(0,1,numMolFrac), 1 - np.linspace(0,1,numMolFrac)]).T + + # Simulate the sensor response for all possible concentrations + arraySimResponse = np.zeros([moleFractionRange.shape[0],sensorID.shape[0]]) + for ii in range(moleFractionRange.shape[0]): + arraySimResponse[ii,:] = simulateSensorArray(sensorID, pressureTotal, + temperature, np.array([moleFractionRange[ii,:]])) + + # Compute the sensitive region for each sensor material in the array + sensitiveRegion = np.zeros([arraySimResponse.shape[1],2]) + # Loop through all the materials + firstDerivative = np.zeros([arraySimResponse.shape[0],arraySimResponse.shape[1]]) + firstDerivativeSimResponse_y1 = np.zeros([moleFractionRange.shape[0],arraySimResponse.shape[1]]) + firstDerivativeSimResponse_y2 = np.zeros([moleFractionRange.shape[0],arraySimResponse.shape[1]]) + for kk in range(arraySimResponse.shape[1]): + firstDerivative[:,kk] = np.gradient(arraySimResponse[:,kk],moleFractionRange[1,0]-moleFractionRange[0,0]) + firstDerivativeSimResponse_y1[:,kk] = np.gradient(moleFractionRange[:,0],arraySimResponse[:,kk]) + firstDerivativeSimResponse_y2[:,kk] = np.gradient(moleFractionRange[:,1],arraySimResponse[:,kk]) + # Get the sign of the first derivative for increasing/decreasing + if all(i >= 0. for i in firstDerivative[:,kk]): + slopeDir = "increasing" + elif all(i < 0. for i in firstDerivative[:,kk]): + slopeDir = "decreasing" + else: + print("Dangerous! I should not be here!!!") + + # Compute the knee/elbow of the curve + kneedle = KneeLocator(moleFractionRange[:,0], arraySimResponse[:,kk], + direction=slopeDir) + elbowPoint = list(kneedle.all_elbows) + + # Obtain coordinates to fill working region + if slopeDir == "increasing": + sensitiveRegion[kk,:] = [0,elbowPoint[0]] + else: + sensitiveRegion[kk,:] = [elbowPoint[0], 1.0] + + # Return the mole fraction, response and sensitive region for each + # material + return moleFractionRange, arraySimResponse, sensitiveRegion + +# fun: concatenateConcEstimate +# Concatenates concentration estimates into a panda dataframe and computes +# the coefficient of variation +def concatenateConcEstimate(loadFileName,sensorText): + import numpy as np + from numpy import load + import os + import pandas as pd + + # Initialize x, y, and type for the plotting + concatenatedX = [] + concatenatedY1 = [] + concatenatedY2 = [] + concatenatedType = [] + + # Loop through the different files to generate the violin plot + for kk in range(len(loadFileName)): + # Initialize x, y, and type for the local loop + xVar = [] + y1Var = [] + y2Var = [] + typeVar = [] + + simResultsFile = os.path.join('..','simulationResults',loadFileName[kk]); + resultOutput = load(simResultsFile)["arrayConcentration"] + moleFrac = load(simResultsFile)["trueMoleFrac"] + + # Loop through all the molefractions + for ii in range(resultOutput.shape[0]): + if resultOutput.shape[2] == 3: + counterInd = -1 + elif resultOutput.shape[2] == 4: + counterInd = 0 + elif resultOutput.shape[2] == 5: + counterInd = 1 + + y1Var = np.concatenate((y1Var,resultOutput[ii,:,counterInd+2])) # y1 + y2Var = np.concatenate((y2Var,resultOutput[ii,:,counterInd+3])) # y2 + xVar = xVar + ([str(moleFrac[ii])] * len(resultOutput[ii,:,counterInd+2])) # x (true mole fraction) + typeVar = typeVar+[sensorText[kk]] * len(resultOutput[ii,:,counterInd+2]) + + # Concatenate all the data to form a data frame with x, y, and type + concatenatedX = concatenatedX + xVar + concatenatedY1 = np.concatenate((concatenatedY1,y1Var)) + concatenatedY2 = np.concatenate((concatenatedY2,y2Var)) + concatenatedType = concatenatedType + typeVar + # Reinitialize all the loaded values to empty variable + simResultsFile = [] + resultOutput = [] + moleFrac = [] + + # Generate panda data frame + # x = molefraction (true) + # y = molefraction (estimated) + # dataType = either sensor id/comparison type + df = pd.DataFrame({'x':concatenatedX, + 'y1':concatenatedY1, + 'y2':concatenatedY2, + 'dataType':concatenatedType}) + + # Compute the mean and standard deviation + meanData = df.groupby(['dataType','x'], as_index=False, sort=False).mean() + stdData = df.groupby(['dataType','x'], as_index=False, sort=False).std() + + # Compute the relative error of the mean (non-negative) + meanErr = stdData.copy() + meanErr['y1'] = abs(meanData['x'].astype(float) - meanData['y1'])/meanData['x'].astype(float) + meanErr['y2'] = abs((1.-meanData['x'].astype(float)) - meanData['y2'])/(1.-meanData['x'].astype(float)) + + # Coefficient of variation + cvData = stdData.copy() + cvData['y1'] = stdData['y1']/meanData['y1'] + cvData['y2'] = stdData['y2']/meanData['y2'] + + # Return the coefficient of variation + return meanErr, cvData \ No newline at end of file diff --git a/plotFunctions/singleColumn.mplstyle b/plotFunctions/singleColumn.mplstyle new file mode 100755 index 0000000..ea278cf --- /dev/null +++ b/plotFunctions/singleColumn.mplstyle @@ -0,0 +1,44 @@ +# matplotlib configuration file for the ERASE project +# https://matplotlib.org/3.3.2/tutorials/introductory/customizing.html + +## Figure property +figure.figsize : 3.3, 3 # width, height in inches +figure.dpi : 600 # dpi +figure.autolayout : true # for labels not being cut out + +## Axes +axes.titlesize : 10 +axes.labelsize : 10 +axes.formatter.limits : -5, 4 + +## Grid +axes.grid : true +grid.color : cccccc +grid.linewidth : 0.5 + +## Lines & Scatter +lines.linewidth : 1.5 +lines.markersize : 4 +scatter.marker: o + +## Ticks +xtick.top : true +ytick.right : true +xtick.labelsize : 8 +ytick.labelsize : 8 + +## Fonts +font.family : arial +font.weight : normal +font.size : 10 + +## Legends +legend.frameon : true +legend.fontsize : 10 +legend.edgecolor : 1 +legend.framealpha : 0.6 + +## Save figure +savefig.dpi : figure +savefig.format : png +savefig.transparent : false \ No newline at end of file diff --git a/screenSensorArray.py b/screenSensorArray.py new file mode 100755 index 0000000..024801c --- /dev/null +++ b/screenSensorArray.py @@ -0,0 +1,128 @@ +############################################################################ +# +# Imperial College London, United Kingdom +# Multifunctional Nanomaterials Laboratory +# +# Project: ERASE +# Year: 2020 +# Python: Python 3.7 +# Authors: Ashwin Kumar Rajagopalan (AK) +# +# Purpose: +# Generates hypothetical sorbents using latin hypercube sampling. The +# sorbents are assumed to exhibit Langmuirian behavior. +# +# Last modified: +# - 2021-02-11, AK: Minor fixes +# - 2021-02-11, AK: Fix for parallel computing +# - 2020-10-29, AK: Add 3 sorbent sensor +# - 2020-10-28, AK: Add auxiliary functions as a module +# - 2020-10-22, AK: Initial creation +# +# Input arguments: +# +# +# Output arguments: +# +# +############################################################################ + +import numpy as np +from numpy import savez +import auxiliaryFunctions +import multiprocessing # For parallel processing +from joblib import Parallel, delayed # For parallel processing +from tqdm import tqdm # To track progress of the loop +from estimateConcentration import estimateConcentration +import os + +# Number of sensors in the array +numSensors = 1 + +# Get the commit ID of the current repository +gitCommitID = auxiliaryFunctions.getCommitID() + +# Find out the total number of cores available for parallel processing +num_cores = multiprocessing.cpu_count() + +# Total number of sensor elements/gases simulated and generated using +# generateHypotheticalAdsorbents.py function +numberOfAdsorbents = 20 +numberOfGases = 2 + +# "True" gas composition that is exposed to the sensor array (0-4) +# Check generateTrueSensorResponse.py for the actual concentrations +moleFracID = 1 + +# Check for number of sensors in the array +if numSensors == 1: + ##### FOR 1 SORBENT SENSOR ARRAY ##### + # Get the current date and time for saving purposes + simulationDT = auxiliaryFunctions.getCurrentDateTime() + + # Loop over all the sorbents for a single material sensor + # Using parallel processing to loop through all the materials + arrayConcentration = np.zeros(numberOfAdsorbents) + arrayConcentration = Parallel(n_jobs=num_cores)(delayed(estimateConcentration) + (numberOfAdsorbents,numberOfGases,moleFracID,[ii]) + for ii in tqdm(range(numberOfAdsorbents))) + + # Convert the output list to a matrix + arrayConcentration = np.array(arrayConcentration) + +elif numSensors == 2: + ##### FOR 2 SORBENT SENSOR ARRAY ##### + # Get the current date and time for saving purposes + simulationDT = auxiliaryFunctions.getCurrentDateTime() + + # Loop over all the sorbents for a single material sensor + # Using parallel processing to loop through all the materials + arrayConcentration = np.zeros(numberOfAdsorbents) + for jj in range(numberOfAdsorbents-1): + arrayConcentrationTemp = Parallel(n_jobs=num_cores)(delayed(estimateConcentration) + (numberOfAdsorbents,numberOfGases,moleFracID,[ii,jj]) + for ii in tqdm(range(jj+1,numberOfAdsorbents))) + # Convert the output list to a matrix + arrayConcentrationTemp = np.array(arrayConcentrationTemp) + if jj == 0: + arrayConcentration = arrayConcentrationTemp + else: + arrayConcentration = np.append(arrayConcentration,arrayConcentrationTemp, axis=0) + +elif numSensors == 3: + ##### FOR 3 SORBENT SENSOR ARRAY ##### + # Get the current date and time for saving purposes + simulationDT = auxiliaryFunctions.getCurrentDateTime() + + # Loop over all the sorbents for a single material sensor + # Using parallel processing to loop through all the materials + arrayConcentration = np.zeros(numberOfAdsorbents) + + for kk in range(numberOfAdsorbents-1): + for jj in range(kk+1,numberOfAdsorbents-1): + arrayConcentrationTemp = Parallel(n_jobs=num_cores)(delayed(estimateConcentration) + (numberOfAdsorbents,numberOfGases,moleFracID,[ii,jj,kk]) + for ii in tqdm(range(jj+1,numberOfAdsorbents))) + # Convert the output list to a matrix + arrayConcentrationTemp = np.array(arrayConcentrationTemp) + if kk == 0 and jj == 1: + arrayConcentration = arrayConcentrationTemp + else: + arrayConcentration = np.append(arrayConcentration,arrayConcentrationTemp, axis=0) + +# Save the array concentration into a native numpy file +# The .npy file is saved in a folder called simulationResults (hardcoded) +filePrefix = "arrayConcentration" +saveFileName = filePrefix + "_" + simulationDT + "_" + gitCommitID; +savePath = os.path.join('simulationResults',saveFileName) + +# Check if inputResources directory exists or not. If not, create the folder +if not os.path.exists('simulationResults'): + os.mkdir('simulationResults') + +# Save the array ceoncentration obtained from estimateConcentration +savez (savePath, arrayConcentration = arrayConcentration, # Estimated Concentration + numberOfAdsorbents = numberOfAdsorbents, # Estimated response + numberOfGases = numberOfGases, # Flag to gases to be sensed + numSensors = numSensors, # Number of sensors in the array + moleFracID = moleFracID) # Mole fraction \ No newline at end of file diff --git a/sensitivityAnalysis.py b/sensitivityAnalysis.py new file mode 100755 index 0000000..a906aa3 --- /dev/null +++ b/sensitivityAnalysis.py @@ -0,0 +1,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) \ No newline at end of file diff --git a/sensorFullModelWrapper.py b/sensorFullModelWrapper.py new file mode 100644 index 0000000..2d9c16b --- /dev/null +++ b/sensorFullModelWrapper.py @@ -0,0 +1,103 @@ +############################################################################ +# +# Imperial College London, United Kingdom +# Multifunctional Nanomaterials Laboratory +# +# Project: ERASE +# Year: 2021 +# Python: Python 3.7 +# Authors: Ashwin Kumar Rajagopalan (AK) +# +# Purpose: +# Wrapper function for sensor full model to check the impact of different +# variables +# +# Last modified: +# - 2021-03-23, AK: Minor fixes +# - 2021-02-03, AK: Change output plot response to absolute values +# - 2021-01-20, AK: Initial creation +# +# Input arguments: +# +# +# Output arguments: +# +# +############################################################################ + +import auxiliaryFunctions +from simulateFullModel import simulateFullModel +from tqdm import tqdm # To track progress of the loop +import numpy as np +import os +from numpy import savez +import socket +import matplotlib.pyplot as plt + +# Save and plot settings +saveFlag = True +saveFileExtension = ".png" +saveFlag = "rateConstant" +colorsForPlot = ["E5383B","CD4448","B55055","9C5D63","846970","6C757D"] + +# Get the commit ID of the current repository +gitCommitID = auxiliaryFunctions.getCommitID() + +# Get the current date and time for saving purposes +currentDT = auxiliaryFunctions.getCurrentDateTime() + +# Define the variable to be looped +# This has to be a be a tuple. For on condition write the values followed by a +# comma to say its a tuple +volTotal = 10e-7 +voidFrac = 0.5 +loopVariable = ([0.0001,0.0001,0.0001], + [0.0005,0.0005,0.0005], + [0.001,0.001,0.001], + [0.005,0.005,0.005], + [0.01,0.01,0.01], + [10000.0,10000.0,10000.0]) + +# Define a dictionary +outputStruct = {} + +# Loop over all the individual elements of the loop variable +for ii in tqdm(range(len(loopVariable))): + # Call the full model with a given rate constant + timeSim, _ , sensorFingerPrint, inputParameters = simulateFullModel(volTotal = volTotal, + voidFrac = voidFrac, + rateConstant = loopVariable[ii]) + outputStruct[ii] = {'timeSim':timeSim, + 'sensorFingerPrint':sensorFingerPrint, + 'inputParameters':inputParameters} + +# Save the array concentration into a native numpy file +# The .npz file is saved in a folder called simulationResults (hardcoded) +# Save the figure +if saveFlag: + filePrefix = "fullModelSensitivity" + saveFileName = filePrefix + "_" + saveFlag + "_" + currentDT + "_" + gitCommitID; + savePath = os.path.join('simulationResults',saveFileName) + savez (savePath, outputStruct = outputStruct, # True response + hostName = socket.gethostname()) # Hostname of the computer + +# Plot the sensor finger print +os.chdir("plotFunctions") +plt.style.use('singleColumn.mplstyle') # Custom matplotlib style file +fig = plt.figure +ax = plt.subplot(1,1,1) +for ii in range(len(loopVariable)): + timeTemp = outputStruct[ii]["timeSim"] + fingerPrintTemp = (outputStruct[ii]["sensorFingerPrint"] + *outputStruct[ii]["inputParameters"][1] + *outputStruct[ii]["inputParameters"][9]) # Compute true response [g] + ax.plot(timeTemp, fingerPrintTemp, + linewidth=1.5,color="#"+colorsForPlot[ii], + label = str(loopVariable[ii])) +ax.set(xlabel='$t$ [s]', + ylabel='$m_i$ [g]', + xlim = [timeSim[0], timeSim[-1]], ylim = [0, None]) +ax.locator_params(axis="x", nbins=4) +ax.locator_params(axis="y", nbins=4) +plt.show() +os.chdir("..") \ No newline at end of file diff --git a/simulateFullModel.py b/simulateFullModel.py new file mode 100644 index 0000000..de7394c --- /dev/null +++ b/simulateFullModel.py @@ -0,0 +1,414 @@ +############################################################################ +# +# Imperial College London, United Kingdom +# Multifunctional Nanomaterials Laboratory +# +# Project: ERASE +# Year: 2021 +# Python: Python 3.7 +# Authors: Ashwin Kumar Rajagopalan (AK) +# +# Purpose: +# Simulates the sensor chamber as a CSTR incorporating kinetic effects +# +# Last modified: +# - 2021-03-20, AK: Bug fix for flow rate calculator +# - 2021-02-19, AK: Add relative tolerances for ode solver +# - 2021-02-03, AK: Add total volume and void fraction +# - 2021-02-02, AK: Add flow rate to output +# - 2021-01-30, AK: Add constant pressure model +# - 2021-01-27, AK: Add volSorbent and volGas to inputs +# - 2021-01-25, AK: Change the time interval definition +# - 2021-01-21, AK: Cosmetic changes +# - 2021-01-20, AK: Change to output time and plot function +# - 2021-01-20, AK: Cosmetic changes +# - 2021-01-19, AK: Initial creation +# +# Input arguments: +# +# +# Output arguments: +# +# +############################################################################ + +def simulateFullModel(**kwargs): + import numpy as np + from scipy.integrate import solve_ivp + from simulateSensorArray import simulateSensorArray + import auxiliaryFunctions + + # Plot flag + plotFlag = False + + # Get the commit ID of the current repository + gitCommitID = auxiliaryFunctions.getCommitID() + + # Get the current date and time for saving purposes + currentDT = auxiliaryFunctions.getCurrentDateTime() + + # Model flag (constant pressure or constant flow rate) + # Default is constant pressure + if 'modelConstF' in kwargs: + modelConstF = kwargs["modelConstF"] + else: + modelConstF = False + + # Sensor ID + if 'sensorID' in kwargs: + sensorID = np.array([kwargs["sensorID"]]) + else: + sensorID = np.array([6]) + + # Kinetic rate constants [/s] + if 'rateConstant' in kwargs: + rateConstant = np.array(kwargs["rateConstant"]) + else: + rateConstant = np.array([.01,.01,.01]) + + # Feed flow rate [m3/s] + if 'flowIn' in kwargs: + flowIn = np.array(kwargs["flowIn"]) + else: + flowIn = np.array([5e-7]) + + # Feed Mole Fraction [-] + if 'feedMoleFrac' in kwargs: + feedMoleFrac = np.array(kwargs["feedMoleFrac"]) + else: + feedMoleFrac = np.array([1.,0.,0.]) + + # Initial Gas Mole Fraction [-] + if 'initMoleFrac' in kwargs: + initMoleFrac = np.array(kwargs["initMoleFrac"]) + else: + # Equilibrium process + initMoleFrac = np.array([0.,0.,1.]) + + # Time span for integration [tuple with t0 and tf] + if 'timeInt' in kwargs: + timeInt = kwargs["timeInt"] + else: + timeInt = (0.0,2000) + + # Time step for output printing + if 'timeStep' in kwargs: + timeStep = kwargs["timeStep"] + else: + timeStep = 5. + + # Volume of sorbent material [m3] + if 'volSorbent' in kwargs: + volSorbent = kwargs["volSorbent"] + else: + volSorbent = 5e-7 + + # Volume of gas chamber (dead volume) [m3] + if 'volGas' in kwargs: + volGas = kwargs["volGas"] + else: + volGas = 5e-7 + + # Total volume of the system [m3] + if 'volTotal' in kwargs: + volTotal = kwargs["volTotal"] + if 'voidFrac' in kwargs: + voidFrac = kwargs["voidFrac"] + else: + raise Exception("You should provide void fraction if you provide total volume!") + volGas = voidFrac * volTotal # Volume of gas chamber (dead volume) [m3] + volSorbent = (1-voidFrac) * volTotal # Volume of solid sorbent [m3] + + if (len(feedMoleFrac) != len(initMoleFrac) + or len(feedMoleFrac) != len(rateConstant)): + raise Exception("The dimensions of the mole fraction or rate constant and the number of gases in the adsorbent is not consistent!") + else: + numberOfGases = len(feedMoleFrac) + + # Total pressure of the gas [Pa] + if 'pressureTotal' in kwargs: + pressureTotal = np.array(kwargs["pressureTotal"]); + else: + pressureTotal = np.array([1.e5]); + + # Temperature of the gas [K] + # Can be a vector of temperatures + if 'temperature' in kwargs: + temperature = np.array(kwargs["temperature"]); + else: + temperature = np.array([298.15]); + + # Compute the initial sensor loading [mol/m3] @ initMoleFrac + sensorLoadingPerGasVol, adsorbentDensity, molecularWeight = simulateSensorArray(sensorID, pressureTotal, + temperature, np.array([initMoleFrac]), + fullModel = True) + + # Prepare tuple of input parameters for the ode solver + inputParameters = (sensorID, adsorbentDensity, rateConstant, numberOfGases, flowIn, + feedMoleFrac, initMoleFrac, pressureTotal, temperature, volSorbent, volGas) + + # Solve the system of ordinary differential equations + # Stiff solver used for the problem: BDF or Radau + # The output is print out every 5 s + # Solves the model assuming constant flow rate at the inlet and outlet + if modelConstF: + # Prepare initial conditions vector + initialConditions = np.zeros([2*numberOfGases]) + initialConditions[0:numberOfGases-1] = initMoleFrac[0:numberOfGases-1] # Gas mole fraction + initialConditions[numberOfGases-1:2*numberOfGases-1] = sensorLoadingPerGasVol # Initial Loading + initialConditions[2*numberOfGases-1] = pressureTotal # Outlet pressure the same as inlet pressure + outputSol = solve_ivp(solveSorptionEquationConstF, timeInt, initialConditions, + method='Radau', t_eval = np.arange(timeInt[0],timeInt[1],timeStep), + rtol = 1e-6, args = inputParameters) + + # Flow out vector in output + flowOutVec = flowIn * np.ones(len(outputSol.t)) # Constant flow rate + + # Parse out the output matrix and add flow rate + resultMat = np.row_stack((outputSol.y,flowOutVec)) + + # Solves the model assuming constant/negligible pressure across the sensor + else: + # Prepare initial conditions vector + initialConditions = np.zeros([2*numberOfGases-1]) + initialConditions[0:numberOfGases-1] = initMoleFrac[0:numberOfGases-1] # Gas mole fraction + initialConditions[numberOfGases-1:2*numberOfGases-1] = sensorLoadingPerGasVol # Initial Loading + outputSol = solve_ivp(solveSorptionEquationConstP, timeInt, initialConditions, + method='Radau', t_eval = np.arange(timeInt[0],timeInt[1],timeStep), + rtol = 1e-6, args = inputParameters) + + # Presure vector in output + pressureVec = pressureTotal * np.ones(len(outputSol.t)) # Constant pressure + + # This is needed to make sure constant pressure model functions well + # If the time of integration and the time step is equal then the code + # will fail because of the flow rate determination step which requires + # a gradient - very first step flow out and in are considered to be the + # same if only one element present + if len(outputSol.t) == 1 : + flowOut = flowIn + else: + # Compute the outlet flow rate + dqdt = np.gradient(outputSol.y[numberOfGases-1:2*numberOfGases-1,:], + outputSol.t, axis=1) # Compute gradient of loading + sum_dqdt = np.sum(dqdt, axis=0) # Time resolved sum of gradient + flowOut = flowIn - ((volSorbent*(8.314*temperature)/pressureTotal)*(sum_dqdt)) + + # Parse out the output matrix and add flow rate + resultMat = np.row_stack((outputSol.y,pressureVec,flowOut)) + + # Parse out the time + timeSim = outputSol.t + + # Compute the time resolved sensor response + sensorFingerPrint = np.zeros([len(timeSim)]) + for ii in range(len(timeSim)): + loadingTemp = resultMat[numberOfGases-1:2*numberOfGases-1,ii] + sensorFingerPrint[ii] = np.dot(loadingTemp,molecularWeight)/adsorbentDensity + + # Call the plotting function + if plotFlag: + plotFullModelResult(timeSim, resultMat, sensorFingerPrint, inputParameters, + gitCommitID, currentDT, modelConstF) + + # Return time and the output matrix + return timeSim, resultMat, sensorFingerPrint, inputParameters + +# func: solveSorptionEquationConstF - Constant flow rate model +# Solves the system of ODEs to evaluate the gas composition, loadings, and pressure +def solveSorptionEquationConstF(t, f, *inputParameters): + import numpy as np + from simulateSensorArray import simulateSensorArray + + # Gas constant + Rg = 8.314; # [J/mol K] + + # Unpack the tuple of input parameters used to solve equations + sensorID, _ , rateConstant, numberOfGases, flowIn, feedMoleFrac, _ , pressureTotal, temperature, volSorbent, volGas = inputParameters + + # Initialize the derivatives to zero + df = np.zeros([2*numberOfGases]) + + # Compute the equilbirium loading at the current gas composition + currentGasComposition = np.concatenate((f[0:numberOfGases-1], + np.array([1.-np.sum(f[0:numberOfGases-1])]))) + sensorLoadingPerGasVol, _ , _ = simulateSensorArray(sensorID, f[2*numberOfGases-1], + temperature, np.array([currentGasComposition]), + fullModel = True) + + # Linear driving force model (derivative of solid phase loadings) + df[numberOfGases-1:2*numberOfGases-1] = np.multiply(rateConstant,(sensorLoadingPerGasVol-f[numberOfGases-1:2*numberOfGases-1])) + + # Total mass balance + # Assumes constant flow rate, so pressure evalauted + term1 = 1/volGas + term2 = ((flowIn*pressureTotal) - (flowIn*f[2*numberOfGases-1]) + - ((volSorbent*(Rg*temperature))*(np.sum(df[numberOfGases-1:2*numberOfGases-1])))) + df[2*numberOfGases-1] = term1*term2 + + # Component mass balance + term1 = 1/volGas + for ii in range(numberOfGases-1): + term2 = (flowIn*(pressureTotal*feedMoleFrac[ii] - f[2*numberOfGases-1]*f[ii]) + - (volSorbent*(Rg*temperature))*df[ii+numberOfGases-1]) + df[ii] = (term1*term2 - f[ii]*df[2*numberOfGases-1])/f[2*numberOfGases-1] + + # Return the derivatives for the solver + return df + +# func: solveSorptionEquationConstP - Constant pressure model +# Solves the system of ODEs to evaluate the gas composition and loadings +def solveSorptionEquationConstP(t, f, *inputParameters): + import numpy as np + from simulateSensorArray import simulateSensorArray + + # Gas constant + Rg = 8.314; # [J/mol K] + + # Unpack the tuple of input parameters used to solve equations + sensorID, _ , rateConstant, numberOfGases, flowIn, feedMoleFrac, _ , pressureTotal, temperature, volSorbent, volGas = inputParameters + + # Initialize the derivatives to zero + df = np.zeros([2*numberOfGases-1]) + + # Compute the equilbirium loading at the current gas composition + currentGasComposition = np.concatenate((f[0:numberOfGases-1], + np.array([1.-np.sum(f[0:numberOfGases-1])]))) + sensorLoadingPerGasVol, _ , _ = simulateSensorArray(sensorID, pressureTotal, + temperature, np.array([currentGasComposition]), + fullModel = True) + + # Linear driving force model (derivative of solid phase loadings) + df[numberOfGases-1:2*numberOfGases-1] = np.multiply(rateConstant,(sensorLoadingPerGasVol-f[numberOfGases-1:2*numberOfGases-1])) + + # Total mass balance + # Assumes constant pressure, so flow rate evalauted + flowOut = flowIn - ((volSorbent*(Rg*temperature)/pressureTotal)*(np.sum(df[numberOfGases-1:2*numberOfGases-1]))) + + # Component mass balance + term1 = 1/volGas + for ii in range(numberOfGases-1): + term2 = ((flowIn*feedMoleFrac[ii] - flowOut*f[ii]) + - (volSorbent*(Rg*temperature)/pressureTotal)*df[ii+numberOfGases-1]) + df[ii] = term1*term2 + + # Return the derivatives for the solver + return df + +# func: plotFullModelResult +# Plots the model output for the conditions simulated locally +def plotFullModelResult(timeSim, resultMat, sensorFingerPrint, inputParameters, + gitCommitID, currentDT, modelConstF): + import numpy as np + import os + import matplotlib.pyplot as plt + + # Save settings + saveFlag = False + saveFileExtension = ".png" + + # Unpack the tuple of input parameters used to solve equations + sensorID, _ , _ , _ , flowIn, _ , _ , _ , temperature, _ , _ = inputParameters + + os.chdir("plotFunctions") + if resultMat.shape[0] == 7: + # Plot the solid phase compositions + plt.style.use('doubleColumn.mplstyle') # Custom matplotlib style file + fig = plt.figure + ax = plt.subplot(1,3,1) + ax.plot(timeSim, resultMat[2,:], + linewidth=1.5,color='r') + ax.set(xlabel='$t$ [s]', + ylabel='$q_1$ [mol m$^{\mathregular{-3}}$]', + xlim = [timeSim[0], timeSim[-1]], ylim = [0, 1.1*np.max(resultMat[2,:])]) + + ax = plt.subplot(1,3,2) + ax.plot(timeSim, resultMat[3,:], + linewidth=1.5,color='b') + ax.set(xlabel='$t$ [s]', + ylabel='$q_2$ [mol m$^{\mathregular{-3}}$]', + xlim = [timeSim[0], timeSim[-1]], ylim = [0, 1.1*np.max(resultMat[3,:])]) + + ax = plt.subplot(1,3,3) + ax.plot(timeSim, resultMat[4,:], + linewidth=1.5,color='g') + ax.set(xlabel='$t$ [s]', + ylabel='$q_3$ [mol m$^{\mathregular{-3}}$]', + xlim = [timeSim[0], timeSim[-1]], ylim = [0, 1.1*np.max(resultMat[4,:])]) + + # Save the figure + if saveFlag: + # FileName: solidLoadingFM___ + saveFileName = "solidLoadingFM_" + str(sensorID) + "_" + currentDT + "_" + gitCommitID + saveFileExtension + savePath = os.path.join('..','simulationFigures',saveFileName.replace('[','').replace(']','')) + # Check if inputResources directory exists or not. If not, create the folder + if not os.path.exists(os.path.join('..','simulationFigures')): + os.mkdir(os.path.join('..','simulationFigures')) + plt.savefig (savePath) + + plt.show() + + # Plot the pressure drop, the flow rate, and the mole fraction + plt.style.use('doubleColumn.mplstyle') # Custom matplotlib style file + fig = plt.figure + ax = plt.subplot(1,3,1) + ax.plot(timeSim, resultMat[5,:], + linewidth=1.5,color='r') + ax.set_xlabel('$t$ [s]') + ax.set_ylabel('$P$ [Pa]', color='r') + ax.tick_params(axis='y', labelcolor='r') + ax.set(xlim = [timeSim[0], timeSim[-1]], + ylim = [0, 1.1*np.max(resultMat[5,:])]) + + ax2 = plt.twinx() + ax2.plot(timeSim, resultMat[6,:], + linewidth=1.5,color='b') + ax2.set_ylabel('$F$ [m$^{\mathregular{3}}$ s$^{\mathregular{-1}}$]', color='b') + ax2.tick_params(axis='y', labelcolor='b') + ax2.set(xlim = [timeSim[0], timeSim[-1]], + ylim = [0, 1.1*np.max(resultMat[6,:])]) + + ax = plt.subplot(1,3,2) + ax.plot(timeSim, resultMat[5,:]*resultMat[6,:]/(temperature*8.314), + linewidth=1.5,color='k') + ax.plot(timeSim, resultMat[5,:]*resultMat[6,:]*resultMat[0,:]/(temperature*8.314), + linewidth=1.5,color='r') + ax.plot(timeSim, resultMat[5,:]*resultMat[6,:]*resultMat[1,:]/(temperature*8.314), + linewidth=1.5,color='b') + ax.plot(timeSim, resultMat[5,:]*resultMat[6,:]*(1-resultMat[0,:]-resultMat[1,:])/(temperature*8.314), + linewidth=1.5,color='g') + ax.set(xlabel='$t$ [s]', + ylabel='$Q$ [mol s$^{\mathregular{-1}}$]', + xlim = [timeSim[0], timeSim[-1]], ylim = [0, 1.1*np.max(resultMat[5,:])*np.max(resultMat[6,:]) + /temperature/8.314]) + + ax = plt.subplot(1,3,3) + ax.plot(timeSim, resultMat[0,:],linewidth=1.5,color='r') + ax.plot(timeSim, resultMat[1,:],linewidth=1.5,color='b') + ax.plot(timeSim, 1-resultMat[0,:]-resultMat[1,:],linewidth=1.5,color='g') + ax.set(xlabel='$t$ [s]', + ylabel='$y$ [-]', + xlim = [timeSim[0], timeSim[-1]], ylim = [0, 1.]) + plt.show() + + # Plot the sensor finger print + plt.style.use('singleColumn.mplstyle') # Custom matplotlib style file + fig = plt.figure + ax = plt.subplot(1,1,1) + ax.plot(timeSim, sensorFingerPrint, + linewidth=1.5,color='k') + ax.set(xlabel='$t$ [s]', + ylabel='$m_i$ [g kg$^{\mathregular{-1}}$]', + xlim = [timeSim[0], timeSim[-1]], ylim = [0, 1.1*np.max(sensorFingerPrint)]) + # Save the figure + if saveFlag: + # FileName: SensorFingerPrintFM___ + saveFileName = "SensorFingerPrintFM_" + str(sensorID) + "_" + currentDT + "_" + gitCommitID + saveFileExtension + savePath = os.path.join('..','simulationFigures',saveFileName.replace('[','').replace(']','')) + # Check if inputResources directory exists or not. If not, create the folder + if not os.path.exists(os.path.join('..','simulationFigures')): + os.mkdir(os.path.join('..','simulationFigures')) + plt.savefig (savePath) + plt.show() + + os.chdir("..") \ No newline at end of file diff --git a/simulateSSL.py b/simulateSSL.py new file mode 100755 index 0000000..8861751 --- /dev/null +++ b/simulateSSL.py @@ -0,0 +1,69 @@ +############################################################################ +# +# Imperial College London, United Kingdom +# Multifunctional Nanomaterials Laboratory +# +# Project: ERASE +# Year: 2020 +# Python: Python 3.7 +# Authors: Ashwin Kumar Rajagopalan (AK) +# +# Purpose: +# Generates the single site Langmuir isotherms for different temperatures +# and concentrations. +# +# Last modified: +# - 2020-10-23, AK: Replace the noncompetitive to competitive isotherm +# - 2020-10-22, AK: Minor cosmetic changes +# - 2020-10-19, AK: Incorporate sorbent density and minor improvements +# - 2020-10-19, AK: Change functionality to work with a single sorbent input +# - 2020-10-16, AK: Initial creation +# +# Input arguments: +# +# +# Output arguments: +# +# +############################################################################ + +def simulateSSL(adsorbentIsotherm, adsorbentDensity, pressureTotal, + temperature, moleFraction): + import numpy as np + + # Gas constant + Rg = 8.314; # [J/mol K] + + # Get the number of gases + numberOfGases = moleFraction.shape[1] + + if moleFraction.shape[1] != adsorbentIsotherm.shape[0]: + raise Exception("The dimensions of the mole fraction and the number of gases in the adsorbent is not consistent!") + + # Initialize the equilibriumLoadings array with zeros + equilibriumLoadings = np.zeros((moleFraction.shape[0], + temperature.shape[0],numberOfGases)) + + # Generate the isotherm + for ii in range(numberOfGases): + for jj in range(moleFraction.shape[0]): + for kk in range(temperature.shape[0]): + # Parse out the saturation capacity, equilibrium constant, and + # heat of adsorption + qsat = adsorbentIsotherm[ii,0] + b0 = adsorbentIsotherm[:,1] + delH = adsorbentIsotherm[:,2] + b = np.multiply(b0,np.exp(-delH/(Rg*temperature[kk]))) + # Compute the concentraiton + conc = pressureTotal*moleFraction[jj,:]/(Rg*temperature[kk]) + # Compute sum(b*c) for the multicomponet competitive eqbm (den.) + bStarConc = np.multiply(b,conc) + sumbStarConc = np.sum(bStarConc) + # Compute the numerator and denominator for SSL + loadingNum = adsorbentDensity*qsat*b[ii]*conc[ii] + loadingDen = 1 + sumbStarConc + # Compute the equilibrium loading + equilibriumLoadings[jj,kk,ii] = loadingNum/loadingDen # [mol/m3] + + # Return the equilibrium loadings + return equilibriumLoadings \ No newline at end of file diff --git a/simulateSensorArray.py b/simulateSensorArray.py new file mode 100755 index 0000000..a3d9ed1 --- /dev/null +++ b/simulateSensorArray.py @@ -0,0 +1,92 @@ +############################################################################ +# +# Imperial College London, United Kingdom +# Multifunctional Nanomaterials Laboratory +# +# Project: ERASE +# Year: 2020 +# Python: Python 3.7 +# Authors: Ashwin Kumar Rajagopalan (AK) +# +# Purpose: +# Generates the response of a sensor array as change in mass due to gas +# sorption at a given pressure, temperature, and mole fraction for +# n sorbents. +# +# Last modified: +# - 2021-01-20, AK: Structure change and add material file for full model +# - 2021-01-19, AK: Add flag for full model +# - 2020-10-30, AK: Fix to find number of gases +# - 2020-10-22, AK: Add two/three gases +# - 2020-10-21, AK: Cosmetic changes and make it a function +# - 2020-10-20, AK: Obtain sensor array finger print +# - 2020-10-19, AK: Initial creation +# +# Input arguments: +# +# +# Output arguments: +# +# +############################################################################ + +def simulateSensorArray(sensorID, pressureTotal, temperature, moleFraction, **kwargs): + import numpy as np + from numpy import load + import os + from simulateSSL import simulateSSL + + # Flag to check if simulation full model or not + if 'fullModel' in kwargs: + if kwargs["fullModel"]: + flagFullModel = True + else: + flagFullModel = False + else: + flagFullModel = False + + # Load a given adsorbent isotherm material file based on full model flag + if flagFullModel: + if moleFraction.shape[1] == 3: + loadFileName = "isothermParameters_20210120_1722_fb57143.npz" # Two gases + Inert + elif moleFraction.shape[1] == 4: + loadFileName = "isothermParameters_20210120_1724_fb57143.npz" # Three gases + Inert + else: + if moleFraction.shape[1] == 2: + loadFileName = "isothermParameters_20201020_1756_5f263af.npz" # Two gases + elif moleFraction.shape[1] == 3: + loadFileName = "isothermParameters_20201022_1056_782efa3.npz" # Three gases + hypoAdsorbentFile = os.path.join('inputResources',loadFileName); + + # Check if the file with the adsorbent properties exist + if os.path.exists(hypoAdsorbentFile): + loadedFileContent = load(hypoAdsorbentFile) + adsorbentIsothermTemp = loadedFileContent['adsIsotherm'] + adsorbentDensityTemp = loadedFileContent['adsDensity'] + molecularWeight = loadedFileContent['molWeight'] + else: + errorString = "Adsorbent property file " + hypoAdsorbentFile + " does not exist." + raise Exception(errorString) + + # Get the equilibrium loading for all the sensors for each gas + # This is a [nxg] matrix where n is the number of sensors and g the number + # of gases + sensorLoadingPerGasVol = np.zeros((sensorID.shape[0],moleFraction.shape[1])) # [mol/m3] + sensorLoadingPerGasMass = np.zeros((sensorID.shape[0],moleFraction.shape[1])) # [mol/kg] + for ii in range(sensorID.shape[0]): + adsorbentID = sensorID[ii] + adsorbentIsotherm = adsorbentIsothermTemp[:,:,adsorbentID] + adsorbentDensity = adsorbentDensityTemp[adsorbentID] + equilibriumLoadings = simulateSSL(adsorbentIsotherm,adsorbentDensity, + pressureTotal,temperature,moleFraction) # [mol/m3] + sensorLoadingPerGasVol[ii,:] = equilibriumLoadings[0,0,:] # [mol/m3] + sensorLoadingPerGasMass[ii,:] = equilibriumLoadings[0,0,:]/adsorbentDensity # [mol/kg] + + # Obtain the sensor finger print # [g of total gas adsorbed/kg of sorbent] + sensorFingerPrint = np.dot(sensorLoadingPerGasMass,molecularWeight) # [g/kg] + + # Flag to check if simulation full model or not + if flagFullModel: + return sensorLoadingPerGasVol, adsorbentDensity, molecularWeight + else: + return sensorFingerPrint \ No newline at end of file