Multiprocessing implemented for FVA (validated on macOS)

modelpruner/core/model_pruner.py

@@ -5,7 +5,6 @@ import os
 import re
 import time
 import math
-import copy
 import numpy as np
 import multiprocessing as mp
 import json
@@ -14,30 +13,72 @@ from modelpruner.models.fba_model import FbaModel
 from modelpruner.core.protected_parts import ProtectedParts
 
 global _fba_model
-global _pf
+global _pfs
 
 
-def init_worker(fba_model, pf):
+def _init_worker(fba_model, pfs):
+    """Multiprocessing initializer for workers.
+
+    Configure global parameters for the worker process
+
+    :param fba_model: Fba model to apply FVA
+    :type fba_model: Class FbaModel
+    :param pfs: Protected Functions
+    :type pfs: dict, key: function id, value: class ProtectedFunction
+    """
     global _fba_model
-    global _pf
+    global _pfs
 
     _fba_model = fba_model.copy()
-    _pf = pf
+    _pfs = pfs
+
+
+def _worker(rids):
+    """Mulitprocessing worker to run FVA for selected reactions over all conditions.
 
+    Checked that best performance achieved when worker performs FVA for
+    all protected functions and a larger range for reactions.
 
-def worker(rid):
+    :param rids: reactions for FVA
+    :type rids: list-like of str
+    """
     global _fba_model
-    global _pf
+    global _pfs
 
-    return _fba_model.fva_pf_flux_ranges(_pf, [rid])
+    flux_min_pfs = np.zeros((len(_pfs), len(rids)))
+    flux_max_pfs = np.zeros((len(_pfs), len(rids)))
+
+    for idx, pf in enumerate(_pfs.values()):
+        res = _fba_model.fva_pf_flux_ranges(pf, rids)
+        if res['success'] is True:
+            flux_min_pfs[idx] = res['min']
+            flux_max_pfs[idx] = res['max']
+        else:
+            print(f'FBA for condition {pf.obj_id} is not feasible')
+            flux_min_pfs[idx] = np.nan
+            flux_max_pfs[idx] = np.nan
+            break
+    # delete linear problem
+    # del _fba_model.delete_fba_lp()
+    res = {'success': True, 'rids': rids, 'min_pfs': flux_min_pfs, 'max_pfs': flux_max_pfs}
+    return res
 
 
 class ModelPruner:
 
     def __init__(self, sbml_fname, protected_parts, resume=False):
-
+        """Init
+
+        :param sbml_fname: pathname of original model (ending with '.xml')
+        :type sbml_fname: str
+        :param protected_parts: pathname of protected parts file or dict
+        :type protected_parts: str or dict with protected parts information
+        :param resume: Indicator if to resume pruning from stored snapshot
+        :type resume: bool, optional (default: False)
+        """
         self.tolerance = 1e-7   # later accessible by set_params() ?
-        self.cpus = int(os.cpu_count() * 0.8)
+        self.cpus = int(os.cpu_count() - 0)
+        self.min_mp_total_flux = 1000
 
         self._snapshot_sbml = re.sub(r'.xml$', '_snapshot.xml', sbml_fname)
         self._snapshot_json = re.sub(r'.xml$', '_snapshot.json', sbml_fname)
@@ -52,7 +93,7 @@ class ModelPruner:
             with open(self._snapshot_json, 'r') as f:
                 snapshot_data = json.load(f)
             self.protected_rids = set(snapshot_data['protected_rids'])
