population examples update

.gitignore

@@ -165,6 +165,10 @@ cython_debug/
 .vscode/settings.json
 psymple/api.py
 examples/test.ipynb
+examples/test*
 psymple/functions.py
 doctest_runner.py
 test_json.json
+examples/mixing_problems/3-state_dependent_rates.py
+examples/bug.py
+examples/population_dynamics/3-reponse_functions.py

docs/examples/population_dynamics/logistic_growth.md

@@ -1,5 +1,9 @@
 # The logistic equation
 
+!!! info "Raw code"
+
+    The raw code for this example without explanations can be found [here](https://github.com/casasglobal-org/psymple/blob/main/examples/population_dynamics/2-logistic_population.py).
+
 This example follows on from [Malthusian growth](malthusian_population.md). The code for a simple Malthusian population produced there is available in the drop-down box below.
 
 ??? example "Malthusian population"
@@ -98,15 +102,15 @@ system_1 = System(logistic_simple)
 system_1.compile()
 
 sim_1 = system_1.create_simulation(initial_values={"x": 1})
-sim_1.simulate(t_end=25)
+sim_1.simulate(t_end=100)
 sim_1.plot_solution()
 
 
 system_2 = System(logistic_pop)
 system_2.compile()
 
 sim_2 = system_2.create_simulation(initial_values={"x": 1})
-sim_2.simulate(t_end=25)
+sim_2.simulate(t_end=100)
 sim_2.plot_solution()
 ```
 

docs/examples/population_dynamics/malthusian_population.md

@@ -1,5 +1,9 @@
 # Malthusian growth
 
+!!! info "Raw code"
+
+    The raw code for this example without explanations can be found [here](https://github.com/casasglobal-org/psymple/blob/main/examples/population_dynamics/1-malthusian_population.py).
+
 A fundamental equation in population dynamics describes a population which reproduces as a constant rate \(r\). Without any other influence, this population grows exponentially according to the differential equation \( \frac{dx}{dt} = rx \). This phenomenon is known as *Malthusian growth*.
 
 This example covers:

docs/examples/population_dynamics/predator_prey.md

@@ -1,5 +1,9 @@
 # Predator-prey systems
 
+!!! info "Raw code"
+
+    The raw code for this example without explanations can be found [here](https://github.com/casasglobal-org/psymple/blob/main/examples/population_dynamics/3-predator_prey.py).
+
 A predator-prey system was introduced in the [defining ODEs](../../components/variable_ported_objects.md#example) tutorial. This example shows how to build this system on top of the single-species population models constructed in [Malthusian growth](malthusian_population.md) and [logistic growth](logistic_growth.md).
 
 ## Lotka-Volterra model
@@ -70,8 +74,8 @@ ecosystem = CompositePortedObject(
     children=[pred, prey, interaction],
     variable_ports=["x", "y"],
     variable_wires=[
-        (["prey.x", "interaction.x"], "x"),
-        (["pred.x", "interaction.y"], "y"),
+        (["prey.x", "pred-prey.x"], "x"),
+        (["pred.x", "pred-prey.y"], "y"),
     ]
 )
 ```

examples/population_dynamics/1-malthusian_population.py

@@ -11,25 +11,7 @@
 from psymple.build import System
 
 """
-The most basic population dynamics model is \( dx/dt = rx \).
-
-This equation can be formed as a VariablePortedObject with an assignment 
-`assg = {"variable": "x", "expression": "r*x"}`. This entry is equivalent to the shorthand
-`assg = ("x", "r*x")`.
-
-The value of parameter 'r' could be specified immediately in the assignment, for example we
-could instead write the assignment ("x", "0.1*x"). 
-
-Alternatively, psymple allows you to specify a value for 'r' through an input port. This
-has two advantages:
-    (1) the variable ported object becomes a general, reusable object
-    (2) the default value can be overwritten from elsewhere in a larger system
-
-We will see an example of (2) shortly. For now, to a default value of 0.1 for the parameter
-'r', we add an input port `port = {"name": "r", "default_value": 0.1}`. As for assignments,
-this is equivalent to the shorthand `port = ("r", 0.1)`.
-
-Putting these together, we define the VariablePortedObject `pop` as follows.
+Exponential population growth
 """
 
 pop = VariablePortedObject(
@@ -39,26 +21,32 @@
 )
 
 """
-The following commands create, run and plot a simulation of our equation with initial condition
-'x = 1'. 
+Running the simulation 
 """
 
 system_1 = System(pop)
 system_1.compile()
 
-sim_1 = system_1.create_simulation("sim", solver="discrete", initial_values={"x": 1})
-sim_1.simulate(10, n_steps=10)
+sim_1 = system_1.create_simulation(initial_values={"x": 1})
+sim_1.simulate(t_end=25)
 sim_1.plot_solution()
 
 """
-Now let's see an example of overriding a default parameter value. Suppose that we individually know
-the birth rate 'b' and death rate 'd' for the population, with values 0.4 and 0.2, respectively. 
-The combined Malthusian rate is 'r = b - d'. We can perform this calculation in a FunctionalPortedObject
-with assignment `assg = {"parameter": "r", "expression": "b-d"}`. As before, this is equivalent to the
-shorthand `assg = ("r", "b-d")`.
+Changing parameters at simulation
+"""
 
-Again, we need to tell psymple what we mean by 'b' and 'd', which we can do as default input ports on
-this object.
+system_1 = System(pop)
+system_1.compile()
+
+sim_2 = system_1.create_simulation(
+    initial_values={"x": 1}, 
+    input_parameters={"r": 0.2}
+)
+sim_2.simulate(t_end=25)
+sim_2.plot_solution()
+
+"""
+Overriding a default values with a function
 """
 
 rate = FunctionalPortedObject(
@@ -68,11 +56,7 @@
 )
 
 """
-Next, we need to tell our variable object 'pop' to read its value of 'r' from the value of 'r' produced
-by the function 'rate'. We do this in a CompositePortedObject containing both 'pop' and 'rate' as 
-children, using a directed wire from 'rate.r' to 'malthusian_pop.r'. This is written as 
-`wire = {"source": "rate.r", "destination": "malthusian_pop.r"}`. There is, of course, an equivalent 
-shorthand `wire = ("rate.r", "malthusian_pop.r")`. 
+Composing the function with the differential equation
 """
 
 pop_system = CompositePortedObject(
@@ -82,14 +66,44 @@
 )
 
 """
-Notice that we did not have to redefine the variable object 'pop'. Since we defined it generally,
-we can reuse its mechanics over and over again. The following commands simulate and plot this
-new system.
+Exposing inputs
+"""
+
+pop_system = CompositePortedObject(
+    name="malthusian_pop_system",
+    children=[pop, rate],
+    input_ports=[("b", 0.4), ("d", 0.2)],
+    directed_wires=[
+        ("rate.r", "malthusian_pop.r"),
+        ("b", "rate.b"),
+        ("d", "rate.d")
+    ],
+)
+
+"""
+Compilation and simulation
 """
 
 system_2 = System(pop_system)
 system_2.compile()
 
-sim_2 = system_2.create_simulation("sim", solver="discrete", initial_values={"malthusian_pop.x": 1})
-sim_2.simulate(10, n_steps=10)
-sim_2.plot_solution()
+sim_3 = system_2.create_simulation(initial_values={"malthusian_pop.x": 1})
+sim_3.simulate(t_end=25)
+sim_3.plot_solution()
+
+"""
+Creating an interface
+"""
+
+pop_system = CompositePortedObject(
+    name="malthusian_pop_system",
+    children=[pop, rate],
+    input_ports=[("b", 0.4), ("d", 0.2)],
+    directed_wires=[
+        ("rate.r", "malthusian_pop.r"),
+        ("b", "rate.b"),
+        ("d", "rate.d")
+    ],
+    variable_ports=["x"],
+    variable_wires=[(["malthusian_pop.x"], "x")],
+)

examples/population_dynamics/2-logistic_population.py

@@ -5,13 +5,12 @@
 
 from psymple.build import (
     CompositePortedObject,
-    FunctionalPortedObject,
     VariablePortedObject,
     System,
 )
 
 """
-Recall the following variable object 'pop' capturing Malthusian growth in a reusable way.
+Basic population model
 """
 
 pop = VariablePortedObject(
@@ -21,131 +20,65 @@
 )
 
 """
-Logistic population growth consists of a Malthusian growth component and a density-dependent limit to some
-carrying capacity K. The differential equation is dx/dt = rx(1 - x/K). 
-
-We could create a new variable object capturing this equation, but since we have already have the object
-'pop' capturing Malthusian growth, it feels as though we are repeating some work. 
-
-Here is an alternative way of viewing the logistic growth equation: its right-hand side is the sum of 
-two components: a Malthusian component rx, and a limit component -x^2/K. Psymple allows us to combine
-two variable objects, each capturing one of these components, to a single dynamic model in a process
-called composition. 
-
-In simple terms, the composition of the equations dx/dt = f(x,t) and dx/dt = g(x,t) is the equation
-dx/dt = f(x,t) + g(x,t). Let's see how to do this in psymple. First we need a variable object 
-capturing the logistic component, with a default value of 'K = 10'.
+Direct logistic implementation
 """
 
-limit = VariablePortedObject(
-    name="logistic_limit",
-    assignments=[("x", "- r* x**2 / K")],
+logistic_simple = VariablePortedObject(
+    name="logistic_simple",
+    assignments=[("x", "r*x*(1-x/K)")],
     input_ports=[("r", 0.1), ("K", 10)],
 )
 
 """
-Next, in a composite object, we use a variable wire to tell psymple to compose the variables
-'malthusian_pop.x' and 'logistic_limit.x' together. To do this, we specify a wire
-`wire = {"child_ports": ["malthusian_pop.x", "logistic_limit.x"], "output_name" = "x"}`. There is
-a shorthand for this, which at first appears unnatural:
-`wire = (["malthusian_pop.x", "logistic_limit.x"], None, "x")`. We will discuss this shortly.
-"""
-
-logistic_pop = CompositePortedObject(
-    name="logistic_pop",
-    children=[pop, limit],
-    variable_wires=[(["malthusian_pop.x", "logistic_limit.x"], None, "x")],
-)
-
-system = System(logistic_pop)
-system.compile()
-
-sim = system.create_simulation("sim", solver="discrete", initial_values={"x": 1})
-sim.simulate(100, n_steps=10)
-sim.plot_solution()
-
+Structural implementation
 """
-You'll that the simulation looks exactly as expected: the population initially undergoes exponential
-growth, before its growth rate slows and it approaches the limit of 10. 
-
-Let's discuss some efficiencies and good practice. Suppose instead that we did implement the logistic
-growth equation as a single variable object:
 
-logistic_pop = VariablePortedObject(
-    name="logistic_pop",
-    assignments=[("x", "r*(x - x**2/K)")],
+limit = VariablePortedObject(
+    name="limit",
+    assignments=[("x", "- r/K* x**2")],
     input_ports=[("r", 0.1), ("K", 10)],
 )
 
-As discussed, this implementation does not exploit the reusability features of psymple, but there are three
-apparent advantages:
-    (1) The default value for "r" only had to be specified once.
-    (2) The variable "x" is available for composition in larger systems.
-    (3) The object as a whole has a clean form: it is an object depending on input parameters "r" and "K",
-        and makes the variable "x" available to the wider system.
-
-Fortunately, composite objects in psymple allow for all of these to be implemented, and this is how
-psymple is intended to be used in best practice to allow for full modularity and reusablity. To achieve this,
-composite objects can be given their own ports: input ports to read parameters or specify default values, and
-variable ports to expose variables and their differential equations for composition. These ports connect to 
-wires in the ways we've already seen.
-
-With this in mind, let's see how we would build up the logistic_pop object to both maximise modularity and
-reusability. We first take the two variable components 'pop' and 'limit'.
-"""
-
-pop = VariablePortedObject(
-    name="malthusian_pop",
-    assignments=[("x", "r*x")],
-)
-
-limit = VariablePortedObject(
-    name="logistic_limit",
-    assignments=[("x", "- r* x**2 / K")],
+logistic_pop = CompositePortedObject(
+    name="logistic_pop",
+    children=[pop, limit],
+    variable_ports = ["x"],
+    variable_wires=[(["malthusian_pop.x", "limit.x"], "x")],
 )
 
 """
-Notice that we have not specified default input ports. Psymple is able to automatically generate the port 
-'r' in 'pop' and the ports 'r' and 'K' in 'limit'. Given we know that we will specify these later, we don't 
-need to give default values at this stage.
-
-Next, we will define the composite object 'logistic_pop' 
+Sharing input parameters
 """
 
 logistic_pop = CompositePortedObject(
     name="logistic_pop",
     children=[pop, limit],
-    input_ports = [("r", 0.1), ("K", 10)],
-    variable_ports=["x"],
-    variable_wires=[(["malthusian_pop.x", "logistic_limit.x"], "x")],
+    input_ports=[("r", 0.1), ("K", 10)],
     directed_wires=[
-        ("r", ["malthusian_pop.r", "logistic_limit.r"]),
-        ("K", "logistic_limit.K"),
+        ("r", ["malthusian_pop.r", "limit.r"]), # (1)!
+        ("K", "limit.K"),
     ],
+    variable_ports = ["x"],
+    variable_wires=[(["malthusian_pop.x", "limit.x"], "x")],
 )
 
 """
-There's a bit to look at here. First, notice that we have given input ports to 'logistic_pop' in
-the same way we did for our previous variable objects. We have also told psymple to form a
-variable port 'x' on 'logistic_pop', so that we can expose internal information.
+System and simulation - checking implementations are the same
+"""
 
-Next, notice that the variable wire has nearly the same format - but this time the given format
-is shorthand for `wire = {"child_ports": ["malthusian_pop.x", "logistic_limit.x"], "parent_port": "x"}`.
-So instead of just assigning the name of the composition of these variables to 'x', we're now saying to
-expose that composition at variable port 'x'.
 
-Finally, we've used two directed wires to connect our ports and objects. The first wire,
-`("r", ["malthusian_pop.r", "logistic_limit.r"])` transfers the default value at port 'r' to *both*
-'malthusian_pop.r' and 'logistic_limit.r'. The second wire transfers the default value of 'K' to
-'logistic_limit.K'. 
+system_1 = System(logistic_simple)
+system_1.compile()
+
+sim_1 = system_1.create_simulation(initial_values={"x": 1})
+sim_1.simulate(t_end=100)
+sim_1.plot_solution()
 
-Let's finish by checking that this new, general object does the same thing as before.
-"""
 
-system = System(logistic_pop)
-system.compile()
+system_2 = System(logistic_pop)
+system_2.compile()
 
-sim = system.create_simulation("sim", solver="discrete", initial_values={"x": 1})
-sim.simulate(100, n_steps=10)
-sim.plot_solution()
+sim_2 = system_2.create_simulation(initial_values={"x": 1})
+sim_2.simulate(t_end=100)
+sim_2.plot_solution()