Implemented extended noise model

src/controlSBML/__init__.py

@@ -6,5 +6,5 @@
 from controlSBML.staircase import Staircase
 from controlSBML.make_roadrunner import makeRoadrunner
 # Lesser used
-from controlSBML.constants import LegendSpec
-import controlSBML.constants as constants
+from controlSBML.constants import LegendSpec, DisturbanceSpec, NoiseSpec
+import controlSBML.constants as constants

src/controlSBML/antimony_builder.py

@@ -269,7 +269,7 @@ def makeAdditionStatement(self, *pargs, is_assignment=True):
                     statement += argument_str
         self.addStatement(statement) 
 
-    def makeSinusoidSignal(self, amplitude, frequency, is_offset_amplitude=True, prefix="sinusoid", suffix=""):
+    def oldmakeSinusoidSignal(self, amplitude, frequency, is_offset_amplitude=True, prefix="sinusoid", suffix=""):
         """
         Makes a sinusoid. The created variable is prefix + suffix + "_ot"
         Prefix is used to scope within a control loop. Suffix is used to scope between control loops.
@@ -291,6 +291,33 @@ def makeSinusoidSignal(self, amplitude, frequency, is_offset_amplitude=True, pre
             statement += " + %f" % (amplitude)
         self.addStatement(statement) 
         return name
+
+    def makeNoiseElement(self, noise:cn.NoiseSpec, prefix:str="sinusoid", suffix:str="")->str:
+        """
+        Makes a sinusoid with random noise and a ramp. The created variable is prefix + suffix + "_ot"
+        Prefix is used to scope within a control loop. Suffix is used to scope between control loops.
+
+        Args:
+            noise: cn.NoiseSpec (noise specification)
+            prefix: str (beginning of the name)
+            suffix: str (ending of the name)
+        Returns:
+            str: name of the sinusoid
+        """
+        self.addStatement("")
+        self.makeComment("Make sinusoid: %s" % str(noise))
+        name = prefix +  suffix +  OT
+        statement = "%s := 0" % (name)
+        if not np.isclose(noise.sine_amp, 0):
+            statement += " + %f*sin(2*pi*%f*time)" % (noise.sine_amp, noise.sine_freq)
+        if not np.isclose(noise.random_mag, 0):
+            statement += " + %f*lognormal(0, %f)" % (noise.random_mag, noise.random_std)
+        if not np.isclose(noise.offset, 0):
+            statement += " + %f" % noise.offset
+        if not np.isclose(noise.slope, 0):
+            statement += " + %f*time" % noise.slope
+        self.addStatement(statement) 
+        return name
 
     def _makeInputOutputNames(self, prefix, suffix):
         base_name = prefix + suffix
@@ -429,7 +456,7 @@ def makePIDControllerElement(self,
         return name_in, name_ot
 
     def makeSISOClosedLoopSystem(self, input_name, output_name, kP=None, kI=None, kD=None, kF=None, setpoint=0,
-                           noise_amplitude=0, noise_frequency=20, disturbance_ampliude=0, disturbance_frequency=20,
+                           noise_spec=cn.NoiseSpec(), disturbance_spec=cn.DisturbanceSpec(),
                            initial_output_value=None, sign=-1):
         """
         Args:
