Eliminated the save_path feature since it conflicted with usage in co…

…lab.

README.md

@@ -20,10 +20,13 @@ ctl.__version__
 ```
 
 ## Version History
-* 1.1.04
+* 1.2.1 5/23/2024
+    * Eliminated the save_path feature since it conflicted with usage in colab.
+* 1.2.0 5/22/2024
     * Implemented differential control
     * Extended noise model to include lognormal distribution, offset, and slope
     * Implemented several examples, some of which are based on student projects
+    * Two algorithms for fitting. gpz fits the transfer function in order by gain (g), poles (p), zeroes (z). poly fits a polynomial in s to the simulations.
 * 1.1.03 1/20/2024
     * Better error checking
     * API uses parameter names kI, kP, kF

docs/source/index.rst

@@ -6,22 +6,30 @@
 ``controlSBML``: Control engineering with SBML models
 =============================================================
 
-``controlSBML`` provides a bridge between SBML models and the
+``controlSBML`` is a python package for doing control analysis and design in
+systems biology
+where the open loop system is a model described by the
+community standard Systems Biology Markup Language
+(SBML).
+SBML models may be specified as strings, files, or URLs. They may be in the
+XML representation or in Antimony.
+
+At present, ``controlSBML`` only supports the development of SISO (single input, single output)
+control loops. The input can be a chemical species or an SMBL parameter (e.g.,
+a kinetic constant). The output can be a chemical species or a reaction flux.
+The package supports the following workflow:
+
+* access the SBML model
+* analyze controllability (the ability of an input to control the output)
+* system identification (find a transfer function that fits the input-output response)
+* control design for proportional-integral-differential (PID) control with a filter
+* analyze the effect of perturbations (either disturbances to the control input to the open loop system or noise in the output measured from the open loop system)
+
+``controlSBML`` leverages the
 `CalTech control systems library <https://python-control.readthedocs.io/en/latest/intro.html>`_.
-At the lowest level, this bridge wraps an SBML model as an object known
-to the ``control`` package so that its tools can be
-applied to SBML models.
-The approach is:
+The control package is used
+to do classical control analysis such as bode plots and root locus plots.
 
-1. Construct a ``controlSBML.ControlSBML`` object for the SBML model.
-
-2. From this object, obtain objects used in the ``control`` package.
-
-   a. ``control.NonlinearIOSystem`` objects wrap the entire SBML model to construct computational testbeds for closed loop systems.
-
-   b. ``control.StateSpace`` objects are linear approximations to the SBML model for control analysis and design using transfer functions.
-
-3. Do control analysis and design using the ``control`` package.
 
 Below, we provide an overview of the features of ``controlSBML``.
 Readers who are less familiar with control engineering

pyproject.toml

@@ -1,6 +1,6 @@
 [project]
 name = "controlSBML"
-version = "1.1.04"
+version = "1.2.1"
 description = "Control analysis of SBML models"
 authors = [{ name = "Joseph Hellerstein", email = "joseph.hellerstein@gmail.com" }]
 license = { file = "LICENSE" }

src/controlSBML/_version.py

@@ -11,5 +11,5 @@
         data = tomli.load(f)
     __version__ = data["project"]["version"]
 except FileNotFoundError:
-    __version__ = "1.1.04"
+    __version__ = "1.2.1"
 

src/controlSBML/control_sbml.py

@@ -132,8 +132,6 @@
 TIMES_KEYS = list(TIMES_DCT.keys())
 OPTION_KEYS = list(OPTION_DCT.keys())
 SIMULATION_KEYS = list(SIMULATION_DCT.keys())
-#
-SAVE_PATH = os.path.join(cn.DATA_DIR, "control_sbml.csv")
 
 # Returned results
 DesignResult = namedtuple("DesignResult", ["timeseries", "antimony_builder", "designs"])
@@ -157,7 +155,6 @@ def __init__(self, model_reference:str,
                  setpoint:Optional[float]=OPTION_DCT[cn.O_SETPOINT],
                  sign:Optional[int]=OPTION_DCT[cn.O_SIGN],
                  times:Optional[np.ndarray[float]]=None,
-                 save_path:Optional[str]=SAVE_PATH,
                  **kwargs):
         """
         model_reference: str
@@ -169,7 +166,6 @@ def __init__(self, model_reference:str,
         is_steady_state: bool
         noise_spec: NoiseSpec
         disturbance_spec: DisturbanceSpec
-        save_path: str (path to file where results are saved after a design)
         kwargs: dict with options listed below. These are the default options used unless overridden.
             Plot options:
                 ax: axis for plot
@@ -231,7 +227,6 @@ def __init__(self, model_reference:str,
         self.is_steady_state = is_steady_state
         self.noise_spec = noise_spec
         self.disturbance_spec = disturbance_spec
-        self.save_path = save_path
         self.fitter_method = OPTION_DCT[cn.O_FITTER_METHOD] if fitter_method is None else fitter_method
         self.setpoint = OPTION_DCT[cn.O_SETPOINT] if setpoint is None else setpoint
         self.sign = OPTION_DCT[cn.O_SIGN] if sign is None else sign
@@ -252,8 +247,7 @@ def copy(self):
                 fitter_method=self.fitter_method,
                 setpoint=self.setpoint,
                 sign=self.sign,
-                times=self.times,
-                save_path=self.save_path)
+                times=self.times)
         for key in cn.CONTROL_PARAMETERS:
             setattr(ctlsb, key, getattr(self, key))
         return ctlsb
@@ -686,8 +680,6 @@ def plotGridDesign(self,
                 timeseries: Timeseries
                 antimony_builder: AntimonyBuilder
         """
