From f896264676a5ecd0c3392d9cd4e5cf0e66f278ea Mon Sep 17 00:00:00 2001 From: Logan Ward Date: Mon, 1 Aug 2022 11:48:41 -0400 Subject: [PATCH 01/29] Updated environment file - Removed duplicate numba channel - Relaxed version for ipopt - Added a newline at the end of the file --- frhodo_environment.yml | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/frhodo_environment.yml b/frhodo_environment.yml index bcb685d..daecfb0 100644 --- a/frhodo_environment.yml +++ b/frhodo_environment.yml @@ -2,12 +2,11 @@ name: frhodo channels: - numba - cantera/label/dev - - numba - conda-forge - defaults dependencies: - cantera=2.6.0a1 - - ipopt=3.14.4 + - ipopt=3.14.* - matplotlib=3.5.1 - nlopt=2.7.0 - numba=0.54.1 @@ -22,4 +21,4 @@ dependencies: - tabulate=0.8.9 - pip: - dtcwt==0.12.0 - - rbfopt==4.2.2 \ No newline at end of file + - rbfopt==4.2.2 From 05bd7899a387054cd46e089b0309b4577088c483 Mon Sep 17 00:00:00 2001 From: Logan Ward Date: Mon, 1 Aug 2022 14:33:57 -0400 Subject: [PATCH 02/29] Refactored to create an importable package - Running by `python src/main.py` is now deprecated. We can have importing as a module or running as a script, both is annoying - Added a `frhodo` command to launch the GUI - Had to change a lot of imports to make it happen - Moved the version to a special file --- {src => frhodo}/CiteSoftLocal.py | 0 {src => frhodo}/UI/error_window.ui | 0 {src => frhodo}/UI/graphics/arrowdown.png | Bin {src => frhodo}/UI/graphics/check_icon.png | Bin {src => frhodo}/UI/graphics/main_icon.png | Bin {src => frhodo}/UI/graphics/x_icon.png | Bin {src => frhodo}/UI/main_window.ui | 0 {src => frhodo}/UI/save_dialog.ui | 0 {src/calculate => frhodo}/__init__.py | 0 {src => frhodo}/appdirs.py | 0 {src => frhodo}/bonmin/bonmin-linux64/bonmin | Bin .../bonmin/bonmin-linux64/coin-license.txt | 0 {src => frhodo}/bonmin/bonmin-osx/bonmin | Bin .../bonmin/bonmin-osx/coin-license.txt | 0 .../bonmin/bonmin-win64/bonmin.exe | Bin .../bonmin/bonmin-win64/coin-license.txt | 0 .../bonmin/bonmin-win64/libipoptfort.dll | Bin .../optimize => frhodo/calculate}/__init__.py | 0 {src => frhodo}/calculate/convert_units.py | 0 {src => frhodo}/calculate/integrate.py | 0 {src => frhodo}/calculate/mech_fcns.py | 5 +- .../optimize/CheKiPEUQ_from_Frhodo.py | 11 ++--- .../optimize/CheKiPEUQ_integration_notes.txt | 0 .../CombinationGeneratorModule.py | 0 .../CheKiPEUQ_local/InverseProblem.py | 13 ++--- .../optimize/CheKiPEUQ_local/LICENSE.txt | 0 .../optimize/CheKiPEUQ_local/MANUAL.txt | 0 .../optimize/CheKiPEUQ_local/UserInput.py | 0 .../optimize/CheKiPEUQ_local/__init__.py | 0 .../mumpce_custom_plotting_example.py | 0 .../nestedObjectsFunctionsLocal.py | 0 .../CheKiPEUQ_local/parallel_processing.py | 0 .../CheKiPEUQ_local/plotting_functions.py | 0 .../optimize/CheKiPEUQ_local/test_pytest.py | 0 .../calculate/optimize}/__init__.py | 0 .../calculate/optimize/fit_coeffs.py | 7 +-- .../calculate/optimize/fit_coeffs_pygmo.py | 9 ++-- {src => frhodo}/calculate/optimize/fit_fcn.py | 10 ++-- .../calculate/optimize/mech_optimize.py | 8 ++-- .../calculate/optimize/misc_fcns.py | 4 +- .../calculate/optimize/optimize_worker.py | 7 +-- {src => frhodo}/calculate/reactors.py | 4 +- {src => frhodo}/calculate/shock_fcns.py | 0 {src => frhodo}/calculate/smooth_data.py | 0 {src => frhodo}/ck2yaml.py | 0 {src => frhodo}/colors.py | 0 {src => frhodo}/config_io.py | 0 .../data/loss_partition_fcn_tck_spline.dat | 0 {src => frhodo}/error_window.py | 0 {src => frhodo}/help_menu.py | 0 .../ipopt/ipopt-linux64/coin-license.txt | 0 {src => frhodo}/ipopt/ipopt-linux64/ipopt | Bin .../ipopt/ipopt-osx/coin-license.txt | 0 {src => frhodo}/ipopt/ipopt-osx/ipopt | Bin .../ipopt/ipopt-win64/coin-license.txt | 0 {src => frhodo}/ipopt/ipopt-win64/ipopt.exe | Bin .../ipopt/ipopt-win64/libipoptfort.dll | Bin {src => frhodo}/main.py | 45 +++++++++++------- {src => frhodo}/mech_widget.py | 8 ++-- {src => frhodo}/misc_widget.py | 2 +- {src => frhodo}/options_panel_widgets.py | 15 +++--- frhodo/plot/__init__.py | 0 {src => frhodo}/plot/base_plot.py | 6 +-- .../plot/custom_mpl_ticker_formatter.py | 0 {src => frhodo}/plot/custom_mplscale.py | 2 +- {src => frhodo}/plot/draggable.py | 0 {src => frhodo}/plot/optimization_plot.py | 4 +- {src => frhodo}/plot/plot_main.py | 2 +- {src => frhodo}/plot/plot_widget.py | 2 +- {src => frhodo}/plot/raw_signal_plot.py | 4 +- {src => frhodo}/plot/signal_plot.py | 6 +-- {src => frhodo}/plot/sim_explorer_plot.py | 4 +- {src => frhodo}/save_output.py | 4 +- {src => frhodo}/save_widget.py | 0 {src => frhodo}/series_viewer_widget.py | 2 +- {src => frhodo}/settings.py | 4 +- {src => frhodo}/sim_explorer_widget.py | 7 ++- {src => frhodo}/soln2ck.py | 0 {src => frhodo}/thermo_widget.py | 8 +--- frhodo/version.py | 4 ++ setup.py | 20 ++++++++ 81 files changed, 130 insertions(+), 97 deletions(-) rename {src => frhodo}/CiteSoftLocal.py (100%) rename {src => frhodo}/UI/error_window.ui (100%) rename {src => frhodo}/UI/graphics/arrowdown.png (100%) rename {src => frhodo}/UI/graphics/check_icon.png (100%) rename {src => frhodo}/UI/graphics/main_icon.png (100%) rename {src => frhodo}/UI/graphics/x_icon.png (100%) rename {src => frhodo}/UI/main_window.ui (100%) rename {src => frhodo}/UI/save_dialog.ui (100%) rename {src/calculate => frhodo}/__init__.py (100%) rename {src => frhodo}/appdirs.py (100%) rename {src => frhodo}/bonmin/bonmin-linux64/bonmin (100%) rename {src => frhodo}/bonmin/bonmin-linux64/coin-license.txt (100%) rename {src => frhodo}/bonmin/bonmin-osx/bonmin (100%) rename {src => frhodo}/bonmin/bonmin-osx/coin-license.txt (100%) rename {src => frhodo}/bonmin/bonmin-win64/bonmin.exe (100%) rename {src => frhodo}/bonmin/bonmin-win64/coin-license.txt (100%) rename {src => frhodo}/bonmin/bonmin-win64/libipoptfort.dll (100%) rename {src/calculate/optimize => frhodo/calculate}/__init__.py (100%) rename {src => frhodo}/calculate/convert_units.py (100%) rename {src => frhodo}/calculate/integrate.py (100%) rename {src => frhodo}/calculate/mech_fcns.py (99%) rename {src => frhodo}/calculate/optimize/CheKiPEUQ_from_Frhodo.py (98%) rename {src => frhodo}/calculate/optimize/CheKiPEUQ_integration_notes.txt (100%) rename {src => frhodo}/calculate/optimize/CheKiPEUQ_local/CombinationGeneratorModule.py (100%) rename {src => frhodo}/calculate/optimize/CheKiPEUQ_local/InverseProblem.py (99%) rename {src => frhodo}/calculate/optimize/CheKiPEUQ_local/LICENSE.txt (100%) rename {src => frhodo}/calculate/optimize/CheKiPEUQ_local/MANUAL.txt (100%) rename {src => frhodo}/calculate/optimize/CheKiPEUQ_local/UserInput.py (100%) rename {src => frhodo}/calculate/optimize/CheKiPEUQ_local/__init__.py (100%) rename {src => frhodo}/calculate/optimize/CheKiPEUQ_local/mumpce_custom_plotting_example.py (100%) rename {src => frhodo}/calculate/optimize/CheKiPEUQ_local/nestedObjectsFunctionsLocal.py (100%) rename {src => frhodo}/calculate/optimize/CheKiPEUQ_local/parallel_processing.py (100%) rename {src => frhodo}/calculate/optimize/CheKiPEUQ_local/plotting_functions.py (100%) rename {src => frhodo}/calculate/optimize/CheKiPEUQ_local/test_pytest.py (100%) rename {src/plot => frhodo/calculate/optimize}/__init__.py (100%) rename {src => frhodo}/calculate/optimize/fit_coeffs.py (99%) rename {src => frhodo}/calculate/optimize/fit_coeffs_pygmo.py (99%) rename {src => frhodo}/calculate/optimize/fit_fcn.py (98%) rename {src => frhodo}/calculate/optimize/mech_optimize.py (99%) rename {src => frhodo}/calculate/optimize/misc_fcns.py (99%) rename {src => frhodo}/calculate/optimize/optimize_worker.py (98%) rename {src => frhodo}/calculate/reactors.py (99%) rename {src => frhodo}/calculate/shock_fcns.py (100%) rename {src => frhodo}/calculate/smooth_data.py (100%) rename {src => frhodo}/ck2yaml.py (100%) rename {src => frhodo}/colors.py (100%) rename {src => frhodo}/config_io.py (100%) rename {src => frhodo}/data/loss_partition_fcn_tck_spline.dat (100%) rename {src => frhodo}/error_window.py (100%) rename {src => frhodo}/help_menu.py (100%) rename {src => frhodo}/ipopt/ipopt-linux64/coin-license.txt (100%) rename {src => frhodo}/ipopt/ipopt-linux64/ipopt (100%) rename {src => frhodo}/ipopt/ipopt-osx/coin-license.txt (100%) rename {src => frhodo}/ipopt/ipopt-osx/ipopt (100%) rename {src => frhodo}/ipopt/ipopt-win64/coin-license.txt (100%) rename {src => frhodo}/ipopt/ipopt-win64/ipopt.exe (100%) rename {src => frhodo}/ipopt/ipopt-win64/libipoptfort.dll (100%) rename {src => frhodo}/main.py (93%) rename {src => frhodo}/mech_widget.py (99%) rename {src => frhodo}/misc_widget.py (99%) rename {src => frhodo}/options_panel_widgets.py (99%) create mode 100644 frhodo/plot/__init__.py rename {src => frhodo}/plot/base_plot.py (99%) rename {src => frhodo}/plot/custom_mpl_ticker_formatter.py (100%) rename {src => frhodo}/plot/custom_mplscale.py (99%) rename {src => frhodo}/plot/draggable.py (100%) rename {src => frhodo}/plot/optimization_plot.py (99%) rename {src => frhodo}/plot/plot_main.py (88%) rename {src => frhodo}/plot/plot_widget.py (99%) rename {src => frhodo}/plot/raw_signal_plot.py (98%) rename {src => frhodo}/plot/signal_plot.py (99%) rename {src => frhodo}/plot/sim_explorer_plot.py (98%) rename {src => frhodo}/save_output.py (99%) rename {src => frhodo}/save_widget.py (100%) rename {src => frhodo}/series_viewer_widget.py (99%) rename {src => frhodo}/settings.py (99%) rename {src => frhodo}/sim_explorer_widget.py (98%) rename {src => frhodo}/soln2ck.py (100%) rename {src => frhodo}/thermo_widget.py (96%) create mode 100644 frhodo/version.py create mode 100644 setup.py diff --git a/src/CiteSoftLocal.py b/frhodo/CiteSoftLocal.py similarity index 100% rename from src/CiteSoftLocal.py rename to frhodo/CiteSoftLocal.py diff --git a/src/UI/error_window.ui b/frhodo/UI/error_window.ui similarity index 100% rename from src/UI/error_window.ui rename to frhodo/UI/error_window.ui diff --git a/src/UI/graphics/arrowdown.png b/frhodo/UI/graphics/arrowdown.png similarity index 100% rename from src/UI/graphics/arrowdown.png rename to frhodo/UI/graphics/arrowdown.png diff --git a/src/UI/graphics/check_icon.png b/frhodo/UI/graphics/check_icon.png similarity index 100% rename from src/UI/graphics/check_icon.png rename to frhodo/UI/graphics/check_icon.png diff --git a/src/UI/graphics/main_icon.png b/frhodo/UI/graphics/main_icon.png similarity index 100% rename from src/UI/graphics/main_icon.png rename to frhodo/UI/graphics/main_icon.png diff --git a/src/UI/graphics/x_icon.png b/frhodo/UI/graphics/x_icon.png similarity index 100% rename from src/UI/graphics/x_icon.png rename to frhodo/UI/graphics/x_icon.png diff --git a/src/UI/main_window.ui b/frhodo/UI/main_window.ui similarity index 100% rename from src/UI/main_window.ui rename to frhodo/UI/main_window.ui diff --git a/src/UI/save_dialog.ui b/frhodo/UI/save_dialog.ui similarity index 100% rename from src/UI/save_dialog.ui rename to frhodo/UI/save_dialog.ui diff --git a/src/calculate/__init__.py b/frhodo/__init__.py similarity index 100% rename from src/calculate/__init__.py rename to frhodo/__init__.py diff --git a/src/appdirs.py b/frhodo/appdirs.py similarity index 100% rename from src/appdirs.py rename to frhodo/appdirs.py diff --git a/src/bonmin/bonmin-linux64/bonmin b/frhodo/bonmin/bonmin-linux64/bonmin similarity index 100% rename from src/bonmin/bonmin-linux64/bonmin rename to frhodo/bonmin/bonmin-linux64/bonmin diff --git a/src/bonmin/bonmin-linux64/coin-license.txt b/frhodo/bonmin/bonmin-linux64/coin-license.txt similarity index 100% rename from src/bonmin/bonmin-linux64/coin-license.txt rename to frhodo/bonmin/bonmin-linux64/coin-license.txt diff --git a/src/bonmin/bonmin-osx/bonmin b/frhodo/bonmin/bonmin-osx/bonmin similarity index 100% rename from src/bonmin/bonmin-osx/bonmin rename to frhodo/bonmin/bonmin-osx/bonmin diff --git a/src/bonmin/bonmin-osx/coin-license.txt b/frhodo/bonmin/bonmin-osx/coin-license.txt similarity index 100% rename from src/bonmin/bonmin-osx/coin-license.txt rename to frhodo/bonmin/bonmin-osx/coin-license.txt diff --git a/src/bonmin/bonmin-win64/bonmin.exe b/frhodo/bonmin/bonmin-win64/bonmin.exe similarity index 100% rename from src/bonmin/bonmin-win64/bonmin.exe rename to frhodo/bonmin/bonmin-win64/bonmin.exe diff --git a/src/bonmin/bonmin-win64/coin-license.txt b/frhodo/bonmin/bonmin-win64/coin-license.txt similarity index 100% rename from src/bonmin/bonmin-win64/coin-license.txt rename to frhodo/bonmin/bonmin-win64/coin-license.txt diff --git a/src/bonmin/bonmin-win64/libipoptfort.dll b/frhodo/bonmin/bonmin-win64/libipoptfort.dll similarity index 100% rename from src/bonmin/bonmin-win64/libipoptfort.dll rename to frhodo/bonmin/bonmin-win64/libipoptfort.dll diff --git a/src/calculate/optimize/__init__.py b/frhodo/calculate/__init__.py similarity index 100% rename from src/calculate/optimize/__init__.py rename to frhodo/calculate/__init__.py diff --git a/src/calculate/convert_units.py b/frhodo/calculate/convert_units.py similarity index 100% rename from src/calculate/convert_units.py rename to frhodo/calculate/convert_units.py diff --git a/src/calculate/integrate.py b/frhodo/calculate/integrate.py similarity index 100% rename from src/calculate/integrate.py rename to frhodo/calculate/integrate.py diff --git a/src/calculate/mech_fcns.py b/frhodo/calculate/mech_fcns.py similarity index 99% rename from src/calculate/mech_fcns.py rename to frhodo/calculate/mech_fcns.py index 65e731b..707d945 100644 --- a/src/calculate/mech_fcns.py +++ b/frhodo/calculate/mech_fcns.py @@ -7,10 +7,11 @@ import cantera as ct from cantera import interrupts, cti2yaml#, ck2yaml, ctml2yaml import numpy as np -from calculate import reactors, shock_fcns, integrate -import ck2yaml from timeit import default_timer as timer +from . import reactors, shock_fcns, integrate +from .. import ck2yaml + class Chemical_Mechanism: def __init__(self): diff --git a/src/calculate/optimize/CheKiPEUQ_from_Frhodo.py b/frhodo/calculate/optimize/CheKiPEUQ_from_Frhodo.py similarity index 98% rename from src/calculate/optimize/CheKiPEUQ_from_Frhodo.py rename to frhodo/calculate/optimize/CheKiPEUQ_from_Frhodo.py index d311152..00351b6 100644 --- a/src/calculate/optimize/CheKiPEUQ_from_Frhodo.py +++ b/frhodo/calculate/optimize/CheKiPEUQ_from_Frhodo.py @@ -6,19 +6,14 @@ except: #import pathlib #sys.path.append(pathlib.Path(__file__).parent.absolute()) # add directory of **this** file to path - import calculate.optimize.CheKiPEUQ_local as CKPQ + import frhodo.calculate.optimize.CheKiPEUQ_local as CKPQ try: import CheKiPEUQ.UserInput as UserInput except: - import calculate.optimize.CheKiPEUQ_local.UserInput as UserInput + from .CheKiPEUQ_local import UserInput -try: - import CiteSoft -except: - #import pathlib #The below lines are to allow CiteSoftLocal to be called regardless of user's working directory. - #sys.path.append(pathlib.Path(__file__).parent.absolute()) # add directory of **this** file to path - import CiteSoftLocal as CiteSoft +import frhodo.CiteSoftLocal as CiteSoft #get things ready for CiteSoft entry... software_name = "CheKiPEUQ Bayesian Parameter Estimation" diff --git a/src/calculate/optimize/CheKiPEUQ_integration_notes.txt b/frhodo/calculate/optimize/CheKiPEUQ_integration_notes.txt similarity index 100% rename from src/calculate/optimize/CheKiPEUQ_integration_notes.txt rename to frhodo/calculate/optimize/CheKiPEUQ_integration_notes.txt diff --git a/src/calculate/optimize/CheKiPEUQ_local/CombinationGeneratorModule.py b/frhodo/calculate/optimize/CheKiPEUQ_local/CombinationGeneratorModule.py similarity index 100% rename from src/calculate/optimize/CheKiPEUQ_local/CombinationGeneratorModule.py rename to frhodo/calculate/optimize/CheKiPEUQ_local/CombinationGeneratorModule.py diff --git a/src/calculate/optimize/CheKiPEUQ_local/InverseProblem.py b/frhodo/calculate/optimize/CheKiPEUQ_local/InverseProblem.py similarity index 99% rename from src/calculate/optimize/CheKiPEUQ_local/InverseProblem.py rename to frhodo/calculate/optimize/CheKiPEUQ_local/InverseProblem.py index 38799cc..2251cff 100644 --- a/src/calculate/optimize/CheKiPEUQ_local/InverseProblem.py +++ b/frhodo/calculate/optimize/CheKiPEUQ_local/InverseProblem.py @@ -9,16 +9,13 @@ from collections.abc import Iterable #import mumce_py.Project as mumce_pyProject #FIXME: Eric to fix plotting/graphing issue described in issue 9 -- https://github.com/AdityaSavara/ODE-KIN-BAYES-SG-EW/issues/9 #import mumce_py.solution mumce_pySolution -try: - import CiteSoft -except: - #import pathlib #The below lines are to allow CiteSoftLocal to be called regardless of user's working directory. - #sys.path.append(pathlib.Path(__file__).parent.absolute()) # add directory of **this** file to path - import CiteSoftLocal as CiteSoft + +import frhodo.CiteSoftLocal as CiteSoft + try: import UnitTesterSG.nestedObjectsFunctions as nestedObjectsFunctions -except: - import calculate.optimize.CheKiPEUQ_local.nestedObjectsFunctionsLocal as nestedObjectsFunctions +except ImportError: + import frhodo.calculate.optimize.CheKiPEUQ_local.nestedObjectsFunctionsLocal as nestedObjectsFunctions class parameter_estimation: #Inside this class, a UserInput namespace is provided. This has dictionaries of UserInput choices. diff --git a/src/calculate/optimize/CheKiPEUQ_local/LICENSE.txt b/frhodo/calculate/optimize/CheKiPEUQ_local/LICENSE.txt similarity index 100% rename from src/calculate/optimize/CheKiPEUQ_local/LICENSE.txt rename to frhodo/calculate/optimize/CheKiPEUQ_local/LICENSE.txt diff --git a/src/calculate/optimize/CheKiPEUQ_local/MANUAL.txt b/frhodo/calculate/optimize/CheKiPEUQ_local/MANUAL.txt similarity index 100% rename from src/calculate/optimize/CheKiPEUQ_local/MANUAL.txt rename to frhodo/calculate/optimize/CheKiPEUQ_local/MANUAL.txt diff --git a/src/calculate/optimize/CheKiPEUQ_local/UserInput.py b/frhodo/calculate/optimize/CheKiPEUQ_local/UserInput.py similarity index 100% rename from src/calculate/optimize/CheKiPEUQ_local/UserInput.py rename to frhodo/calculate/optimize/CheKiPEUQ_local/UserInput.py diff --git a/src/calculate/optimize/CheKiPEUQ_local/__init__.py b/frhodo/calculate/optimize/CheKiPEUQ_local/__init__.py similarity index 100% rename from src/calculate/optimize/CheKiPEUQ_local/__init__.py rename to frhodo/calculate/optimize/CheKiPEUQ_local/__init__.py diff --git a/src/calculate/optimize/CheKiPEUQ_local/mumpce_custom_plotting_example.py b/frhodo/calculate/optimize/CheKiPEUQ_local/mumpce_custom_plotting_example.py similarity index 100% rename from src/calculate/optimize/CheKiPEUQ_local/mumpce_custom_plotting_example.py rename to frhodo/calculate/optimize/CheKiPEUQ_local/mumpce_custom_plotting_example.py diff --git a/src/calculate/optimize/CheKiPEUQ_local/nestedObjectsFunctionsLocal.py b/frhodo/calculate/optimize/CheKiPEUQ_local/nestedObjectsFunctionsLocal.py similarity index 100% rename from src/calculate/optimize/CheKiPEUQ_local/nestedObjectsFunctionsLocal.py rename to frhodo/calculate/optimize/CheKiPEUQ_local/nestedObjectsFunctionsLocal.py diff --git a/src/calculate/optimize/CheKiPEUQ_local/parallel_processing.py b/frhodo/calculate/optimize/CheKiPEUQ_local/parallel_processing.py similarity index 100% rename from src/calculate/optimize/CheKiPEUQ_local/parallel_processing.py rename to frhodo/calculate/optimize/CheKiPEUQ_local/parallel_processing.py diff --git a/src/calculate/optimize/CheKiPEUQ_local/plotting_functions.py b/frhodo/calculate/optimize/CheKiPEUQ_local/plotting_functions.py similarity index 100% rename from src/calculate/optimize/CheKiPEUQ_local/plotting_functions.py rename to frhodo/calculate/optimize/CheKiPEUQ_local/plotting_functions.py diff --git a/src/calculate/optimize/CheKiPEUQ_local/test_pytest.py b/frhodo/calculate/optimize/CheKiPEUQ_local/test_pytest.py similarity index 100% rename from src/calculate/optimize/CheKiPEUQ_local/test_pytest.py rename to frhodo/calculate/optimize/CheKiPEUQ_local/test_pytest.py diff --git a/src/plot/__init__.py b/frhodo/calculate/optimize/__init__.py similarity index 100% rename from src/plot/__init__.py rename to frhodo/calculate/optimize/__init__.py diff --git a/src/calculate/optimize/fit_coeffs.py b/frhodo/calculate/optimize/fit_coeffs.py similarity index 99% rename from src/calculate/optimize/fit_coeffs.py rename to frhodo/calculate/optimize/fit_coeffs.py index 4b5af1c..e53827c 100644 --- a/src/calculate/optimize/fit_coeffs.py +++ b/frhodo/calculate/optimize/fit_coeffs.py @@ -12,8 +12,9 @@ from timeit import default_timer as timer import itertools -from calculate.convert_units import OoM, Bisymlog -from calculate.optimize.misc_fcns import penalized_loss_fcn, set_arrhenius_bnds, min_pos_system_value, max_pos_system_value +from ..convert_units import OoM, Bisymlog +from ..optimize.misc_fcns import penalized_loss_fcn, set_arrhenius_bnds, min_pos_system_value, max_pos_system_value +import frhodo Ru = ct.gas_constant # Ru = 1.98720425864083 @@ -134,7 +135,7 @@ def ln_arrhenius_jac(T, *args): return popt -path = {'main': pathlib.Path(sys.argv[0]).parents[0].resolve()} +path = {'main': pathlib.Path(frhodo.__file__).parents[0].resolve()} # TODO (wardlt) - Make a globally-accessible Path OS_type = platform.system() if OS_type == 'Windows': path['bonmin'] = path['main'] / 'bonmin/bonmin-win64/bonmin.exe' diff --git a/src/calculate/optimize/fit_coeffs_pygmo.py b/frhodo/calculate/optimize/fit_coeffs_pygmo.py similarity index 99% rename from src/calculate/optimize/fit_coeffs_pygmo.py rename to frhodo/calculate/optimize/fit_coeffs_pygmo.py index 1fc8716..fbaf8d6 100644 --- a/src/calculate/optimize/fit_coeffs_pygmo.py +++ b/frhodo/calculate/optimize/fit_coeffs_pygmo.py @@ -12,8 +12,9 @@ from timeit import default_timer as timer import itertools -from calculate.convert_units import OoM -from calculate.optimize.misc_fcns import penalized_loss_fcn, set_arrhenius_bnds +import frhodo +from ..convert_units import OoM +from ..optimize.misc_fcns import penalized_loss_fcn, set_arrhenius_bnds Ru = ct.gas_constant # Ru = 1.98720425864083 @@ -50,7 +51,7 @@ def ln_arrhenius_k(T, Ea, ln_A, n): # LPL, HPL, Fcent def fit_arrhenius(rates, T, x0=[], coefNames=default_arrhenius_coefNames, bnds=[], loss='linear'): def fit_fcn_decorator(x0, alter_idx, jac=False): def set_coeffs(*args): - coeffs = x0 + coeffs = x0git for n, idx in enumerate(alter_idx): coeffs[idx] = args[n] return coeffs @@ -118,7 +119,7 @@ def ln_arrhenius_jac(T, *args): return popt -path = {'main': pathlib.Path(sys.argv[0]).parents[0].resolve()} +path = {'main': pathlib.Path(frhodo.__file__).parents[0].resolve()} # TODO (wardlt) - Make a globally-accessible Path OS_type = platform.system() if OS_type == 'Windows': path['bonmin'] = path['main'] / 'bonmin/bonmin-win64/bonmin.exe' diff --git a/src/calculate/optimize/fit_fcn.py b/frhodo/calculate/optimize/fit_fcn.py similarity index 98% rename from src/calculate/optimize/fit_fcn.py rename to frhodo/calculate/optimize/fit_fcn.py index e32db8d..dd40365 100644 --- a/src/calculate/optimize/fit_fcn.py +++ b/frhodo/calculate/optimize/fit_fcn.py @@ -10,11 +10,11 @@ from scipy import stats from copy import deepcopy -from calculate.mech_fcns import Chemical_Mechanism -from calculate.convert_units import OoM, Bisymlog -from calculate.optimize.misc_fcns import weighted_quantile, outlier, penalized_loss_fcn -from calculate.optimize.fit_coeffs import fit_coeffs -from calculate.optimize.CheKiPEUQ_from_Frhodo import CheKiPEUQ_Frhodo_interface +from ..mech_fcns import Chemical_Mechanism +from ..convert_units import OoM, Bisymlog +from ..optimize.misc_fcns import weighted_quantile, outlier, penalized_loss_fcn +from ..optimize.fit_coeffs import fit_coeffs +from ..optimize.CheKiPEUQ_from_Frhodo import CheKiPEUQ_Frhodo_interface mpMech = {} def initialize_parallel_worker(mech_dict, species_dict, coeffs, coeffs_bnds, rate_bnds): diff --git a/src/calculate/optimize/mech_optimize.py b/frhodo/calculate/optimize/mech_optimize.py similarity index 99% rename from src/calculate/optimize/mech_optimize.py rename to frhodo/calculate/optimize/mech_optimize.py index f15a6d0..7cca45c 100644 --- a/src/calculate/optimize/mech_optimize.py +++ b/frhodo/calculate/optimize/mech_optimize.py @@ -12,10 +12,10 @@ from scipy import stats -from calculate.optimize.optimize_worker import Worker -from calculate.optimize.fit_fcn import update_mech_coef_opt -from calculate.optimize.misc_fcns import rates, set_bnds -from calculate.optimize.fit_coeffs import fit_generic as Troe_fit +from ..optimize.optimize_worker import Worker +from ..optimize.fit_fcn import update_mech_coef_opt +from ..optimize.misc_fcns import rates, set_bnds +from ..optimize.fit_coeffs import fit_generic as Troe_fit Ru = ct.gas_constant default_arrhenius_coefNames = ['activation_energy', 'pre_exponential_factor', 'temperature_exponent'] diff --git a/src/calculate/optimize/misc_fcns.py b/frhodo/calculate/optimize/misc_fcns.py similarity index 99% rename from src/calculate/optimize/misc_fcns.py rename to frhodo/calculate/optimize/misc_fcns.py index 8e38a9d..6653c27 100644 --- a/src/calculate/optimize/misc_fcns.py +++ b/frhodo/calculate/optimize/misc_fcns.py @@ -8,6 +8,8 @@ import cantera as ct import pathlib, sys +import frhodo + Ru = ct.gas_constant min_pos_system_value = (np.finfo(float).tiny*(1E20))**(1/2) @@ -20,7 +22,7 @@ # interpolation function for Z from loss function -path = {'main': pathlib.Path(sys.argv[0]).parents[0].resolve()} +path = {'main': pathlib.Path(frhodo.__file__).parents[0].resolve()} path['Z_tck_spline.dat'] = path['main'] / 'data/loss_partition_fcn_tck_spline.dat' tck = [] diff --git a/src/calculate/optimize/optimize_worker.py b/frhodo/calculate/optimize/optimize_worker.py similarity index 98% rename from src/calculate/optimize/optimize_worker.py rename to frhodo/calculate/optimize/optimize_worker.py index 5fb4019..4f2d900 100644 --- a/src/calculate/optimize/optimize_worker.py +++ b/frhodo/calculate/optimize/optimize_worker.py @@ -12,8 +12,9 @@ from timeit import default_timer as timer -from calculate.optimize.fit_fcn import initialize_parallel_worker, Fit_Fun -from calculate.optimize.misc_fcns import rates +from ..optimize.fit_fcn import initialize_parallel_worker, Fit_Fun +from ..optimize.misc_fcns import rates +import frhodo debug = True @@ -210,7 +211,7 @@ class WorkerSignals(QObject): 'Optimization failed: Out of memory', 'Optimization failed: Roundoff errors limited progress', 'Optimization failed: Forced termination'] -path = {'main': pathlib.Path(sys.argv[0]).parents[0].resolve()} +path = {'main': pathlib.Path(frhodo.__file__).parents[0].resolve()} # TODO (wardlt) - Make a globally-accessible Path OS_type = platform.system() if OS_type == 'Windows': path['bonmin'] = path['main'] / 'bonmin/bonmin-win64/bonmin.exe' diff --git a/src/calculate/reactors.py b/frhodo/calculate/reactors.py similarity index 99% rename from src/calculate/reactors.py rename to frhodo/calculate/reactors.py index 2c9bdbe..bd4480c 100644 --- a/src/calculate/reactors.py +++ b/frhodo/calculate/reactors.py @@ -7,8 +7,8 @@ import cantera as ct from cantera import interrupts, cti2yaml#, ck2yaml, ctml2yaml import numpy as np -from calculate import shock_fcns, integrate -import ck2yaml +from . import shock_fcns, integrate +import frhodo.ck2yaml from timeit import default_timer as timer diff --git a/src/calculate/shock_fcns.py b/frhodo/calculate/shock_fcns.py similarity index 100% rename from src/calculate/shock_fcns.py rename to frhodo/calculate/shock_fcns.py diff --git a/src/calculate/smooth_data.py b/frhodo/calculate/smooth_data.py similarity index 100% rename from src/calculate/smooth_data.py rename to frhodo/calculate/smooth_data.py diff --git a/src/ck2yaml.py b/frhodo/ck2yaml.py similarity index 100% rename from src/ck2yaml.py rename to frhodo/ck2yaml.py diff --git a/src/colors.py b/frhodo/colors.py similarity index 100% rename from src/colors.py rename to frhodo/colors.py diff --git a/src/config_io.py b/frhodo/config_io.py similarity index 100% rename from src/config_io.py rename to frhodo/config_io.py diff --git a/src/data/loss_partition_fcn_tck_spline.dat b/frhodo/data/loss_partition_fcn_tck_spline.dat similarity index 100% rename from src/data/loss_partition_fcn_tck_spline.dat rename to frhodo/data/loss_partition_fcn_tck_spline.dat diff --git a/src/error_window.py b/frhodo/error_window.py similarity index 100% rename from src/error_window.py rename to frhodo/error_window.py diff --git a/src/help_menu.py b/frhodo/help_menu.py similarity index 100% rename from src/help_menu.py rename to frhodo/help_menu.py diff --git a/src/ipopt/ipopt-linux64/coin-license.txt b/frhodo/ipopt/ipopt-linux64/coin-license.txt similarity index 100% rename from src/ipopt/ipopt-linux64/coin-license.txt rename to frhodo/ipopt/ipopt-linux64/coin-license.txt diff --git a/src/ipopt/ipopt-linux64/ipopt b/frhodo/ipopt/ipopt-linux64/ipopt similarity index 100% rename from src/ipopt/ipopt-linux64/ipopt rename to frhodo/ipopt/ipopt-linux64/ipopt diff --git a/src/ipopt/ipopt-osx/coin-license.txt b/frhodo/ipopt/ipopt-osx/coin-license.txt similarity index 100% rename from src/ipopt/ipopt-osx/coin-license.txt rename to frhodo/ipopt/ipopt-osx/coin-license.txt diff --git a/src/ipopt/ipopt-osx/ipopt b/frhodo/ipopt/ipopt-osx/ipopt similarity index 100% rename from src/ipopt/ipopt-osx/ipopt rename to frhodo/ipopt/ipopt-osx/ipopt diff --git a/src/ipopt/ipopt-win64/coin-license.txt b/frhodo/ipopt/ipopt-win64/coin-license.txt similarity index 100% rename from src/ipopt/ipopt-win64/coin-license.txt rename to frhodo/ipopt/ipopt-win64/coin-license.txt diff --git a/src/ipopt/ipopt-win64/ipopt.exe b/frhodo/ipopt/ipopt-win64/ipopt.exe similarity index 100% rename from src/ipopt/ipopt-win64/ipopt.exe rename to frhodo/ipopt/ipopt-win64/ipopt.exe diff --git a/src/ipopt/ipopt-win64/libipoptfort.dll b/frhodo/ipopt/ipopt-win64/libipoptfort.dll similarity index 100% rename from src/ipopt/ipopt-win64/libipoptfort.dll rename to frhodo/ipopt/ipopt-win64/libipoptfort.dll diff --git a/src/main.py b/frhodo/main.py similarity index 93% rename from src/main.py rename to frhodo/main.py index e3712c9..29c8f57 100644 --- a/src/main.py +++ b/frhodo/main.py @@ -5,21 +5,20 @@ # and licensed under BSD-3-Clause. See License.txt in the top-level # directory for license and copyright information. -version = '1.3.1' import os, sys, platform, multiprocessing, pathlib, ctypes # os.environ['QT_API'] = 'pyside2' # forces pyside2 +from typing import Tuple from qtpy.QtWidgets import QMainWindow, QApplication, QMessageBox from qtpy import uic, QtCore, QtGui - import numpy as np -# from timeit import default_timer as timer -from plot.plot_main import All_Plots as plot -from misc_widget import MessageWindow -from calculate import mech_fcns, reactors, convert_units -import appdirs, options_panel_widgets, sim_explorer_widget -import settings, config_io, save_widget, error_window, help_menu +from .plot.plot_main import All_Plots as plot +from .misc_widget import MessageWindow +from .calculate import mech_fcns, reactors, convert_units +from .version import __version__ +from . import appdirs, options_panel_widgets, sim_explorer_widget +from . import settings, config_io, save_widget, error_window, help_menu if os.environ['QT_API'] == 'pyside2': # Silence warning: "Qt WebEngine seems to be initialized from a plugin." QApplication.setAttribute(QtCore.Qt.AA_ShareOpenGLContexts) @@ -31,7 +30,7 @@ QApplication.setAttribute(QtCore.Qt.AA_UseHighDpiPixmaps, True) # set main folder -path = {'main': pathlib.Path(sys.argv[0]).parents[0].resolve()} +path = {'main': pathlib.Path(__file__).parents[0].resolve()} # set appdata folder using AppDirs library (but just using the source code file) dirs = appdirs.AppDirs(appname='Frhodo', roaming=True, appauthor=False) @@ -86,7 +85,7 @@ def __init__(self, app, path): self.initialize_settings() # ~ 4 sec # Setup help menu - self.version = version + self.version = __version__ help_menu.HelpMenu(self) def initialize_settings(self): # TODO: Solving for loaded shock twice @@ -290,17 +289,31 @@ def run_single(self, event=None, t_save=None, rxn_changed=False): # def raise_error(self): # assert False - -if __name__ == '__main__': + +def launch_gui() -> Tuple[QApplication, Main]: + """Launch the GUI + + Returns: + - The QApplication instance + - Link to the main window + """ if platform.system() == 'Windows': # this is required for pyinstaller on windows multiprocessing.freeze_support() - if getattr(sys, 'frozen', False): # if frozen minimize console immediately + if getattr(sys, 'frozen', False): # if frozen minimize console immediately ctypes.windll.user32.ShowWindow(ctypes.windll.kernel32.GetConsoleWindow(), 0) - + app = QApplication(sys.argv) sys.excepthook = error_window.excepthookDecorator(app, path, shut_down) - main = Main(app, path) + main_window = Main(app, path) + return app, main_window + + +def main(): + """Launch the GUI and then block until it finishes""" + # Launch the application + app, main_window = launch_gui() + + # Pass the exit code forward sys.exit(app.exec_()) - diff --git a/src/mech_widget.py b/frhodo/mech_widget.py similarity index 99% rename from src/mech_widget.py rename to frhodo/mech_widget.py index 1963f32..052858c 100644 --- a/src/mech_widget.py +++ b/frhodo/mech_widget.py @@ -2,17 +2,17 @@ # and licensed under BSD-3-Clause. See License.txt in the top-level # directory for license and copyright information. +from copy import deepcopy import sys, ast, re -import misc_widget + import cantera as ct import numpy as np -from copy import deepcopy -from functools import partial from scipy.optimize import root_scalar from qtpy.QtWidgets import * from qtpy import QtWidgets, QtGui, QtCore -from timeit import default_timer as timer +from . import misc_widget + def silentSetValue(obj, value): obj.blockSignals(True) # stop changing text from signaling diff --git a/src/misc_widget.py b/frhodo/misc_widget.py similarity index 99% rename from src/misc_widget.py rename to frhodo/misc_widget.py index 36f3e59..ccaebcf 100644 --- a/src/misc_widget.py +++ b/frhodo/misc_widget.py @@ -6,7 +6,7 @@ import numpy as np from qtpy.QtWidgets import * from qtpy import QtWidgets, QtGui, QtCore -from calculate.convert_units import OoM +from .calculate.convert_units import OoM # Regular expression to find floats. Match groups are the whole string, the # whole coefficient, the decimal part of the coefficient, and the exponent diff --git a/src/options_panel_widgets.py b/frhodo/options_panel_widgets.py similarity index 99% rename from src/options_panel_widgets.py rename to frhodo/options_panel_widgets.py index a88b59e..f5bb4a4 100644 --- a/src/options_panel_widgets.py +++ b/frhodo/options_panel_widgets.py @@ -3,17 +3,20 @@ # directory for license and copyright information. import pathlib, os, sys +from copy import deepcopy + import numpy as np from scipy.optimize import minimize import nlopt -import mech_widget, misc_widget, thermo_widget, series_viewer_widget, save_output -from calculate import shock_fcns -from calculate.optimize.mech_optimize import Multithread_Optimize -from calculate.convert_units import OoM -from settings import double_sigmoid from qtpy.QtWidgets import * from qtpy import QtWidgets, QtGui, QtCore -from copy import deepcopy + +from . import mech_widget, misc_widget, thermo_widget, series_viewer_widget, save_output +from .calculate import shock_fcns +from .calculate.optimize.mech_optimize import Multithread_Optimize +from .calculate.convert_units import OoM +from .settings import double_sigmoid + class Initialize(QtCore.QObject): diff --git a/frhodo/plot/__init__.py b/frhodo/plot/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/plot/base_plot.py b/frhodo/plot/base_plot.py similarity index 99% rename from src/plot/base_plot.py rename to frhodo/plot/base_plot.py index 57f01cd..946fab9 100644 --- a/src/plot/base_plot.py +++ b/frhodo/plot/base_plot.py @@ -24,9 +24,9 @@ #mpl.use("module://mplcairo.qt") # This implements mplcairo, faster/more accurate. Issues with other OSes? import numpy as np -from calculate.convert_units import Bisymlog -from plot.custom_mplscale import * -from plot.custom_mpl_ticker_formatter import * +from ..calculate.convert_units import Bisymlog +from .custom_mplscale import * +from .custom_mpl_ticker_formatter import * from timeit import default_timer as timer diff --git a/src/plot/custom_mpl_ticker_formatter.py b/frhodo/plot/custom_mpl_ticker_formatter.py similarity index 100% rename from src/plot/custom_mpl_ticker_formatter.py rename to frhodo/plot/custom_mpl_ticker_formatter.py diff --git a/src/plot/custom_mplscale.py b/frhodo/plot/custom_mplscale.py similarity index 99% rename from src/plot/custom_mplscale.py rename to frhodo/plot/custom_mplscale.py index 935aaeb..d1c949d 100644 --- a/src/plot/custom_mplscale.py +++ b/frhodo/plot/custom_mplscale.py @@ -6,7 +6,7 @@ from matplotlib import scale as mplscale import numpy as np -from calculate.convert_units import Bisymlog +from frhodo.calculate.convert_units import Bisymlog class AbsoluteLogScale(mplscale.LogScale): name = 'abslog' diff --git a/src/plot/draggable.py b/frhodo/plot/draggable.py similarity index 100% rename from src/plot/draggable.py rename to frhodo/plot/draggable.py diff --git a/src/plot/optimization_plot.py b/frhodo/plot/optimization_plot.py similarity index 99% rename from src/plot/optimization_plot.py rename to frhodo/plot/optimization_plot.py index 1bc5a36..d30a072 100644 --- a/src/plot/optimization_plot.py +++ b/frhodo/plot/optimization_plot.py @@ -2,12 +2,12 @@ # and licensed under BSD-3-Clause. See License.txt in the top-level # directory for license and copyright information. -import colors +from frhodo import colors import matplotlib as mpl import numpy as np -from plot.base_plot import Base_Plot +from .base_plot import Base_Plot colormap = colors.colormap(reorder_from=1, num_shift=4) diff --git a/src/plot/plot_main.py b/frhodo/plot/plot_main.py similarity index 88% rename from src/plot/plot_main.py rename to frhodo/plot/plot_main.py index b0cd415..771b59b 100644 --- a/src/plot/plot_main.py +++ b/frhodo/plot/plot_main.py @@ -2,7 +2,7 @@ # and licensed under BSD-3-Clause. See License.txt in the top-level # directory for license and copyright information. -from plot import raw_signal_plot, signal_plot, sim_explorer_plot, optimization_plot, plot_widget +from . import raw_signal_plot, signal_plot, sim_explorer_plot, optimization_plot, plot_widget class All_Plots: # container to hold all plots diff --git a/src/plot/plot_widget.py b/frhodo/plot/plot_widget.py similarity index 99% rename from src/plot/plot_widget.py rename to frhodo/plot/plot_widget.py index b668e57..034d928 100644 --- a/src/plot/plot_widget.py +++ b/frhodo/plot/plot_widget.py @@ -7,7 +7,7 @@ from qtpy import QtWidgets, QtGui, QtCore from copy import deepcopy -import misc_widget +from .. import misc_widget class Observable_Widgets(QtCore.QObject): def __init__(self, parent): diff --git a/src/plot/raw_signal_plot.py b/frhodo/plot/raw_signal_plot.py similarity index 98% rename from src/plot/raw_signal_plot.py rename to frhodo/plot/raw_signal_plot.py index dd7b7e1..1139caa 100644 --- a/src/plot/raw_signal_plot.py +++ b/frhodo/plot/raw_signal_plot.py @@ -9,8 +9,8 @@ import numpy as np from scipy import stats -from plot.base_plot import Base_Plot -from plot.draggable import Draggable +from .base_plot import Base_Plot +from .draggable import Draggable class Plot(Base_Plot): diff --git a/src/plot/signal_plot.py b/frhodo/plot/signal_plot.py similarity index 99% rename from src/plot/signal_plot.py rename to frhodo/plot/signal_plot.py index f543c12..31d1468 100644 --- a/src/plot/signal_plot.py +++ b/frhodo/plot/signal_plot.py @@ -8,9 +8,9 @@ import numpy as np from scipy import stats -from calculate.convert_units import OoM -from plot.base_plot import Base_Plot -from plot.draggable import Draggable +from frhodo.calculate.convert_units import OoM +from .base_plot import Base_Plot +from .draggable import Draggable def shape_data(x, y): diff --git a/src/plot/sim_explorer_plot.py b/frhodo/plot/sim_explorer_plot.py similarity index 98% rename from src/plot/sim_explorer_plot.py rename to frhodo/plot/sim_explorer_plot.py index cbb41c4..a420b87 100644 --- a/src/plot/sim_explorer_plot.py +++ b/frhodo/plot/sim_explorer_plot.py @@ -5,8 +5,8 @@ import matplotlib as mpl import numpy as np -from plot.base_plot import Base_Plot -from plot.draggable import DraggableLegend +from .base_plot import Base_Plot +from .draggable import DraggableLegend class Plot(Base_Plot): diff --git a/src/save_output.py b/frhodo/save_output.py similarity index 99% rename from src/save_output.py rename to frhodo/save_output.py index 9f30668..2162e63 100644 --- a/src/save_output.py +++ b/frhodo/save_output.py @@ -4,10 +4,10 @@ import numpy as np from tabulate import tabulate -import soln2ck -import pathlib # from cantera import ck2cti # Maybe later use ck2cti.Parser.writeCTI to write cti file +from . import soln2ck + class Save: def __init__(self, parent): self.parent = parent diff --git a/src/save_widget.py b/frhodo/save_widget.py similarity index 100% rename from src/save_widget.py rename to frhodo/save_widget.py diff --git a/src/series_viewer_widget.py b/frhodo/series_viewer_widget.py similarity index 99% rename from src/series_viewer_widget.py rename to frhodo/series_viewer_widget.py index d7e6c77..a3dc9f8 100644 --- a/src/series_viewer_widget.py +++ b/frhodo/series_viewer_widget.py @@ -2,7 +2,7 @@ # and licensed under BSD-3-Clause. See License.txt in the top-level # directory for license and copyright information. -from calculate import shock_fcns +from .calculate import shock_fcns import numpy as np from qtpy.QtWidgets import * from qtpy import QtWidgets, QtGui, QtCore diff --git a/src/settings.py b/frhodo/settings.py similarity index 99% rename from src/settings.py rename to frhodo/settings.py index 471727c..40738bd 100644 --- a/src/settings.py +++ b/frhodo/settings.py @@ -10,8 +10,8 @@ from scipy.interpolate import CubicSpline from qtpy import QtCore -from calculate.convert_units import OoM, RoundToSigFigs -from calculate.smooth_data import dual_tree_complex_wavelet_filter +from .calculate.convert_units import OoM, RoundToSigFigs +from .calculate.smooth_data import dual_tree_complex_wavelet_filter min_pos_system_value = (np.finfo(float).tiny*(1E20))**(1/2) max_pos_system_value = (np.finfo(float).max*(1E-20))**(1/2) diff --git a/src/sim_explorer_widget.py b/frhodo/sim_explorer_widget.py similarity index 98% rename from src/sim_explorer_widget.py rename to frhodo/sim_explorer_widget.py index 21648dc..5754be5 100644 --- a/src/sim_explorer_widget.py +++ b/frhodo/sim_explorer_widget.py @@ -3,13 +3,12 @@ # directory for license and copyright information. import numpy as np -from calculate import shock_fcns -import mech_widget, misc_widget, thermo_widget, series_viewer_widget, save_output from qtpy.QtWidgets import * from qtpy import QtWidgets, QtGui, QtCore from copy import deepcopy -import re -from timeit import default_timer as timer + +from .calculate import shock_fcns +from . import mech_widget, misc_widget, thermo_widget, series_viewer_widget, save_output class SIM_Explorer_Widgets(QtCore.QObject): def __init__(self, parent): diff --git a/src/soln2ck.py b/frhodo/soln2ck.py similarity index 100% rename from src/soln2ck.py rename to frhodo/soln2ck.py diff --git a/src/thermo_widget.py b/frhodo/thermo_widget.py similarity index 96% rename from src/thermo_widget.py rename to frhodo/thermo_widget.py index 023607b..c92159d 100644 --- a/src/thermo_widget.py +++ b/frhodo/thermo_widget.py @@ -2,15 +2,11 @@ # and licensed under BSD-3-Clause. See License.txt in the top-level # directory for license and copyright information. -import sys, ast, re -import misc_widget -import cantera as ct -import numpy as np -from copy import deepcopy -from scipy.optimize import root_scalar from qtpy.QtWidgets import * from qtpy import QtWidgets, QtGui, QtCore +from . import misc_widget + def silentSetValue(obj, value): obj.blockSignals(True) # stop changing text from signaling obj.setValue(value) diff --git a/frhodo/version.py b/frhodo/version.py new file mode 100644 index 0000000..3b949b4 --- /dev/null +++ b/frhodo/version.py @@ -0,0 +1,4 @@ +"""Single source of truth for package version""" +# https://packaging.python.org/en/latest/guides/single-sourcing-package-version/ + +__version__ = '1.3.1' diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..15d522c --- /dev/null +++ b/setup.py @@ -0,0 +1,20 @@ +from setuptools import setup, find_packages +from pathlib import Path + +# single source of truth for package version +version_ns = {} +with (Path("frhodo") / "version.py").open() as f: + exec(f.read(), version_ns) +version = version_ns['__version__'] + +setup( + name='frhodo', + version='0.0.1', + packages=find_packages(), + description='Simulate experimental data and optimize chemical kinetics mechanisms with this GUI-based application', + entry_points={ + 'console_scripts': [ + 'frhodo=frhodo.main:main' + ] + } +) From 682b28fad3536530afcc1cf0fc13ae0fbaa7536d Mon Sep 17 00:00:00 2001 From: Logan Ward Date: Fri, 19 Aug 2022 13:10:46 -0400 Subject: [PATCH 03/29] Install Frhodo with the environment --- frhodo_environment.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/frhodo_environment.yml b/frhodo_environment.yml index daecfb0..73381d0 100644 --- a/frhodo_environment.yml +++ b/frhodo_environment.yml @@ -22,3 +22,4 @@ dependencies: - pip: - dtcwt==0.12.0 - rbfopt==4.2.2 + - -e . From 2d9de9349280c3ddc3ed7768b07237e62d95340b Mon Sep 17 00:00:00 2001 From: Logan Ward Date: Mon, 22 Aug 2022 11:51:59 -0400 Subject: [PATCH 04/29] Add testing to env, version # to setup --- frhodo_environment.yml | 1 + setup.py | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/frhodo_environment.yml b/frhodo_environment.yml index 73381d0..06c53e1 100644 --- a/frhodo_environment.yml +++ b/frhodo_environment.yml @@ -19,6 +19,7 @@ dependencies: - requests=2.27.1 - scipy=1.7.3 - tabulate=0.8.9 + - pytest - pip: - dtcwt==0.12.0 - rbfopt==4.2.2 diff --git a/setup.py b/setup.py index 15d522c..fd3d99d 100644 --- a/setup.py +++ b/setup.py @@ -9,7 +9,7 @@ setup( name='frhodo', - version='0.0.1', + version=version, packages=find_packages(), description='Simulate experimental data and optimize chemical kinetics mechanisms with this GUI-based application', entry_points={ From cc560d4fbd93bd263e15f4a56968fedd9856d10d Mon Sep 17 00:00:00 2001 From: Logan Ward Date: Mon, 22 Aug 2022 11:54:11 -0400 Subject: [PATCH 05/29] Mostly adding documentation --- frhodo/bonmin/bonmin-win64/bonmin.exe | Bin frhodo/calculate/mech_fcns.py | 46 ++++++- frhodo/calculate/optimize/fit_fcn.py | 54 +++++++- frhodo/calculate/optimize/mech_optimize.py | 24 +++- frhodo/calculate/optimize/optimize_worker.py | 48 +++++-- frhodo/calculate/reactors.py | 2 + frhodo/config_io.py | 7 +- frhodo/ipopt/ipopt-win64/ipopt.exe | Bin frhodo/main.py | 127 +++++++++++-------- frhodo/mech_widget.py | 30 ++++- frhodo/options_panel_widgets.py | 18 ++- frhodo/settings.py | 24 +++- frhodo/version.py | 8 +- 13 files changed, 290 insertions(+), 98 deletions(-) mode change 100644 => 100755 frhodo/bonmin/bonmin-win64/bonmin.exe mode change 100644 => 100755 frhodo/ipopt/ipopt-win64/ipopt.exe diff --git a/frhodo/bonmin/bonmin-win64/bonmin.exe b/frhodo/bonmin/bonmin-win64/bonmin.exe old mode 100644 new mode 100755 diff --git a/frhodo/calculate/mech_fcns.py b/frhodo/calculate/mech_fcns.py index 707d945..d767fa4 100644 --- a/frhodo/calculate/mech_fcns.py +++ b/frhodo/calculate/mech_fcns.py @@ -4,6 +4,8 @@ import os, io, stat, contextlib, pathlib, time from copy import deepcopy +from typing import Tuple + import cantera as ct from cantera import interrupts, cti2yaml#, ck2yaml, ctml2yaml import numpy as np @@ -295,6 +297,7 @@ def copy_bnds(new_bnds, bnds, rxnIdx, bnds_type, keys=[]): raise(f'{rxn} is a {rxn.__class__.__name__} and is currently unsupported in Frhodo, but this error should never be seen') def get_coeffs_keys(self, rxn, coefAbbr, rxnIdx=None): + """Get the keys in a reaction dictionary associated that define the par""" if type(rxn) in [ct.ElementaryReaction, ct.ThreeBodyReaction]: bnds_key = 'rate' coef_key = 0 @@ -320,6 +323,9 @@ def get_coeffs_keys(self, rxn, coefAbbr, rxnIdx=None): else: coef_key = bnds_key = 'falloff_parameters' + else: + raise ValueError(f'rxn type {rxn} not yet supported') + return coef_key, bnds_key def set_thermo_expression_coeffs(self): # TODO Doesn't work with NASA 9 @@ -525,18 +531,48 @@ def get_M(rxn): return M - def run(self, reactor_choice, t_end, T_reac, P_reac, mix, **kwargs): + def run(self, reactor_choice, t_end, T_reac, P_reac, mix, **kwargs) -> Tuple[dict, dict]: + """Perform a simulation of the reaction output + + Args: + reactor_choice: Choice of the reactor + t_end: How long of a simulation to run + T_reac: Temperature at which to perform the reaction + P_reac: Pressure at which to perform the simulation + mix: Initial concentrations of the reactants + kwargs: Options that are passed forward to the specific type of reactor being used and ODE solver + Returns: + - Results of the simulation + - Extensive details of the model performance + """ return self.reactor.run(reactor_choice, t_end, T_reac, P_reac, mix, **kwargs) class Uncertainty: # alternate name: why I hate pickle part 10 + """Computes the uncertainty bounds for parameters""" + def __init__(self, unc_type, rxnIdx, **kwargs): + """ + Args: + unc_type: Type of the + rxnIdx: Index of the reaction + **kwargs: + """ # self.gas = gas self.unc_type = unc_type self.rxnIdx = rxnIdx self.unc_dict = kwargs - - def _unc_fcn(self, x, uncVal, uncType): # uncertainty function + + def _unc_fcn(self, x, uncVal, uncType): # TODO (wardlt): Can be static or + """Compute the bounds given the value of a coefficient + + Args: + x: Assembled value of the coefficient + uncVal: Magnitude of the uncertainty + uncType: Type of the uncertainty (e.g., %, +, F) + Returns: + Bounds for the value + """ if np.isnan(uncVal): return [np.nan, np.nan] elif uncType == 'F': @@ -549,6 +585,8 @@ def _unc_fcn(self, x, uncVal, uncType): # uncertainty function return np.sort([x, x + uncVal], axis=0) elif uncType == '-': return np.sort([x - uncVal, x], axis=0) + else: + raise ValueError(f'Unsupported type {uncType}') def __call__(self, x=None): if self.unc_type == 'rate': @@ -570,4 +608,4 @@ def __call__(self, x=None): coef_val = coef_dict['resetVal'] unc_value = coef_dict['value'] unc_type = coef_dict['type'] - return self._unc_fcn(coef_val, unc_value, unc_type) \ No newline at end of file + return self._unc_fcn(coef_val, unc_value, unc_type) diff --git a/frhodo/calculate/optimize/fit_fcn.py b/frhodo/calculate/optimize/fit_fcn.py index dd40365..5fec164 100644 --- a/frhodo/calculate/optimize/fit_fcn.py +++ b/frhodo/calculate/optimize/fit_fcn.py @@ -72,7 +72,18 @@ def update_mech_coef_opt(mech, coef_opt, x): if mech_changed: mech.modify_reactions(mech.coeffs) # Update mechanism with new coefficients -def calculate_residuals(args_list): +def calculate_residuals(args_list): + """Evaluate the different between a shock experiment and the simulation given a list of parameters + + Args: + args_list: Tuple with arguments in the following order: + var: Variables being optimized + coef_opt: Description of the coefficients being optimized + x: Guess for these new parameters + shock: Description of the shock experiment + Returns: + ??? + """ def resid_func(t_offset, t_adjust, t_sim, obs_sim, t_exp, obs_exp, weights, obs_bounds=[], loss_alpha=2, loss_c=1, loss_penalty=True, scale='Linear', bisymlog_scaling_factor=1.0, DoF=1, opt_type='Residual', verbose=False): @@ -249,7 +260,20 @@ def calc_density(x, data, dim=1): # Using optimization vs least squares curve fit because y_range's change if time_offset != 0 class Fit_Fun: + """Fitness function to be minimized during optimization""" + def __init__(self, input_dict): + """ + + Args: + input_dict: Dictionary containing all input data + - parent: Link to the host GUI. Accesses optimization settings from here + - shocks2run: List of the shock experiments to compare against + - coef_opt: Dictionary describing which coefficients are being optimized + - rxn_rate_opttheir initial values, and the bounds + - mech: Link to the mechanism object being tuned + - signals: Link to the signals in the main application + """ self.parent = input_dict['parent'] self.shocks2run = input_dict['shocks2run'] self.data = self.parent.series.shock @@ -288,10 +312,20 @@ def __init__(self, input_dict): input_dict['opt_settings'] = self.opt_settings self.CheKiPEUQ_Frhodo_interface = CheKiPEUQ_Frhodo_interface(input_dict) - def __call__(self, s, optimizing=True): + def __call__(self, s, optimizing=True): + """Evaluate the fitness function + + Args: + s: + optimizing: Whether this function is being used by an optimizer. If not, return extra data + + Returns: + The objective function value if optimizing. Objective function value plus + + """ def append_output(output_dict, calc_resid_output): for key in calc_resid_output: - if key not in output_dict: + if key not in output_dict: # TODO (wardlt) replace with defaultdict output_dict[key] = [] output_dict[key].append(calc_resid_output[key]) @@ -318,7 +352,8 @@ def append_output(output_dict, calc_resid_output): display_ind_var = None display_observable = None - + + # Evaluate the simulations in parallel if we have > 1 if self.multiprocessing and len(self.shocks2run) > 1: args_list = ((var_dict, self.coef_opt, x, shock) for shock in self.shocks2run) calc_resid_outputs = self.pool.map(calculate_residuals, args_list) @@ -434,7 +469,14 @@ def calculate_obj_fcn(self, x, loss_resid, alpha, log_opt_rates, output_dict, ob return obj_fcn - def fit_all_coeffs(self, all_rates): + def fit_all_coeffs(self, all_rates): + """Adjust the input rates such that they are all consistent and within user-defined bounds + + Args: + all_rates: List of all the rates being optimized + Returns: + Updated set of reaction coefficients + """ coeffs = [] i = 0 for rxn_coef in self.rxn_coef_opt: @@ -456,4 +498,4 @@ def fit_all_coeffs(self, all_rates): i += len(T) - return coeffs \ No newline at end of file + return coeffs diff --git a/frhodo/calculate/optimize/mech_optimize.py b/frhodo/calculate/optimize/mech_optimize.py index 7cca45c..fb788f8 100644 --- a/frhodo/calculate/optimize/mech_optimize.py +++ b/frhodo/calculate/optimize/mech_optimize.py @@ -41,17 +41,22 @@ def __init__(self, parent): def start_threads(self): parent = self.parent + + # Prepare the name of the output file parent.path_set.optimized_mech() + # Controls for the display self.last_plot_timer = 0.0 self.time_between_plots = 0.0 # maximum update rate, updated based on time to plot + # Check whether we should be using multiprocessing parent.multiprocessing = parent.multiprocessing_box.isChecked() ## Check fit_coeffs #from optimize.fit_coeffs import debug #debug(parent.mech) + # Check if the optimization cannot start for some reason if parent.directory.invalid: parent.log.append('Invalid directory found\n') return @@ -72,8 +77,9 @@ def start_threads(self): continue self.shocks2run.append(shock) - - if len(self.shocks2run) == 0: # optimize current shock if nothing selected + + # Optimize current shock if nothing selected + if len(self.shocks2run) == 0: self.shocks2run = [parent.display_shock] else: if not parent.load_full_series_box.isChecked(): # TODO: May want to remove this limitation in future @@ -150,7 +156,8 @@ def start_threads(self): # Create Progress Bar # parent.create_progress_bar() - + + # Invoke the optimization to begin in a separate thread if not parent.abort: s = 'Optimization starting\n\n Iteration\t\t Objective Func\tBest Objetive Func' parent.log.append(s, alert=False) @@ -163,15 +170,16 @@ def _has_opt_rxns(self): return True return False - def _initialize_opt(self): # initialize various dictionaries for optimization + def _initialize_opt(self): + """Initialize various dictionaries for optimization""" self.coef_opt = self._set_coef_opt() self.rxn_coef_opt = self._set_rxn_coef_opt() self.rxn_rate_opt = self._set_rxn_rate_opt() - def _set_coef_opt(self): + def _set_coef_opt(self): + """Find which coefficients of kinetic models should be optimized""" mech = self.parent.mech coef_opt = [] - #for rxnIdx in range(mech.gas.n_reactions): # searches all rxns for rxnIdx, rxn in enumerate(mech.gas.reactions()): # searches all rxns if not mech.rate_bnds[rxnIdx]['opt']: continue # ignore fixed reactions @@ -186,8 +194,10 @@ def _set_coef_opt(self): return coef_opt def _set_rxn_coef_opt(self, min_T_range=500, min_P_range_factor=2): + """Get the initial value and bounds for the parameters being optimized""" mech = self.parent.mech rxn_coef_opt = [] + print(self.coef_opt) for coef in self.coef_opt: if len(rxn_coef_opt) == 0 or coef['rxnIdx'] != rxn_coef_opt[-1]['rxnIdx']: rxn_coef_opt.append(deepcopy(coef)) @@ -266,7 +276,7 @@ def _set_rxn_coef_opt(self, min_T_range=500, min_P_range_factor=2): rxn_coef['invT'] = [] rxn_coef['P'] = [] - # set condtions for upper and lower rates + # set conditions for upper and lower rates for coef_type in ['low_rate', 'high_rate']: n_coef = 0 for coef in rxn_coef['key']: diff --git a/frhodo/calculate/optimize/optimize_worker.py b/frhodo/calculate/optimize/optimize_worker.py index 4f2d900..145f4ef 100644 --- a/frhodo/calculate/optimize/optimize_worker.py +++ b/frhodo/calculate/optimize/optimize_worker.py @@ -14,6 +14,7 @@ from ..optimize.fit_fcn import initialize_parallel_worker, Fit_Fun from ..optimize.misc_fcns import rates +from frhodo.calculate.mech_fcns import Chemical_Mechanism import frhodo @@ -35,7 +36,19 @@ class Worker(QRunnable): ''' - def __init__(self, parent, shocks2run, mech, coef_opt, rxn_coef_opt, rxn_rate_opt, *args, **kwargs): + def __init__(self, parent, shocks2run, mech: Chemical_Mechanism, + coef_opt, rxn_coef_opt, rxn_rate_opt, *args, **kwargs): + """ + Args: + parent: Attachment to parent QApplication + shocks2run: Data of the shocks to be compared against + mech: Connection to the :class:`~Chemical_Mechanism` holding reaction data + coef_opt: List of the coefficients being optimized, path to them in `mech` + rxn_coef_opt: The extracted coefficients and bounds for each parameter in `ceof_opt` + rxn_rate_opt: A parsimonious form of the coefficients and bounds, usable by optimization codes + *args: Unused? + **kwargs: Unused + """ super(Worker, self).__init__() # Store constructor arguments (re-used for processing) self.parent = parent @@ -57,7 +70,7 @@ def __init__(self, parent, shocks2run, mech, coef_opt, rxn_coef_opt, rxn_rate_op def _initialize(self): mech = self.mech - # Calculate initial rate scalers + # Calculate initial rate scalars lb, ub = self.rxn_rate_opt['bnds'].values() self.s = rates(self.rxn_coef_opt, mech) - self.rxn_rate_opt['x0'] # this initializes from current GUI settings @@ -76,7 +89,10 @@ def trim_shocks(self): # trim shocks from zero weighted data shock['abs_uncertainties_trim'] = shock['abs_uncertainties'][exp_bounds,:] def optimize_coeffs(self): + """Perform the optimization""" parent = self.parent + + # Create a pool where each worker has a copy of the full reaction mechanism dictionary pool = mp.Pool(processes=parent.max_processors, initializer=initialize_parallel_worker, initargs=(parent.mech.reset_mech, parent.mech.thermo_coeffs, parent.mech.coeffs, parent.mech.coeffs_bnds, @@ -87,9 +103,11 @@ def optimize_coeffs(self): input_dict = {'parent': parent, 'pool': pool, 'mech': self.mech, 'shocks2run': self.shocks2run, 'coef_opt': self.coef_opt, 'rxn_coef_opt': self.rxn_coef_opt, 'rxn_rate_opt': self.rxn_rate_opt, 'multiprocessing': parent.multiprocessing, 'signals': self.signals} - + + # Create the function that will be called during optimization Scaled_Fit_Fun = Fit_Fun(input_dict) - def eval_fun(s, grad=None): + def eval_fun(s, grad=None): + """Wrapper for evaluation function. Handles aborting when desired""" if self.__abort: parent.optimize_running = False self.signals.log.emit('\nOptimization aborted') @@ -225,6 +243,14 @@ class WorkerSignals(QObject): class Optimize: def __init__(self, obj_fcn, x0, bnds, opt_options, Scaled_Fit_Fun): + """ + Args: + obj_fcn: Objection function wrapper that invokes `Scaled_Fit_Fun` and includes an "abort" feature + x0: Initial guess + bnds: + opt_options: + Scaled_Fit_Fun: Fitness function used by `obj_fun` + """ self.obj_fcn = obj_fcn self.x0 = x0 self.bnds = bnds @@ -295,7 +321,8 @@ def nlopt(self, x0, bnds, options): x = opt.optimize(x0) # optimize! #s = parent.optimize.HoF['s'] - + + # Invoke the code one last time with the proposed solution obj_fcn, x, shock_output = self.Scaled_Fit_Fun(x, optimizing=False) if nlopt.SUCCESS > 0: @@ -374,23 +401,26 @@ def rbfopt(self, x0, bnds, options): # noisy, cheap function option. supports d timer_start = timer() + # Determine how to stop the optimization if options['stop_criteria_type'] == 'Iteration Maximum': max_eval = int(options['stop_criteria_val']) max_time = 1E30 elif options['stop_criteria_type'] == 'Maximum Time [min]': max_eval = 10000 # will need to check if rbfopt changes based on iteration max_time = options['stop_criteria_val']*60 + else: + raise ValueError(f'Stopping condition "{options["stop_criteria_type"]}" not yet implemented') - var_type = ['R']*np.size(x0) # specifies that all variables are continious + var_type = ['R']*np.size(x0) # specifies that all variables are continuous output = {'success': False, 'message': []} - # Intialize and report any problems to log, not to console window + # Initialize and report any problems to log, not to console window stdout = io.StringIO() stderr = io.StringIO() with contextlib.redirect_stderr(stderr): with contextlib.redirect_stdout(stdout): bb = rbfopt.RbfoptUserBlackBox(np.size(x0), np.array(bnds[0]), np.array(bnds[1]), - np.array(var_type), self.obj_fcn) + np.array(var_type), self.obj_fcn) settings = rbfopt.RbfoptSettings(max_iterations=max_eval, max_evaluations=max_eval, max_cycles=1E30, @@ -414,4 +444,4 @@ def rbfopt(self, x0, bnds, options): # noisy, cheap function option. supports d res = {'x': x, 'shock': shock_output, 'fval': obj_fcn, 'nfev': evalcount + fast_evalcount, 'success': output['success'], 'message': output['message'], 'time': timer() - timer_start} - return res \ No newline at end of file + return res diff --git a/frhodo/calculate/reactors.py b/frhodo/calculate/reactors.py index bd4480c..07c79c4 100644 --- a/frhodo/calculate/reactors.py +++ b/frhodo/calculate/reactors.py @@ -114,6 +114,7 @@ def __call__(self, idx=None, units='CGS'): # units must be 'CGS' or 'SI' class Simulation_Result: + """Storage for results of a reactor simulation""" def __init__(self, num=None, states=None, reactor_vars=[]): self.states = states self.all_var = all_var @@ -175,6 +176,7 @@ def finalize(self, success, ind_var, observable, units='CGS'): class Reactor: + """Simulates all types of reactors""" def __init__(self, mech): self.mech = mech self.ODE_success = False diff --git a/frhodo/config_io.py b/frhodo/config_io.py index 070b38e..b0a4cb2 100644 --- a/frhodo/config_io.py +++ b/frhodo/config_io.py @@ -202,6 +202,8 @@ def from_yaml(self, src=None): class GUI_settings: + """Link between the UI display and the underlying :class:`GUI_config` data""" + def __init__(self, parent): self.parent = parent # Need to find a better solution than passing parent self.cfg_io = GUI_Config() @@ -317,7 +319,8 @@ def set_box(box, val): def save(self, save_all=False): parent = self.parent - + + # Accesses sepcific settings = {'directory': self.cfg['Directory Settings'], 'exp': self.cfg['Experiment Settings'], 'reactor': self.cfg['Reactor Settings'], @@ -410,4 +413,4 @@ def save(self, save_all=False): # gui_cfg.dump(gui_cfg.settings, sys.stdout) gui_cfg.to_yaml(path['default_config']) - gui_cfg.from_yaml(path['default_config']) \ No newline at end of file + gui_cfg.from_yaml(path['default_config']) diff --git a/frhodo/ipopt/ipopt-win64/ipopt.exe b/frhodo/ipopt/ipopt-win64/ipopt.exe old mode 100644 new mode 100755 diff --git a/frhodo/main.py b/frhodo/main.py index 29c8f57..09b988e 100644 --- a/frhodo/main.py +++ b/frhodo/main.py @@ -19,10 +19,10 @@ from .version import __version__ from . import appdirs, options_panel_widgets, sim_explorer_widget from . import settings, config_io, save_widget, error_window, help_menu - + if os.environ['QT_API'] == 'pyside2': # Silence warning: "Qt WebEngine seems to be initialized from a plugin." QApplication.setAttribute(QtCore.Qt.AA_ShareOpenGLContexts) - + # Handle high resolution displays: Minimum recommended resolution 1280 x 960 if hasattr(QtCore.Qt, 'AA_EnableHighDpiScaling'): QApplication.setAttribute(QtCore.Qt.AA_EnableHighDpiScaling, True) @@ -35,39 +35,48 @@ # set appdata folder using AppDirs library (but just using the source code file) dirs = appdirs.AppDirs(appname='Frhodo', roaming=True, appauthor=False) path['appdata'] = pathlib.Path(dirs.user_config_dir) -path['appdata'].mkdir(parents=True, exist_ok=True) # Make path if it doesn't exist +path['appdata'].mkdir(parents=True, exist_ok=True) # Make path if it doesn't exist shut_down = {'bool': False} + class Main(QMainWindow): def __init__(self, app, path): super().__init__() self.app = app - self.path_set = settings.Path(self, path) + + # Create an object that will facilitate loading data from disk + self.path_set = settings.Path(self, path) # Also loads in mechanism data according to default configuration uic.loadUi(str(self.path['main']/'UI'/'main_window.ui'), self) # ~0.4 sec self.splitter.moveSplitter(0, 1) # moves splitter 0 as close to 1 as possible self.setWindowIcon(QtGui.QIcon(str(self.path['main']/'UI'/'graphics'/'main_icon.png'))) - + # Start threadpools self.threadpool = QtCore.QThreadPool() self.threadpool.setMaxThreadCount(2) # Sets thread count to 1 (1 for gui - this is implicit, 1 for calc) - + # Set selected tabs for tab_widget in [self.option_tab_widget, self.plot_tab_widget]: tab_widget.setCurrentIndex(0) - + # Set Clipboard self.clipboard = QApplication.clipboard() - - self.var = {'reactor': {'t_unit_conv': 1}} + + # Initialize the panels that manage experimental data + self.var = {'reactor': {'t_unit_conv': 1}} # `var`, like `parent`, is a catch-all variable storage self.SIM = reactors.Simulation_Result() self.mech_loaded = False self.run_block = True self.convert_units = convert_units.Convert_Units(self) self.series = settings.series(self) - + + # Initialize the panels that display optimization results self.sim_explorer = sim_explorer_widget.SIM_Explorer_Widgets(self) self.plot = plot(self) + + # Create the options panel(s) and loads in the experimental data options_panel_widgets.Initialize(self) + + # Create storage for the chemical mechanism data self.mech = mech_fcns.Chemical_Mechanism() # Setup save sim @@ -79,38 +88,42 @@ def __init__(self, app, path): sys.exit() else: self.show() - self.app.processEvents() # allow everything to draw properly - + self.app.processEvents() # allow everything to draw properly + # Initialize Settings self.initialize_settings() # ~ 4 sec # Setup help menu self.version = __version__ help_menu.HelpMenu(self) - + def initialize_settings(self): # TODO: Solving for loaded shock twice + """Load in the user-defined settings from configuration files on disk""" msgBox = MessageWindow(self, 'Loading...') self.app.processEvents() self.var['old_shock_choice'] = self.var['shock_choice'] = 1 - + + # Load in the UI selections for the last box self.user_settings = config_io.GUI_settings(self) self.user_settings.load() - + + # Whether to load >1 files from the self.load_full_series = self.load_full_series_box.isChecked() # TODO: Move to somewhere else? - + # load previous paths if file in path, can be accessed, and is a file - if ('path_file' in self.path and os.access(self.path['path_file'], os.R_OK) and + if ('path_file' in self.path and os.access(self.path['path_file'], os.R_OK) and self.path['path_file'].is_file()): - + self.path_set.load_dir_file(self.path['path_file']) # ~3.9 sec - + self.update_user_settings() self.run_block = False # Block multiple simulations from running during initialization self.run_single() # Attempt simulation after initialization completed msgBox.close() - + def load_mech(self, event = None): + """Load in the mechanism data from disk""" def mechhasthermo(mech_path): f = open(mech_path, 'r') while True: @@ -119,21 +132,21 @@ def mechhasthermo(mech_path): continue if 'ther' in line.split('!')[0].strip().lower(): return True - + if not line: break - + f.close() return False - + if self.mech_select_comboBox.count() == 0: return # if no items return, unsure if this is needed now - + # Specify mech file path self.path['mech'] = self.path['mech_main'] / str(self.mech_select_comboBox.currentText()) if not self.path['mech'].is_file(): # if it's not a file, then it was deleted self.path_set.mech() # update mech pulldown choices return - + # Check use thermo box viability if mechhasthermo(self.path['mech']): if self.thermo_select_comboBox.count() == 0: @@ -143,7 +156,7 @@ def mechhasthermo(mech_path): # Autoselect checkbox off if thermo exists in mech if self.sender() is None or 'use_thermo_file_box' not in self.sender().objectName(): self.use_thermo_file_box.blockSignals(True) # stop set from sending signal, causing double load - self.use_thermo_file_box.setChecked(False) + self.use_thermo_file_box.setChecked(False) self.use_thermo_file_box.blockSignals(False) # allow signals again else: self.use_thermo_file_box.blockSignals(True) # stop set from sending signal, causing double load @@ -156,7 +169,7 @@ def mechhasthermo(mech_path): self.thermo_select_comboBox.setEnabled(True) else: self.thermo_select_comboBox.setDisabled(True) - + # Specify thermo file path if self.use_thermo_file_box.isChecked(): if self.thermo_select_comboBox.count() > 0: @@ -166,22 +179,24 @@ def mechhasthermo(mech_path): return else: self.path['thermo'] = None - + # Initialize Mechanism self.log.clear([]) # Clear log when mechanism changes to avoid log errors about prior mech mech_load_output = self.mech.load_mechanism(self.path) self.log.append(mech_load_output['message'], alert=not mech_load_output['success']) + self.mech_loaded = mech_load_output['success'] - + if not mech_load_output['success']: # if error: update species and return self.mix.update_species() self.log._blink(True) # updating_species is causing blink to stop due to successful shock calculation return - + # Initialize tables and trees self.tree.set_trees(self.mech) self.mix.update_species() # this was commented out, could be because multiple calls to solver from update_mix / setItems - + + # Update the appropriate display tab tabIdx = self.plot_tab_widget.currentIndex() tabText = self.plot_tab_widget.tabText(tabIdx) if tabText == 'Signal/Sim': @@ -190,34 +205,34 @@ def mechhasthermo(mech_path): self.plot.observable_widget.widget['main_parameter'].currentIndexChanged[str].emit(observable) elif tabText == 'Sim Explorer': # TODO: This gets called twice? self.sim_explorer.update_all_main_parameter() - + def shock_choice_changed(self, event): if 'exp_main' in self.directory.invalid: # don't allow shock change if problem with exp directory return - + self.var['old_shock_choice'] = self.var['shock_choice'] self.var['shock_choice'] = event - + self.shockRollingList = ['P1', 'u1'] # reset rolling list self.rxn_change_history = [] # reset tracking of rxn numbers changed - + if not self.optimize_running: self.log.clear([]) self.series.change_shock() # link display_shock to correct set and - + def update_user_settings(self, event = None): # This is one is located on the Files tab shock = self.display_shock self.series.set('series_name', self.exp_series_name_box.text()) - + t_unit_conv = self.var['reactor']['t_unit_conv'] if self.time_offset_box.value()*t_unit_conv != shock['time_offset']: # if values are different self.series.set('time_offset', self.time_offset_box.value()*t_unit_conv) if hasattr(self.mech_tree, 'rxn'): # checked if copy valid in function self.tree._copy_expanded_tab_rates() # copy rates and time offset - + self.var['time_unc'] = self.time_unc_box.value()*t_unit_conv - + # self.user_settings.save() # saves settings everytime a variable is changed if event is not None: @@ -235,45 +250,49 @@ def update_user_settings(self, event = None): for i in self.var: print('key: {:<14s} value: {:<16s}'.format(i, str(self.var[i]))) ''' - + def keyPressEvent(self, event): pass # THIS IS NOT FULLY FUNCTIONING # http://ftp.ics.uci.edu/pub/centos0/ics-custom-build/BUILD/PyQt-x11-gpl-4.7.2/doc/html/qkeyevent.html # print(event.modifiers(),event.text()) - + def run_single(self, event=None, t_save=None, rxn_changed=False): if self.run_block: return if not self.mech_loaded: return # if mech isn't loaded successfully, exit if not hasattr(self.mech_tree, 'rxn'): return # if mech tree not set up, exit - + shock = self.display_shock - + + # Get the conditions of the current reactor T_reac, P_reac, mix = shock['T_reactor'], shock['P_reactor'], shock['thermo_mix'] + + # Make sure the rate constants are update-to-date with the conditions on the table and specified shock self.tree.update_rates() - + # calculate all properties or observable by sending t_save tabIdx = self.plot_tab_widget.currentIndex() tabText = self.plot_tab_widget.tabText(tabIdx) if tabText == 'Sim Explorer': t_save = np.array([0]) - - SIM_kwargs = {'u_reac': shock['u2'], 'rho1': shock['rho1'], 'observable': self.display_shock['observable'], - 't_lab_save': t_save, 'sim_int_f': self.var['reactor']['sim_interp_factor'], - 'ODE_solver': self.var['reactor']['ode_solver'], + + # Formulate the output arguments for the + SIM_kwargs = {'u_reac': shock['u2'], 'rho1': shock['rho1'], 'observable': self.display_shock['observable'], + 't_lab_save': t_save, 'sim_int_f': self.var['reactor']['sim_interp_factor'], + 'ODE_solver': self.var['reactor']['ode_solver'], 'rtol': self.var['reactor']['ode_rtol'], 'atol': self.var['reactor']['ode_atol']} - + if '0d Reactor' in self.var['reactor']['name']: SIM_kwargs['solve_energy'] = self.var['reactor']['solve_energy'] SIM_kwargs['frozen_comp'] = self.var['reactor']['frozen_comp'] - - self.SIM, verbose = self.mech.run(self.var['reactor']['name'], self.var['reactor']['t_end'], + + self.SIM, verbose = self.mech.run(self.var['reactor']['name'], self.var['reactor']['t_end'], T_reac, P_reac, mix, **SIM_kwargs) - + if verbose['success']: - self.log._blink(False) + self.blink = self.log._blink(False) else: self.log.append(verbose['message']) - + if self.SIM is not None: self.plot.signal.update_sim(self.SIM.independent_var, self.SIM.observable, rxn_changed) if tabText == 'Sim Explorer': diff --git a/frhodo/mech_widget.py b/frhodo/mech_widget.py index 052858c..549a1aa 100644 --- a/frhodo/mech_widget.py +++ b/frhodo/mech_widget.py @@ -12,6 +12,7 @@ from qtpy import QtWidgets, QtGui, QtCore from . import misc_widget +from .calculate.mech_fcns import Chemical_Mechanism def silentSetValue(obj, value): @@ -74,18 +75,29 @@ def item_clicked(self, event): tree.expand(event) item.info['isExpanded'] = True - def set_trees(self, mech): + def set_trees(self, mech: Chemical_Mechanism): + """Update the widget display with a new chemical mechanism""" parent = self.parent() #parent.mech_tree.reset() - self.model.removeRows(0, self.model.rowCount()) + self.model.removeRows(0, self.model.rowCount()) # Clear the old version + + # Determine the display mode if 'Chemkin' in parent.tab_select_comboBox.currentText(): self.mech_tree_type = 'Chemkin' else: self.mech_tree_type = 'Bilbo' + + # Update the self.mech_tree_data = self._set_mech_tree_data(self.mech_tree_type, mech) self._set_mech_tree(self.mech_tree_data) def _set_mech_tree_data(self, selection, mech): + """Update the mechanism tree + + Args: + selection: Which type of mechanism file to use + mech: Mechanism data + """ def get_coef_abbreviation(coefName): if 'activation_energy' == coefName: return 'Ea' @@ -367,6 +379,11 @@ def getRateConst(parent, rxnNum, coef_key, coefName, value): parent.run_single(rxn_changed=True) def update_rates(self, rxnNum=None): + """Update the reaction rates _for the user-selected T/P/concentrations_? + + Args: + rxnNum: Specific reaction number to update + """ parent = self.parent() shock = parent.display_shock @@ -398,6 +415,7 @@ def update_rates(self, rxnNum=None): self._copy_expanded_tab_rates() def update_uncertainties(self, event=None, sender=None): + """Update the uncertainty box based on an event from the UI""" parent = self.parent() mech = parent.mech @@ -410,7 +428,7 @@ def update_uncertainties(self, event=None, sender=None): if 'coefName' in sender.info: # this means the coef unc was changed coefNum, coefName, coefAbbr = sender.info['coefNum'], sender.info['coefName'], sender.info['coefAbbr'] - # get correct uncertainty diction based on reaction type + # get correct uncertainty dictionary based on reaction type coef_key, bnds_key = keysFromBox(sender, mech) coefUncDict = mech.coeffs_bnds[rxnNum][bnds_key][coefName] @@ -466,6 +484,12 @@ def update_uncertainties(self, event=None, sender=None): parent.series.rate_bnds(parent.display_shock) def update_coef_rate_from_opt(self, coef_opt, x): + """Updated all coefficients given the coefficient dictionary and vector used in optimization + + Args: + coef_opt: Dictionary describing the parameters in `x` + x: Coefficients set by the optimizer + """ parent = self.parent() conv_type = 'Cantera2' + self.mech_tree_type diff --git a/frhodo/options_panel_widgets.py b/frhodo/options_panel_widgets.py index f5bb4a4..cb0d42a 100644 --- a/frhodo/options_panel_widgets.py +++ b/frhodo/options_panel_widgets.py @@ -4,6 +4,7 @@ import pathlib, os, sys from copy import deepcopy +from typing import Union, List import numpy as np from scipy.optimize import minimize @@ -144,7 +145,7 @@ def preset(self, selection): def select(self): parent = self.parent() - + key = '_'.join(self.sender().objectName().split("_")[:-1]) if 'path_file_load' in key: key = 'path_file' @@ -153,6 +154,7 @@ def select(self): dialog = 'select' type = self.sender().objectName().split("_")[-1] + if 'button' in type: description_text = eval('parent.' + key + '_box.placeholderText()') initial_dir = pathlib.Path.home() # set user as initial folder @@ -189,8 +191,8 @@ def fn(parent, text): return self.sender().setPlainText(text) self.QTextEdit_function(self.sender(), fn, parent, text) - parent.path[key] = pathlib.Path(text) - + parent.path[key] = pathlib.Path(text) + # select will modify box, this section is under if box to prevent double calling self.update_icons() if 'mech_main' in key and 'mech_main' not in self.invalid: # Mech path changed: update mech combobox @@ -973,6 +975,8 @@ def switch_unc_type(self): class Tables_Tab(QtCore.QObject): + """Table holding the thermodynamic and kinetic models""" + def __init__(self, parent): super().__init__(parent) parent = self.parent() @@ -1025,7 +1029,13 @@ def __init__(self, tab_widget, log_box, clear_log_button, copy_log_button): clear_log_button.clicked.connect(self.clear) copy_log_button.clicked.connect(self.copy) - def append(self, message, alert=True): + def append(self, message: Union[str, List[str]], alert: bool = True) -> None: + """ + + Args: + message: Message to be displayed + alert: Whether to blink the log tab in the GUI window + """ if isinstance(message, list): message = '\n'.join(message) diff --git a/frhodo/settings.py b/frhodo/settings.py index 40738bd..7c650f7 100644 --- a/frhodo/settings.py +++ b/frhodo/settings.py @@ -16,7 +16,10 @@ min_pos_system_value = (np.finfo(float).tiny*(1E20))**(1/2) max_pos_system_value = (np.finfo(float).max*(1E-20))**(1/2) + class Path: + """Facilitates loading data from disk""" + def __init__(self, parent, path): self.parent = parent self.loading_dir_file = False @@ -38,9 +41,11 @@ def __init__(self, parent, path): self.fs_watcher.directoryChanged.connect(self.mech) def mech(self): + """Load in the chemical mechanism files from the "mech_main" directory""" parent = self.parent - # Test for existance of mech folder and return if it doesn't + # Test for existence of mech folder and return if it doesn't + # TODO (wardlt): Remove after confirming this value is never used if parent.path['mech_main'].is_dir(): self.mech_main_exists = True else: @@ -255,6 +260,7 @@ def sim_output(self, var_name): # takes variable name and creates path for it return self.parent.path[var_name] def optimized_mech(self, file_out='opt_mech'): + """Define the filename of the optimized mechanism file""" parent = self.parent mech_name = parent.path['mech'].stem @@ -281,6 +287,7 @@ def optimized_mech(self, file_out='opt_mech'): return parent.path['Optimized_Mech_recast.mech'] def load_dir_file(self, file_path): + """Load the directories specified by the user""" parent = self.parent self.loading_dir_file = True self.config.read(file_path) @@ -988,9 +995,16 @@ def thermo_mix(self, shock=[]): # don't run if a species isn't in the mech or mech doesn't exist # elif not hasattr(parent.mech.gas, 'species_names'): return # elif species not in parent.mech.gas.species_names: return - - def rates(self, shock, rxnIdx=[]): # This resets and updates all rates in shock - if not self.parent.mech.isLoaded: return + + def rates(self, shock, rxnIdx=()): + """Resets and updates all rates in shock + + Args: + shock: Parameters of the shock to update rates for + rxnIdx: List of reactions to update + """ + if not self.parent.mech.isLoaded: + return mech = self.parent.mech mech_out = mech.set_TPX(shock['T_reactor'], shock['P_reactor'], shock['thermo_mix']) @@ -1000,7 +1014,7 @@ def rates(self, shock, rxnIdx=[]): # This resets and updates all rates in shock if not rxnIdx: rxnIdxRange = range(mech.gas.n_reactions) - elif not isinstance(rxns, list): # does not function right now + elif not isinstance(rxn, list): # does not function right now rxnIdxRange = [rxnIdx] else: rxnIdxRange = rxnIdx diff --git a/frhodo/version.py b/frhodo/version.py index 3b949b4..087b58e 100644 --- a/frhodo/version.py +++ b/frhodo/version.py @@ -1,4 +1,4 @@ -"""Single source of truth for package version""" -# https://packaging.python.org/en/latest/guides/single-sourcing-package-version/ - -__version__ = '1.3.1' +"""Single source of truth for package version""" +# https://packaging.python.org/en/latest/guides/single-sourcing-package-version/ + +__version__ = '1.3.1' From 35c5e4ebf0d06f5e61e4536390b89f87b84d2b7d Mon Sep 17 00:00:00 2001 From: Logan Ward Date: Mon, 22 Aug 2022 12:01:35 -0400 Subject: [PATCH 06/29] Added some initial unit tests --- .../optimize/CheKiPEUQ_local/plotting_functions.py | 2 +- .../optimize/CheKiPEUQ_local/test_pytest.py | 10 ++++------ tests/conftest.py | 13 +++++++++++++ tests/test_api.py | 7 +++++++ 4 files changed, 25 insertions(+), 7 deletions(-) create mode 100644 tests/conftest.py create mode 100644 tests/test_api.py diff --git a/frhodo/calculate/optimize/CheKiPEUQ_local/plotting_functions.py b/frhodo/calculate/optimize/CheKiPEUQ_local/plotting_functions.py index 17e2af3..c0b3899 100644 --- a/frhodo/calculate/optimize/CheKiPEUQ_local/plotting_functions.py +++ b/frhodo/calculate/optimize/CheKiPEUQ_local/plotting_functions.py @@ -1,6 +1,6 @@ import sys #sys.path.insert(0, '/mumpce/') -import CheKiPEUQ.mumpce.Project as mumpceProject +import CheKiPEUQ.mumpce.Project as mumpceProject import CheKiPEUQ.mumpce.solution as mumpceSolution import numpy as np import matplotlib diff --git a/frhodo/calculate/optimize/CheKiPEUQ_local/test_pytest.py b/frhodo/calculate/optimize/CheKiPEUQ_local/test_pytest.py index 05b3f6c..68e0e06 100644 --- a/frhodo/calculate/optimize/CheKiPEUQ_local/test_pytest.py +++ b/frhodo/calculate/optimize/CheKiPEUQ_local/test_pytest.py @@ -1,13 +1,11 @@ # Type 'pytest' in the current directory to run tests. # EAW 2020/01/17 -import plotting_functions -from plotting_functions import plotting_functions -import UserInput_ODE_KIN_BAYES_SG_EW as UserInput -import run_me -from run_me import ip +# LWard - Broken as of 22Aug22 +# from .plotting_functions import plotting_functions_class +# from . import UserInput def test_mumpce_plots(): - plot_object = plotting_functions() + plot_object = plotting_functions_class() assert plot_object.mumpce_plots(model_parameter_info = UserInput.model_parameter_info, active_parameters = UserInput.active_parameters, pairs_of_parameter_indices = UserInput.pairs_of_parameter_indices, posterior_mu_vector = UserInput.posterior_mu_vector, posterior_cov_matrix = UserInput.posterior_cov_matrix, prior_mu_vector = UserInput.prior_mu_vector, prior_cov_matrix = UserInput.prior_cov_matrix, contour_settings_custom = UserInput.contour_settings_custom) == None assert exec(open("mumpce_custom_plotting_example.py").read()) == None diff --git a/tests/conftest.py b/tests/conftest.py new file mode 100644 index 0000000..302bf3a --- /dev/null +++ b/tests/conftest.py @@ -0,0 +1,13 @@ +from typing import Tuple + +from PyQt5.QtWidgets import QApplication +from pytest import fixture + +from frhodo.main import launch_gui, Main + + +@fixture() +def frhodo_app() -> Tuple[QApplication, Main]: + app, window = launch_gui() + yield app, window + app.quit() diff --git a/tests/test_api.py b/tests/test_api.py new file mode 100644 index 0000000..781f198 --- /dev/null +++ b/tests/test_api.py @@ -0,0 +1,7 @@ +"""Testing the API components of Frhodo""" + + +def test_launch(frhodo_app): + """Make sure we can launch Frhodo""" + app, window = frhodo_app + assert window.isVisible() From 1d1e97c0e061cd1652b2e25a05337fd825aa6b8d Mon Sep 17 00:00:00 2001 From: Logan Ward Date: Mon, 22 Aug 2022 13:11:28 -0400 Subject: [PATCH 07/29] Marked the old tests as xfail Not yet sure when these worked last --- frhodo/calculate/optimize/CheKiPEUQ_local/test_pytest.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/frhodo/calculate/optimize/CheKiPEUQ_local/test_pytest.py b/frhodo/calculate/optimize/CheKiPEUQ_local/test_pytest.py index 68e0e06..ddcdcb8 100644 --- a/frhodo/calculate/optimize/CheKiPEUQ_local/test_pytest.py +++ b/frhodo/calculate/optimize/CheKiPEUQ_local/test_pytest.py @@ -4,11 +4,16 @@ # from .plotting_functions import plotting_functions_class # from . import UserInput +from pytest import mark + + +@mark.xfail def test_mumpce_plots(): plot_object = plotting_functions_class() assert plot_object.mumpce_plots(model_parameter_info = UserInput.model_parameter_info, active_parameters = UserInput.active_parameters, pairs_of_parameter_indices = UserInput.pairs_of_parameter_indices, posterior_mu_vector = UserInput.posterior_mu_vector, posterior_cov_matrix = UserInput.posterior_cov_matrix, prior_mu_vector = UserInput.prior_mu_vector, prior_cov_matrix = UserInput.prior_cov_matrix, contour_settings_custom = UserInput.contour_settings_custom) == None assert exec(open("mumpce_custom_plotting_example.py").read()) == None +@mark.xfail def MH_then_scatterplot(): parseUserInputParameters() ip_object = ip() From 985baaf38895abc55af8a2ef5ad87b0bd82f52f1 Mon Sep 17 00:00:00 2001 From: Logan Ward Date: Mon, 22 Aug 2022 13:16:36 -0400 Subject: [PATCH 08/29] Have tests launch Frhodo headless --- frhodo/main.py | 12 +++++++++--- tests/conftest.py | 4 +++- 2 files changed, 12 insertions(+), 4 deletions(-) diff --git a/frhodo/main.py b/frhodo/main.py index 09b988e..bf966c2 100644 --- a/frhodo/main.py +++ b/frhodo/main.py @@ -7,7 +7,7 @@ import os, sys, platform, multiprocessing, pathlib, ctypes # os.environ['QT_API'] = 'pyside2' # forces pyside2 -from typing import Tuple +from typing import Tuple, Optional, List from qtpy.QtWidgets import QMainWindow, QApplication, QMessageBox from qtpy import uic, QtCore, QtGui @@ -309,20 +309,26 @@ def run_single(self, event=None, t_save=None, rxn_changed=False): # assert False -def launch_gui() -> Tuple[QApplication, Main]: +def launch_gui(args: Optional[List[str]] = None) -> Tuple[QApplication, Main]: """Launch the GUI + Args: + args: Arguments to pass to the QApplication. Use ``sys.argv`` by default Returns: - The QApplication instance - Link to the main window """ + # Get the default argument if none specified + if args is None: + args = sys.argv.copy() + if platform.system() == 'Windows': # this is required for pyinstaller on windows multiprocessing.freeze_support() if getattr(sys, 'frozen', False): # if frozen minimize console immediately ctypes.windll.user32.ShowWindow(ctypes.windll.kernel32.GetConsoleWindow(), 0) - app = QApplication(sys.argv) + app = QApplication(args) sys.excepthook = error_window.excepthookDecorator(app, path, shut_down) main_window = Main(app, path) diff --git a/tests/conftest.py b/tests/conftest.py index 302bf3a..569a639 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,3 +1,4 @@ +import sys from typing import Tuple from PyQt5.QtWidgets import QApplication @@ -8,6 +9,7 @@ @fixture() def frhodo_app() -> Tuple[QApplication, Main]: - app, window = launch_gui() + # Launch the application headless + app, window = launch_gui(['frhodo', '-platform', 'offscreen']) yield app, window app.quit() From 89bb77a783711387ab8f5735db7c4250b2d035c5 Mon Sep 17 00:00:00 2001 From: Logan Ward Date: Mon, 22 Aug 2022 13:49:16 -0400 Subject: [PATCH 09/29] First GUI function: load data from specified path Also adds ability to launch without load past files --- frhodo/api/__init__.py | 36 ++++++++++++++++++++++++++++++++ frhodo/config_io.py | 5 +++++ frhodo/main.py | 37 +++++++++++++++++++++++++-------- frhodo/options_panel_widgets.py | 14 ++++++------- frhodo/settings.py | 29 ++++++++++++++++++++------ tests/conftest.py | 10 +++++++-- tests/test_api.py | 11 ++++++++++ 7 files changed, 118 insertions(+), 24 deletions(-) create mode 100644 frhodo/api/__init__.py diff --git a/frhodo/api/__init__.py b/frhodo/api/__init__.py new file mode 100644 index 0000000..176464b --- /dev/null +++ b/frhodo/api/__init__.py @@ -0,0 +1,36 @@ +"""Functions related to using Frhodo from a Python interface + +Many of these functions require launching the Frhodo GUI first. +using :meth:`~fhrodo.main.launch_gui` +""" + +from pathlib import Path + +from PyQt5.QtWidgets import QApplication + +from ..main import Main + + +def load_files(window: Main, + app: QApplication, + experiment_path: Path, + mechanism_path: Path, + output_path: Path): + """Load the input files for Frhodo + + Args: + window: Link to the main window + app: Link to the Qt application + experiment_path: Path to the experimental data + mechanism_path: Path to the mechanism files + output_path: Path to which simulation results should be stored + """ + + window.exp_main_box.setPlainText(str(experiment_path.resolve().absolute())) + window.mech_main_box.setPlainText(str(mechanism_path.resolve().absolute())) + window.sim_main_box.setPlainText(str(output_path.resolve().absolute())) + + # Trigger Frhodo to process these files + app.processEvents() + + return diff --git a/frhodo/config_io.py b/frhodo/config_io.py index b0a4cb2..449a70e 100644 --- a/frhodo/config_io.py +++ b/frhodo/config_io.py @@ -192,6 +192,11 @@ def to_yaml(self, dest=None): self.dump(out, configFile) def from_yaml(self, src=None): + """Read settings from a YAML file, if present + + Args: + src: Path to where the settings file should be + """ if src is None: return if not src.exists(): return diff --git a/frhodo/main.py b/frhodo/main.py index bf966c2..3fc1fe7 100644 --- a/frhodo/main.py +++ b/frhodo/main.py @@ -40,7 +40,14 @@ class Main(QMainWindow): - def __init__(self, app, path): + def __init__(self, app: QApplication, path: dict, skip_config: bool = True): + """Launch the Frhodo GUI application + + Args: + app: Link the QTApplication + path: Collection of useful paths + skip_config: Skip loading in the default configuration. Primarily useful if running Frhodo from API + """ super().__init__() self.app = app @@ -52,7 +59,7 @@ def __init__(self, app, path): # Start threadpools self.threadpool = QtCore.QThreadPool() - self.threadpool.setMaxThreadCount(2) # Sets thread count to 1 (1 for gui - this is implicit, 1 for calc) + self.threadpool.setMaxThreadCount(2) # Sets thread count to 1 (1 for gui - this is implicit, 1 for calc) # Set selected tabs for tab_widget in [self.option_tab_widget, self.plot_tab_widget]: @@ -91,14 +98,18 @@ def __init__(self, app, path): self.app.processEvents() # allow everything to draw properly # Initialize Settings - self.initialize_settings() # ~ 4 sec + self.initialize_settings(skip_config) # ~ 4 sec # Setup help menu self.version = __version__ help_menu.HelpMenu(self) - def initialize_settings(self): # TODO: Solving for loaded shock twice - """Load in the user-defined settings from configuration files on disk""" + def initialize_settings(self, skip_config: bool): # TODO: Solving for loaded shock twice + """Load in the user-defined settings from configuration files on disk + + Args: + skip_config: Skip loading the default configuration + """ msgBox = MessageWindow(self, 'Loading...') self.app.processEvents() @@ -106,7 +117,8 @@ def initialize_settings(self): # TODO: Solving for loaded shock twice # Load in the UI selections for the last box self.user_settings = config_io.GUI_settings(self) - self.user_settings.load() + if not skip_config: + self.user_settings.load() # Whether to load >1 files from the self.load_full_series = self.load_full_series_box.isChecked() # TODO: Move to somewhere else? @@ -309,11 +321,12 @@ def run_single(self, event=None, t_save=None, rxn_changed=False): # assert False -def launch_gui(args: Optional[List[str]] = None) -> Tuple[QApplication, Main]: +def launch_gui(args: Optional[List[str]] = None, fresh_gui: bool = False) -> Tuple[QApplication, Main]: """Launch the GUI Args: args: Arguments to pass to the QApplication. Use ``sys.argv`` by default + fresh_gui: Whether to skip loading the previous configurations Returns: - The QApplication instance - Link to the main window @@ -328,10 +341,16 @@ def launch_gui(args: Optional[List[str]] = None) -> Tuple[QApplication, Main]: if getattr(sys, 'frozen', False): # if frozen minimize console immediately ctypes.windll.user32.ShowWindow(ctypes.windll.kernel32.GetConsoleWindow(), 0) + # Make a copy of the default paths, so that multiple launches of the Frhodo do not interfere + # WardLT: I added this while building unit tests, which involve launching Frhodo many times + my_path = path.copy() + + # Launch the QT application app = QApplication(args) - sys.excepthook = error_window.excepthookDecorator(app, path, shut_down) + sys.excepthook = error_window.excepthookDecorator(app, my_path, shut_down) - main_window = Main(app, path) + # Create the Frhodo main window + main_window = Main(app, my_path, fresh_gui) return app, main_window diff --git a/frhodo/options_panel_widgets.py b/frhodo/options_panel_widgets.py index cb0d42a..20b610a 100644 --- a/frhodo/options_panel_widgets.py +++ b/frhodo/options_panel_widgets.py @@ -132,8 +132,8 @@ def __init__(self, parent): parent.load_full_series_box.stateChanged.connect(self.set_load_full_set) self.set_load_full_set() - self.x_icon = QtGui.QPixmap(str(parent.path['graphics']/'x_icon.png')) - self.check_icon = QtGui.QPixmap(str(parent.path['graphics']/'check_icon.png')) + self.x_icon = QtGui.QPixmap(str(parent.path['graphics'] / 'x_icon.png')) + self.check_icon = QtGui.QPixmap(str(parent.path['graphics'] / 'check_icon.png')) self.update_icons() def preset(self, selection): @@ -169,7 +169,7 @@ def select(self): initial_dir = parent.path[key].parents[0] # set path_file as initial folder # if initial_dir doesn't exist or can't be accessed, choose source folder - if not os.access(parent.path[key], os.R_OK) or not initial_dir.is_dir(): + if not os.access(parent.path[key], os.R_OK) or not initial_dir.is_dir(): initial_dir = parent.path['main'] path = QFileDialog.getOpenFileName(parent, description_text, str(initial_dir), 'ini (*.ini)') @@ -274,8 +274,8 @@ def update_icons(self, invalid=[]): # This also checks if paths are valid if key in self.invalid: eval('parent.' + key + '_label.setPixmap(self.x_icon)') - elif (key in parent.path and os.access(parent.path[key], os.R_OK) - and parent.path[key].is_dir() and str(parent.path[key]) != '.'): + elif (key in parent.path and os.access(parent.path[key], os.R_OK) + and parent.path[key].is_dir() and str(parent.path[key]) != '.'): eval('parent.' + key + '_label.setPixmap(self.check_icon)') else: @@ -563,7 +563,7 @@ def create_thermo_boxes(self, species=[]): species.insert(0, '') self.thermoSpecies_box = [] # create down_arrow_path with forward slashes as required by QT stylesheet url - down_arrow_path = '"' + str((self.parent().path['graphics']/'arrowdown.png').as_posix()) + '"' + down_arrow_path = '"' + str((self.parent().path['graphics'] / 'arrowdown.png').as_posix()) + '"' for row in range(self.table.rowCount()): self.thermoSpecies_box.append(misc_widget.SearchComboBox()) self.thermoSpecies_box[-1].addItems(species) @@ -632,7 +632,7 @@ def isValidRow(table, row): # if path_file exists and species_aliases exist and not loading preset, save aliases if save_species_alias or len(parent.series.current['species_alias']) > 0: - if parent.path['path_file'].is_file() and not parent.path_set.loading_dir_file: + if parent.path['path_file'].is_file() and not parent.path_set.loading_dir_file: parent.path_set.save_aliases(parent.path['path_file']) parent.series.thermo_mix() diff --git a/frhodo/settings.py b/frhodo/settings.py index 7c650f7..c2e9476 100644 --- a/frhodo/settings.py +++ b/frhodo/settings.py @@ -25,12 +25,12 @@ def __init__(self, parent, path): self.loading_dir_file = False parent.path = path - parent.path['graphics'] = parent.path['main']/'UI'/'graphics' + parent.path['graphics'] = parent.path['main'] / 'UI' / 'graphics' self.config = configparser.RawConfigParser() # Specify yaml files - parent.path['default_config'] = parent.path['appdata']/'default_config.yaml' - parent.path['Cantera_Mech'] = parent.path['appdata']/'generated_mech.yaml' + parent.path['default_config'] = parent.path['appdata'] / 'default_config.yaml' + parent.path['Cantera_Mech'] = parent.path['appdata'] / 'generated_mech.yaml' for key in ['default_config', 'Cantera_Mech']: if parent.path[key].exists(): # Check that file is readable and writable if not os.access(parent.path[key], os.R_OK) or not os.access(parent.path[key], os.W_OK): @@ -344,7 +344,7 @@ def set_watch_dir(self): if self.fs_watcher.directories(): self.fs_watcher.removePaths(self.fs_watcher.directories()) - if self.parent.path['mech_main'].is_dir(): + if self.parent.path['mech_main'].is_dir(): self.fs_watcher.addPath(str(self.parent.path['mech_main'])) @@ -626,6 +626,9 @@ def b_eval(x, k, x0): return res class series: + """Holds data from a collection of related shock experiments. + + Also provides functionality for switching between different shocks in displays""" def __init__(self, parent): self.parent = parent self.exp = experiment(parent) @@ -666,7 +669,9 @@ def initialize_shock(self): self.shock.append(shock) def _create_shock(self, num, shock_path): + """Create a placeholder for data from a shock experiment""" parent = self.parent + # TODO (wardlt): This complex dictionary is ripe for implementation as a dataclass shock = {'num': num, 'path': deepcopy(shock_path), 'include': False, 'series_name': self.name[-1], 'run_SIM': True, # TODO: make run_SIM a trigger to calculate or not @@ -735,7 +740,7 @@ def add_series(self): # need to think about what to do when mech ch self.path.append(deepcopy(parent.path['exp_main'])) self.name.append(parent.exp_series_name_box.text()) - self.shock_num.append(list(parent.path['shock'][:,0].astype(int))) + self.shock_num.append(list(parent.path['shock'][:, 0].astype(int))) self.species_alias.append({}) self.in_table.append(False) @@ -753,10 +758,16 @@ def added_to_table(self, n): # update if in table self.in_table[n] = True def clear_series(self, n): + """Delete a series from the collection + + Args: + n: Index of series to be deleted + """ del self.path[n], self.name[n], self.shock_num[n], self.species_alias[n] del self.shock[n], self.in_table[n] def clear_shocks(self): + """Remove all shocks from a certain series""" if self.parent.load_full_series: return self.update_idx() # in case series are changed and load full set not selected @@ -1027,7 +1038,12 @@ def rates(self, shock, rxnIdx=()): return shock['rate_val'] - def rate_bnds(self, shock): + def rate_bnds(self, shock: dict): + """Update the bonds on the reaction rates based on a certain shock experiment + + Args: + shock: Dictionary contain the shock information + """ if not self.parent.mech.isLoaded: return mech = self.parent.mech @@ -1064,6 +1080,7 @@ def load_full_series(self): self.parent.series_viewer.update() def change_shock(self): + """Charge the selected shock based on the choice in the GUI""" parent = self.parent # parent.user_settings.save(save_all = False) diff --git a/tests/conftest.py b/tests/conftest.py index 569a639..f0c0170 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,4 +1,4 @@ -import sys +from pathlib import Path from typing import Tuple from PyQt5.QtWidgets import QApplication @@ -7,9 +7,15 @@ from frhodo.main import launch_gui, Main +@fixture() +def example_dir() -> Path: + """Directory containing example input files""" + return Path(__file__).parent.parent / 'Example' + + @fixture() def frhodo_app() -> Tuple[QApplication, Main]: # Launch the application headless - app, window = launch_gui(['frhodo', '-platform', 'offscreen']) + app, window = launch_gui(['frhodo', '-platform', 'offscreen'], fresh_gui=True) yield app, window app.quit() diff --git a/tests/test_api.py b/tests/test_api.py index 781f198..5ba08f2 100644 --- a/tests/test_api.py +++ b/tests/test_api.py @@ -1,7 +1,18 @@ """Testing the API components of Frhodo""" +from frhodo.api import load_files def test_launch(frhodo_app): """Make sure we can launch Frhodo""" app, window = frhodo_app assert window.isVisible() + + +def test_load(frhodo_app, example_dir, tmp_path): + """Load loading in desired data""" + app, window = frhodo_app + load_files(window, app, + example_dir / 'Experiment', + example_dir / 'Mechanism', + tmp_path) + assert len(window.series.shock) == 1 From a764a9acc99dc64e358d53820d3fefc04262cd7a Mon Sep 17 00:00:00 2001 From: Logan Ward Date: Mon, 22 Aug 2022 15:36:30 -0400 Subject: [PATCH 10/29] Make an API driver class Avoid having to pass window and app to every function --- frhodo/api/__init__.py | 53 +++++++++++++++++++++++++----------------- tests/conftest.py | 5 ++-- tests/test_api.py | 20 +++++++--------- 3 files changed, 44 insertions(+), 34 deletions(-) diff --git a/frhodo/api/__init__.py b/frhodo/api/__init__.py index 176464b..287c179 100644 --- a/frhodo/api/__init__.py +++ b/frhodo/api/__init__.py @@ -11,26 +11,37 @@ from ..main import Main -def load_files(window: Main, - app: QApplication, - experiment_path: Path, - mechanism_path: Path, - output_path: Path): - """Load the input files for Frhodo - - Args: - window: Link to the main window - app: Link to the Qt application - experiment_path: Path to the experimental data - mechanism_path: Path to the mechanism files - output_path: Path to which simulation results should be stored - """ - - window.exp_main_box.setPlainText(str(experiment_path.resolve().absolute())) - window.mech_main_box.setPlainText(str(mechanism_path.resolve().absolute())) - window.sim_main_box.setPlainText(str(output_path.resolve().absolute())) +class FrhodoDriver: + """Driver the Frhodo GUI application - # Trigger Frhodo to process these files - app.processEvents() + Implements simple interfaces for common tasks, such as loading differnet datasets + or evaluating changes to different mechanisms + """ - return + def __init__(self, window: Main, app: QApplication): + """Create the driver given connection to a running instance of Frhodo + + Args: + window: Link to the main Frhodo window + app: Link to the Qt backend + """ + self.window = window + self.app = app + + def load_files(self, experiment_path: Path, + mechanism_path: Path, + output_path: Path): + """Load the input files for Frhodo + + Args: + experiment_path: Path to the experimental data + mechanism_path: Path to the mechanism files + output_path: Path to which simulation results should be stored + """ + + self.window.exp_main_box.setPlainText(str(experiment_path.resolve().absolute())) + self.window.mech_main_box.setPlainText(str(mechanism_path.resolve().absolute())) + self.window.sim_main_box.setPlainText(str(output_path.resolve().absolute())) + + # Trigger Frhodo to process these files + self.app.processEvents() diff --git a/tests/conftest.py b/tests/conftest.py index f0c0170..5f4aa62 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -4,6 +4,7 @@ from PyQt5.QtWidgets import QApplication from pytest import fixture +from frhodo.api import FrhodoDriver from frhodo.main import launch_gui, Main @@ -14,8 +15,8 @@ def example_dir() -> Path: @fixture() -def frhodo_app() -> Tuple[QApplication, Main]: +def frhodo_driver() -> FrhodoDriver: # Launch the application headless app, window = launch_gui(['frhodo', '-platform', 'offscreen'], fresh_gui=True) - yield app, window + yield FrhodoDriver(window, app) app.quit() diff --git a/tests/test_api.py b/tests/test_api.py index 5ba08f2..a7cbae7 100644 --- a/tests/test_api.py +++ b/tests/test_api.py @@ -1,18 +1,16 @@ """Testing the API components of Frhodo""" -from frhodo.api import load_files -def test_launch(frhodo_app): +def test_launch(frhodo_driver): """Make sure we can launch Frhodo""" - app, window = frhodo_app - assert window.isVisible() + assert frhodo_driver.window.isVisible() -def test_load(frhodo_app, example_dir, tmp_path): +def test_load(frhodo_driver, example_dir, tmp_path): """Load loading in desired data""" - app, window = frhodo_app - load_files(window, app, - example_dir / 'Experiment', - example_dir / 'Mechanism', - tmp_path) - assert len(window.series.shock) == 1 + frhodo_driver.load_files( + example_dir / 'Experiment', + example_dir / 'Mechanism', + tmp_path + ) + assert len(frhodo_driver.window.series.shock) == 1 From 6e9a3e440379802f7a62a858d5c9e0349beb9432 Mon Sep 17 00:00:00 2001 From: Logan Ward Date: Tue, 23 Aug 2022 11:50:30 -0400 Subject: [PATCH 11/29] Add function for returning the experimental data --- frhodo/api/__init__.py | 36 ++++++++++++++++++++++++++++++++++++ tests/conftest.py | 13 ++++++------- tests/test_api.py | 14 +++++++++++++- 3 files changed, 55 insertions(+), 8 deletions(-) diff --git a/frhodo/api/__init__.py b/frhodo/api/__init__.py index 287c179..a1b4738 100644 --- a/frhodo/api/__init__.py +++ b/frhodo/api/__init__.py @@ -5,7 +5,9 @@ """ from pathlib import Path +from typing import List +import numpy as np from PyQt5.QtWidgets import QApplication from ..main import Main @@ -45,3 +47,37 @@ def load_files(self, experiment_path: Path, # Trigger Frhodo to process these files self.app.processEvents() + + @property + def n_shocks(self): + """Number of shock experiments loaded""" + n_series = len(self.window.series.shock) + assert n_series <= 1, "We only support one series" + if n_series == 0: + return 0 + return len(self.window.series.shock[0]) + + def _select_shock(self, n: int): + """Change which shock experiment is running + + Args: + n: Which shock experiment to evaluate + """ + + self.window.shock_choice_box.setValue(n + 1) + self.app.processEvents() + + def get_observables(self) -> List[np.ndarray]: + """Get the observable data from each shock experiment + + Returns: + List of experimental data arrays where each is a 2D + array with the first column is the time and second + is the observable + """ + # Loop over each shock + output = [] + for i in range(self.n_shocks): + self._select_shock(i) + output.append(self.window.display_shock['exp_data']) + return output diff --git a/tests/conftest.py b/tests/conftest.py index 5f4aa62..d4c92a8 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,11 +1,13 @@ from pathlib import Path -from typing import Tuple -from PyQt5.QtWidgets import QApplication from pytest import fixture from frhodo.api import FrhodoDriver -from frhodo.main import launch_gui, Main +from frhodo.main import launch_gui + +# Launch the application headless and without reading the configuration from a past run +app, window = launch_gui(['frhodo', '-platform', 'offscreen'], fresh_gui=True) +driver = FrhodoDriver(window, app) @fixture() @@ -16,7 +18,4 @@ def example_dir() -> Path: @fixture() def frhodo_driver() -> FrhodoDriver: - # Launch the application headless - app, window = launch_gui(['frhodo', '-platform', 'offscreen'], fresh_gui=True) - yield FrhodoDriver(window, app) - app.quit() + return driver diff --git a/tests/test_api.py b/tests/test_api.py index a7cbae7..a563873 100644 --- a/tests/test_api.py +++ b/tests/test_api.py @@ -13,4 +13,16 @@ def test_load(frhodo_driver, example_dir, tmp_path): example_dir / 'Mechanism', tmp_path ) - assert len(frhodo_driver.window.series.shock) == 1 + assert frhodo_driver.n_shocks == 1 + + +def test_observables(frhodo_driver, example_dir, tmp_path): + frhodo_driver.load_files( + example_dir / 'Experiment', + example_dir / 'Mechanism', + tmp_path + ) + runs = frhodo_driver.get_observables() + assert len(runs) == 1 + assert runs[0].ndim == 2 + assert runs[0].shape[1] == 2 From d30e1831b8fd6ec18e3a69c614a1c75939e308a5 Mon Sep 17 00:00:00 2001 From: Logan Ward Date: Tue, 23 Aug 2022 13:29:38 -0400 Subject: [PATCH 12/29] Add a utility for running simulations Uses a side-effect in the GUI where a simulation is ran whenever you change the simulation being displayed --- frhodo/api/__init__.py | 36 ++++++++++++++++++++++++++++++++---- tests/test_api.py | 28 ++++++++++++++++++---------- 2 files changed, 50 insertions(+), 14 deletions(-) diff --git a/frhodo/api/__init__.py b/frhodo/api/__init__.py index a1b4738..2441f7d 100644 --- a/frhodo/api/__init__.py +++ b/frhodo/api/__init__.py @@ -58,12 +58,15 @@ def n_shocks(self): return len(self.window.series.shock[0]) def _select_shock(self, n: int): - """Change which shock experiment is running + """Change which shock experiment is being displayed and simulated Args: n: Which shock experiment to evaluate """ + # Check if it is in the list + assert self.n_shocks > 0 + assert any(n == x['num'] for x in self.window.series.shock[0]) self.window.shock_choice_box.setValue(n + 1) self.app.processEvents() @@ -75,9 +78,34 @@ def get_observables(self) -> List[np.ndarray]: array with the first column is the time and second is the observable """ + + if self.n_shocks == 0: + return [] + # Loop over each shock output = [] - for i in range(self.n_shocks): - self._select_shock(i) - output.append(self.window.display_shock['exp_data']) + for shock in self.window.series.shock[0]: + output.append(shock['exp_data']) + return output + + def run_simulations(self) -> List[np.ndarray]: + """Run the simulation for each of the observed + + Returns: + List of simulated data arrays where each is a 2D + array with the first column is the time and second + is the simulated observable + """ + + # Loop over all shocks + output = [] + for shock in self.window.series.shock[0]: + # We force the simulation by changing the shock index in the GUI + # That could be an issue if the behavior of the GUI changes + self._select_shock(shock['num']) + assert self.window.SIM.success, 'Simulation failed' + output.append(np.stack([ + self.window.SIM.independent_var, + self.window.SIM.observable + ], axis=1)) return output diff --git a/tests/test_api.py b/tests/test_api.py index a563873..cf1c1ef 100644 --- a/tests/test_api.py +++ b/tests/test_api.py @@ -1,5 +1,13 @@ """Testing the API components of Frhodo""" +def _load_example(frhodo_driver, example_dir, tmp_path): + """Set up the driver with a specific problem case""" + frhodo_driver.load_files( + example_dir / 'Experiment', + example_dir / 'Mechanism', + tmp_path + ) + def test_launch(frhodo_driver): """Make sure we can launch Frhodo""" @@ -8,21 +16,21 @@ def test_launch(frhodo_driver): def test_load(frhodo_driver, example_dir, tmp_path): """Load loading in desired data""" - frhodo_driver.load_files( - example_dir / 'Experiment', - example_dir / 'Mechanism', - tmp_path - ) + _load_example(frhodo_driver, example_dir, tmp_path) assert frhodo_driver.n_shocks == 1 def test_observables(frhodo_driver, example_dir, tmp_path): - frhodo_driver.load_files( - example_dir / 'Experiment', - example_dir / 'Mechanism', - tmp_path - ) + _load_example(frhodo_driver, example_dir, tmp_path) runs = frhodo_driver.get_observables() assert len(runs) == 1 assert runs[0].ndim == 2 assert runs[0].shape[1] == 2 + + +def test_simulate(frhodo_driver, example_dir, tmp_path): + _load_example(frhodo_driver, example_dir, tmp_path) + runs = frhodo_driver.run_simulations() + assert len(runs) == 1 + assert runs[0].ndim == 2 + assert runs[0].shape[1] == 2 From fa6f4fa653f6cd2b4ab6a66d5b141a19e85ba00a Mon Sep 17 00:00:00 2001 From: Logan Ward Date: Tue, 23 Aug 2022 15:25:44 -0400 Subject: [PATCH 13/29] Bug fix: We must call always call `.load()` Even if we're avoiding loading the default configuration --- frhodo/main.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/frhodo/main.py b/frhodo/main.py index 3fc1fe7..b0b4c3e 100644 --- a/frhodo/main.py +++ b/frhodo/main.py @@ -117,8 +117,9 @@ def initialize_settings(self, skip_config: bool): # TODO: Solving for loaded sh # Load in the UI selections for the last box self.user_settings = config_io.GUI_settings(self) - if not skip_config: - self.user_settings.load() + if skip_config: + self.path['default_config'] = None # Will result in configuration loading to fail + self.user_settings.load() # Whether to load >1 files from the self.load_full_series = self.load_full_series_box.isChecked() # TODO: Move to somewhere else? From 8982ca6ae57cc6fe74a0609e628fe9860aed758b Mon Sep 17 00:00:00 2001 From: Logan Ward Date: Tue, 23 Aug 2022 15:29:45 -0400 Subject: [PATCH 14/29] Add ability to set aliases from API Sets them for the current series and each shock within the series --- frhodo/api/__init__.py | 27 ++++++++++++++++++++++----- tests/test_api.py | 8 +++++++- 2 files changed, 29 insertions(+), 6 deletions(-) diff --git a/frhodo/api/__init__.py b/frhodo/api/__init__.py index 2441f7d..82741e5 100644 --- a/frhodo/api/__init__.py +++ b/frhodo/api/__init__.py @@ -5,7 +5,7 @@ """ from pathlib import Path -from typing import List +from typing import List, Optional, Dict import numpy as np from PyQt5.QtWidgets import QApplication @@ -16,8 +16,13 @@ class FrhodoDriver: """Driver the Frhodo GUI application - Implements simple interfaces for common tasks, such as loading differnet datasets + Implements simple interfaces for common tasks, such as loading different datasets or evaluating changes to different mechanisms + + **Limitations** + + The driver only supports a subset of + We only provide support for experiments that use experimental data from a single series """ def __init__(self, window: Main, app: QApplication): @@ -32,13 +37,16 @@ def __init__(self, window: Main, app: QApplication): def load_files(self, experiment_path: Path, mechanism_path: Path, - output_path: Path): + output_path: Path, + aliases: Optional[Dict[str, str]] = None): """Load the input files for Frhodo Args: experiment_path: Path to the experimental data mechanism_path: Path to the mechanism files output_path: Path to which simulation results should be stored + aliases: List of aliases of chemical species as mapping of the name + of a chemical in the experiment data to name of the same chemical in the mechanism. """ self.window.exp_main_box.setPlainText(str(experiment_path.resolve().absolute())) @@ -48,6 +56,13 @@ def load_files(self, experiment_path: Path, # Trigger Frhodo to process these files self.app.processEvents() + # Add in the aliases + if aliases is not None: + aliases = aliases.copy() # Ensure that later changes to aliases don't propagate here + self.window.series.species_alias[0] = aliases + for shock in self.window.series.shock[0]: + shock['species_alias'] = aliases + @property def n_shocks(self): """Number of shock experiments loaded""" @@ -65,8 +80,10 @@ def _select_shock(self, n: int): """ # Check if it is in the list - assert self.n_shocks > 0 - assert any(n == x['num'] for x in self.window.series.shock[0]) + assert self.n_shocks > 0, 'No shocks are loaded' + assert n in self.window.series.current['shock_num'], f'Shock number {n} is not in the current series' + + # Trigger an update of the displayed system self.window.shock_choice_box.setValue(n + 1) self.app.processEvents() diff --git a/tests/test_api.py b/tests/test_api.py index cf1c1ef..794dd69 100644 --- a/tests/test_api.py +++ b/tests/test_api.py @@ -1,11 +1,13 @@ """Testing the API components of Frhodo""" + def _load_example(frhodo_driver, example_dir, tmp_path): """Set up the driver with a specific problem case""" frhodo_driver.load_files( example_dir / 'Experiment', example_dir / 'Mechanism', - tmp_path + tmp_path, + {'A': 'B'} # A fake alias ) @@ -34,3 +36,7 @@ def test_simulate(frhodo_driver, example_dir, tmp_path): assert len(runs) == 1 assert runs[0].ndim == 2 assert runs[0].shape[1] == 2 + + # Make sure the aliases propagated through + # They are reloaded from configuration variables when we re-run a simulation + assert frhodo_driver.window.display_shock['species_alias']['A'] == 'B' From a910f4b6520ba0559c3f7c4f1b9e256f6833f18e Mon Sep 17 00:00:00 2001 From: Logan Ward Date: Tue, 23 Aug 2022 18:08:37 -0400 Subject: [PATCH 15/29] Add function for modifying reaction coefficients Still need to figure out how to modify parameters from pressure-dependent reactions --- frhodo/api/__init__.py | 52 +++++++++++++++++++++++++++++++++++++++--- frhodo/mech_widget.py | 1 + frhodo/settings.py | 4 +++- tests/test_api.py | 49 ++++++++++++++++++++++++++++++--------- 4 files changed, 91 insertions(+), 15 deletions(-) diff --git a/frhodo/api/__init__.py b/frhodo/api/__init__.py index 82741e5..05be762 100644 --- a/frhodo/api/__init__.py +++ b/frhodo/api/__init__.py @@ -5,7 +5,7 @@ """ from pathlib import Path -from typing import List, Optional, Dict +from typing import List, Optional, Dict, Tuple import numpy as np from PyQt5.QtWidgets import QApplication @@ -21,8 +21,10 @@ class FrhodoDriver: **Limitations** - The driver only supports a subset of - We only provide support for experiments that use experimental data from a single series + The driver only supports a subset of features of Frhodo + + - Only supports experiments that use experimental data from a single series + - Does not support changing parameters to reactions which are pressure-dependent """ def __init__(self, window: Main, app: QApplication): @@ -72,6 +74,11 @@ def n_shocks(self): return 0 return len(self.window.series.shock[0]) + @property + def rxn_coeffs(self) -> List[List[Dict[str, float]]]: + """Reaction rate coefficients""" + return self.window.mech.coeffs + def _select_shock(self, n: int): """Change which shock experiment is being displayed and simulated @@ -126,3 +133,42 @@ def run_simulations(self) -> List[np.ndarray]: self.window.SIM.observable ], axis=1)) return output + + def get_reaction_rates(self) -> np.ndarray: + """Get the reaction rates for each shock experiment + + Returns: + Array where each row is a different reaction rate and each column is + a different shock tube experiment + """ + + # Run them directly through the `series` interface rather than changing the display + # to avoid re-running the + output = [] + for shock in self.window.series.shock[0]: + output.append(self.window.series.rates(shock)) + + # Call with the current shock to ensure `mech` hasn't been changed + self.window.series.rates(self.window.display_shock) + return np.stack(output, axis=1) + + def change_coefficient(self, new_values: Dict[Tuple[int, int, str], float]): + """Update the parameters of a reaction parameter + + Args: + new_values: A dictionary where the key is a tuple defining which + coefficient is being altered: (, , ) + The value is the new value for that coefficient + """ + + # Update the value in the Chemical_Mechanism dictionary + for key, new_value in new_values.items(): + rxn_id, prs_id, coef_name = key + rxn_model = self.window.mech.coeffs[rxn_id][prs_id] + assert coef_name in rxn_model, f'Key {coef_name} not present for reaction #{rxn_id} at pressure #{prs_id}' + rxn_model[coef_name] = new_value + + # Update the reaction rates for the current shock + self.window.mech.modify_reactions(self.window.mech.coeffs) + self.window.tree.update_rates() + self.app.processEvents() diff --git a/frhodo/mech_widget.py b/frhodo/mech_widget.py index 549a1aa..8c46819 100644 --- a/frhodo/mech_widget.py +++ b/frhodo/mech_widget.py @@ -295,6 +295,7 @@ def currentRxn(self): return None def update_value(self, event): + """Update the value from when a user interacts with the mechanism widget""" def getRateConst(parent, rxnNum, coef_key, coefName, value): shock = parent.display_shock parent.mech.coeffs[rxnNum][coef_key][coefName] = value diff --git a/frhodo/settings.py b/frhodo/settings.py index c2e9476..0b453c6 100644 --- a/frhodo/settings.py +++ b/frhodo/settings.py @@ -1008,7 +1008,9 @@ def thermo_mix(self, shock=[]): # elif species not in parent.mech.gas.species_names: return def rates(self, shock, rxnIdx=()): - """Resets and updates all rates in shock + """Resets and updates all rates in a single shock experiment + + Modifies the `mech` object of the Frhodo main window Args: shock: Parameters of the shock to update rates for diff --git a/tests/test_api.py b/tests/test_api.py index 794dd69..6929639 100644 --- a/tests/test_api.py +++ b/tests/test_api.py @@ -1,7 +1,10 @@ """Testing the API components of Frhodo""" +import numpy as np +from pytest import fixture -def _load_example(frhodo_driver, example_dir, tmp_path): +@fixture +def loaded_frhodo(frhodo_driver, example_dir, tmp_path): """Set up the driver with a specific problem case""" frhodo_driver.load_files( example_dir / 'Experiment', @@ -9,6 +12,7 @@ def _load_example(frhodo_driver, example_dir, tmp_path): tmp_path, {'A': 'B'} # A fake alias ) + return frhodo_driver def test_launch(frhodo_driver): @@ -16,27 +20,50 @@ def test_launch(frhodo_driver): assert frhodo_driver.window.isVisible() -def test_load(frhodo_driver, example_dir, tmp_path): +def test_load(loaded_frhodo): """Load loading in desired data""" - _load_example(frhodo_driver, example_dir, tmp_path) - assert frhodo_driver.n_shocks == 1 + assert loaded_frhodo.n_shocks == 1 -def test_observables(frhodo_driver, example_dir, tmp_path): - _load_example(frhodo_driver, example_dir, tmp_path) - runs = frhodo_driver.get_observables() +def test_observables(loaded_frhodo): + runs = loaded_frhodo.get_observables() assert len(runs) == 1 assert runs[0].ndim == 2 assert runs[0].shape[1] == 2 -def test_simulate(frhodo_driver, example_dir, tmp_path): - _load_example(frhodo_driver, example_dir, tmp_path) - runs = frhodo_driver.run_simulations() +def test_simulate(loaded_frhodo): + runs = loaded_frhodo.run_simulations() assert len(runs) == 1 assert runs[0].ndim == 2 assert runs[0].shape[1] == 2 # Make sure the aliases propagated through # They are reloaded from configuration variables when we re-run a simulation - assert frhodo_driver.window.display_shock['species_alias']['A'] == 'B' + assert loaded_frhodo.window.display_shock['species_alias']['A'] == 'B' + + +def test_update(loaded_frhodo): + # Get the initial reaction rates + rates = loaded_frhodo.get_reaction_rates() + assert rates.shape == (66, 1) + + # Get the simulation before + sim = loaded_frhodo.run_simulations()[0] + + # Update one of the rates + rxn_id = 3 # A simple reaction + loaded_frhodo.change_coefficient({(rxn_id, 0, 'pre_exponential_factor'): 200}) + assert loaded_frhodo.rxn_coeffs[rxn_id][0]['pre_exponential_factor'] == 200 + + # Make sure that only one rate changes + new_rates = loaded_frhodo.get_reaction_rates() + changes = np.abs(rates - new_rates) + assert changes[rxn_id, 0] > 0 + mask = np.ones_like(rates, bool) + mask[rxn_id, 0] = False + assert np.isclose(changes[mask], 0).all() + + # Make sure that this changes the simulation + new_sim = loaded_frhodo.run_simulations()[0] + assert not np.isclose(sim[-1, 1], new_sim[-1, 1], rtol=1e-4) # Look just at the last point From 3f5053eaa649668872f2854cd2953b812130c335 Mon Sep 17 00:00:00 2001 From: Logan Ward Date: Thu, 25 Aug 2022 10:32:54 -0400 Subject: [PATCH 16/29] Add tests for pressure-dependent reactions Fails for now --- tests/test_api.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/tests/test_api.py b/tests/test_api.py index 6929639..441b818 100644 --- a/tests/test_api.py +++ b/tests/test_api.py @@ -1,6 +1,6 @@ """Testing the API components of Frhodo""" import numpy as np -from pytest import fixture +from pytest import fixture, mark @fixture @@ -43,7 +43,10 @@ def test_simulate(loaded_frhodo): assert loaded_frhodo.window.display_shock['species_alias']['A'] == 'B' -def test_update(loaded_frhodo): +@mark.parametrize( + 'rxn_id', [1, 3] # TODO (wardlt): Reaction 1 fails because we don't yet support pressure-dependent reactions +) +def test_update(loaded_frhodo, rxn_id): # Get the initial reaction rates rates = loaded_frhodo.get_reaction_rates() assert rates.shape == (66, 1) @@ -52,7 +55,6 @@ def test_update(loaded_frhodo): sim = loaded_frhodo.run_simulations()[0] # Update one of the rates - rxn_id = 3 # A simple reaction loaded_frhodo.change_coefficient({(rxn_id, 0, 'pre_exponential_factor'): 200}) assert loaded_frhodo.rxn_coeffs[rxn_id][0]['pre_exponential_factor'] == 200 From bfe4b9aaffc2f6da37fac8a69b979ad4ff865b97 Mon Sep 17 00:00:00 2001 From: Logan Ward Date: Thu, 25 Aug 2022 10:38:45 -0400 Subject: [PATCH 17/29] Add shortcut for launching the driver --- frhodo/api/__init__.py | 16 +++++++++++++++- tests/test_api.py | 7 +++++-- 2 files changed, 20 insertions(+), 3 deletions(-) diff --git a/frhodo/api/__init__.py b/frhodo/api/__init__.py index 05be762..4bfe8cb 100644 --- a/frhodo/api/__init__.py +++ b/frhodo/api/__init__.py @@ -10,7 +10,7 @@ import numpy as np from PyQt5.QtWidgets import QApplication -from ..main import Main +from ..main import Main, launch_gui class FrhodoDriver: @@ -172,3 +172,17 @@ def change_coefficient(self, new_values: Dict[Tuple[int, int, str], float]): self.window.mech.modify_reactions(self.window.mech.coeffs) self.window.tree.update_rates() self.app.processEvents() + + @classmethod + def create_driver(cls, headless: bool = True, fresh: bool = True) -> 'FrhodoDriver': + """Create a Frhodo driver + + Args: + headless: Whether to suppress the display + fresh: Whether to skip loading in previous configuration + Returns: + Driver connected to a new Frhodo instance + """ + qt_args = ['frhodo', '-platform', 'offscreen'] if headless else ['frhodo'] + app, window = launch_gui(qt_args, fresh_gui=fresh) + return cls(window, app) diff --git a/tests/test_api.py b/tests/test_api.py index 441b818..e24f3b6 100644 --- a/tests/test_api.py +++ b/tests/test_api.py @@ -2,6 +2,8 @@ import numpy as np from pytest import fixture, mark +from frhodo.api import FrhodoDriver + @fixture def loaded_frhodo(frhodo_driver, example_dir, tmp_path): @@ -15,9 +17,10 @@ def loaded_frhodo(frhodo_driver, example_dir, tmp_path): return frhodo_driver -def test_launch(frhodo_driver): +def test_launch(): """Make sure we can launch Frhodo""" - assert frhodo_driver.window.isVisible() + driver = FrhodoDriver.create_driver() + assert driver.window.isVisible() def test_load(loaded_frhodo): From 8f86e484877da861b8c55c5377f7cde3c4fbd68b Mon Sep 17 00:00:00 2001 From: Logan Ward Date: Thu, 25 Aug 2022 14:59:06 -0400 Subject: [PATCH 18/29] Retrieve weights too, re-load data before extracting Also test with Falloff reactions --- frhodo/api/__init__.py | 26 +++++++++++++++++--------- frhodo/calculate/mech_fcns.py | 20 ++++++++++++++------ tests/test_api.py | 10 ++++++---- 3 files changed, 37 insertions(+), 19 deletions(-) diff --git a/frhodo/api/__init__.py b/frhodo/api/__init__.py index 4bfe8cb..7e21961 100644 --- a/frhodo/api/__init__.py +++ b/frhodo/api/__init__.py @@ -5,7 +5,7 @@ """ from pathlib import Path -from typing import List, Optional, Dict, Tuple +from typing import List, Optional, Dict, Tuple, Union import numpy as np from PyQt5.QtWidgets import QApplication @@ -51,6 +51,7 @@ def load_files(self, experiment_path: Path, of a chemical in the experiment data to name of the same chemical in the mechanism. """ + # Set the files in the text boxes self.window.exp_main_box.setPlainText(str(experiment_path.resolve().absolute())) self.window.mech_main_box.setPlainText(str(mechanism_path.resolve().absolute())) self.window.sim_main_box.setPlainText(str(output_path.resolve().absolute())) @@ -94,26 +95,33 @@ def _select_shock(self, n: int): self.window.shock_choice_box.setValue(n + 1) self.app.processEvents() - def get_observables(self) -> List[np.ndarray]: - """Get the observable data from each shock experiment + def get_observables(self) -> Tuple[List[np.ndarray], List[np.ndarray]]: + """Get the observable data and associated weights from each shock experiment Returns: - List of experimental data arrays where each is a 2D + - List of experimental data arrays where each is a 2D array with the first column is the time and second is the observable + - 1D array of the weights for each point """ if self.n_shocks == 0: - return [] + return [], [] # Loop over each shock output = [] + output_weights = [] for shock in self.window.series.shock[0]: - output.append(shock['exp_data']) - return output + # Changing the index in the GUI forces data to be loaded + self._select_shock(shock['num']) + + # Get the data from the display shock + output.append(self.window.display_shock['exp_data']) + output_weights.append(self.window.display_shock['normalized_weights']) + return output, output_weights def run_simulations(self) -> List[np.ndarray]: - """Run the simulation for each of the observed + """Run the simulation for each of the shocks Returns: List of simulated data arrays where each is a 2D @@ -152,7 +160,7 @@ def get_reaction_rates(self) -> np.ndarray: self.window.series.rates(self.window.display_shock) return np.stack(output, axis=1) - def change_coefficient(self, new_values: Dict[Tuple[int, int, str], float]): + def change_coefficient(self, new_values: Dict[Tuple[int, Union[str, int], str], float]): """Update the parameters of a reaction parameter Args: diff --git a/frhodo/calculate/mech_fcns.py b/frhodo/calculate/mech_fcns.py index d767fa4..ea46930 100644 --- a/frhodo/calculate/mech_fcns.py +++ b/frhodo/calculate/mech_fcns.py @@ -4,7 +4,7 @@ import os, io, stat, contextlib, pathlib, time from copy import deepcopy -from typing import Tuple +from typing import Tuple, Optional, Sequence, Union import cantera as ct from cantera import interrupts, cti2yaml#, ck2yaml, ctml2yaml @@ -59,7 +59,7 @@ def loader(self, path): print(e) output = {'success': False, 'message': []} - # Intialize and report any problems to log, not to console window + # Initialize and report any problems to log, not to console window stdout = io.StringIO() stderr = io.StringIO() with contextlib.redirect_stderr(stderr): @@ -338,13 +338,21 @@ def set_thermo_expression_coeffs(self): # TODO Doesn't work with NASA 9 'coeffs': np.array(S.thermo.coeffs), 'h_scaler': 1, 's_scaler': 1,} - self.thermo_coeffs.append(thermo_dict) - - def modify_reactions(self, coeffs, rxnIdxs=[]): # Only works for Arrhenius equations currently + self.thermo_coeffs.append(thermo_dict) + + def modify_reactions(self, coeffs, rxnIdxs: Optional[Union[int, float, Sequence[int]]] = ()): + """Update the reaction models in the underlying Cantera reaction models associated with this class + + Args: + coeffs: New coefficient values + rxnIdxs: IDs of reactions to be updated + """ + # TODO (wardlt): This function is only called with `self.mech` as inputs + # TODO (sikes): Only works for Arrhenius equations currently if not rxnIdxs: # if rxnNums does not exist, modify all rxnIdxs = range(len(coeffs)) else: - if isinstance(rxnIdxs, (float, int)): # if single reaction given, run that one + if isinstance(rxnIdxs, (float, int)): # if single reaction given, run that onee rxnIdxs = [rxnIdxs] for rxnIdx in rxnIdxs: diff --git a/tests/test_api.py b/tests/test_api.py index e24f3b6..5e919ea 100644 --- a/tests/test_api.py +++ b/tests/test_api.py @@ -47,9 +47,11 @@ def test_simulate(loaded_frhodo): @mark.parametrize( - 'rxn_id', [1, 3] # TODO (wardlt): Reaction 1 fails because we don't yet support pressure-dependent reactions + 'rxn_id, prs_id', [(1, 0), # PLog TODO (wardlt): Reaction 1 fails because we don't yet support PLog reactions + (3, 0), # Elementary + (36, 'low_rate')] # Three-body ) -def test_update(loaded_frhodo, rxn_id): +def test_update(loaded_frhodo, rxn_id, prs_id): # Get the initial reaction rates rates = loaded_frhodo.get_reaction_rates() assert rates.shape == (66, 1) @@ -58,8 +60,8 @@ def test_update(loaded_frhodo, rxn_id): sim = loaded_frhodo.run_simulations()[0] # Update one of the rates - loaded_frhodo.change_coefficient({(rxn_id, 0, 'pre_exponential_factor'): 200}) - assert loaded_frhodo.rxn_coeffs[rxn_id][0]['pre_exponential_factor'] == 200 + loaded_frhodo.change_coefficient({(rxn_id, prs_id, 'pre_exponential_factor'): 200}) + assert loaded_frhodo.rxn_coeffs[rxn_id][prs_id]['pre_exponential_factor'] == 200 # Make sure that only one rate changes new_rates = loaded_frhodo.get_reaction_rates() From 6fa12b4012c4e47ed613a58e8d9fada20667f602 Mon Sep 17 00:00:00 2001 From: Logan Ward Date: Thu, 25 Aug 2022 15:31:25 -0400 Subject: [PATCH 19/29] Add commands for listing and retreiving parameter values --- frhodo/api/__init__.py | 57 +++++++++++++++++++++++++++++++++++++++--- tests/test_api.py | 24 ++++++++++++++++-- 2 files changed, 76 insertions(+), 5 deletions(-) diff --git a/frhodo/api/__init__.py b/frhodo/api/__init__.py index 7e21961..b8db875 100644 --- a/frhodo/api/__init__.py +++ b/frhodo/api/__init__.py @@ -5,13 +5,16 @@ """ from pathlib import Path -from typing import List, Optional, Dict, Tuple, Union +from typing import List, Optional, Dict, Tuple, Union, Any import numpy as np +import cantera as ct from PyQt5.QtWidgets import QApplication from ..main import Main, launch_gui +CoefIndex = Tuple[int, Union[str, int], str] +"""How to specify the index of a reaction parameter: reaction index, pressure index, name of the parameter""" class FrhodoDriver: """Driver the Frhodo GUI application @@ -76,7 +79,7 @@ def n_shocks(self): return len(self.window.series.shock[0]) @property - def rxn_coeffs(self) -> List[List[Dict[str, float]]]: + def rxn_coeffs(self) -> List[Union[List[Dict[str, float]], Dict[str, Any]]]: """Reaction rate coefficients""" return self.window.mech.coeffs @@ -160,7 +163,55 @@ def get_reaction_rates(self) -> np.ndarray: self.window.series.rates(self.window.display_shock) return np.stack(output, axis=1) - def change_coefficient(self, new_values: Dict[Tuple[int, Union[str, int], str], float]): + def get_fittable_parameters(self, chosen_reactions: Optional[List[int]] = None) -> List[CoefIndex]: + """Get a list of all parameters than can be fit during optimization + + Args: + chosen_reactions: Which reactions we're allowed to optimize + Returns: + Index of each parameter that could be modified. Indices are defined by: + - Reaction number + - Index of the pressure level + - Name of the coefficient + """ + + output = [] + # Loop over all reactions, getting both the coefficients (stored in Fhrodo) + # and reaction models (stored in a Cantera reaction object) + for rxn_id, (coeffs, rxn_obj) in enumerate(zip(self.rxn_coeffs, self.window.mech.gas.reactions())): + # Determine whether to skip this reaction + if chosen_reactions is not None and rxn_id not in chosen_reactions: + continue + + # We get the high and low pressure for the falloff reactions + if type(rxn_obj) is ct.FalloffReaction: + for p_id in ['low_rate', 'high_rate']: + for coeff in coeffs[p_id].keys(): + output.append((rxn_id, p_id, coeff)) + else: + # Loop over pressure levels + for p_id, p_coeffs in enumerate(coeffs): + for coeff in p_coeffs.keys(): + output.append((rxn_id, p_id, coeff)) + return output + + def get_coefficients(self, indices: List[CoefIndex]) -> List[float]: + """Get the values of specific coefficients + + Args: + indices: List of coefficients to extract + Returns: + List of their current values + """ + + output = [] + for rxn_id, prs_id, coef_name in indices: + rxn_model = self.window.mech.coeffs[rxn_id][prs_id] + assert coef_name in rxn_model, f'Key {coef_name} not present for reaction #{rxn_id} at pressure #{prs_id}' + output.append(rxn_model[coef_name]) + return output + + def change_coefficient(self, new_values: Dict[CoefIndex, float]): """Update the parameters of a reaction parameter Args: diff --git a/tests/test_api.py b/tests/test_api.py index 5e919ea..465a95e 100644 --- a/tests/test_api.py +++ b/tests/test_api.py @@ -29,11 +29,14 @@ def test_load(loaded_frhodo): def test_observables(loaded_frhodo): - runs = loaded_frhodo.get_observables() + runs, weights = loaded_frhodo.get_observables() assert len(runs) == 1 assert runs[0].ndim == 2 assert runs[0].shape[1] == 2 + assert len(weights) == 1 + assert weights[0].size == runs[0].shape[0] + def test_simulate(loaded_frhodo): runs = loaded_frhodo.run_simulations() @@ -47,11 +50,12 @@ def test_simulate(loaded_frhodo): @mark.parametrize( - 'rxn_id, prs_id', [(1, 0), # PLog TODO (wardlt): Reaction 1 fails because we don't yet support PLog reactions + 'rxn_id, prs_id', [(0, 0), # PLog TODO (wardlt): Reaction 0 fails because we don't yet support PLog reactions (3, 0), # Elementary (36, 'low_rate')] # Three-body ) def test_update(loaded_frhodo, rxn_id, prs_id): + """Test that we can update individual parameters""" # Get the initial reaction rates rates = loaded_frhodo.get_reaction_rates() assert rates.shape == (66, 1) @@ -74,3 +78,19 @@ def test_update(loaded_frhodo, rxn_id, prs_id): # Make sure that this changes the simulation new_sim = loaded_frhodo.run_simulations()[0] assert not np.isclose(sim[-1, 1], new_sim[-1, 1], rtol=1e-4) # Look just at the last point + + +def test_fittable_parameters(loaded_frhodo): + """Test getting lists of parameters to update""" + + # Make sure we can get all parameters + total = loaded_frhodo.get_fittable_parameters() + assert len(total) == 434 + + # Make sure it works with each reaction type + assert len(loaded_frhodo.get_fittable_parameters([0])) == 32 # plog + assert len(loaded_frhodo.get_fittable_parameters([3])) == 3 # elementary + assert len(loaded_frhodo.get_fittable_parameters([36])) == 6 # Troe + + # Test getting a parameter + assert np.isclose(loaded_frhodo.get_coefficients([(1, 0, 'pre_exponential_factor')]), 5.9102033e+93) From 91c699e4bed86c4166a619fa76f5b3bee7fadf3e Mon Sep 17 00:00:00 2001 From: Logan Ward Date: Thu, 25 Aug 2022 16:46:12 -0400 Subject: [PATCH 20/29] Pressure and efficiencies are not adjustable parameters --- frhodo/api/__init__.py | 7 +++++-- tests/test_api.py | 4 ++-- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/frhodo/api/__init__.py b/frhodo/api/__init__.py index b8db875..3873392 100644 --- a/frhodo/api/__init__.py +++ b/frhodo/api/__init__.py @@ -174,6 +174,7 @@ def get_fittable_parameters(self, chosen_reactions: Optional[List[int]] = None) - Index of the pressure level - Name of the coefficient """ + _not_parameter = ['efficiencies', 'Pressure'] output = [] # Loop over all reactions, getting both the coefficients (stored in Fhrodo) @@ -187,12 +188,14 @@ def get_fittable_parameters(self, chosen_reactions: Optional[List[int]] = None) if type(rxn_obj) is ct.FalloffReaction: for p_id in ['low_rate', 'high_rate']: for coeff in coeffs[p_id].keys(): - output.append((rxn_id, p_id, coeff)) + if coeff not in _not_parameter: + output.append((rxn_id, p_id, coeff)) else: # Loop over pressure levels for p_id, p_coeffs in enumerate(coeffs): for coeff in p_coeffs.keys(): - output.append((rxn_id, p_id, coeff)) + if coeff not in _not_parameter: + output.append((rxn_id, p_id, coeff)) return output def get_coefficients(self, indices: List[CoefIndex]) -> List[float]: diff --git a/tests/test_api.py b/tests/test_api.py index 465a95e..c96009c 100644 --- a/tests/test_api.py +++ b/tests/test_api.py @@ -85,10 +85,10 @@ def test_fittable_parameters(loaded_frhodo): # Make sure we can get all parameters total = loaded_frhodo.get_fittable_parameters() - assert len(total) == 434 + assert len(total) == 360 # Make sure it works with each reaction type - assert len(loaded_frhodo.get_fittable_parameters([0])) == 32 # plog + assert len(loaded_frhodo.get_fittable_parameters([0])) == 24 # plog assert len(loaded_frhodo.get_fittable_parameters([3])) == 3 # elementary assert len(loaded_frhodo.get_fittable_parameters([36])) == 6 # Troe From 972ae9742ff4331a74cb2964fffdb3073aef20cd Mon Sep 17 00:00:00 2001 From: Logan Ward Date: Thu, 25 Aug 2022 16:50:35 -0400 Subject: [PATCH 21/29] Broke FrhodoDriver off into its own module --- frhodo/api/__init__.py | 250 +---------------------------------------- frhodo/api/driver.py | 247 ++++++++++++++++++++++++++++++++++++++++ tests/conftest.py | 2 +- tests/test_api.py | 2 +- 4 files changed, 250 insertions(+), 251 deletions(-) create mode 100644 frhodo/api/driver.py diff --git a/frhodo/api/__init__.py b/frhodo/api/__init__.py index 3873392..50785c0 100644 --- a/frhodo/api/__init__.py +++ b/frhodo/api/__init__.py @@ -1,250 +1,2 @@ -"""Functions related to using Frhodo from a Python interface +"""Functions related to using Frhodo from a Python interface""" -Many of these functions require launching the Frhodo GUI first. -using :meth:`~fhrodo.main.launch_gui` -""" - -from pathlib import Path -from typing import List, Optional, Dict, Tuple, Union, Any - -import numpy as np -import cantera as ct -from PyQt5.QtWidgets import QApplication - -from ..main import Main, launch_gui - -CoefIndex = Tuple[int, Union[str, int], str] -"""How to specify the index of a reaction parameter: reaction index, pressure index, name of the parameter""" - -class FrhodoDriver: - """Driver the Frhodo GUI application - - Implements simple interfaces for common tasks, such as loading different datasets - or evaluating changes to different mechanisms - - **Limitations** - - The driver only supports a subset of features of Frhodo - - - Only supports experiments that use experimental data from a single series - - Does not support changing parameters to reactions which are pressure-dependent - """ - - def __init__(self, window: Main, app: QApplication): - """Create the driver given connection to a running instance of Frhodo - - Args: - window: Link to the main Frhodo window - app: Link to the Qt backend - """ - self.window = window - self.app = app - - def load_files(self, experiment_path: Path, - mechanism_path: Path, - output_path: Path, - aliases: Optional[Dict[str, str]] = None): - """Load the input files for Frhodo - - Args: - experiment_path: Path to the experimental data - mechanism_path: Path to the mechanism files - output_path: Path to which simulation results should be stored - aliases: List of aliases of chemical species as mapping of the name - of a chemical in the experiment data to name of the same chemical in the mechanism. - """ - - # Set the files in the text boxes - self.window.exp_main_box.setPlainText(str(experiment_path.resolve().absolute())) - self.window.mech_main_box.setPlainText(str(mechanism_path.resolve().absolute())) - self.window.sim_main_box.setPlainText(str(output_path.resolve().absolute())) - - # Trigger Frhodo to process these files - self.app.processEvents() - - # Add in the aliases - if aliases is not None: - aliases = aliases.copy() # Ensure that later changes to aliases don't propagate here - self.window.series.species_alias[0] = aliases - for shock in self.window.series.shock[0]: - shock['species_alias'] = aliases - - @property - def n_shocks(self): - """Number of shock experiments loaded""" - n_series = len(self.window.series.shock) - assert n_series <= 1, "We only support one series" - if n_series == 0: - return 0 - return len(self.window.series.shock[0]) - - @property - def rxn_coeffs(self) -> List[Union[List[Dict[str, float]], Dict[str, Any]]]: - """Reaction rate coefficients""" - return self.window.mech.coeffs - - def _select_shock(self, n: int): - """Change which shock experiment is being displayed and simulated - - Args: - n: Which shock experiment to evaluate - """ - - # Check if it is in the list - assert self.n_shocks > 0, 'No shocks are loaded' - assert n in self.window.series.current['shock_num'], f'Shock number {n} is not in the current series' - - # Trigger an update of the displayed system - self.window.shock_choice_box.setValue(n + 1) - self.app.processEvents() - - def get_observables(self) -> Tuple[List[np.ndarray], List[np.ndarray]]: - """Get the observable data and associated weights from each shock experiment - - Returns: - - List of experimental data arrays where each is a 2D - array with the first column is the time and second - is the observable - - 1D array of the weights for each point - """ - - if self.n_shocks == 0: - return [], [] - - # Loop over each shock - output = [] - output_weights = [] - for shock in self.window.series.shock[0]: - # Changing the index in the GUI forces data to be loaded - self._select_shock(shock['num']) - - # Get the data from the display shock - output.append(self.window.display_shock['exp_data']) - output_weights.append(self.window.display_shock['normalized_weights']) - return output, output_weights - - def run_simulations(self) -> List[np.ndarray]: - """Run the simulation for each of the shocks - - Returns: - List of simulated data arrays where each is a 2D - array with the first column is the time and second - is the simulated observable - """ - - # Loop over all shocks - output = [] - for shock in self.window.series.shock[0]: - # We force the simulation by changing the shock index in the GUI - # That could be an issue if the behavior of the GUI changes - self._select_shock(shock['num']) - assert self.window.SIM.success, 'Simulation failed' - output.append(np.stack([ - self.window.SIM.independent_var, - self.window.SIM.observable - ], axis=1)) - return output - - def get_reaction_rates(self) -> np.ndarray: - """Get the reaction rates for each shock experiment - - Returns: - Array where each row is a different reaction rate and each column is - a different shock tube experiment - """ - - # Run them directly through the `series` interface rather than changing the display - # to avoid re-running the - output = [] - for shock in self.window.series.shock[0]: - output.append(self.window.series.rates(shock)) - - # Call with the current shock to ensure `mech` hasn't been changed - self.window.series.rates(self.window.display_shock) - return np.stack(output, axis=1) - - def get_fittable_parameters(self, chosen_reactions: Optional[List[int]] = None) -> List[CoefIndex]: - """Get a list of all parameters than can be fit during optimization - - Args: - chosen_reactions: Which reactions we're allowed to optimize - Returns: - Index of each parameter that could be modified. Indices are defined by: - - Reaction number - - Index of the pressure level - - Name of the coefficient - """ - _not_parameter = ['efficiencies', 'Pressure'] - - output = [] - # Loop over all reactions, getting both the coefficients (stored in Fhrodo) - # and reaction models (stored in a Cantera reaction object) - for rxn_id, (coeffs, rxn_obj) in enumerate(zip(self.rxn_coeffs, self.window.mech.gas.reactions())): - # Determine whether to skip this reaction - if chosen_reactions is not None and rxn_id not in chosen_reactions: - continue - - # We get the high and low pressure for the falloff reactions - if type(rxn_obj) is ct.FalloffReaction: - for p_id in ['low_rate', 'high_rate']: - for coeff in coeffs[p_id].keys(): - if coeff not in _not_parameter: - output.append((rxn_id, p_id, coeff)) - else: - # Loop over pressure levels - for p_id, p_coeffs in enumerate(coeffs): - for coeff in p_coeffs.keys(): - if coeff not in _not_parameter: - output.append((rxn_id, p_id, coeff)) - return output - - def get_coefficients(self, indices: List[CoefIndex]) -> List[float]: - """Get the values of specific coefficients - - Args: - indices: List of coefficients to extract - Returns: - List of their current values - """ - - output = [] - for rxn_id, prs_id, coef_name in indices: - rxn_model = self.window.mech.coeffs[rxn_id][prs_id] - assert coef_name in rxn_model, f'Key {coef_name} not present for reaction #{rxn_id} at pressure #{prs_id}' - output.append(rxn_model[coef_name]) - return output - - def change_coefficient(self, new_values: Dict[CoefIndex, float]): - """Update the parameters of a reaction parameter - - Args: - new_values: A dictionary where the key is a tuple defining which - coefficient is being altered: (, , ) - The value is the new value for that coefficient - """ - - # Update the value in the Chemical_Mechanism dictionary - for key, new_value in new_values.items(): - rxn_id, prs_id, coef_name = key - rxn_model = self.window.mech.coeffs[rxn_id][prs_id] - assert coef_name in rxn_model, f'Key {coef_name} not present for reaction #{rxn_id} at pressure #{prs_id}' - rxn_model[coef_name] = new_value - - # Update the reaction rates for the current shock - self.window.mech.modify_reactions(self.window.mech.coeffs) - self.window.tree.update_rates() - self.app.processEvents() - - @classmethod - def create_driver(cls, headless: bool = True, fresh: bool = True) -> 'FrhodoDriver': - """Create a Frhodo driver - - Args: - headless: Whether to suppress the display - fresh: Whether to skip loading in previous configuration - Returns: - Driver connected to a new Frhodo instance - """ - qt_args = ['frhodo', '-platform', 'offscreen'] if headless else ['frhodo'] - app, window = launch_gui(qt_args, fresh_gui=fresh) - return cls(window, app) diff --git a/frhodo/api/driver.py b/frhodo/api/driver.py new file mode 100644 index 0000000..d2985a1 --- /dev/null +++ b/frhodo/api/driver.py @@ -0,0 +1,247 @@ +"""Low-level driver for controlling a Frhodo GUI session""" + +from pathlib import Path +from typing import Optional, Dict, List, Union, Any, Tuple + +import cantera as ct +import numpy as np +from PyQt5.QtWidgets import QApplication + +from frhodo.main import Main, launch_gui + +CoefIndex = Tuple[int, Union[str, int], str] +"""How to specify the index of a reaction parameter: reaction index, pressure index, name of the parameter""" + + +class FrhodoDriver: + """Driver the Frhodo GUI application + + Implements simple interfaces for common tasks, such as loading different datasets + or evaluating changes to different mechanisms + + **Limitations** + + The driver only supports a subset of features of Frhodo + + - Only supports experiments that use experimental data from a single series + - Does not support changing parameters to reactions which are pressure-dependent + """ + + def __init__(self, window: Main, app: QApplication): + """Create the driver given connection to a running instance of Frhodo + + Args: + window: Link to the main Frhodo window + app: Link to the Qt backend + """ + self.window = window + self.app = app + + def load_files(self, experiment_path: Path, + mechanism_path: Path, + output_path: Path, + aliases: Optional[Dict[str, str]] = None): + """Load the input files for Frhodo + + Args: + experiment_path: Path to the experimental data + mechanism_path: Path to the mechanism files + output_path: Path to which simulation results should be stored + aliases: List of aliases of chemical species as mapping of the name + of a chemical in the experiment data to name of the same chemical in the mechanism. + """ + + # Set the files in the text boxes + self.window.exp_main_box.setPlainText(str(experiment_path.resolve().absolute())) + self.window.mech_main_box.setPlainText(str(mechanism_path.resolve().absolute())) + self.window.sim_main_box.setPlainText(str(output_path.resolve().absolute())) + + # Trigger Frhodo to process these files + self.app.processEvents() + + # Add in the aliases + if aliases is not None: + aliases = aliases.copy() # Ensure that later changes to aliases don't propagate here + self.window.series.species_alias[0] = aliases + for shock in self.window.series.shock[0]: + shock['species_alias'] = aliases + + @property + def n_shocks(self): + """Number of shock experiments loaded""" + n_series = len(self.window.series.shock) + assert n_series <= 1, "We only support one series" + if n_series == 0: + return 0 + return len(self.window.series.shock[0]) + + @property + def rxn_coeffs(self) -> List[Union[List[Dict[str, float]], Dict[str, Any]]]: + """Reaction rate coefficients""" + return self.window.mech.coeffs + + def _select_shock(self, n: int): + """Change which shock experiment is being displayed and simulated + + Args: + n: Which shock experiment to evaluate + """ + + # Check if it is in the list + assert self.n_shocks > 0, 'No shocks are loaded' + assert n in self.window.series.current['shock_num'], f'Shock number {n} is not in the current series' + + # Trigger an update of the displayed system + self.window.shock_choice_box.setValue(n + 1) + self.app.processEvents() + + def get_observables(self) -> Tuple[List[np.ndarray], List[np.ndarray]]: + """Get the observable data and associated weights from each shock experiment + + Returns: + - List of experimental data arrays where each is a 2D + array with the first column is the time and second + is the observable + - 1D array of the weights for each point + """ + + if self.n_shocks == 0: + return [], [] + + # Loop over each shock + output = [] + output_weights = [] + for shock in self.window.series.shock[0]: + # Changing the index in the GUI forces data to be loaded + self._select_shock(shock['num']) + + # Get the data from the display shock + output.append(self.window.display_shock['exp_data']) + output_weights.append(self.window.display_shock['normalized_weights']) + return output, output_weights + + def run_simulations(self) -> List[np.ndarray]: + """Run the simulation for each of the shocks + + Returns: + List of simulated data arrays where each is a 2D + array with the first column is the time and second + is the simulated observable + """ + + # Loop over all shocks + output = [] + for shock in self.window.series.shock[0]: + # We force the simulation by changing the shock index in the GUI + # That could be an issue if the behavior of the GUI changes + self._select_shock(shock['num']) + assert self.window.SIM.success, 'Simulation failed' + output.append(np.stack([ + self.window.SIM.independent_var, + self.window.SIM.observable + ], axis=1)) + return output + + def get_reaction_rates(self) -> np.ndarray: + """Get the reaction rates for each shock experiment + + Returns: + Array where each row is a different reaction rate and each column is + a different shock tube experiment + """ + + # Run them directly through the `series` interface rather than changing the display + # to avoid re-running the + output = [] + for shock in self.window.series.shock[0]: + output.append(self.window.series.rates(shock)) + + # Call with the current shock to ensure `mech` hasn't been changed + self.window.series.rates(self.window.display_shock) + return np.stack(output, axis=1) + + def get_fittable_parameters(self, chosen_reactions: Optional[List[int]] = None) -> List[CoefIndex]: + """Get a list of all parameters than can be fit during optimization + + Args: + chosen_reactions: Which reactions we're allowed to optimize + Returns: + Index of each parameter that could be modified. Indices are defined by: + - Reaction number + - Index of the pressure level + - Name of the coefficient + """ + _not_parameter = ['efficiencies', 'Pressure'] + + output = [] + # Loop over all reactions, getting both the coefficients (stored in Fhrodo) + # and reaction models (stored in a Cantera reaction object) + for rxn_id, (coeffs, rxn_obj) in enumerate(zip(self.rxn_coeffs, self.window.mech.gas.reactions())): + # Determine whether to skip this reaction + if chosen_reactions is not None and rxn_id not in chosen_reactions: + continue + + # We get the high and low pressure for the falloff reactions + if type(rxn_obj) is ct.FalloffReaction: + for p_id in ['low_rate', 'high_rate']: + for coeff in coeffs[p_id].keys(): + if coeff not in _not_parameter: + output.append((rxn_id, p_id, coeff)) + else: + # Loop over pressure levels + for p_id, p_coeffs in enumerate(coeffs): + for coeff in p_coeffs.keys(): + if coeff not in _not_parameter: + output.append((rxn_id, p_id, coeff)) + return output + + def get_coefficients(self, indices: List[CoefIndex]) -> List[float]: + """Get the values of specific coefficients + + Args: + indices: List of coefficients to extract + Returns: + List of their current values + """ + + output = [] + for rxn_id, prs_id, coef_name in indices: + rxn_model = self.window.mech.coeffs[rxn_id][prs_id] + assert coef_name in rxn_model, f'Key {coef_name} not present for reaction #{rxn_id} at pressure #{prs_id}' + output.append(rxn_model[coef_name]) + return output + + def change_coefficient(self, new_values: Dict[CoefIndex, float]): + """Update the parameters of a reaction parameter + + Args: + new_values: A dictionary where the key is a tuple defining which + coefficient is being altered: (, , ) + The value is the new value for that coefficient + """ + + # Update the value in the Chemical_Mechanism dictionary + for key, new_value in new_values.items(): + rxn_id, prs_id, coef_name = key + rxn_model = self.window.mech.coeffs[rxn_id][prs_id] + assert coef_name in rxn_model, f'Key {coef_name} not present for reaction #{rxn_id} at pressure #{prs_id}' + rxn_model[coef_name] = new_value + + # Update the reaction rates for the current shock + self.window.mech.modify_reactions(self.window.mech.coeffs) + self.window.tree.update_rates() + self.app.processEvents() + + @classmethod + def create_driver(cls, headless: bool = True, fresh: bool = True) -> 'FrhodoDriver': + """Create a Frhodo driver + + Args: + headless: Whether to suppress the display + fresh: Whether to skip loading in previous configuration + Returns: + Driver connected to a new Frhodo instance + """ + qt_args = ['frhodo', '-platform', 'offscreen'] if headless else ['frhodo'] + app, window = launch_gui(qt_args, fresh_gui=fresh) + return cls(window, app) diff --git a/tests/conftest.py b/tests/conftest.py index d4c92a8..efb6070 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -2,7 +2,7 @@ from pytest import fixture -from frhodo.api import FrhodoDriver +from frhodo.api.driver import FrhodoDriver from frhodo.main import launch_gui # Launch the application headless and without reading the configuration from a past run diff --git a/tests/test_api.py b/tests/test_api.py index c96009c..d3dc43c 100644 --- a/tests/test_api.py +++ b/tests/test_api.py @@ -2,7 +2,7 @@ import numpy as np from pytest import fixture, mark -from frhodo.api import FrhodoDriver +from frhodo.api.driver import FrhodoDriver @fixture From 0abcc5c09196c8b0ec9fa45ab883e36706751e96 Mon Sep 17 00:00:00 2001 From: Logan Ward Date: Fri, 26 Aug 2022 12:38:34 -0400 Subject: [PATCH 22/29] Introduce multiprocessing-ready optimization function --- frhodo/api/driver.py | 3 +- frhodo/api/optimize.py | 201 +++++++++++++++++++++++++++++++++++++++++ tests/test_api.py | 57 +++++++++++- 3 files changed, 258 insertions(+), 3 deletions(-) create mode 100644 frhodo/api/optimize.py diff --git a/frhodo/api/driver.py b/frhodo/api/driver.py index d2985a1..6361d34 100644 --- a/frhodo/api/driver.py +++ b/frhodo/api/driver.py @@ -23,6 +23,7 @@ class FrhodoDriver: The driver only supports a subset of features of Frhodo + - We can only run one driver per process. You will even need to "spawn" new processes rather than fork on Linux - Only supports experiments that use experimental data from a single series - Does not support changing parameters to reactions which are pressure-dependent """ @@ -211,7 +212,7 @@ def get_coefficients(self, indices: List[CoefIndex]) -> List[float]: output.append(rxn_model[coef_name]) return output - def change_coefficient(self, new_values: Dict[CoefIndex, float]): + def set_coefficients(self, new_values: Dict[CoefIndex, float]): """Update the parameters of a reaction parameter Args: diff --git a/frhodo/api/optimize.py b/frhodo/api/optimize.py new file mode 100644 index 0000000..0ad205e --- /dev/null +++ b/frhodo/api/optimize.py @@ -0,0 +1,201 @@ +"""Utilities useful for using Frhodo from external optimizers""" +import multiprocessing +import warnings +from pathlib import Path +from typing import List, Optional, Dict, Sequence + +import numpy as np +from scipy import stats as ss +from scipy.interpolate import interp1d + +from frhodo.api.driver import CoefIndex, FrhodoDriver + +_frhodo: Optional[FrhodoDriver] = None + + +class BaseObjectiveFunction: + """Base class for a Frhodo-based objective function + + Launches its own Frhodo GUI instance the first time you invoke a prediction function. + """ + + def __init__( + self, + exp_directory: Path, + mech_directory: Path, + parameters: List[CoefIndex], + aliases: Optional[Dict[str, str]] = None, + ): + """ + + Args: + exp_directory: Path to the experimental data + mech_directory: Path to the mechanism file(s) + parameters: Parameters to be adjusted + aliases: Aliases that map species in the experiment to the mechanism description + """ + self.parameters = parameters.copy() + + # Make a placeholder for the Frhodo instance + self._exp_directory = exp_directory + self._mech_directory = mech_directory + self._aliases = aliases + self._frhodo: Optional[FrhodoDriver] = None + + # Make a placeholder for experimental data + self._observations: Optional[List[np.ndarray]] = None + self._weights: Optional[List[np.ndarray]] = None + + def __getstate__(self): + if multiprocessing.get_start_method() != "spawn": + warnings.warn('You must set the multiprocessing start method to spawn so we can run >1 multiprocessing instance') + + @property + def frhodo(self) -> FrhodoDriver: + """Link to the underlying Frhodo instance""" + global _frhodo # Use the module-level instance of Frhodo + # Launch if needed + if _frhodo is None: + _frhodo = FrhodoDriver.create_driver() + + # Load in the data and mechanism + _frhodo.load_files( + self._exp_directory, + self._mech_directory, + self._mech_directory / 'outputs', + aliases=self._aliases + ) + return _frhodo + + @property + def observations(self) -> List[np.ndarray]: + """Experimental data less any regions that are masked out""" + # Load the data in if needed + if self._observations is None: + self._load_experiments() + return self._observations + + @property + def weights(self) -> List[np.ndarray]: + """Weights for each observation""" + if self._observations is None: + self._load_experiments() + return self._weights + + @property + def x(self) -> np.ndarray: + """Get the current state of the coefficient being optimized + + This is either the initial values from loading the mechanism or the last values ran. + """ + return np.array(self.frhodo.get_coefficients(self.parameters)) + + def _load_experiments(self): + """Load observations and weights from disk + + Sets the values in `self._observations` and `self._weights` + """ + self._observations = [] + self._weights = [] + for obs, weights in zip(*self.frhodo.get_observables()): + # Find portions that are weighted sufficiently + mask = weights > 1e-6 + + # Store the masked weights and observations + self._observations.append(obs[mask, :]) + self._weights.append(weights[mask]) + + def run_simulations(self, x: Sequence[float]) -> List[np.ndarray]: + """Run the simulations with a new set of parameters + + Args: + x: New set of reaction coefficients + Returns: + Simulated for each experiment interpolated at the same time increments + as :meth:`observations`. + """ + + # Update the parameters + assert len(x) == len(self.parameters), f"Expected {len(self.parameters)} parameters but got {len(x)}" + self.frhodo.set_coefficients(dict(zip(self.parameters, x))) + + # Run each simulation to each set of experimental data + sims = self.frhodo.run_simulations() + output = [] + for sim, obs in zip(sims, self.observations): + # Interpolate simulation data over the same steps as the experiments + sim_func = interp1d(sim[:, 0], sim[:, 1], kind='cubic') + output.append(sim_func(obs[:, 0])) + return output + + def compute_residuals(self, x: Sequence[float]) -> List[np.ndarray]: + """Compute the residual between simulation and experiment + + Args: + x: New set of reaction coefficients + Returns: + Residuals for each shock experiment at each point + """ + + # Run the simulations with the new parameters + sims = self.run_simulations(x) + + # Compute residuals + output = [] + for sim, obs in zip(sims, self.observations): + output.append(np.subtract(sim, obs[:, 1])) + return output + + def __call__(self, x: np.ndarray, **kwargs): + """Invoke the objective function""" + raise NotImplementedError() + + +class BayesianObjectiveFunction(BaseObjectiveFunction): + """Computes the log-probability of observing experimental data given the simulated results + + Users must also pass an estimated "uncertainty" of experimental measurements along with the + other input parameters to :meth:`__call__`. We place it as the first argument in the list. + + Uses a t-Distribution error model for the data to be robust against noise in the data and + weighs data differently by scaling the width of the t-distribution so that it is wider + for data with smaller weights. + These approaches are modelled after the techniques demonstrated by + `Paulson et al. 2019 `_. + """ + + def compute_log_probs(self, x: np.ndarray, uncertainty: float) -> np.ndarray: + """Compute the log probability of observing each shock experiment + + Args: + x: New values for each coefficient + uncertainty: Size of the uncertainty for the observables + """ + + resids = self.compute_residuals(x) + output = [] + for resid, weights in zip(resids, self.weights): + # Adjust uncertainty so that the error tolerance of less-important data is larger + # See doi:10.1016/j.ijengsci.2019.05.011 + std = uncertainty / weights + output.append(ss.t(loc=0, scale=std, df=2.1).logpdf(resid).sum()) + return np.array(output) + + def _load_experiments(self): + super()._load_experiments() + + # Adjust the weights such that the largest is 1 + for i, weight in enumerate(self.weights): + self.weights[i] /= weight.max() + + def __call__(self, x: np.ndarray, **kwargs): + """Invoke the objective function + + Args: + x: Uncertainty of the observations followed by the new coefficients + """ + assert len(kwargs) == 0, "This function does not take keyword arguments" + + uncertainty, coeffs = x[0], x[1:] + logprobs = self.compute_log_probs(coeffs, uncertainty) + return logprobs.sum() diff --git a/tests/test_api.py b/tests/test_api.py index d3dc43c..682992e 100644 --- a/tests/test_api.py +++ b/tests/test_api.py @@ -1,8 +1,11 @@ """Testing the API components of Frhodo""" import numpy as np from pytest import fixture, mark +from multiprocessing import Pool, set_start_method from frhodo.api.driver import FrhodoDriver +from frhodo.api.optimize import BayesianObjectiveFunction +from frhodo.api import optimize @fixture @@ -63,9 +66,14 @@ def test_update(loaded_frhodo, rxn_id, prs_id): # Get the simulation before sim = loaded_frhodo.run_simulations()[0] + # Get the original value + coef_ind = (rxn_id, prs_id, 'pre_exponential_factor') + orig_value = loaded_frhodo.get_coefficients([coef_ind])[0] + # Update one of the rates - loaded_frhodo.change_coefficient({(rxn_id, prs_id, 'pre_exponential_factor'): 200}) + loaded_frhodo.set_coefficients({coef_ind: 200}) assert loaded_frhodo.rxn_coeffs[rxn_id][prs_id]['pre_exponential_factor'] == 200 + assert loaded_frhodo.get_coefficients([coef_ind]) == [200] # Make sure that only one rate changes new_rates = loaded_frhodo.get_reaction_rates() @@ -77,7 +85,11 @@ def test_update(loaded_frhodo, rxn_id, prs_id): # Make sure that this changes the simulation new_sim = loaded_frhodo.run_simulations()[0] - assert not np.isclose(sim[-1, 1], new_sim[-1, 1], rtol=1e-4) # Look just at the last point + try: + assert not np.isclose(sim[-1, 1], new_sim[-1, 1], rtol=1e-4) # Look just at the last point + finally: + # Set it back to not mess-up our other tasks + loaded_frhodo.set_coefficients({coef_ind: orig_value}) def test_fittable_parameters(loaded_frhodo): @@ -94,3 +106,44 @@ def test_fittable_parameters(loaded_frhodo): # Test getting a parameter assert np.isclose(loaded_frhodo.get_coefficients([(1, 0, 'pre_exponential_factor')]), 5.9102033e+93) + + +def test_optimizer(loaded_frhodo, example_dir, tmp_path): + set_start_method("spawn") # Allows us to run >1 Frhodo instance + opt = BayesianObjectiveFunction( + exp_directory=example_dir / 'Experiment', + mech_directory=example_dir / 'Mechanism', + parameters=[(3, 0, 'pre_exponential_factor')] + ) + + # Set the frhodo executable for that module (we can have only 1 per process) + optimize._frhodo = loaded_frhodo + + # Test the state + assert opt.weights[0].max() == 1 + assert len(opt.x) == 1 + + # Make sure the optimizer produces different results with different inputs + x0 = opt.x.tolist() + x0.insert(0, 1e-4) + y0 = opt(x0) + + x1 = list(x0) + x1[1] = x0[1] * 100 + y1 = opt(x1) + assert y0 != y1 + + # Re-running the initial guess should produce the same result + y0_repeat = opt(x0) + assert y0 == y0_repeat + + # Make sure the serialization works + with Pool(2) as p: + # Make sure we still get repeatability + y = p.map(opt, [x0] * 32) + assert np.isclose(y, y0).all() + + # Run with different inputs and make sure they match up + y = p.map(opt, [x0, x1]) + assert np.isclose(y, [y0, y1]).all() + From 3cf0ff10ff5794839aca211b7478f4e8750797a4 Mon Sep 17 00:00:00 2001 From: Logan Ward Date: Fri, 23 Sep 2022 11:26:14 -0400 Subject: [PATCH 23/29] Add function for running simulator without invoking GUI --- frhodo/api/driver.py | 44 +++++++++++++++++++++ frhodo/api/optimize.py | 72 +++++++++++++++++++++++++++++------ frhodo/calculate/mech_fcns.py | 3 +- frhodo/main.py | 28 ++++++++------ tests/test_api.py | 14 ++----- 5 files changed, 128 insertions(+), 33 deletions(-) diff --git a/frhodo/api/driver.py b/frhodo/api/driver.py index 6361d34..62f4419 100644 --- a/frhodo/api/driver.py +++ b/frhodo/api/driver.py @@ -7,6 +7,7 @@ import numpy as np from PyQt5.QtWidgets import QApplication +from frhodo.calculate.reactors import Simulation_Result from frhodo.main import Main, launch_gui CoefIndex = Tuple[int, Union[str, int], str] @@ -96,6 +97,23 @@ def _select_shock(self, n: int): self.window.shock_choice_box.setValue(n + 1) self.app.processEvents() + def get_simulator_inputs(self) -> Tuple[List[dict], List[Tuple[float, float, Dict[str, float]]]]: + """Get the list of variables that serve as input for the simulator for each shock that runs + + Returns: + - Keyword arguments passed to the simulator + - Temperature, pressure and concentration of different guasses + """ + var_outputs = [] + rxn_outputs = [] + + for shock in self.window.series.shock[0]: + self._select_shock(shock['num']) + var_outputs.append(self.window.get_simulation_kwargs(shock, np.array([0]))) + rxn_outputs.append((shock['T_reactor'], shock['P_reactor'], shock['thermo_mix'])) + + return var_outputs, rxn_outputs + def get_observables(self) -> Tuple[List[np.ndarray], List[np.ndarray]]: """Get the observable data and associated weights from each shock experiment @@ -143,6 +161,32 @@ def run_simulations(self) -> List[np.ndarray]: ], axis=1)) return output + def run_simulation_from_kwargs(self, sim_kwargs, rxn_conditions) -> np.ndarray: + """Perform a simulation for an arbitrary reactor and reaction conditions. + + Args: + sim_kwargs: Keywords defining the reactor + rxn_conditions: Conditions at which the reaction occurs + + Returns: + Simulated data arrays where each is a 2D + array with the first column is the time and second + is the simulated observable + """ + + # Update the mechanism for the specified reaction conditions + self.window.mech.set_TPX(*rxn_conditions) + + # Run the simulation + # TODO (wardlt): Do not hardcode the runtime or reactor conditions + sim, _ = self.window.mech.run('Incident Shock Reactor', 1.2e-05, *rxn_conditions, **sim_kwargs) + + assert sim.success, "Simulation failed" + return np.stack([ + sim.independent_var, + sim.observable + ], axis=1) + def get_reaction_rates(self) -> np.ndarray: """Get the reaction rates for each shock experiment diff --git a/frhodo/api/optimize.py b/frhodo/api/optimize.py index 0ad205e..7fe7a8e 100644 --- a/frhodo/api/optimize.py +++ b/frhodo/api/optimize.py @@ -2,7 +2,7 @@ import multiprocessing import warnings from pathlib import Path -from typing import List, Optional, Dict, Sequence +from typing import List, Optional, Dict, Sequence, Tuple import numpy as np from scipy import stats as ss @@ -45,10 +45,28 @@ def __init__( # Make a placeholder for experimental data self._observations: Optional[List[np.ndarray]] = None self._weights: Optional[List[np.ndarray]] = None + self._sim_kwargs: Optional[List[dict]] = None + self._rxn_conditions: Optional[List[Tuple[float, float, dict]]] = None + + # Track whether the data has been loaded or not + self._loaded = False def __getstate__(self): if multiprocessing.get_start_method() != "spawn": warnings.warn('You must set the multiprocessing start method to spawn so we can run >1 multiprocessing instance') + state = self.__dict__.copy() + return state + + def set_driver(self, driver: FrhodoDriver): + """Set the driver to be used for this class + + Use this if you have already launched Frhodo. + + Args: + driver: Driver to use in place of creating a new one + """ + global _frhodo + _frhodo = driver @property def frhodo(self) -> FrhodoDriver: @@ -59,12 +77,14 @@ def frhodo(self) -> FrhodoDriver: _frhodo = FrhodoDriver.create_driver() # Load in the data and mechanism - _frhodo.load_files( - self._exp_directory, - self._mech_directory, - self._mech_directory / 'outputs', - aliases=self._aliases - ) + if not self._loaded: + _frhodo.load_files( + self._exp_directory, + self._mech_directory, + self._mech_directory / 'outputs', + aliases=self._aliases + ) + assert self.observations is not None return _frhodo @property @@ -105,6 +125,8 @@ def _load_experiments(self): self._observations.append(obs[mask, :]) self._weights.append(weights[mask]) + self._sim_kwargs, self._rxn_conditions = self.frhodo.get_simulator_inputs() + def run_simulations(self, x: Sequence[float]) -> List[np.ndarray]: """Run the simulations with a new set of parameters @@ -120,10 +142,13 @@ def run_simulations(self, x: Sequence[float]) -> List[np.ndarray]: self.frhodo.set_coefficients(dict(zip(self.parameters, x))) # Run each simulation to each set of experimental data - sims = self.frhodo.run_simulations() + sims = [] + for sim_kwargs, rxn_cond in zip(self._sim_kwargs, self._rxn_conditions): + sims.append(self.frhodo.run_simulation_from_kwargs(sim_kwargs, rxn_cond)) + + # Interpolate simulation data over the same steps as the experiments output = [] for sim, obs in zip(sims, self.observations): - # Interpolate simulation data over the same steps as the experiments sim_func = interp1d(sim[:, 0], sim[:, 1], kind='cubic') output.append(sim_func(obs[:, 0])) return output @@ -164,6 +189,25 @@ class BayesianObjectiveFunction(BaseObjectiveFunction): `Paulson et al. 2019 `_. """ + def __init__( + self, + exp_directory: Path, + mech_directory: Path, + parameters: List[CoefIndex], + priors: Optional[List[ss.rv_continuous]] = None, + aliases: Optional[Dict[str, str]] = None, + ): + """ + + Args: + exp_directory: Path to the experimental data + mech_directory: Path to the mechanism file(s) + parameters: Parameters to be adjusted + aliases: Aliases that map species in the experiment to the mechanism description + """ + super().__init__(exp_directory, mech_directory, parameters, aliases) + self.priors = priors + def compute_log_probs(self, x: np.ndarray, uncertainty: float) -> np.ndarray: """Compute the log probability of observing each shock experiment @@ -196,6 +240,12 @@ def __call__(self, x: np.ndarray, **kwargs): """ assert len(kwargs) == 0, "This function does not take keyword arguments" + # Compute the log probability of the data uncertainty, coeffs = x[0], x[1:] - logprobs = self.compute_log_probs(coeffs, uncertainty) - return logprobs.sum() + logprob = self.compute_log_probs(coeffs, uncertainty).sum() + + # Compute the log probability from the priors + if self.priors is not None: + logprob += sum(p.logpdf(c) for p, c in zip(self.priors, x)) + + return logprob diff --git a/frhodo/calculate/mech_fcns.py b/frhodo/calculate/mech_fcns.py index ea46930..20436b7 100644 --- a/frhodo/calculate/mech_fcns.py +++ b/frhodo/calculate/mech_fcns.py @@ -12,6 +12,7 @@ from timeit import default_timer as timer from . import reactors, shock_fcns, integrate +from .reactors import Simulation_Result from .. import ck2yaml @@ -539,7 +540,7 @@ def get_M(rxn): return M - def run(self, reactor_choice, t_end, T_reac, P_reac, mix, **kwargs) -> Tuple[dict, dict]: + def run(self, reactor_choice, t_end, T_reac, P_reac, mix, **kwargs) -> Tuple[Simulation_Result, dict]: """Perform a simulation of the reaction output Args: diff --git a/frhodo/main.py b/frhodo/main.py index b0b4c3e..3e0030b 100644 --- a/frhodo/main.py +++ b/frhodo/main.py @@ -279,24 +279,17 @@ def run_single(self, event=None, t_save=None, rxn_changed=False): # Get the conditions of the current reactor T_reac, P_reac, mix = shock['T_reactor'], shock['P_reactor'], shock['thermo_mix'] - # Make sure the rate constants are update-to-date with the conditions on the table and specified shock - self.tree.update_rates() - # calculate all properties or observable by sending t_save tabIdx = self.plot_tab_widget.currentIndex() tabText = self.plot_tab_widget.tabText(tabIdx) if tabText == 'Sim Explorer': t_save = np.array([0]) - # Formulate the output arguments for the - SIM_kwargs = {'u_reac': shock['u2'], 'rho1': shock['rho1'], 'observable': self.display_shock['observable'], - 't_lab_save': t_save, 'sim_int_f': self.var['reactor']['sim_interp_factor'], - 'ODE_solver': self.var['reactor']['ode_solver'], - 'rtol': self.var['reactor']['ode_rtol'], 'atol': self.var['reactor']['ode_atol']} + # Make sure the rate constants are update-to-date with the conditions on the table and specified shock + self.tree.update_rates() - if '0d Reactor' in self.var['reactor']['name']: - SIM_kwargs['solve_energy'] = self.var['reactor']['solve_energy'] - SIM_kwargs['frozen_comp'] = self.var['reactor']['frozen_comp'] + # Get the keyword arguments to being passed to the simulator + SIM_kwargs = self.get_simulation_kwargs(shock, t_save) self.SIM, verbose = self.mech.run(self.var['reactor']['name'], self.var['reactor']['t_end'], T_reac, P_reac, mix, **SIM_kwargs) @@ -318,6 +311,19 @@ def run_single(self, event=None, t_save=None, rxn_changed=False): self.sim_explorer.update_plot(None) return # If mech error exit function + def get_simulation_kwargs(self, shock, t_save) -> dict: + """Get the keyword arguments passed to a simulator for a certain shock""" + + # Formulate the output arguments for the + SIM_kwargs = {'u_reac': shock['u2'], 'rho1': shock['rho1'], 'observable': self.display_shock['observable'], + 't_lab_save': t_save, 'sim_int_f': self.var['reactor']['sim_interp_factor'], + 'ODE_solver': self.var['reactor']['ode_solver'], + 'rtol': self.var['reactor']['ode_rtol'], 'atol': self.var['reactor']['ode_atol']} + if '0d Reactor' in self.var['reactor']['name']: + SIM_kwargs['solve_energy'] = self.var['reactor']['solve_energy'] + SIM_kwargs['frozen_comp'] = self.var['reactor']['frozen_comp'] + return SIM_kwargs + # def raise_error(self): # assert False diff --git a/tests/test_api.py b/tests/test_api.py index 682992e..e5c077c 100644 --- a/tests/test_api.py +++ b/tests/test_api.py @@ -51,6 +51,10 @@ def test_simulate(loaded_frhodo): # They are reloaded from configuration variables when we re-run a simulation assert loaded_frhodo.window.display_shock['species_alias']['A'] == 'B' + # Test running the simulation from the keyword arguments + kwargs, rxn_cond = loaded_frhodo.get_simulator_inputs() + manual_sim = loaded_frhodo.run_simulation_from_kwargs(kwargs[0], rxn_cond[0]) + assert np.isclose(runs[0], manual_sim).all() @mark.parametrize( 'rxn_id, prs_id', [(0, 0), # PLog TODO (wardlt): Reaction 0 fails because we don't yet support PLog reactions @@ -137,13 +141,3 @@ def test_optimizer(loaded_frhodo, example_dir, tmp_path): y0_repeat = opt(x0) assert y0 == y0_repeat - # Make sure the serialization works - with Pool(2) as p: - # Make sure we still get repeatability - y = p.map(opt, [x0] * 32) - assert np.isclose(y, y0).all() - - # Run with different inputs and make sure they match up - y = p.map(opt, [x0, x1]) - assert np.isclose(y, [y0, y1]).all() - From 6876938cc9d60ae232f24a780cab0802924ac188 Mon Sep 17 00:00:00 2001 From: Logan Ward Date: Mon, 26 Sep 2022 08:54:59 -0400 Subject: [PATCH 24/29] No longer require modifying GUI to run simulation --- frhodo/api/driver.py | 88 +++++++---- frhodo/api/optimize.py | 75 ++++----- frhodo/calculate/mech_fcns.py | 276 ++++++++++++++++++++-------------- 3 files changed, 254 insertions(+), 185 deletions(-) diff --git a/frhodo/api/driver.py b/frhodo/api/driver.py index 62f4419..a5776fe 100644 --- a/frhodo/api/driver.py +++ b/frhodo/api/driver.py @@ -7,6 +7,7 @@ import numpy as np from PyQt5.QtWidgets import QApplication +from frhodo.calculate.mech_fcns import Chemical_Mechanism from frhodo.calculate.reactors import Simulation_Result from frhodo.main import Main, launch_gui @@ -14,6 +15,61 @@ """How to specify the index of a reaction parameter: reaction index, pressure index, name of the parameter""" +def set_coefficients(mech: Chemical_Mechanism, new_values: Dict[CoefIndex, float]): + """Get the coefficients of a chemical mechanism object + + Args: + mech: Mechanism to modify + new_values: New values for different coefficients + """ + for key, new_value in new_values.items(): + rxn_id, prs_id, coef_name = key + rxn_model = mech.coeffs[rxn_id][prs_id] + assert coef_name in rxn_model, f'Key {coef_name} not present for reaction #{rxn_id} at pressure #{prs_id}' + rxn_model[coef_name] = new_value + + # Update the reaction rates for the current shock + mech.modify_reactions(mech.coeffs) + + +def run_simulation(mech: Chemical_Mechanism, rxn_conditions: Tuple[float, float, dict], sim_kwargs: dict) -> np.ndarray: + """Run a simulation for a single reaction condition + + Args: + mech: Mechanism describing the chemical kinetics + rxn_conditions: Conditions of the reaction (temperature, pressure, composition) + sim_kwargs: Keywords to the simulator + Returns: + Array with time as first column and observable as the second + """ + mech.set_TPX(*rxn_conditions) + # Run the simulation + # TODO (wardlt): Do not hardcode the runtime or reactor conditions + sim, _ = mech.run('Incident Shock Reactor', 1.2e-05, *rxn_conditions, **sim_kwargs) + assert sim.success, "Simulation failed" + return np.stack([ + sim.independent_var, + sim.observable + ], axis=1) + + +def get_coefficients(mech: Chemical_Mechanism, indices: List[CoefIndex]) -> List[float]: + """Get specific coefficients from the reaction model + + Args: + mech: Mechanism to interrogate + indices: List of coefficients to retrieve + Returns: + Desired coefficents + """ + output = [] + for rxn_id, prs_id, coef_name in indices: + rxn_model = mech.coeffs[rxn_id][prs_id] + assert coef_name in rxn_model, f'Key {coef_name} not present for reaction #{rxn_id} at pressure #{prs_id}' + output.append(rxn_model[coef_name]) + return output + + class FrhodoDriver: """Driver the Frhodo GUI application @@ -175,17 +231,8 @@ def run_simulation_from_kwargs(self, sim_kwargs, rxn_conditions) -> np.ndarray: """ # Update the mechanism for the specified reaction conditions - self.window.mech.set_TPX(*rxn_conditions) - - # Run the simulation - # TODO (wardlt): Do not hardcode the runtime or reactor conditions - sim, _ = self.window.mech.run('Incident Shock Reactor', 1.2e-05, *rxn_conditions, **sim_kwargs) - - assert sim.success, "Simulation failed" - return np.stack([ - sim.independent_var, - sim.observable - ], axis=1) + mech = self.window.mech + return run_simulation(mech, rxn_conditions, sim_kwargs) def get_reaction_rates(self) -> np.ndarray: """Get the reaction rates for each shock experiment @@ -249,11 +296,8 @@ def get_coefficients(self, indices: List[CoefIndex]) -> List[float]: List of their current values """ - output = [] - for rxn_id, prs_id, coef_name in indices: - rxn_model = self.window.mech.coeffs[rxn_id][prs_id] - assert coef_name in rxn_model, f'Key {coef_name} not present for reaction #{rxn_id} at pressure #{prs_id}' - output.append(rxn_model[coef_name]) + mech = self.window.mech + output = get_coefficients(mech, indices) return output def set_coefficients(self, new_values: Dict[CoefIndex, float]): @@ -266,16 +310,8 @@ def set_coefficients(self, new_values: Dict[CoefIndex, float]): """ # Update the value in the Chemical_Mechanism dictionary - for key, new_value in new_values.items(): - rxn_id, prs_id, coef_name = key - rxn_model = self.window.mech.coeffs[rxn_id][prs_id] - assert coef_name in rxn_model, f'Key {coef_name} not present for reaction #{rxn_id} at pressure #{prs_id}' - rxn_model[coef_name] = new_value - - # Update the reaction rates for the current shock - self.window.mech.modify_reactions(self.window.mech.coeffs) - self.window.tree.update_rates() - self.app.processEvents() + mech = self.window.mech + set_coefficients(mech, new_values) @classmethod def create_driver(cls, headless: bool = True, fresh: bool = True) -> 'FrhodoDriver': diff --git a/frhodo/api/optimize.py b/frhodo/api/optimize.py index 7fe7a8e..6b5c1b0 100644 --- a/frhodo/api/optimize.py +++ b/frhodo/api/optimize.py @@ -8,9 +8,8 @@ from scipy import stats as ss from scipy.interpolate import interp1d -from frhodo.api.driver import CoefIndex, FrhodoDriver - -_frhodo: Optional[FrhodoDriver] = None +from frhodo.api.driver import CoefIndex, FrhodoDriver, set_coefficients, get_coefficients, run_simulation +from frhodo.calculate.mech_fcns import Chemical_Mechanism class BaseObjectiveFunction: @@ -43,63 +42,31 @@ def __init__( self._frhodo: Optional[FrhodoDriver] = None # Make a placeholder for experimental data + self.mech: Optional[Chemical_Mechanism] = None self._observations: Optional[List[np.ndarray]] = None self._weights: Optional[List[np.ndarray]] = None self._sim_kwargs: Optional[List[dict]] = None self._rxn_conditions: Optional[List[Tuple[float, float, dict]]] = None - # Track whether the data has been loaded or not - self._loaded = False - def __getstate__(self): if multiprocessing.get_start_method() != "spawn": warnings.warn('You must set the multiprocessing start method to spawn so we can run >1 multiprocessing instance') state = self.__dict__.copy() return state - def set_driver(self, driver: FrhodoDriver): - """Set the driver to be used for this class - - Use this if you have already launched Frhodo. - - Args: - driver: Driver to use in place of creating a new one - """ - global _frhodo - _frhodo = driver - - @property - def frhodo(self) -> FrhodoDriver: - """Link to the underlying Frhodo instance""" - global _frhodo # Use the module-level instance of Frhodo - # Launch if needed - if _frhodo is None: - _frhodo = FrhodoDriver.create_driver() - - # Load in the data and mechanism - if not self._loaded: - _frhodo.load_files( - self._exp_directory, - self._mech_directory, - self._mech_directory / 'outputs', - aliases=self._aliases - ) - assert self.observations is not None - return _frhodo - @property def observations(self) -> List[np.ndarray]: """Experimental data less any regions that are masked out""" # Load the data in if needed if self._observations is None: - self._load_experiments() + self.load_experiments() return self._observations @property def weights(self) -> List[np.ndarray]: """Weights for each observation""" if self._observations is None: - self._load_experiments() + self.load_experiments() return self._weights @property @@ -108,16 +75,29 @@ def x(self) -> np.ndarray: This is either the initial values from loading the mechanism or the last values ran. """ - return np.array(self.frhodo.get_coefficients(self.parameters)) + return np.array(get_coefficients(self.mech, self.parameters)) - def _load_experiments(self): + def load_experiments(self, frhodo: Optional[FrhodoDriver] = None): """Load observations and weights from disk Sets the values in `self._observations` and `self._weights` """ + + if frhodo is None: + frhodo = FrhodoDriver.create_driver() + + # Use Frhodo to load the data + frhodo.load_files( + self._exp_directory, + self._mech_directory, + self._mech_directory / 'outputs', + aliases=self._aliases + ) + + # Extract the data self._observations = [] self._weights = [] - for obs, weights in zip(*self.frhodo.get_observables()): + for obs, weights in zip(*frhodo.get_observables()): # Find portions that are weighted sufficiently mask = weights > 1e-6 @@ -125,7 +105,8 @@ def _load_experiments(self): self._observations.append(obs[mask, :]) self._weights.append(weights[mask]) - self._sim_kwargs, self._rxn_conditions = self.frhodo.get_simulator_inputs() + self._sim_kwargs, self._rxn_conditions = frhodo.get_simulator_inputs() + self.mech = frhodo.window.mech def run_simulations(self, x: Sequence[float]) -> List[np.ndarray]: """Run the simulations with a new set of parameters @@ -139,17 +120,17 @@ def run_simulations(self, x: Sequence[float]) -> List[np.ndarray]: # Update the parameters assert len(x) == len(self.parameters), f"Expected {len(self.parameters)} parameters but got {len(x)}" - self.frhodo.set_coefficients(dict(zip(self.parameters, x))) + set_coefficients(self.mech, dict(zip(self.parameters, x))) # Run each simulation to each set of experimental data sims = [] for sim_kwargs, rxn_cond in zip(self._sim_kwargs, self._rxn_conditions): - sims.append(self.frhodo.run_simulation_from_kwargs(sim_kwargs, rxn_cond)) + sims.append(run_simulation(self.mech, rxn_cond, sim_kwargs)) # Interpolate simulation data over the same steps as the experiments output = [] for sim, obs in zip(sims, self.observations): - sim_func = interp1d(sim[:, 0], sim[:, 1], kind='cubic') + sim_func = interp1d(sim[:, 0], sim[:, 1], kind='cubic', fill_value='extrapolate') output.append(sim_func(obs[:, 0])) return output @@ -225,8 +206,8 @@ def compute_log_probs(self, x: np.ndarray, uncertainty: float) -> np.ndarray: output.append(ss.t(loc=0, scale=std, df=2.1).logpdf(resid).sum()) return np.array(output) - def _load_experiments(self): - super()._load_experiments() + def load_experiments(self, frhodo: Optional[FrhodoDriver] = None): + super().load_experiments(frhodo) # Adjust the weights such that the largest is 1 for i, weight in enumerate(self.weights): diff --git a/frhodo/calculate/mech_fcns.py b/frhodo/calculate/mech_fcns.py index 20436b7..195f2ad 100644 --- a/frhodo/calculate/mech_fcns.py +++ b/frhodo/calculate/mech_fcns.py @@ -4,16 +4,18 @@ import os, io, stat, contextlib, pathlib, time from copy import deepcopy +from pathlib import Path +from tempfile import NamedTemporaryFile, TemporaryDirectory from typing import Tuple, Optional, Sequence, Union import cantera as ct -from cantera import interrupts, cti2yaml#, ck2yaml, ctml2yaml +from cantera import interrupts, cti2yaml # , ck2yaml, ctml2yaml import numpy as np from timeit import default_timer as timer from . import reactors, shock_fcns, integrate from .reactors import Simulation_Result -from .. import ck2yaml +from .. import ck2yaml, soln2ck class Chemical_Mechanism: @@ -21,35 +23,68 @@ def __init__(self): self.isLoaded = False self.reactor = reactors.Reactor(self) + def __getstate__(self): + state = self.__dict__.copy() + # Save the current state of the gas as Chemkin Mech file + # TODO: Figure out why we cannot save it to YML + with NamedTemporaryFile() as fp: + fp.close() + soln2ck.write(self.gas, fp.name) + with open(fp.name) as fp: + state['gas'] = fp.read() + return state + + def __setstate__(self, state): + # Replace "gas" with a Cantera object + with TemporaryDirectory() as tmp: + # Convert the ChemKin output to YAML + # First write the mechanism to disk as a chemkin file + tmp = Path(tmp) + ck_file = tmp / 'gas.mech' + with ck_file.open('w') as fp: + print(state.pop('gas'), file=fp) + + # Now convert it to a YAML + ct_file = tmp / 'gas.yaml' + ck2yaml.convert_mech(ck_file, thermo_file=None, transport_file=None, surface_file=None, + phase_name='gas', out_name=ct_file, quiet=False, + permissive=True) + state['gas'] = ct.Solution(yaml=ct_file.read_text()) + + self.__dict__.update(state) + def load_mechanism(self, path, silent=False): def chemkin2cantera(path): if path['thermo'] is not None: - surfaces = ck2yaml.convert_mech(path['mech'], thermo_file=path['thermo'], transport_file=None, surface_file=None, - phase_name='gas', out_name=path['Cantera_Mech'], quiet=False, permissive=True) + surfaces = ck2yaml.convert_mech(path['mech'], thermo_file=path['thermo'], transport_file=None, + surface_file=None, + phase_name='gas', out_name=path['Cantera_Mech'], quiet=False, + permissive=True) else: surfaces = ck2yaml.convert_mech(path['mech'], thermo_file=None, transport_file=None, surface_file=None, - phase_name='gas', out_name=path['Cantera_Mech'], quiet=False, permissive=True) - + phase_name='gas', out_name=path['Cantera_Mech'], quiet=False, + permissive=True) + return surfaces def loader(self, path): # path is assumed to be the path dictionary surfaces = [] - if path['mech'].suffix in ['.yaml', '.yml']: # check if it's a yaml cantera file + if path['mech'].suffix in ['.yaml', '.yml']: # check if it's a yaml cantera file mech_path = str(path['mech']) - else: # if not convert into yaml cantera file + else: # if not convert into yaml cantera file mech_path = str(path['Cantera_Mech']) - + if path['mech'].suffix == '.cti': cti2yaml.convert(path['mech'], path['Cantera_Mech']) elif path['mech'].suffix in ['.ctml', '.xml']: raise Exception('not implemented') - #ctml2yaml.convert(path['mech'], path['Cantera_Mech']) - else: # if not a cantera file, assume chemkin + # ctml2yaml.convert(path['mech'], path['Cantera_Mech']) + else: # if not a cantera file, assume chemkin surfaces = chemkin2cantera(path) - - print('Validating mechanism...', end='') - try: # This test taken from ck2cti + + print('Validating mechanism...', end='') + try: # This test taken from ck2cti yaml_txt = path['Cantera_Mech'].read_text() self.gas = ct.Solution(yaml=yaml_txt) for surfname in surfaces: @@ -58,7 +93,7 @@ def loader(self, path): except RuntimeError as e: print('FAILED.') print(e) - + output = {'success': False, 'message': []} # Initialize and report any problems to log, not to console window stdout = io.StringIO() @@ -73,31 +108,32 @@ def loader(self, path): except: pass # output['message'].append('Error when loading mech:\n') - + ct_out = stdout.getvalue() ct_err = stderr.getvalue().replace('INFO:root:', 'Warning: ') - + if 'FAILED' in ct_out: output['success'] = False self.isLoaded = False elif 'PASSED' in ct_out: output['success'] = True self.isLoaded = True - + for log_str in [ct_out, ct_err]: if log_str != '' and not silent: - if (path['Cantera_Mech'], pathlib.WindowsPath): # reformat string to remove \\ making it unable to be copy paste + if (path['Cantera_Mech'], + pathlib.WindowsPath): # reformat string to remove \\ making it unable to be copy paste cantera_path = str(path['Cantera_Mech']).replace('\\', '\\\\') log_str = log_str.replace(cantera_path, str(path['Cantera_Mech'])) output['message'].append(log_str) output['message'].append('\n') - + if self.isLoaded: - self.set_rate_expression_coeffs() # set copy of coeffs - self.set_thermo_expression_coeffs() # set copy of thermo coeffs + self.set_rate_expression_coeffs() # set copy of coeffs + self.set_thermo_expression_coeffs() # set copy of thermo coeffs return output - + def set_mechanism(self, mech_dict, species_dict={}, bnds=[]): def get_Arrhenius_parameters(entry): A = entry['pre_exponential_factor'] @@ -113,12 +149,12 @@ def get_Arrhenius_parameters(entry): for n in range(len(species_dict)): s_dict = species_dict[n] s = ct.Species(name=s_dict['name'], composition=s_dict['composition'], - charge=s_dict['charge'], size=s_dict['size']) + charge=s_dict['charge'], size=s_dict['size']) thermo = s_dict['type'](s_dict['T_low'], s_dict['T_high'], s_dict['P_ref'], s_dict['coeffs']) s.thermo = thermo species.append(s) - + # Set kinetics data rxns = [] for rxnIdx in range(len(mech_dict)): @@ -172,25 +208,26 @@ def get_Arrhenius_parameters(entry): elif 'ChebyshevReaction' == mech_dict[rxnIdx]['rxnType']: rxn = ct.ChebyshevReaction(mech_dict[rxnIdx]['reactants'], mech_dict[rxnIdx]['products']) - rxn.set_parameters(Tmin=mech_dict['Tmin'], Tmax=mech_dict['Tmax'], + rxn.set_parameters(Tmin=mech_dict['Tmin'], Tmax=mech_dict['Tmax'], Pmin=mech_dict['Pmin'], Pmax=mech_dict['Pmax'], coeffs=mech_dict['coeffs']) - + rxn.duplicate = mech_dict[rxnIdx]['duplicate'] rxn.reversible = mech_dict[rxnIdx]['reversible'] rxn.allow_negative_orders = True rxn.allow_nonreactant_orders = True rxns.append(rxn) - + self.gas = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', species=species, reactions=rxns) - - self.set_rate_expression_coeffs(bnds) # set copy of coeffs - self.set_thermo_expression_coeffs() # set copy of thermo coeffs - - def gas(self): return self.gas - + + self.set_rate_expression_coeffs(bnds) # set copy of coeffs + self.set_thermo_expression_coeffs() # set copy of thermo coeffs + + def gas(self): + return self.gas + def set_rate_expression_coeffs(self, bnds=[]): def copy_bnds(new_bnds, bnds, rxnIdx, bnds_type, keys=[]): if len(bnds) == 0: return new_bnds @@ -198,7 +235,7 @@ def copy_bnds(new_bnds, bnds, rxnIdx, bnds_type, keys=[]): if bnds_type == 'rate': for key in ['value', 'type', 'opt']: new_bnds[rxnIdx][key] = bnds['rate_bnds'][rxnIdx][key] - + else: bndsKey, attrs = keys for coefName in attrs: @@ -213,89 +250,103 @@ def copy_bnds(new_bnds, bnds, rxnIdx, bnds_type, keys=[]): self.reset_mech = reset_mech = [] for rxnIdx, rxn in enumerate(self.gas.reactions()): - rate_bnds.append({'value': np.nan, 'limits': Uncertainty('rate', rxnIdx, rate_bnds=rate_bnds), 'type': 'F', 'opt': False}) + rate_bnds.append({'value': np.nan, 'limits': Uncertainty('rate', rxnIdx, rate_bnds=rate_bnds), 'type': 'F', + 'opt': False}) rate_bnds = copy_bnds(rate_bnds, bnds, rxnIdx, 'rate') if type(rxn) in [ct.ElementaryReaction, ct.ThreeBodyReaction]: - attrs = [p for p in dir(rxn.rate) if not p.startswith('_')] # attributes not including __ + attrs = [p for p in dir(rxn.rate) if not p.startswith('_')] # attributes not including __ coeffs.append([{attr: getattr(rxn.rate, attr) for attr in attrs}]) if type(rxn) is ct.ThreeBodyReaction: coeffs[-1][0]['efficiencies'] = rxn.efficiencies - coeffs_bnds.append({'rate': {attr: {'resetVal': coeffs[-1][0][attr], 'value': np.nan, - 'limits': Uncertainty('coef', rxnIdx, key='rate', coef_name=attr, coeffs_bnds=coeffs_bnds), - 'type': 'F', 'opt': False} for attr in attrs}}) - + coeffs_bnds.append({'rate': {attr: {'resetVal': coeffs[-1][0][attr], 'value': np.nan, + 'limits': Uncertainty('coef', rxnIdx, key='rate', coef_name=attr, + coeffs_bnds=coeffs_bnds), + 'type': 'F', 'opt': False} for attr in attrs}}) + coeffs_bnds = copy_bnds(coeffs_bnds, bnds, rxnIdx, 'coeffs', ['rate', attrs]) - reset_mech.append({'reactants': rxn.reactants, 'products': rxn.products, 'rxnType': rxn.__class__.__name__, - 'duplicate': rxn.duplicate, 'reversible': rxn.reversible, 'orders': rxn.orders, - 'rxnCoeffs': deepcopy(coeffs[-1])}) - + reset_mech.append( + {'reactants': rxn.reactants, 'products': rxn.products, 'rxnType': rxn.__class__.__name__, + 'duplicate': rxn.duplicate, 'reversible': rxn.reversible, 'orders': rxn.orders, + 'rxnCoeffs': deepcopy(coeffs[-1])}) + elif type(rxn) is ct.PlogReaction: coeffs.append([]) coeffs_bnds.append({}) for n, rate in enumerate(rxn.rates): - attrs = [p for p in dir(rate[1]) if not p.startswith('_')] # attributes not including __ + attrs = [p for p in dir(rate[1]) if not p.startswith('_')] # attributes not including __ coeffs[-1].append({'Pressure': rate[0]}) coeffs[-1][-1].update({attr: getattr(rate[1], attr) for attr in attrs}) - if n == 0 or n == len(rxn.rates)-1: # only going to allow coefficient uncertainties to be placed on upper and lower pressures + if n == 0 or n == len( + rxn.rates) - 1: # only going to allow coefficient uncertainties to be placed on upper and lower pressures if n == 0: key = 'low_rate' else: key = 'high_rate' - coeffs_bnds[-1][key] = {attr: {'resetVal': coeffs[-1][-1][attr], 'value': np.nan, - 'limits': Uncertainty('coef', rxnIdx, key=key, coef_name=attr, coeffs_bnds=coeffs_bnds), - 'type': 'F', 'opt': False} for attr in attrs} + coeffs_bnds[-1][key] = {attr: {'resetVal': coeffs[-1][-1][attr], 'value': np.nan, + 'limits': Uncertainty('coef', rxnIdx, key=key, coef_name=attr, + coeffs_bnds=coeffs_bnds), + 'type': 'F', 'opt': False} for attr in attrs} coeffs_bnds = copy_bnds(coeffs_bnds, bnds, rxnIdx, 'coeffs', [key, attrs]) - reset_mech.append({'reactants': rxn.reactants, 'products': rxn.products, 'rxnType': rxn.__class__.__name__, - 'duplicate': rxn.duplicate, 'reversible': rxn.reversible, 'orders': rxn.orders, - 'rxnCoeffs': deepcopy(coeffs[-1])}) + reset_mech.append( + {'reactants': rxn.reactants, 'products': rxn.products, 'rxnType': rxn.__class__.__name__, + 'duplicate': rxn.duplicate, 'reversible': rxn.reversible, 'orders': rxn.orders, + 'rxnCoeffs': deepcopy(coeffs[-1])}) elif type(rxn) is ct.FalloffReaction: coeffs_bnds.append({}) - coeffs.append({'falloff_type': rxn.falloff.type, 'high_rate': [], 'low_rate': [], 'falloff_parameters': list(rxn.falloff.parameters), + coeffs.append({'falloff_type': rxn.falloff.type, 'high_rate': [], 'low_rate': [], + 'falloff_parameters': list(rxn.falloff.parameters), 'default_efficiency': rxn.default_efficiency, 'efficiencies': rxn.efficiencies}) for key in ['low_rate', 'high_rate']: rate = getattr(rxn, key) - attrs = [p for p in dir(rate) if not p.startswith('_')] # attributes not including __ + attrs = [p for p in dir(rate) if not p.startswith('_')] # attributes not including __ coeffs[-1][key] = {attr: getattr(rate, attr) for attr in attrs} - coeffs_bnds[-1][key] = {attr: {'resetVal': coeffs[-1][key][attr], 'value': np.nan, - 'limits': Uncertainty('coef', rxnIdx, key=key, coef_name=attr, coeffs_bnds=coeffs_bnds), - 'type': 'F', 'opt': False} for attr in attrs} + coeffs_bnds[-1][key] = {attr: {'resetVal': coeffs[-1][key][attr], 'value': np.nan, + 'limits': Uncertainty('coef', rxnIdx, key=key, coef_name=attr, + coeffs_bnds=coeffs_bnds), + 'type': 'F', 'opt': False} for attr in attrs} coeffs_bnds = copy_bnds(coeffs_bnds, bnds, rxnIdx, 'coeffs', [key, attrs]) key = 'falloff_parameters' n_coef = len(rxn.falloff.parameters) - coeffs_bnds[-1][key] = {n: {'resetVal': coeffs[-1][key][n], 'value': np.nan, - 'limits': Uncertainty('coef', rxnIdx, key=key, coef_name=n, coeffs_bnds=coeffs_bnds), - 'type': 'F', 'opt': True} for n in range(0,n_coef)} - - reset_mech.append({'reactants': rxn.reactants, 'products': rxn.products, 'rxnType': rxn.__class__.__name__, - 'duplicate': rxn.duplicate, 'reversible': rxn.reversible, 'orders': rxn.orders, - 'falloffType': rxn.falloff.type, 'rxnCoeffs': deepcopy(coeffs[-1])}) - + coeffs_bnds[-1][key] = {n: {'resetVal': coeffs[-1][key][n], 'value': np.nan, + 'limits': Uncertainty('coef', rxnIdx, key=key, coef_name=n, + coeffs_bnds=coeffs_bnds), + 'type': 'F', 'opt': True} for n in range(0, n_coef)} + + reset_mech.append( + {'reactants': rxn.reactants, 'products': rxn.products, 'rxnType': rxn.__class__.__name__, + 'duplicate': rxn.duplicate, 'reversible': rxn.reversible, 'orders': rxn.orders, + 'falloffType': rxn.falloff.type, 'rxnCoeffs': deepcopy(coeffs[-1])}) + elif type(rxn) is ct.ChebyshevReaction: coeffs.append({}) coeffs_bnds.append({}) if len(bnds) == 0: rate_bnds.append({}) - - reset_coeffs = {'Pmin': rxn.Pmin, 'Pmax': rxn.Pmax, 'Tmin': rxn.Tmin, 'Tmax': rxn.Tmax, 'coeffs': rxn.coeffs} - reset_mech.append({'reactants': rxn.reactants, 'products': rxn.products, 'rxnType': rxn.__class__.__name__, - 'duplicate': rxn.duplicate, 'reversible': rxn.reversible, 'orders': rxn.orders, - 'rxnCoeffs': reset_coeffs}) + + reset_coeffs = {'Pmin': rxn.Pmin, 'Pmax': rxn.Pmax, 'Tmin': rxn.Tmin, 'Tmax': rxn.Tmax, + 'coeffs': rxn.coeffs} + reset_mech.append( + {'reactants': rxn.reactants, 'products': rxn.products, 'rxnType': rxn.__class__.__name__, + 'duplicate': rxn.duplicate, 'reversible': rxn.reversible, 'orders': rxn.orders, + 'rxnCoeffs': reset_coeffs}) else: coeffs.append({}) coeffs_bnds.append({}) if len(bnds) == 0: rate_bnds.append({}) - reset_mech.append({'reactants': rxn.reactants, 'products': rxn.products, 'rxnType': rxn.__class__.__name__}) - raise(f'{rxn} is a {rxn.__class__.__name__} and is currently unsupported in Frhodo, but this error should never be seen') + reset_mech.append( + {'reactants': rxn.reactants, 'products': rxn.products, 'rxnType': rxn.__class__.__name__}) + raise ( + f'{rxn} is a {rxn.__class__.__name__} and is currently unsupported in Frhodo, but this error should never be seen') def get_coeffs_keys(self, rxn, coefAbbr, rxnIdx=None): """Get the keys in a reaction dictionary associated that define the par""" @@ -309,13 +360,13 @@ def get_coeffs_keys(self, rxn, coefAbbr, rxnIdx=None): for rxnIdx, mechRxn in enumerate(self.gas.reactions()): if rxn is mechRxn: break - + bnds_key = 'high_rate' coef_key = len(self.coeffs[rxnIdx]) - 1 elif 'low' in coefAbbr: bnds_key = 'low_rate' coef_key = 0 - + elif type(rxn) is ct.FalloffReaction: if 'high' in coefAbbr: coef_key = bnds_key = 'high_rate' @@ -329,7 +380,7 @@ def get_coeffs_keys(self, rxn, coefAbbr, rxnIdx=None): return coef_key, bnds_key - def set_thermo_expression_coeffs(self): # TODO Doesn't work with NASA 9 + def set_thermo_expression_coeffs(self): # TODO Doesn't work with NASA 9 self.thermo_coeffs = [] for i in range(self.gas.n_species): S = self.gas.species(i) @@ -337,8 +388,8 @@ def set_thermo_expression_coeffs(self): # TODO Doesn't work with NASA 9 'size': S.size, 'type': type(S.thermo), 'P_ref': S.thermo.reference_pressure, 'T_low': S.thermo.min_temp, 'T_high': S.thermo.max_temp, 'coeffs': np.array(S.thermo.coeffs), - 'h_scaler': 1, 's_scaler': 1,} - + 'h_scaler': 1, 's_scaler': 1, } + self.thermo_coeffs.append(thermo_dict) def modify_reactions(self, coeffs, rxnIdxs: Optional[Union[int, float, Sequence[int]]] = ()): @@ -350,12 +401,12 @@ def modify_reactions(self, coeffs, rxnIdxs: Optional[Union[int, float, Sequence[ """ # TODO (wardlt): This function is only called with `self.mech` as inputs # TODO (sikes): Only works for Arrhenius equations currently - if not rxnIdxs: # if rxnNums does not exist, modify all + if not rxnIdxs: # if rxnNums does not exist, modify all rxnIdxs = range(len(coeffs)) else: if isinstance(rxnIdxs, (float, int)): # if single reaction given, run that onee rxnIdxs = [rxnIdxs] - + for rxnIdx in rxnIdxs: rxn = self.gas.reaction(rxnIdx) rxnChanged = False @@ -363,8 +414,8 @@ def modify_reactions(self, coeffs, rxnIdxs: Optional[Union[int, float, Sequence[ for coefName in ['activation_energy', 'pre_exponential_factor', 'temperature_exponent']: if coeffs[rxnIdx][0][coefName] != eval(f'rxn.rate.{coefName}'): rxnChanged = True - - if rxnChanged: # Update reaction rate + + if rxnChanged: # Update reaction rate A = coeffs[rxnIdx][0]['pre_exponential_factor'] b = coeffs[rxnIdx][0]['temperature_exponent'] Ea = coeffs[rxnIdx][0]['activation_energy'] @@ -392,36 +443,36 @@ def modify_reactions(self, coeffs, rxnIdxs: Optional[Union[int, float, Sequence[ rxn.falloff = ct.TroeFalloff(coeffs[rxnIdx][key][:-1]) else: rxn.falloff = ct.TroeFalloff(coeffs[rxnIdx][key]) - else: # could also be SRI. For optimization this would need to be cast as Troe + else: # could also be SRI. For optimization this would need to be cast as Troe rxn.falloff = ct.SriFalloff(coeffs[rxnIdx][key]) - elif type(rxn) is ct.ChebyshevReaction: - pass + elif type(rxn) is ct.ChebyshevReaction: + pass else: continue - + if rxnChanged: self.gas.modify_reaction(rxnIdx, rxn) # Not sure if this is necessary, but it reduces strange behavior in incident shock reactor - #time.sleep(10E-3) # TODO: if incident shock reactor is written in C++, this can likely be removed - + # time.sleep(10E-3) # TODO: if incident shock reactor is written in C++, this can likely be removed + def rxn2Troe(self, rxnIdx, HPL, LPL, eff={}): reactants = self.gas.reaction(rxnIdx).reactants products = self.gas.reaction(rxnIdx).products r = ct.FalloffReaction(reactants, products) print(r) - #r.high_rate = ct.Arrhenius(7.4e10, -0.37, 0.0) - #r.low_rate = ct.Arrhenius(2.3e12, -0.9, -1700*1000*4.184) - #r.falloff = ct.TroeFalloff((0.7346, 94, 1756, 5182)) - #r.efficiencies = {'AR':0.7, 'H2':2.0, 'H2O':6.0} + # r.high_rate = ct.Arrhenius(7.4e10, -0.37, 0.0) + # r.low_rate = ct.Arrhenius(2.3e12, -0.9, -1700*1000*4.184) + # r.falloff = ct.TroeFalloff((0.7346, 94, 1756, 5182)) + # r.efficiencies = {'AR':0.7, 'H2':2.0, 'H2O':6.0} print(dir(self.gas)) print(self.gas.thermo_model) print(self.gas.kinetics_model) start = timer() - #self.gas.thermo_model - #self.gas.kinetics_model + # self.gas.thermo_model + # self.gas.kinetics_model self.gas = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', species=self.gas.species(), reactions=self.gas.reactions()) print(timer() - start) @@ -440,7 +491,7 @@ def modify_thermo(self, multipliers): # Only works for NasaPoly2 (NASA 7) curre T_high = S_initial.thermo.max_temp P_ref = S_initial.thermo.reference_pressure coeffs = S_initial.thermo.coeffs - + # Update thermo properties coeffs[1:] *= multipliers[i] S.thermo = ct.NasaPoly2(T_low, T_high, P_ref, coeffs) @@ -454,13 +505,13 @@ def modify_thermo(self, multipliers): # Only works for NasaPoly2 (NASA 7) curre continue self.gas.modify_species(i, S) - + def reset(self, rxnIdxs=None, coefNames=None): if rxnIdxs is None: rxnIdxs = range(self.gas.n_reactions) - elif type(rxnIdxs) is not list: # if not list then assume given single rxnIdx + elif type(rxnIdxs) is not list: # if not list then assume given single rxnIdx rxnIdxs = [rxnIdxs] - + # if not list then assume given single coefficient # if Arrhenius, expects list of coefficients. If Plog or Falloff expected [['high_rate', 'activation_energy'], ...] if coefNames is not None and type(coefNames) is not list: @@ -468,9 +519,9 @@ def reset(self, rxnIdxs=None, coefNames=None): prior_coeffs = deepcopy(self.coeffs) for rxnIdx in rxnIdxs: - if coefNames is None: # resets all coefficients in rxn + if coefNames is None: # resets all coefficients in rxn self.coeffs[rxnIdx] = self.reset_mech[rxnIdx]['rxnCoeffs'] - + elif self.reset_mech[rxnIdx]['rxnType'] in ['ElementaryReaction', 'ThreeBodyReaction']: for coefName in coefNames: self.coeffs[rxnIdx][coefName] = self.reset_mech[rxnIdx]['rxnCoeffs'][coefName] @@ -485,7 +536,8 @@ def reset(self, rxnIdxs=None, coefNames=None): elif 'FalloffReaction' == self.reset_mech[rxnIdx]['rxnType']: self.coeffs[rxnIdx]['falloff_type'] = self.reset_mech[rxnIdx]['falloffType'] for [limit_type, coefName] in coefNames: - self.coeffs[rxnIdx][limit_type][coefName] = self.reset_mech[rxnIdx]['rxnCoeffs'][limit_type][coefName] + self.coeffs[rxnIdx][limit_type][coefName] = self.reset_mech[rxnIdx]['rxnCoeffs'][limit_type][ + coefName] self.modify_reactions(self.coeffs) @@ -506,23 +558,23 @@ def set_TPX(self, T, P, X=[]): if species not in self.gas.species_names: output['message'].append('Species: {:s} is not in the mechanism'.format(species)) return output - + self.gas.TPX = T, P, X else: self.gas.TP = T, P - + output['success'] = True return output - def M(self, rxn, TPX=[]): # kmol/m^3 + def M(self, rxn, TPX=[]): # kmol/m^3 def get_M(rxn): M = self.gas.density_mole if hasattr(rxn, 'efficiencies') and rxn.efficiencies: M *= rxn.default_efficiency for (s, conc) in zip(self.gas.species_names, self.gas.concentrations): if s in rxn.efficiencies: - M += conc*(rxn.efficiencies[s] - 1.0) + M += conc * (rxn.efficiencies[s] - 1.0) else: M += conc return M @@ -533,9 +585,9 @@ def get_M(rxn): [T, P, X] = TPX M = [] for i in range(0, len(T)): - self.set_TPX(T[i], P[i], X) # IF MIXTURE CHANGES THEN THIS NEEDS TO BE VARIABLE + self.set_TPX(T[i], P[i], X) # IF MIXTURE CHANGES THEN THIS NEEDS TO BE VARIABLE M.append(get_M(rxn)) - + M = np.array(M) return M @@ -557,7 +609,7 @@ def run(self, reactor_choice, t_end, T_reac, P_reac, mix, **kwargs) -> Tuple[Sim return self.reactor.run(reactor_choice, t_end, T_reac, P_reac, mix, **kwargs) -class Uncertainty: # alternate name: why I hate pickle part 10 +class Uncertainty: # alternate name: why I hate pickle part 10 """Computes the uncertainty bounds for parameters""" def __init__(self, unc_type, rxnIdx, **kwargs): @@ -585,9 +637,9 @@ def _unc_fcn(self, x, uncVal, uncType): # TODO (wardlt): Can be static or if np.isnan(uncVal): return [np.nan, np.nan] elif uncType == 'F': - return np.sort([x/uncVal, x*uncVal], axis=0) + return np.sort([x / uncVal, x * uncVal], axis=0) elif uncType == '%': - return np.sort([x/(1 + uncVal), x*(1 + uncVal)], axis=0) + return np.sort([x / (1 + uncVal), x * (1 + uncVal)], axis=0) elif uncType == '±': return np.sort([x - uncVal, x + uncVal], axis=0) elif uncType == '+': @@ -599,7 +651,7 @@ def _unc_fcn(self, x, uncVal, uncType): # TODO (wardlt): Can be static or def __call__(self, x=None): if self.unc_type == 'rate': - #if x is None: # defaults to giving current rate bounds + # if x is None: # defaults to giving current rate bounds # x = self.gas.forward_rate_constants[self.rxnIdx] rate_bnds = self.unc_dict['rate_bnds'] unc_value = rate_bnds[self.rxnIdx]['value'] @@ -610,7 +662,7 @@ def __call__(self, x=None): key = self.unc_dict['key'] coefName = self.unc_dict['coef_name'] - if key == 'falloff_parameters': # falloff parameters have no limits + if key == 'falloff_parameters': # falloff parameters have no limits return [np.nan, np.nan] coef_dict = coeffs_bnds[self.rxnIdx][key][coefName] From a58c3ddf9cb36d96941099e37711c08644f7afe3 Mon Sep 17 00:00:00 2001 From: Logan Ward Date: Mon, 26 Sep 2022 09:30:14 -0400 Subject: [PATCH 25/29] Detect when cantera segfaults by running in separate process --- frhodo/api/optimize.py | 29 +++++++++++++++++++---------- frhodo/calculate/mech_fcns.py | 2 +- tests/test_api.py | 11 +++++++---- 3 files changed, 27 insertions(+), 15 deletions(-) diff --git a/frhodo/api/optimize.py b/frhodo/api/optimize.py index 6b5c1b0..bc255ea 100644 --- a/frhodo/api/optimize.py +++ b/frhodo/api/optimize.py @@ -1,6 +1,8 @@ """Utilities useful for using Frhodo from external optimizers""" import multiprocessing import warnings +from concurrent.futures import ProcessPoolExecutor +from concurrent.futures.process import BrokenProcessPool from pathlib import Path from typing import List, Optional, Dict, Sequence, Tuple @@ -48,12 +50,6 @@ def __init__( self._sim_kwargs: Optional[List[dict]] = None self._rxn_conditions: Optional[List[Tuple[float, float, dict]]] = None - def __getstate__(self): - if multiprocessing.get_start_method() != "spawn": - warnings.warn('You must set the multiprocessing start method to spawn so we can run >1 multiprocessing instance') - state = self.__dict__.copy() - return state - @property def observations(self) -> List[np.ndarray]: """Experimental data less any regions that are masked out""" @@ -124,8 +120,14 @@ def run_simulations(self, x: Sequence[float]) -> List[np.ndarray]: # Run each simulation to each set of experimental data sims = [] - for sim_kwargs, rxn_cond in zip(self._sim_kwargs, self._rxn_conditions): - sims.append(run_simulation(self.mech, rxn_cond, sim_kwargs)) + with ProcessPoolExecutor(1) as exec: + for sim_kwargs, rxn_cond in zip(self._sim_kwargs, self._rxn_conditions): + exc = exec.submit(run_simulation, self.mech, rxn_cond, sim_kwargs) + try: + res = exc.result() + except BrokenProcessPool: + raise ValueError('Running the simulation failed') + sims.append(res) # Interpolate simulation data over the same steps as the experiments output = [] @@ -152,10 +154,17 @@ def compute_residuals(self, x: Sequence[float]) -> List[np.ndarray]: output.append(np.subtract(sim, obs[:, 1])) return output - def __call__(self, x: np.ndarray, **kwargs): + def _call(self, x: np.ndarray, **kwargs): """Invoke the objective function""" raise NotImplementedError() + def __call__(self, x: np.ndarray, **kwargs): + """Invoke the objective function""" + try: + return self._call(x) + except ValueError: + return np.inf + class BayesianObjectiveFunction(BaseObjectiveFunction): """Computes the log-probability of observing experimental data given the simulated results @@ -213,7 +222,7 @@ def load_experiments(self, frhodo: Optional[FrhodoDriver] = None): for i, weight in enumerate(self.weights): self.weights[i] /= weight.max() - def __call__(self, x: np.ndarray, **kwargs): + def _call(self, x: np.ndarray, **kwargs): """Invoke the objective function Args: diff --git a/frhodo/calculate/mech_fcns.py b/frhodo/calculate/mech_fcns.py index 195f2ad..0b4c438 100644 --- a/frhodo/calculate/mech_fcns.py +++ b/frhodo/calculate/mech_fcns.py @@ -47,7 +47,7 @@ def __setstate__(self, state): # Now convert it to a YAML ct_file = tmp / 'gas.yaml' ck2yaml.convert_mech(ck_file, thermo_file=None, transport_file=None, surface_file=None, - phase_name='gas', out_name=ct_file, quiet=False, + phase_name='gas', out_name=ct_file, quiet=True, permissive=True) state['gas'] = ct.Solution(yaml=ct_file.read_text()) diff --git a/tests/test_api.py b/tests/test_api.py index e5c077c..8d19c64 100644 --- a/tests/test_api.py +++ b/tests/test_api.py @@ -1,11 +1,11 @@ """Testing the API components of Frhodo""" import numpy as np +import pickle as pkl from pytest import fixture, mark -from multiprocessing import Pool, set_start_method +from multiprocessing import set_start_method from frhodo.api.driver import FrhodoDriver from frhodo.api.optimize import BayesianObjectiveFunction -from frhodo.api import optimize @fixture @@ -113,7 +113,6 @@ def test_fittable_parameters(loaded_frhodo): def test_optimizer(loaded_frhodo, example_dir, tmp_path): - set_start_method("spawn") # Allows us to run >1 Frhodo instance opt = BayesianObjectiveFunction( exp_directory=example_dir / 'Experiment', mech_directory=example_dir / 'Mechanism', @@ -121,7 +120,7 @@ def test_optimizer(loaded_frhodo, example_dir, tmp_path): ) # Set the frhodo executable for that module (we can have only 1 per process) - optimize._frhodo = loaded_frhodo + opt.load_experiments(loaded_frhodo) # Test the state assert opt.weights[0].max() == 1 @@ -141,3 +140,7 @@ def test_optimizer(loaded_frhodo, example_dir, tmp_path): y0_repeat = opt(x0) assert y0 == y0_repeat + # Make sure it is serializable + opt2 = pkl.loads(pkl.dumps(opt)) + y0_opt2 = opt2(x0) + assert np.isclose(y0, y0_opt2).all() From e5df03facbca81d903f30b458e8747eb7346620c Mon Sep 17 00:00:00 2001 From: Logan Ward Date: Mon, 26 Sep 2022 10:52:58 -0400 Subject: [PATCH 26/29] Avoid re-constituting the process pool between invocations --- frhodo/api/optimize.py | 41 ++++++++++++++++++++++++++++++----------- 1 file changed, 30 insertions(+), 11 deletions(-) diff --git a/frhodo/api/optimize.py b/frhodo/api/optimize.py index bc255ea..d38592a 100644 --- a/frhodo/api/optimize.py +++ b/frhodo/api/optimize.py @@ -3,6 +3,7 @@ import warnings from concurrent.futures import ProcessPoolExecutor from concurrent.futures.process import BrokenProcessPool +from copy import deepcopy from pathlib import Path from typing import List, Optional, Dict, Sequence, Tuple @@ -26,6 +27,7 @@ def __init__( mech_directory: Path, parameters: List[CoefIndex], aliases: Optional[Dict[str, str]] = None, + num_workers: Optional[int] = None ): """ @@ -34,10 +36,11 @@ def __init__( mech_directory: Path to the mechanism file(s) parameters: Parameters to be adjusted aliases: Aliases that map species in the experiment to the mechanism description + num_workers: Maximum number of parallel workers to allow """ self.parameters = parameters.copy() - # Make a placeholder for the Frhodo instance + # Store where to find the experimental data self._exp_directory = exp_directory self._mech_directory = mech_directory self._aliases = aliases @@ -50,6 +53,20 @@ def __init__( self._sim_kwargs: Optional[List[dict]] = None self._rxn_conditions: Optional[List[Tuple[float, float, dict]]] = None + # Holder for the process pool, which is not serialized + self._exec: Optional[ProcessPoolExecutor] = ProcessPoolExecutor(num_workers) + self._num_workers = num_workers + + def __getstate__(self): + state = self.__dict__.copy() + del state['_exec'] + return state + + def __setstate__(self, state: dict): + state = state.copy() + state['_exec'] = ProcessPoolExecutor(state['_num_workers']) + self.__dict__.update(state) + @property def observations(self) -> List[np.ndarray]: """Experimental data less any regions that are masked out""" @@ -113,21 +130,22 @@ def run_simulations(self, x: Sequence[float]) -> List[np.ndarray]: Simulated for each experiment interpolated at the same time increments as :meth:`observations`. """ + # Make a copy of the mech so we're threadsafe + mech = deepcopy(self.mech) # Update the parameters assert len(x) == len(self.parameters), f"Expected {len(self.parameters)} parameters but got {len(x)}" - set_coefficients(self.mech, dict(zip(self.parameters, x))) + set_coefficients(mech, dict(zip(self.parameters, x))) # Run each simulation to each set of experimental data sims = [] - with ProcessPoolExecutor(1) as exec: - for sim_kwargs, rxn_cond in zip(self._sim_kwargs, self._rxn_conditions): - exc = exec.submit(run_simulation, self.mech, rxn_cond, sim_kwargs) - try: - res = exc.result() - except BrokenProcessPool: - raise ValueError('Running the simulation failed') - sims.append(res) + for sim_kwargs, rxn_cond in zip(self._sim_kwargs, self._rxn_conditions): + future = self._exec.submit(run_simulation, mech, rxn_cond, sim_kwargs) + try: + res = future.result() + except BrokenProcessPool: + raise ValueError('Running the simulation failed') + sims.append(res) # Interpolate simulation data over the same steps as the experiments output = [] @@ -186,6 +204,7 @@ def __init__( parameters: List[CoefIndex], priors: Optional[List[ss.rv_continuous]] = None, aliases: Optional[Dict[str, str]] = None, + **kwargs ): """ @@ -195,7 +214,7 @@ def __init__( parameters: Parameters to be adjusted aliases: Aliases that map species in the experiment to the mechanism description """ - super().__init__(exp_directory, mech_directory, parameters, aliases) + super().__init__(exp_directory, mech_directory, parameters, aliases, **kwargs) self.priors = priors def compute_log_probs(self, x: np.ndarray, uncertainty: float) -> np.ndarray: From 66ef87d5a911bdf80797da06a12d3af2161b6533 Mon Sep 17 00:00:00 2001 From: Logan Ward Date: Mon, 26 Sep 2022 11:35:32 -0400 Subject: [PATCH 27/29] One process per parameters, not per simulation Speeds us up significantly --- frhodo/api/optimize.py | 68 ++++++++++++++++++++++++++---------------- 1 file changed, 42 insertions(+), 26 deletions(-) diff --git a/frhodo/api/optimize.py b/frhodo/api/optimize.py index d38592a..2931210 100644 --- a/frhodo/api/optimize.py +++ b/frhodo/api/optimize.py @@ -1,9 +1,6 @@ """Utilities useful for using Frhodo from external optimizers""" -import multiprocessing -import warnings from concurrent.futures import ProcessPoolExecutor from concurrent.futures.process import BrokenProcessPool -from copy import deepcopy from pathlib import Path from typing import List, Optional, Dict, Sequence, Tuple @@ -15,6 +12,26 @@ from frhodo.calculate.mech_fcns import Chemical_Mechanism +def _run_simulation(x, mech, parameters, observations, rxn_conditions, sim_kwargs): + """Private method to run the simulations. Intended to be run in a Process as the simulator can throw seg faults""" + + # Update the parameters + set_coefficients(mech, dict(zip(parameters, x))) + + # Run each simulation to each set of experimental data + sims = [] + for sim_kwargs, rxn_cond in zip(sim_kwargs, rxn_conditions): + res = run_simulation(mech, rxn_cond, sim_kwargs) + sims.append(res) + + # Interpolate simulation data over the same steps as the experiments + output = [] + for sim, obs in zip(sims, observations): + sim_func = interp1d(sim[:, 0], sim[:, 1], kind='cubic', fill_value='extrapolate') + output.append(sim_func(obs[:, 0])) + return output + + class BaseObjectiveFunction: """Base class for a Frhodo-based objective function @@ -82,6 +99,19 @@ def weights(self) -> List[np.ndarray]: self.load_experiments() return self._weights + @property + def sim_kwargs(self) -> List[dict]: + """Keyword arguments used to simulate each experiment""" + if self._sim_kwargs is None: + self.load_experiments() + return self._sim_kwargs.copy() + + def rxn_conditions(self) -> List[Tuple[float, float, dict]]: + """Starting conditions for each experiment""" + if self._rxn_conditions is None: + self.load_experiments() + return self._rxn_conditions.copy() + @property def x(self) -> np.ndarray: """Get the current state of the coefficient being optimized @@ -130,29 +160,15 @@ def run_simulations(self, x: Sequence[float]) -> List[np.ndarray]: Simulated for each experiment interpolated at the same time increments as :meth:`observations`. """ - # Make a copy of the mech so we're threadsafe - mech = deepcopy(self.mech) - - # Update the parameters - assert len(x) == len(self.parameters), f"Expected {len(self.parameters)} parameters but got {len(x)}" - set_coefficients(mech, dict(zip(self.parameters, x))) - - # Run each simulation to each set of experimental data - sims = [] - for sim_kwargs, rxn_cond in zip(self._sim_kwargs, self._rxn_conditions): - future = self._exec.submit(run_simulation, mech, rxn_cond, sim_kwargs) - try: - res = future.result() - except BrokenProcessPool: - raise ValueError('Running the simulation failed') - sims.append(res) - - # Interpolate simulation data over the same steps as the experiments - output = [] - for sim, obs in zip(sims, self.observations): - sim_func = interp1d(sim[:, 0], sim[:, 1], kind='cubic', fill_value='extrapolate') - output.append(sim_func(obs[:, 0])) - return output + + # Submit the simulation as a subprocess + future = self._exec.submit(_run_simulation, + x, self.mech, self.parameters, self.observations, + self._rxn_conditions, self._sim_kwargs) + try: + return future.result() + except BrokenProcessPool: + raise ValueError(f'Process pool failed for: {x}') def compute_residuals(self, x: Sequence[float]) -> List[np.ndarray]: """Compute the residual between simulation and experiment From 40184a065018adb2d5d007aab1e436c9cc33475f Mon Sep 17 00:00:00 2001 From: Logan Ward Date: Mon, 26 Sep 2022 15:54:44 -0400 Subject: [PATCH 28/29] Split off interface functions, add one for k's --- frhodo/api/driver.py | 63 +------------------------------ frhodo/api/interface.py | 83 +++++++++++++++++++++++++++++++++++++++++ frhodo/api/optimize.py | 6 ++- tests/test_api.py | 6 +++ 4 files changed, 95 insertions(+), 63 deletions(-) create mode 100644 frhodo/api/interface.py diff --git a/frhodo/api/driver.py b/frhodo/api/driver.py index a5776fe..39aecf8 100644 --- a/frhodo/api/driver.py +++ b/frhodo/api/driver.py @@ -7,68 +7,9 @@ import numpy as np from PyQt5.QtWidgets import QApplication -from frhodo.calculate.mech_fcns import Chemical_Mechanism -from frhodo.calculate.reactors import Simulation_Result +from frhodo.api.interface import set_coefficients, run_simulation, get_coefficients, CoefIndex from frhodo.main import Main, launch_gui -CoefIndex = Tuple[int, Union[str, int], str] -"""How to specify the index of a reaction parameter: reaction index, pressure index, name of the parameter""" - - -def set_coefficients(mech: Chemical_Mechanism, new_values: Dict[CoefIndex, float]): - """Get the coefficients of a chemical mechanism object - - Args: - mech: Mechanism to modify - new_values: New values for different coefficients - """ - for key, new_value in new_values.items(): - rxn_id, prs_id, coef_name = key - rxn_model = mech.coeffs[rxn_id][prs_id] - assert coef_name in rxn_model, f'Key {coef_name} not present for reaction #{rxn_id} at pressure #{prs_id}' - rxn_model[coef_name] = new_value - - # Update the reaction rates for the current shock - mech.modify_reactions(mech.coeffs) - - -def run_simulation(mech: Chemical_Mechanism, rxn_conditions: Tuple[float, float, dict], sim_kwargs: dict) -> np.ndarray: - """Run a simulation for a single reaction condition - - Args: - mech: Mechanism describing the chemical kinetics - rxn_conditions: Conditions of the reaction (temperature, pressure, composition) - sim_kwargs: Keywords to the simulator - Returns: - Array with time as first column and observable as the second - """ - mech.set_TPX(*rxn_conditions) - # Run the simulation - # TODO (wardlt): Do not hardcode the runtime or reactor conditions - sim, _ = mech.run('Incident Shock Reactor', 1.2e-05, *rxn_conditions, **sim_kwargs) - assert sim.success, "Simulation failed" - return np.stack([ - sim.independent_var, - sim.observable - ], axis=1) - - -def get_coefficients(mech: Chemical_Mechanism, indices: List[CoefIndex]) -> List[float]: - """Get specific coefficients from the reaction model - - Args: - mech: Mechanism to interrogate - indices: List of coefficients to retrieve - Returns: - Desired coefficents - """ - output = [] - for rxn_id, prs_id, coef_name in indices: - rxn_model = mech.coeffs[rxn_id][prs_id] - assert coef_name in rxn_model, f'Key {coef_name} not present for reaction #{rxn_id} at pressure #{prs_id}' - output.append(rxn_model[coef_name]) - return output - class FrhodoDriver: """Driver the Frhodo GUI application @@ -158,7 +99,7 @@ def get_simulator_inputs(self) -> Tuple[List[dict], List[Tuple[float, float, Dic Returns: - Keyword arguments passed to the simulator - - Temperature, pressure and concentration of different guasses + - Temperature, pressure and concentration of different experiments """ var_outputs = [] rxn_outputs = [] diff --git a/frhodo/api/interface.py b/frhodo/api/interface.py new file mode 100644 index 0000000..c5350bc --- /dev/null +++ b/frhodo/api/interface.py @@ -0,0 +1,83 @@ +"""Simple functions for common tasks in Frhodo. + +The goal is to eventually distribute alongside the GUI components""" +from typing import Dict, Tuple, List, Union + +import numpy as np + +from frhodo.calculate.mech_fcns import Chemical_Mechanism + +CoefIndex = Tuple[int, Union[str, int], str] +"""How to specify the index of a reaction parameter: reaction index, pressure index, name of the parameter""" + +RxnConditions = Tuple[float, float, dict] +"""Specification for the temperature, pressure and concentration for a reaction""" + + +def set_coefficients(mech: Chemical_Mechanism, new_values: Dict[CoefIndex, float]): + """Get the coefficients of a chemical mechanism object + + Args: + mech: Mechanism to modify + new_values: New values for different coefficients + """ + for key, new_value in new_values.items(): + rxn_id, prs_id, coef_name = key + rxn_model = mech.coeffs[rxn_id][prs_id] + assert coef_name in rxn_model, f'Key {coef_name} not present for reaction #{rxn_id} at pressure #{prs_id}' + rxn_model[coef_name] = new_value + + # Update the reaction rates for the current shock + mech.modify_reactions(mech.coeffs) + + +def run_simulation(mech: Chemical_Mechanism, rxn_conditions: RxnConditions, sim_kwargs: dict) -> np.ndarray: + """Run a simulation for a single reaction condition + + Args: + mech: Mechanism describing the chemical kinetics + rxn_conditions: Conditions of the reaction (temperature, pressure, composition) + sim_kwargs: Keywords to the simulator + Returns: + Array with time as first column and observable as the second + """ + mech.set_TPX(*rxn_conditions) + # Run the simulation + # TODO (wardlt): Do not hardcode the runtime or reactor conditions + sim, _ = mech.run('Incident Shock Reactor', 1.2e-05, *rxn_conditions, **sim_kwargs) + assert sim.success, "Simulation failed" + return np.stack([ + sim.independent_var, + sim.observable + ], axis=1) + + +def get_coefficients(mech: Chemical_Mechanism, indices: List[CoefIndex]) -> List[float]: + """Get specific coefficients from the reaction model + + Args: + mech: Mechanism to interrogate + indices: List of coefficients to retrieve + Returns: + Desired coefficents + """ + output = [] + for rxn_id, prs_id, coef_name in indices: + rxn_model = mech.coeffs[rxn_id][prs_id] + assert coef_name in rxn_model, f'Key {coef_name} not present for reaction #{rxn_id} at pressure #{prs_id}' + output.append(rxn_model[coef_name]) + return output + + +def compute_kinetic_coefficients(mech: Chemical_Mechanism, rxn_conditions: RxnConditions) -> np.ndarray: + """Get the kinetic coefficients at a certain conditions + + Args: + mech: Mechanism object to use + rxn_conditions: Reaction conditions (T, P, X) + """ + + mech.set_TPX(*rxn_conditions) + return np.array([ + mech.gas.forward_rate_constants[i] for i in range(mech.gas.n_reactions) + ]) diff --git a/frhodo/api/optimize.py b/frhodo/api/optimize.py index 2931210..30e9682 100644 --- a/frhodo/api/optimize.py +++ b/frhodo/api/optimize.py @@ -8,7 +8,8 @@ from scipy import stats as ss from scipy.interpolate import interp1d -from frhodo.api.driver import CoefIndex, FrhodoDriver, set_coefficients, get_coefficients, run_simulation +from frhodo.api.driver import FrhodoDriver +from frhodo.api.interface import set_coefficients, run_simulation, get_coefficients, CoefIndex from frhodo.calculate.mech_fcns import Chemical_Mechanism @@ -31,7 +32,7 @@ def _run_simulation(x, mech, parameters, observations, rxn_conditions, sim_kwarg output.append(sim_func(obs[:, 0])) return output - +# TODO (wardlt): Add in the class BaseObjectiveFunction: """Base class for a Frhodo-based objective function @@ -106,6 +107,7 @@ def sim_kwargs(self) -> List[dict]: self.load_experiments() return self._sim_kwargs.copy() + @property def rxn_conditions(self) -> List[Tuple[float, float, dict]]: """Starting conditions for each experiment""" if self._rxn_conditions is None: diff --git a/tests/test_api.py b/tests/test_api.py index 8d19c64..7d298c9 100644 --- a/tests/test_api.py +++ b/tests/test_api.py @@ -5,6 +5,7 @@ from multiprocessing import set_start_method from frhodo.api.driver import FrhodoDriver +from frhodo.api.interface import compute_kinetic_coefficients from frhodo.api.optimize import BayesianObjectiveFunction @@ -67,6 +68,11 @@ def test_update(loaded_frhodo, rxn_id, prs_id): rates = loaded_frhodo.get_reaction_rates() assert rates.shape == (66, 1) + # Get them at a specific TPX + rxn_conditions = loaded_frhodo.get_simulator_inputs()[1][0] + rates_single = compute_kinetic_coefficients(loaded_frhodo.window.mech, rxn_conditions) + assert np.isclose(rates_single, rates[:, 0]).all() + # Get the simulation before sim = loaded_frhodo.run_simulations()[0] From 11e05f6e491df202061102996eed002484bc6a3e Mon Sep 17 00:00:00 2001 From: Logan Ward Date: Thu, 29 Sep 2022 17:35:56 -0400 Subject: [PATCH 29/29] Add ability to optimize log coefficients --- frhodo/api/optimize.py | 31 +++++++++++++++++++++---------- tests/test_api.py | 15 ++++++++++++--- 2 files changed, 33 insertions(+), 13 deletions(-) diff --git a/frhodo/api/optimize.py b/frhodo/api/optimize.py index 30e9682..1f64d2a 100644 --- a/frhodo/api/optimize.py +++ b/frhodo/api/optimize.py @@ -32,7 +32,7 @@ def _run_simulation(x, mech, parameters, observations, rxn_conditions, sim_kwarg output.append(sim_func(obs[:, 0])) return output -# TODO (wardlt): Add in the + class BaseObjectiveFunction: """Base class for a Frhodo-based objective function @@ -44,6 +44,7 @@ def __init__( exp_directory: Path, mech_directory: Path, parameters: List[CoefIndex], + parameter_is_log: Optional[List[bool]] = None, aliases: Optional[Dict[str, str]] = None, num_workers: Optional[int] = None ): @@ -52,7 +53,8 @@ def __init__( Args: exp_directory: Path to the experimental data mech_directory: Path to the mechanism file(s) - parameters: Parameters to be adjusted + parameters: Kinetic coefficients to be adjusted + parameter_is_log: Whether the input coefficient is a log value. Optional: All false aliases: Aliases that map species in the experiment to the mechanism description num_workers: Maximum number of parallel workers to allow """ @@ -64,6 +66,11 @@ def __init__( self._aliases = aliases self._frhodo: Optional[FrhodoDriver] = None + # Determine whether the inputs are logs + if parameter_is_log is None: + parameter_is_log = [False] * len(parameters) + self.parameter_is_log = parameter_is_log.copy() + # Make a placeholder for experimental data self.mech: Optional[Chemical_Mechanism] = None self._observations: Optional[List[np.ndarray]] = None @@ -114,13 +121,12 @@ def rxn_conditions(self) -> List[Tuple[float, float, dict]]: self.load_experiments() return self._rxn_conditions.copy() - @property - def x(self) -> np.ndarray: - """Get the current state of the coefficient being optimized - - This is either the initial values from loading the mechanism or the last values ran. - """ - return np.array(get_coefficients(self.mech, self.parameters)) + def get_initial_values(self) -> np.ndarray: + """Get the initial values from loading the mechanism""" + return np.array([ + np.log(x) if l else x for x, l in + zip(get_coefficients(self.mech, self.parameters), + self.parameter_is_log)]) def load_experiments(self, frhodo: Optional[FrhodoDriver] = None): """Load observations and weights from disk @@ -163,6 +169,9 @@ def run_simulations(self, x: Sequence[float]) -> List[np.ndarray]: as :meth:`observations`. """ + # Convert the inputs to real values, if needed + x = [np.exp(xi) if li else xi for xi, li in zip(x, self.parameter_is_log)] + # Submit the simulation as a subprocess future = self._exec.submit(_run_simulation, x, self.mech, self.parameters, self.observations, @@ -220,6 +229,7 @@ def __init__( exp_directory: Path, mech_directory: Path, parameters: List[CoefIndex], + parameter_is_log: Optional[List[bool]] = None, priors: Optional[List[ss.rv_continuous]] = None, aliases: Optional[Dict[str, str]] = None, **kwargs @@ -230,9 +240,10 @@ def __init__( exp_directory: Path to the experimental data mech_directory: Path to the mechanism file(s) parameters: Parameters to be adjusted + parameter_is_log: Whether the input coefficient is a log value. Optional: All false aliases: Aliases that map species in the experiment to the mechanism description """ - super().__init__(exp_directory, mech_directory, parameters, aliases, **kwargs) + super().__init__(exp_directory, mech_directory, parameters, parameter_is_log, aliases, **kwargs) self.priors = priors def compute_log_probs(self, x: np.ndarray, uncertainty: float) -> np.ndarray: diff --git a/tests/test_api.py b/tests/test_api.py index 7d298c9..5400f2a 100644 --- a/tests/test_api.py +++ b/tests/test_api.py @@ -130,10 +130,10 @@ def test_optimizer(loaded_frhodo, example_dir, tmp_path): # Test the state assert opt.weights[0].max() == 1 - assert len(opt.x) == 1 + assert len(opt.get_initial_values()) == 1 # Make sure the optimizer produces different results with different inputs - x0 = opt.x.tolist() + x0 = opt.get_initial_values().tolist() x0.insert(0, 1e-4) y0 = opt(x0) @@ -147,6 +147,15 @@ def test_optimizer(loaded_frhodo, example_dir, tmp_path): assert y0 == y0_repeat # Make sure it is serializable - opt2 = pkl.loads(pkl.dumps(opt)) + opt2: BayesianObjectiveFunction = pkl.loads(pkl.dumps(opt)) y0_opt2 = opt2(x0) assert np.isclose(y0, y0_opt2).all() + + # Create another version where we test the log of the initial parameter + assert opt2.parameter_is_log == [False] + opt2.parameter_is_log = [True] + assert np.isclose(opt.get_initial_values(), np.exp(opt2.get_initial_values())).all() + + x2 = [1e-4] + np.log(x0[1:]).tolist() + y2 = opt2(x2) + assert np.isclose(y0, y2).all()