-            print(f'snapshot restored at {len(self.fba_model.rids):4d} remaining reactions '
+            print(f'Snapshot restored at {len(self.fba_model.rids):4d} remaining reactions '
                   f'from: {self._snapshot_json}')
         else:
             self.protected_rids = self.nrp.protected_rids
@@ -62,19 +103,24 @@ class ModelPruner:
         self.fba_model.update_model_flux_bounds(self.nrp.overwrite_bounds)
         self.fba_model.update_model_reversible(self.get_pp_reversible_reactions())
 
-        print('determine optimal values for protected functions')
+        print('Determine optimal values for protected functions')
         self.set_function_optimum()
         print('Check consistency of protected parts')
         self.check_protected_parts()
 
     def get_total_protected_rids(self):
+        """Determine protected reactions, incl. those leading to protected metabolites.
+
+        :return: total protected reaction ids
+        :rtype: set of str
+
+        """
         # identify reactions that need protection due to protected metabolites
         #  over and above the protected rids provided in protected parts
         sid_mask = np.array([True if sid in self.protected_sids else False
                              for sid in self.fba_model.sids])
-        psf_mask = np.any(self.fba_model.s_mat[sid_mask] != 0, axis=0)
+        psf_mask = self.fba_model.s_mat_coo.tocsr()[sid_mask].getnnz(axis=0)
         protected_sid_rids = set(self.fba_model.rids[psf_mask])
-
         protected_rids = self.nrp.protected_rids.union(protected_sid_rids)
         return protected_rids
 
@@ -109,7 +155,10 @@ class ModelPruner:
                     pf.set_target_ub(pf.target_frac * pf.optimal)
 
     def check_protected_parts(self):
+        """Report on consisteny issues wrt to protected parts.
 
+        E.g. report on reactions and metabolties that do not exist in the model
+        """
         extra_p_rids = self.protected_rids.difference(set(self.fba_model.rids))
         extra_p_sids = self.protected_sids.difference(set(self.fba_model.sids))
         pf_rids = set()
@@ -117,27 +166,27 @@ class ModelPruner:
             pf_rids |= {rid for rid in pf.objective}
             pf_rids |= {rid for rid in pf.overwrite_bounds}
             if pf.fba_success is False:
-                print(f'protected function {pf.obj_id} has unsuccessful FBA '
+                print(f'Protected function {pf.obj_id} has unsuccessful FBA '
                       f' with optimal value {pf.optimal}')
 
         extra_pf_rids = pf_rids.difference(set(self.fba_model.rids))
         pf_unprotected_rids = pf_rids.difference(set(self.protected_rids))
         if len(extra_p_rids) > 0:
-            print('protected reactions that do not exist in model:', extra_p_rids)
+            print('Protected reactions that do not exist in model:', extra_p_rids)
         if len(extra_p_sids) > 0:
-            print('protected metabolites that do not exist in model:', extra_p_sids)
+            print('Protected metabolites that do not exist in model:', extra_p_sids)
         if len(extra_pf_rids) > 0:
-            print('reactions in protected functions that do not exist in model:', extra_pf_rids)
+            print('Reactions in protected functions that do not exist in model:', extra_pf_rids)
         if len(pf_unprotected_rids) > 0:
-            print('reactions in protected functions that are not protected:', pf_unprotected_rids)
+            print('Reactions in protected functions that are not protected:', pf_unprotected_rids)
 
         blocked_rids = self.reaction_types_fva(self.fba_model.rids)['blocked_rids']
         blocked_p_rids = self.nrp.protected_sids.intersection(blocked_rids)
         blocked_pf_rids = pf_rids.intersection(blocked_rids)
         if len(blocked_p_rids) > 0:
-            print('protected reactions that are blocked (no flux):', blocked_p_rids)
+            print('Protected reactions that are blocked (no flux):', blocked_p_rids)
         if len(blocked_pf_rids) > 0:
-            print('reactions in protected functions that are blocked (no flux):', blocked_pf_rids)
+            print('Reactions in protected functions that are blocked (no flux):', blocked_pf_rids)
 
     def identify_parallel_reactions(self):
         """Identify parallel reactions having same stoichiometry (considering reverse).
@@ -183,17 +232,24 @@ class ModelPruner:
         flux_max_pfs = np.zeros((n_pf, len(free_rids)))
 
         res = {}
-        processes = min(self.cpus, len(free_rids))
-        for idx, pf in enumerate(self.nrp.protected_functions.values()):
-            if processes > 100:
+        # setting up multiprocessing takes quite some resources
+        #    there is a tradeoff wrt running on one processor
+        processes = min(self.cpus, int(len(free_rids)*n_pf/self.min_mp_total_flux))
+        # processes = min(self.cpus, len(free_rids))
+
+        if processes > 1:
+            print(f'Multiprocessing with {processes} processes')
             frid2idx = {rid: idx for idx, rid in enumerate(free_rids)}
-                # place to split fva for free_rids over several
-                with mp.Pool(processes, initializer=init_worker, initargs=(self.fba_model, pf)) as pool:
-                    for res in pool.imap_unordered(worker, free_rids):
-                        rid = res['rids'][0]
-                        flux_min_pfs[idx, frid2idx[rid]] = res['min'][0]
-                        flux_max_pfs[idx, frid2idx[rid]] = res['max'][0]
+            # delete lp problem prior to pool creation (so FbaProblem can be pickled)
+            self.fba_model.delete_fba_lp()
+            with mp.Pool(processes, initializer=_init_worker,
+                         initargs=(self.fba_model, self.nrp.protected_functions)) as pool:
+                for res in pool.imap_unordered(_worker, np.array_split(free_rids, processes)):
+                    idxs = [frid2idx[rid] for rid in res['rids']]
+                    flux_min_pfs[:, idxs] = res['min_pfs']
+                    flux_max_pfs[:, idxs] = res['max_pfs']
         else:
+            for idx, pf in enumerate(self.nrp.protected_functions.values()):
                 # print(f'\nCheck FVA for pf {pf.obj_id} with {len(free_rids)} free variables')
                 res = self.fba_model.fva_pf_flux_ranges(pf, free_rids)
                 if res['success'] is True:
@@ -262,16 +318,27 @@ class ModelPruner:
 
         return reaction_types
 
-    def init_worker(self):
-        global _fba_model
-        _fba_model = copy.deepcopy(self.fba_model)
+    def prune(self, snapshot_freq=0):
+        """Prune the metabolic network
+
+        Once the Metabolic network has been loaded it can be pruned.
+        With snapshot_freq one can deterimine if (snapshot_freq > 0) and
+        how often intermediate results shall be store, from which one can resume in future
+
+        Step during Pruning
+        1. identification and removal of parallel reactions
+        2. check feasibility of the network wrt to all protected functions
+        3. in a loop till all free reactions have been consumed
+        3.1 identify types of remaining free reactions using FVA across all protected functions
+        3.2 blocked reactions (zero flux across all functions) get removed
+        3.3 essential reactions (strict positive/negative flux across conditions) get protected
+        3.4 a single candidate reaction (supporting zero flux) gets identify and removed
 
-    @staticmethod
-    def fva_single_rids(rid):
-        return _fba_model.fva_single_rid(rid)
+        :param snapshot_freq: frequency to store intermediate results (every x reactions)
+        :type snapshot_freq: int, optional (default: 0, no snapshots)
+        """
 
-    def prune(self, snapshot_freq=None):
-        print(f'initial S-matrix: {self.fba_model.s_mat.shape}; '
+        print(f'Initial S-matrix: {self.fba_model.s_mat_coo.shape}; '
               f' protected reactions: {len(self.protected_rids)}:', self.protected_rids)
         drop_rids = self.identify_parallel_reactions().difference(self.protected_rids)
         if len(drop_rids) > 0:
@@ -279,13 +346,13 @@ class ModelPruner:
             print(f'{len(drop_rids)} parallel reaction(s) dropped:', drop_rids)
 
         feasible = np.array([self.fba_model.check_pf_feasibility(pf) for pf in self.nrp.protected_functions.values()])
-        print('protected functions feasibility:', np.all(feasible))
+        print('Protected functions feasibility:', np.all(feasible))
 
         next_snapshot = 0
-        if type(snapshot_freq) is int:
+        if snapshot_freq > 0:
             next_snapshot = max(0, len(self.fba_model.rids) - snapshot_freq)
 
-        while len(self.protected_rids) < self.fba_model.s_mat.shape[1]:
+        while len(self.protected_rids) < self.fba_model.s_mat_coo.shape[1]:
             free_rids = list(set(self.fba_model.rids).difference(self.protected_rids))
             reaction_types = self.reaction_types_fva(free_rids)
             print()
@@ -299,11 +366,11 @@ class ModelPruner:
 
             # try to remove one reaction while still satisfying all protected functions
             for rid in candidate_rids:
-                # print('check feasibility of', rid)
+                # print('Check feasibility of', rid)
                 feasible = True
                 for pf in self.nrp.protected_functions.values():
                     feasible = self.fba_model.check_pf_feasibility(pf, zero_flux_rid=rid)
-                    # print(f'feasibility of {rid} is {feasible}')
+                    # print(f'Feasibility of {rid} is {feasible}')
                     if feasible is False:
                         essential_rids.add(rid)
                         break
@@ -317,7 +384,7 @@ class ModelPruner:
                 print(f"{len(essential_rids)} reaction(s) added to protected reactions:", essential_rids)
 
             if feasible is False:
-                print('no more reactions to remove')
+                print('No more reactions to remove')
 
             # store intermediate results (snapshots)
             if len(self.fba_model.rids) < next_snapshot:
@@ -328,7 +395,7 @@ class ModelPruner:
                 next_snapshot = max(0, len(self.fba_model.rids) - snapshot_freq)
                 print('snapshot data saved.')
 
-            print(f'S-matrix: {self.fba_model.s_mat.shape}; '
+            print(f'S-matrix: {self.fba_model.s_mat_coo.shape}; '
                   f'protected reactions: {len(self.protected_rids)}')
 
         # check if blocked reactions exist in the finally pruned model
@@ -336,10 +403,12 @@ class ModelPruner:
         # if len(blocked_rids) > 0:
         #     print(f'{len(blocked_rids)} blocked reaction(s) in pruned model', blocked_rids)
 
-        # remove any snapshot files
-        # for fname in [self._snapshot_sbml, self._snapshot_json]:
-        #    if os.path.exists(fname):
-        #        os.remove(fname)
+    def remove_snapshot_files(self):
+        """Remove any snapshot files still existing"""
+        for pname in [self._snapshot_sbml, self._snapshot_json]:
+            if os.path.exists(pname):
+                os.remove(pname)
+                print(f'{pname} removed')
 
     def get_reduction_status(self):
         """ Return status in model reduction process
@@ -355,6 +424,7 @@ class ModelPruner:
         return status
 
     def print_status(self):
+        """Print status of pruning process"""
 
         status = self.get_reduction_status()
         print("reactions removed/remaining/protected: "
@@ -362,11 +432,21 @@ class ModelPruner:
               f"/{len(self.protected_rids):4d}; dof: {status['dof']:3d}")
 
     def export_pruned_model(self, pruned_sbml):
+        """Export the pruned network (or snapshot) to SBML.
+
+        Component information of initial network will be carried over,
+        e.g. annotation data, groups, gene products
+
+        :param pruned_sbml: pathname of pruned network, ending with '.xml'
+        :type pruned_sbml: str
+        """
         success = self.fba_model.export_pruned_model(pruned_sbml)
         if success is True:
-            print('pruned model exported to', pruned_sbml)
+            print('Pruned model exported to', pruned_sbml)
+            # self.remove_snapshot_files()
         else:
-            print('pruned model could not be exported')
+            print('Pruned model could not be exported')
 
     def __del__(self):
+        """Delete corresponding LpProblem"""
         del self.fba_model

modelpruner/models/fba_model.py

@@ -5,6 +5,7 @@ import time
 import re
 import math
 import numpy as np
+from scipy import sparse
 import sbmlxdf
 
 from modelpruner.problems.lp_problem import LpProblem
@@ -41,12 +42,13 @@ class FbaModel:
                 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)}
