OptResults added and small adjustments

xbanalysis/_version.py

 """
 Definition of version string.
 """
-__version__ = "0.1"
+__version__ = "0.2.0"
 program_name = 'xbanalysis'

xbanalysis/model/xba_parameter.py

@@ -3,6 +3,7 @@
 Peter Schubert, HHU Duesseldorf, May 2022
 """
 
+
 class XbaParameter:
 
     def __init__(self, s_parameter):

xbanalysis/model/xba_reaction.py

@@ -31,14 +31,14 @@ class XbaReaction:
 
     def set_enzyme(self, species):
         for sid in self.modifiers:
-            if species[sid].sboterm == 'SBO:0000250' or species[sid].sboterm == 'SBO:0000252':
+            if species[sid].sboterm in ('SBO:0000252', 'SBO:0000250'):
                 self.enzyme = sid
                 break
 
     def add_species_usage(self, species):
         # for metabolic reactions add to the species consumed/produced the metabolic reaction id
         # !!! we could also enter the Reaction Object references instead !!!!!!!!!!
-        if (self.sboterm == 'SBO:0000167') or (self.sboterm == 'SBO:0000179') or (self.sboterm == 'SBO:0000182'):
+        if self.sboterm in ('SBO:0000167', 'SBO:0000179', 'SBO:0000182'):
             for sid, val in self.products.items():
                 species[sid].add_m_reaction(self.id, val)
             for sid, val in self.reactants.items():
@@ -46,35 +46,31 @@ class XbaReaction:
         # protein degradation reaction, assume that first protein in reactants is getting degraded
         # a protein degradation reaction is still a metabolic reaction
         if self.sboterm == 'SBO:0000179':
-            for sid in self.reactants:
-                if (species[sid].sboterm == 'SBO:0000252') or (species[sid].sboterm == 'SBO:0000250'):
-                    species[sid].set_p_degrad_rid(self.id)
+            for sid, val in self.reactants.items():
+                if species[sid].sboterm in ('SBO:0000252', 'SBO:0000250'):
+                    species[sid].add_p_degrad_rid(self.id, -val)
                     break
-        # protein interconversion reactions
+        # protein transitions
         if self.sboterm == 'SBO:0000182':
-            for sid in self.reactants:
-                if (species[sid].sboterm == 'SBO:0000252') or (species[sid].sboterm == 'SBO:0000250'):
-                    species[sid].set_p_conversion_rid(self.id)
+            for sid, val in self.reactants.items():
+                if species[sid].sboterm in ('SBO:0000252', 'SBO:0000250'):
+                    species[sid].add_p_transition_rid(self.id, -val)
                     break
-            for sid in self.products:
-                if (species[sid].sboterm == 'SBO:0000252') or (species[sid].sboterm == 'SBO:0000250'):
-                    species[sid].set_p_conversion_rid(self.id)
+            for sid, val in self.products.items():
+                if species[sid].sboterm in ('SBO:0000252', 'SBO:0000250'):
+                    species[sid].add_p_transition_rid(self.id, val)
                     break
 
         # for protein synthesis reactions identify the enzyme being produced from reaction.products
         #     to species being consumed add the enzyme product with corresponding stoichiometry
         #     to species being produced add the enzyme product with corresponding stoichiometry, unless
         #        the species is itself the enzyme product, in which case we add the ps reaction producing it.
-        if (self.sboterm == 'SBO:0000589') or (self.sboterm == 'SBO:0000205') :
-            enz_sid = ''
-            for sid in self.products:
-                if (species[sid].sboterm == 'SBO:0000252') or (species[sid].sboterm == 'SBO:0000250'):
-                    enz_sid = sid
-                    self.ps_product = enz_sid
-                    species[enz_sid].set_ps_rid(self.id)
+        if self.sboterm in ('SBO:0000589', 'SBO:0000205'):
+            for sid, val in self.products.items():
+                if species[sid].sboterm in ('SBO:0000252', 'SBO:0000250'):
+                    self.ps_product = sid
                     break
             for sid, val in self.products.items():
-                if sid != enz_sid:
-                    species[sid].add_ps_reaction(enz_sid, val)
+                species[sid].add_ps_reaction(self.id, val)
             for sid, val in self.reactants.items():
-                species[sid].add_ps_reaction(enz_sid, -val)
+                species[sid].add_ps_reaction(self.id, -val)

xbanalysis/model/xba_species.py

@@ -15,14 +15,18 @@ class XbaSpecies:
         self.compartment = s_species['compartment']
         self.initial_conc = s_species['initialConcentration']
         self.constant = s_species['constant']
+        self.boundary = s_species['boundaryCondition']
         self.units = s_species['substanceUnits']
         attrs = sbmlxdf.misc.extract_xml_attrs(s_species['xmlAnnotation'], token='molecule')
-        self.weight = float(attrs['weight_Da'])
+        self.mw = float(attrs['weight_Da'])
+        self.ps_rid = None
         self.m_reactions = {}
         self.ps_reactions = {}
-        self.ps_rid = None
-        self.p_degrad_rid = None
-        self.p_convert_rid = None
+        self.p_degrad_rids = {}
+        self.p_trans_rids = {}
+
+    def set_ps_rid(self, rid, stoic):
+        self.ps_rid = {'rid': rid, 'stoic': stoic}
 
     def add_m_reaction(self, rid, stoic):
         self.m_reactions[rid] = stoic
@@ -30,11 +34,8 @@ class XbaSpecies:
     def add_ps_reaction(self, rid, stoic):
         self.ps_reactions[rid] = stoic
 
-    def set_ps_rid(self, rid):
-        self.ps_rid = rid
-
-    def set_p_degrad_rid(self, rid):
-        self.p_degrad_rid = rid
+    def add_p_degrad_rid(self, rid, stoic):
+        self.p_degrad_rids[rid] = stoic
 
-    def set_p_conversion_rid(self, rid):
-        self.p_convert_rid = rid
+    def add_p_transition_rid(self, rid, stoic):
+        self.p_trans_rids[rid] = stoic

xbanalysis/utils/__init__.py

 """Subpackage with utility functions """
 
-#rom .xba_model import XBAmodel
-
-__all__ = []
+from .opt_results import OptResults
 
+__all__ = ['OptResults']

xbanalysis/utils/hessians.py

+"""Implementation of Hessians calculations.
+
+Peter Schubert, HHU Duesseldorf, May 2022
+"""
+
+import numpy as np
+
+
+def hess_from_grad_grad(grad_l, grad_r):
+    """Calculate Hessian from two gradiends wrt x
+
+    :param grad_l: left scalar gradient wrt x variables
+    :param grad_l: numpy.ndarray (1D) with length of x
+    :param grad_r: rigth scalar gradient wrt x variables
+    :param grad_r: numpy.ndarray (1D) with length of x
+    :returns: result_hess Hessian wrt to x variables, lower triangle matrix
+    :rtype: numpy.ndarray (2D) with shape [x,x]
+    """
+    assert grad_l.shape[0] == grad_r.shape[0]
+    size = grad_l.shape[0]
+    result_hess = np.zeros((size, size))
+    for i in range(size):
+        for j in range(i+1):
+            result_hess[i, j] = grad_l[i] * grad_r[j]
+    return result_hess
+
+
+def hess_from_v_hesss(v, hesss):
+    """Calculate Hessian from a vector and Hessians
+
+    :param v: vector of length n
+    :param v: numpy.ndarray (1D)
+    :param hesss: Hessians of n vectors wrt to x variables
+    :param hesss: numpy.ndarray (3D) with shape [n,x,x]
+    :returns: result_hess Hessian wrt to x variables, lower triangle matrix
+    :rtype: numpy.ndarray (2D) with shape [x,x]
+    """
+    assert hesss.shape[1] == hesss.shape[2]
+    size = hesss.shape[1]
+    result_hess = np.zeros((size, size))
+    for i in range(size):
+        for j in range(i + 1):
+            result_hess[i, j] = v.dot(hesss[:, i, j])
+    return result_hess
+
+
+def hess_from_jac_jac(jac_l, jac_r):
+    """Calculate Hessian from two Jacobians
+
+    :param jac_l: left Jacobian for vector of size n wrt to x variables
+    :param jac_l: numpy.ndarray (2D) with shape [n,x]
+    :param jac_r: right Jacobian for vector of size n wrt to x variables
+    :param jac_r: numpy.ndarray (2D) with shape [n,x]
+    :returns: result_hess Hessian wrt to x variables, lower triangle matrix
+    :rtype: numpy.ndarray (2D) with shape [x,x]
+    """
+    assert jac_l.shape == jac_r.shape
+    size = jac_l.shape[1]
+    result_hess = np.zeros((size, size))
+    for i in range(size):
+        for j in range(i + 1):
+            result_hess[i, j] = jac_l[:, i].dot(jac_r[:, j])
+    return result_hess
+
+
+def hesss_from_matrix_hesss(matrix, hesss):
+    """Calculate dot product of a stoichiometric matrix with Hessians
+
+    :param matrix: stoichiometric sub-matrix with shape (m, n)
+    :param matrix: numpy.ndarray (2D) with shape [m,n]
+    :param hesss: Hessians of n vectors wrt to x variables
+    :param hesss: numpy.ndarray (3D) with shape [n,x,x]
+    :returns: result_hess Hessians wrt to x variables, lower triangle matrix
+    :rtype: numpy.ndarray (3D) with shape [n,x,x]
+    """
+    assert matrix.shape[1] == hesss.shape[0], "expect as many Hessians as there are matrix columns"
+    assert hesss.shape[1] == hesss.shape[2]
+    size = hesss.shape[1]
+    n_rows, n_cols = matrix.shape
+    result_hesss = np.zeros((n_rows, size, size))
+    for i in range(n_rows):
+        for j in range(n_cols):
+            if matrix[i, j] != 0:
+                result_hesss[i] += matrix[i, j] * hesss[j]
+    return result_hesss
+
+
+def hesss_from_v_hess(v, hess):
+    """Calculate Hessians from vector and a Hessian
+
+    :param v: vector of length m
+    :param v: numpy.ndarray (1D)
+    :param hess: Hessian of scalar wrt to x variables
+    :param hess: numpy.ndarray (2D) with shape [x,x]
+    :returns: result_hesss Hessians wrt to x variables, lower triangle matrix
+    :rtype: numpy.ndarray (3D) with shape [n,x,x]
+    """
+    assert hess.shape[0] == hess.shape[1]
+    size = hess.shape[1]
+    n_rows = v.shape[0]
+
+    result_hesss = np.zeros((n_rows, size, size))
+    for i in range(n_rows):
+        result_hesss[i] = v[i] * hess
+    return result_hesss
+
+
+def hesss_from_grad_jac(grad, jac):
+    """Calculate Hessians from gradient and a Jacobian
+
+    :param grad: gradient of scalar wrt to x variables
+    :param grad: numpy.ndarray (1D) with shape [x]
+    :param jac: Jacobian of a vector size m wrt to x variables
+    :param jac: numpy.ndarray (2D) with shape [m,x]
+    :returns: result_hesss Hessians wrt to x variables, lower triangle matrix
+    :rtype: numpy.ndarray (3D) with shape [m,x,x]
+    """
+    assert grad.shape[0] == jac.shape[1]
+    size = jac.shape[1]
+    n_rows = jac.shape[0]
+
+    result_hesss = np.zeros((n_rows, size, size))
+    for i in range(n_rows):
+        result_hesss[i] = hess_from_grad_grad(jac[i], grad)
+    return result_hesss

xbanalysis/utils/opt_results.py

+"""Implementation of OptResults class.
+
+Peter Schubert, HHU Duesseldorf, June 2022
+"""
+
+import numpy as np
+
+
+class OptResults:
+
+    def __init__(self, problem, n_points):
+        """Initialize OptResults.
+
+        :param problem: GbaProblem from where to extract relevant parameters
+        :type problem: GbaProblem
+        :param n_points: number of points to record
+        :type n_points: integer
+        """
+        self.problem = problem
+        self._n_points = n_points
+        self._idx = 0
+        self._n_m = len(self.problem.var_metab_sids)
+        self._n_e = len(self.problem.var_enz_sids)
+
+        self.ress = {}
+        self.ress['metabolites'] = self.problem.var_metab_sids
+        self.ress['enzymes'] = self.problem.var_enz_sids
+        self.ress['mr_reactions'] = self.problem.mr_rids
+        self.ress['enz_sid2idx'] = {sid: idx for idx, sid in enumerate(self.problem.var_enz_sids)}
+        self.ress['var_mask_metab'] = self.problem.var_mask_metab
+        self.ress['var_mask_enz'] = self.problem.var_mask_enz
+
+        self.ress['grs_per_h'] = np.zeros(n_points)
+        self.ress['status'] = np.zeros(n_points)
+        self.ress['x_opt'] = np.zeros((n_points, self.problem.n_vars))
+        self.ress['cm_M'] = np.zeros((n_points, self._n_m))
+        self.ress['rho_m'] = np.zeros((n_points, self._n_m))
+        self.ress['ce_M'] = np.zeros((n_points, self._n_e))
+        self.ress['rho_e'] = np.zeros((n_points, self._n_e))
+        self.ress['alphas'] = np.zeros((n_points, self._n_e))
+        self.ress['mr_fluxes'] = np.zeros((n_points, len(self.problem.mr_rids)))
+        self.ress['psm_fluxes'] = np.zeros((n_points, self._n_e))
+
+    def add(self, info):
+        if self._idx < self._n_points:
+            i = self._idx
+            cx = info['x']
+            self.ress['x_opt'][i] = cx
+            self.ress['status'][i] = info['status']
+            self.ress['grs_per_h'][i] = -3600 * info['obj_val']
+            self.ress['cm_M'][i] = cx[self.problem.var_mask_metab]
+            self.ress['ce_M'][i] = cx[self.problem.var_mask_enz]
+            rho = self.problem.model_params['mws'] * cx
+            self.ress['rho_m'][i] = rho[self.problem.var_mask_metab]
+            self.ress['rho_e'][i] = rho[self.problem.var_mask_enz]
+
+            # !!! what volume to consider, possibly the volume of the ribosome ??
+            scale_factor = 3600 * 1e3 / self.problem.model_params['rho']
+            self.ress['mr_fluxes'][i] = scale_factor * self.problem.mras(cx) / self.problem.var_vols[-1]
+            self.ress['psm_fluxes'][i] = scale_factor / (self.problem.pstm_scale * self.problem.pstmas(cx))
+
+            self.ress['alphas'][i] = self.problem.get_alphas(cx)
+            self._idx += 1
+
+    def get(self):
+        return self.ress