-        save_path = None   # Disable "save_path" feature
-        #
         new_keys = list(PLOT_KEYS)
         self._checkKwargs(new_keys, **kwargs)
         option_dct = self.getOptions(sign=sign, setpoint=setpoint, times=times, selections=selections,
@@ -699,11 +691,6 @@ def plotGridDesign(self,
         selections = option_dct[cn.O_SELECTIONS]
         num_process = option_dct[cn.O_NUM_PROCESS]
         num_restart = option_dct[cn.O_NUM_RESTART]
-        if save_path is not None:
-            if len(save_path) == 0:
-                save_path = self.save_path
-            if os.path.isfile(save_path):
-                os.remove(save_path)
         plot_dct = util.subsetDct(option_dct, PLOT_KEYS)
         # Process the request
         designer = SISOClosedLoopDesigner(self._sbml_system, self.getOpenLoopTransferFunction(),
@@ -712,8 +699,7 @@ def plotGridDesign(self,
                 sign=sign,
                 setpoint=setpoint,
                 input_name=self.input_name,
-                output_name=self.output_name,
-                save_path=self.save_path)
+                output_name=self.output_name)
         # Translate axis names
         designs = designer.designAlongGrid(grid, is_greedy=self.is_greedy, num_process=num_process,  # type: ignore
                                  num_restart=num_restart)    # type: ignore
@@ -847,8 +833,7 @@ def setSpec(key, val):
                 sign=sign,
                 setpoint=setpoint,
                 input_name=self.input_name,
-                output_name=self.output_name,
-                save_path=self.save_path)
+                output_name=self.output_name)
         designs = designer.design(
             kP_spec=kP_spec,
             kI_spec=kI_spec,
@@ -893,21 +878,16 @@ def setSpec(key, val):
                 **plot_dct)
         return DesignResult(timeseries=response_ts, antimony_builder=antimony_builder, designs=designs)
 
-    def _plotDesignResult(self, save_path:Optional[str]=None, **kwargs):
+    def _plotDesignResult(self, **kwargs):
         """
         Plots the mean squared error of all points searched in the last design.
 
         Args:
-            save_path: str (path to CSV file where design results are saved)
             kwargs: dict (plot options)
             AntimonyBuilder
         """
-        valids = ["save_path"]
-        valids = list(set(valids).union(PLOT_KEYS))
         self._checkKwargs(**kwargs)
-        if save_path is None:
-            save_path = self.save_path
-        designer = SISOClosedLoopDesigner(self._sbml_system, self.getOpenLoopTransferFunction(), save_path=save_path)
+        designer = SISOClosedLoopDesigner(self._sbml_system, self.getOpenLoopTransferFunction())
         if designer.design_result_df is None:
             msgs.warn("No design found!")
             return None, None

src/controlSBML/siso_closed_loop_designer.py

@@ -70,20 +70,17 @@ def _calculateClosedLoopTransferFunction(open_loop_transfer_function=None, kP=No
 class SISOClosedLoopDesigner(object):
 
     def __init__(self, system, open_loop_transfer_function, times=None, setpoint=SETPOINT, is_steady_state=False,
-                 is_history=True, sign=-1, input_name=None, output_name=None,
-                 save_path=None):
+                 sign=-1, input_name=None, output_name=None):
         """
         Constructs a SISOClosedLoopDesigner object. If the system has more than one input (output) and not input (output) is specified),
         then the first input (output) is used.
         Args:
             sbml_system: SBMLSystem
             sys_tf: control.TransferFunction (open loop system)
             is_steady_state: bool (if True, then the steady state is used)
-            is_history: bool (if True, then history is maintained)
             sign: int (-1: negative feedback; +1: positive feedback)
             input_name: str (name of input species)
             output_name: str (name of output species)
-            save_path: str (path to save the results or to an existing result)
         """
         self.system = system
         self.open_loop_transfer_function = open_loop_transfer_function