-        self.s_mat = sbml_model.get_s_matrix().to_numpy()
+        # using a sparse matrix speeds up significantly pickling (factor > 350 for S-matrix)
+        # pickling required to support multiprocessing
+        # also data exchange from sbmlxdf.Model and towards LpProblem improves
+        self.s_mat_coo = sparse.coo_matrix(sbml_model.get_s_matrix(sparse=True))
 
-        # TODO: create LP Problem only once (but what is impact on deepcopy()
-        #   function to assign new LP problem
-
-        # lp problem will only be created once required. This gives chances to update e.g. flux bounds
+        # lp problem will only be created once required.
+        # this gives chances to update e.g. flux bounds and to support pickling (multiprocessing)
         self.lp = None
 
     def update_reversible_rids(self):
@@ -87,11 +89,13 @@ class FbaModel:
         :return: lists of reaction ids that are parallel
         :rtype: list with lists of str
         """
-        fwd_rev_s_mat = np.hstack((self.s_mat, -self.s_mat[:, self.reversible]))
+
+        fwd_rev_s_mat_coo = sparse.hstack([self.s_mat_coo,
+                                           -self.s_mat_coo.tocsc()[:, self.reversible]])
         fwd_rev_rids = np.hstack((self.rids, self.rids_rev))
 
         # identify groups of parallel reactions
-        sorted_idx = np.lexsort(fwd_rev_s_mat)
+        sorted_idx = np.lexsort(fwd_rev_s_mat_coo.todense())[0]
         groups = []
         consumed_idxs = set()
         for pos in range(len(sorted_idx)):
@@ -100,7 +104,7 @@ class FbaModel:
                 group = [idx]
                 for adjacent_pos in range(pos + 1, len(sorted_idx)):
                     candiate_idx = sorted_idx[adjacent_pos]
-                    if np.all(fwd_rev_s_mat[:, idx] == fwd_rev_s_mat[:, candiate_idx]):
+                    if (fwd_rev_s_mat_coo.getcol(idx) != fwd_rev_s_mat_coo.getcol(candiate_idx)).getnnz() == 0:
                         consumed_idxs.add(candiate_idx)
                         group.append(candiate_idx)
                     else:
@@ -131,14 +135,24 @@ class FbaModel:
         :returns: status information in a dictionary
         :rtype: dict
         """
