Fix bug with times; handline simulation of transfer functions with in…

…itial state

examples/cancer-vitamins.ipynb

@@ -410,7 +410,7 @@
    ],
    "source": [
     "_ = CTLSB.plotDesign(setpoint=2, kP_spec=0.2, kI_spec=0.1, min_parameter_value=1,\n",
-    "                max_parameter_value=100, num_restart=1 )"
+    "                max_parameter_value=100)"
    ]
   },
   {

src/controlSBML/point_evaluator.py

@@ -9,7 +9,6 @@
 
 
 ##################################################################
-# FIXME: Need to specify the sign 
 class PointEvaluator(Evaluator):
     # Evaluates a point in the design space
     def __init__(self, sbml_system:SBMLSystem, input_name:str, output_name:str, setpoint:float, 

src/controlSBML/siso_transfer_function_builder.py

@@ -274,8 +274,9 @@ def adjustArray(arr):
         #
         _, y_arr = fitter.simulateTransferFunction(fitter.transfer_function)
         df = new_data_ts.copy()
+        times = new_data_ts.index/cn.MS_IN_SEC
         df[cn.O_PREDICTED] = y_arr
-        output_ts = ctl.Timeseries(df)
+        output_ts = ctl.Timeseries(df, times=times)
         residuals = output_ts[self.output_name] - output_ts[cn.O_PREDICTED].values
         rms_residuals = np.sqrt(np.mean(residuals**2))
         fitter_result = cn.FitterResult(

src/controlSBML/siso_transfer_function_fitter.py

@@ -10,8 +10,8 @@
 import control # type: ignore
 import numpy as np
 import pandas as pd  # type: ignore
-import scipy  # type: ignore
-from typing import Optional, Tuple, Union, Iterator
+from typing import Optional, Tuple, Union
+import warnings
 
 
 MIN_PARAMETER_VALUE = -1e6
@@ -48,7 +48,7 @@ def __init__(self, timeseries:Timeseries, num_zero:Optional[int]=0, num_pole:Opt
         max_zero_value: float
         """
         if num_zero > num_pole:  # type: ignore
-            raise ValueError("num_zero must be less than or equal to num_pole.")
+            raise ValueError("num_zero cannot be larger thatn num_pole.")
         self.timeseries = timeseries
         self.input_name, self.output_name = self._extractNames(timeseries, input_name, output_name)
         self.in_arr = timeseries[self.input_name].values
@@ -124,6 +124,11 @@ def simulateTransferFunction(self,
         transfer function assumes 0 initial conditions, a fit is first done by adjusting for the mean of
         the input and the output.
 
+        The simulation takes the mean of the input and output as the initial conditions for the state space model.
+        The state representation is the order of derviatives of the output, with the highest order derivative
+        being the 0th value of the initial state vector, and the lowest order derivative being the last value of
+        the initial state vector.
+
         Parameters
         ----------
         transfer_function: control.TransferFunction
@@ -134,12 +139,21 @@ def simulateTransferFunction(self,
         np.ndarray - times
         np.ndarray - output
         """
+        # FIXME: Do I need to consider initial state? The output starts at 0. Calculate state in terms of input?
         if transfer_function is None:
             transfer_function = self.transfer_function
         in_arr = self.in_arr - np.mean(self.in_arr)
-        result_input = control.forced_response(transfer_function, T=self.times, U=in_arr)
+        out_mean = np.mean(self.out_arr)
+        out_arr = self.out_arr - out_mean
+        # Specify initial state for the transfer function, just the initial output
+        X0 = np.zeros(len(transfer_function.den[0][0])-1)  # type: ignore
+        X0[-1] = out_arr[0]  # state position for undifferentiated output
+        # Do the simulation
+        warnings.filterwarnings("ignore")
+        result_input = control.forced_response(transfer_function, T=self.times, U=in_arr, X0=X0)
+        warnings.filterwarnings("default")
         y_arr_output = np.reshape(result_input.y, (self.length,)) 
-        y_arr_output = y_arr_output + np.mean(self.out_arr)
+        y_arr_output = y_arr_output + out_mean
         return self.times, y_arr_output
 
     def _calculateNormalizedMSE(self, transfer_function:control.TransferFunction)->float:
@@ -153,7 +167,7 @@ def _calculateNormalizedMSE(self, transfer_function:control.TransferFunction)->f
 
         Returns
         -------
-        float (root of the mean square of residuals)
+        float (root of the mean square of residuals divided by the standard deviation of the output)
         """
         _, y_arr = self.simulateTransferFunction(transfer_function)
         residuals = self.out_arr - y_arr

src/controlSBML/timeseries.py

@@ -17,7 +17,7 @@
 import controlSBML.constants as cn
 
 import numpy as np
-import pandas as pd
+import pandas as pd # type: ignore
 
 ############# FUNCTIONS ###############
 def findCommonIndices(index1, index2):

tests/test_control_sbml.py

@@ -13,8 +13,8 @@
 import unittest
 
 
-IGNORE_TEST = False
-IS_PLOT = False
+IGNORE_TEST = True
+IS_PLOT = True
 TIMES = cn.TIMES
 FIGSIZE = (5, 5)
 SAVE1_PATH = os.path.join(cn.TEST_DIR, "control_sbml_save_path.csv")
@@ -414,6 +414,17 @@ def testBug10(self):
                                   is_plot=IS_PLOT)
         self.assertEqual(result.designs.dataframe.loc[0, cn.REASON], cn.DESIGN_RESULT_SUCCESS)
 
+    def testBug11(self):
+        # Bogus initial transient on fit
+        #if IGNORE_TEST:
+        #    return
+        URL = "https://www.ebi.ac.uk/biomodels/services/download/get-files/MODEL1809060006/5/Tsai2014.xml"
+        CTLSB = ControlSBML(URL, times=np.linspace(0, 100, 1000))
+        CTLSB.setSystem(input_name="Plx1_active", output_name="APC_C_active")
+        result = CTLSB.plotTransferFunctionFit(fit_start_time=20, final_value=1.0,
+                                               is_plot=IS_PLOT)
+        import pdb; pdb.set_trace()
+