From 4908ba43bb647ece712f23e19e2e799fab3e721c Mon Sep 17 00:00:00 2001
From: Peter Schubert <Peter.Schubert@hhu.de>
Date: Tue, 7 Jun 2022 19:43:00 +0200
Subject: [PATCH] GbaIpoptProblem added (version 0.3.0)
---
xbanalysis/__init__.py | 2 +
xbanalysis/_version.py | 2 +-
xbanalysis/problems/gba_problem.py | 32 +++++++----
xbanalysis/solvers/__init__.py | 5 ++
xbanalysis/solvers/gba_ipopt_problem.py | 70 +++++++++++++++++++++++++
xbanalysis/utils/opt_results.py | 2 +-
6 files changed, 102 insertions(+), 11 deletions(-)
create mode 100644 xbanalysis/solvers/__init__.py
create mode 100644 xbanalysis/solvers/gba_ipopt_problem.py
diff --git a/xbanalysis/__init__.py b/xbanalysis/__init__.py
index 57b2abe..c2e384f 100644
--- a/xbanalysis/__init__.py
+++ b/xbanalysis/__init__.py
@@ -7,6 +7,7 @@ Package generating XBA problems from SBML.
from .model import *
from .problems import *
+from .solvers import *
from .utils import *
from . import model, problems, utils
@@ -15,3 +16,4 @@ __all__ = []
__all__ += model.__all__
__all__ += problems.__all__
__all__ += utils.__all__
+__all__ += solvers.__all__
diff --git a/xbanalysis/_version.py b/xbanalysis/_version.py
index 2a1d8ea..120eb33 100644
--- a/xbanalysis/_version.py
+++ b/xbanalysis/_version.py
@@ -1,5 +1,5 @@
"""
Definition of version string.
"""
-__version__ = "0.2.0"
+__version__ = "0.3.0"
program_name = 'xbanalysis'
diff --git a/xbanalysis/problems/gba_problem.py b/xbanalysis/problems/gba_problem.py
index 3f73fc7..eaeb62b 100644
--- a/xbanalysis/problems/gba_problem.py
+++ b/xbanalysis/problems/gba_problem.py
@@ -15,10 +15,10 @@ from xbanalysis.utils.hessians import (hess_from_grad_grad, hess_from_v_hesss, h
# Decorator to cache last result
def cache_last(func):
- def wrapper(self, x):
+ def wrapper(self, *args):
fname = wrapper.__name__
- if fname not in GbaProblem.last_values or np.array_equal(GbaProblem.last_values[fname][0], x) is False:
- GbaProblem.last_values[fname] = [x, func(self, x)]
+ if fname not in GbaProblem.last_values or np.array_equal(GbaProblem.last_values[fname][0], args[0]) is False:
+ GbaProblem.last_values[fname] = [args[0], func(self, args[0])]
return GbaProblem.last_values[fname][1]
wrapper.__name__ = func.__name__
return wrapper
@@ -298,7 +298,7 @@ class GbaProblem:
enz_trans_hesss = types.FunctionType(hesss_code.co_consts[0], self.model_params)
return enz_trans, enz_trans_jac, enz_trans_hesss
-# @cache_last
+ @cache_last
def get_growth_rate_0(self, x):
"""Calculate growth rate at x, assuming no protein degradation
@@ -323,6 +323,7 @@ class GbaProblem:
# self.last_gr0 = gr0
return gr0
+ @cache_last
def get_growth_rate_0_grad(self, x):
"""Calculate growth rate gradient at x, assuming no protein degradation
@@ -345,6 +346,7 @@ class GbaProblem:
gr0_grad = - gr0 ** 2 * tau_grad
return np.nan_to_num(gr0_grad, nan=0.0)
+ @cache_last
def get_growth_rate_0_hess(self, x):
"""Calculate growth rate Hessian at x, assuming no protein degradation
@@ -389,6 +391,7 @@ class GbaProblem:
gr0_hess = - 2 * gr0 * hess_a - gr0 ** 2 * (hess_b1 + hess_b2 + hess_b3)
return np.nan_to_num(gr0_hess, nan=0.0)
+ @cache_last
def get_growth_rate(self, x):
"""Calculate growth rate at x, considering protein degradation
@@ -413,6 +416,7 @@ class GbaProblem:
gr = gr0 * (1.0 + enz_to_mras.dot(enz_to_pstmas))
return gr
+ @cache_last
def get_growth_rate_grad(self, x):
"""Calculate growth rate gradient at x, considering protein degradation
@@ -440,6 +444,7 @@ class GbaProblem:
+ gr0 * (enz_to_mras.dot(enz_to_pstmas_jac) + enz_to_pstmas.dot(enz_to_mras_jac)))
return np.nan_to_num(gr_grad, nan=0.0)
+ @cache_last
def get_growth_rate_hess(self, x):
"""Calculate growth rate Hessian at x, considering protein degradation
@@ -490,6 +495,7 @@ class GbaProblem:
gr_hess = hess_a + 2 * hess_b + gr0 * (hess_c1 + hess_c2 + hess_c3 + hess_c4)
return np.nan_to_num(gr_hess, nan=0.0)
+ @cache_last
def get_alphas(self, x):
"""Calculate ribosomal allocations at x, considering protein degradation
@@ -514,6 +520,7 @@ class GbaProblem:
alphas = alphas0 - enz_to_mras * enz_to_pstmas
return alphas
+ @cache_last
def heq_mass_balance(self, x):
"""Calculate mass balance equations at x for metabolites
@@ -548,6 +555,7 @@ class GbaProblem:
metab_mb = metabolic_mb + gr * (ps_drain - cm) + synth_enz_mb
return metab_mb
+ @cache_last
def heq_mass_balance_jac(self, x):
"""Calculate Jacobian of mass balance equations at x for metabolites
@@ -601,6 +609,7 @@ class GbaProblem:
metab_mb_jac = metabolic_jac + dilute_1_jac + dilute_2_jac + synth_enz_mb_jac
return np.nan_to_num(metab_mb_jac, nan=0.0)
+ @cache_last
def heq_mass_balance_hess(self, x):
"""Calculate Hessians of mass balance equations at x
@@ -681,15 +690,19 @@ class GbaProblem:
metab_mb_hesss = metabolic_hesss + dilute_1_hesss + 2 * dilute_2_hesss + synth_enz_mb_hesss
return np.nan_to_num(metab_mb_hesss, nan=0.0)
+ @cache_last
def heq_enz_trans(self, x):
return self.enz_trans(x)
+ @cache_last
def heq_enz_trans_jac(self, x):
return self.enz_trans_jac(x)
+ @cache_last
def heq_enz_trans_hess(self, x):
return self.enz_trans_hesss(x)
+ @cache_last
def heq_density(self, x):
"""Calculate density (g/l) constraint at x
@@ -712,10 +725,10 @@ class GbaProblem:
density = self.model_params['mws'][self.var_mask_enz].dot(ce)
else:
density = self.model_params['mws'].dot(x)
- # !!! when using Ipopt, the rho could be specified as a bound!
- return np.array(self.model_params['rho'] - density)
+ # TODO check if rho should be specified as a bound!
+ return np.array([self.model_params['rho'] - density])
- # noinspection PyUnusedLocal
+ @cache_last
def heq_density_jac(self, x):
"""Calculate gradient of density constraint at x
@@ -732,13 +745,14 @@ class GbaProblem:
return density_grad
# noinspection PyUnusedLocal
+ @cache_last
def heq_density_hess(self, x):
"""Calculate hessian of density constraint at x
:param x: optimization variables
:type x: numpy.ndarray
:returns: mass balance Hessian at x
- :rtype: numpy.ndarray (2D), shape [x,x]
+ :rtype: numpy.ndarray (3D), shape [1,x,x]
"""
- hess = np.zeros((self.n_vars, self.n_vars))
+ hess = np.zeros((1, self.n_vars, self.n_vars))
return np.tril(hess)
diff --git a/xbanalysis/solvers/__init__.py b/xbanalysis/solvers/__init__.py
new file mode 100644
index 0000000..37b8472
--- /dev/null
+++ b/xbanalysis/solvers/__init__.py
@@ -0,0 +1,5 @@
+"""Subpackage with XBA model classes """
+
+from .gba_ipopt_problem import GbaIpoptProblem
+
+__all__ = ['GbaIpoptProblem']
diff --git a/xbanalysis/solvers/gba_ipopt_problem.py b/xbanalysis/solvers/gba_ipopt_problem.py
new file mode 100644
index 0000000..480c309
--- /dev/null
+++ b/xbanalysis/solvers/gba_ipopt_problem.py
@@ -0,0 +1,70 @@
+"""Implementation of GbaIpoptProblem class.
+
+Peter Schubert, HHU Duesseldorf, June 2022
+"""
+
+import numpy as np
+
+
+# based on cyipopt examples
+class GbaIpoptProblem:
+
+ def __init__(self, gba_problem):
+ """Initialize GbaIpoptProblem.
+
+ :param gba_problem: GbaProblem from where to extract relevant parameters
+ :type gba_problem: GbaProblem
+ """
+ self.gba_problem = gba_problem
+ self._n_vars = len(gba_problem.var_initial)
+ self.report_freq = 0
+
+ self.nfev = 0
+ self.njev = 0
+ self.nit = 0
+ tmp_x = np.ones(self._n_vars) * .01
+ self.constraint_dims = np.array([len(gba_problem.heq_mass_balance(tmp_x)),
+ len(gba_problem.heq_enz_trans(tmp_x)),
+ len(gba_problem.heq_density(tmp_x))])
+
+ def objective(self, x):
+ self.nfev += 1
+ return -self.gba_problem.get_growth_rate(x)
+
+ def gradient(self, x):
+ self.njev += 1
+ return -self.gba_problem.get_growth_rate_grad(x)
+
+ def constraints(self, x):
+ return np.concatenate((self.gba_problem.heq_mass_balance(x),
+ self.gba_problem.heq_enz_trans(x),
+ self.gba_problem.heq_density(x)))
+
+ def jacobian(self, x):
+ return (np.vstack((self.gba_problem.heq_mass_balance_jac(x),
+ self.gba_problem.heq_enz_trans_jac(x),
+ self.gba_problem.heq_density_jac(x)))
+ ).flatten()
+
+ def hessianstructure(self):
+ return np.tril_indices(self._n_vars)
+
+ def hessian(self, x, lagrange, obj_factor):
+ hess = obj_factor * -self.gba_problem.get_growth_rate_hess(x)
+ hesss = np.vstack((self.gba_problem.heq_mass_balance_hess(x),
+ self.gba_problem.heq_enz_trans_hess(x),
+ self.gba_problem.heq_density_hess(x)))
+ hess += (lagrange.reshape((-1, 1, 1)) * hesss).sum(axis=0)
+ return hess[np.tril_indices(self._n_vars)]
+
+ def set_report_freq(self, mod_val=0):
+ self.report_freq = mod_val
+
+ # noinspection PyUnusedLocal
+ def intermediate(self, alg_mod, iter_count, obj_value, inf_pr, inf_du, mu,
+ d_norm, regularization_size, alpha_du, alpha_pr,
+ ls_trials):
+ self.nit = iter_count
+ if self.report_freq > 0:
+ if iter_count % self.report_freq == 0:
+ print(f'[{iter_count:5d}] growth rate: {-3600 *obj_value:.5f}')
diff --git a/xbanalysis/utils/opt_results.py b/xbanalysis/utils/opt_results.py
index df73f03..ec6b9af 100644
--- a/xbanalysis/utils/opt_results.py
+++ b/xbanalysis/utils/opt_results.py
@@ -54,7 +54,7 @@ class OptResults:
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 ??
+ # TODO extract reaction volume from kinetics to scale fluxes here
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))
--
GitLab