-        status = {'dof': self.s_mat.shape[1] - np.linalg.matrix_rank(self.s_mat),
+        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)}
         status.update(self.full_model_shape)
         return status
 
     def export_pruned_model(self, pruned_sbml):
+        """Export the pruned network (or snapshot) to SBML.
+
+        Component information of initial network will be carried over,
+        e.g. annotation data, groups, gene products
+        History information, if available, will be updated.
+
+        Model will also be validated with respect to SBML compliance.
 
+        :param pruned_sbml: pathname of pruned network, ending with '.xml'
+        :type pruned_sbml: str
+        """
         sbml_model = sbmlxdf.Model(self.sbml_fname)
         model_dict = sbml_model.to_df()
         pruned_mdict = model_dict.copy()
@@ -147,8 +161,11 @@ class FbaModel:
 
         # update model attributes (ids, names and modification date)
         pruned_mdict['modelAttrs']['id'] += '_pruned'
+        if 'name' in pruned_mdict['modelAttrs']:
             pruned_mdict['modelAttrs']['name'] += '_pruned'
+        if 'metaid' in pruned_mdict['modelAttrs']:
             pruned_mdict['modelAttrs']['metaid'] += '_pruned'
+        if 'modifiedHistory' in pruned_mdict['modelAttrs']:
             pruned_mdict['modelAttrs']['modifiedHistory'] += '; localtime'
 
         # clean groups table
