Skip to content

Commit

Permalink
Plum through noise, disturbance.
Browse files Browse the repository at this point in the history
  • Loading branch information
@joseph-hellerstein
joseph-hellerstein committed May 14, 2024
1 parent 1389821 commit 67008dd
Show file tree
Hide file tree
Showing 10 changed files with 251 additions and 51 deletions.
46 changes: 35 additions & 11 deletions examples/Analysis-of-Immune-Response-Pneumocci.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -166,7 +166,7 @@
},
{
"cell_type": "code",
"execution_count": 9,
"execution_count": 6,
"metadata": {},
"outputs": [
{
Expand All @@ -188,7 +188,7 @@
},
{
"cell_type": "code",
"execution_count": 10,
"execution_count": 7,
"metadata": {},
"outputs": [
{
Expand All @@ -200,7 +200,7 @@
"TransferFunction(array([0.22887252]), array([0.31748021, 1.78740105, 1. ]))"
]
},
"execution_count": 10,
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
Expand All @@ -212,7 +212,7 @@
},
{
"cell_type": "code",
"execution_count": 11,
"execution_count": 8,
"metadata": {},
"outputs": [
{
Expand Down Expand Up @@ -254,14 +254,14 @@
},
{
"cell_type": "code",
"execution_count": 18,
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|█████████████████████████████████████████████| 1/1 [00:01<00:00, 1.47s/it]\n"
"100%|█████████████████████████████████████████████| 1/1 [00:01<00:00, 1.33s/it]\n"
]
},
{
Expand Down Expand Up @@ -304,7 +304,7 @@
},
{
"cell_type": "code",
"execution_count": 15,
"execution_count": 10,
"metadata": {},
"outputs": [
{
Expand Down Expand Up @@ -339,14 +339,14 @@
},
{
"cell_type": "code",
"execution_count": 16,
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|██████████| 100/100 [00:02<00:00, 34.76it/s]\n"
"100%|██████████| 100/100 [00:03<00:00, 25.23it/s]\n"
]
},
{
Expand Down Expand Up @@ -383,14 +383,14 @@
},
{
"cell_type": "code",
"execution_count": 17,
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|█████████████████████████████████████████████| 1/1 [00:01<00:00, 1.73s/it]\n"
"100%|█████████████████████████████████████████████| 1/1 [00:01<00:00, 1.29s/it]\n"
]
},
{
Expand Down Expand Up @@ -422,6 +422,30 @@
"# Effect of Disturbance and Noise"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"ename": "ValueError",
"evalue": "Invalid keyword arguments: {'noise_spec'}",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)",
"Cell \u001b[0;32mIn[13], line 2\u001b[0m\n\u001b[1;32m 1\u001b[0m noise_spec \u001b[38;5;241m=\u001b[39m ctl\u001b[38;5;241m.\u001b[39mNoiseSpec(sine_amp\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m1\u001b[39m)\n\u001b[0;32m----> 2\u001b[0m _ \u001b[38;5;241m=\u001b[39m \u001b[43mCTLSB\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mplotDesign\u001b[49m\u001b[43m(\u001b[49m\u001b[43mkP_spec\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mCTLSB\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mkP\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mkI_spec\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mCTLSB\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mkI\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtimes\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mTIMES\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mnoise_spec\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mnoise_spec\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43msetpoint\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mSETPOINT\u001b[49m\u001b[43m)\u001b[49m\n",
"File \u001b[0;32m~/home/Technical/repos/controlSBML/src/controlSBML/control_sbml.py:795\u001b[0m, in \u001b[0;36mControlSBML.plotDesign\u001b[0;34m(self, kP_spec, kI_spec, kD_spec, kF_spec, setpoint, sign, times, num_process, num_restart, selections, min_parameter_value, max_parameter_value, num_coordinate, is_report, **kwargs)\u001b[0m\n\u001b[1;32m 793\u001b[0m kD_spec\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39msetSpec(kD_spec)\n\u001b[1;32m 794\u001b[0m kF_spec\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39msetSpec(kF_spec)\n\u001b[0;32m--> 795\u001b[0m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_checkKwargs\u001b[49m\u001b[43m(\u001b[49m\u001b[43mPLOT_KEYS\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 796\u001b[0m option_dct \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mgetOptions(kP_spec\u001b[38;5;241m=\u001b[39mkP_spec,\n\u001b[1;32m 797\u001b[0m kI_spec\u001b[38;5;241m=\u001b[39mkI_spec,\n\u001b[1;32m 798\u001b[0m kD_spec\u001b[38;5;241m=\u001b[39mkD_spec,\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 805\u001b[0m times\u001b[38;5;241m=\u001b[39mtimes,\n\u001b[1;32m 806\u001b[0m \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs)\n\u001b[1;32m 807\u001b[0m kP_spec \u001b[38;5;241m=\u001b[39m option_dct[cn\u001b[38;5;241m.\u001b[39mO_KP_SPEC]\n",
"File \u001b[0;32m~/home/Technical/repos/controlSBML/src/controlSBML/control_sbml.py:260\u001b[0m, in \u001b[0;36mControlSBML._checkKwargs\u001b[0;34m(self, valids, invalids, **kwargs)\u001b[0m\n\u001b[1;32m 258\u001b[0m not_valids \u001b[38;5;241m=\u001b[39m keys\u001b[38;5;241m.\u001b[39mdifference(\u001b[38;5;28mset\u001b[39m(valids)) \u001b[38;5;66;03m# type: ignore\u001b[39;00m\n\u001b[1;32m 259\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mlen\u001b[39m(not_valids) \u001b[38;5;241m>\u001b[39m \u001b[38;5;241m0\u001b[39m:\n\u001b[0;32m--> 260\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mInvalid keyword arguments: \u001b[39m\u001b[38;5;132;01m%s\u001b[39;00m\u001b[38;5;124m\"\u001b[39m \u001b[38;5;241m%\u001b[39m \u001b[38;5;28mstr\u001b[39m(not_valids))\n\u001b[1;32m 261\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m invalids \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m 262\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m\n",
"\u001b[0;31mValueError\u001b[0m: Invalid keyword arguments: {'noise_spec'}"
]
}
],
"source": [
"noise_spec = ctl.NoiseSpec(sine_amp=1)\n",
"_ = CTLSB.plotDesign(kP_spec=CTLSB.kP, kI_spec=CTLSB.kI, times=TIMES, noise_spec=noise_spec, setpoint=SETPOINT)"
]
},
{
"cell_type": "markdown",
"metadata": {},
Expand Down
45 changes: 37 additions & 8 deletions src/controlSBML/antimony_builder.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,6 +73,7 @@ def __init__(self, antimony, symbol_dct=None):
if symbol_dct is None:
symbol_dct = util.makeRoadrunnerSymbolDct(self.roadrunner)
self.symbol_dct = dict(symbol_dct)
self.closed_loop_symbols = []

def copy(self):
"""
Expand Down Expand Up @@ -310,9 +311,10 @@ def makeNoiseElement(self, noise:cn.NoiseSpec, prefix:str="sinusoid", suffix:str
"""
self.addStatement("")
self.makeComment("Make sinusoid: %s" % str(noise))
name = prefix + suffix + OT
new_suffix = suffix + OT
name = self._makeScopedName(prefix, new_suffix)
statement = "%s := 0" % (name)
if not np.isclose(noise.sine_amp, 0):
if (not np.isclose(noise.sine_amp, 0)) and (not np.isclose(noise.sine_freq, 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)
Expand All @@ -323,14 +325,18 @@ def makeNoiseElement(self, noise:cn.NoiseSpec, prefix:str="sinusoid", suffix:str
self.addStatement(statement)
return name

@staticmethod
def _makeScopedName(prefix, suffix):
return prefix + suffix
def _makeScopedName(self, prefix, suffix, is_symbol=True):
name = prefix + suffix
if is_symbol:
self.closed_loop_symbols.append(name)
return name

def _makeInputOutputName(self, prefix, suffix):
base_name = self._makeScopedName(prefix, suffix)
base_name = self._makeScopedName(prefix, suffix, is_symbol=False)
name_in = base_name + IN
self.closed_loop_symbols.append(name_in)
name_ot = base_name + OT
self.closed_loop_symbols.append(name_ot)
return name_in, name_ot

def makeFilterElement(self, kF:float, prefix:str="filter", suffix:str=""):
Expand Down Expand Up @@ -388,7 +394,8 @@ def makeControlErrorSignal(self, setpoint, forward_output_name, sign, prefix="co
"""
self.addStatement("")
self.makeComment("Make the control error")
name_ot = prefix + suffix + OT
new_suffix = suffix + OT
name_ot = self._makeScopedName(prefix, new_suffix)
if sign == 1:
statement = "%s := %s - %s" % (name_ot, forward_output_name, str(setpoint))
elif sign == -1:
Expand Down Expand Up @@ -435,13 +442,17 @@ def makePIDControllerElement(self,
kD_name = base_name % "kD"
if kP is not None:
self.makeAdditionStatement(kP_name, str(kP), is_assignment=False)
self.closed_loop_symbols.append(kP_name)
if kI is not None:
self.makeAdditionStatement(kI_name, str(kI), is_assignment=False)
self.closed_loop_symbols.append(kI_name)
if kD is not None:
self.makeAdditionStatement(kD_name, str(kD), is_assignment=False)
self.closed_loop_symbols.append(kD_name)
# Make the derivative of the control error
if kD is not None:
derivative_error_name = base_name % "derivative_error"
self.closed_loop_symbols.append(derivative_error_name)
if filter_calculation is None:
filter_calculation = "rateOf(%s)" % (output_name) # type: ignore
sign_filter_calculation = "%d*(%s)" % (sign, filter_calculation) # type: ignore
Expand All @@ -450,6 +461,7 @@ def makePIDControllerElement(self,
# Make the integral of the control error
if kI is not None:
integral_error_name = base_name % "integral_error"
self.closed_loop_symbols.append(integral_error_name)
statement = "%s' = %s" % (integral_error_name, name_in)
self.addStatement(statement)
self.makeAdditionStatement(integral_error_name, 0, is_assignment=False) # integral_error_name = 0
Expand All @@ -469,6 +481,8 @@ def makeSISOClosedLoopSystem(self, input_name, output_name, kP=None, kI=None, kD
noise_spec=cn.NoiseSpec(), disturbance_spec=cn.DisturbanceSpec(),
initial_output_value=None, sign=-1):
"""
Creates a closed loop system with a single input and a single output.
Args:
input_name: str (input to system)
output_name: str (output from system)
Expand All @@ -485,7 +499,7 @@ def makeSISOClosedLoopSystem(self, input_name, output_name, kP=None, kI=None, kD
self.makeComment("**CREATING CLOSED LOOP SYSTEM**")
suffix = self.makeClosedLoopSuffix(input_name, output_name)
# Symbol for setpoint
setpoint_name = "setpoint" + suffix
setpoint_name = self._makeScopedName("setpoint", suffix)
self.makeAdditionStatement(setpoint_name, setpoint, is_assignment=False)
# Handle the initial value of the input
if initial_output_value is None:
Expand Down Expand Up @@ -522,6 +536,21 @@ def getInputManipulationName(self, input_name):
name = input_name
return name

def _getClosedLoopSymbols(self, input_name, output_name)->list[str]:
"""
Finds the names of closed loop symbols used for the system defined with the current input and output.
This method is used for validation purposes only.
Returns:
input_name: str
output_name: str
list[str]: list of symbols
"""
search_string = self.makeClosedLoopSuffix(input_name, output_name)
rr = te.loada(str(self))
symbols = [n for n in rr.keys() if (search_string in n) and (not "(" in n) and (not "[" in n)]
return symbols

def makeStaircase(self, input_name, times=cn.TIMES, initial_value=cn.DEFAULT_INITIAL_VALUE,
num_step=cn.DEFAULT_NUM_STEP, final_value=cn.DEFAULT_FINAL_VALUE):
"""
Expand Down
14 changes: 11 additions & 3 deletions src/controlSBML/constants.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,14 +28,20 @@ def __init__(self, sine_amp=0, sine_freq=0, random_mag=0, random_std=0, offset=N
self.random_mag = random_mag
self.random_std = random_std
if offset is None:
offset = sine_amp
if (np.isclose(sine_amp, 0) and np.isclose(slope, 0)):
offset = sine_amp
else:
offset = 0
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

def __eq__(self, other):
return self.__dict__ == other.__dict__

class DisturbanceSpec(NoiseSpec):
# Specifications for a disturbance model
Expand Down Expand Up @@ -126,7 +132,7 @@ def equals(self, other):
# Keyword options
O_AX = "ax"
O_AX2 = "ax2"
O_SELECTIONS = "selections" # Selections for the plot
O_DISTURBANCE_SPEC = "disturbance_spec"
O_END_TIME = "end_time"
O_FIGURE = "figure"
O_FIGSIZE = "figsize"
Expand All @@ -136,7 +142,6 @@ def equals(self, other):
O_KI_SPEC = "kI_spec"
O_KD_SPEC = "kD_spec"
O_KF_SPEC = "kF_spec"
O_STEP_VAL = "step_val"
O_INITIAL_VALUE = "initial_value"
O_INPUT_NAME = "input_name"
O_INPUT_NAMES = "input_names"
Expand All @@ -149,6 +154,7 @@ def equals(self, other):
O_LEGEND_SPEC = "legend_spec"
O_MARKERS = "markers"
O_MARKERSIZE = "markersize"
O_NOISE_SPEC = "noise_spec"
O_NUM_POINT = "num_point"
O_NUM_PROCESS = "num_process"
O_NUM_RESTART = "num_restart"
Expand All @@ -157,10 +163,12 @@ def equals(self, other):
O_OUTPUT_NAMES = "output_names"
O_POINTS_PER_TIME = "points_per_time"
O_PREDICTED = "predicted"
O_SELECTIONS = "selections" # Selections for the plot
O_SETPOINT = "setpoint"
O_SIGN = "sign"
O_STAIRCASE = "staircase"
O_START_TIME = "start_time"
O_STEP_VAL = "step_val"
O_SUPTITLE = "suptitle"
O_TIMES = "times"
O_TITLE = "title"
Expand Down
Loading

0 comments on commit 67008dd

Please sign in to comment.