@@ -440,10 +467,8 @@ def makeSISOClosedLoopSystem(self, input_name, output_name, kP=None, kI=None, kD
             kD: float
             kF: float
             setpoint: float (setpoint)
-            noise_amplitude: float (Amplitude of the additions to the output)
-            noise_frequency: float (Frequency of the additions to the output)
-            disturbance_amplitude: float (Amplitude of the disturbance)
-            disturbance_frequency: float (Frequency of the disturbance)
+            noise_spec: cn.NoiseSpec
+            disturbance_spec: cn.DisturbanceSpec
             initial_input_value: float (initial value of the input)
         """
         self.addStatement("")
@@ -457,8 +482,8 @@ def makeSISOClosedLoopSystem(self, input_name, output_name, kP=None, kI=None, kD
             initial_output_value = self.roadrunner[output_name]
         self.makeAdditionStatement(output_name, initial_output_value, is_assignment=False)
         # Make the elements of the closed loop
-        noise_ot = self.makeSinusoidSignal(noise_amplitude, noise_frequency, prefix="noise", suffix=suffix)
-        disturbance_ot = self.makeSinusoidSignal(disturbance_ampliude, disturbance_frequency, prefix="disturbance", suffix=suffix)
+        noise_ot = self.makeNoiseElement(noise_spec, prefix="noise", suffix=suffix)
+        disturbance_ot = self.makeNoiseElement(disturbance_spec, prefix="disturbance", suffix=suffix)
         filter_in, filter_ot, filter_calculation = self.makeFilterElement(kF, prefix="filter", suffix=suffix)
         controller_in, controller_ot = self.makePIDControllerElement(output_name,
               filter_calculation=filter_calculation,

src/controlSBML/constants.py

@@ -7,9 +7,41 @@
 import os
 
 
-DEFAULT_NUM_STEP = 5
-
+##########################
 # Classes
+##########################
+class NoiseSpec(object):
+    # Specifications for a noise model
+    def __init__(self, sine_amp=0, sine_freq=0, random_mag=0, random_std=0, offset=None, slope=0):
+        """_summary_
+
+        Args:
+            sine_amp (int, optional): Amplitude of the sine wave. Defaults to 0.
+            sine_freq (int, optional): Frequency of the sine wave. Defaults to 0.
+            random_mag (int, optional): Magnitude scalinging of the log-normal random variable. Defaults to 0.
+            random_std (int, optional): Standard deviation of the normal for the log-normal random variable added to the function. Defaults to 0.
+            offset (int, optional): Constant offset. Defaults to  sine_amp.
+            slope (int, optional): Slope w.r.t. time. Defaults to 0.
+        """
+        self.sine_amp = sine_amp
+        self.sine_freq = sine_freq
+        self.random_mag = random_mag
+        self.random_std = random_std
+        if offset is None:
+            offset = sine_amp
+        self.offset = offset
+        self.slope = slope
+
+    def __str__(self):
+        result = f"NoiseSpec(sine_amp={self.sine_amp}, sine_freq={self.sine_freq}, random_mag={self.random_mag}"
+        result += f" random_std={self.random_std}, dc_gain={self.offset}, slope={self.slope})"
+        return result
+
+class DisturbanceSpec(NoiseSpec):
+    # Specifications for a disturbance model
+    pass
+
+
 class LegendSpec():
 
     def __init__(self, names, crd=(1.15, 1), loc="upper right"):

src/controlSBML/control_sbml.py

@@ -198,6 +198,10 @@ def __init__(self, model_reference:str,
                 num_step: int (number of steps in staircase; default: 5)
             Miscellaneous options
                 times: list-float (times of simulation; default: np.linspace(0, 5, 50))
+                disturbance: NoiseSpec (specification of disturbance added to control input)
+                       sine_amp, sine_freq, rand_mag, rand_std, dc_gain, slope
+                noise: NoiseSpec (specification of noise added to measured output)
+                       sine_amp, sine_freq, rand_mag, rand_std, dc_gain, slope
 
         """
         self._checkKwargs(OPTION_KEYS, **kwargs)

tests/test_antimony_builder.py

@@ -1,7 +1,6 @@
 import controlSBML.constants as cn # type: ignore
 from controlSBML import antimony_builder as ab
 
-import matplotlib.pyplot as plt
 import numpy as np
 import re # type: ignore
 import unittest
@@ -114,11 +113,13 @@ def testMakeAdditionStatement(self):
         result = re.search("S2.* =.*S3", self.getStatement())
         self.assertTrue(result)
 
-    def testMakeSinusoidSignal(self):
+    def testMakeNoiseElement(self):
         if IGNORE_TEST:
             return
-        ot_name = self.builder.makeSinusoidSignal(1, 2, suffix="_S1_S2")
-        result = re.search("%s.*=.*1.*sin.*2" % ot_name, self.getStatement())
+        self.init()
+        noise_spec = cn.NoiseSpec(sine_amp=1, sine_freq=2, random_mag=3, random_std=4, offset=5, slope=0.0006)
+        ot_name = self.builder.makeNoiseElement(noise_spec, suffix="_S1_S2")
+        result = re.search("%s.*=.*1.*sin.*2.*3*lognormal.*4.*5.*6.*time" % ot_name, self.getStatement())
         self.assertTrue(result)
         self.builder.makeBoundarySpecies("S1")
         self.builder.makeAdditionStatement("S1", ot_name)
@@ -129,7 +130,8 @@ def testMakeFilterElement(self):
             return
         self.init()
         filter_in, filter_ot, calculation = self.builder.makeFilterElement(1.0, suffix="_S1_S3")
-        sin_ot = self.builder.makeSinusoidSignal(1, 2, suffix="_S1_S2")
+        noise_spec = cn.NoiseSpec(sine_amp=1, sine_freq=2)
+        sin_ot = self.builder.makeNoiseElement(noise_spec, suffix="_S1_S2")
         self.builder.makeAdditionStatement(filter_in, sin_ot)
         self.check()
 
@@ -196,26 +198,27 @@ def testMakeSISOClosedLoopSystem(self):
             return
         self.init()
         self.builder.makeBoundarySpecies("S1")
-        self.builder.makeSISOClosedLoopSystem("S1", "S3", kP=10, kI=1, setpoint=5,
-                           noise_amplitude=1, noise_frequency=2, disturbance_ampliude=0, disturbance_frequency=3)
+        noise_spec = cn.NoiseSpec(sine_amp=1, sine_freq=2)
+        self.builder.makeSISOClosedLoopSystem("S1", "S3", kP=10, kI=1, setpoint=5, noise_spec=noise_spec,
+                                              disturbance_spec=cn.NoiseSpec(sine_amp=2, sine_freq=3))
         self.check()
         #
         self.builder.initializeOutput()
         self.builder.makeBoundarySpecies("S1")
+        noise_spec = cn.NoiseSpec(sine_amp=1, sine_freq=2)
         self.builder.makeSISOClosedLoopSystem("S1", "S3", kP=1, kI=0.1, kF=0.1, kD=2, setpoint=5,
-                           noise_amplitude=1, noise_frequency=2, disturbance_ampliude=2, disturbance_frequency=3)
+                                              noise_spec=noise_spec)
         self.check()
         #
         self.builder.initializeOutput()
         self.builder.makeBoundarySpecies("S1")
         self.builder.makeSISOClosedLoopSystem("S1", "S3", kP=10, setpoint=5, kI=5, kD=2,
-                           noise_amplitude=1, noise_frequency=20, disturbance_ampliude=2, disturbance_frequency=20)
+                                              noise_spec=noise_spec)
         self.check()
         #
         self.builder.initializeOutput()
         self.builder.makeBoundarySpecies("S1")
-        self.builder.makeSISOClosedLoopSystem("S1", "S3", kP=10, setpoint=5,
-                           noise_amplitude=1, noise_frequency=20, disturbance_ampliude=2, disturbance_frequency=20,
+        self.builder.makeSISOClosedLoopSystem("S1", "S3", kP=10, setpoint=5, noise_spec=noise_spec,
                            initial_output_value=33)
         self.check()