Skip to content
Snippets Groups Projects
Commit 4908ba43 authored by Peter Schubert's avatar Peter Schubert
Browse files

GbaIpoptProblem added (version 0.3.0)

parent 7c80ca04
Branches
No related tags found
No related merge requests found
...@@ -7,6 +7,7 @@ Package generating XBA problems from SBML. ...@@ -7,6 +7,7 @@ Package generating XBA problems from SBML.
from .model import * from .model import *
from .problems import * from .problems import *
from .solvers import *
from .utils import * from .utils import *
from . import model, problems, utils from . import model, problems, utils
...@@ -15,3 +16,4 @@ __all__ = [] ...@@ -15,3 +16,4 @@ __all__ = []
__all__ += model.__all__ __all__ += model.__all__
__all__ += problems.__all__ __all__ += problems.__all__
__all__ += utils.__all__ __all__ += utils.__all__
__all__ += solvers.__all__
""" """
Definition of version string. Definition of version string.
""" """
__version__ = "0.2.0" __version__ = "0.3.0"
program_name = 'xbanalysis' program_name = 'xbanalysis'
...@@ -15,10 +15,10 @@ from xbanalysis.utils.hessians import (hess_from_grad_grad, hess_from_v_hesss, h ...@@ -15,10 +15,10 @@ from xbanalysis.utils.hessians import (hess_from_grad_grad, hess_from_v_hesss, h
# Decorator to cache last result # Decorator to cache last result
def cache_last(func): def cache_last(func):
def wrapper(self, x): def wrapper(self, *args):
fname = wrapper.__name__ fname = wrapper.__name__
if fname not in GbaProblem.last_values or np.array_equal(GbaProblem.last_values[fname][0], x) is False: if fname not in GbaProblem.last_values or np.array_equal(GbaProblem.last_values[fname][0], args[0]) is False:
GbaProblem.last_values[fname] = [x, func(self, x)] GbaProblem.last_values[fname] = [args[0], func(self, args[0])]
return GbaProblem.last_values[fname][1] return GbaProblem.last_values[fname][1]
wrapper.__name__ = func.__name__ wrapper.__name__ = func.__name__
return wrapper return wrapper
...@@ -298,7 +298,7 @@ class GbaProblem: ...@@ -298,7 +298,7 @@ class GbaProblem:
enz_trans_hesss = types.FunctionType(hesss_code.co_consts[0], self.model_params) enz_trans_hesss = types.FunctionType(hesss_code.co_consts[0], self.model_params)
return enz_trans, enz_trans_jac, enz_trans_hesss return enz_trans, enz_trans_jac, enz_trans_hesss
# @cache_last @cache_last
def get_growth_rate_0(self, x): def get_growth_rate_0(self, x):
"""Calculate growth rate at x, assuming no protein degradation """Calculate growth rate at x, assuming no protein degradation
...@@ -323,6 +323,7 @@ class GbaProblem: ...@@ -323,6 +323,7 @@ class GbaProblem:
# self.last_gr0 = gr0 # self.last_gr0 = gr0
return gr0 return gr0
@cache_last
def get_growth_rate_0_grad(self, x): def get_growth_rate_0_grad(self, x):
"""Calculate growth rate gradient at x, assuming no protein degradation """Calculate growth rate gradient at x, assuming no protein degradation
...@@ -345,6 +346,7 @@ class GbaProblem: ...@@ -345,6 +346,7 @@ class GbaProblem:
gr0_grad = - gr0 ** 2 * tau_grad gr0_grad = - gr0 ** 2 * tau_grad
return np.nan_to_num(gr0_grad, nan=0.0) return np.nan_to_num(gr0_grad, nan=0.0)
@cache_last
def get_growth_rate_0_hess(self, x): def get_growth_rate_0_hess(self, x):
"""Calculate growth rate Hessian at x, assuming no protein degradation """Calculate growth rate Hessian at x, assuming no protein degradation
...@@ -389,6 +391,7 @@ class GbaProblem: ...@@ -389,6 +391,7 @@ class GbaProblem:
gr0_hess = - 2 * gr0 * hess_a - gr0 ** 2 * (hess_b1 + hess_b2 + hess_b3) gr0_hess = - 2 * gr0 * hess_a - gr0 ** 2 * (hess_b1 + hess_b2 + hess_b3)
return np.nan_to_num(gr0_hess, nan=0.0) return np.nan_to_num(gr0_hess, nan=0.0)
@cache_last
def get_growth_rate(self, x): def get_growth_rate(self, x):
"""Calculate growth rate at x, considering protein degradation """Calculate growth rate at x, considering protein degradation
...@@ -413,6 +416,7 @@ class GbaProblem: ...@@ -413,6 +416,7 @@ class GbaProblem:
gr = gr0 * (1.0 + enz_to_mras.dot(enz_to_pstmas)) gr = gr0 * (1.0 + enz_to_mras.dot(enz_to_pstmas))
return gr return gr
@cache_last
def get_growth_rate_grad(self, x): def get_growth_rate_grad(self, x):
"""Calculate growth rate gradient at x, considering protein degradation """Calculate growth rate gradient at x, considering protein degradation
...@@ -440,6 +444,7 @@ class GbaProblem: ...@@ -440,6 +444,7 @@ class GbaProblem:
+ gr0 * (enz_to_mras.dot(enz_to_pstmas_jac) + enz_to_pstmas.dot(enz_to_mras_jac))) + 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) return np.nan_to_num(gr_grad, nan=0.0)
@cache_last
def get_growth_rate_hess(self, x): def get_growth_rate_hess(self, x):
"""Calculate growth rate Hessian at x, considering protein degradation """Calculate growth rate Hessian at x, considering protein degradation
...@@ -490,6 +495,7 @@ class GbaProblem: ...@@ -490,6 +495,7 @@ class GbaProblem:
gr_hess = hess_a + 2 * hess_b + gr0 * (hess_c1 + hess_c2 + hess_c3 + hess_c4) 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) return np.nan_to_num(gr_hess, nan=0.0)
@cache_last
def get_alphas(self, x): def get_alphas(self, x):
"""Calculate ribosomal allocations at x, considering protein degradation """Calculate ribosomal allocations at x, considering protein degradation
...@@ -514,6 +520,7 @@ class GbaProblem: ...@@ -514,6 +520,7 @@ class GbaProblem:
alphas = alphas0 - enz_to_mras * enz_to_pstmas alphas = alphas0 - enz_to_mras * enz_to_pstmas
return alphas return alphas
@cache_last
def heq_mass_balance(self, x): def heq_mass_balance(self, x):
"""Calculate mass balance equations at x for metabolites """Calculate mass balance equations at x for metabolites
...@@ -548,6 +555,7 @@ class GbaProblem: ...@@ -548,6 +555,7 @@ class GbaProblem:
metab_mb = metabolic_mb + gr * (ps_drain - cm) + synth_enz_mb metab_mb = metabolic_mb + gr * (ps_drain - cm) + synth_enz_mb
return metab_mb return metab_mb
@cache_last
def heq_mass_balance_jac(self, x): def heq_mass_balance_jac(self, x):
"""Calculate Jacobian of mass balance equations at x for metabolites """Calculate Jacobian of mass balance equations at x for metabolites
...@@ -601,6 +609,7 @@ class GbaProblem: ...@@ -601,6 +609,7 @@ class GbaProblem:
metab_mb_jac = metabolic_jac + dilute_1_jac + dilute_2_jac + synth_enz_mb_jac 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) return np.nan_to_num(metab_mb_jac, nan=0.0)
@cache_last
def heq_mass_balance_hess(self, x): def heq_mass_balance_hess(self, x):
"""Calculate Hessians of mass balance equations at x """Calculate Hessians of mass balance equations at x
...@@ -681,15 +690,19 @@ class GbaProblem: ...@@ -681,15 +690,19 @@ class GbaProblem:
metab_mb_hesss = metabolic_hesss + dilute_1_hesss + 2 * dilute_2_hesss + synth_enz_mb_hesss 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) return np.nan_to_num(metab_mb_hesss, nan=0.0)
@cache_last
def heq_enz_trans(self, x): def heq_enz_trans(self, x):
return self.enz_trans(x) return self.enz_trans(x)
@cache_last
def heq_enz_trans_jac(self, x): def heq_enz_trans_jac(self, x):
return self.enz_trans_jac(x) return self.enz_trans_jac(x)
@cache_last
def heq_enz_trans_hess(self, x): def heq_enz_trans_hess(self, x):
return self.enz_trans_hesss(x) return self.enz_trans_hesss(x)
@cache_last
def heq_density(self, x): def heq_density(self, x):
"""Calculate density (g/l) constraint at x """Calculate density (g/l) constraint at x
...@@ -712,10 +725,10 @@ class GbaProblem: ...@@ -712,10 +725,10 @@ class GbaProblem:
density = self.model_params['mws'][self.var_mask_enz].dot(ce) density = self.model_params['mws'][self.var_mask_enz].dot(ce)
else: else:
density = self.model_params['mws'].dot(x) density = self.model_params['mws'].dot(x)
# !!! when using Ipopt, the rho could be specified as a bound! # TODO check if rho should be specified as a bound!
return np.array(self.model_params['rho'] - density) return np.array([self.model_params['rho'] - density])
# noinspection PyUnusedLocal @cache_last
def heq_density_jac(self, x): def heq_density_jac(self, x):
"""Calculate gradient of density constraint at x """Calculate gradient of density constraint at x
...@@ -732,13 +745,14 @@ class GbaProblem: ...@@ -732,13 +745,14 @@ class GbaProblem:
return density_grad return density_grad
# noinspection PyUnusedLocal # noinspection PyUnusedLocal
@cache_last
def heq_density_hess(self, x): def heq_density_hess(self, x):
"""Calculate hessian of density constraint at x """Calculate hessian of density constraint at x
:param x: optimization variables :param x: optimization variables
:type x: numpy.ndarray :type x: numpy.ndarray
:returns: mass balance Hessian at x :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) return np.tril(hess)
"""Subpackage with XBA model classes """
from .gba_ipopt_problem import GbaIpoptProblem
__all__ = ['GbaIpoptProblem']
"""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}')
...@@ -54,7 +54,7 @@ class OptResults: ...@@ -54,7 +54,7 @@ class OptResults:
self.ress['rho_m'][i] = rho[self.problem.var_mask_metab] self.ress['rho_m'][i] = rho[self.problem.var_mask_metab]
self.ress['rho_e'][i] = rho[self.problem.var_mask_enz] 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'] 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['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['psm_fluxes'][i] = scale_factor / (self.problem.pstm_scale * self.problem.pstmas(cx))
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment