Changes motivated by pharmapy-opti work

InstallOnMac.sh

@@ -5,9 +5,13 @@ echo --------------------------------------------------
 echo Enter environment name:
 read env_name 
 echo ------------------------------
-echo Creating pandas environment...
+echo Creating conda environment...
 echo ------------------------------
 conda create -n $env_name python=3.9 --file requirements.txt -c conda-forge
+echo ------------------------------
+echo Activating conda environment...
+echo ------------------------------
+conda init
 conda activate $env_name
 echo ----------------------
 echo Installing PharmaPy...

PharmaPy/Kinetics.py

@@ -239,6 +239,8 @@ def __init__(self, path, k_params, ea_params, rxn_list=None,
             stoich_matrix = get_stoich(di, partic_species)
         else:
             stoich_matrix = np.atleast_2d(stoich_matrix)
+            if partic_species is None:
+                raise PharmaPyTypeError('Please provide a participating species list when using a stoichiometric matrix.')
 
         perm_idx = get_permutation_indexes(name_species, partic_species)
         stoich_matrix = stoich_matrix[:, perm_idx]

PharmaPy/ParamEstim.py

@@ -22,7 +22,11 @@
 from PharmaPy.LevMarq import levenberg_marquardt
 from PharmaPy.Commons import plot_sens, reorder_sens
 
-from cyipopt import minimize_ipopt
+try:
+    from cyipopt import minimize_ipopt
+    have_cyipopt = True
+except ImportError:
+    have_cyipopt = False
 
 linestyles = cycle(['-', '--', '-.', ':'])
 
@@ -615,6 +619,8 @@ def optimize_fn(self, optim_options=None, simulate=False, verbose=True,
                 **optim_options)
 
         elif method == 'IPOPT':
+            if not have_cyipopt:
+                raise ImportError('cyipopt is an optional import. Please install cyipopt to use IPOPT as a solver for parameter estimation. conda install -c conda-forge cyipopt')
             if optim_options is None:
                 optim_options = {'print_level': int(verbose) * 5}
             else:

PharmaPy/Reactors.py

@@ -19,7 +19,7 @@
 from PharmaPy.CheckModule import check_modeling_objects
 
 import numpy as np
-from numpy.core.umath_tests import inner1d
+# from numpy.core.umath_tests import inner1d
 import matplotlib.pyplot as plt
 from matplotlib.ticker import AutoMinorLocator
 from matplotlib.animation import FuncAnimation
@@ -125,10 +125,10 @@ class _BaseReactor:
         f(time) = state_value
     return_sens : bool (optional, default = True)
         whether or not the paramest_wrapper method should return
-        the sensitivity system along with the concentratio profiles.
+        the sensitivity system along with the concentration profiles.
         Use False if you want the parameter estimation platform to
         estimate the sensitivity system using finite differences
-    state_events : lsit of dict(s)
+    state_events : list of dict(s)
         list of dictionaries, each one containing the specification of a
         state event
     """
@@ -506,7 +506,7 @@ def plot_profiles(self, pick_comp=None, **fig_kwargs):
             can be either the name of a species (str) or the index of the
             species (int). The default is None.
         **fig_kwargs : keyword arguments to plt.subplots()
-            named arguments passed to the plotting functions. A yypical field
+            named arguments passed to the plotting functions. A typical field
             is 'figsize', passed as a (width, height) tuple.
 
         Returns
@@ -605,7 +605,7 @@ class BatchReactor(_BaseReactor):
         'coil', 'bath']
     return_sens : bool (optional, default = True)
         whether or not the paramest_wrapper method should return
-        the sensitivity system along with the concentratio profiles.
+        the sensitivity system along with the concentration profiles.
         Use False if you want the parameter estimation platform to
         estimate the sensitivity system using finite differences
     """
@@ -688,7 +688,10 @@ def energy_balances(self, time, mole_conc, vol, temp, temp_ht, inputs,
                                             delta_hrxn=deltah_rxn)
 
         # Balance terms (W)
-        source_term = -inner1d(deltah_rxn, rates) * vol*1000  # vol in L
+        #source_term = -inner1d(deltah_rxn, rates) * vol * 1000  # vol in L
+        # TODO: Check if this is correct
+        # source_term = -np.inner(deltah_rxn, rates) * vol * 1000  # vol in L
+        source_term = -(deltah_rxn * rates).sum(axis=1) * vol * 1000  # vol in L
 
         if heat_prof:
             if 'temp' in self.controls.keys():
@@ -942,7 +945,7 @@ class CSTR(_BaseReactor):
         'coil', 'bath']
     return_sens : bool (optional, default = True)
         whether or not the paramest_wrapper method should return