@@ -99,7 +96,6 @@ def __init__(self, system, open_loop_transfer_function, times=None, setpoint=SET
             output_name = self.system.output_names[0]
         self.input_name = input_name
         self.output_name = output_name
-        self.save_path = save_path
         # Calculated properties
         self.start_time = self.times[0]
         self.end_time = self.times[-1]
@@ -113,7 +109,7 @@ def __init__(self, system, open_loop_transfer_function, times=None, setpoint=SET
         self.closed_loop_system = None
         self.residual_mse = None
         self.minimizer_result = None
-        self._design_result_df = None   # Calculated in property
+        self.design_result_df = None   # Calculated in property
         #
         self._initializeDesigner()
         self.siso = None # SISOClosedLoopSystem
@@ -358,26 +354,9 @@ def getValue(val):
             self.kD = getValue(sorted_mean_df.loc[0, CP_kD])
         if CP_kF in sorted_mean_df.columns:
             self.kF = getValue(sorted_mean_df.loc[0, CP_kF])
-        # Save the results
-        if self.save_path is not None:
-            sorted_mean_df.to_csv(self.save_path, index=False)
-        self._design_result_df = sorted_mean_df
+        self.design_result_df = sorted_mean_df
         return DesignResult(dataframe=result_df, max_count=max_count)
 
-    @property
-    def design_result_df(self):
-        """
-        Returns:
-            pd.DataFrame
-                columns: kP, kI, kD, kF, mse
-        """
-        if self._design_result_df is None:
-            if (self.save_path is not None) and (os.path.isfile(self.save_path)):
-                self._design_result_df = pd.read_csv(self.save_path)
-            else:
-                return None
-        return self._design_result_df
-
     def simulateTransferFunction(self, transfer_function=None, period=None)->tuple[np.array, np.array]:
         """
         Simulates the closed loop transfer function based on the parameters of the object.

tests/test_control_sbml.py

@@ -11,14 +11,14 @@
 import matplotlib.pyplot as plt
 import numpy as np
 import os
+from typing import List
 import unittest
 
 
 IGNORE_TEST = False
 IS_PLOT = False
 TIMES = cn.TIMES
 FIGSIZE = (5, 5)
-SAVE1_PATH = os.path.join(cn.TEST_DIR, "control_sbml_save_path.csv")
 LINEAR_MDL = """
 model *main_model()
 // Illustrate Antimony File
@@ -43,9 +43,7 @@
 """
 SPECIES_NAMES = ["S1", "S2", "S3", "S4"]
 CTLSB = ControlSBML(LINEAR_MDL, input_name="S1", output_name="S3")
-CSV_FILE1 = os.path.join(cn.TEST_DIR, "test_control_sbml1.csv")
-CSV_FILE2 = os.path.join(cn.TEST_DIR, "test_control_sbml2.csv")
-REMOVE_FILES = [CSV_FILE1, CSV_FILE2]
+REMOVE_FILES:List[str] = []
 
 
 #############################
@@ -182,16 +180,6 @@ def testPlotGridDesign(self):
         grid.addAxis("kD", num_coordinate=3)
         _ = ctlsb.plotGridDesign(grid, is_plot=IS_PLOT)
 
-    def testPlotDesign1(self):
-        if IGNORE_TEST:
-            return
-        setpoint = 5
-        ctlsb = ControlSBML(LINEAR_MDL, final_value=10, input_name="S1", output_name="S3", save_path=CSV_FILE1)
-        _ = ctlsb.plotDesign(setpoint=setpoint, kP_spec=True, kI_spec=False, is_plot=IS_PLOT,
-                                            min_parameter_value=0.001, max_parameter_value=10, num_restart=1,
-                                            num_coordinate=2)
-        self.assertTrue(os.path.isfile(CSV_FILE1))
-
     def testEqualsCopy(self):
         if IGNORE_TEST:
             return
@@ -242,7 +230,7 @@ def testPlotDesignResult(self):
         if IGNORE_TEST:
             return
         setpoint = 5
-        ctlsb = ControlSBML(LINEAR_MDL, final_value=10, input_name="S1", output_name="S3", save_path=CSV_FILE2)
+        ctlsb = ControlSBML(LINEAR_MDL, final_value=10, input_name="S1", output_name="S3")
         _ = ctlsb.plotDesign(setpoint=setpoint, kP_spec=True, kI_spec=True, kD_spec=True, is_plot=False,
                                             min_parameter_value=0.001, max_parameter_value=10, num_restart=1,
                                             num_coordinate=4)
@@ -314,7 +302,7 @@ def testBug3(self):
         path = util.getModelPath("Varusai2018")
         INPUT_NAME = "pIRS"
         OUTPUT_NAME = "pmTORC1"
-        ctlsb = ControlSBML(path, save_path=SAVE1_PATH, input_name=INPUT_NAME, output_name=OUTPUT_NAME)
+        ctlsb = ControlSBML(path, input_name=INPUT_NAME, output_name=OUTPUT_NAME)
         #
         grid = ctlsb.getGrid(kP_spec=True, kI_spec=False, kD_spec=True, num_coordinate=10, is_random=False)
         axis = grid.getAxis("kP")
@@ -330,7 +318,7 @@ def testBug4(self):
         path = util.getModelPath("Varusai2018")
         INPUT_NAME = "pIRS"
         OUTPUT_NAME = "pmTORC1"
-        ctlsb = ControlSBML(path, save_path=SAVE1_PATH, input_name=INPUT_NAME, output_name=OUTPUT_NAME)
+        ctlsb = ControlSBML(path, input_name=INPUT_NAME, output_name=OUTPUT_NAME)
         #
         grid = ctlsb.getGrid(kP_spec=True, kI_spec=False, num_coordinate=40, is_random=False)
         axis = grid.getAxis("kP")
@@ -464,7 +452,7 @@ def testBug13(self):
         INPUT_NAME = "E"
         OUTPUT_NAME = "R"
         CTLSB = ControlSBML(path, figsize=(5, 5), times=np.linspace(0, 10, 100), markers=False, is_fixed_input_species=True,
-                   save_path="data.csv", input_name=INPUT_NAME, output_name=OUTPUT_NAME)  # Specify default value of options
+                   input_name=INPUT_NAME, output_name=OUTPUT_NAME)  # Specify default value of options
         TIMES = np.linspace(0, 10**5, 10**6)
         _ = CTLSB.plotDesign(setpoint=0.1, kP_spec=1, kI_spec=0.1, times=TIMES, num_restart=1, is_plot=IS_PLOT,
                              num_coordinate=2)

tests/test_siso_closed_loop_designer.py

@@ -14,6 +14,7 @@
 import os
 import pandas as pd  # type: ignore
 import sympy # type: ignore
+from typing import List
 import unittest
 
 IGNORE_TEST = False
@@ -73,8 +74,7 @@
 TIMES = np.linspace(0, 20, 200)
 PARAMETER_DCT = {p: n+1 for n, p in enumerate(CONTROL_PARAMETERS)}
 SETPOINT = 3
-SAVE_PATH = os.path.join(cn.TEST_DIR, "siso_closed_loop_designer.csv")
-REMOVE_FILES = [SAVE_PATH]
+REMOVE_FILES:List[str] = []
 CONTROL_PARAMETER_SPECS = ["kP_spec", "kI_spec", "kD_spec", "kF_spec"]
 #if False:
 #    # Required to construct the transfer function
@@ -113,8 +113,8 @@ def init(self):
             return
         self.sys_tf = copy.deepcopy(TRANSFER_FUNCTION)
         self.system = copy.deepcopy(SYSTEM)
-        self.designer = cld.SISOClosedLoopDesigner(self.system, self.sys_tf, times=TIMES, setpoint=SETPOINT,
-                                                   save_path=SAVE_PATH)
+        self.designer = cld.SISOClosedLoopDesigner(self.system, self.sys_tf, times=TIMES, setpoint=SETPOINT)
+
 
     def testGetSet(self):
         if IGNORE_TEST:
@@ -183,15 +183,11 @@ def testPlot(self):
         self.designer.set(kP=10, kI=5, kD=2)
         self.designer.plot(is_plot=IS_PLOT)
 
-    def makeDesigner(self, end_time=200, is_save_path=True):
+    def makeDesigner(self, end_time=200):
         times = np.linspace(0, end_time, 10*end_time)
         system = copy.deepcopy(SYSTEM)
         transfer_function = copy.deepcopy(TRANSFER_FUNCTION)
-        if is_save_path:
-            save_path = SAVE_PATH
-        else:
-            save_path = None
-        designer = cld.SISOClosedLoopDesigner(system, transfer_function, times=times, save_path=save_path)
+        designer = cld.SISOClosedLoopDesigner(system, transfer_function, times=times)
         return designer
 
     def testPlot2(self):
@@ -245,7 +241,7 @@ def testPlotDesignResult(self):
         # Plots a previously computed result
         if IGNORE_TEST:
             return
-        designer = self.makeDesigner(is_save_path=False)
+        designer = self.makeDesigner()
         def test(parameter_names):
             dct = {}
             for spec in CONTROL_PARAMETER_SPECS:
@@ -266,7 +262,7 @@ def testPlotDesignResult2(self):
         if IGNORE_TEST:
             return
         def test(parameter_names):
-            designer = self.makeDesigner(is_save_path=False)
+            designer = self.makeDesigner()
             dct = {}
             for spec in CONTROL_PARAMETER_SPECS:
                 if spec in parameter_names: