Header menu logo bristlecone

Composing model parts together

Bristlecone has been designed to enable swift model composition, whereby smaller models may be composed into larger systems of models.

To get started, open the Bristlecone libraries.

open Bristlecone // Opens Bristlecone core library and estimation engine
open Bristlecone.Language // Open the language for writing Bristlecone models
open FSharp.Data.UnitSystems.SI.UnitNames

Nested model systems

One key use of compositional modelling is to understand which (ecological) theories best represent one or many particular time-series data. A nested model is one where one or more processes may be substituted for different theoretical forms and competed to determine which combination of model components best represents the given data.

Using Bristlecone, we can simply define a nested model as a Bristlecone model but where the definition takes some additional model components as function parameters.

We can define a nestable model, in which a number of components may be swapped for alternative mathematical forms based on (ecological) theory. Here, we will apply a simple model of non-linear plant growth. To do this, we first define a base model as a ModelSystem that takes a number of interchangable components:

let r = parameter "r" noConstraints 0.01< / Time.month> 1.00< / Time.month>
let m = state<kilogram> "m"

let baseModel growthLimit lossRate =
    Model.empty<Time.month>
    |> Model.addRateEquation m (P r * growthLimit This<kilogram> - lossRate This<kilogram>)
    |> Model.estimateParameter r
    |> Model.useLikelihoodFunction (ModelLibrary.Likelihood.sumOfSquares [ Require.state m ])

You can see from the type signatures that growthLimit returns a kilogram value, whereas the loss rate returns a rate (kilogram per month).

Importantly, unlike when normally building models, we do not call Model.compile at the end. This step will be done later.

We can then design the growth limit and loss rate components that will fit into this model. These are made using the modelComponent and subComponent functions as follows:

open Bristlecone.Language.Components

let K = parameter "K" notNegative 10.0<kilogram> 50.0<kilogram>

let growthLimit =
    modelComponent
        "Limitation to growth"
        [

          subComponent "Linear" (fun _ -> Constant 1.<kilogram>)

          subComponent "Monomolecular" (fun (m: ModelExpression<kilogram>) -> P K - m)
          |> estimateParameter K

          subComponent "Gompertz" (fun (m: ModelExpression<kilogram>) -> m * Logarithm(P K / m))
          |> estimateParameter K

          subComponent "Logistic (three parameter)" (fun (m: ModelExpression<kilogram>) -> m * (Constant 1. - m / P K))
          |> estimateParameter K ]

If a component requires additional parameters over the base model, these may be added by piping into the estimateParameter function as above.

We can similarly define the other required component for the base model - the loss rate:

let alpha = parameter "alpha" notNegative 0.001< / Time.month> 0.10< / Time.month>

let lossRate =
    modelComponent
        "Biomass loss rate"
        [

          subComponent "None" (fun _ -> Constant 0.<kilogram / Time.month>)

          subComponent "Density-dependent" (fun m -> m * P alpha)
          |> estimateParameter alpha ]

Once the components are set up, we can compile the nested model by making all possible combinations of the model components. We do this using the Hypotheses module, which will compile a list of model hypotheses for us to test.

let hypotheses =
    baseModel
    |> Hypotheses.createFromModel
    |> Hypotheses.apply growthLimit
    |> Hypotheses.apply lossRate
    |> Hypotheses.compile

The resultant constructed hypotheses are shown in the following printout:

