diff --git a/setup.py b/setup.py
index b4799408e7640c880eb3a5142ca9d18ba81acd97..4837556d61afef394108028c8f7e3ffde32869b3 100755
--- a/setup.py
+++ b/setup.py
@@ -39,7 +39,7 @@ setup(
'numpy >= 1.20.0',
'sympy >= 1.9.0',
'scipy >= 1.7.0',
- 'swiglpk >= 5.0.3',
+ 'swiglpk >= 5.0.5',
'sbmlxdf >= 0.2.7'],
python_requires=">=3.7",
keywords=['modeling', 'standardization', 'SBML'],
diff --git a/xbanalysis/model/xba_model.py b/xbanalysis/model/xba_model.py
index 29bcbccbd4b531cd1374efb9c81b6753d9f595eb..056a5c828e081ff12f1950a91b2135d93d3c4939 100644
--- a/xbanalysis/model/xba_model.py
+++ b/xbanalysis/model/xba_model.py
@@ -7,6 +7,7 @@ import os
import time
import numpy as np
import sbmlxdf
+from scipy.sparse import coo_matrix
from .xba_compartment import XbaCompartment
from .xba_function import XbaFunction
@@ -23,7 +24,8 @@ class XbaModel:
def __init__(self, sbml_file):
if os.path.exists(sbml_file) is False:
print(f'{sbml_file} does not exist')
- return
+ raise FileNotFoundError
+
print(f'loading: {sbml_file} (last modified: {time.ctime(os.path.getmtime(sbml_file))})')
sbml_model = sbmlxdf.Model(sbml_file)
if sbml_model.isModel is False:
@@ -76,3 +78,28 @@ class XbaModel:
if sid in map_metab_row:
stoic_mat[map_metab_row[sid], col] = stoic
return stoic_mat
+
+ def get_stoic_matrix_coo(self, sids, rids):
+ """retrieve stoichiometric sub-matrix [sids x rids].
+
+ :param sids: species ids - rows of sub-matrix
+ :type sids: list of strings
+ :param rids: reaction ids - columns of sub-matrix
+ :type rids: list of strings
+ :returns: stoichmetric sub-matrix in sparse format
+ :rtype: scipy.sparse.coo_matrix
+ """
+ # TODO: sparce matrix manipulation
+ sid2idx = {sid: idx for idx, sid in enumerate(sids)}
+ stoic_data = []
+ for col_idx, rid in enumerate(rids):
+ for sid, stoic in self.reactions[rid].reactants.items():
+ if sid in sids:
+ stoic_data.append([sid2idx[sid], col_idx, -stoic])
+ for sid, stoic in self.reactions[rid].products.items():
+ if sid in sids:
+ stoic_data.append([sid2idx[sid], col_idx, stoic])
+ coo_data = np.array(stoic_data).copy() # ???
+ stoic_mat_coo = coo_matrix((coo_data[:, 2], (coo_data[:, 0], coo_data[:, 1])),
+ shape=(len(sids), len(rids)))
+ return stoic_mat_coo
diff --git a/xbanalysis/problems/__init__.py b/xbanalysis/problems/__init__.py
index 3adc55401bc3b973d4b49b99d8f153dcef000be2..a7326aaac46b6415ea4bd89bc8c8e210d9b8256c 100644
--- a/xbanalysis/problems/__init__.py
+++ b/xbanalysis/problems/__init__.py
@@ -3,6 +3,5 @@
from .gba_problem import GbaProblem
from .rba_problem import RbaProblem
from .fba_problem import FbaProblem
-from .lp_problem import LpProblem
-__all__ = ['GbaProblem', 'FbaProblem', 'LpProblem', 'RbaProblem']
+__all__ = ['GbaProblem', 'FbaProblem', 'RbaProblem']
diff --git a/xbanalysis/problems/fba_problem.py b/xbanalysis/problems/fba_problem.py
index 539a47a580e295352d6523897ae79884404c0538..7cd2a55034302bad5d15f7d8839feddadb2a45e0 100644
--- a/xbanalysis/problems/fba_problem.py
+++ b/xbanalysis/problems/fba_problem.py
@@ -1,11 +1,7 @@
-"""Implementation of FBA Problem class.
-
-Peter Schubert, HHU Duesseldorf, June 2022
-"""
-
import numpy as np
+from scipy import sparse
-from xbanalysis.problems.lp_problem import LpProblem
+from xbanalysis.solvers.glpk_linear_problem import GlpkLinearProblem
class FbaProblem:
@@ -38,37 +34,47 @@ class FbaProblem:
self.xba_model = xba_model
# identify reaction excluding protein synthesis and protein degradation
- self.mr_rids = [r.id for r in xba_model.reactions.values()
- if r.sboterm.is_in_branch('SBO:0000205') is False and
- r.sboterm.is_in_branch('SBO:0000179') is False]
- self.mrid2idx = {mrid: idx for idx, mrid in enumerate(self.mr_rids)}
- self.reversible = np.array([self.xba_model.reactions[rid].reversible for rid in self.mr_rids])
-
- sids = set()
- for rid in self.mr_rids:
+ self.rids = np.array([r.id for r in xba_model.reactions.values()
+ if r.sboterm.is_in_branch('SBO:0000205') is False and
+ r.sboterm.is_in_branch('SBO:0000179') is False])
+ self.rid2idx = {mrid: idx for idx, mrid in enumerate(self.rids)}
+ self.flux_bounds = np.array([[xba_model.reactions[rid].fbc_lb for rid in self.rids],
+ [xba_model.reactions[rid].fbc_ub for rid in self.rids]])
+ self.reversible = np.array([self.xba_model.reactions[rid].reversible for rid in self.rids])
+ for idx in range(len(self.reversible)):
+ if self.flux_bounds[0, idx] < 0 and np.logical_not(self.reversible[idx]):
+ self.reversible[idx] = True
+ print(f'reaction {idx} set to reversible')
+ self.rids_rev = self.rids[self.reversible]
+ self.rid2idx_rev = {rid: idx + len(self.rids) for idx, rid in enumerate(self.rids_rev)}
+
+ used_sids = set()
+ for rid in self.rids:
for sid in xba_model.reactions[rid].reactants:
- sids.add(sid)
+ used_sids.add(sid)
for sid in xba_model.reactions[rid].products:
- sids.add(sid)
- self.sids = list(sids)
+ used_sids.add(sid)
+ self.sids = np.array(list(used_sids))
+ self.sid2idx = {sid: idx for idx, sid in enumerate(self.sids)}
- self.s_mat = xba_model.get_stoic_matrix(self.sids, self.mr_rids)
+ # sparse matrices faster and more efficient to handle
+ self.s_mat_coo = xba_model.get_stoic_matrix_coo(self.sids, self.rids)
self.is_fba_model = self.check_fba_model(xba_model)
if self.is_fba_model is False:
print('Error: not a FBA model - missing parameters')
return
- self.flux_bounds = np.array([[xba_model.reactions[rid].fbc_lb for rid in self.mr_rids],
- [xba_model.reactions[rid].fbc_ub for rid in self.mr_rids]])
- self.obj_coefs = None
for oid, fbc_obj in xba_model.fbc_objectives.items():
if fbc_obj.active is True:
- self.obj_id = fbc_obj.id
+ # self.obj_id = fbc_obj.id
self.obj_dir = fbc_obj.direction
- self.set_obj_coefs(fbc_obj.coefficiants)
+ self.objective = fbc_obj.coefficiants
break
+ # lp problem will only be created once required.
+ self.lp = None
+
def check_fba_model(self, xba_model):
"""Check if fba related data available in the XBA model
@@ -81,7 +87,7 @@ class FbaProblem:
r_check = {'flux lower bound': 'fbc_lb', 'flux upper bound': 'fbc_ub'}
for check, param in r_check.items():
- missing = [rid for rid in self.mr_rids if hasattr(xba_model.reactions[rid], param) is False]
+ missing = [rid for rid in self.rids if hasattr(xba_model.reactions[rid], param) is False]
if len(missing) > 0:
print(f'missing {check} in:', missing)
is_fba_model = False
@@ -97,7 +103,11 @@ class FbaProblem:
return is_fba_model
- def combine_fwd_bwd(self, vector):
+ def update_reversible_rids(self):
+ self.rids_rev = self.rids[self.reversible]
+ self.rid2idx_rev = {rid: idx + len(self.rids) for idx, rid in enumerate(self.rids_rev)}
+
+ def combine_fwd_rev(self, vector):
"""combine forward and backward flux information.
:param vector: contains data on all reactions + backward information for reversible reactions
@@ -105,179 +115,206 @@ class FbaProblem:
:return: combined flux information
:rtype: 1D ndarray of shape (len(mr_rids))
"""
- n_cols = len(self.mr_rids)
+ n_cols = len(self.rids)
combined = vector[:n_cols]
combined[self.reversible] -= vector[n_cols:]
return combined
- def create_fba_lp(self):
- """create Linear Programming Problem for FBA with split reactions
+ def get_model_status(self):
+ """ Return status on model
+
+ :returns: status information in a dictionary
+ :rtype: dict
+ """
+ status = {'dof': self.s_mat_coo.shape[1] - np.linalg.matrix_rank(self.s_mat_coo.todense()),
+ 'n_rids': len(self.rids), 'n_sids': len(self.sids)}
+ return status
+
+ def create_core_fba_lp(self):
+ """create a core Linear Programming Problem for FBA with split reactions
+
+ Core Problem contains coefficient matrix, variable and row bounds
+ Objective function yet needs to be added prior to any optimization
Splitting reactions in forward/backward directions reduces high values during FBA
+ configured bounds are >= 0
Splitting required for pFBA (minimize absolute flux values)
:return: LpProblem configured for FBA
:rtype: LpProblem
"""
- if self.is_fba_model is False:
- return None
-
- lp = LpProblem()
- # self.reversible = np.array([self.xba_model.reactions[rid].reversible for rid in self.mr_rids])
- fwd_bwd_s_mat = np.hstack((self.s_mat, -self.s_mat[:, self.reversible]))
- fwd_bwd_obj_coefs = np.hstack((self.obj_coefs, -self.obj_coefs[self.reversible]))
- fwd_bwd_bnds = np.hstack((np.fmax(self.flux_bounds, 0),
- np.vstack((np.fmax(-self.flux_bounds[1, self.reversible], 0),
- np.fmax(-self.flux_bounds[0, self.reversible], 0)))))
- row_bnds = np.zeros((fwd_bwd_s_mat.shape[0], 2))
- lp.configure(fwd_bwd_obj_coefs, fwd_bwd_s_mat, row_bnds, fwd_bwd_bnds, self.obj_dir)
+ fwd_rev_s_mat_coo = sparse.hstack([self.s_mat_coo,
+ -self.s_mat_coo.tocsc()[:, self.reversible]])
+ fwd_rev_col_bnds = np.hstack((np.fmax(self.flux_bounds[:], 0),
+ np.fmax(np.flip(-self.flux_bounds[:, self.reversible], axis=0), 0)))
+ row_bnds = np.zeros((fwd_rev_s_mat_coo.shape[0], 2))
+ lp = GlpkLinearProblem(fwd_rev_s_mat_coo, row_bnds, fwd_rev_col_bnds)
return lp
- def fba_optimize(self):
+ def delete_fba_lp(self):
+ """Delete the corresponding linear problem, when required"""
+ if self.lp is not None:
+ del self.lp
+ self.lp = None
+
+ def get_coefficients(self, objective):
+ """get coefficients for objective or constraint function
+
+ :param objective: objective/constraint coefficients
+ :type objective: dict with key: reaction id, value: stoichiometric coefficient
+ :return: coefficients aligned with the linear problem
+ :rtype: 1D ndarray with float
+ """
+ coefs = np.zeros(len(self.rids) + len(self.rids_rev))
+ for rid, coef in objective.items():
+ coefs[self.rid2idx[rid]] = coef
+ if rid in self.rids_rev:
+ coefs[self.rid2idx_rev[rid]] = -coef
+ return coefs
+
+ def set_lp_objective(self, objective, direction='maximize'):
+ """Set objective coefficients and direction in linear problem
+
+ :param objective: objective coefficients
+ :type objective: dict with key: reaction id, value: stoichiometric coefficient
+ :param direction: direction of optimization
+ :type direction: str (optional, default: 'maximize')
+ """
+ assert self.lp is not None
+ self.lp.set_obj_coefs(self.get_coefficients(objective))
+ self.lp.set_obj_dir(direction)
+
+ def remove_lp_objective(self, objective):
+ """remove objective coefficients
+
+ :param objective: objective coefficients
+ :type objective: dict with key: reaction id, value: stoichiometric coefficient
+ """
+ self.lp.remove_obj_coefs(self.get_coefficients(objective))
+
+ def fba_optimize(self, short_result=False):
"""FBA optimization of current problem
+ :param short_result: short result (only objective value and status) or long result
+ :type short_result: bool (optional, default: False)
:return: results information from the optimization
:rtype: dict
"""
- lp = self.create_fba_lp()
- if lp is None:
- return {'success': False, 'message': 'not an FBA model'}
- res = lp.solve()
- res['x'] = self.combine_fwd_bwd(res['x'])
- res['reduced_costs'] = self.combine_fwd_bwd(res['reduced_costs'])
+ if self.lp is None:
+ self.lp = self.create_core_fba_lp()
+ # print('Core LP Problem created')
+ assert self.lp is not None
+
+ self.set_lp_objective(self.objective, self.obj_dir)
+ res = self.lp.solve(short_result)
+ # print(res['success'], res['fun'])
+ self.remove_lp_objective(self.objective)
+
+ if short_result is False:
+ res['x'] = self.combine_fwd_rev(res['x'])
+ res['reduced_costs'] = self.combine_fwd_rev(res['reduced_costs'])
return res
- def pfba_optimize(self, fract_of_optimum=1.0):
- """ pfba analysis
+ def pfba_optimize(self, frac_of_opt=1.0):
+ """pFBA parsimoneous flux balance analysis.
- Reactions need to be split in forward/backard (to minimize absolute flux value)
+ 0 ≤ frac_of_opt ≤ 1.0 (for minimization problems. frac_of_opt > 0)
+ for maximizing objective, this sets lower bound to (frac_of_opt * optimum)
+ for minimization objective, this sets upper bound to (optimum / frac_of_opt)
- :param fract_of_optimum: value between 0.0 and 1.0.
- :type fract_of_optimum: float (default 1.0)
+ :param frac_of_opt: objective value within
+ :type frac_of_opt: float [0 ... 1] , optional (default: 1.0)`
:return:
+ :rtype: dict
"""
- lp = self.create_fba_lp()
- if lp is None:
- return {'success': False, 'message': 'not an FBA model'}
- res = lp.solve(short_result=False)
- if res['success']:
- # fix FBA objective as constraint
- fwd_bwd_obj_coefs = np.hstack((self.obj_coefs, -self.obj_coefs[self.reversible]))
- if self.obj_dir == 'maximize':
- row_bds = np.array([res['fun'] * fract_of_optimum, np.nan])
- else:
- row_bds = np.array([np.nan, res['fun'] / fract_of_optimum])
- lp.add_constraint(fwd_bwd_obj_coefs, row_bds)
- # new LP objective is to minimize sum of fluxes
- n_cols = len(self.mr_rids)
- lp.set_obj_coefs(np.ones(n_cols + sum(self.reversible)))
- lp.set_obj_dir('minimize')
- res = lp.solve()
- res['x'] = self.combine_fwd_bwd(res['x'])
- res['reduced_costs'] = self.combine_fwd_bwd(res['reduced_costs'])
+ # get optimum value for model objective function
+ res = self.fba_optimize(short_result=True)
+ if res['success'] is False:
+ return res
+
+ # implement model main objective as a constraint
+ if self.obj_dir == 'maximize':
+ row_bnds = np.array([res['fun'] * frac_of_opt, np.nan])
+ else:
+ row_bnds = np.array([np.nan, res['fun'] / frac_of_opt])
+ row_idx0 = self.lp.add_constraint(self.get_coefficients(self.objective), row_bnds)
+
+ # create pFBA objective (all ones)
+ pfba_obj_coefs = np.ones(len(self.rids) + len(self.rids_rev))
+ self.lp.set_obj_coefs(pfba_obj_coefs)
+ self.lp.set_obj_dir('minimize')
+ res = self.lp.solve(short_result=False)
+ # remove pFBA objective
+ self.lp.remove_obj_coefs(pfba_obj_coefs)
+ # remove main model objective from constraints matrix
+ self.lp.del_constraint(row_idx0)
+ res['x'] = self.combine_fwd_rev(res['x'])
+ res['reduced_costs'] = self.combine_fwd_rev(res['reduced_costs'])
return res
- def fva_optimize(self, rids=None, fract_of_optimum=1.0):
- """FVA - flux variability analysis for all or selected reactions.
+ def fva_optimize(self, rids=None, frac_of_opt=1.0):
+ """Flux variability analyzis on reactions within fraction of optimum
+
+ 0 ≤ frac_of_opt ≤ 1.0 (for minimization problems. frac_of_opt > 0)
+ for maximizing objective, this sets lower bound to (frac_of_opt * optimum)
+ for minimization objective, this sets upper bound to (optimum / frac_of_opt)
- :param rids: selected reaction ids to run FVA (alternatively, all)
- :type rids: list of strings (default: None - all reactions)
- :param fract_of_optimum: value between 0.0 and 1.0.
- :type fract_of_optimum: float (default 1.0)
- :return: Flux ranges and objective values
- :rtype: pandas DataFrame
+ :param rids: reaction ids to determine flux ranges
+ :type rids: list-like of str, optional (default: None, all reactions)
+ :param frac_of_opt: objective value within
+ :type frac_of_opt: float [0 ... 1] , optional (default: 1.0)
+ :return: rids, min and max flux values
+ :rtype: dict
"""
- lp = self.create_fba_lp()
- if lp is None:
- return {'success': False, 'message': 'not an FBA model'}
- res = lp.solve(short_result=True)
+ # get optimum value for model objective function
+ res = self.fba_optimize(short_result=True)
+ if res['success'] is False:
+ return res
+
+ # implement model main objective as a constraint
+ if self.obj_dir == 'maximize':
+ row_bnds = np.array([res['fun'] * frac_of_opt, np.nan])
+ else:
+ row_bnds = np.array([np.nan, res['fun'] / frac_of_opt])
+ row_idx0 = self.lp.add_constraint(self.get_coefficients(self.objective), row_bnds)
if rids is None:
- rids = self.mr_rids
+ rids = self.rids
flux_ranges = np.zeros((4, len(rids)))
- if res['success']:
- # fix FBA objective as constraint
- fwd_bwd_obj_coefs = np.hstack((self.obj_coefs, -self.obj_coefs[self.reversible]))
- if self.obj_dir == 'maximize':
- row_bds = np.array([res['fun'] * fract_of_optimum, np.nan])
- else:
- row_bds = np.array([np.nan, res['fun'] / fract_of_optimum])
- fba_constr_row = lp.add_constraint(fwd_bwd_obj_coefs, row_bds)
-
- n_cols = len(self.mr_rids)
- for idx, rid in enumerate(rids):
- # select reaction maximize/minimize flux
- new_obj_coefs = np.zeros(n_cols)
- new_obj_coefs[self.mrid2idx[rid]] = 1
- new_fwd_bwd_obj_coefs = np.hstack((new_obj_coefs, -new_obj_coefs[self.reversible]))
- lp.set_obj_coefs(new_fwd_bwd_obj_coefs)
-
- lp.set_obj_dir('minimize')
- res = lp.solve(short_result=True)
- if res['success']:
- flux_ranges[0, idx] = res['fun']
- flux_ranges[2, idx] = lp.get_row_value(fba_constr_row)
- else:
- print(f"FVA min failed for {rid}")
- flux_ranges[0, idx] = np.nan
-
- lp.set_obj_dir('maximize')
- res = lp.solve(short_result=True)
- if res['success']:
- flux_ranges[1, idx] = res['fun']
- flux_ranges[3, idx] = lp.get_row_value(fba_constr_row)
- else:
- print(f"FVA max failed for {rid}")
- flux_ranges[1, idx] = np.nan
+ for idx, rid in enumerate(rids):
+ # determine minimal flux for reaction rid
+ self.set_lp_objective({rid: 1}, 'minimize')
+ res = self.lp.solve(short_result=True)
+ if res['success']:
+ flux_ranges[0, idx] = res['fun']
+ flux_ranges[2, idx] = self.lp.get_row_value(row_idx0)
+ # determine maximal flux for reaction rid
+ self.lp.set_obj_dir('maximize')
+ res = self.lp.solve(short_result=True)
+ if res['success']:
+ flux_ranges[1, idx] = res['fun']
+ flux_ranges[3, idx] = self.lp.get_row_value(row_idx0)
+ # remove objective for reaction rid
+ self.remove_lp_objective({rid: 1})
+ # remove main model objective from constraints matrix
+ self.lp.del_constraint(row_idx0)
res = {'rids': rids, 'min': flux_ranges[0], 'max': flux_ranges[1],
'obj_min': flux_ranges[2], 'obj_max': flux_ranges[3]}
return res
- @property
- def obj_dir(self):
- return self.__obj_dir
-
- @obj_dir.setter
- def obj_dir(self, dir_string):
- if 'min' in dir_string:
- self.__obj_dir = 'minimize'
- else:
- self.__obj_dir = 'maximize'
-
- def get_flux_bounds(self):
- """retrieve flux bounds for all reactions in a dict.
-
- :returns: flux bounds configured in the FBA problem
- :rtype: dict (key: reaction id, value: list with two float values for lower/upper bound)
- """
- return {mrid: [bds[0], bds[1]] for mrid, bds in zip(self.mr_rids, self.flux_bounds.T)}
-
- def set_flux_bounds(self, flux_bounds):
- """selectivly set lower and upper flux bounds for given reactions.
-
- :param flux_bounds: new lower and upper flux bounds selected reactions
- :type flux_bounds: dict (key: reaction id, value: list with two float values for lower/upper bound)
- """
- for rid, (lb, ub) in flux_bounds.items():
- self.flux_bounds[0, self.mrid2idx[rid]] = lb
- self.flux_bounds[1, self.mrid2idx[rid]] = ub
-
- def get_obj_coefs(self):
- """retrieve the coefficients of the FBA objective function.
-
- :returns: objective coefficients
- :rtype: dict (key: reaction id, value: float)
- """
- return {rid: self.obj_coefs[idx] for rid, idx in self.mrid2idx.items()
- if self.obj_coefs[idx] != 0.0}
-
- def set_obj_coefs(self, coefficiants):
- """set the coefficients of the FBA objective function.
+ def copy(self):
+ """Copy the FbaModel instance (while removing linear problem).
- :param coefficiants: objective coefficients of selected reactions
- :type coefficiants: dict (key: reaction id, value: float)
+ A shallow copy, to support multiprocessing of FVA,
+ which do not impact FBA model attributes, except for linear problem
"""
- self.obj_coefs = np.zeros(len(self.mr_rids))
- for rid, coef in coefficiants.items():
- self.obj_coefs[self.mrid2idx[rid]] = coef
+ cls = self.__class__
+ result = cls.__new__(cls)
+ result.__dict__.update(self.__dict__)
+ result.__dict__['lp'] = None
+ return result
+
+ def __del__(self):
+ print('deletion of fba_model')
+ del self.lp
diff --git a/xbanalysis/problems/rba_problem.py b/xbanalysis/problems/rba_problem.py
index 11c6dd7ecfd7686a1cc95e736fcbbb5d522e5b0b..071369ae4f8ac2f90f4127901068ef9eab486bcf 100644
--- a/xbanalysis/problems/rba_problem.py
+++ b/xbanalysis/problems/rba_problem.py
@@ -14,7 +14,7 @@ Peter Schubert, HHU Duesseldorf, June 2022
import numpy as np
-from xbanalysis.problems.lp_problem import LpProblem
+from xbanalysis.solvers.lp_problem import LpProblem
class RbaProblem:
@@ -138,7 +138,6 @@ class RbaProblem:
res = {'success': False, 'fun': 0.0, 'message': 'not an RBA model'}
return res
- res = {}
# Preanalysis: find range of growth rates
gr = {'lb': 0.0, 'ub': np.nan, 'try': 1.0}
for i in range(self.max_iter):
@@ -147,6 +146,7 @@ class RbaProblem:
lp = LpProblem()
lp.configure(self.obj_coefs, constr_coefs, self.fix_row_bnds, self.col_bnds, self.obj_dir)
res = lp.solve(scaling=self.scaling)
+ del lp
if res['success'] and res['fun'] > 1e-10:
gr['lb'] = gr['try']
gr['try'] *= 10
@@ -157,7 +157,6 @@ class RbaProblem:
gr['try'] /= 10
if gr['lb'] > 0.0 and np.isfinite(gr['ub']):
break
-
if gr['lb'] == 0.0 or np.isnan(gr['ub']):
res = {'success': False, 'fun': 0.0, 'message': 'no solution found'}
return res
@@ -170,6 +169,7 @@ class RbaProblem:
lp = LpProblem()
lp.configure(self.obj_coefs, constr_coefs, self.fix_row_bnds, self.col_bnds, self.obj_dir)
res = lp.solve(scaling=self.scaling)
+ del lp
if res['success'] and res['fun'] > 1e-10:
gr['lb'] = gr['try']
else:
@@ -183,6 +183,7 @@ class RbaProblem:
lp = LpProblem()
lp.configure(self.obj_coefs, constr_coefs, self.fix_row_bnds, self.col_bnds, self.obj_dir)
res = lp.solve(scaling=self.scaling)
+ del lp
res['fun'] = gr['try']
return res
diff --git a/xbanalysis/problems/rba_problem_coo.py b/xbanalysis/problems/rba_problem_coo.py
new file mode 100644
index 0000000000000000000000000000000000000000..071369ae4f8ac2f90f4127901068ef9eab486bcf
--- /dev/null
+++ b/xbanalysis/problems/rba_problem_coo.py
@@ -0,0 +1,286 @@
+"""Implementation of GBA Problem class.
+
+ - several process machines supported (yet to be validated)
+ - non-enzymatic metabolic reactions supported (yet to be validated)
+
+ Yet to implement:
+ - enzyme degradation
+ - enzyme transitions
+ - several compartments constraints
+ - target concentrations and target fluxes
+
+Peter Schubert, HHU Duesseldorf, June 2022
+"""
+
+import numpy as np
+
+from xbanalysis.solvers.lp_problem import LpProblem
+
+
+class RbaProblem:
+
+ def __init__(self, xba_model):
+ """
+
+ :param xba_model:
+ """
+ self.xba_model = xba_model
+
+ # metabolic reactions and enzyme transitions (no FBA reactions, no protein synthesis, no degradation)
+ self.mr_rids = [r.id for r in xba_model.reactions.values() if r.sboterm.is_in_branch('SBO:0000167')
+ and r.sboterm.is_in_branch('SBO:0000179') is False]
+ self.mrid2enz = {rid: xba_model.reactions[rid].enzyme for rid in self.mr_rids}
+ self.enz_sids = [enz_sid for enz_sid in self.mrid2enz.values() if enz_sid is not None]
+ self.rev_enz_sids = [enz_sid for mr_rid, enz_sid in self.mrid2enz.items()
+ if enz_sid is not None and xba_model.reactions[mr_rid].reversible]
+ self.processes = {s.rba_process_machine: sid for sid, s in xba_model.species.items()
+ if hasattr(s, 'rba_process_machine')}
+ self.psenz_rids = [xba_model.species[enz_sid].ps_rid for enz_sid in self.enz_sids]
+ self.pspm_rids = [xba_model.species[enz_sid].ps_rid for enz_sid in self.processes.values()]
+ self.densities = {cid: c.rba_frac_gcdw for cid, c in xba_model.compartments.items()
+ if hasattr(c, 'rba_frac_gcdw')}
+
+ # metabolites and enzymes used in metabolic reactions and enzyme transitions
+ sids = set()
+ for rid in self.mr_rids:
+ for sid in xba_model.reactions[rid].reactants:
+ sids.add(sid)
+ for sid in xba_model.reactions[rid].products:
+ sids.add(sid)
+ used_sids = list(sids)
+ self.sids = [sid for sid in used_sids if xba_model.species[sid].constant is False]
+ self.ext_conc = {sid: xba_model.species[sid].initial_conc
+ for sid in used_sids if xba_model.species[sid].constant is True}
+
+ self.rba_cols = self.mr_rids + self.enz_sids + list(self.processes)
+ self.rba_rows = (self.sids + list(self.processes)
+ + [sid + '_feff' for sid in self.enz_sids]
+ + [sid + '_reff' for sid in self.rev_enz_sids]
+ + list(self.densities))
+
+ self.is_rba_model = self.check_rba_model(xba_model)
+ if self.is_rba_model is False:
+ print('Error: not a GBA model - missing parameters')
+ return
+
+ fix_m, var_m, row_bnds_m = self.construct_mass_balance(xba_model)
+ fix_pc, var_pc, row_bnds_pc = self.construct_process_capacities(xba_model)
+ fix_eff, var_eff, row_bnds_eff = self.construct_efficencies(xba_model)
+ fix_cd, var_cd, row_bnds_cd = self.construct_densities(xba_model)
+ self.fix_cd = fix_cd
+ self.fix_mat = np.vstack((fix_m, fix_pc, fix_eff, fix_cd))
+ self.var_mat = np.vstack((var_m, var_pc, var_eff, var_cd))
+ self.fix_row_bnds = np.vstack((np.zeros_like(row_bnds_m), np.zeros_like(row_bnds_pc),
+ row_bnds_eff, row_bnds_cd))
+ self.var_row_bnds = np.vstack((row_bnds_m, row_bnds_pc,
+ np.zeros_like(row_bnds_eff), np.zeros_like(row_bnds_cd)))
+
+ # in RBA we only check feasibility.
+ # Here we maximize enzyme amount and check that objective value e.g. > 1e-10 mmol/gcdw
+ self.obj_dir = 'maximize'
+ self.obj_coefs = np.hstack((np.zeros(len(self.mr_rids)),
+ np.ones(len(self.enz_sids) + len(self.processes))))
+ col_lbs = np.hstack(([xba_model.reactions[rid].fbc_lb for rid in self.mr_rids],
+ np.zeros(len(self.enz_sids) + len(self.processes))))
+ col_ubs = np.hstack(([xba_model.reactions[rid].fbc_ub for rid in self.mr_rids],
+ np.nan * np.ones(len(self.enz_sids) + len(self.processes))))
+ self.col_bnds = np.vstack((col_lbs, col_ubs))
+ self.scaling = '2N'
+ self.max_iter = 20
+
+ def check_rba_model(self, xba_model):
+ """Check if rba related data available in the XBA model
+
+ :param xba_model: Xba model with data extracted from SBML file
+ :type xba_model: XbaModel
+ :return: indicator if model contains data required for RBA problem formulation
+ :rtype: bool
+ """
+ is_rba_model = True
+
+ r_check = {'flux lower bound': 'fbc_lb', 'flux upper bound': 'fbc_ub', 'kapp forward': 'rba_kapp_f'}
+ for check, param in r_check.items():
+ missing = [rid for rid in self.mr_rids if hasattr(xba_model.reactions[rid], param) is False]
+ if len(missing) > 0:
+ print(f'missing {check} in:', missing)
+ is_rba_model = False
+
+ missing = [rid for rid in self.mr_rids if xba_model.reactions[rid].reversible
+ and hasattr(xba_model.reactions[rid], 'rba_kapp_r') is False]
+ if len(missing) > 0:
+ print('missing kapp reverse in:', missing)
+ is_rba_model = False
+
+ missing = [sid for sid in self.sids + self.enz_sids + list(self.processes.values())
+ if hasattr(xba_model.species[sid], 'mw') is False]
+ if len(missing) > 0:
+ print('missing molecular weight in:', missing)
+ is_rba_model = False
+
+ missing = [sid for sid in list(self.processes.values())
+ if hasattr(xba_model.species[sid], 'rba_process_capacity') is False]
+ if len(missing) > 0:
+ print('missing process capacity in:', missing)
+ is_rba_model = False
+
+ if len(self.densities) == 0:
+ print('missing a density constraint')
+ is_rba_model = False
+ return is_rba_model
+
+ def rba_optimize(self):
+ """
+
+ :return:
+ """
+ if self.is_rba_model is False:
+ print('Error: not a RBA model - missing parameters')
+ res = {'success': False, 'fun': 0.0, 'message': 'not an RBA model'}
+ return res
+
+ # Preanalysis: find range of growth rates
+ gr = {'lb': 0.0, 'ub': np.nan, 'try': 1.0}
+ for i in range(self.max_iter):
+ constr_coefs = self.fix_mat + gr['try'] * self.var_mat
+ # row_bnds = self.fix_row_bnds + gr['try'] * self.var_row_bnds
+ lp = LpProblem()
+ lp.configure(self.obj_coefs, constr_coefs, self.fix_row_bnds, self.col_bnds, self.obj_dir)
+ res = lp.solve(scaling=self.scaling)
+ del lp
+ if res['success'] and res['fun'] > 1e-10:
+ gr['lb'] = gr['try']
+ gr['try'] *= 10
+ else:
+ if gr['try'] <= 1e-4:
+ break
+ gr['ub'] = gr['try']
+ gr['try'] /= 10
+ if gr['lb'] > 0.0 and np.isfinite(gr['ub']):
+ break
+ if gr['lb'] == 0.0 or np.isnan(gr['ub']):
+ res = {'success': False, 'fun': 0.0, 'message': 'no solution found'}
+ return res
+
+ # find exact growth rate
+ for i in range(self.max_iter):
+ gr['try'] = 0.5 * (gr['lb'] + gr['ub'])
+ constr_coefs = self.fix_mat + gr['try'] * self.var_mat
+ # row_bnds = self.fix_row_bnds + gr['try'] * self.var_row_bnds
+ lp = LpProblem()
+ lp.configure(self.obj_coefs, constr_coefs, self.fix_row_bnds, self.col_bnds, self.obj_dir)
+ res = lp.solve(scaling=self.scaling)
+ del lp
+ if res['success'] and res['fun'] > 1e-10:
+ gr['lb'] = gr['try']
+ else:
+ gr['ub'] = gr['try']
+ if (gr['ub'] - gr['lb']) < gr['try'] * 1e-5:
+ break
+
+ # get optimum values for the final solution
+ gr['try'] = gr['lb']
+ constr_coefs = self.fix_mat + gr['try'] * self.var_mat
+ lp = LpProblem()
+ lp.configure(self.obj_coefs, constr_coefs, self.fix_row_bnds, self.col_bnds, self.obj_dir)
+ res = lp.solve(scaling=self.scaling)
+ del lp
+ res['fun'] = gr['try']
+ return res
+
+ def construct_mass_balance(self, xba_model):
+
+ fix_m_x_mr = xba_model.get_stoic_matrix(self.sids, self.mr_rids)
+ var_m_x_enz = xba_model.get_stoic_matrix(self.sids, self.psenz_rids)
+ var_m_x_pm = xba_model.get_stoic_matrix(self.sids, self.pspm_rids)
+ fix_m = np.hstack((fix_m_x_mr,
+ np.zeros((len(self.sids), len(self.enz_sids)+len(self.processes)))))
+ var_m = np.hstack((np.zeros((len(self.sids), len(self.mr_rids))),
+ var_m_x_enz,
+ var_m_x_pm))
+ targets = np.array([getattr(xba_model.species[sid], 'rba_target', 0.0) for sid in self.sids])
+ row_bnds_m = np.vstack((targets, targets)).T
+ return fix_m, var_m, row_bnds_m
+
+ def construct_process_capacities(self, xba_model):
+
+ var_pc_x_enz = np.zeros((len(self.processes), len(self.psenz_rids)))
+ for idx, pm in enumerate(self.processes):
+ var_pc_x_enz[idx] = [xba_model.species[enz_sid].rba_process_costs.get(pm, 0.0)
+ for enz_sid in self.enz_sids]
+ var_pc_x_pm = np.zeros((len(self.processes), len(self.processes)))
+ for idx, pm in enumerate(self.processes):
+ var_pc_x_pm[idx] = [xba_model.species[enz_sid].rba_process_costs.get(pm, 0.0)
+ for enz_sid in self.processes.values()]
+ fix_pc_x_pm = -np.diag([xba_model.species[enz_sid].rba_process_capacity
+ for enz_sid in self.processes.values()]) * 3600
+ fix_pc = np.hstack((np.zeros((len(self.processes), len(self.mr_rids) + len(self.enz_sids))),
+ fix_pc_x_pm))
+ var_pc = np.hstack((np.zeros((len(self.processes), len(self.mr_rids))),
+ var_pc_x_enz, var_pc_x_pm))
+ row_bnds_pc = np.zeros((len(self.processes), 2))
+ return fix_pc, var_pc, row_bnds_pc
+
+ def get_efficiencies(self, xba_model):
+
+ efficiencies = np.zeros((2, len(self.enz_sids)))
+ idx = 0
+ for mr_rid, enz_sid in self.mrid2enz.items():
+ if enz_sid is not None:
+ r = xba_model.reactions[mr_rid]
+ saturation = 1.0
+ if hasattr(r, 'rba_km_ext'):
+ for sid in r.reactants:
+ if sid in self.ext_conc:
+ saturation *= self.ext_conc[sid] / (self.ext_conc[sid] + r.rba_km_ext)
+ efficiencies[0, idx] = r.rba_kapp_f * 3600 * saturation
+ saturation = 1.0
+ if r.reversible:
+ if hasattr(r, 'rba_km_ext'):
+ for sid in r.products:
+ if sid in self.ext_conc:
+ saturation *= self.ext_conc[sid] / (self.ext_conc[sid] + r.rba_km_ext)
+ efficiencies[1, idx] = r.rba_kapp_r * 3600 * saturation
+ idx += 1
+ return efficiencies
+
+ def construct_efficencies(self, xba_model):
+
+ efficiencies = self.get_efficiencies(xba_model)
+ mask = np.array([enz_sid is not None for enz_sid in self.mrid2enz.values()])
+ fix_feff_x_mr = np.diag(np.ones(len(self.mrid2enz)))[mask]
+ fix_feff_x_enz = np.diag(-efficiencies[0])
+ fix_feff_x_pm = np.zeros((len(self.enz_sids), len(self.processes)))
+ fix_feff = np.hstack((fix_feff_x_mr, fix_feff_x_enz, fix_feff_x_pm))
+ row_bnds_feff = np.vstack((np.ones(len(self.enz_sids)) * np.nan, np.zeros(len(self.enz_sids)))).T
+
+ mask = [enz_sid is not None and xba_model.reactions[mr_rid].reversible
+ for mr_rid, enz_sid in self.mrid2enz.items()]
+ fix_reff_x_mr = np.diag(-np.ones(len(self.mrid2enz)))[mask]
+ mask = [enz_sid in self.rev_enz_sids for enz_sid in self.enz_sids]
+ fix_reff_x_enz = np.diag(-efficiencies[1])[mask]
+ fix_reff_x_pm = np.zeros((len(self.rev_enz_sids), len(self.processes)))
+ fix_reff = np.hstack((fix_reff_x_mr, fix_reff_x_enz, fix_reff_x_pm))
+ row_bnds_reff = np.vstack((np.ones(len(self.rev_enz_sids)) * np.nan, np.zeros(len(self.rev_enz_sids)))).T
+ fix_eff = np.vstack((fix_feff, fix_reff))
+ var_eff = np.zeros_like(fix_eff)
+ row_bnds_eff = np.vstack((row_bnds_feff, row_bnds_reff))
+ return fix_eff, var_eff, row_bnds_eff
+
+ def construct_densities(self, xba_model):
+ fix_cd_x_mr = np.zeros((len(self.densities), len(self.mr_rids)))
+ fix_cd_x_enz = np.zeros((len(self.densities), len(self.enz_sids)))
+ for row_idx, cid in enumerate(self.densities):
+ for col_idx, sid in enumerate(self.enz_sids):
+ if xba_model.species[sid].compartment == cid:
+ fix_cd_x_enz[row_idx, col_idx] = xba_model.species[sid].mw / 1000
+ fix_cd_x_pm = np.zeros((len(self.densities), len(self.processes)))
+ for row_idx, cid in enumerate(self.densities):
+ for col_idx, sid in enumerate(self.processes.values()):
+ if xba_model.species[sid].compartment == cid:
+ fix_cd_x_pm[row_idx, col_idx] = xba_model.species[sid].mw / 1000
+ fix_cd = np.hstack((fix_cd_x_mr, fix_cd_x_enz, fix_cd_x_pm))
+
+ var_cd = np.zeros_like(fix_cd)
+ row_bnds_cd = np.vstack((np.zeros(len(self.densities)),
+ [val for val in self.densities.values()])).T
+ return fix_cd, var_cd, row_bnds_cd
diff --git a/xbanalysis/solvers/__init__.py b/xbanalysis/solvers/__init__.py
index 37b8472d758df8fde32209a4a5daa2c4def29a96..bee7a56080a80855f0258e4ef44072aa3f7a4f14 100644
--- a/xbanalysis/solvers/__init__.py
+++ b/xbanalysis/solvers/__init__.py
@@ -1,5 +1,7 @@
"""Subpackage with XBA model classes """
+from .glpk_linear_problem import GlpkLinearProblem
+from .lp_problem import LpProblem
from .gba_ipopt_problem import GbaIpoptProblem
-__all__ = ['GbaIpoptProblem']
+__all__ = ['GlpkLinearProblem', 'GbaIpoptProblem', 'LpProblem']
diff --git a/xbanalysis/solvers/glpk_linear_problem.py b/xbanalysis/solvers/glpk_linear_problem.py
new file mode 100644
index 0000000000000000000000000000000000000000..a1c83b8f2e02f475ea14a19e445b404a15441d54
--- /dev/null
+++ b/xbanalysis/solvers/glpk_linear_problem.py
@@ -0,0 +1,315 @@
+"""Implementation of LpProblem (linear programming problem) class.
+
+built upon swiglpk
+
+Peter Schubert, HHU Duesseldorf, June 2022
+"""
+
+import numpy as np
+import swiglpk as glpk
+
+
+# status reported
+simplex_status = {
+ 0: 'LP problem instance successfully solved',
+ glpk.GLP_EBADB: 'Unable to start search, initial basis specified in the '
+ 'problem object is invalid. Number basic variables is not same as the number rows.',
+ glpk.GLP_ESING: 'Unable to start search, basis matrix corresponding to initial basis is singular.',
+ glpk.GLP_ECOND: 'Unable to start search, basis matrix corresponding to initial basis is ill-conditioned.',
+ glpk.GLP_EBOUND: 'Unable to start search, some double-bounded variables have incorrect bounds.',
+ glpk.GLP_EFAIL: 'Search prematurely terminated: solver failure.',
+ glpk.GLP_EOBJLL: 'Search prematurely terminated: objective function being '
+ 'maximized reached its lower limit and continues decreasing.',
+ glpk.GLP_EOBJUL: 'Search prematurely terminated: objective function being '
+ 'minimized reached its upper limit and continues increasing.',
+ glpk.GLP_EITLIM: 'Search prematurely terminated: simplex iteration limit exceeded.',
+ glpk.GLP_ETMLIM: 'Search prematurely terminated: time limit exceeded.',
+ glpk.GLP_ENOPFS: 'LP problem instance has no primal feasible solution.',
+ glpk.GLP_ENODFS: 'LP problem instance has no dual feasible solution.'}
+
+generic_status = {glpk.GLP_UNDEF: 'solution is undefined',
+ glpk.GLP_FEAS: 'solution is feasible',
+ glpk.GLP_INFEAS: 'solution is infeasible',
+ glpk.GLP_NOFEAS: 'problem has no feasible solution',
+ glpk.GLP_OPT: 'solution is optimal',
+ glpk.GLP_UNBND: 'problem has unbounded solution'}
+
+dual_status = {glpk.GLP_UNDEF: 'dual solution is undefined',
+ glpk.GLP_FEAS: 'dual solution is feasible',
+ glpk.GLP_INFEAS: 'dual solution is infeasible',
+ glpk.GLP_NOFEAS: 'no dual feasible solution exists'}
+
+prim_status = {glpk.GLP_UNDEF: 'primal solution is undefined',
+ glpk.GLP_FEAS: 'primal solution is feasible',
+ glpk.GLP_INFEAS: 'primal solution is infeasible',
+ glpk.GLP_NOFEAS: 'no primal feasible solution exists'}
+
+var_status = {glpk.GLP_BS: 'basic variable',
+ glpk.GLP_NL: 'non-basic variable on its lower bound',
+ glpk.GLP_NU: 'non-basic variable on its upper bound',
+ glpk.GLP_NF: 'non-basic free (unbounded) variable',
+ glpk.GLP_NS: 'non-basic fixed variable'}
+
+scaling_options = {'GM': glpk.GLP_SF_GM,
+ 'EQ': glpk.GLP_SF_EQ,
+ '2N': glpk.GLP_SF_2N,
+ 'SKIP': glpk.GLP_SF_SKIP,
+ 'AUTO': glpk.GLP_SF_AUTO}
+
+
+class GlpkLinearProblem:
+
+ def __init__(self, constr_coefs_coo, row_bnds, col_bnds):
+ """create core linear problem (constraints, variable and row bounds)
+
+ Problem needs to be properly deleted usin del(LpProblem)
+ Note: glpk related indices start at 1 (not zero)
+ :param constr_coefs_coo: matrix with constraint coefficients
+ :type constr_coefs_coo: scipy.sparce.coo_matrix with shape(nrows, ncols)
+ :param row_bnds: lower and upper bounds of auxhiliary (row) variables
+ :type row_bnds: 2D ndarray of shape(nrows, 2) (1st col: lower bounds, 2nd col: upper bounds)
+ :param col_bnds: lower and upper bounds of structural (col) variables
+ :type col_bnds: 2D ndarray of shape(2, ncols) (1st row: lower bounds, 2nd row: upper bounds)
+ """
+ glp_msg_lev = glpk.GLP_MSG_OFF # GLP_MSG_OFF, GLP_MSG_ERR, GLP_MSG_ON, GLP_MSG_ALL
+ glp_tm_lim = 10 * 1000
+ glp_presolve = glpk.GLP_OFF # or glpk.GLP_OFF (default) or GLP_ON
+ self.res = {}
+
+ self.nrows, self.ncols = constr_coefs_coo.shape
+ assert col_bnds.shape == (2, self.ncols)
+ assert row_bnds.shape == (self.nrows, 2)
+
+ self.lp = glpk.glp_create_prob()
+
+ # configure variables/columns and related bounds
+ glpk.glp_add_cols(self.lp, self.ncols)
+ for i, lb_ub in enumerate(col_bnds.T):
+ var_type, lb, ub = self._get_variable_type(lb_ub)
+ glpk.glp_set_col_bnds(self.lp, 1 + i, var_type, lb, ub)
+
+ # configure rows/constraints
+ glpk.glp_add_rows(self.lp, self.nrows)
+ for i, lb_ub in enumerate(row_bnds):
+ var_type, lb, ub = self._get_variable_type(lb_ub)
+ glpk.glp_set_row_bnds(self.lp, 1 + i, var_type, lb, ub)
+ n_coefs = len(constr_coefs_coo.row)
+ ia = glpk.intArray(1 + n_coefs)
+ ja = glpk.intArray(1 + n_coefs)
+ ar = glpk.doubleArray(1 + n_coefs)
+ for i in range(n_coefs):
+ ia[1 + i] = int(1 + constr_coefs_coo.row[i])
+ ja[1 + i] = int(1 + constr_coefs_coo.col[i])
+ ar[1 + i] = constr_coefs_coo.data[i]
+ glpk.glp_load_matrix(self.lp, n_coefs, ia, ja, ar)
+
+ # configure options
+ self._smcp = glpk.glp_smcp()
+ glpk.glp_init_smcp(self._smcp)
+ self._smcp.msg_lev = glp_msg_lev
+ self._smcp.tm_lim = glp_tm_lim
+ self._smcp.presolve = glp_presolve
+
+ @staticmethod
+ def _get_variable_type(bounds):
+ """Determine variable type for upper/lower bounds.
+
+ :param bounds: upper and lower bounds
+ :type bounds: 1D ndarray of length 2
+ :return var_type: glpk variable type to configure
+ :rtype: integer constant
+ :return lb: value of lower bound
+ :rtype float
+ :return ub: value of upper bound
+ :rtype float
+ """
+ lb, ub = bounds
+ if np.isfinite(lb) and np.isfinite(ub):
+ if lb < ub:
+ var_type = glpk.GLP_DB
+ else:
+ ub = 0.0
+ var_type = glpk.GLP_FX
+ elif not np.isfinite(lb) and not np.isfinite(ub):
+ lb = 0.0
+ ub = 0.0
+ var_type = glpk.GLP_FR
+ elif np.isfinite(ub):
+ lb = 0.0
+ var_type = glpk.GLP_UP
+ else:
+ ub = 0.0
+ var_type = glpk.GLP_LO
+ return var_type, lb, ub
+
+ def set_obj_coefs(self, obj_coefs):
+ """set coefficients for objective function
+
+ configure coefficients that are not zero
+
+ Note: glpk related indices start at 1 (not zero)
+ :param obj_coefs: coefficients of objective function
+ :type obj_coefs: 1D ndarray of shape (1, ncols)
+ """
+ assert len(obj_coefs) == self.ncols
+ for i in range(self.ncols):
+ if obj_coefs[i] != 0.0:
+ glpk.glp_set_obj_coef(self.lp, int(1 + i), obj_coefs[i])
+
+ def remove_obj_coefs(self, obj_coefs):
+ """remove coefficients from objective function
+
+ set coefficients that are not zero to zero in the liner problem
+
+ Note: glpk related indices start at 1 (not zero)
+ :param obj_coefs: coefficients of objective function
+ :type obj_coefs: list-like of floats
+ """
+ assert len(obj_coefs) == self.ncols
+ for i in range(self.ncols):
+ if obj_coefs[i] != 0.0:
+ glpk.glp_set_obj_coef(self.lp, int(1 + i), 0.0)
+
+ def set_obj_dir(self, direction):
+ """Configure the optimization direction
+
+ :param direction: direction for the optimization ('maximize' or 'minimize')
+ :type direction: string
+ """
+ obj_dir = glpk.GLP_MAX if 'max' in direction else glpk.GLP_MIN
+ glpk.glp_set_obj_dir(self.lp, obj_dir)
+
+ def add_constraint(self, constr_coefs, row_bnds):
+ """add a single constraint to the lp problem
+
+ :param constr_coefs: constraint coefficients for the new row
+ :type constr_coefs: 1d ndarray of length ncols
+ :param row_bnds: upper and lower bounds of the auxiliary (row) variable
+ :type row_bnds: 1D ndarray of length 2
+ :return: row index (start from 0)
+ :rtype: int
+ """
+ row_idx1 = glpk.glp_add_rows(self.lp, 1)
+ glpk.glp_set_row_bnds(self.lp, row_idx1, *self._get_variable_type(row_bnds))
+
+ n_coef = len(np.nonzero(constr_coefs)[0])
+ ind = glpk.intArray(n_coef + 1)
+ val = glpk.doubleArray(n_coef + 1)
+ for num_idx0, col_idx0 in enumerate(np.nonzero(constr_coefs)[0]):
+ ind[num_idx0 + 1] = int(col_idx0 + 1)
+ val[num_idx0 + 1] = constr_coefs[col_idx0]
+ glpk.glp_set_mat_row(self.lp, row_idx1, n_coef, ind, val)
+ self.nrows += 1
+ row_idx0 = row_idx1 - 1
+ return row_idx0
+
+ def del_constraint(self, row_idx0):
+ """remove single constraint from the lp problem
+
+ :param row_idx0: row index (start from 0) to be removed
+ :type row_idx0: int
+ """
+ assert row_idx0 < self.nrows
+ num = glpk.intArray(1 + 1)
+ num[1] = row_idx0 + 1
+ glpk.glp_del_rows(self.lp, 1, num)
+ self.nrows -= 1
+
+ def set_col_bnd(self, col_idx0, bounds):
+ """Set a single column (variable) bound
+ :param col_idx0: column index (start at 0) of variable
+ :type col_idx0: int
+ :param bounds: lower and upper bound to set
+ :type bounds: list with two float values
+ """
+ assert col_idx0 < self.ncols
+ glpk.glp_set_col_bnds(self.lp, col_idx0 + 1, *self._get_variable_type(bounds))
+
+ def remove_cols(self, col_idxs0):
+ """ remove variables (columns) from the linear problem
+
+ :param col_idxs0: col indices (start at 0) to be removed
+ :type col_idxs0: list or 1D ndarray of col indices
+ """
+ num = glpk.intArray(len(col_idxs0) + 1)
+ for num_idx0, col_idx0 in enumerate(col_idxs0):
+ assert col_idx0 < self.ncols
+ num[num_idx0 + 1] = col_idx0 + 1
+ glpk.glp_del_cols(self.lp, len(col_idxs0), num)
+ self.ncols -= len(col_idxs0)
+
+ def remove_rows(self, row_idxs0):
+ """ remove constraints (rows) from the linear problem.
+
+ :param row_idxs0: row indices (start at 0) to be removed
+ :type row_idxs0: list or 1D ndarray of row indices
+ """
+ num = glpk.intArray(len(row_idxs0) + 1)
+ for num_idx0, row_idx0 in enumerate(row_idxs0):
+ assert row_idx0 < self.nrows
+ num[num_idx0 + 1] = row_idx0 + 1
+ glpk.glp_del_rows(self.lp, len(row_idxs0), num)
+ self.nrows -= len(row_idxs0)
+
+ def construct_inital_basis(self):
+ """Re-construct an initial basis, e.g. after problem update
+ """
+ glpk.glp_adv_basis(self.lp, 0)
+
+ def solve(self, short_result=True, scaling=None):
+ """Solve the Linear Programming Problem.
+
+ :param short_result: short result (only objective value and status) or long result
+ :type short_result: bool (optional, default: True)
+ :param scaling: problem scaling options ('GM', 'EQ', '2N', 'SKIP', 'AUTO'
+ :type scaling: string (optional, default: None)
+ :return: result from linear optimization
+ :rtype: dict
+ """
+ if scaling is not None:
+ glpk.glp_scale_prob(self.lp, scaling_options.get(scaling, glpk.GLP_SF_AUTO))
+
+ simplex_result = glpk.glp_simplex(self.lp, self._smcp)
+
+ # in case of unsuccessful result, set basis and redo optimization
+ if glpk.glp_get_status(self.lp) != glpk.GLP_OPT:
+ glpk.glp_adv_basis(self.lp, 0)
+ # print('GLP ADV BASIS')
+ # glpk.glp_cpx_basis(self.lp)
+ simplex_result = glpk.glp_simplex(self.lp, self._smcp)
+
+ self.res = {
+ 'fun': glpk.glp_get_obj_val(self.lp),
+ 'success': glpk.glp_get_status(self.lp) == glpk.GLP_OPT,
+ 'generic_status': generic_status[glpk.glp_get_status(self.lp)],
+ }
+ if short_result is False:
+ self.res['x'] = np.array([glpk.glp_get_col_prim(self.lp, 1 + i) for i in range(self.ncols)])
+ self.res['reduced_costs'] = np.array([glpk.glp_get_col_dual(self.lp, 1 + i)
+ for i in range(self.ncols)])
+ self.res['shadow_prices'] = np.array([glpk.glp_get_row_dual(self.lp, 1 + i)
+ for i in range(self.nrows)])
+ self.res['simplex_result'] = simplex_result
+ self.res['simplex_message'] = simplex_status[simplex_result],
+ self.res['primal_status'] = prim_status[glpk.glp_get_prim_stat(self.lp)]
+ self.res['dual_status'] = dual_status[glpk.glp_get_dual_stat(self.lp)]
+ return self.res
+
+ def get_row_value(self, row_idx0):
+ """Retrieve right-hand side of a constraint (auxiliary variable)
+
+ :param row_idx0: row index into constraints matrix (starting with 0)
+ :type row_idx0: int
+ :return: value of auxiliary varibale (value of rhs)
+ :rtype: float
+ """
+ assert row_idx0 < glpk.glp_get_num_rows(self.lp)
+ return glpk.glp_get_row_prim(self.lp, row_idx0+1)
+
+ def __del__(self):
+ """Explicitly delete the LpProblem. (can we do without?)
+ """
+ # print('deletion of LpProblem')
+ assert self.lp is not None
+ if getattr(self, 'lp'):
+ glpk.glp_delete_prob(self.lp)
diff --git a/xbanalysis/problems/lp_problem.py b/xbanalysis/solvers/lp_problem.py
similarity index 81%
rename from xbanalysis/problems/lp_problem.py
rename to xbanalysis/solvers/lp_problem.py
index b52d80344ed2597b5fc11393fe0a2d9db43281ad..2e1479c67d33e82d09cf79858f26984e193243d5 100644
--- a/xbanalysis/problems/lp_problem.py
+++ b/xbanalysis/solvers/lp_problem.py
@@ -14,43 +14,43 @@ import swiglpk as glpk
# status reported
simplex_status = {
0: 'LP problem instance successfully solved',
- 1: 'Unable to start search, initial basis specified in the '
+ glpk.GLP_EBADB: 'Unable to start search, initial basis specified in the '
'problem object is invalid. Number basic variables is not same as the number rows.',
- 2: 'Unable to start search, basis matrix corresponding to initial basis is singular.',
- 3: 'Unable to start search, basis matrix corresponding to initial basis is ill-conditioned.',
- 4: 'Unable to start search, some double-bounded variables have incorrect bounds.',
- 5: 'Search prematurely terminated: solver failure.',
- 6: 'Search prematurely terminated: objective function being '
+ glpk.GLP_ESING: 'Unable to start search, basis matrix corresponding to initial basis is singular.',
+ glpk.GLP_ECOND: 'Unable to start search, basis matrix corresponding to initial basis is ill-conditioned.',
+ glpk.GLP_EBOUND: 'Unable to start search, some double-bounded variables have incorrect bounds.',
+ glpk.GLP_EFAIL: 'Search prematurely terminated: solver failure.',
+ glpk.GLP_EOBJLL: 'Search prematurely terminated: objective function being '
'maximized reached its lower limit and continues decreasing.',
- 7: 'Search prematurely terminated: objective function being '
+ glpk.GLP_EOBJUL: 'Search prematurely terminated: objective function being '
'minimized reached its upper limit and continues increasing.',
- 8: 'Search prematurely terminated: simplex iteration limit exceeded.',
- 9: 'Search prematurely terminated: time limit exceeded.',
- 10: 'LP problem instance has no primal feasible solution.',
- 11: 'LP problem instance has no dual feasible solution.'}
-
-generic_status = {1: 'solution is undefined',
- 2: 'solution is feasible',
- 3: 'solution is infeasible',
- 4: 'problem has no feasible solution',
- 5: 'solution is optimal',
- 6: 'problem has unbounded solution'}
-
-dual_status = {1: 'dual solution is undefined',
- 2: 'dual solution is feasible',
- 3: 'dual solution is infeasible',
- 4: 'no dual feasible solution exists'}
-
-prim_status = {1: 'primal solution is undefined',
- 2: 'primal solution is feasible',
- 3: 'primal solution is infeasible',
- 4: 'no primal feasible solution exists'}
-
-var_status = {1: 'basic variable',
- 2: 'non-basic variable on its lower bound',
- 3: 'non-basic variable on its upper bound',
- 4: 'non-basic free (unbounded) variable',
- 5: 'non-basic fixed variable'}
+ glpk.GLP_EITLIM: 'Search prematurely terminated: simplex iteration limit exceeded.',
+ glpk.GLP_ETMLIM: 'Search prematurely terminated: time limit exceeded.',
+ glpk.GLP_ENOPFS: 'LP problem instance has no primal feasible solution.',
+ glpk.GLP_ENODFS: 'LP problem instance has no dual feasible solution.'}
+
+generic_status = {glpk.GLP_UNDEF: 'solution is undefined',
+ glpk.GLP_FEAS: 'solution is feasible',
+ glpk.GLP_INFEAS: 'solution is infeasible',
+ glpk.GLP_NOFEAS: 'problem has no feasible solution',
+ glpk.GLP_OPT: 'solution is optimal',
+ glpk.GLP_UNBND: 'problem has unbounded solution'}
+
+dual_status = {glpk.GLP_UNDEF: 'dual solution is undefined',
+ glpk.GLP_FEAS: 'dual solution is feasible',
+ glpk.GLP_INFEAS: 'dual solution is infeasible',
+ glpk.GLP_NOFEAS: 'no dual feasible solution exists'}
+
+prim_status = {glpk.GLP_UNDEF: 'primal solution is undefined',
+ glpk.GLP_FEAS: 'primal solution is feasible',
+ glpk.GLP_INFEAS: 'primal solution is infeasible',
+ glpk.GLP_NOFEAS: 'no primal feasible solution exists'}
+
+var_status = {glpk.GLP_BS: 'basic variable',
+ glpk.GLP_NL: 'non-basic variable on its lower bound',
+ glpk.GLP_NU: 'non-basic variable on its upper bound',
+ glpk.GLP_NF: 'non-basic free (unbounded) variable',
+ glpk.GLP_NS: 'non-basic fixed variable'}
scaling_options = {'GM': glpk.GLP_SF_GM,
'EQ': glpk.GLP_SF_EQ,
@@ -69,6 +69,8 @@ class LpProblem:
self.lp = glpk.glp_create_prob()
self.ncols = 0
self.nrows = 0
+ self.msg_lev = glpk.GLP_MSG_ERR
+ self.tm_lim = 10 * 1000
self.res = {}
@staticmethod
@@ -245,12 +247,15 @@ class LpProblem:
:return: result information from the optimization stored in a dictionary
:rtype: dict
"""
+ smcp = glpk.glp_smcp()
+ glpk.glp_init_smcp(smcp)
+ smcp.msg_lev = self.msg_lev
+ smcp.tm_lim = self.tm_lim
+
if scaling is not None:
opt = scaling_options.get(scaling, glpk.GLP_SF_AUTO)
glpk.glp_scale_prob(self.lp, opt)
- # sjj = [glpk.glp_get_sjj(self.lp, 1+i) for i in range(glpk.glp_get_num_cols(self.lp))]
- # rii = [glpk.glp_get_rii(self.lp, 1+i) for i in range(glpk.glp_get_num_rows(self.lp))]
- simplex_result = glpk.glp_simplex(self.lp, None)
+ simplex_result = glpk.glp_simplex(self.lp, smcp)
return self.results(simplex_result, short_result)
def get_row_value(self, row):