@@ -203,13 +220,23 @@ class FbaModel:
         :return: LpProblem configured for FBA
         :rtype: LpProblem
         """
-        fwd_rev_s_mat = np.hstack((self.s_mat[:], -self.s_mat[:, self.reversible]))
+        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.shape[0], 2))
-        lp = LpProblem(fwd_rev_s_mat, row_bnds, fwd_rev_col_bnds)
+        row_bnds = np.zeros((fwd_rev_s_mat_coo.shape[0], 2))
+        lp = LpProblem(fwd_rev_s_mat_coo, row_bnds, fwd_rev_col_bnds)
         return lp
 
+    def delete_fba_lp(self):
+        """Delete the corresponding linear problem.
+
+        E.g. required for multiprocessing so FbaModel can be pickled
+        """
+        if self.lp is not None:
+            del self.lp
+            self.lp = None
+
     def set_lp_objective(self, objective, direction='maximize'):
         """Set objective coefficients and direction in linear problem
 
@@ -241,7 +268,7 @@ class FbaModel:
         """
         if self.lp is None:
             self.lp = self.create_core_fba_lp()
-            print('Core LP Problem created')
+            # print('Core LP Problem created')
         for rid, bound_upd in flux_bounds.items():
             if rid in self.rids:
                 temp_bounds = self.flux_bounds[:, self.rid2idx[rid]].copy()
@@ -291,6 +318,9 @@ class FbaModel:
         :param drop_rids: reaction ids to be removed from the model
         :type drop_rids: set of str
         """
