- implement write method and initialization from init file

 - add more tests

pyenzyme/thinlayers/TL_Copasi.py

@@ -11,8 +11,9 @@
 Copyright (c) 2021 Institute of Biochemistry and Technical Biochemistry Stuttgart
 '''
 import logging
-from typing import Union
+from typing import Union, Optional
 
+from pyenzyme import EnzymeMLDocument
 from pyenzyme.thinlayers import BaseThinLayer
 import os
 
@@ -23,18 +24,21 @@
 
 class ThinLayerCopasi(BaseThinLayer):
 
-    def __init__(self, path, outdir, measurement_ids: Union[str, list] = "all"):
+    def __init__(self, path, outdir,
+                 measurement_ids: Union[str, list] = "all",
+                 init_file: Optional[str] = None):
         """
         Initializes a new instance of the COPASI thin layer, by loading the EnzymeML file
         specified in `path` and creating a COPASI file (+ data) in `outdir`.
 
         :param path: the enzyme ml document to load
         :param outdir: the output dir
         :param measurement_ids: the measurement ids or all
+        :param init_file: optional initialization file for fit items
         """
 
         # initialize base class, let it do the reading
-        BaseThinLayer.__init__(self, path, measurement_ids)
+        BaseThinLayer.__init__(self, path, measurement_ids, init_file=init_file)
         self.name = self.enzmldoc.name
         self.working_dir = os.path.join(os.path.abspath(outdir), self.name)
         self.cps_file = os.path.join(self.working_dir, self.name + '.cps')
@@ -65,7 +69,10 @@ def __init__(self, path, outdir, measurement_ids: Union[str, list] = "all"):
         self._import_experiments()
 
         # set all parameters as fit items
-        self._set_default_items()
+        if init_file is None:
+            self._set_default_items()
+        else:
+            self._set_default_items_from_init_file()
 
         self.dm.saveModel(self.cps_file, True)
 
@@ -183,12 +190,15 @@ def get_fit_parameters(self):
                 continue
             name = obj.getObjectName() if obj.getObjectType() != "Reference" else obj.getObjectParent().getObjectName()
 
+            r = obj.getObjectAncestor('Reaction')
+            reaction_id = r.getSBMLId() if r else None
+
             result.append({
                 'name': name,
                 'start': p.getStartValue(),
                 'lower': p.getLowerBoundValue(),
                 'upper': p.getUpperBoundValue(),
-                'reaction_id': obj.getObjectAncestor('Reaction').getSBMLId()
+                'reaction_id': reaction_id
             })
         return result
 
@@ -231,6 +241,39 @@ def optimize(self):
         self.task.process(True)
         self.task.restore()
 
+    def write(self):
+        """Writes the estimated parameters to a copy of the EnzymeMLDocument"""
+
+        nu_enzmldoc = self.enzmldoc.copy(deep=True)
+
+        assert (isinstance(self.problem, COPASI.CFitProblem))
+        results = self.problem.getSolutionVariables()
+
+        logging.debug('OBJ: {0}'.format(self.problem.getSolutionValue()))
+        logging.debug('RMS: {0}'.format(self.problem.getRMS()))
+
+        for i in range(self.problem.getOptItemSize()):
+            item = self.problem.getOptItem(i)
+            obj = self.dm.getObject(item.getObjectCN())
+            if obj is None:
+                continue
+
+            name = obj.getObjectName() if obj.getObjectType() != 'Reference' else obj.getObjectParent().getObjectName()
+            value = results.get(i)
+            logging.debug(name, value)
+
+            reaction = obj.getObjectAncestor('Reaction')
+            if reaction is not None:
+                enz_reaction = nu_enzmldoc.getReaction(reaction.getSBMLId())
+                if enz_reaction:
+                    parameter = enz_reaction.model.getParameter(name)
+                    parameter.value = value
+            else:
+                p = nu_enzmldoc.global_parameters.get(name)
+                p.value = value
+
+        return nu_enzmldoc
+
     def update_enzymeml_doc(self):
         """ Updates the enzyme ml document with the values from last optimization run
 
@@ -277,6 +320,69 @@ def _set_default_items(self):
                 fit_item.setUpperBound(COPASI.CCommonName(str(1e6)))
                 fit_item.setStartValue(float(p.value))
 
+    def _set_default_items_from_init_file(self):
+        """ Use this to create a default template, when an init file was passed to the thin layer
+            and it has been already appied
+
+        :return: None
+        """
+        assert (isinstance(self.problem, COPASI.CFitProblem))
+
+        # remove old items
+        while self.problem.getOptItemSize() > 0:
+            self.problem.removeOptItem(0)
+
+        for global_param in self.global_parameters.values():
+
+            if not global_param.lower or not global_param.upper:
+                # nan values used to indicate that this should not be fitted
+                continue
+
+            mv = self.dm.getModel().getModelValue(global_param.name)
+            if not mv:
+                logging.warning("No global parameter {0} in the model".format(global_param.name))
+                continue
+
+            value = global_param.value if global_param.value else global_param.initial_value
+
+            if not value:
+                raise ValueError(
+                    f"Neither initial_value nor value given for parameter {global_param.name} in global parameters"
+                )
+
+            cn = mv.getInitialValueReference().getCN()
+            fit_item = self.problem.addFitItem(cn)
+            assert (isinstance(fit_item, COPASI.CFitItem))
+            fit_item.setLowerBound(COPASI.CCommonName(str(global_param.lower)))
+            fit_item.setUpperBound(COPASI.CCommonName(str(global_param.upper)))
+            fit_item.setStartValue(float(value))
+
+        for reaction_id, (model, _) in self.reaction_data.items():
+            r = self.sbml_id_map[reaction_id]
+            assert (isinstance(r, COPASI.CReaction))
+            for p in model.parameters:
+                if p.is_global:
+                    continue
+
+                if not p.lower or not p.upper:
+                    # nan values used to indicate that this should not be fitted
+                    continue
+
+                value = p.value if p.value else p.initial_value
+
+                if not value:
+                    raise ValueError(
+                        f"Neither initial_value nor value given for parameter {p.name} in reaction {reaction_id}"
+                    )
+
+                obj = r.getParameterObjects(p.name)[0].getObject(COPASI.CCommonName('Reference=Value'))
+                cn = obj.getCN()
+                fit_item = self.problem.addFitItem(cn)
+                assert (isinstance(fit_item, COPASI.CFitItem))
+                fit_item.setLowerBound(COPASI.CCommonName(str(p.lower)))
+                fit_item.setUpperBound(COPASI.CCommonName(str(p.upper)))
+                fit_item.setStartValue(float(value))
+
 
 if __name__ == '__main__':
     this_dir = os.path.dirname(__file__)

tests/test_tl_copasi.py

@@ -9,19 +9,19 @@
 if not os.path.exists(temp_dir):
     os.mkdir(temp_dir)
 
-example_file = os.path.join(this_dir, '..', 'examples', 'ThinLayers', 'COPASI', '3IZNOK_Simulated.omex')
-
 
 class TestTlCopasi(unittest.TestCase):
 
+    example_file = os.path.join(this_dir, '..', 'examples', 'ThinLayers', 'COPASI', '3IZNOK_Simulated.omex')
+
     def test_copasi_version(self):
         self.assertGreaterEqual(int(COPASI.CVersion.VERSION.getVersionDevel()), 214, "Need newer COPASI version")
 
     def test_example_file_exists(self):
-        self.assertTrue(os.path.exists(example_file))
+        self.assertTrue(os.path.exists(self.example_file))
 
     def test_example(self):
-        thin_layer = ThinLayerCopasi(path=example_file, outdir=temp_dir)
+        thin_layer = ThinLayerCopasi(path=self.example_file, outdir=temp_dir)
         self.assertEqual(thin_layer.reaction_data['r0'][0].parameters[0].name, 'k_cat')
         self.assertEqual(thin_layer.reaction_data['r0'][0].parameters[0].value, 0.015)
         self.assertEqual(thin_layer.reaction_data['r0'][0].parameters[1].name, 'k_m')
@@ -32,3 +32,54 @@ def test_example(self):
         self.assertNotEqual(thin_layer.reaction_data['r0'][0].parameters[0].value, 0.015)
         self.assertEqual(thin_layer.reaction_data['r0'][0].parameters[1].name, 'k_m')
         self.assertNotEqual(thin_layer.reaction_data['r0'][0].parameters[1].value, 0.01)
+
+
+class TestTlCopasiENO(unittest.TestCase):
+    example_file = os.path.join(this_dir, '..', 'examples', 'ThinLayers', 'COPASI', 'PGM-ENO.omex')
+
+    def test_copasi_version(self):
+        self.assertGreaterEqual(int(COPASI.CVersion.VERSION.getVersionDevel()), 214, "Need newer COPASI version")
+
+    def test_example_file_exists(self):
+        self.assertTrue(os.path.exists(self.example_file))
+
+    def test_example(self):
+        thin_layer = ThinLayerCopasi(path=self.example_file, outdir=temp_dir)
+        fit_items = thin_layer.get_fit_parameters()
+        initial_values = [val['start'] for val in fit_items]
+        thin_layer.optimize()
+        new_values = [val for val in thin_layer.problem.getSolutionVariables()]
+        self.assertNotEqual(initial_values, new_values)
+        thin_layer.update_enzymeml_doc()
+
+
+class TestModel4(unittest.TestCase):
+    example_file = os.path.join(this_dir, '..', 'examples', 'ThinLayers', 'COPASI', 'Model_4.omex')
+    init_file = os.path.join(this_dir, '..', 'examples', 'ThinLayers', 'COPASI', 'Model_4_init.yaml')
+
+    def test_copasi_version(self):
+        self.assertGreaterEqual(int(COPASI.CVersion.VERSION.getVersionDevel()), 214, "Need newer COPASI version")
+
+    def test_example_file_exists(self):
+        self.assertTrue(os.path.exists(self.example_file))
+        self.assertTrue(os.path.exists(self.init_file))
+
+    def test_example(self):
+        thin_layer = ThinLayerCopasi(path=self.example_file, outdir=temp_dir, init_file=self.init_file)
+        fit_items = thin_layer.get_fit_parameters()
+        initial_values = [val['start'] for val in fit_items]
+        thin_layer.optimize()
+        new_values = [val for val in thin_layer.problem.getSolutionVariables()]
+        self.assertNotEqual(initial_values, new_values)
+        new_doc = thin_layer.write()
+
+        self.assertIsNotNone(new_doc)
+
+        #assert (isinstance(thin_layer.task, COPASI.CFitTask))
+        #thin_layer.task.setMethodType(COPASI.CTaskEnum.Method_LevenbergMarquardt)
+        #thin_layer.optimize()
+        #new_values2 = [val for val in thin_layer.problem.getSolutionVariables()]
+        #thin_layer.task.setMethodType(COPASI.CTaskEnum.Method_LevenbergMarquardt)
+        #thin_layer.optimize()
+        #new_values3 = [val for val in thin_layer.problem.getSolutionVariables()]
+        #self.assertNotEqual(initial_values, new_values)