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/frhodo/api/__init__.py b/frhodo/api/__init__.py new file mode 100644 index 0000000..50785c0 --- /dev/null +++ b/frhodo/api/__init__.py @@ -0,0 +1,2 @@ +"""Functions related to using Frhodo from a Python interface""" + diff --git a/frhodo/api/driver.py b/frhodo/api/driver.py new file mode 100644 index 0000000..39aecf8 --- /dev/null +++ b/frhodo/api/driver.py @@ -0,0 +1,269 @@ +"""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.api.interface import set_coefficients, run_simulation, get_coefficients, CoefIndex +from frhodo.main import Main, launch_gui + + +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 + + - 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 + """ + + 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_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 experiments + """ + 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 + + 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 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 + 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 + + 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 + """ + + mech = self.window.mech + output = get_coefficients(mech, indices) + return output + + def set_coefficients(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 + mech = self.window.mech + set_coefficients(mech, new_values) + + @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/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 new file mode 100644 index 0000000..1f64d2a --- /dev/null +++ b/frhodo/api/optimize.py @@ -0,0 +1,289 @@ +"""Utilities useful for using Frhodo from external optimizers""" +from concurrent.futures import ProcessPoolExecutor +from concurrent.futures.process import BrokenProcessPool +from pathlib import Path +from typing import List, Optional, Dict, Sequence, Tuple + +import numpy as np +from scipy import stats as ss +from scipy.interpolate import interp1d + +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 + + +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 + + 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], + parameter_is_log: Optional[List[bool]] = None, + aliases: Optional[Dict[str, str]] = None, + num_workers: Optional[int] = None + ): + """ + + Args: + exp_directory: Path to the experimental data + mech_directory: Path to the mechanism file(s) + 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 + """ + self.parameters = parameters.copy() + + # Store where to find the experimental data + self._exp_directory = exp_directory + self._mech_directory = mech_directory + 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 + self._weights: Optional[List[np.ndarray]] = None + 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""" + # 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 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() + + @property + 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() + + 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 + + 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(*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]) + + 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 + + Args: + x: New set of reaction coefficients + Returns: + Simulated for each experiment interpolated at the same time increments + 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, + 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 + + 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() + + 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 + + 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 __init__( + self, + 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 + ): + """ + + Args: + 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, parameter_is_log, aliases, **kwargs) + self.priors = priors + + 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, frhodo: Optional[FrhodoDriver] = None): + super().load_experiments(frhodo) + + # 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" + + # Compute the log probability of the data + uncertainty, coeffs = x[0], x[1:] + 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/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 old mode 100644 new mode 100755 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 67% rename from src/calculate/mech_fcns.py rename to frhodo/calculate/mech_fcns.py index 65e731b..0b4c438 100644 --- a/src/calculate/mech_fcns.py +++ b/frhodo/calculate/mech_fcns.py @@ -4,48 +4,87 @@ 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 calculate import reactors, shock_fcns, integrate -import ck2yaml from timeit import default_timer as timer +from . import reactors, shock_fcns, integrate +from .reactors import Simulation_Result +from .. import ck2yaml, soln2ck + class Chemical_Mechanism: 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=True, + 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: @@ -54,9 +93,9 @@ def loader(self, path): except RuntimeError as e: print('FAILED.') 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): @@ -69,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'] @@ -109,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)): @@ -168,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 @@ -194,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: @@ -209,91 +250,106 @@ 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""" if type(rxn) in [ct.ElementaryReaction, ct.ThreeBodyReaction]: bnds_key = 'rate' coef_key = 0 @@ -304,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' @@ -319,9 +375,12 @@ 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 + 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) @@ -329,17 +388,25 @@ 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,} - - self.thermo_coeffs.append(thermo_dict) - - def modify_reactions(self, coeffs, rxnIdxs=[]): # Only works for Arrhenius equations currently - if not rxnIdxs: # if rxnNums does not exist, modify all + 'h_scaler': 1, 's_scaler': 1, } + + 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: rxn = self.gas.reaction(rxnIdx) rxnChanged = False @@ -347,8 +414,8 @@ def modify_reactions(self, coeffs, rxnIdxs=[]): # Only works for Arrhenius e 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'] @@ -376,36 +443,36 @@ def modify_reactions(self, coeffs, rxnIdxs=[]): # Only works for Arrhenius e 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) @@ -424,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) @@ -438,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: @@ -452,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] @@ -469,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) @@ -490,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 @@ -517,41 +585,73 @@ 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 - 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[Simulation_Result, 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 +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': - 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 == '+': 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': - #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'] @@ -562,11 +662,11 @@ 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] 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/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 99% rename from src/calculate/optimize/CheKiPEUQ_local/plotting_functions.py rename to frhodo/calculate/optimize/CheKiPEUQ_local/plotting_functions.py index 17e2af3..c0b3899 100644 --- a/src/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/src/calculate/optimize/CheKiPEUQ_local/test_pytest.py b/frhodo/calculate/optimize/CheKiPEUQ_local/test_pytest.py similarity index 81% rename from src/calculate/optimize/CheKiPEUQ_local/test_pytest.py rename to frhodo/calculate/optimize/CheKiPEUQ_local/test_pytest.py index 05b3f6c..ddcdcb8 100644 --- a/src/calculate/optimize/CheKiPEUQ_local/test_pytest.py +++ b/frhodo/calculate/optimize/CheKiPEUQ_local/test_pytest.py @@ -1,16 +1,19 @@ # 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 +from pytest import mark + + +@mark.xfail 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 +@mark.xfail def MH_then_scatterplot(): parseUserInputParameters() ip_object = ip() 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 90% rename from src/calculate/optimize/fit_fcn.py rename to frhodo/calculate/optimize/fit_fcn.py index e32db8d..5fec164 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): @@ -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/src/calculate/optimize/mech_optimize.py b/frhodo/calculate/optimize/mech_optimize.py similarity index 96% rename from src/calculate/optimize/mech_optimize.py rename to frhodo/calculate/optimize/mech_optimize.py index f15a6d0..fb788f8 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'] @@ -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/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 88% rename from src/calculate/optimize/optimize_worker.py rename to frhodo/calculate/optimize/optimize_worker.py index 5fb4019..145f4ef 100644 --- a/src/calculate/optimize/optimize_worker.py +++ b/frhodo/calculate/optimize/optimize_worker.py @@ -12,8 +12,10 @@ 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 +from frhodo.calculate.mech_fcns import Chemical_Mechanism +import frhodo debug = True @@ -34,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 @@ -56,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 @@ -75,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, @@ -86,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') @@ -210,7 +229,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' @@ -224,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 @@ -294,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: @@ -373,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, @@ -413,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/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..07c79c4 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 @@ -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/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 98% rename from src/config_io.py rename to frhodo/config_io.py index 070b38e..449a70e 100644 --- a/src/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 @@ -202,6 +207,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 +324,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 +418,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/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 old mode 100644 new mode 100755 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 72% rename from src/main.py rename to frhodo/main.py index e3712c9..3e0030b 100644 --- a/src/main.py +++ b/frhodo/main.py @@ -5,25 +5,24 @@ # 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, Optional, List 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) - + # Handle high resolution displays: Minimum recommended resolution 1280 x 960 if hasattr(QtCore.Qt, 'AA_EnableHighDpiScaling'): QApplication.setAttribute(QtCore.Qt.AA_EnableHighDpiScaling, True) @@ -31,44 +30,60 @@ 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) 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): + 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 - 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) - + 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 @@ -80,38 +95,48 @@ 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 + self.initialize_settings(skip_config) # ~ 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 + + 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() 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) + 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? - + # 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: @@ -120,21 +145,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: @@ -144,7 +169,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 @@ -157,7 +182,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: @@ -167,22 +192,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': @@ -191,34 +218,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: @@ -236,45 +263,42 @@ 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'] - 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'], - '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'], + + # Make sure the rate constants are update-to-date with the conditions on the table and specified shock + self.tree.update_rates() + + # 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) - + 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': @@ -287,20 +311,60 @@ 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 - -if __name__ == '__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 + """ + # 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 + 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) + # 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, my_path, shut_down) + + # Create the Frhodo main window + main_window = Main(app, my_path, fresh_gui) + 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 97% rename from src/mech_widget.py rename to frhodo/mech_widget.py index 1963f32..8c46819 100644 --- a/src/mech_widget.py +++ b/frhodo/mech_widget.py @@ -2,17 +2,18 @@ # 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 +from .calculate.mech_fcns import Chemical_Mechanism + def silentSetValue(obj, value): obj.blockSignals(True) # stop changing text from signaling @@ -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' @@ -283,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 @@ -367,6 +380,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 +416,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 +429,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 +485,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/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 98% rename from src/options_panel_widgets.py rename to frhodo/options_panel_widgets.py index a88b59e..20b610a 100644 --- a/src/options_panel_widgets.py +++ b/frhodo/options_panel_widgets.py @@ -3,17 +3,21 @@ # directory for license and copyright information. import pathlib, os, sys +from copy import deepcopy +from typing import Union, List + 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): @@ -128,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): @@ -141,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' @@ -150,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 @@ -164,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)') @@ -186,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 @@ -269,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: @@ -558,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) @@ -627,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() @@ -970,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() @@ -1022,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/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 95% rename from src/settings.py rename to frhodo/settings.py index 471727c..0b453c6 100644 --- a/src/settings.py +++ b/frhodo/settings.py @@ -10,24 +10,27 @@ 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) + class Path: + """Facilitates loading data from disk""" + def __init__(self, parent, path): self.parent = parent 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): @@ -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) @@ -337,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'])) @@ -619,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) @@ -659,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 @@ -728,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) @@ -746,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 @@ -988,9 +1006,18 @@ 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 a single shock experiment + + Modifies the `mech` object of the Frhodo main window + + 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 +1027,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 @@ -1013,7 +1040,12 @@ def rates(self, shock, rxnIdx=[]): # This resets and updates all rates in shock 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 @@ -1050,6 +1082,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/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..087b58e --- /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/frhodo_environment.yml b/frhodo_environment.yml index bcb685d..06c53e1 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 @@ -20,6 +19,8 @@ dependencies: - requests=2.27.1 - scipy=1.7.3 - tabulate=0.8.9 + - pytest - pip: - dtcwt==0.12.0 - - rbfopt==4.2.2 \ No newline at end of file + - rbfopt==4.2.2 + - -e . diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..fd3d99d --- /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=version, + 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' + ] + } +) diff --git a/tests/conftest.py b/tests/conftest.py new file mode 100644 index 0000000..efb6070 --- /dev/null +++ b/tests/conftest.py @@ -0,0 +1,21 @@ +from pathlib import Path + +from pytest import fixture + +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 +app, window = launch_gui(['frhodo', '-platform', 'offscreen'], fresh_gui=True) +driver = FrhodoDriver(window, app) + + +@fixture() +def example_dir() -> Path: + """Directory containing example input files""" + return Path(__file__).parent.parent / 'Example' + + +@fixture() +def frhodo_driver() -> FrhodoDriver: + return driver diff --git a/tests/test_api.py b/tests/test_api.py new file mode 100644 index 0000000..5400f2a --- /dev/null +++ b/tests/test_api.py @@ -0,0 +1,161 @@ +"""Testing the API components of Frhodo""" +import numpy as np +import pickle as pkl +from pytest import fixture, mark +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 + + +@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', + example_dir / 'Mechanism', + tmp_path, + {'A': 'B'} # A fake alias + ) + return frhodo_driver + + +def test_launch(): + """Make sure we can launch Frhodo""" + driver = FrhodoDriver.create_driver() + assert driver.window.isVisible() + + +def test_load(loaded_frhodo): + """Load loading in desired data""" + assert loaded_frhodo.n_shocks == 1 + + +def test_observables(loaded_frhodo): + 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() + 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 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 + (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) + + # 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] + + # 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.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() + 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] + 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): + """Test getting lists of parameters to update""" + + # Make sure we can get all parameters + total = loaded_frhodo.get_fittable_parameters() + assert len(total) == 360 + + # Make sure it works with each reaction type + 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 + + # 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): + 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) + opt.load_experiments(loaded_frhodo) + + # Test the state + assert opt.weights[0].max() == 1 + assert len(opt.get_initial_values()) == 1 + + # Make sure the optimizer produces different results with different inputs + x0 = opt.get_initial_values().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 it is serializable + 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()