+        if self.lp is None:
+            self.lp = self.create_core_fba_lp()
+            # print('Core LP Problem created')
         assert self.lp is not None
 
         # remove variables (columns) from linear problem
@@ -304,16 +334,16 @@ class FbaModel:
         self.reversible = self.reversible[keep_rids_mask]
         self.update_reversible_rids()
         self.flux_bounds = self.flux_bounds[:, keep_rids_mask]
-        self.s_mat = self.s_mat[:, keep_rids_mask]
+        self.s_mat_coo = self.s_mat_coo.tocsc()[:, keep_rids_mask].tocoo()
 
         # drop blocked metabolites due to deleted reactions
-        keep_sids_mask = np.any(self.s_mat != 0, axis=1)
+        keep_sids_mask = (self.s_mat_coo.getnnz(axis=1) > 0)
         row_idxs0 = [idx for idx, val in enumerate(keep_sids_mask) if val == np.bool_(False)]
         if len(row_idxs0) > 0:
             self.lp.remove_rows(row_idxs0)
             self.sids = self.sids[keep_sids_mask]
             self.sid2idx = {sid: idx for idx, sid in enumerate(self.sids)}
-            self.s_mat = self.s_mat[keep_sids_mask, :]
+            self.s_mat_coo = self.s_mat_coo.tocsr()[keep_sids_mask, :].tocoo()
 
     def fba_pf_optimize(self, pf, zero_flux_rid=None):
         """FBA optimization of current problem
@@ -329,7 +359,7 @@ class FbaModel:
         """
         if self.lp is None:
             self.lp = self.create_core_fba_lp()
-            print('Core LP Problem created')
+            # print('Core LP Problem created')
         assert self.lp is not None
 
         if zero_flux_rid is not None:
@@ -430,9 +460,6 @@ class FbaModel:
         result.__dict__['lp'] = None
         return result
 
-    def __reduce__(self):
-        return self.__class__, (self.sbml_fname, )
-
     def __del__(self):
         print('deletion of fba_model')
         del self.lp

modelpruner/problems/lp_problem.py

@@ -6,8 +6,6 @@ Peter Schubert, HHU Duesseldorf, June 2022
 """
 
 import numpy as np
-import scipy
-import scipy.sparse
 import swiglpk as glpk
 
 
@@ -61,13 +59,13 @@ scaling_options = {'GM': glpk.GLP_SF_GM,
 
 class LpProblem:
 
-    def __init__(self, constr_coefs, row_bnds, col_bnds):
+    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: matrix with constraint coefficients
-        :type constr_coefs: 2D nparray with shape(nrows, ncols)
+        :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
@@ -78,7 +76,7 @@ class LpProblem:
         glp_presolve = glpk.GLP_OFF   # or glpk.GLP_OFF (default) or GLP_ON
         self.res = {}
 
-        self.nrows, self.ncols = constr_coefs.shape
+        self.nrows, self.ncols = constr_coefs_coo.shape
         assert col_bnds.shape == (2, self.ncols)
         assert row_bnds.shape == (self.nrows, 2)
 
@@ -95,15 +93,14 @@ class LpProblem:
         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)
-        constr_coefs_sparse = scipy.sparse.coo_matrix(constr_coefs)
-        n_coefs = len(constr_coefs_sparse.row)
+        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_sparse.row[i])
-            ja[1 + i] = int(1 + constr_coefs_sparse.col[i])
-            ar[1 + i] = constr_coefs_sparse.data[i]
+            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
@@ -307,7 +304,7 @@ class LpProblem:
     def __del__(self):
         """Explicitly delete the LpProblem. (can we do without?)
         """
-        print('deletion of LpProblem')
+        # print('deletion of LpProblem')
         assert self.lp is not None
         if getattr(self, 'lp'):
             glpk.glp_delete_prob(self.lp)