-        the sensitivity system along with the concentratio profiles.
+        the sensitivity system along with the concentration profiles.
         Use False if you want the parameter estimation platform to
         estimate the sensitivity system using finite differences
     """
@@ -1045,7 +1048,10 @@ def energy_balances(self, time, mole_conc, vol, temp, temp_ht, inputs,
         flow_term = inlet_flow * (h_in - h_temp)  # W
 
         # Balance terms (W) - convert vol to L
-        source_term = -inner1d(deltah_rxn, rates) * vol * 1000
+        # source_term = -inner1d(deltah_rxn, rates) * vol * 1000
+        # TODO: Check if this is correct
+        # source_term = -np.dot(deltah_rxn, rates) * vol * 1000  # vol in L
+        source_term = -(deltah_rxn * rates).sum(axis=1) * vol * 1000  # vol in L
 
         if heat_prof:
             if self.isothermal:
@@ -1242,7 +1248,7 @@ class SemibatchReactor(CSTR):
         'coil', 'bath']
     return_sens : bool (optional, default = True)
         whether or not the paramest_wrapper method should return
-        the sensitivity system along with the concentratio profiles.
+        the sensitivity system along with the concentration profiles.
         Use False if you want the parameter estimation platform to
         estimate the sensitivity system using finite differences
     """
@@ -1442,7 +1448,7 @@ class PlugFlowReactor(_BaseReactor):
         'coil', 'bath']
     return_sens : bool (optional, default = True)
         whether or not the paramest_wrapper method should return
-        the sensitivity system along with the concentratio profiles.
+        the sensitivity system along with the concentration profiles.
         Use False if you want the parameter estimation platform to
         estimate the sensitivity system using finite differences
     """
@@ -1550,7 +1556,10 @@ def energy_steady(self, conc, temp):
                                             delta_hrxn=deltah_rxn)
 
         # ---------- Balance terms (W)
-        source_term = -inner1d(deltah_rxn, rates) * 1000  # W/m**3
+        # source_term = -inner1d(deltah_rxn, rates) * 1000  # W/m**3
+        # TODO: Check if this is correct
+        # source_term = -np.dot(deltah_rxn, rates) * 1000  # W / m**3
+        source_term = -(deltah_rxn * rates).sum(axis=1) * 1000  # W / m**3
 
         if self.adiabatic:
             heat_transfer = 0
@@ -1645,7 +1654,11 @@ def energy_balances(self, time, mole_conc, vol_diff, temp, flow_in, rate_i,
         _, cp_j = self.Liquid_1.getCpPure(temp)
 
         # Volumetric heat capacity
-        cp_vol = inner1d(cp_j, mole_conc) * 1000  # J/m**3/K
+        # cp_vol = inner1d(cp_j, mole_conc) * 1000  # J/m**3/K
+        # TODO: Check if this is correct
+        # cp_vol = -np.dot(cp_j, mole_conc) * 1000  # J/m**3/K
+        cp_vol = (cp_j * mole_conc).sum(axis=1) * 1000  # vol in L
+
 
         # Heat of reaction
         delta_href = self.Kinetics.delta_hrxn
@@ -1656,7 +1669,10 @@ def energy_balances(self, time, mole_conc, vol_diff, temp, flow_in, rate_i,
             stoich, temp, self.mask_species, delta_href, tref_hrxn)  # J/mol
 
         # ---------- Balance terms (W)
-        source_term = -inner1d(deltah_rxn, rate_i * 1000)  # W/m**3
+        # source_term = -inner1d(deltah_rxn, rate_i * 1000)  # W/m**3
+        # TODO: Check if this is correct
+        # source_term = -np.dot(deltah_rxn, rate_i * 1000)  # W/m**3
+        source_term = -(deltah_rxn * rate_i * 1000).sum(axis=1)  # W/m**3
 
         temp_diff = np.diff(temp)
         flow_term = -flow_in * temp_diff / vol_diff  # K/s
@@ -1910,7 +1926,7 @@ def retrieve_results(self, time, states):
     def plot_steady(self, fig_size=None, title=None):
         fig, axes = plt.subplots(1, 2, figsize=fig_size)
 
-        # Concentration
+        # Concentrationn
         axes[0].plot(self.volPosition, self.concProfSteady)
         axes[0].set_ylabel('$C_j$ (mol/L)')
         axes[0].legend(self.name_species)

PharmaPy/Results.py

@@ -103,11 +103,11 @@ def pprint(di, name_items, fields, str_out=True):
 
     items = list(di.keys())
 
-    # ---------- Lenghts
-    # Lenght of first column
+    # ---------- Lengths
+    # Length of first column
     max_lens = {name_items: max([len(name) for name in items])}
 
-    # Lenght of remaining columns
+    # Length of remaining columns
     for field in fields:
         field_vals = [di[name].get(field, '') for name in items]