Goal:
MUQ's modelling module provides tools for constructing, combining, and manipulating model components. In this setting, a model component can be any deterministic or stochastic inputoutput relationship. This example creates a single model component to evaluate the right hand side of an ODE describing a classic predatorprey system.
Mathematically, the system of interest is given by \begin{eqnarray} \frac{d P}{dt} &=& P \left[r\left(1  \frac{P}{K}\right)  s \frac{Q}{a + P}\right] \ \frac{d Q}{dt} &=& Q\left[ \frac{vP}{a + P}  u\right] \end{eqnarray} where $P$ is the prey population; $Q$ is the predator population; and $r$, $K$, $s$, $a$, $v$, and $u$ are model parameters. Thie example will define an object that takes $[P,Q]$ as inputs, and uses fixed values of $[r,K,s,a,v,u]$ to compute $dP/dt$ and $dQ/dt$.
ModPieces
The ModPiece class provides the basic building block for constructing sophisticated models in MUQ. To define a model component, MUQ users need to create a child of the ModPiece class with an EvaluateImpl function. Here we will create a child class that evaluates the right hand side of our ODE system in the EvaluateImpl function. For simplicity, we inherit from another class called OneInputNoDerivModPiece, which is simple a ModPiece with one input that does not provide any derivative information.
Import MUQ's modelling library
First we import the modelling library, which is called libmuqModelling in Python
import libmuqModelling
Define the model component
Now we are ready to define our model component. Our new class, PredPreyModel, will evaluate the right hand side of the ODE system above. The constructor takes the 6 model parameters as inputs and initialized the OneInputNoDerivModPiece with the input and output dimensions (both 2). Note that the input and output dimensions of all models in MUQ are constant, which means that these dimensions must be set in the ModPiece constructor.
The EvaluateImpl function must be defined in all children of the ModPiece class (including children of OneInputNoDerivModPiece). This function is where all the magic happens; it is where the model computation is performed.
class PredPreyModel(libmuqModelling.OneInputNoDerivModPiece): def __init__(self, r, k, s, a, u, v): libmuqModelling.OneInputNoDerivModPiece.__init__(self, 2, 2) self.preyGrowth = r self.preyCapacity = k self.predationRate = s self.predationHunger = a self.predatorLoss = u self.predatorGrowth = v def EvaluateImpl(self, ins): preyPop = ins[0] predPop = ins[1] output = [None]*2 output[0] = preyPop * (self.preyGrowth * (1  preyPop / self.preyCapacity)  self.predationRate * predPop / (self.predationHunger + preyPop)) output[1] = predPop * (self.predatorGrowth * preyPop / (self.predationHunger + preyPop)  self.predatorLoss) return output
Create the model instance
Now that we've defined the model in our PredPreyModel class, lets put it to work. In the following cell, we create an instance of the PredPreyModel class.
# Before constructing the model, we will set the fixed model values. preyGrowth = 0.8 preyCapacity = 100 predationRate = 4.0 predationHunger = 15 predatorLoss = 0.6 predatorGrowth = 0.8 # create an instance of the ModPiece defining the predator prey model. predatorPreyModel = PredPreyModel(preyGrowth, preyCapacity, predationRate, predationHunger, predatorLoss, predatorGrowth)
Evaluate the model
To evaluate our model, we can simply call the Evaluate function of the PredPreyModel. Note that this function is inherited from the ModPiece and is distinct from the EvaluateImpl function defined above. Calling Evaluate instead of EvaluateImpl allows MUQ to perform some additional error correcting and caching that is not present in EvaluateImpl.
populations = [50.0, 5.0] growthRates = predatorPreyModel.Evaluate(populations)
Display the model output
Print the results.
print "The prey growth rate evaluated at (", populations[0], ",", populations[1], ") is ", growthRates[0] print "The predator growth rate evaluated at (", populations[0], ",", populations[1], ") is ", growthRates[1]
The prey growth rate evaluated at ( 50.0 , 5.0 ) is 4.61538461538 The predator growth rate evaluated at ( 50.0 , 5.0 ) is 0.0769230769231
Get information about the model
We can also ask the model component for it's input and output sizes. Every Model in MUQ has a constant vector, "inputSizes", that contains the length of each vector valued input. Since our Model only has one input, inputSizes only has one entry. The output dimension of the model is stored in the "outputSize" variable. Since MUQ models can only have one output, "outputSize" is simply an integer.
print "The input size is ", predatorPreyModel.inputSizes[0] print "The output size is ", predatorPreyModel.outputSize
The input size is 2 The output size is 2
Completed code:
import libmuqModelling class PredPreyModel(libmuqModelling.OneInputNoDerivModPiece): def __init__(self, r, k, s, a, u, v): libmuqModelling.OneInputNoDerivModPiece.__init__(self, 2, 2) self.preyGrowth = r self.preyCapacity = k self.predationRate = s self.predationHunger = a self.predatorLoss = u self.predatorGrowth = v def EvaluateImpl(self, ins): preyPop = ins[0] predPop = ins[1] output = [None]*2 output[0] = preyPop * (self.preyGrowth * (1  preyPop / self.preyCapacity)  self.predationRate * predPop / (self.predationHunger + preyPop)) output[1] = predPop * (self.predatorGrowth * preyPop / (self.predationHunger + preyPop)  self.predatorLoss) return output # Before constructing the model, we will set the fixed model values. preyGrowth = 0.8 preyCapacity = 100 predationRate = 4.0 predationHunger = 15 predatorLoss = 0.6 predatorGrowth = 0.8 # create an instance of the ModPiece defining the predator prey model. predatorPreyModel = PredPreyModel(preyGrowth, preyCapacity, predationRate, predationHunger, predatorLoss, predatorGrowth) populations = [50.0, 5.0] growthRates = predatorPreyModel.Evaluate(populations) print "The prey growth rate evaluated at (", populations[0], ",", populations[1], ") is ", growthRates[0] print "The predator growth rate evaluated at (", populations[0], ",", populations[1], ") is ", growthRates[1] print "The input size is ", predatorPreyModel.inputSizes[0] print "The output size is ", predatorPreyModel.outputSize

Join our slack workspace to connect with other users, get help, and discuss new features.
This work was supported by the DOE Office of Science through the QUEST and FastMath SciDAC Institutes.
Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.