[Hypothesis
   ({ Parameters =
       Pool
         (map
            [(ShortCode "r",
              { Name = "r"
  ToTensorRealIO = <fun:make@119-1>
  FromTensorRealIO = <fun:make@124-2>
  GetConstraint = <fun:make@130-3>
  TryGetReal = <fun:make@131-4>
  TryGetBounds = <fun:make@133-6> })])
      EnvironmentKeys = []
      Equations =
       DifferentialEqs (map [(ShortCode "m", <fun:ToFSharpFunc@3936-2>)])
      Measures = map []
      Initialisers = map []
      NegLogLikelihood = { RequiredCodes = [State ShortCode "m"]
                           Evaluate = <fun:Pipe #1 input at line 27@27> } },
    [{ Component = "Limitation to growth"
       Implementation = "Linear" }; { Component = "Biomass loss rate"
                                      Implementation = "None" }]);
 Hypothesis
   ({ Parameters =
       Pool
         (map
            [(ShortCode "K",
              { Name = "K"
  ToTensorRealIO = <fun:make@119-1>
  FromTensorRealIO = <fun:make@124-2>
  GetConstraint = <fun:make@130-3>
  TryGetReal = <fun:make@131-4>
  TryGetBounds = <fun:make@133-6> });
             (ShortCode "r",
              { Name = "r"
  ToTensorRealIO = <fun:make@119-1>
  FromTensorRealIO = <fun:make@124-2>
  GetConstraint = <fun:make@130-3>
  TryGetReal = <fun:make@131-4>
  TryGetBounds = <fun:make@133-6> })])
      EnvironmentKeys = []
      Equations =
       DifferentialEqs (map [(ShortCode "m", <fun:ToFSharpFunc@3936-2>)])
      Measures = map []
      Initialisers = map []
      NegLogLikelihood = { RequiredCodes = [State ShortCode "m"]
                           Evaluate = <fun:Pipe #1 input at line 27@27> } },
    [{ Component = "Limitation to growth"
       Implementation = "Monomolecular" }; { Component = "Biomass loss rate"
                                             Implementation = "None" }]);
 Hypothesis
   ({ Parameters =
       Pool
         (map
            [(ShortCode "K",
              { Name = "K"
  ToTensorRealIO = <fun:make@119-1>
  FromTensorRealIO = <fun:make@124-2>
  GetConstraint = <fun:make@130-3>
  TryGetReal = <fun:make@131-4>
  TryGetBounds = <fun:make@133-6> });
             (ShortCode "r",
              { Name = "r"
  ToTensorRealIO = <fun:make@119-1>
  FromTensorRealIO = <fun:make@124-2>
  GetConstraint = <fun:make@130-3>
  TryGetReal = <fun:make@131-4>
  TryGetBounds = <fun:make@133-6> })])
      EnvironmentKeys = []
      Equations =
       DifferentialEqs (map [(ShortCode "m", <fun:ToFSharpFunc@3936-2>)])
      Measures = map []
      Initialisers = map []
      NegLogLikelihood = { RequiredCodes = [State ShortCode "m"]
                           Evaluate = <fun:Pipe #1 input at line 27@27> } },
    [{ Component = "Limitation to growth"
       Implementation = "Gompertz" }; { Component = "Biomass loss rate"
                                        Implementation = "None" }]);
 Hypothesis
   ({ Parameters =
       Pool
         (map
            [(ShortCode "K",
              { Name = "K"
  ToTensorRealIO = <fun:make@119-1>
  FromTensorRealIO = <fun:make@124-2>
  GetConstraint = <fun:make@130-3>
  TryGetReal = <fun:make@131-4>
  TryGetBounds = <fun:make@133-6> });
             (ShortCode "r",
              { Name = "r"
  ToTensorRealIO = <fun:make@119-1>
  FromTensorRealIO = <fun:make@124-2>
  GetConstraint = <fun:make@130-3>
  TryGetReal = <fun:make@131-4>
  TryGetBounds = <fun:make@133-6> })])
      EnvironmentKeys = []
      Equations =
       DifferentialEqs (map [(ShortCode "m", <fun:ToFSharpFunc@3936-2>)])
      Measures = map []
      Initialisers = map []
      NegLogLikelihood = { RequiredCodes = [State ShortCode "m"]
                           Evaluate = <fun:Pipe #1 input at line 27@27> } },
    [{ Component = "Limitation to growth"
       Implementation = "Logistic (three parameter)" };
     { Component = "Biomass loss rate"
       Implementation = "None" }]);
 Hypothesis
   ({ Parameters =
       Pool
         (map
            [(ShortCode "alpha",
              { Name = "alpha"
  ToTensorRealIO = <fun:make@119-1>
  FromTensorRealIO = <fun:make@124-2>
  GetConstraint = <fun:make@130-3>
  TryGetReal = <fun:make@131-4>
  TryGetBounds = <fun:make@133-6> });
             (ShortCode "r",
              { Name = "r"
  ToTensorRealIO = <fun:make@119-1>
  FromTensorRealIO = <fun:make@124-2>
  GetConstraint = <fun:make@130-3>
  TryGetReal = <fun:make@131-4>
  TryGetBounds = <fun:make@133-6> })])
      EnvironmentKeys = []
      Equations =
       DifferentialEqs (map [(ShortCode "m", <fun:ToFSharpFunc@3936-2>)])
      Measures = map []
      Initialisers = map []
      NegLogLikelihood = { RequiredCodes = [State ShortCode "m"]
                           Evaluate = <fun:Pipe #1 input at line 27@27> } },
    [{ Component = "Limitation to growth"
       Implementation = "Linear" }; { Component = "Biomass loss rate"
                                      Implementation = "Density-dependent" }]);
 Hypothesis
   ({ Parameters =
       Pool
         (map
            [(ShortCode "K",
              { Name = "K"
  ToTensorRealIO = <fun:make@119-1>
  FromTensorRealIO = <fun:make@124-2>
  GetConstraint = <fun:make@130-3>
  TryGetReal = <fun:make@131-4>
  TryGetBounds = <fun:make@133-6> });
             (ShortCode "alpha",
              { Name = "alpha"
  ToTensorRealIO = <fun:make@119-1>
  FromTensorRealIO = <fun:make@124-2>
  GetConstraint = <fun:make@130-3>
  TryGetReal = <fun:make@131-4>
  TryGetBounds = <fun:make@133-6> });
             (ShortCode "r",
              { Name = "r"
  ToTensorRealIO = <fun:make@119-1>
  FromTensorRealIO = <fun:make@124-2>
  GetConstraint = <fun:make@130-3>
  TryGetReal = <fun:make@131-4>
  TryGetBounds = <fun:make@133-6> })])
      EnvironmentKeys = []
      Equations =
       DifferentialEqs (map [(ShortCode "m", <fun:ToFSharpFunc@3936-2>)])
      Measures = map []
      Initialisers = map []
      NegLogLikelihood = { RequiredCodes = [State ShortCode "m"]
                           Evaluate = <fun:Pipe #1 input at line 27@27> } },
    [{ Component = "Limitation to growth"
       Implementation = "Monomolecular" };
     { Component = "Biomass loss rate"
       Implementation = "Density-dependent" }]);
 Hypothesis
   ({ Parameters =
       Pool
         (map
            [(ShortCode "K",
              { Name = "K"
  ToTensorRealIO = <fun:make@119-1>
  FromTensorRealIO = <fun:make@124-2>
  GetConstraint = <fun:make@130-3>
  TryGetReal = <fun:make@131-4>
  TryGetBounds = <fun:make@133-6> });
             (ShortCode "alpha",
              { Name = "alpha"
  ToTensorRealIO = <fun:make@119-1>
  FromTensorRealIO = <fun:make@124-2>
  GetConstraint = <fun:make@130-3>
  TryGetReal = <fun:make@131-4>
  TryGetBounds = <fun:make@133-6> });
             (ShortCode "r",
              { Name = "r"
  ToTensorRealIO = <fun:make@119-1>
  FromTensorRealIO = <fun:make@124-2>
  GetConstraint = <fun:make@130-3>
  TryGetReal = <fun:make@131-4>
  TryGetBounds = <fun:make@133-6> })])
      EnvironmentKeys = []
      Equations =
       DifferentialEqs (map [(ShortCode "m", <fun:ToFSharpFunc@3936-2>)])
      Measures = map []
      Initialisers = map []
      NegLogLikelihood = { RequiredCodes = [State ShortCode "m"]
                           Evaluate = <fun:Pipe #1 input at line 27@27> } },
    [{ Component = "Limitation to growth"
       Implementation = "Gompertz" }; { Component = "Biomass loss rate"
                                        Implementation = "Density-dependent" }]);
 Hypothesis
   ({ Parameters =
       Pool
         (map
            [(ShortCode "K",
              { Name = "K"
  ToTensorRealIO = <fun:make@119-1>
  FromTensorRealIO = <fun:make@124-2>
  GetConstraint = <fun:make@130-3>
  TryGetReal = <fun:make@131-4>
  TryGetBounds = <fun:make@133-6> });
             (ShortCode "alpha",
              { Name = "alpha"
  ToTensorRealIO = <fun:make@119-1>
  FromTensorRealIO = <fun:make@124-2>
  GetConstraint = <fun:make@130-3>
  TryGetReal = <fun:make@131-4>
  TryGetBounds = <fun:make@133-6> });
             (ShortCode "r",
              { Name = "r"
  ToTensorRealIO = <fun:make@119-1>
  FromTensorRealIO = <fun:make@124-2>
  GetConstraint = <fun:make@130-3>
  TryGetReal = <fun:make@131-4>
  TryGetBounds = <fun:make@133-6> })])
      EnvironmentKeys = []
      Equations =
       DifferentialEqs (map [(ShortCode "m", <fun:ToFSharpFunc@3936-2>)])
      Measures = map []
      Initialisers = map []
      NegLogLikelihood = { RequiredCodes = [State ShortCode "m"]
                           Evaluate = <fun:Pipe #1 input at line 27@27> } },
    [{ Component = "Limitation to growth"
       Implementation = "Logistic (three parameter)" };
     { Component = "Biomass loss rate"
       Implementation = "Density-dependent" }])]

For further analysis, you may choose to use the orchestration functionality to run many hypotheses and replicates (multiple fits per hypothesis) in parallel.

Multiple items
module Bristlecone from Bristlecone
<namespacedoc><summary>The core library of Bristlecone, containing model-fitting functions.</summary></namespacedoc>
Main functionality of Bristlecone, including functions to scaffold `ModelSystem`s and for model-fitting (tests and real fits).


--------------------
namespace Bristlecone
module Language from Bristlecone
<summary> An F# Domain Specific Language (DSL) for scripting with Bristlecone. </summary>
namespace Microsoft.FSharp
namespace Microsoft.FSharp.Data
namespace Microsoft.FSharp.Data.UnitSystems
namespace Microsoft.FSharp.Data.UnitSystems.SI
namespace Microsoft.FSharp.Data.UnitSystems.SI.UnitNames
val r: IncludedParameter</Time.month>
Multiple items
val parameter: code: string -> con: Parameter.Constraint -> lower: float<'u> -> upper: float<'u> -> IncludedParameter<'u>
<summary> Define an estimatable parameter for a Bristlecone model. </summary>

--------------------
[<Measure>] type parameter
val noConstraints: Parameter.Constraint
Multiple items
val Time<'t> : ModelExpression<'t>

--------------------
module Time from Bristlecone
[<Measure>] type month
val m: StateId<kilogram>
val state: name: string -> StateId<'u>
[<Measure>] type kilogram
val baseModel: growthLimit: (ModelExpression<kilogram> -> ModelExpression<kilogram>) -> lossRate: (ModelExpression<kilogram> -> ModelExpression<kilogram/Time.month>) -> ModelBuilder.ModelBuilder<Time.month>
val growthLimit: (ModelExpression<kilogram> -> ModelExpression<kilogram>)
val lossRate: (ModelExpression<kilogram> -> ModelExpression<kilogram/Time.month>)
module Model from Bristlecone.Language
<summary> Terms for scaffolding a model system for use with Bristlecone. </summary>
val empty<'time> : ModelBuilder.ModelBuilder<'time>
val addRateEquation: name: StateId<'state> -> expr: ModelExpression<'state/'time> -> mb: ModelBuilder.ModelBuilder<'time> -> ModelBuilder.ModelBuilder<'time>
val P: name: IncludedParameter<'u> -> ModelExpression<'u>
val This<'u> : ModelExpression<'u>
val estimateParameter: p: IncludedParameter<'u> -> builder: ModelBuilder.ModelBuilder<'time> -> ModelBuilder.ModelBuilder<'time>
val useLikelihoodFunction: likelihoodFn: ModelSystem.Likelihood<ModelSystem.state> -> builder: ModelBuilder.ModelBuilder<'u> -> ModelBuilder.ModelBuilder<'u>
namespace Bristlecone.ModelLibrary
module Likelihood from Bristlecone.ModelLibrary
<summary>Likelihood functions to represent a variety of distributions and data types.</summary>
<namespacedoc><summary>Pre-built model parts for use in Bristlecone</summary></namespacedoc>
val sumOfSquares: keys: ModelSystem.LikelihoodRequirement list -> ModelSystem.Likelihood<'u>
<summary> Residual sum of squares. Provides a simple metric of distance between observed data and model predictions. </summary>
module Require from Bristlecone.Language
val state: s: StateId<'u> -> ModelSystem.LikelihoodRequirement
module Components from Bristlecone.Language
val K: IncludedParameter<kilogram>
val notNegative: Parameter.Constraint
val growthLimit: ModelComponent<(ModelExpression<kilogram> -> ModelExpression<kilogram>)>
val modelComponent: name: string -> list: SubComponent<'a> list -> ModelComponent<'a>
val subComponent: name: string -> transform: 'a -> SubComponent<'a>
<summary> Creates a nested component that can be inserted into a base model. </summary>
val Constant: x: float<'u> -> ModelExpression<'u>
val m: ModelExpression<kilogram>
type ModelExpression<'u> = private | ME of ModelExpressionUntyped static member (%) : ModelExpression<'u> * ModelExpression<'u> -> ModelExpression<'u> static member ( * ) : ModelExpression<'u1> * ModelExpression<'u2> -> ModelExpression<'u1 'u2> static member (+) : ModelExpression<'u> * ModelExpression<'u> -> ModelExpression<'u> static member (-) : ModelExpression<'u> * ModelExpression<'u> -> ModelExpression<'u> static member (.<) : ModelExpression<'u> * ModelExpression<'u> -> BoolExpression static member (.>) : ModelExpression<'u> * ModelExpression<'u> -> BoolExpression static member (/) : ModelExpression<'u1> * ModelExpression<'u2> -> ModelExpression<'u1/'u2> static member (=.) : ModelExpression<'u> * ModelExpression<'u> -> BoolExpression static member Environment: sid: StateId<'u> -> ModelExpression<'u> static member Exp: ModelExpression<1> -> ModelExpression<1> ...
val estimateParameter: p: IncludedParameter<'u> -> comp: SubComponent<'a> -> SubComponent<'a>
val Logarithm: expr: ModelExpression<1> -> ModelExpression<1>
val alpha: IncludedParameter</Time.month>
val lossRate: ModelComponent<(ModelExpression<kilogram> -> ModelExpression<kilogram/Time.month>)>
val hypotheses: Hypotheses.Hypothesis<Time.month> list
module Hypotheses from Bristlecone.Language
<summary>Types to represent a hypothesis, given that a hypothesis is a model system that contains some alternate formulations of certain components.</summary>
val createFromModel: baseModel: ('a -> 'rest) -> Writer.Writer<('a -> 'rest),(Hypotheses.ComponentName * CodedMap<Parameter.Pool.AnyParameter>)> list
val apply: comp: ModelComponent<'a> -> builders: Writer.Writer<('a -> 'rest),(Hypotheses.ComponentName * CodedMap<Parameter.Pool.AnyParameter>)> list -> Writer.Writer<'rest,(Hypotheses.ComponentName * CodedMap<Parameter.Pool.AnyParameter>)> list
val compile: builders: Writer.Writer<ModelBuilder.ModelBuilder<'u>,(Hypotheses.ComponentName * CodedMap<Parameter.Pool.AnyParameter>)> list -> Hypotheses.Hypothesis<'u> list
<summary> Compile: run the writer(s), add parameters into the model builder, and wrap in Hypothesis </summary>

Type something to start searching.