Use parameter value accessor in likelihood functions

src/Bristlecone/Equations.fs

@@ -1,5 +1,5 @@
 /// Pre-built model parts for use in Bristlecone.
-module ModelLibrary
+namespace Bristlecone.ModelLibrary
 
 /// Likelihood functions to represent a variety of distributions
 /// and data types.
@@ -18,11 +18,14 @@ module Likelihood =
     let internal sumOfSquares' exp obs =
         Array.zip obs exp |> Array.sumBy (fun d -> ((fst d) - (snd d)) ** 2.)
 
-    let sumOfSquares keys _ (data: CodedMap<PredictedSeries>) =
-        keys
-        |> List.sumBy (fun k ->
-            let d = data |> getData k
-            sumOfSquares' d.Expected d.Observed)
+    /// Residual sum of squares. Provides a simple metric of distance between
+    /// observed data and model predictions.
+    let sumOfSquares keys : LikelihoodFn =
+        fun _ (data: CodedMap<PredictedSeries>) ->
+            keys
+            |> List.sumBy (fun k ->
+                let d = data |> getData k
+                sumOfSquares' d.Expected d.Observed)
 
     /// Negative log likelihood for a bivariate normal distribution.
     /// For two random variables with bivariate normal N(u1,u2,sigma1,sigma2,rho).
@@ -32,13 +35,15 @@ module Likelihood =
 
     /// <summary>
     /// Log likelihood function for single equation system, assuming Gaussian error for x.
+    /// Requires a parameter 'σ[x]' to be included in any `ModelSystem` that uses it.
     /// </summary>
-    let gaussian key pool data =
-        let x = data |> getData key
-        let sigmax = pool |> Parameter.Pool.tryGetTransformedValue "σ[x]" |> Option.get
+    let gaussian key : LikelihoodFn =
+        fun paramAccessor data ->
+            let x = data |> getData key
+            let sigmax = paramAccessor.Get "σ[x]"
 
-        [ 1 .. (Array.length x.Observed) - 1 ]
-        |> List.sumBy (fun i -> (gaussian' sigmax x.Observed.[i] x.Expected.[i]))
+            [ 1 .. (Array.length x.Observed) - 1 ]
+            |> List.sumBy (fun i -> (gaussian' sigmax x.Observed.[i] x.Expected.[i]))
 
     /// Negative log likelihood for a bivariate normal distribution.
     /// For two random variables with bivariate normal N(u1,u2,sigma1,sigma2,rho).
@@ -54,14 +59,16 @@ module Likelihood =
 
     /// <summary>
     /// Log likelihood function for dual simultaneous system, assuming Gaussian error for both x and y.
+    /// Requires parameters 'σ[x]', 'σ[y]' and 'ρ' to be included in any `ModelSystem` that uses it.
     /// </summary>
-    let bivariateGaussian key1 key2 pool data =
-        let x = data |> getData key1
-        let y = data |> getData key2
-        let sigmax = pool |> Parameter.Pool.tryGetTransformedValue "σ[x]" |> Option.get
-        let sigmay = pool |> Parameter.Pool.tryGetTransformedValue "σ[y]" |> Option.get
-        let rho = pool |> Parameter.Pool.tryGetTransformedValue "ρ" |> Option.get
+    let bivariateGaussian key1 key2 : LikelihoodFn =
+        fun paramAccessor data ->
+            let x = data |> getData key1
+            let y = data |> getData key2
+            let sigmax = paramAccessor.Get "σ[x]"
+            let sigmay = paramAccessor.Get "σ[y]"
+            let rho = paramAccessor.Get "ρ"
 
-        [ 1 .. (Array.length x.Observed) - 1 ]
-        |> List.sumBy (fun i ->
-            (bivariateGaussian' sigmax sigmay rho x.Observed.[i] y.Observed.[i] x.Expected.[i] y.Expected.[i]))
+            [ 1 .. (Array.length x.Observed) - 1 ]
+            |> List.sumBy (fun i ->
+                (bivariateGaussian' sigmax sigmay rho x.Observed.[i] y.Observed.[i] x.Expected.[i] y.Expected.[i]))

src/Bristlecone/EstimationEngine.fs

@@ -19,19 +19,22 @@ module ModelSystem =
     type PredictedSeries =
         { Expected: float[]; Observed: float[] }
 
+    /// A function that returns a parameter's current value by its name.
+    type ParameterValueAccessor = ParameterValueAccessor of (string -> float)
+        with member this.Get name = let (ParameterValueAccessor v) = this in v name
+
     /// A function that computes the likelihood of a set of parameters.
-    type Likelihood = Parameter.Pool -> CodedMap<PredictedSeries> -> float
+    type LikelihoodFn = ParameterValueAccessor -> CodedMap<PredictedSeries> -> float
 
     /// A function that computes a measured system property given a
     /// current (time t) and previous (time t-1) system state.
     type MeasureEquation = float -> Environment -> Environment -> float
 
-    /// TODO Constrain model system so that codes cannot be duplicated
     type ModelSystem =
         { Parameters: Parameter.Pool
           Equations: CodedMap<ModelEquation>
           Measures: CodedMap<MeasureEquation>
-          Likelihood: Likelihood }
+          Likelihood: LikelihoodFn }
 
     type FitValue = { Fit: float; Obs: float }
     type FitSeries = TimeSeries<FitValue>

src/Bristlecone/Language.fs

@@ -181,7 +181,7 @@ module Language =
         type ModelFragment =
             | EquationFragment of ModelExpression
             | ParameterFragment of Parameter.Parameter
-            | LikelihoodFragment of ModelSystem.Likelihood
+            | LikelihoodFragment of ModelSystem.LikelihoodFn
             | MeasureFragment of ModelSystem.MeasureEquation
 
         type ModelBuilder = private ModelBuilder of Map<ShortCode.ShortCode, ModelFragment>

src/Bristlecone/Objective.fs

@@ -47,4 +47,7 @@ module Objective =
             point
             |> predict system integrate solveDiscrete
             |> pairObservationsToExpected observed
-            |> system.Likelihood(point |> Parameter.Pool.fromPointInTransformedSpace system.Parameters)
+            |> system.Likelihood(
+                    let pool = point |> Parameter.Pool.fromPointInTransformedSpace system.Parameters
+                    ParameterValueAccessor
+                    <| fun name -> Parameter.Pool.tryGetRealValue name pool |> Option.get)

src/Bristlecone/Parameter.fs

@@ -137,18 +137,13 @@ module Parameter =
             |> Option.map snd
             |> Option.map getTransformedValue
 
-        let private resultToOption r =
-            match r with
-            | Ok x -> Some x
-            | Error _ -> None
-
         /// Gets the 'real' / non-transformed value for use in model
         /// calculation.
         let internal tryGetRealValue key (pool: ParameterPool) : float option =
             pool
             |> unwrap
             |> Map.tryFindBy (fun k -> k.Value = key)
-            |> Option.bind (getEstimate >> resultToOption)
+            |> Option.bind (getEstimate >> Result.toOption)
 
         /// Returns the starting bounds in transformed parameter space if
         /// the parameter has not been estimated. If the parameter has already
@@ -188,7 +183,8 @@ module Parameter =
             if count pool = (point |> Array.length) then
                 pool
                 |> toList
-                |> List.mapi (fun i (sc, p) -> sc, setTransformedValue p point.[i])
+                |> List.mapi (fun i (sc, p) -> 
+                    sc, setTransformedValue p point.[i])
                 |> List.choose (fun (c, r) ->
                     match r with
                     | Ok x -> Some(c, x)

src/Bristlecone/Types.fs

@@ -295,6 +295,13 @@ module Result =
         | Some o -> Ok o
         | None -> Error errorMessage
 
+    /// If OK, return Some, otherwise None.
+    let toOption r =
+        match r with
+        | Ok x -> Some x
+        | Error _ -> None
+
+
 [<AutoOpen>]
 module ListExtensions =
 

tests/Bristlecone.Benchmark/Bristlecone.Benchmark.fsproj

@@ -4,6 +4,7 @@
     <OutputType>Exe</OutputType>
     <TargetFramework>net6.0</TargetFramework>
     <IsPackable>false</IsPackable>
+    <IsTestProject>true</IsTestProject>
   </PropertyGroup>
   <ItemGroup>
     <ProjectReference Include="../../src/Bristlecone/Bristlecone.fsproj" />

tests/Bristlecone.Benchmark/TestFunctions.fs

@@ -97,13 +97,13 @@ module Timeseries =
 
     let ``predator-prey`` =
         predatorPreyBase
-        |> Model.useLikelihoodFunction (ModelLibrary.Likelihood.sumOfSquares [ "hare"; "lynx" ])
+        |> Model.useLikelihoodFunction (Bristlecone.ModelLibrary.Likelihood.sumOfSquares [ "hare"; "lynx" ])
         |> Model.compile
 
     let ``predator-prey [with noise]`` =
         predatorPreyBase
         |> Model.estimateParameter "σ[x]" notNegative 0.01 0.5
         |> Model.estimateParameter "σ[y]" notNegative 0.01 0.5
         |> Model.estimateParameter "ρ" noConstraints -0.5 0.5
-        |> Model.useLikelihoodFunction (ModelLibrary.Likelihood.bivariateGaussian "hare" "lynx")
+        |> Model.useLikelihoodFunction (Bristlecone.ModelLibrary.Likelihood.bivariateGaussian "hare" "lynx")
         |> Model.compile

tests/Bristlecone.Tests/Bristlecone.Tests.fsproj

@@ -4,6 +4,7 @@
     <OutputType>Exe</OutputType>
     <TargetFramework>net6.0</TargetFramework>
     <IsPackable>false</IsPackable>
+    <IsTestProject>true</IsTestProject>
   </PropertyGroup>
   <ItemGroup>
     <ProjectReference Include="../../src/Bristlecone/Bristlecone.fsproj" />

tests/Bristlecone.Tests/Language.fs

@@ -144,7 +144,7 @@ let modelBuilder =
         [
 
           testProperty "Does not compile when more than one likelihood function"
-          <| fun (likelihoodFns: ModelSystem.Likelihood list) ->
+          <| fun (likelihoodFns: ModelSystem.LikelihoodFn list) ->
               let mb =
                   likelihoodFns
                   |> Seq.fold (fun mb l -> mb |> Model.useLikelihoodFunction l) Model.empty