diff --git a/build.fsx b/build.fsx
index e145c59..7ed2668 100644
--- a/build.fsx
+++ b/build.fsx
@@ -31,13 +31,16 @@ open Fake.DotNet
let projectName = "Bristlecone"
let projectSummary = "Time-series modelling in F#"
-let projectDescription = """
+
+let projectDescription =
+ """
An F# library for model-fitting model-selection (MFMS) of ecological
models to observational data. The library was developed for
individual-based plant modelling using tree ring time-series,
but can be used for any ecological models."""
+
let authors = "Andrew Martin"
-let companyName = ""
+let companyName = "University of Oxford; University of Cambridge"
let tags = "modelling time-series F# fsharp R data-science ecology model-fitting"
let license = "MIT"
let iconUrl = ""
@@ -55,40 +58,55 @@ let repositoryContentUrl = "https://raw.githubusercontent.com/AndrewIOM/bristlec
// Read release notes & version info from RELEASE_NOTES.md
System.Environment.CurrentDirectory <- __SOURCE_DIRECTORY__
let binDir = __SOURCE_DIRECTORY__ @@ "bin"
-let release = System.IO.File.ReadLines "RELEASE_NOTES.MD" |> Fake.Core.ReleaseNotes.parse
+
+let release =
+ System.IO.File.ReadLines "RELEASE_NOTES.MD" |> Fake.Core.ReleaseNotes.parse
// Generate assembly info files with the right version & up-to-date information
Target.create "AssemblyInfo" (fun _ ->
- let fileName = "src/Bristlecone/AssemblyInfo.fs"
- AssemblyInfoFile.createFSharpWithConfig fileName
- [ Fake.DotNet.AssemblyInfo.Title projectName
- Fake.DotNet.AssemblyInfo.Company companyName
- Fake.DotNet.AssemblyInfo.Product projectName
- Fake.DotNet.AssemblyInfo.Description projectSummary
- Fake.DotNet.AssemblyInfo.Version release.AssemblyVersion
- Fake.DotNet.AssemblyInfo.FileVersion release.AssemblyVersion ]
- (AssemblyInfoFileConfig(false))
-)
+ let fileName = "src/Bristlecone/AssemblyInfo.fs"
+
+ AssemblyInfoFile.createFSharpWithConfig
+ fileName
+ [ Fake.DotNet.AssemblyInfo.Title projectName
+ Fake.DotNet.AssemblyInfo.Company companyName
+ Fake.DotNet.AssemblyInfo.Product projectName
+ Fake.DotNet.AssemblyInfo.Description projectSummary
+ Fake.DotNet.AssemblyInfo.Version release.AssemblyVersion
+ Fake.DotNet.AssemblyInfo.FileVersion release.AssemblyVersion ]
+ (AssemblyInfoFileConfig(false)))
// --------------------------------------------------------------------------------------
// Clean build results & restore NuGet packages
Target.create "Clean" (fun _ ->
- Fake.IO.Shell.cleanDirs ["bin"; "output" ]
- Fake.IO.Shell.cleanDirs ["tests/Bristlecone.Tests/bin"; "tests/Bristlecone.Tests/obj" ]
-)
+ Fake.IO.Shell.cleanDirs [ "bin"; "output" ]
+ Fake.IO.Shell.cleanDirs [ "tests/Bristlecone.Tests/bin"; "tests/Bristlecone.Tests/obj" ])
+
+Target.create "CleanDocs" (fun _ -> Fake.IO.Shell.cleanDirs [ ".fsdocs" ])
+
+// --------------------------------------------------------------------------------------
+// Check formatting with Fantomas
+
+Target.create "CheckFormat" (fun _ ->
+ let result = DotNet.exec id "fantomas" "./src --check"
+ let resultDocs = DotNet.exec id "fantomas" "./docs --check"
-Target.create "CleanDocs" (fun _ ->
- Fake.IO.Shell.cleanDirs [".fsdocs"]
-)
+ if result.ExitCode = 0 && resultDocs.ExitCode = 0 then
+ Trace.log "No files need formatting"
+ elif result.ExitCode = 99 then
+ failwith "Some files need formatting, run \"dotnet fantomas ./src\" to resolve this."
+ elif result.ExitCode = 99 then
+ failwith "Some files need formatting, run \"dotnet fantomas ./docs\" to resolve this."
+ else
+ Trace.logf "Errors while formatting: %A" result.Errors)
// --------------------------------------------------------------------------------------
// Build library & test project
Target.create "Build" (fun _ ->
Trace.log " --- Building the app --- "
- Fake.DotNet.DotNet.build id ("bristlecone.sln")
-)
+ Fake.DotNet.DotNet.build id ("bristlecone.sln"))
// --------------------------------------------------------------------------------------
// Run the unit tests using test runner & kill test runner when complete
@@ -96,51 +114,60 @@ Target.create "Build" (fun _ ->
Target.create "RunTests" (fun _ ->
let rHome = Environment.environVarOrFail "R_HOME"
Trace.logf "R_HOME is set as %s" rHome
- let result = Fake.DotNet.DotNet.exec (fun args ->
- { args with Verbosity = Some Fake.DotNet.DotNet.Verbosity.Normal}) "run" ("--project tests/Bristlecone.Tests")
- if result.ExitCode <> 0 then failwith "Tests failed"
-)
+
+ let result =
+ Fake.DotNet.DotNet.exec
+ (fun args ->
+ { args with
+ Verbosity = Some Fake.DotNet.DotNet.Verbosity.Normal })
+ "run"
+ ("--project tests/Bristlecone.Tests")
+
+ if result.ExitCode <> 0 then
+ failwith "Tests failed")
// --------------------------------------------------------------------------------------
// Build a NuGet package
Target.create "NuGet" (fun _ ->
// Format the description to fit on a single line (remove \r\n and double-spaces)
- let projectDescription = projectDescription.Replace("\r", "").Replace("\n", "").Replace(" ", " ")
-
+ let projectDescription =
+ projectDescription.Replace("\r", "").Replace("\n", "").Replace(" ", " ")
+
// Format the release notes
let releaseNotes = release.Notes |> String.concat "\n"
- let properties = [
- ("Version", release.NugetVersion)
- ("Authors", authors)
- ("PackageProjectUrl", packageProjectUrl)
- ("PackageTags", tags)
- ("RepositoryType", repositoryType)
- ("RepositoryUrl", repositoryUrl)
- ("PackageLicenseExpression", license)
- ("PackageRequireLicenseAcceptance", "false")
- ("PackageReleaseNotes", releaseNotes)
- ("Summary", projectSummary)
- ("PackageDescription", projectDescription)
- ("PackageIcon", "logo.png")
- ("PackageIconUrl", iconUrl)
- ("EnableSourceLink", "true")
- ("PublishRepositoryUrl", "true")
- ("EmbedUntrackedSources", "true")
- ("IncludeSymbols", "true")
- ("IncludeSymbols", "false")
- ("SymbolPackageFormat", "snupkg")
- ("Copyright", copyright)
- ]
-
- DotNet.pack (fun p ->
- { p with
- Configuration = DotNet.BuildConfiguration.Release
- OutputPath = Some "bin"
- MSBuildParams = { p.MSBuildParams with Properties = properties}
- }
- ) ("bristlecone.sln"))
+ let properties =
+ [ ("Version", release.NugetVersion)
+ ("Authors", authors)
+ ("PackageProjectUrl", packageProjectUrl)
+ ("PackageTags", tags)
+ ("RepositoryType", repositoryType)
+ ("RepositoryUrl", repositoryUrl)
+ ("PackageLicenseExpression", license)
+ ("PackageRequireLicenseAcceptance", "false")
+ ("PackageReleaseNotes", releaseNotes)
+ ("Summary", projectSummary)
+ ("PackageDescription", projectDescription)
+ ("PackageIcon", "logo.png")
+ ("PackageIconUrl", iconUrl)
+ ("EnableSourceLink", "true")
+ ("PublishRepositoryUrl", "true")
+ ("EmbedUntrackedSources", "true")
+ ("IncludeSymbols", "true")
+ ("IncludeSymbols", "false")
+ ("SymbolPackageFormat", "snupkg")
+ ("Copyright", copyright) ]
+
+ DotNet.pack
+ (fun p ->
+ { p with
+ Configuration = DotNet.BuildConfiguration.Release
+ OutputPath = Some "bin"
+ MSBuildParams =
+ { p.MSBuildParams with
+ Properties = properties } })
+ ("bristlecone.sln"))
//--------------------------------------------------------------------------------------
//Generate the documentation
@@ -163,21 +190,19 @@ Target.create "DocsMeta" (fun _ ->
"true"
"default"
""
- ""]
- |> Fake.IO.File.write false "Directory.Build.props"
-)
+ "" ]
+ |> Fake.IO.File.write false "Directory.Build.props")
Target.create "GenerateDocs" (fun _ ->
- Fake.IO.Shell.cleanDir ".fsdocs"
- DotNet.exec id "fsdocs" "build --clean --eval --strict" |> ignore
-)
+ Fake.IO.Shell.cleanDir ".fsdocs"
+ DotNet.exec id "fsdocs" "build --clean --eval --strict" |> ignore)
// --------------------------------------------------------------------------------------
// Run all targets by default. Invoke 'build ' to override
Target.create "All" ignore
-"Clean" ==> "AssemblyInfo" ==> "Build"
+"Clean" ==> "CheckFormat" ==> "AssemblyInfo" ==> "Build"
"Build" ==> "CleanDocs" ==> "DocsMeta" ==> "GenerateDocs" ==> "All"
"Build" ==> "NuGet" ==> "All"
"Build" ==> "RunTests" ==> "All"
diff --git a/components/allometric.fsx b/components/allometric.fsx
deleted file mode 100644
index 0568638..0000000
--- a/components/allometric.fsx
+++ /dev/null
@@ -1,125 +0,0 @@
-module Allometric
-
-#load "../bristlecone.fsx"
-
-open Bristlecone
-
-let pi = System.Math.PI
-
-module Constants =
-
- // Empirically-derived parameters:
- let k5 = 19.98239 // Allometric fit to Yamal shrub BD-length data #1 (in centimetres)
- let k6 = 0.42092 // Allometric fit to Yamal shrub BD-length data #2 (in centimetres)
-
- // Constants from the literature:
- let a = 2. // the number of child branches added to previous branches (including the tops of the stems) for a shrub
- let p = 0.5 // the length of a child branch as a proportion of its parent branch/stem
- let lmin = 20. //cm. the length at which a stem or branch gets child branches
- let rtip = 0.1 //cm. the radius of the outermost tip of a stem or branch
- let b = 0.0075 // the ratio of the basal radius of a stem or branch and its length
- let salixWoodDensity = 0.5 // g / cm3 (from internet)
- let numberOfStems = 2.2
-
-module NiklasAndSpatz_Allometry =
-
- let nthroot n A =
- let rec f x =
- let m = n - 1.
- let x' = (m * x + A/x**m) / n
- match abs(x' - x) with
- | t when t < abs(x * 1e-9) -> x'
- | _ -> f x'
- f (A / double n)
-
- /// Gives the basal radius in centimetres of a stem/branch given its length in centimetres. Function from Niklas and Spatz (2004).
- /// The basal radius is always positive.
- let basalRadius k5 k6 stemLength =
- max (100. * (( 0.01 * stemLength + k6) / k5) ** (3. / 2.) / 2.) 1e-06
-
- /// Inverse equation of basalRadius.
- let stemLength k5 k6 radius =
- max (2. * ((nthroot 3. 2.) * 5. ** (2./3.) * k5 * radius ** (2./3.) - 50. * k6)) 1e-06
-
-module Götmark2016_ShrubModel =
-
- /// Total shrub volume given height and number of stems
- let shrubVolume b a rtip p lmin k5 k6 n h =
-
- let radius = NiklasAndSpatz_Allometry.basalRadius k5 k6
- let mainStemVolume =
- match radius h with
- | r when r > rtip -> n * pi * h * ((radius h) ** 2. + (radius h) * rtip + rtip ** 2.) / 3.
- | _ -> n * pi * h * rtip ** 2.
-
- let mutable volume = mainStemVolume
- let mutable k = 0.
-
- while (p ** k * h > lmin * 2./3.) do
- let volToAdd =
- match (p ** k * h < lmin) with
- | true ->
- match (b * 3. * p * (p ** k * h - 2. * lmin / 3.) > rtip) with
- | true ->
- n * a * (a + 1.) ** (float k) * pi * 3. * p * (p ** k * h - 2. * lmin / 3.) * ((radius (3. * p * (p ** k * h - 2. * lmin / 3.))) * (radius (3. * p * (p ** k * h - 2. * lmin / 3.)) * rtip + rtip ** 2.)) / 3.
- | false ->
- n * a * (a + 1.) ** (float k) * 3. * p * (p ** k * h - 2. * lmin / 3.) * pi * rtip ** 2.
- | false ->
- match (radius (p ** (k + 1.) * h) > rtip) with
- | true ->
- n * a * (a + 1.) ** (float k) * pi * p ** (k + 1.) * h * ((radius (p ** (k+1.) * h)) ** 2. + (radius (p ** (k + 1.) * h)) * rtip + rtip ** 2.) / 3.
- | false ->
- n * a * (a + 1.) ** (float k) * p ** (k + 1.) * h * pi * rtip ** 2.
-
- volume <- volume + volToAdd
- k <- k + 1.
-
- k, volume
-
-
-module Allometrics =
-
- open Götmark2016_ShrubModel
- open Bristlecone.Dendro
-
- let private removeUnit (x:float<_>) = float x
-
- let mass woodDensity volume =
- volume * woodDensity
-
- let massToVolume woodDensity mass =
- mass / woodDensity
-
- let shrubBiomass b a rtip p lmin k5 k6 n woodDensity (radius:float) =
- radius
- |> removeUnit
- |> NiklasAndSpatz_Allometry.stemLength k5 k6
- |> shrubVolume b a rtip p lmin k5 k6 n |> snd
- |> mass woodDensity
-
- let shrubRadius b a rtip p lmin k5 k6 n woodDensity mass =
- let findRadius volume =
- let v x = x |> NiklasAndSpatz_Allometry.stemLength k5 k6 |> shrubVolume b a rtip p lmin k5 k6 n |> snd
- let f = (fun x -> (v x) - volume )
- Statistics.RootFinding.bisect 0 200 f 0.01 100.00 1e-8 // Assumption that shrub radius is between 0.01 and 100.0cm.
- mass
- |> massToVolume woodDensity
- |> findRadius
-
- let shrubHeight k5 k6 radius =
- radius |> NiklasAndSpatz_Allometry.stemLength k5 k6
-
-module Proxies =
-
- open Bristlecone.Dendro
-
- /// Radius in millimetres
- let toBiomassMM (radius:float) =
- radius / 10. |> Allometrics.shrubBiomass Constants.b Constants.a Constants.rtip Constants.p Constants.lmin Constants.k5 Constants.k6 Constants.numberOfStems Constants.salixWoodDensity
-
- /// Biomass in grams.
- let toRadiusMM biomassGrams =
- if System.Double.IsNaN biomassGrams || System.Double.IsInfinity biomassGrams || System.Double.IsNegativeInfinity biomassGrams then nan
- else
- let radiusCm = biomassGrams |> Allometrics.shrubRadius Constants.b Constants.a Constants.rtip Constants.p Constants.lmin Constants.k5 Constants.k6 Constants.numberOfStems Constants.salixWoodDensity
- radiusCm * 10.
diff --git a/docs/composing-models.fsx b/docs/composing-models.fsx
new file mode 100644
index 0000000..861280c
--- /dev/null
+++ b/docs/composing-models.fsx
@@ -0,0 +1,133 @@
+(**
+---
+title: Composing models into many hypotheses
+category: The Bristlecone Language
+categoryindex: 1
+index: 4
+---
+*)
+
+(*** condition: prepare ***)
+#nowarn "211"
+#r "../src/Bristlecone/bin/Debug/netstandard2.0/Bristlecone.dll"
+(*** condition: fsx ***)
+#if FSX
+#r "nuget: Bristlecone,{{fsdocs-package-version}}"
+#endif // FSX
+(*** condition: ipynb ***)
+#if IPYNB
+#r "nuget: Bristlecone,{{fsdocs-package-version}}"
+#endif // IPYNB
+
+(**
+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
+
+(**
+
+
+
+### 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. For example,
+take the following model:
+*)
+
+let someModel someFn =
+ Parameter "K" + someFn This - Environment "T" * Parameter "J"
+
+(**
+In the above expression, `someFn` is a function of signature `ModelExpression -> ModelExpression`.
+Looking at the equation, you can see that `someFn` takes one argument, the model expression
+`This`, which represents the current state of the model equation.
+
+We can take this concept further to 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 baseModel growthLimit lossRate =
+ Model.empty
+ |> Model.addEquation "m" (Parameter "r" * (growthLimit This) - lossRate This)
+ |> Model.estimateParameter "r" noConstraints 0.01 1.00
+ |> Model.useLikelihoodFunction (ModelLibrary.Likelihood.sumOfSquares [ "x" ])
+
+(**
+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:
+*)
+
+let growthLimit =
+ modelComponent
+ "Limitation to growth"
+ [
+
+ subComponent "Linear" (fun _ -> Constant 1.)
+
+ subComponent "Monomolecular" (fun m -> Parameter "K" - m)
+ |> estimateParameter "K" notNegative 10.0 50.0
+
+ subComponent "Gompertz" (fun m -> m * log (Parameter "K" / m))
+ |> estimateParameter "K" notNegative 10.0 50.0
+
+ subComponent "Logisitc (three parameter)" (fun m -> m * (Constant 1. - (m / Parameter "K")))
+ |> estimateParameter "K" notNegative 10.0 50.0 ]
+
+(**
+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 lossRate =
+ modelComponent
+ "Biomass loss rate"
+ [
+
+ subComponent "None" (fun _ -> Constant 0.)
+
+ subComponent "Density-dependent" (fun m -> m * Parameter "alpha")
+ |> estimateParameter "alpha" notNegative 0.001 0.10 ]
+
+(**
+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.createFromComponent growthLimit
+ |> Hypotheses.useAnother lossRate
+ |> Hypotheses.compile
+
+(**
+The resultant constructed hypotheses are shown in the following printout:
+*)
+
+(***include-value: hypotheses ***)
+
+(**
+For further analysis, you may choose to use the orchestration functionality to
+run many hypotheses and replicates (multiple fits per hypothesis) in parallel.
+*)
diff --git a/docs/confidence-intervals.fsx b/docs/confidence-intervals.fsx
index 55f67c1..8b8d6a5 100644
--- a/docs/confidence-intervals.fsx
+++ b/docs/confidence-intervals.fsx
@@ -37,20 +37,20 @@ open Bristlecone
open Bristlecone.Optimisation
open Bristlecone.Data
-// // fitFn =
-// fun engine dataset hypothesis result ->
+fun engine dataset hypothesis result ->
-// // The function used to fit the model, which unless an
-// // advanced scenario is usually Bristlecone.fit
-// let fitFn = Bristlecone.fit
+ // The function used to fit the model, which unless an
+ // advanced scenario is usually Bristlecone.fit
+ let fitFn = Bristlecone.fit
-// // The number of jumps to perform in parameter space
-// let n = 10000
+ // The number of jumps to perform in parameter space
+ let n = 10000
-// let ci = ConfidenceInterval.ProfileLikelihood.profile fitFn engine dataset hypothesis n result
+ let ci =
+ Confidence.ProfileLikelihood.profile fitFn engine dataset hypothesis n result
-// // Optionally, save the result
-// let saveDir = "/some/save/dir"
-// let subjectName = "some subject"
-// let modelId = "some model hypothesis"
-// Confidence.save saveDir subjectName modelId result.ResultId ci
\ No newline at end of file
+ // Optionally, save the result
+ let saveDir = "/some/save/dir"
+ let subjectName = "some subject"
+ let modelId = "some model hypothesis"
+ Confidence.save saveDir subjectName modelId result.ResultId ci
diff --git a/docs/data.fsx b/docs/data.fsx
index 842d8b6..65ae137 100644
--- a/docs/data.fsx
+++ b/docs/data.fsx
@@ -5,6 +5,9 @@ category: Techniques
categoryindex: 2
index: 6
---
+
+[]({{root}}/{{fsdocs-source-basename}}.fsx)
+[]({{root}}/{{fsdocs-source-basename}}.ipynb)
*)
(*** condition: prepare ***)
@@ -27,38 +30,87 @@ The `Bristlecone.Data` namespace includes methods for saving Bristlecone
results and diagnostics, and loading saved results for further analysis
at a later date.
-### Estimation Results
+### Individual estimations results
+
+#### The estimation result itself
Saving and loading estimation results is conducted as follows:
*)
open Bristlecone
open Bristlecone.Data
+open Bristlecone.Language
let resultsDirectory = "/some/data/dir"
-let thinTraceBy = 1000 // Only trace every n iterations to save disk space
+let thinTraceBy = Some 1000 // Only trace every n iterations to save disk space
let subject = "Mosquito population" // The name of the 'subject' of analysis
let modelName = "Logistic growth"
-/// Output three files: the maximum likelihood estimate,
-/// a trace of the optimisation process, and the estimated vs
-/// observed time series.
-fun result ->
- EstimationResult.saveAll resultsDirectory subject modelName thinTraceBy result
+fun result -> EstimationResult.saveAll resultsDirectory subject modelName thinTraceBy result
(**
-### Other statistics with save functions
+This save function outputs three files: one each for the maximum likelihood estimate,
+a trace of the optimisation process, and the estimated vs
+observed time series.
-`Bristlecone.Data` includes methods for saving results from the following
-statistical functions that Bristlecone performs:
-*)
+#### Other statistics and estimates
+
+Some other Bristlecone functions have functions to load and save their outputs to CSV files.
-// Save confidence intervals:
+Save confidence intervals:
+*)
Bristlecone.Data.Confidence.save
-// Save convergence statistics based on multiple 'chains'
-// of analysis on the same subject and model:
-Bristlecone.Data.Convergence.save
+(**
+Save n-step ahead predictions made using `Bristlecone.oneStepAhead` and similar functions:
+*)
+
+Bristlecone.Data.NStepAhead.save
+
+(**
+### Ensemble-level results
+
+Some processes, for example model-selection, work across many model results. Bristlecone
+includes loading and saving functions (to and from CSV files) for many of these procedures.
+
+#### Model-selection
+
+For model selection, you may calculate Akaike weights either from a sequence of `EstimationResult` (for simpler tasks)
+or from a sequence of `ResultSet` (for working with many hypotheses):
+*)
+
+fun modelResults subjectIds hypothesisIds ->
+ modelResults
+ |> ModelSelection.Akaike.akaikeWeights
+ |> Seq.zip3 subjectIds hypothesisIds
+ |> Seq.map (fun (s, h, (a, b)) -> (s, h, a, b))
+ |> Bristlecone.Data.ModelSelection.save "/some/dir"
+
+fun (hypothesisResults: ModelSelection.ResultSet.ResultSet seq) ->
+ hypothesisResults
+ |> ModelSelection.Akaike.akaikeWeightsForSet (fun h -> h.ReferenceCode)
+ |> Bristlecone.Data.ModelSelection.save "/some/dir"
+
+(**
+#### Convergence of Monte-Carlo based optimisation approaches
+
+You can assess the convergence of multiple 'chains' between MCMC or Monte-Carlo based
+optimisation methods using the functions in the `Bristlecone.Diagnostics` namespace, such as
+the per-parameter R-hat values based on the optimisation trace.
+*)
+
+// Given a set of results on a per-subject and per-hypothesis basis:
+fun resultSet ->
+ resultSet
+ |> Bristlecone.Diagnostics.Convergence.gelmanRubin 10000 "some subject" "some hypothesis"
+ |> Option.iter (Bristlecone.Data.Convergence.save "/some/dir")
+
+(**
+#### Summarise Root Mean Squared Error for many hypotheses / subjects
+
+Given individual n-step ahead predictions made using `Bristlecone.oneStepAhead`
+and similar functions, you can save the resultant RMSE statistic (an indication
+of model performance) in a summary table across all hypotheses and subjects:
+*)
-// Save Akaike weights calculated for model selection:
-Bristlecone.Data.ModelSelection.save
\ No newline at end of file
+fun nStepFits -> Bristlecone.Data.NStepAhead.saveAllRMSE "/some/dir" nStepFits
diff --git a/docs/dendro.fsx b/docs/dendro.fsx
index 92e1aee..577f537 100644
--- a/docs/dendro.fsx
+++ b/docs/dendro.fsx
@@ -5,10 +5,14 @@ category: Bristlecone.Dendro
categoryindex: 3
index: 1
---
+
+[]({{root}}/{{fsdocs-source-basename}}.fsx)
+[]({{root}}/{{fsdocs-source-basename}}.ipynb)
*)
(*** condition: prepare ***)
#nowarn "211"
#r "../src/Bristlecone/bin/Debug/netstandard2.0/Bristlecone.dll"
+#r "../src/Bristlecone.Dendro/bin/Debug/netstandard2.0/Bristlecone.Dendro.dll"
(*** condition: fsx ***)
#if FSX
#r "nuget: Bristlecone,{{fsdocs-package-version}}"
@@ -19,20 +23,46 @@ index: 1
#endif // IPYNB
(**
-Bristlecone
-======================
-
-Full documentation is forthcoming.
-
-
-
-*)
\ No newline at end of file
+## Bristlecone as a tool in dendroecology
+
+Bristlecone has been developed and used for working with
+long-term ecological time-series, which includes wood ring
+proxy records.
+
+The `Bristlecone.Dendro` package contains additional functions
+that can be usesful when conducting top-down mechanistic modelling
+with tree- and shrub-ring data, such as ring width and stable isotope
+records. You can reference it in an F# script as follows:
+*)
+
+(*** do-not-eval ***)
+#r "nuget: Bristlecone.Dendro"
+
+(**
+### Seasonal Cycles and Day Length
+
+In high and mid-latitudes, seasonal cycles are the dominant control
+on plant productivity. In Bristlecone, we can represent seasonal variability
+in day length and environmental conditions - such as soil temperature
+and moisture - within model systems using some features of `Bristlecone.Dendro`.
+
+#### Sunrise, Sunset and Hours of Light
+
+To calculate whether there is total sun or darkness, or partial light, at any
+location on any day of the year, we can use the `Sunrise.calculate` function.
+*)
+
+open Bristlecone.Dendro
+
+let latitude = 54.2
+let longitude = -4.2
+
+Sunrise.calculate 2024 04 20 latitude longitude "Europe/London"
+(*** include-it ***)
+
+(**
+As an illustration, if we looked at Tromsø in January we would see far less available light:
+*)
+
+Sunrise.calculate 2024 01 10 69.64961 18.95702 "Europe/London"
+(*** include-it ***)
diff --git a/docs/diagnostics.fsx b/docs/diagnostics.fsx
index 6d8bf61..703bdc6 100644
--- a/docs/diagnostics.fsx
+++ b/docs/diagnostics.fsx
@@ -5,6 +5,9 @@ category: Techniques
categoryindex: 2
index: 5
---
+
+[]({{root}}/{{fsdocs-source-basename}}.fsx)
+[]({{root}}/{{fsdocs-source-basename}}.ipynb)
*)
(*** condition: prepare ***)
@@ -46,9 +49,10 @@ let saveDirectory = "/some/save/dir"
fun listOfResults ->
- // Calculate Gelman-Rubin statistic for all parameters
+ // Calculate Gelman-Rubin statistic for all parameters
// across multiple model runs
- let stat = Convergence.gelmanRubinAll 1000 listOfResults
+ let stat =
+ Convergence.gelmanRubinAll 1000 (fun _ -> "some subject") (fun _ -> "some hypothesis") listOfResults
// Save results to a single file
Convergence.save saveDirectory stat
@@ -66,7 +70,8 @@ obtain computed values across the time-series for components within
a model system. First, you must set up a model that has a parameter
of the `IComponentLogger)` type:
-[ NB TODO: This function must be refactored to work with the new Bristlecone Language ]
+[ NB TODO: This function must be refactored to work with the new v2 Bristlecone Language,
+ and as such is currently disabled. ]
*)
// open Bristlecone
@@ -75,12 +80,12 @@ of the `IComponentLogger)` type:
// let hypothesis (cLog:IComponentLogger) =
-// let ``dN/dt`` =
+// let ``dN/dt`` =
// Parameter "r" * This * cLog.StoreValue "log this thing" (Constant 1. - (This / Parameter "K"))
// Model.empty
// |> Model.addEquation "N" ``dN/dt``
-// |> Model.estimateParameter "r" noConstraints 0.10 5.00
+// |> Model.estimateParameter "r" noConstraints 0.10 5.00
// |> Model.estimateParameter "K" noConstraints 50.0 150.
// |> Model.useLikelihoodFunction (ModelLibrary.Likelihood.sumOfSquares [ "N" ])
// |> Model.compile
@@ -90,4 +95,4 @@ of the `IComponentLogger)` type:
// // A wrapper around the Bristlecone.fit function
// let fitFn ts m e = Bristlecone.fit e (Optimisation.EndConditions.afterIteration 1) ts m
// let logs = ModelComponents.calculateComponents fitFn engine singleResult
-// ()
\ No newline at end of file
+// ()
diff --git a/docs/estimation-engine.fsx b/docs/estimation-engine.fsx
index 90f8f57..ae5bfb6 100644
--- a/docs/estimation-engine.fsx
+++ b/docs/estimation-engine.fsx
@@ -5,6 +5,9 @@ category: The Bristlecone Language
categoryindex: 1
index: 3
---
+
+[]({{root}}/{{fsdocs-source-basename}}.fsx)
+[]({{root}}/{{fsdocs-source-basename}}.ipynb)
*)
(*** condition: prepare ***)
@@ -76,4 +79,4 @@ To use Bristlecone functions requires a configured ``EstimationEngine``. The eas
| ``withContinuousTime`` | t : Integrate<'a,'b> -> engine: EstimationEngine<'a,'b> -> EstimationEngine<'a,'b> | Transforms engine to continuous time mode, using the given integrator function. |
| ``withConditioning`` | c: Conditioning -> engine: EstimationEngine<'a,'b> -> EstimationEngine<'a,'b> | Choose how the start point is chosen when solving the model system. |
-**)
\ No newline at end of file
+**)
diff --git a/docs/examples/dendroecology.fsx b/docs/examples/dendroecology.fsx.todo
similarity index 98%
rename from docs/examples/dendroecology.fsx
rename to docs/examples/dendroecology.fsx.todo
index 985e9a5..4428909 100644
--- a/docs/examples/dendroecology.fsx
+++ b/docs/examples/dendroecology.fsx.todo
@@ -1,9 +1,9 @@
(**
---
-title: Dendroecology: investigating plant-environment interactions
+title: Investigating plant-environment interactions using wood rings
category: Examples
-categoryindex: 2
-index: 1
+categoryindex: 3
+index: 2
---
*)
diff --git a/docs/examples/predator-prey.fsx b/docs/examples/predator-prey.fsx
index b924cbd..5a59428 100644
--- a/docs/examples/predator-prey.fsx
+++ b/docs/examples/predator-prey.fsx
@@ -2,7 +2,7 @@
---
title: Predator-Prey Dynamics
category: Examples
-categoryindex: 1
+categoryindex: 3
index: 1
---
@@ -31,8 +31,8 @@ To get started, we first load and open the Bristlecone library in
an F# script file (.fsx):
*)
-open Bristlecone // Opens Bristlecone core library and estimation engine
-open Bristlecone.Language // Open the language for writing Bristlecone models
+open Bristlecone // Opens Bristlecone core library and estimation engine
+open Bristlecone.Language // Open the language for writing Bristlecone models
(**### Defining the ecological model
@@ -55,26 +55,24 @@ in lynx data, the variability in hare data, and their covariance.
let ``predator-prey`` =
let ``dh/dt`` = Parameter "α" * This - Parameter "β" * This * Environment "lynx"
- let ``dl/dt`` = - Parameter "γ" * This + Parameter "Δ" * Environment "hare" * This
+ let ``dl/dt`` = -Parameter "δ" * This + Parameter "γ" * Environment "hare" * This
Model.empty
- |> Model.addEquation "hare" ``dh/dt``
- |> Model.addEquation "lynx" ``dl/dt``
+ |> Model.addEquation "hare" ``dh/dt``
+ |> Model.addEquation "lynx" ``dl/dt``
- |> Model.estimateParameter "α" noConstraints 0.01 1.00 // Natural growth rate of hares in absence of predation
- |> Model.estimateParameter "β" noConstraints 0.01 1.00 // Death rate per encounter of hares due to predation
- |> Model.estimateParameter "Δ" noConstraints 0.01 0.20 // Efficiency of turning predated hares into lynx
- |> Model.estimateParameter "γ" noConstraints 0.01 0.20 // Natural death rate of lynx in the absence of food
+ |> Model.estimateParameter "α" noConstraints 0.75 1.25 // Natural growth rate of hares in absence of predation
+ |> Model.estimateParameter "β" noConstraints 0.01 0.20 // Death rate per encounter of hares due to predation
+ |> Model.estimateParameter "δ" noConstraints 0.75 1.25 // Natural death rate of lynx in the absence of food
+ |> Model.estimateParameter "γ" noConstraints 0.01 0.20 // Efficiency of turning predated hares into lynx
- |> Model.useLikelihoodFunction (ModelLibrary.Likelihood.sumOfSquares ["hare"; "lynx"])
- |> Model.estimateParameter "ρ" noConstraints -0.500 0.500
- |> Model.estimateParameter "σ[x]" notNegative 0.001 0.100
- |> Model.estimateParameter "σ[y]" notNegative 0.001 0.100
+ |> Model.useLikelihoodFunction (ModelLibrary.Likelihood.bivariateGaussian "hare" "lynx")
+ |> Model.estimateParameter "ρ" noConstraints -0.500 0.500
+ |> Model.estimateParameter "σ[x]" notNegative 0.001 0.100
+ |> Model.estimateParameter "σ[y]" notNegative 0.001 0.100
|> Model.compile
- // TODO Error when setting constraints. Optim allows them to go negative!
-
(**
### Setting up the *Bristlecone Engine*
@@ -84,9 +82,14 @@ This engine uses a gradident descent method (Nelder Mead simplex), and a basic
Runge-Kutta 4 integration method provided by MathNet Numerics.
*)
-let engine =
+let engine =
Bristlecone.mkContinuous
- |> Bristlecone.withTunedMCMC []
+ // |> Bristlecone.withCustomOptimisation (Optimisation.Amoeba.swarm 5 20 Optimisation.Amoeba.Solver.Default)
+ |> Bristlecone.withCustomOptimisation (
+ Optimisation.MonteCarlo.Filzbach.filzbach
+ { Optimisation.MonteCarlo.Filzbach.FilzbachSettings.Default with
+ BurnLength = Optimisation.EndConditions.afterIteration 10000 }
+ )
|> Bristlecone.withContinuousTime Integration.MathNet.integrate
|> Bristlecone.withConditioning Conditioning.RepeatFirstDataPoint
|> Bristlecone.withSeed 1000 // We are setting a seed for this example - see below
@@ -115,16 +118,16 @@ but here we will configure some additional settings:
let testSettings =
Test.create
- |> Test.addStartValues [
- "hare", 50.
- "lynx", 75. ]
+ |> Test.addStartValues [ "hare", 50.; "lynx", 75. ]
|> Test.addNoise (Test.Noise.tryAddNormal "σ[y]" "lynx")
|> Test.addNoise (Test.Noise.tryAddNormal "σ[x]" "hare")
- |> Test.addGenerationRules [
- Test.GenerationRules.alwaysLessThan 10000. "lynx"
- Test.GenerationRules.alwaysLessThan 10000. "hare" ]
+ |> Test.addGenerationRules
+ [ Test.GenerationRules.alwaysLessThan 100000. "lynx"
+ Test.GenerationRules.alwaysMoreThan 10. "lynx"
+ Test.GenerationRules.alwaysLessThan 100000. "hare"
+ Test.GenerationRules.alwaysMoreThan 10. "hare" ]
|> Test.withTimeSeriesLength 30
- |> Test.endWhen (Optimisation.EndConditions.afterIteration 1)
+ |> Test.endWhen (Optimisation.EndConditions.afterIteration 100)
(**
In our `TestSettings`, we have specified the initial time point (t = 0)
@@ -135,9 +138,7 @@ should be 30 years in length.
With these test settings, we can now run the test.
*)
-let testResult =
- ``predator-prey``
- |> Bristlecone.testModel engine testSettings
+let testResult = ``predator-prey`` |> Bristlecone.tryTestModel engine testSettings
(*** include-output ***)
(**
@@ -149,53 +150,53 @@ module Graphing =
open Plotly.NET
- let pairedFits (series:Map) =
+ let pairedFits (series: Map) =
match testResult with
| Ok r ->
series
- |> Seq.map(fun kv ->
- let lines =
+ |> Seq.map (fun kv ->
+ let lines =
kv.Value
|> Bristlecone.Time.TimeSeries.toObservations
- |> Seq.collect(fun (d,v) ->
- [
- v, "Modelled", d.Fit
- v, "Observed", d.Obs
- ])
- |> Seq.groupBy(fun (_,x,_) -> x)
- |> Seq.map(fun (_,s) -> s |> Seq.map(fun (x,_,y) -> x,y))
+ |> Seq.collect (fun (d, v) -> [ v, "Modelled", d.Fit; v, "Observed", d.Obs ])
+ |> Seq.groupBy (fun (_, x, _) -> x)
+ |> Seq.map (fun (_, s) -> s |> Seq.map (fun (x, _, y) -> x, y))
|> Seq.toList
- // Each chart has the modelled and observed series
- Chart.combine [ Chart.Line(xy = lines.[0], Name = "Modelled"); Chart.Line(xy = lines.[1], Name = "Observed") ]
- |> Chart.withTitle kv.Key
- )
- |> Chart.Grid(2,1)
- |> fun x -> printfn "%A" x; x
+ // Each chart has the modelled and observed series
+ Chart.combine
+ [ Chart.Line(xy = lines.[0], Name = "Modelled")
+ Chart.Line(xy = lines.[1], Name = "Observed") ]
+ |> Chart.withTitle kv.Key)
+ |> Chart.Grid(2, 1)
+ |> fun x ->
+ printfn "%A" x
+ x
|> GenericChart.toChartHTML
- |> fun x -> printfn "%s" x; x
+ |> fun x ->
+ printfn "%s" x
+ x
| Error e -> sprintf "Cannot display data, as model fit did not run successfully (%s)" e
- let pairedFitsForTestResult (testResult:Result) =
+ let pairedFitsForTestResult (testResult: Result) =
match testResult with
| Ok r -> pairedFits r.Series
| Error e -> sprintf "Cannot display data, as model fit did not run successfully (%s)" e
- let pairedFitsForResult (testResult:Result) =
+ let pairedFitsForResult (testResult: Result) =
match testResult with
- | Ok r -> pairedFits (r.Series |> Seq.map(fun kv -> kv.Key.Value, kv.Value) |> Map.ofSeq)
+ | Ok r -> pairedFits (r.Series |> Seq.map (fun kv -> kv.Key.Value, kv.Value) |> Map.ofSeq)
| Error e -> sprintf "Cannot display data, as model fit did not run successfully (%s)" e
- let parameterTrace (result:Result) =
+ let parameterTrace (result: Result) =
match result with
- | Ok r ->
+ | Ok r ->
r.Trace
|> Seq.map snd
|> Seq.map Seq.toList
|> Seq.toList
|> List.flip
- |> List.map(fun values ->
- Chart.Line(y = values, x = [ 1 .. values.Length ]))
- |> Chart.Grid(3,3)
+ |> List.map (fun values -> Chart.Line(y = values, x = [ 1 .. values.Length ]))
+ |> Chart.Grid(3, 3)
|> GenericChart.toChartHTML
| Error _ -> "Model did not fit successfully"
@@ -217,12 +218,14 @@ a Bristlecone `TimeSeries` type using `TimeSeries.fromObservations`:
[]
let ResolutionFolder = __SOURCE_DIRECTORY__
-type PopulationData = FSharp.Data.CsvProvider<"data/lynx-hare.csv", ResolutionFolder = ResolutionFolder>
+type PopulationData = FSharp.Data.CsvProvider<"data/lynx-hare.csv", ResolutionFolder=ResolutionFolder>
-let data =
- let csv = PopulationData.Load (__SOURCE_DIRECTORY__ + "/data/lynx-hare.csv")
- [ (code "hare").Value, Time.TimeSeries.fromObservations (csv.Rows |> Seq.map(fun r -> float r.Hare, r.Year))
- (code "lynx").Value, Time.TimeSeries.fromObservations (csv.Rows |> Seq.map(fun r -> float r.Lynx, r.Year)) ] |> Map.ofList
+let data =
+ let csv = PopulationData.Load(__SOURCE_DIRECTORY__ + "/data/lynx-hare.csv")
+
+ [ (code "hare").Value, Time.TimeSeries.fromObservations (csv.Rows |> Seq.map (fun r -> float r.Hare, r.Year))
+ (code "lynx").Value, Time.TimeSeries.fromObservations (csv.Rows |> Seq.map (fun r -> float r.Lynx, r.Year)) ]
+ |> Map.ofList
(*** include-value: data ***)
@@ -233,9 +236,7 @@ the main fitting function of the Bristlecone library.
let endCondition = Optimisation.EndConditions.afterIteration 10000
-let result =
- ``predator-prey``
- |> Bristlecone.fit engine endCondition data
+let result = ``predator-prey`` |> Bristlecone.tryFit engine endCondition data
(*** include-value: result ***)
@@ -267,4 +268,4 @@ the optimisation routine:
*)
Graphing.parameterTrace result
-(*** include-it-raw ***)
\ No newline at end of file
+(*** include-it-raw ***)
diff --git a/docs/examples/shrub-resource.fsx b/docs/examples/shrub-resource.fsx
index 4e48196..08be732 100644
--- a/docs/examples/shrub-resource.fsx
+++ b/docs/examples/shrub-resource.fsx
@@ -3,13 +3,14 @@
title: Shrub response to a single limiting resource
category: Examples
categoryindex: 3
-index: 1
+index: 3
---
*)
(*** condition: prepare ***)
#nowarn "211"
#r "../../src/Bristlecone/bin/Debug/netstandard2.0/Bristlecone.dll"
+#r "../../src/Bristlecone.Dendro/bin/Debug/netstandard2.0/Bristlecone.Dendro.dll"
#r "nuget: MathNet.Numerics.FSharp,5.0.0"
(*** condition: fsx ***)
#if FSX
@@ -20,461 +21,722 @@ index: 1
#r "nuget: Bristlecone,{{fsdocs-package-version}}"
#endif // IPYNB
-////////////////////////////////////////////////////
-// Plant nitrogen limitation using wood rings
-// and nitrogen isotopes
-////////////////////////////////////////////////////
+(*** hide ***)
+module Allometric =
+
+ open Bristlecone
+
+ let pi = System.Math.PI
+
+ module Constants =
+
+ // Empirically-derived parameters:
+ let k5 = 19.98239 // Allometric fit to Yamal shrub BD-length data #1 (in centimetres)
+ let k6 = 0.42092 // Allometric fit to Yamal shrub BD-length data #2 (in centimetres)
+
+ // Constants from the literature:
+ let a = 2. // the number of child branches added to previous branches (including the tops of the stems) for a shrub
+ let p = 0.5 // the length of a child branch as a proportion of its parent branch/stem
+ let lmin = 20. //cm. the length at which a stem or branch gets child branches
+ let rtip = 0.1 //cm. the radius of the outermost tip of a stem or branch
+ let b = 0.0075 // the ratio of the basal radius of a stem or branch and its length
+ let salixWoodDensity = 0.5 // g / cm3 (from internet)
+ let numberOfStems = 2.2
+
+ module NiklasAndSpatz_Allometry =
+
+ let nthroot n A =
+ let rec f x =
+ let m = n - 1.
+ let x' = (m * x + A / x ** m) / n
+
+ match abs (x' - x) with
+ | t when t < abs (x * 1e-9) -> x'
+ | _ -> f x'
+
+ f (A / double n)
+
+ /// Gives the basal radius in centimetres of a stem/branch given its length in centimetres. Function from Niklas and Spatz (2004).
+ /// The basal radius is always positive.
+ let basalRadius k5 k6 stemLength =
+ max (100. * ((0.01 * stemLength + k6) / k5) ** (3. / 2.) / 2.) 1e-06
+
+ /// Inverse equation of basalRadius.
+ let stemLength k5 k6 radius =
+ max (2. * ((nthroot 3. 2.) * 5. ** (2. / 3.) * k5 * radius ** (2. / 3.) - 50. * k6)) 1e-06
+
+ module Götmark2016_ShrubModel =
+
+ /// Total shrub volume given height and number of stems
+ let shrubVolume b a rtip p lmin k5 k6 n h =
+
+ let radius = NiklasAndSpatz_Allometry.basalRadius k5 k6
+
+ let mainStemVolume =
+ match radius h with
+ | r when r > rtip -> n * pi * h * ((radius h) ** 2. + (radius h) * rtip + rtip ** 2.) / 3.
+ | _ -> n * pi * h * rtip ** 2.
+
+ let mutable volume = mainStemVolume
+ let mutable k = 0.
+
+ while (p ** k * h > lmin * 2. / 3.) do
+ let volToAdd =
+ match (p ** k * h < lmin) with
+ | true ->
+ match (b * 3. * p * (p ** k * h - 2. * lmin / 3.) > rtip) with
+ | true ->
+ n
+ * a
+ * (a + 1.) ** (float k)
+ * pi
+ * 3.
+ * p
+ * (p ** k * h - 2. * lmin / 3.)
+ * ((radius (3. * p * (p ** k * h - 2. * lmin / 3.)))
+ * (radius (3. * p * (p ** k * h - 2. * lmin / 3.)) * rtip + rtip ** 2.))
+ / 3.
+ | false ->
+ n
+ * a
+ * (a + 1.) ** (float k)
+ * 3.
+ * p
+ * (p ** k * h - 2. * lmin / 3.)
+ * pi
+ * rtip ** 2.
+ | false ->
+ match (radius (p ** (k + 1.) * h) > rtip) with
+ | true ->
+ n
+ * a
+ * (a + 1.) ** (float k)
+ * pi
+ * p ** (k + 1.)
+ * h
+ * ((radius (p ** (k + 1.) * h)) ** 2.
+ + (radius (p ** (k + 1.) * h)) * rtip
+ + rtip ** 2.)
+ / 3.
+ | false -> n * a * (a + 1.) ** (float k) * p ** (k + 1.) * h * pi * rtip ** 2.
+
+ volume <- volume + volToAdd
+ k <- k + 1.
+
+ k, volume
+
+
+ module Allometrics =
+
+ open Götmark2016_ShrubModel
+ open Bristlecone.Dendro
+
+ let private removeUnit (x: float<_>) = float x
+
+ let mass woodDensity volume = volume * woodDensity
+
+ let massToVolume woodDensity mass = mass / woodDensity
+
+ let shrubBiomass b a rtip p lmin k5 k6 n woodDensity (radius: float) =
+ radius
+ |> removeUnit
+ |> NiklasAndSpatz_Allometry.stemLength k5 k6
+ |> shrubVolume b a rtip p lmin k5 k6 n
+ |> snd
+ |> mass woodDensity
+
+ let shrubRadius b a rtip p lmin k5 k6 n woodDensity mass =
+ let findRadius volume =
+ let v x =
+ x
+ |> NiklasAndSpatz_Allometry.stemLength k5 k6
+ |> shrubVolume b a rtip p lmin k5 k6 n
+ |> snd
+
+ let f = (fun x -> (v x) - volume)
+ Statistics.RootFinding.bisect 0 200 f 0.01 100.00 1e-8 // Assumption that shrub radius is between 0.01 and 100.0cm.
+
+ mass |> massToVolume woodDensity |> findRadius
+
+ let shrubHeight k5 k6 radius =
+ radius |> NiklasAndSpatz_Allometry.stemLength k5 k6
+
+ module Proxies =
+
+ open Bristlecone.Dendro
+
+ /// Radius in millimetres
+ let toBiomassMM (radius: float) =
+ radius / 10.
+ |> Allometrics.shrubBiomass
+ Constants.b
+ Constants.a
+ Constants.rtip
+ Constants.p
+ Constants.lmin
+ Constants.k5
+ Constants.k6
+ Constants.numberOfStems
+ Constants.salixWoodDensity
+
+ /// Biomass in grams.
+ let toRadiusMM biomassGrams =
+ if
+ System.Double.IsNaN biomassGrams
+ || System.Double.IsInfinity biomassGrams
+ || System.Double.IsNegativeInfinity biomassGrams
+ then
+ nan
+ else
+ let radiusCm =
+ biomassGrams
+ |> Allometrics.shrubRadius
+ Constants.b
+ Constants.a
+ Constants.rtip
+ Constants.p
+ Constants.lmin
+ Constants.k5
+ Constants.k6
+ Constants.numberOfStems
+ Constants.salixWoodDensity
+
+ radiusCm * 10.
-(* An example Bristlecone script for working with
- wood ring datasets. *)
-open Bristlecone // Opens Bristlecone core library and estimation engine
-open Bristlecone.Language // Open the language for writing Bristlecone models
+(**
+## A fully-worked example: soil nutrients and plant growth
+
+The following code demonstrates the key functionality of
+Bristlecone, using a real-world example. Specifically, it
+demonstrates model-fitting and model-selection for the
+mechanisms controlling nutrient limitation to shrub individuals
+in the Arctic tundra.
+
+The hypotheses are tested on a per-shrub basis. For each shrub,
+there are four components that may vary in the model, which
+when multiplied together results in 12 possible combinations
+of mechanisms for each shrub.
+
+First, we must load Bristlecone.
+*)
+
+open Bristlecone // Opens Bristlecone core library and estimation engine
+open Bristlecone.Language // Open the language for writing Bristlecone models
open Bristlecone.Time
-// // 1. Define basic model system
-// // ----------------------------
-// // First, we define a 'base model' into which we can insert
-// // components that represent different hypotheses.
-
-// let baseModel =
-
-// /// Transform δ15N to N availability.
-// let ``δ15N -> N availability`` =
-// (Constant 100. * Environment "N" + Constant 309.) / Constant 359.
-
-// /// Plant uptake of N from soil, which may be turned on or off
-// let uptake f geom =
-// (geom This) * This * (f ``δ15N -> N availability``)
-
-// /// 1. Cumulative stem biomass
-// let ``db/dt`` geom nLimitation =
-// Parameter "r" * (uptake nLimitation geom) - Parameter "γ[b]" * This
-
-// /// 2. Soil nitrogen availability
-// let ``dN/dt`` geom feedback limitationName nLimitation =
-// if limitationName = "None"
-// then Parameter "λ" - Parameter "γ[N]" * ``δ15N -> N availability`` + feedback This
-// else
-// Parameter "λ" - Parameter "γ[N]" * ``δ15N -> N availability``
-// + feedback This
-// - (geom This) * This * (nLimitation ``δ15N -> N availability``)
-
-// /// 3. Stem radius (a 'Measurement' variable)
-// let stemRadius lastRadius lastEnv env =
-// let oldCumulativeMass = lastEnv |> lookup "bs"
-// let newCumulativeMass = env |> lookup "bs"
-// if (newCumulativeMass - oldCumulativeMass) > 0.
-// then newCumulativeMass |> Allometric.Proxies.toRadiusMM
-// else lastRadius
-
-// fun geom feedback (nLimitMode,nLimitation) ->
-// Model.empty
-// |> Model.addEquation "bs" (``db/dt`` geom nLimitation)
-// |> Model.addEquation "N" (``dN/dt`` geom feedback nLimitMode nLimitation)
-// |> Model.includeMeasure "x" stemRadius
-// |> Model.estimateParameter "λ" notNegative 0.001 0.500
-// |> Model.estimateParameter "γ[N]" notNegative 0.001 0.200
-// |> Model.estimateParameter "γ[b]" notNegative 0.001 0.200
-// |> Model.useLikelihoodFunction (ModelLibrary.Likelihood.bivariateGaussian "x" "N")
-// |> Model.estimateParameter "ρ" noConstraints -0.500 0.500
-// |> Model.estimateParameter "σ[x]" notNegative 0.001 0.100
-// |> Model.estimateParameter "σ[y]" notNegative 0.001 0.100
-
-
-// // 2. Define competing hypotheses
-// // ----------------------------
-// /// Define 12 alternative hypotheses by defining three interchangeable components:
-// /// - Asymptotic plant size (2 types);
-// /// - Plant-soil feedback presence / absence (2 types); and
-// /// - Nitrogen limitation form (3 types).
-// /// The product of each of the three components with the base model forms 12
-// /// alternative hypotheses, each represented as a `ModelSystem`.
-// let hypotheses =
-
-// // 1. Setup two alternatives for geometric mode.
-// let chapmanRichards mass = Constant 1. - (mass / (Parameter "k" * Constant 1000.))
-// let geometricModes = modelComponent "Geometric constraint" [
-// subComponent "None" (Constant 1. |> (*))
-// subComponent "Chapman-Richards" chapmanRichards
-// |> estimateParameter "k" notNegative 3.00 5.00 // Asymptotic biomass (in kilograms)
-// ]
-
-// // 2. Setup two alternatives for plant-soil feedbacks.
-// let biomassLoss biomass = (Parameter "ɑ" / Constant 100.) * biomass * Parameter "γ[b]"
-// let feedbackModes = modelComponent "Plant-Soil Feedback" [
-// subComponent "None" (Constant 1. |> (*))
-// subComponent "Biomass Loss" biomassLoss
-// |> estimateParameter "ɑ" notNegative 0.01 1.00 // N-recycling efficiency
-// ]
-
-// // 3. Setup three alternatives for the form of plant N limitation.
-
-// let saturating minimumNutrient nutrient =
-// let hollingModel n = (Parameter "a" * n) / (Constant 1. + (Parameter "a" * Parameter "b" * Parameter "h" * n))
-// Conditional(fun compute ->
-// if compute (hollingModel minimumNutrient) < 1e-12
-// then Invalid
-// else hollingModel nutrient )
-
-// let linear min resource =
-// Conditional(fun compute ->
-// if compute (Parameter "a" * min) < 1e-12 then Invalid else Parameter "a" * resource)
-
-// let limitationModes = modelComponent "N-limitation" [
-// subComponent "Saturating" (saturating (Constant 5.))
-// |> estimateParameter "a" notNegative 0.100 0.400
-// |> estimateParameter "h" notNegative 0.100 0.400
-// |> estimateParameter "r" notNegative 0.500 1.000
-// subComponent "Linear" (linear (Constant 5.))
-// |> estimateParameter "a" notNegative 0.100 0.400
-// |> estimateParameter "r" notNegative 0.500 1.000
-// subComponent "None" (Constant 1. |> (*))
-// |> estimateParameter "r" notNegative 0.500 1.000
-// ]
-
-// baseModel
-// |> Hypotheses.createFromComponent geometricModes
-// |> Hypotheses.useAnother feedbackModes
-// |> Hypotheses.useAnotherWithName limitationModes
-// |> Hypotheses.compile
-
-
-// // 3. Setup Bristlecone Engine
-// // ----------------------------
-// // A bristlecone engine provides a fixed setup for estimating parameters from data.
-// // Use the same engine for all model fits within a single study.
-
-// let engine =
-// Bristlecone.mkContinuous
-// |> Bristlecone.withContinuousTime Integration.MathNet.integrate
-// |> Bristlecone.withConditioning Conditioning.RepeatFirstDataPoint
-// |> Bristlecone.withTunedMCMC [ Optimisation.MonteCarlo.TuneMethod.CovarianceWithScale 0.200, 250, Optimisation.EndConditions.afterIteration 20000 ]
-
-
-// // 4. Test Engine and Model
-// // ----------------------------
-// // Running a full test is strongly recommended. The test will demonstrate if the current
-// // configuration can find known parameters for a model. If this step fails, there is an
-// // issue with either your model, or the Bristlecone configuration.
-
-// let testSettings =
-// Test.create
-// |> Test.addNoise (Test.Noise.tryAddNormal "sigma[y]" "N")
-// |> Test.addNoise (Test.Noise.tryAddNormal "sigma[x]" "bs")
-// |> Test.addGenerationRules [
-// Test.GenerationRules.alwaysMoreThan -3. "N"
-// Test.GenerationRules.alwaysLessThan 20. "N"
-// Test.GenerationRules.alwaysMoreThan 0. "bs"
-// Test.GenerationRules.monotonicallyIncreasing "x" ] // There must be at least 10mm of wood production
-// |> Test.addStartValues [
-// "x", 5.0
-// "bs", 5.0 |> Allometric.Proxies.toBiomassMM
-// "N", 3.64 ]
-// |> Test.withTimeSeriesLength 30
-// |> Test.endWhen (Optimisation.EndConditions.afterIteration 1000)
-
-// let testResult =
-// hypotheses
-// |> List.map(fst >> Bristlecone.testModel engine testSettings)
-
-
-// // 5. Load Real Data
-// // ----------------------------
-// // Here, we are using the Bristlecone.Dendro package to
-// // read in dendroecological data.
-
-// open Bristlecone.Dendro
-
-// let shrubData =
-// let plants = Data.PlantIndividual.loadRingWidths (__SOURCE_DIRECTORY__ + "/../data/yamal-rw.csv")
-// let isotopeData = Data.PlantIndividual.loadLocalEnvironmentVariable (__SOURCE_DIRECTORY__ + "/../data/yuribei-d15N-imputed.csv")
-// plants |> PlantIndividual.zipEnvMany "N" isotopeData
-
-
-// // 6. Fit Models to Real Data
-// // -----------------------------------
-
-
-
-// // 7. Calculate model comparison statistics
-// // -----------------------------------
-
-
-// // 3. Load Real Data and Estimate
-// // ----------------------------
-
-// // Define the start values for a model system, as follows:
-// // initial radius =
-// let startValues (startDate:System.DateTime) (plant:PlantIndividual.PlantIndividual) =
-// let initialRadius =
-// match plant.Growth with
-// | PlantIndividual.PlantGrowth.RingWidth s ->
-// match s with
-// | GrowthSeries.Absolute c -> c.Head |> fst
-// | _ -> invalidOp "Not applicable"
-// | _ -> invalidOp "Not applicable"
-// let initialMass = initialRadius |> ModelComponents.Proxies.toBiomassMM
-// let initialNitrogen = plant.Environment.[code "N"].Head |> fst
-// [ ("x", initialRadius)
-// ("N", initialNitrogen)
-// ("bs", initialMass) ] |> Map.ofList
-
-// let workPackages shrubs hypotheses engine saveDirectory =
-// seq {
-// for s in shrubs do
-
-// // 1. Arrange the subject and settings
-// let shrub = s |> PlantIndividual.toCumulativeGrowth
-// let common = shrub |> PlantIndividual.keepCommonYears
-// let startDate = (common.Environment.["N"]).StartDate |> snd
-// let startConditions = getStartValues startDate shrub
-// let e = engine |> Bristlecone.withConditioning (Custom startConditions)
-
-// // 2. Setup batches of dependent analyses
-// for h in [ 1 .. hypotheses |> List.length ] do
-// for _ in [ 1 .. Options.chains ] do
-// if h = 1 || h = 2 || h = 7 || h = 8 || h = 9 || h = 10 then yield async {
-// // A. Compute result
-// let result = Bristlecone.PlantIndividual.fit e Options.endWhen hypotheses.[h-1] common |> fst
-// // B. Save to file
-// Bristlecone.Data.EstimationResult.saveAll saveDirectory s.Identifier.Value h 1 result
-// return result }
-// }
-
-// // Orchestrate the analyses
-// let work = workPackages shrubs hypotheses Options.engine Options.resultsDirectory
-// let run() = work |> Seq.iter (OrchestrationMessage.StartWorkPackage >> Options.orchestrator.Post)
-
-
-// let saveDiagnostics () =
-
-// // 1. Get all results sliced by plant and hypothesis
-// let results =
-// let get subject model modelId = Bristlecone.Data.EstimationResult.loadAll Options.resultsDirectory subject.Identifier.Value model modelId
-// Bristlecone.ModelSelection.ResultSet.arrangeResultSets shrubs hypotheses get
-
-// // 2. Save convergence statistics to file
-// // NB include only results within 5 likelihood of the minimum (to remove outliers)
-// results
-// // |> Seq.map(fun (x,a,b,c) -> x.Identifier.Value,a,b,c)
-// // |> Seq.toList
-// |> Diagnostics.Convergence.gelmanRubinAll 10000
-// |> Data.Convergence.save Options.resultsDirectory
-
-// // 3. Save Akaike weights to file
-// results
-// |> ModelSelection.Select.weights
-// |> Seq.map(fun (x,a,b,c) -> x.Identifier.Value,a,b,c)
-// |> Data.ModelSelection.save Options.resultsDirectory
-
-// // 4. Save out logged components
-// // results
-// // |> Seq.map(fun r -> calculateComponents fit Options.engine r)
-
-
-
-
-
-// // What about one-step-ahead predictions?
-
-// open Bristlecone.Diagnostics.ModelComponents
-// open FSharp.Data
-
-// // 1. Load in all MLE results and pick the best for each shrub x hypothesis.
-// // - Must convert from old bristlecone format to new one.
-
-// module LoadOldTrace =
-
-// open Bristlecone.Data.Trace
-// open Bristlecone.Data.Config
-
-// type OldTrace = CsvProvider<"/Users/andrewmartin/Desktop/Bristlecone Results/Paper1-Yuribei-Annual/dphil-shrub-output-YUSL03A-1-92aadb09-f035-4373-aacd-a16ed7dec822.csv">
-
-// let fileMatch directory subject modelId =
-// let path = System.IO.DirectoryInfo(directory)
-// if path.Exists then
-// let files = path.GetFiles(sprintf "dphil-shrub-output-%s-%i-*.csv" subject modelId)
-// let regex = sprintf "dphil-shrub-output-%s-%i-(%s).csv" subject modelId regexGuid
-// files |> Seq.choose(fun f ->
-// let m = System.Text.RegularExpressions.Regex.Match(f.Name, regex)
-// if m.Success
-// then
-// let guid = m.Groups.[1].Value |> System.Guid.Parse
-// (guid, f) |> Some
-// else None )
-// else invalidArg "directory" "The specified directory does not exist"
-
-// let toTrace (data:OldTrace) : (float * float []) list =
-// data.Rows
-// |> Seq.groupBy(fun r -> r.Iteration)
-// |> Seq.map(fun (i,r) ->
-// i,
-// (r |> Seq.head).NegativeLogLikelihood,
-// r |> Seq.map(fun r -> r.ParameterValue) |> Seq.toArray)
-// |> Seq.sortByDescending(fun (i,_,_) -> i)
-// |> Seq.map(fun (_,x,y) -> x,y)
-// |> Seq.toList
-
-// let loadOldTrace directory subject modelId =
-// let traceFiles = fileMatch directory subject modelId
-// traceFiles
-// |> Seq.choose(fun (i,f) ->
-// printfn "File loading is %s" f.FullName
-// let data = OldTrace.Load f.FullName
-// match data.Rows |> Seq.length with
-// | 0 -> None
-// | _ -> data |> toTrace |> Some )
-
-// let oldShrubDir = "/Users/andrewmartin/Documents/DPhil Zoology/Bristlecone Results/YuribeiAnnual/"
-
-// let oldMleResults =
-// seq {
-// for s in shrubs do
-// for hi in [ 1 .. 12 ] do
-// let oldResult = LoadOldTrace.loadOldTrace oldShrubDir s.Identifier.Value hi |> Seq.concat |> Seq.toList
-// let mle = oldResult |> Seq.minBy(fun (p,l) -> p)
-// yield s, hi, mle
-// }
-// |> Seq.toList
-
-
-// // 2. Setup each so that
-
-// module Cool =
-
-// open Bristlecone.Bristlecone
-// open Bristlecone.Optimisation
-
-// /// "How good am I at predicting the next data point"?
-// ///
-// let oneStepAhead engine hypothesis (preTransform:CodedMap>->CodedMap>) (timeSeries) (estimatedTheta:ParameterPool) =
-// let mleToBounds mlePool = mlePool |> Map.map(fun k v -> Parameter.create (Parameter.detatchConstraint v |> snd) (v |> Parameter.getEstimate) (v |> Parameter.getEstimate))
-// let hypothesisMle : ModelSystem = { hypothesis with Parameters = mleToBounds estimatedTheta }
-// let pairedDataFrames =
-// timeSeries
-// |> Map.map(fun _ fitSeries ->
-// fitSeries
-// |> TimeSeries.toObservations
-// |> Seq.pairwise
-// |> Seq.map (fun (t1,t2) -> TimeSeries.fromObservations [t1; t2] |> TimeSeries.map(fun (x,y) -> x )))
-// printfn "Paired data frames = %A" pairedDataFrames
-// let timeParcelCount = (pairedDataFrames |> Seq.head).Value |> Seq.length
-// printfn "Time parcels: %i" timeParcelCount
-// let data =
-// seq { 1 .. timeParcelCount }
-// |> Seq.map(fun i -> pairedDataFrames |> Map.map(fun _ v -> v |> Seq.item (i-1)) |> preTransform)
-// printfn "Data: %A" data
-
-// // TODO Remove this hack:
-// // It is predicting with a repeated first point... so...
-// // The next point estimate is at t1
-// // The next point observation is at t2
-// data
-// |> Seq.map (fun d ->
-// let est = Bristlecone.fit (engine |> withCustomOptimisation Optimisation.None.passThrough |> withConditioning RepeatFirstDataPoint) (EndConditions.afterIteration 0) d hypothesisMle
-// printfn "Estimated time series is %A" (est |> fst).Series
-
-// let nextObservation = d |> Map.map(fun c ts -> ts |> TimeSeries.toObservations |> Seq.skip 1 |> Seq.head)
-// let paired =
-// nextObservation
-// |> Map.map(fun code obs ->
-// let nextEstimate = ((est |> fst).Series.[code].Values |> Seq.head).Fit
-// obs |> snd, obs |> fst, nextEstimate
-// )
-
-// paired
-// )
-// |> Seq.toList
-
-// /// Perform n-step-ahead computation on the hypothesis and plant.
-// let predictAhead (engine:EstimationEngine.EstimationEngine) system (plant:PlantIndividual) preTransform estimate =
-// let g =
-// match plant.Growth |> growthSeries with
-// | Absolute g -> g
-// | Cumulative g -> g
-// | Relative g -> g
-// let predictors = plant.Environment |> Map.add (ShortCode.create "x") g
-// oneStepAhead engine system preTransform predictors estimate
-
-
-// let oneStepPredictions =
-// oldMleResults
-// |> List.map(fun (s,hi,mle) ->
-
-// // 0. Convert x into biomass
-// let preTransform (data:CodedMap>) =
-// data
-// |> Map.toList
-// |> List.collect(fun (k,v) ->
-// if k = code "x"
-// then [ (k, v); (code "bs", v |> TimeSeries.map(fun (x,_) -> x |> ModelComponents.Proxies.toBiomassMM)) ]
-// else [ (k, v)] )
-// |> Map.ofList
-
-// // 1. Arrange the subject and settings
-// let shrub = s |> PlantIndividual.toCumulativeGrowth
-// let common = shrub |> PlantIndividual.keepCommonYears
-// let startDate = (common.Environment.[code "N"]).StartDate |> snd
-// let startConditions = getStartValues startDate shrub
-// let e = Options.engine |> Bristlecone.withConditioning (Custom startConditions)
-
-// let mlePool =
-// hypotheses.[hi-1].Parameters
-// |> Map.toList
-// |> List.mapi(fun i (sc,p) ->
-// let est = (mle |> snd).[i]
-// sc, Parameter.setEstimate p est)
-// |> ParameterPool
-
-// Cool.predictAhead e hypotheses.[hi-1] common preTransform mlePool
-// |> List.map(fun r -> s.Identifier.Value, hi, r)
-// )
-
-
-// type SOSRow = {
-// Year: int
-// Variable: string
-// Shrub: string
-// Hypothesis: int
-// Observation: float
-// OneStepPrediction: float
-// }
-
-// let predictions =
-// oneStepPredictions
-// |> List.collect(fun m ->
-// m |> Seq.collect(fun (s,h,m) ->
-// m |> Seq.map(fun kv ->
-// let y,a,b = kv.Value
-// let x =
-// {
-// Year = y.Year
-// Variable = kv.Key.Value
-// Observation = a
-// Shrub = s
-// Hypothesis = h
-// OneStepPrediction = b }
-// x
-// )) |> Seq.toList )
-
-// predictions
-// |> List.map(fun x -> sprintf "%s,%i,%i,%s,%f,%f" x.Shrub x.Hypothesis x.Year x.Variable x.Observation x.OneStepPrediction)
-// |> List.append ["Shrub,Hypothesis,Year,Variable,Observed,OneStepPrediction"]
-// |> String.concat "\n"
-// |> fun x -> System.IO.File.WriteAllText("/Users/andrewmartin/Desktop/onestepahead.csv", x)
-
-// // Root mean squared error is the average squared differences, then square rooted.
-// let rmse =
-// predictions
-// |> List.groupBy(fun x -> x.Shrub, x.Hypothesis, x.Variable)
-// |> List.map(fun ((s,h,v),x) ->
-// let sos = x |> Seq.averageBy(fun x -> (x.Observation - x.OneStepPrediction) ** 2.)
-// s, h, v, sqrt sos)
-
-// rmse
-// |> List.map(fun (s,h,v,sos) -> sprintf "%s,%i,%s,%f" s h v sos)
-// |> List.append ["Shrub,Hypothesis,Variable,RMSE"]
-// |> String.concat "\n"
-// |> fun x -> System.IO.File.WriteAllText("/Users/andrewmartin/Desktop/onestepahead-2.csv", x)
+(**
+### Defining a 'base model'
+
+Then, we define a `base model` system, into which we can nest
+the components that vary (which represent different hypotheses).
+
+Our base model consists of the following parts:
+
+* **An empirical transform.** As this model uses long-term ecological data,
+ there is a function that translates from an ecological proxy (stable nitrogen isotope)
+ into the desired variable - soil nitrogen availability.
+* **Two ordinary differential equations (ODEs).** These represent the soil N availability
+ and plant biomass through time.
+* **A *measurement* variable.** Bristlecone 'measurements' are effectively
+ derived time-series that are calulated from the integrated ODE output. Here,
+ we use a measurement variable to transform biomass into wood ring production,
+ where rings are only produced where plant biomass has increased during the year.
+
+The code for the parts of the base model is shown below.
+*)
+
+// Empirical transform from δ15N to N availability.
+let ``δ15N -> N availability`` =
+ (Constant 100. * Environment "N" + Constant 309.) / Constant 359.
+
+// ODE1. Cumulative stem biomass
+let ``db/dt`` geom nLimitation =
+ Parameter "r" * (geom This) * This * (nLimitation ``δ15N -> N availability``)
+ - Parameter "γ[b]" * This
+
+// ODE2. Soil nitrogen availability
+let ``dN/dt`` geom feedback limitationName nLimitation =
+ if limitationName = "None" then
+ Parameter "λ" - Parameter "γ[N]" * ``δ15N -> N availability`` + feedback This
+ else
+ Parameter "λ" - Parameter "γ[N]" * ``δ15N -> N availability`` + feedback This
+ - (geom This) * This * (nLimitation ``δ15N -> N availability``)
+
+// Measurement variable: stem radius
+let stemRadius lastRadius lastEnv env =
+ let oldCumulativeMass = lastEnv |> lookup "bs"
+ let newCumulativeMass = env |> lookup "bs"
+
+ if (newCumulativeMass - oldCumulativeMass) > 0. then
+ newCumulativeMass |> Allometric.Proxies.toRadiusMM // NB Allometrics is given in full in the downloadable script.
+ else
+ lastRadius
+
+(**
+Once we have defined the components, we can scaffold them into a model system.
+We can plug in the nestable hypotheses (defined further below) by defining
+the base model as a function that takes parameters representing the
+alternative hypotheses.
+*)
+
+let ``base model`` geometricConstraint plantSoilFeedback (nLimitMode, nLimitation) =
+ Model.empty
+ |> Model.addEquation "bs" (``db/dt`` geometricConstraint nLimitation)
+ |> Model.addEquation "N" (``dN/dt`` geometricConstraint plantSoilFeedback nLimitMode nLimitation)
+ |> Model.includeMeasure "x" stemRadius
+ |> Model.estimateParameter "λ" notNegative 0.001 0.500
+ |> Model.estimateParameter "γ[N]" notNegative 0.001 0.200
+ |> Model.estimateParameter "γ[b]" notNegative 0.001 0.200
+ |> Model.useLikelihoodFunction (ModelLibrary.Likelihood.bivariateGaussian "x" "N")
+ |> Model.estimateParameter "ρ" noConstraints -0.500 0.500
+ |> Model.estimateParameter "σ[x]" notNegative 0.001 0.100
+ |> Model.estimateParameter "σ[y]" notNegative 0.001 0.100
+
+(**
+### Defining the competing hypotheses
+
+Here, we define 12 alternative hypotheses by defining three interchangeable
+components:
+
+ - Asymptotic plant size (2 types);
+ - Plant-soil feedback presence / absence (2 types); and
+ - Nitrogen limitation form (3 types).
+
+Once we have defined each of the three components, we can take the
+product of them with the base model which forms 12
+alternative hypotheses, each represented as a `ModelSystem`.
+
+#### Geometric constraint
+
+Plants do not grow indefinitely, but eventually reach an asymptotic mass
+owing either to geometric or resource constraints. Here, we define two
+competing hypotheses: that a shrub does not show evidence of nearing its
+asymptote, or that it does (based on a Chapman-Richards growth function).
+
+In Bristlecone, we use `modelComponent` to construct a pluggable component
+into a model system. We pass `modelComponent` a list of `subComponent`s, which
+each have a name and an equation. In its equation, a model component can
+take one or many parameters, but these must match the signiture required
+by the hole in the base model. For example, the geometric constraint here
+takes the current `mass` only.
+
+In addition, a `subComponent` may require additional estimatable parameters
+over the base model. In this case, the Chapman-Richards model requires an
+extra parameter *K*, which represents the asymptotic biomass. These may be
+added to a `subComponent` by using `|> estimateParameter` afterwards, as below.
+*)
+
+let ``geometric constraint`` =
+ modelComponent
+ "Geometric constraint"
+ [ subComponent "None" (Constant 1. |> (*))
+ subComponent "Chapman-Richards" (fun mass -> Constant 1. - (mass / (Parameter "k" * Constant 1000.)))
+ |> estimateParameter "k" notNegative 3.00 5.00 ] // Asymptotic biomass (in kilograms)
+
+(**
+#### Plant-soil feedback
+
+The plant-soil feedback is the flow of nutrients from plant biomass into the
+soil available nitrogen pool. Here, we effectively turn on or off N input into
+the soil pool on biomass loss. In the base model, density-dependent biomass
+loss occurs. Turning on the feedback here creates an N input into soil based
+on the biomass lost multiplied by a conversion factor `ɑ`.
+*)
+
+let ``plant-soil feedback`` =
+
+ let biomassLoss biomass =
+ (Parameter "ɑ" / Constant 100.) * biomass * Parameter "γ[b]"
+
+ modelComponent
+ "Plant-Soil Feedback"
+ [ subComponent "None" (Constant 1. |> (*))
+ subComponent "Biomass Loss" biomassLoss
+ |> estimateParameter "ɑ" notNegative 0.01 1.00 ] // N-recycling efficiency
+
+(**
+#### Nitrogen limitation
+
+We specify three plausable mechanisms for nutrient limitation to shrub
+growth: (1) that growth is independent on soil N availability; (2) that
+growth is dependent on soil N availability in a linear way; or (3) that
+a mechanistic model of root-foraging (saturating) best represents
+N-limitation of shrub growth.
+
+A new concept here is the `Conditional` element in an equation. This
+term exposes a `ModelExpression -> float` (compute), allowing a calculation
+to be conditional on the state of parameters or values. In this example,
+we use it to restrict the models such that the N-limiting effect cannot
+be zero.
+*)
+
+let ``N-limitation to growth`` =
+
+ let saturating minimumNutrient nutrient =
+ let hollingModel n =
+ (Parameter "a" * n)
+ / (Constant 1. + (Parameter "a" (** Parameter "b"*) * Parameter "h" * n))
+
+ Conditional
+ <| fun compute ->
+ if compute (hollingModel minimumNutrient) < 1e-12 then
+ Invalid
+ else
+ hollingModel nutrient
+
+ let linear min resource =
+ Conditional
+ <| fun compute ->
+ if compute (Parameter "a" * min) < 1e-12 then
+ Invalid
+ else
+ Parameter "a" * resource
+
+ modelComponent
+ "N-limitation"
+ [ subComponent "Saturating" (saturating (Constant 5.))
+ |> estimateParameter "a" notNegative 0.100 0.400
+ |> estimateParameter "h" notNegative 0.100 0.400
+ |> estimateParameter "r" notNegative 0.500 1.000
+ subComponent "Linear" (linear (Constant 5.))
+ |> estimateParameter "a" notNegative 0.100 0.400
+ |> estimateParameter "r" notNegative 0.500 1.000
+ subComponent "None" (Constant 1. |> (*))
+ |> estimateParameter "r" notNegative 0.500 1.000 ]
+
+(**
+#### Putting the hypotheses together
+
+
+
+*)
+
+let hypotheses =
+ ``base model``
+ |> Hypotheses.createFromComponent ``geometric constraint``
+ |> Hypotheses.useAnother ``plant-soil feedback``
+ |> Hypotheses.useAnotherWithName ``N-limitation to growth``
+ |> Hypotheses.compile
+
+(**
+The resultant hypotheses (each of which is a `Hypotheses.Hypothesis`) are:
+*)
+
+(*** hide ***)
+printfn "There are %i hypotheses" hypotheses.Length
+(*** hide ***)
+hypotheses
+|> Seq.iteri (fun i h ->
+ printfn
+ "Hypothesis %i [%s]: %s"
+ (i + 1)
+ h.ReferenceCode
+ (h.Components |> Seq.map (fun x -> x.Implementation) |> String.concat " + "))
+(*** include-output ***)
+
+(**
+### Setting up a *Bristlecone engine*
+
+A bristlecone engine provides a fixed setup for estimating parameters from data.
+We use the same engine for all model fits within a single study.
+
+Here, we scaffold an engine from `Bristlecone.mkContinuous`, as we are working
+with continuous-time models.
+
+**Note.** In a non-docs environment, we may choose to use a logger for this
+example that outputs real-time traces in graphical format using `ggplot` in `R`.
+To do this, we can call `withOutput` with `Logging.RealTimeTrace.graphWithConsole 60. 10000`
+to use the graphics functions from the `Bristlecone.Charts.R` package.
+
+*)
+
+let output = Logging.Console.logger 1000
+
+let engine =
+ Bristlecone.mkContinuous
+ |> Bristlecone.withContinuousTime Integration.MathNet.integrate
+ |> Bristlecone.withOutput output
+ |> Bristlecone.withConditioning Conditioning.RepeatFirstDataPoint
+ |> Bristlecone.withCustomOptimisation (
+ Optimisation.MonteCarlo.SimulatedAnnealing.fastSimulatedAnnealing
+ 0.01
+ true
+ { Optimisation.MonteCarlo.SimulatedAnnealing.AnnealSettings.Default with
+ InitialTemperature = 100.
+ TemperatureCeiling = Some 100.
+ HeatRamp = (fun t -> t + 5.00)
+ BoilingAcceptanceRate = 0.85
+ HeatStepLength = Optimisation.EndConditions.afterIteration 1000
+ TuneLength = 1000
+ AnnealStepLength =
+ (fun x n ->
+ Optimisation.MonteCarlo.SimulatedAnnealing.EndConditions.improvementCount 5000 250 x n
+ || Optimisation.EndConditions.afterIteration 10000 x n) }
+ )
+
+(**
+### Testing the engine and model
+
+Running a full test is strongly recommended. The test will demonstrate if the current
+configuration can find known parameters for a model. If this step fails, there is an
+issue with either your model, or the Bristlecone configuration.
+*)
+
+let testSettings =
+ Test.create
+ |> Test.addNoise (Test.Noise.tryAddNormal "σ[y]" "N")
+ |> Test.addNoise (Test.Noise.tryAddNormal "σ[x]" "bs")
+ |> Test.addGenerationRules
+ [ Test.GenerationRules.alwaysMoreThan -3. "N"
+ Test.GenerationRules.alwaysLessThan 20. "N"
+ Test.GenerationRules.alwaysMoreThan 0. "bs"
+ Test.GenerationRules.monotonicallyIncreasing "x" ] // There must be at least 10mm of wood production
+ |> Test.addStartValues [ "x", 5.0; "bs", 5.0 |> Allometric.Proxies.toBiomassMM; "N", 3.64 ]
+ |> Test.withTimeSeriesLength 30
+ |> Test.endWhen (Optimisation.EndConditions.afterIteration 1000)
+
+(*** do-not-eval ***)
+let testResult =
+ hypotheses
+ |> List.map ((fun h -> h.Model) >> Bristlecone.testModel engine testSettings)
+
+
+(**
+### Load in the real dataset
+
+Here, we are using the Bristlecone.Dendro package to
+read in dendroecological data. For other problems, any
+method to wrangle the data into a `TimeSeries` is acceptable.
+
+*)
+
+open Bristlecone.Dendro
+
+let dataset =
+ let plants =
+ Data.PlantIndividual.loadRingWidths (__SOURCE_DIRECTORY__ + "/../data/yamal-rw.csv")
+
+ let isotopeData =
+ Data.PlantIndividual.loadLocalEnvironmentVariable (__SOURCE_DIRECTORY__ + "/../data/yuribei-d15N-imputed.csv")
+
+ plants |> PlantIndividual.zipEnvMany "N" isotopeData
+
+(**
+### Model-fitting to real data
+
+In this step, we will apply our `EstimationEngine` to the real
+shrub data and model hypotheses.
+
+Because there are so many hypotheses (x12) and individuals (x10),
+and we want to run replicate analyses (x3), it makes sense to
+run these 360 workloads in parallel. To do this, we can make use of
+the workflow tools in the `Bristlecone.Workflow` namespace.
+
+An orchestration agent is an agent that queues and runs work items
+in parallel. Below, we demonstrate setting up an agent:
+*)
+
+open Bristlecone.Workflow
+
+let orchestrator =
+ Orchestration.OrchestrationAgent(
+ writeOut = output,
+ maxSimultaneous = System.Environment.ProcessorCount,
+ retainResults = false
+ )
+
+(**
+
+Before we can run the models, there are some specific considerations
+required for this problem.
+
+#### Setting the start values
+
+Given time-series data, a common modelling problem is the question of
+what to set t = 0 as. One strategy is to repeat the first data point (t1)
+as t0. In this instance, out isotope time-series are shorter than the
+ring-width time-series, so we generally know the real value of one
+time-series but not the other. Because of this, we use a custom start
+point for each shrub individual.
+*)
+
+let startValues (startDate: System.DateTime) (plant: PlantIndividual.PlantIndividual) =
+ let removeUnit (x: float<_>) = float x
+
+ let initialRadius =
+ match plant.Growth with
+ | PlantIndividual.PlantGrowth.RingWidth s ->
+ match s with
+ | GrowthSeries.Cumulative c ->
+ let trimmed = c |> TimeSeries.trimStart (startDate - System.TimeSpan.FromDays(366.))
+
+ match trimmed with
+ | Some t -> t.Values |> Seq.head
+ | None -> failwith "Could not get t0 from ring-width series"
+ | _ -> invalidOp "Not applicable"
+ | _ -> invalidOp "Not applicable"
+
+ let initialMass = initialRadius |> Allometric.Proxies.toBiomassMM
+ let initialNitrogen = plant.Environment.[(code "N").Value].Head |> fst
+
+ [ ((code "x").Value, initialRadius |> removeUnit)
+ ((code "N").Value, initialNitrogen)
+ ((code "bs").Value, initialMass) ]
+ |> Map.ofList
+
+(**
+Next, we set up some configuration settings, then wrap up our hypotheses and shrubs
+as work packages for the orchestrator, where a work package is a `Async`.
+
+The function to make the work packages first sets up the common time period and custom
+start values for each shrub, then creates per-hypothesis work packages for that shrub.
+*)
+
+module Config =
+
+ let numberOfReplicates = 3
+ let resultsDirectory = "results/test/"
+ let thinTrace = Some 100
+ let endWhen = Optimisation.EndConditions.afterIteration 100000
+
+
+// Function to scaffold work packages
+let workPackages shrubs (hypotheses: Hypotheses.Hypothesis list) engine saveDirectory =
+ seq {
+ for s in shrubs do
+
+ // 1. Arrange the subject and settings
+ let shrub = s |> PlantIndividual.toCumulativeGrowth
+ let common = shrub |> PlantIndividual.keepCommonYears
+ let startDate = (common.Environment.[(code "N").Value]).StartDate |> snd
+ let startConditions = startValues startDate shrub
+ let e = engine |> Bristlecone.withConditioning (Conditioning.Custom startConditions)
+
+ // 2. Setup batches of dependent analyses
+ for h in hypotheses do
+ for _ in [ 1 .. Config.numberOfReplicates ] do
+ yield
+ async {
+ // A. Compute result
+ let result =
+ Bristlecone.tryFit e Config.endWhen (Bristlecone.fromDendro common) h.Model
+ // B. Save to file
+ match result with
+ | Ok r ->
+ Bristlecone.Data.EstimationResult.saveAll
+ saveDirectory
+ s.Identifier.Value
+ h.ReferenceCode
+ Config.thinTrace
+ r
+
+ return r
+ | Error e -> return failwithf "Error in package: %s" e
+ }
+ }
+
+let work = workPackages dataset hypotheses engine Config.resultsDirectory
+
+// Calling this function will run all of the analyses, but for this
+// example script we won't run them here.
+let run () =
+ work
+ |> Seq.iter (Orchestration.OrchestrationMessage.StartWorkPackage >> orchestrator.Post)
+
+(**
+### Model selection
+
+After calling `run ()`, we should have a folder containing csv files, three
+for each `EstimationResult`. We can load all of these in at once to
+calculate model comparison statistics.
+
+Functions for loading data are in the `Bristlecone.Data` namespace. You
+may notice that in the work package definition above we used functions from
+this namespace to save the results.
+*)
+let saveDiagnostics () =
+
+ // 1. Get all results sliced by plant and hypothesis
+ let results =
+ let get (subject: PlantIndividual.PlantIndividual) (hypothesis: Hypotheses.Hypothesis) =
+ Bristlecone.Data.EstimationResult.loadAll
+ Config.resultsDirectory
+ subject.Identifier.Value
+ hypothesis.Model
+ hypothesis.ReferenceCode
+
+ Bristlecone.ModelSelection.ResultSet.arrangeResultSets dataset hypotheses get
+
+ // 2. Save convergence statistics to file
+ results
+ |> Diagnostics.Convergence.gelmanRubinAll
+ 10000
+ (fun (s: PlantIndividual.PlantIndividual) -> s.Identifier.Value)
+ (fun (h: Hypotheses.Hypothesis) -> h.ReferenceCode)
+ |> Data.Convergence.save Config.resultsDirectory
+
+ // 3. Save Akaike weights to file
+ results
+ |> ModelSelection.Akaike.akaikeWeightsForSet (fun (h: Hypotheses.Hypothesis) -> h.ReferenceCode)
+ |> Seq.map (fun (x, a, b, c) -> x.Identifier.Value, a, b, c)
+ |> Data.ModelSelection.save Config.resultsDirectory
+
+
+ // // 4. Save out logged components
+ // results
+ // |> Seq.map(fun r ->
+ // Diagnostics.ModelComponents.calculateComponents fit engine r)
+
+ // 5. One-step ahead predictions
+
+ let bestFits =
+ Seq.allPairs dataset hypotheses
+ |> Seq.map (fun (s, h) ->
+ s, h, Bristlecone.Data.MLE.loadBest Config.resultsDirectory s.Identifier.Value h.Model h.ReferenceCode)
+
+ let oneStepPredictions =
+ bestFits
+ |> Seq.map (fun (s, h, mle) ->
+
+ // 0. Convert x into biomass
+ let preTransform (data: CodedMap>) =
+ data
+ |> Map.toList
+ |> List.collect (fun (k, v) ->
+ if k.Value = "x" then
+ [ (k, v)
+ ((code "bs").Value,
+ v |> TimeSeries.map (fun (x, _) -> x * 1. |> Allometric.Proxies.toBiomassMM)) ]
+ else
+ [ (k, v) ])
+ |> Map.ofList
+
+ // 1. Arrange the subject and settings (same as in model-fitting)
+ let shrub = s |> PlantIndividual.toCumulativeGrowth
+ let common = shrub |> PlantIndividual.keepCommonYears
+ let startDate = (common.Environment.[(code "N").Value]).StartDate |> snd
+ let startConditions = startValues startDate shrub
+
+ let e =
+ engine
+ |> Bristlecone.withConditioning (Bristlecone.Conditioning.Custom startConditions)
+
+ let result =
+ Bristlecone.oneStepAhead e h.Model preTransform (Bristlecone.fromDendro common) (mle |> snd |> snd)
+
+ // Save each n-step ahead result to a csv file
+ Bristlecone.Data.NStepAhead.save
+ Config.resultsDirectory
+ s.Identifier.Value
+ h.ReferenceCode
+ (mle |> fst)
+ 1
+ result
+
+ s.Identifier.Value, h.ReferenceCode, result)
+
+ Bristlecone.Data.NStepAhead.saveAllRMSE Config.resultsDirectory oneStepPredictions
+
+ ()
diff --git a/docs/examples/tutorial.fsx b/docs/examples/tutorial.fsx
deleted file mode 100644
index 0dbbe55..0000000
--- a/docs/examples/tutorial.fsx
+++ /dev/null
@@ -1,209 +0,0 @@
-(**
----
-title: Quick Tour
-category: Tutorials
-categoryindex: 5
-index: 1
----
-*)
-
-(*** condition: prepare ***)
-#nowarn "211"
-#r "../src/Bristlecone/bin/Debug/netstandard2.0/Bristlecone.dll"
-(*** condition: fsx ***)
-#if FSX
-#r "nuget: Bristlecone,{{fsdocs-package-version}}"
-#endif // FSX
-(*** condition: ipynb ***)
-#if IPYNB
-#r "nuget: Bristlecone,{{fsdocs-package-version}}"
-#endif // IPYNB
-
-// (**
-// Tutorial
-// ========================
-
-// * [Quick Start: Predator-Prey Dynamics](#quick-start-predator-prey-dynamics)
-// * [Defining the model](#defining-the-model)
-// * [The likelihood function](#The-likelihood-function)
-// * [Putting it all together](#Putting-it-all-together)
-// * [Does my model and method work?](#Does-my-model-and-method-work?)
-// * [Fitting to observation data](#Fitting-to-observation-data)
-// * [Time Modes (+ Integration)](#time-modes)
-// * [Optimisation](#optimisation)
-// * Working with Time Series
-// * Setting a Likelihood Function
-// * Measurement Variables
-// * Environmental Forcing
-
-// ## Quick Start: Predator-Prey Dynamics
-
-// Here we use the classic example of snowshoe hare and lynx predator-prey dynamics, to demonstrate the basic functions of Bristlecone. The full example is in ``/samples/1. lynx-hare.fsx``.
-
-// In an F# script or interactive session, first load and open the Bristlecone library.
-
-// *)
-
-// #load "bristlecone.fsx"
-// open Bristlecone
-// open Bristlecone.ModelSystem
-
-// (**### Defining the model
-
-// A model system is defined as a series of differential equations, parameters, and a likelihood function. The model system is here defined as a pair of ordinary differential equations.
-
-// *)
-// /// Number of snowshoe hares
-// let dhdt' hare lynx alpha beta =
-// alpha * hare - beta * hare * lynx
-
-// /// Number of lynx
-// let dldt' lynx hare delta gamma =
-// - gamma * lynx + delta * hare * lynx
-
-
-// (**The parameters required for the model are defined as a ``CodedMap``. For the predator-prey model, they are defined as:*)
-
-// [ // Natural growth rate of hares in absence of predation
-// ShortCode.create "alpha", Parameter.create Unconstrained 0.10 0.60
-// // Death rate per encounter of hares due to predation
-// ShortCode.create "beta", Parameter.create Unconstrained 0.001 0.0135
-// // Efficiency of turning predated hares into lynx
-// ShortCode.create "delta", Parameter.create Unconstrained 0.001 0.0135
-// // Natural death rate of lynx in the absence of food
-// ShortCode.create "gamma", Parameter.create Unconstrained 0.10 0.60
-// ] |> Map.ofList
-
-
-// (**### The likelihood function
-
-// For this example, we are simply using sum of squares as our measure of goodness of fit. Bristlecone includes sum of squares, and some neagtive log likelihood functions within the ``Bristlecone.Likelihood`` module.
-// *)
-
-// ModelLibrary.Likelihood.sumOfSquares ["hare"; "lynx"]
-
-// (**### Putting it all together
-
-// The above models, parameters, and likelihood function must be tied together in a ``ModelSystem`` type. A complete model system for this problem is defined here:
-// *)
-
-// let ``predator-prey`` =
-
-// let dhdt p _ x (e:Environment) =
-// dhdt' x
-// (e.[ShortCode.create "lynx"])
-// (p |> Pool.getEstimate "alpha")
-// (p |> Pool.getEstimate "beta")
-
-// let dldt p _ y (e:Environment) =
-// dldt' y
-// (e.[ShortCode.create "hare"])
-// (p |> Pool.getEstimate "delta")
-// (p |> Pool.getEstimate "gamma")
-
-// { Equations = [ code "hare", dhdt
-// code "lynx", dldt ] |> Map.ofList
-// Measures = [] |> Map.ofList
-// Parameters = [ code "alpha", parameter PositiveOnly 0.10 0.60
-// code "beta", parameter PositiveOnly 0.001 0.0135
-// code "delta", parameter PositiveOnly 0.001 0.0135
-// code "gamma", parameter PositiveOnly 0.10 0.60
-// ] |> Map.ofList
-// Likelihood = ModelLibrary.Likelihood.sumOfSquares ["hare"; "lynx"] }
-// (**
-
-// The intermediate functions ``dhdt`` and ``dldt`` connect the raw ODEs to the functional signature required by Bristlecone. This feeds the current parameter estimates and variable values into the model.
-
-// The modelling protocol is defined by an ``EstimationEngine``, which can be reused for any number of analyses.
-// *)
-
-// let engine =
-// Bristlecone.mkContinuous
-// |> Bristlecone.withConditioning RepeatFirstDataPoint
-
-
-// (**### Does my model and method work?
-
-// It is recommended that any new model and method combination is tested, to verify that - given a known set of parameters - the same parameters are estimated. Bristlecone can draw many sets of random parameters and test that these can be correctly estimated:
-// *)
-
-// let startValues =
-// [ code "lynx", 30.09
-// code "hare", 19.58 ] |> Map.ofList
-
-// //``predator-prey`` |> Bristlecone.testModel engine Options.testSeriesLength startValues Options.iterations Options.burn
-
-// (**
-
-// ### Fitting to observation data
-
-// In the example script, we load data using FSharp.Data, from a CSV file. You can choose to load data using any method. The raw data is parsed into a Map of TimeSeries using the ``TimeSeries.createVarying`` function.
-// *)
-
-// type PopulationData = FSharp.Data.CsvProvider<"../samples/data/lynx-hare.csv">
-
-// let data =
-// let csv = PopulationData.Load "../samples/data/lynx-hare.csv"
-// [ code "hare",
-// TimeSeries.fromObservations (csv.Rows |> Seq.map(fun r -> r.Year, float r.Hare))
-// code "lynx",
-// TimeSeries.fromObservations (csv.Rows |> Seq.map(fun r -> r.Year, float r.Lynx))
-// ] |> Map.ofList
-
-
-// (**To get model fitting results:*)
-
-// let result =
-// ``predator-prey``
-// //|> Bristlecone.fit engine Options.iterations Options.burn data
-
-
-// (**## Time Modes
-
-// Bristlecone can run models in either discrete time, or continuous time. When running models in continuous time, an integration function is required:
-// *)
-
-// type TimeMode<'data, 'time> =
-// | Discrete
-// | Continuous of Integrate<'data, 'time>
-
-// and Integrate<'data,'time> = 'time -> 'time -> 'time -> CodedMap<'data> -> CodedMap -> CodedMap<'data[]>
-
-// (**
-// Currently only fixed timesteps are supported, but variable timestep support (e.g. for sediment core data) is planned.
-
-// Two integration functions are included:
-
-// | Solver | Function | Description |
-// | - | - | - |
-// | Runge-Kutta 4 (MathNet Numerics) | ``Integration.MathNet.integrate`` | A fourth-order Runge Kutta method to provide approximate solutions to ODE systems. |
-// | Runge-Kutta 547M (Open Solving Library for ODEs - Microsoft Research) | ``Integration.MsftOslo.integrate`` | A method based on classic Runge-Kutta, but with automatic error and step size control. [See the documentation](https://www.microsoft.com/en-us/research/wp-content/uploads/2014/07/osloUserGuide.pdf).|
-
-// ## Optimisation
-
-// Bristlecone supports optimisation functions that have the following type signature:
-
-// ```fsharp
-// type Optimise<'data> = int -> int -> Domain -> ('data[] -> 'data) -> ('data * 'data []) list
-// ```
-
-// There are two optimsation techniques currently built-in:
-
-// | Method | Function | Description |
-// | - | - | - |
-// | Amoeba (Nelder-Mead) | ``Optimisation.Amoeba.solve`` | A gradient-descent method. |
-// | MCMC Random Walk | ``Optimisation.MonteCarlo.randomWalk`` | A method based on classic Runge-Kutta, but with automatic error and step size control. [See the documentation](https://www.microsoft.com/en-us/research/wp-content/uploads/2014/07/osloUserGuide.pdf).|
-
-
-// ## Estimation Engines
-
-// To use Bristlecone functions requires a configured ``EstimationEngine``. The easiest way is with the helper functions within the ``Bristlecone`` module:
-
-// | Function | Type | Description |
-// | - | - | - |
-// | ``mkContinuous`` | EstimationEngine | Default continuous engine |
-// | ``mkDiscrete`` | EstimationEngine | Default discrete model engine |
-// | ``withContinuousTime`` | t : Integrate<'a,'b> -> engine: EstimationEngine<'a,'b> -> EstimationEngine<'a,'b> | Transforms engine to continuous time mode, using the given integrator function. |
-// | ``withConditioning`` | c: Conditioning -> engine: EstimationEngine<'a,'b> -> EstimationEngine<'a,'b> | Choose how the start point is chosen when solving the model system. |
-
-// **)
\ No newline at end of file
diff --git a/docs/index.fsx b/docs/index.fsx
index f371d58..39a9094 100644
--- a/docs/index.fsx
+++ b/docs/index.fsx
@@ -5,6 +5,9 @@ category: Getting Started
categoryindex: 1
index: 1
---
+
+[]({{root}}/{{fsdocs-source-basename}}.fsx)
+[]({{root}}/{{fsdocs-source-basename}}.ipynb)
*)
(*** condition: prepare ***)
@@ -24,7 +27,8 @@ index: 1
Bristlecone
======================
-Bristlecone is a library for conducting model-fitting and model-selection analyses on time-series data.
+Bristlecone is a library for easily *composing* theoretical models into larger systems, and
+subsequently conducting model-fitting and model-selection analyses on time-series data.
Although originally designed for the investigation of non-linear dynamics within ecological
and environmental sciences, the library can be used across finance, econometrics and
other fields that apply non-linear modelling techniques.
@@ -45,23 +49,22 @@ The nuget package [is available here](https://nuget.org/packages/Bristlecone).
This example demonstrates the layout of a model when defined in Bristlecone.
*)
-open Bristlecone // Opens Bristlecone core library and estimation engine
-open Bristlecone.Language // Open the language for writing Bristlecone models
+open Bristlecone // Opens Bristlecone core library and estimation engine
+open Bristlecone.Language // Open the language for writing Bristlecone models
let hypothesis =
- let vonBertalanffy =
- Parameter "η" * This ** Parameter "β" - Parameter "κ" * This
+ let vonBertalanffy = Parameter "η" * This ** Parameter "β" - Parameter "κ" * This
Model.empty
- |> Model.addEquation "mass" vonBertalanffy
- |> Model.estimateParameter "η" noConstraints 0.50 1.50
- |> Model.estimateParameter "β" noConstraints 0.01 1.00
- |> Model.estimateParameter "κ" noConstraints 0.01 1.00
+ |> Model.addEquation "mass" vonBertalanffy
+ |> Model.estimateParameter "η" noConstraints 0.50 1.50
+ |> Model.estimateParameter "β" noConstraints 0.01 1.00
+ |> Model.estimateParameter "κ" noConstraints 0.01 1.00
|> Model.useLikelihoodFunction (ModelLibrary.Likelihood.sumOfSquares [ "mass" ])
|> Model.compile
-let engine =
+let engine =
Bristlecone.mkContinuous
|> Bristlecone.withCustomOptimisation (Optimisation.Amoeba.swarm 5 20 Optimisation.Amoeba.Solver.Default)
@@ -84,7 +87,7 @@ as an F# script or Jupyter notebook using the buttons at the top of each example
* [The predator-prey example](examples/predator-prey.html) covers basic model-fitting with Bristlecone.
- * [The shrub-resource example](example/shrub-resource.html) is a more
+ * [The shrub-resource example](examples/shrub-resource.html) is a more
comprehensive example that covers model-fitting and model-selection (MFMS) with Bristlecone.
* The [API Reference](reference/index.html) contains automatically generated documentation for all types, modules
diff --git a/docs/integration.fsx b/docs/integration.fsx
index a6f7eb2..4a7af1b 100644
--- a/docs/integration.fsx
+++ b/docs/integration.fsx
@@ -5,6 +5,9 @@ category: Components
categoryindex: 4
index: 3
---
+
+[]({{root}}/{{fsdocs-source-basename}}.fsx)
+[]({{root}}/{{fsdocs-source-basename}}.ipynb)
*)
(*** condition: prepare ***)
@@ -50,13 +53,12 @@ function to match the `EstimationEngine.Integration` type signature as follows:
open Bristlecone
-let myCustomIntegrator : EstimationEngine.Integrate<'data,'time> =
+let myCustomIntegrator: EstimationEngine.Integrate<'data, 'time> =
fun writeOut tInitial tEnd tStep initialConditions externalEnvironment model ->
invalidOp "Doesn't actually do anything!"
let engine =
- Bristlecone.mkContinuous
- |> Bristlecone.withContinuousTime myCustomIntegrator
+ Bristlecone.mkContinuous |> Bristlecone.withContinuousTime myCustomIntegrator
(**
@@ -64,4 +66,4 @@ When defining a custom integration routine, you may wish to use the
`Base.solve` function from the `Bristlecone.Integration` namespace.
This is used to arrange the equations into a form suitable for a
simple integration routine.
-*)
\ No newline at end of file
+*)
diff --git a/docs/language.fsx b/docs/language.fsx
index 8d80527..4810cd5 100644
--- a/docs/language.fsx
+++ b/docs/language.fsx
@@ -5,6 +5,9 @@ category: The Bristlecone Language
categoryindex: 1
index: 1
---
+
+[]({{root}}/{{fsdocs-source-basename}}.fsx)
+[]({{root}}/{{fsdocs-source-basename}}.ipynb)
*)
(*** condition: prepare ***)
@@ -29,8 +32,8 @@ namespace:
*)
-open Bristlecone // Opens Bristlecone core library and estimation engine
-open Bristlecone.Language // Open the language for writing Bristlecone models
+open Bristlecone // Opens Bristlecone core library and estimation engine
+open Bristlecone.Language // Open the language for writing Bristlecone models
(**
@@ -93,9 +96,10 @@ Please see the below example:
let linear =
Conditional(fun compute ->
- if compute (Parameter "a" * Constant 5.) < 1e-12
- then Invalid
- else Parameter "a" * This)
+ if compute (Parameter "a" * Constant 5.) < 1e-12 then
+ Invalid
+ else
+ Parameter "a" * This)
(**
In this example, the function is limited such that
@@ -125,8 +129,7 @@ to simplify this process within the `Arbitrary` module:
/// x: population size
/// t: current time
/// z: some environmental value
-let someArbitraryFn x t z =
- x + t * 1.0 - z
+let someArbitraryFn x t z = x + t * 1.0 - z
// TODO Insert example here
@@ -136,4 +139,4 @@ let someArbitraryFn x t z =
The `Time` keyword simply retrives the current time index; as such, it
can be used to formulate time-dependent models.
-*)
\ No newline at end of file
+*)
diff --git a/docs/logging.fsx b/docs/logging.fsx
deleted file mode 100644
index 9bb0f2e..0000000
--- a/docs/logging.fsx
+++ /dev/null
@@ -1,28 +0,0 @@
-(**
----
-title: Logging
-category: Components
-categoryindex: 4
-index: 3
----
-*)
-
-(*** condition: prepare ***)
-#nowarn "211"
-#r "../src/Bristlecone/bin/Debug/netstandard2.0/Bristlecone.dll"
-(*** condition: fsx ***)
-#if FSX
-#r "nuget: Bristlecone,{{fsdocs-package-version}}"
-#endif // FSX
-(*** condition: ipynb ***)
-#if IPYNB
-#r "nuget: Bristlecone,{{fsdocs-package-version}}"
-#endif // IPYNB
-
-(**
-Logging in Bristlecone
-========================
-
-Forthcoming
-
-*)
\ No newline at end of file
diff --git a/docs/model-fitting.fsx b/docs/model-fitting.fsx
index edae468..b63695a 100644
--- a/docs/model-fitting.fsx
+++ b/docs/model-fitting.fsx
@@ -5,6 +5,9 @@ category: Techniques
categoryindex: 1
index: 2
---
+
+[]({{root}}/{{fsdocs-source-basename}}.fsx)
+[]({{root}}/{{fsdocs-source-basename}}.ipynb)
*)
(*** condition: prepare ***)
@@ -25,4 +28,4 @@ Model Fitting
-*)
\ No newline at end of file
+*)
diff --git a/docs/model-selection.fsx b/docs/model-selection.fsx
index 0627721..9b74010 100644
--- a/docs/model-selection.fsx
+++ b/docs/model-selection.fsx
@@ -5,6 +5,9 @@ category: Techniques
categoryindex: 1
index: 3
---
+
+[]({{root}}/{{fsdocs-source-basename}}.fsx)
+[]({{root}}/{{fsdocs-source-basename}}.ipynb)
*)
(*** condition: prepare ***)
@@ -33,19 +36,28 @@ and limitations of alternative model selection statistics.
To calculate Akaike weights for a set of hypotheses, you must first obtain
your results by either loading in saved result files, or running models directly.
Once you have obtained your results, weights can be saved after calculation
-by using the functions within the `Bristlecone.Data` namespace as below:
+by using the functions within the `Bristlecone.Data` namespace.
+
+The `ModelSelection.save` function takes a sequence of tuples of four values: the subject's ID,
+the model or hypothesis ID, the estimated result, and the Akaike weight for that result. In the
+below example, we do some wrangling to get the weights into the correct format for saving:
*)
open Bristlecone
-fun results ->
-
- let resultsDirectory = "some/results/directory/"
+fun modelResults subjectIds hypothesisIds ->
+ modelResults
+ |> ModelSelection.Akaike.akaikeWeights
+ |> Seq.zip3 subjectIds hypothesisIds
+ |> Seq.map (fun (s, h, (a, b)) -> (s, h, a, b))
+ |> Bristlecone.Data.ModelSelection.save "/some/dir"
- let weights =
- results
- |> ModelSelection.weights
+(**
+If you are working with `ResultSet`s, calculating and saving the weights is easier, and can
+be completed as below:
+*)
- // Save the weights into a csv file
- weights
- |> Data.ModelSelection.save resultsDirectory
+fun (hypothesisResults: ModelSelection.ResultSet.ResultSet seq) ->
+ hypothesisResults
+ |> ModelSelection.Akaike.akaikeWeightsForSet (fun h -> h.ReferenceCode)
+ |> Bristlecone.Data.ModelSelection.save "/some/dir"
diff --git a/docs/model.fsx b/docs/model.fsx
index 3ebcf7d..61bec1e 100644
--- a/docs/model.fsx
+++ b/docs/model.fsx
@@ -5,6 +5,9 @@ category: The Bristlecone Language
categoryindex: 1
index: 2
---
+
+[]({{root}}/{{fsdocs-source-basename}}.fsx)
+[]({{root}}/{{fsdocs-source-basename}}.ipynb)
*)
(*** condition: prepare ***)
@@ -36,8 +39,8 @@ A `ModelSystem` is defined by:
To build a `ModelSystem`, first open the Bristlecone language:
*)
-open Bristlecone // Opens Bristlecone core library and estimation engine
-open Bristlecone.Language // Open the language for writing Bristlecone models
+open Bristlecone // Opens Bristlecone core library and estimation engine
+open Bristlecone.Language // Open the language for writing Bristlecone models
(**
Second, define the parts of your hypothesis. In this example,
@@ -72,10 +75,10 @@ The `ModelSystem` can be created with forward pipes (`|>`) as follows:
let hypothesis =
Model.empty
- |> Model.addEquation "N" ``dN/dt``
- |> Model.estimateParameter "r" noConstraints 0.10 5.00
- |> Model.estimateParameter "K" noConstraints 50.0 150.
- |> Model.useLikelihoodFunction (ModelLibrary.Likelihood.sumOfSquares [ "N" ])
+ |> Model.addEquation "N" ``dN/dt``
+ |> Model.estimateParameter "r" noConstraints 0.10 5.00
+ |> Model.estimateParameter "K" noConstraints 50.0 150.
+ |> Model.useLikelihoodFunction (ModelLibrary.Likelihood.sumOfSquares [ "N" ])
|> Model.compile
(**
@@ -89,7 +92,7 @@ required parts are in place.
Let's look at defining the estimatable parameters in more detail:
*)
-Model.estimateParameter "r" noConstraints 0.10 5.00
+Model.estimateParameter "r" noConstraints 0.10 5.00
(**
For the `Model.estimateParameter` function, first pass
@@ -105,4 +108,4 @@ from which to draw an initial starting value for the parameter.
The starting value is drawn from a continuous uniform distribution,
using a `Random` that you define in the `EstimationEngine` (see [estimating](estimation-engine.html)).
-*)
\ No newline at end of file
+*)
diff --git a/docs/optimisation.fsx b/docs/optimisation.fsx
index 1598f89..5c1292c 100644
--- a/docs/optimisation.fsx
+++ b/docs/optimisation.fsx
@@ -5,6 +5,9 @@ category: Components
categoryindex: 4
index: 1
---
+
+[]({{root}}/{{fsdocs-source-basename}}.fsx)
+[]({{root}}/{{fsdocs-source-basename}}.ipynb)
*)
(*** condition: prepare ***)
@@ -36,13 +39,12 @@ procedure can be used with Bristlecone by plugging in as follows:
open Bristlecone
-let myCustomOptimiser : EstimationEngine.Optimiser =
- EstimationEngine.InDetachedSpace <|
- fun writeOut n domain f -> invalidOp "Doesn't actually do anything!"
+let myCustomOptimiser: EstimationEngine.Optimiser =
+ EstimationEngine.InDetachedSpace
+ <| fun writeOut n domain f -> invalidOp "Doesn't actually do anything!"
let engine =
- Bristlecone.mkContinuous
- |> Bristlecone.withCustomOptimisation myCustomOptimiser
+ Bristlecone.mkContinuous |> Bristlecone.withCustomOptimisation myCustomOptimiser
(**
## Included optimisation methods
@@ -103,8 +105,7 @@ It can be used as follows:
let settingsNM = Amoeba.Solver.Default
(*** include-value: settingsNM ***)
-let single : EstimationEngine.Optimise =
- Amoeba.Solver.solve settingsNM
+let single: EstimationEngine.Optimise = Amoeba.Solver.solve settingsNM
(**
A single Nelder-Mead solver is highly subject to local minima. To reduce
@@ -179,10 +180,11 @@ tuning.
open Bristlecone.Optimisation.MonteCarlo.Filzbach
-let settingsFB = { TuneAfterChanges = 20
- MaxScaleChange = 100.
- MinScaleChange = 0.01
- BurnLength = EndConditions.afterIteration 100000 }
+let settingsFB =
+ { TuneAfterChanges = 20
+ MaxScaleChange = 100.
+ MinScaleChange = 0.01
+ BurnLength = EndConditions.afterIteration 100000 }
filzbach settingsFB
diff --git a/docs/workflow.fsx b/docs/workflow.fsx
index 44712cf..cd21981 100644
--- a/docs/workflow.fsx
+++ b/docs/workflow.fsx
@@ -5,6 +5,9 @@ category: Techniques
categoryindex: 2
index: 4
---
+
+[]({{root}}/{{fsdocs-source-basename}}.fsx)
+[]({{root}}/{{fsdocs-source-basename}}.ipynb)
*)
(*** condition: prepare ***)
@@ -53,14 +56,8 @@ let workPackages datasets hypotheses engine =
seq {
for d in datasets do
for h in [ 1 .. hypotheses |> List.length ] do
- for _ in [ 1 .. replicates ] do
- yield async {
- return Bristlecone.fit
- engine
- endCondition
- d
- hypotheses.[h-1]
- }
+ for _ in [ 1..replicates ] do
+ yield async { return Bristlecone.fit engine endCondition d hypotheses.[h - 1] }
}
(**
@@ -114,21 +111,21 @@ iteration, for each chain (along with process IDs). Next, let's create and setup
the orchestration agent:
*)
-let orchestrator = Orchestration.OrchestrationAgent(logger, System.Environment.ProcessorCount, false)
+let orchestrator =
+ Orchestration.OrchestrationAgent(logger, System.Environment.ProcessorCount, false)
+
+fun datasets hypotheses engine ->
-// fun datasets hypotheses engine ->
+ // Orchestrate the analyses
+ let work = workPackages datasets hypotheses engine
-// // Orchestrate the analyses
-// let work = workPackages datasets hypotheses engine
-// let run() =
-// work
-// |> Seq.iter (
-// Orchestration.OrchestrationMessage.StartWorkPackage
-// >> orchestrator.Post)
+ let run () =
+ work
+ |> Seq.iter (Orchestration.OrchestrationMessage.StartWorkPackage >> orchestrator.Post)
-// run()
+ run ()
(**
If the above code is supplied with datasets, hypotheses, and an
`EstimationEngine`, it will schedule, queue and run the jobs until complete.
-*)
\ No newline at end of file
+*)
diff --git a/samples/bristlecone.fsx b/samples/bristlecone.fsx
deleted file mode 100644
index b51e012..0000000
--- a/samples/bristlecone.fsx
+++ /dev/null
@@ -1,6 +0,0 @@
-#r "nuget: MathNet.Numerics.FSharp"
-#r "nuget: FSharp.Data"
-
-#r "../src/Bristlecone/bin/Debug/netstandard2.0/Microsoft.Research.Oslo.dll"
-#r "../src/Bristlecone/bin/Debug/netstandard2.0/Bristlecone.dll"
-#r "../src/Bristlecone.Dendro/bin/Debug/netstandard2.0/Bristlecone.Dendro.dll"
\ No newline at end of file
diff --git a/samples/readme.md b/samples/readme.md
deleted file mode 100644
index 95b1ff7..0000000
--- a/samples/readme.md
+++ /dev/null
@@ -1,5 +0,0 @@
-To run these samples:
-
-1. Install the latest .NET Core SDK. The minimum required version is 3.1.
-2. In the samples directory, run `dotnet fsi` from the terminal / command line.
-3. Send script linewise to F# interactive.
\ No newline at end of file
diff --git a/src/Bristlecone.Dendro/AssemblyInfo.fs b/src/Bristlecone.Dendro/AssemblyInfo.fs
index 7c1f75d..a568e78 100644
--- a/src/Bristlecone.Dendro/AssemblyInfo.fs
+++ b/src/Bristlecone.Dendro/AssemblyInfo.fs
@@ -1,5 +1,6 @@
// Auto-Generated by FAKE; do not edit
namespace System
+
open System.Reflection
[]
@@ -11,9 +12,20 @@ open System.Reflection
do ()
module internal AssemblyVersionInformation =
- let [] AssemblyTitle = "Bristlecone.Dendro"
- let [] AssemblyProduct = "Bristlecone"
- let [] AssemblyDescription = "Time-series modelling in F#"
- let [] AssemblyVersion = "1.0.0"
- let [] AssemblyFileVersion = "1.0.0"
- let [] AssemblyConfiguration = "Release"
+ []
+ let AssemblyTitle = "Bristlecone.Dendro"
+
+ []
+ let AssemblyProduct = "Bristlecone"
+
+ []
+ let AssemblyDescription = "Time-series modelling in F#"
+
+ []
+ let AssemblyVersion = "1.0.0"
+
+ []
+ let AssemblyFileVersion = "1.0.0"
+
+ []
+ let AssemblyConfiguration = "Release"
diff --git a/src/Bristlecone.Dendro/Library.fs b/src/Bristlecone.Dendro/Library.fs
index 5b1c205..648af48 100644
--- a/src/Bristlecone.Dendro/Library.fs
+++ b/src/Bristlecone.Dendro/Library.fs
@@ -5,94 +5,41 @@ open Bristlecone.Dendro
module Bristlecone =
- /// Fit a model to plant growth time-series. The plant individual's
- /// growth data is always labelled as `x`.
- let fitDendro engine iterations system (plant: PlantIndividual.PlantIndividual) =
+ /// Lower time-series data from a plant individual type into
+ /// basic time-series that may be used for model-fitting
+ /// A plant individual
+ /// A coded map of time-series data for model-fitting
+ let fromDendro (plant: PlantIndividual.PlantIndividual) =
let g =
match plant.Growth |> PlantIndividual.growthSeries with
| GrowthSeries.Absolute g -> g
| GrowthSeries.Cumulative g -> g
| GrowthSeries.Relative g -> g
- let predictors = plant.Environment |> Map.add (ShortCode.create "x" |> Option.get) g
- Bristlecone.fit engine iterations predictors system
+ plant.Environment |> Map.add (ShortCode.create "x" |> Option.get) g
- /// Perform n-step-ahead computation on the hypothesis and plant.
- /// The plant individual's growth data is always labelled as `x`.
- let predictAheadDendro
+ /// Fit a model to plant growth time-series. The plant individual's growth data is always labelled as `x`.
+ /// A configured estimation engine
+ /// The end condition to apply to optimisation
+ /// The compiled model system for analysis
+ /// A plant individual
+ /// An estimation result for the given plant, model and fitting method (engine)
+ let fitDendro engine endCondition system (plant: PlantIndividual.PlantIndividual) =
+ Bristlecone.fit engine endCondition (fromDendro plant) system
+
+ /// Perform n-step-ahead computation on the hypothesis and plant. The plant individual's growth data is always labelled as `x`.
+ ///
+ ///
+ ///
+ ///
+ ///
+ ///
+ let oneStepAheadDendro
(engine: EstimationEngine.EstimationEngine)
system
(plant: PlantIndividual.PlantIndividual)
preTransform
estimate
=
- let g =
- match plant.Growth |> PlantIndividual.growthSeries with
- | GrowthSeries.Absolute g -> g
- | GrowthSeries.Cumulative g -> g
- | GrowthSeries.Relative g -> g
-
- let predictors = plant.Environment |> Map.add (ShortCode.create "x" |> Option.get) g
+ let predictors = fromDendro plant
Bristlecone.oneStepAhead engine system preTransform predictors estimate
-
-//
-//module Test =
-//
-// // How to process plants?
-// // 1. Where time-series have variable start and end dates, bound to common time.
-// // 2. Determine conditioning based on real data at new t-1, or can be calculated by allometrics.
-// // 3.
-//
-// type GrowthModellingBasis =
-// | Biomass
-// | AbsoluteGrowthRate
-// | RelativeGrowthRate
-//
-// let startValues (startDate:System.DateTime) (plant:PlantIndividual.PlantIndividual) =
-// let initialRadius =
-// match plant.Growth with
-// | PlantIndividual.PlantGrowth.RingWidth s ->
-// match s with
-// | GrowthSeries.Absolute c -> c.Head |> fst |> removeUnit
-// | GrowthSeries.Cumulative c ->
-// let start = (c |> TimeSeries.trimStart (startDate - System.TimeSpan.FromDays(366.))).Values |> Seq.head |> removeUnit
-// printfn "Start cumulative growth = %f" start
-// start
-// | GrowthSeries.Relative _ -> invalidOp "Not implemented"
-// | _ -> invalidOp "Not implemented 2"
-// let initialMass = initialRadius |> removeUnit |> ModelComponents.Proxies.toBiomassMM
-// let initialNitrogen = plant.Environment.[ShortCode.create "N"].Head |> fst
-// let initialT = plant.Environment.[code "T[max]"] |> TimeSeries.findExact startDate |> fst
-// [ (ShortCode.create "x", initialRadius)
-// (ShortCode.create "N", initialNitrogen)
-// (ShortCode.create "T[max]", initialT)
-// (ShortCode.create "bs", initialMass) ] |> Map.ofList
-//
-// let preparePlantAnalysis growthBasis localEnvironment plant =
-//
-// // 1. Add environment to plant
-// let plant2 =
-// localEnvironment
-// |> List.fold(fun plant (c,ts) -> plant |> PlantIndividual.zipEnv c ts) plant
-//
-// // 2. Transform growth into required basis
-// let plant3 =
-// match growthBasis with
-// | Biomass -> plant |> PlantIndividual.toCumulativeGrowth
-// | _ -> failwith "Not finished"
-//
-// // 3.
-//
-// let shrub = s |> PlantIndividual.toCumulativeGrowth |> PlantIndividual.zipEnv (code "T[max]" |> Option.get) (TemperatureData.monthly)
-// let common = shrub |> PlantIndividual.keepCommonYears |> PlantIndividual.zipEnv (code "T[max]" |> Option.get) (TemperatureData.monthly)
-// let startDate = (common.Environment.[code "N" |> Option.get]).StartDate |> snd
-// let startConditions = getStartValues startDate shrub
-// let e = engine |> Bristlecone.withConditioning (Conditioning.Custom startConditions)
-//
-//
-// failwith "Cool"
-//
-//
-// type StartStrategy =
-// | CommonTime
-// | FirstGrowthYear
diff --git a/src/Bristlecone.Dendro/Sunrise.fs b/src/Bristlecone.Dendro/Sunrise.fs
index 73884c8..ff84232 100644
--- a/src/Bristlecone.Dendro/Sunrise.fs
+++ b/src/Bristlecone.Dendro/Sunrise.fs
@@ -2,6 +2,7 @@ namespace Bristlecone.Dendro
open System
+/// Provides functions for converting to and from dates in the Julian calendar
module JulianDate =
let minutesInDay = 24 * 60
@@ -47,12 +48,12 @@ module JulianDate =
(jDate - (float J2000)) / daysInCentury
-/// The sunrise equation can be used to calculate the
+/// The sunrise equation can be used to calculate the
/// time of sunrise and sunset for any latitude, longitude,
-/// and date.
-/// See: https://en.wikipedia.org/wiki/Sunrise_equation#Complete_calculation_on_Earth
+/// and date.
+/// See: https://en.wikipedia.org/wiki/Sunrise_equation#Complete_calculation_on_Earth
/// Adapted from the SolarCalc by NOAA.
-/// Source: https://dotnetfiddle.net/N3j5th
+/// Source: https://dotnetfiddle.net/N3j5th
module Sunrise =
type DayLength =
@@ -168,10 +169,16 @@ module Sunrise =
let sunsetOffset = JulianDate.toDate timeZone gDate sunset
(sunriseOffset, sunsetOffset) |> PartialLight
+ /// A cache that may be used in model computation where day lengths
+ /// are continually accessed for the same dates over and over again. Stashes
+ /// each calculated day length per latitude/longitude and time.
type DayLengthCache(latitude, longitude, timeZone) =
let mutable (lightData: Map) = [] |> Map.ofList
+ /// Daylight hours
+ ///
+ ///
member __.GetLight(date) =
match lightData |> Map.tryFind date with
| Some l -> l
diff --git a/src/Bristlecone/AssemblyInfo.fs b/src/Bristlecone/AssemblyInfo.fs
index 3184497..436ddbf 100644
--- a/src/Bristlecone/AssemblyInfo.fs
+++ b/src/Bristlecone/AssemblyInfo.fs
@@ -1,12 +1,12 @@
// Auto-Generated by FAKE; do not edit
namespace System
+
open System.Reflection
[]
-[]
+[]
[]
[]
[]
[]
do ()
-
diff --git a/src/Bristlecone/Bristlecone.fsproj b/src/Bristlecone/Bristlecone.fsproj
index ba9fd95..9c70795 100755
--- a/src/Bristlecone/Bristlecone.fsproj
+++ b/src/Bristlecone/Bristlecone.fsproj
@@ -12,10 +12,10 @@
-
+
-
+
diff --git a/src/Bristlecone/ConfidenceInterval.fs b/src/Bristlecone/Confidence.fs
similarity index 91%
rename from src/Bristlecone/ConfidenceInterval.fs
rename to src/Bristlecone/Confidence.fs
index 5d99caf..e26d4ea 100755
--- a/src/Bristlecone/ConfidenceInterval.fs
+++ b/src/Bristlecone/Confidence.fs
@@ -1,11 +1,11 @@
-namespace Bristlecone.Optimisation.ConfidenceInterval
+namespace Bristlecone.Confidence
open Bristlecone.Logging
///
/// Contains functions for calculating confidence intervals on model fits.
///
-///
+///
/// The 95% and 68% confidence interval around a point estimate.
type ConfidenceInterval =
{ ``95%``: Interval
@@ -14,8 +14,8 @@ type ConfidenceInterval =
and Interval = { Lower: float; Upper: float }
-/// Keeps a list of all optimisation events that occur for further analysis.
-type OptimisationEventStash() =
+/// Keeps a list of all optimisation events that occur for further analysis.
+type internal OptimisationEventStash() =
let mutable events = []
member __.SaveEvent(e) =
@@ -26,21 +26,21 @@ type OptimisationEventStash() =
member __.GetAll() = events
-/// Differences in likelihood for different confidence intervals
-/// based on a chi squared distribution.
+/// Differences in likelihood for different confidence intervals
+/// based on a chi squared distribution.
module Bounds =
- /// The difference in likelihood at 68% confidence
+ /// The difference in likelihood at 68% confidence
let lowerBound = 0.49447324 // qchisq(0.68,1)/2
- /// The difference in likelihood at 95% confidence
+ /// The difference in likelihood at 95% confidence
let upperBound = 1.92072941 // qchisq(0.95,1)/2
-/// Given a maximum likelihood estimate (MLE), the profile likelihood method
-/// runs a Monte Carlo algorithm that samples around the MLE. The range for
-/// each parameter is discovered at 95% and 68% confidence based on a chi squared
-/// distribution.
+/// Given a maximum likelihood estimate (MLE), the profile likelihood method
+/// runs a Monte Carlo algorithm that samples around the MLE.
+/// The range for each parameter is discovered at 95% and 68%
+/// confidence based on a chi squared distribution.
module ProfileLikelihood =
open Bristlecone
@@ -51,7 +51,6 @@ module ProfileLikelihood =
type EstimateFunction<'data, 'time, 'subject> =
EstimationEngine<'data, 'time> -> ModelSystem -> 'subject -> EstimationResult
-
module CustomOptimisationMethod =
type TuneSettings =
@@ -72,7 +71,7 @@ module ProfileLikelihood =
// 1. Initial conditions
let draw' = jump random
- let theta1 = Bristlecone.Optimisation.MonteCarlo.initialise domain random
+ let theta1 = Bristlecone.Optimisation.Initialise.initialise domain random
let l1 = f theta1
// 2. Tune individual parameter scales (for many univariate distributions)
diff --git a/src/Bristlecone/DataAccess.fs b/src/Bristlecone/Data.fs
similarity index 62%
rename from src/Bristlecone/DataAccess.fs
rename to src/Bristlecone/Data.fs
index e9743fb..86462f5 100755
--- a/src/Bristlecone/DataAccess.fs
+++ b/src/Bristlecone/Data.fs
@@ -1,19 +1,18 @@
namespace Bristlecone.Data
-// Bristlecone loads and saves CSV data in these formats:
-// 1. [Invididual] MLE: best-fitting parameter steps and their likelihood
-// 2. [Individual] Trace: record of optimisation steps
-// 3. [Invididual] Series: best-fitting parameter steps and their likelihood
-// 4. [Ensemble] Model-selection (weights)
-
open Bristlecone
open Bristlecone.ModelSystem
open FSharp.Data
///
-/// Contains functions that enable loading and saving of core Bristlecone types to file storage.
+///
+/// Contains functions that enable loading and saving of core Bristlecone
+/// types to file storage. Bristlecone loads and saves CSV data for both
+/// individual model results (MLEs, traces, predicted vs observed series),
+/// and for ensemble-level statistics (model selection / weights).
+///
///
-///
+///
/// Specifies how file paths are constructed for different data types.
module Config =
@@ -22,11 +21,13 @@ module Config =
| Trace
| Series
| Intervals
+ | StepAhead of steps: int
| Components
type EnsembleType =
| Weights
| Convergence
+ | RMSE
let typeAsLabel dataType =
match dataType with
@@ -35,11 +36,13 @@ module Config =
| Series -> "series"
| Intervals -> "ci"
| Components -> "components"
+ | StepAhead(steps) -> sprintf "%i-step-ahead-prediction" steps
let ensembleAsLabel dataType =
match dataType with
| Weights -> "ensemble-weights"
| Convergence -> "ensemble-convergence"
+ | RMSE -> "ensemble-rmse"
let filePath directory subject modelId (resultId: System.Guid) dataType =
let path = System.IO.DirectoryInfo(directory)
@@ -75,6 +78,7 @@ module Config =
| Series -> "series"
| Intervals -> "ci"
| Components -> "components"
+ | StepAhead(steps) -> sprintf "%istepahead" steps
let files = path.GetFiles(sprintf "bristlecone-%s-%s-%s-*.csv" subject modelId t)
let regex = sprintf "bristlecone-%s-%s-%s-(%s).csv" subject modelId t regexGuid
@@ -91,14 +95,15 @@ module Config =
else
invalidArg "directory" "The specified directory does not exist"
+/// Functions for loading and saving optimisation traces
[]
module Trace =
type BristleconeTrace = CsvProvider<"templates/individual-trace.csv", IgnoreErrors=true>
- module Row =
+ module internal Row =
- let internal fromEstimate thinBy subject modelId (result: EstimationResult) : seq =
+ let fromEstimate thinBy subject modelId (result: EstimationResult) : seq =
result.Trace
|> Seq.rev
|> Seq.mapi (fun iterationNumber (likelihood, values) ->
@@ -107,10 +112,13 @@ module Trace =
|> Seq.mapi (fun i (name, _) ->
(subject, modelId, iterationNumber + 1, result.ResultId, name.Value, likelihood, values.[i])
|> BristleconeTrace.Row))
- |> Seq.everyNth thinBy
+ |> if (thinBy |> Option.isSome) then
+ Seq.everyNth thinBy.Value
+ else
+ id
|> Seq.concat
- let internal toTrace (data: BristleconeTrace) : (float * float[]) list =
+ let toTrace (data: BristleconeTrace) : (float * float[]) list =
data.Rows
|> Seq.groupBy (fun r -> r.Iteration)
|> Seq.map (fun (i, r) ->
@@ -140,12 +148,13 @@ module Trace =
| 0 -> None
| _ -> (i, data |> Row.toTrace) |> Some)
+/// Functions of loading and saving Maximum Likelihood Estimates from model fits
[]
module MLE =
type IndividualMLE = CsvProvider<"templates/individual-mle.csv", IgnoreErrors=true>
- module Row =
+ module internal Row =
let fromResult subject hypothesisId (result: EstimationResult) =
result.Parameters
@@ -159,18 +168,28 @@ module MLE =
v |> Parameter.getTransformedValue)
|> IndividualMLE.Row)
- let toResult (data: IndividualMLE) =
+ let toResult (modelSystem: ModelSystem) (data: IndividualMLE) =
if data.Rows |> Seq.isEmpty then
Error "An MLE file is corrupt"
else
let mle = (data.Rows |> Seq.head).NegativeLogLikelihood
- let pool =
+ let estimatedTheta =
data.Rows
|> Seq.choose (fun r ->
ShortCode.create r.ParameterCode |> Option.map (fun o -> o, r.ParameterValue))
|> Map.ofSeq
+ let pool =
+ modelSystem.Parameters
+ |> Parameter.Pool.toList
+ |> List.map (fun (k, v) -> k, Parameter.setTransformedValue v (estimatedTheta |> Map.find k))
+ |> List.choose (fun (c, r) ->
+ match r with
+ | Ok x -> Some(c, x)
+ | Error _ -> None)
+ |> Parameter.Pool.fromList
+
(mle, pool) |> Ok
let save directory subject modelId result =
@@ -181,22 +200,41 @@ module MLE =
csv.Save(filePath)
- let load directory subject modelId =
- let traceFiles = Config.fileMatch directory subject modelId Config.DataType.MLE
-
- traceFiles
+ /// Load the Maximum Likelihood Estimates
+ /// from the given directory for a given subject and model combination.
+ /// The results folder where files are saved
+ /// An identifier for the subject of the fit
+ /// An identifier for the model that was fit
+ /// A sequence of tuples which contain the analysis ID followed by another tuple
+ /// that contains the likelihood and theta (parameter set)
+ let load directory subject (modelSystem: ModelSystem) modelId =
+ let files = Config.fileMatch directory subject modelId Config.DataType.MLE
+
+ files
|> Seq.choose (fun (i, f) ->
let data = IndividualMLE.Load f.FullName
match data.Rows |> Seq.length with
| 0 -> None
| _ ->
- let parsed = data |> Row.toResult
+ let parsed = data |> Row.toResult modelSystem
match parsed with
| Ok mle -> (i, mle) |> Some
| Error _ -> None)
+ /// Load the Maximum Likelihood Estimate (with the lowest -log Liklihood)
+ /// from the given directory for a given subject and model combination.
+ /// The results folder where files are saved
+ /// An identifier for the subject of the fit
+ /// An identifier for the model that was fit
+ /// A tuple containing the analysis ID followed by another tuple
+ /// that contains the likelihood and theta (parameter set)
+ let loadBest directory subject modelSystem modelId =
+ load directory subject modelSystem modelId
+ |> Seq.minBy (fun (_, (mle, _)) -> mle)
+
+/// Functions for loading and saving predicted vs observed time-series for model fits
[]
module Series =
@@ -204,7 +242,7 @@ module Series =
type IndividualSeries = CsvProvider<"templates/individual-series.csv", IgnoreErrors=true>
- module Row =
+ module internal Row =
let fromResult subject hypothesisId (result: EstimationResult) =
result.Series
@@ -249,8 +287,14 @@ module Series =
[]
module EstimationResult =
- /// Save the Maximum Likelihood Estimate, trace of the optimisation
- /// procedure, and time-series.
+ /// Save the Maximum Likelihood Estimate, trace of the optimisation
+ /// procedure, and time-series.
+ /// Relative or absolute directory to save files to
+ /// An identifier for the subject of the test
+ /// An identifier for the model used
+ /// If Some, an integer representing the nth traces to keep from the optimisation trace.
+ /// If None, does not thin the trace.
+ /// The estimation result to save
let saveAll directory subject modelId thinTraceBy result =
Trace.save directory subject modelId thinTraceBy result
MLE.save directory subject modelId result
@@ -261,7 +305,8 @@ module EstimationResult =
/// when file names and formats are in original Bristlecone format.
let loadAll directory subject (modelSystem: ModelSystem) modelId =
let mles =
- MLE.load directory subject modelId |> Seq.map (fun (k, v) -> k.ToString(), v)
+ MLE.load directory subject modelSystem modelId
+ |> Seq.map (fun (k, v) -> k.ToString(), v)
let series =
Series.load directory subject modelId |> Seq.map (fun (k, v) -> k.ToString(), v)
@@ -276,15 +321,7 @@ module EstimationResult =
|> Seq.map (fun (k, t, (s, (l, p))) ->
{ ResultId = k |> System.Guid.Parse
Likelihood = l
- Parameters =
- modelSystem.Parameters
- |> Parameter.Pool.toList
- |> List.map (fun (k, v) -> k, Parameter.setTransformedValue v (p |> Map.find k))
- |> List.choose (fun (c, r) ->
- match r with
- | Ok x -> Some(c, x)
- | Error _ -> None)
- |> Parameter.Pool.fromList
+ Parameters = p
Series = s
InternalDynamics = None
Trace = t })
@@ -292,7 +329,7 @@ module EstimationResult =
[]
module Confidence =
- open Bristlecone.Optimisation.ConfidenceInterval
+ open Bristlecone.Confidence
type IndividualCI = CsvProvider<"templates/individual-mle-ci.csv">
@@ -340,9 +377,8 @@ module Convergence =
type ConvergenceStats = CsvProvider<"templates/ensemble-convergence.csv">
- let toCsvRows (result: seq) =
+ let toCsvRows (result: seq) =
result
- |> Seq.concat
|> Seq.map (fun r ->
ConvergenceStats.Row(r.Subject, r.HypothesisId, r.Parameter.Value, r.StatisticName, r.StatisticValue))
@@ -350,3 +386,60 @@ module Convergence =
let csv = new ConvergenceStats(result |> toCsvRows)
let filePath = Config.filePathEnsemble directory Config.EnsembleType.Convergence
csv.Save(filePath)
+
+/// Loads and saves one-step-ahead and n-step-ahead prediction
+/// results into and out of simple CSV files
+[]
+module NStepAhead =
+
+ type NStepAhead = CsvProvider<"templates/individual-n-step-ahead.csv">
+ type NStepAheadStats = CsvProvider<"templates/ensemble-rmse.csv">
+
+ let internal toCsvRows
+ subjectId
+ hypothesisId
+ analysisId
+ nSteps
+ (result: CodedMap)
+ =
+ result
+ |> Seq.collect (fun r ->
+ r.Value
+ |> fst
+ |> Time.TimeSeries.toObservations
+ |> Seq.map (fun (fit, t) ->
+ NStepAhead.Row(subjectId, hypothesisId, analysisId, t, fit.Obs, nSteps, fit.Fit)))
+
+ let internal toStatCsvRows
+ (results: seq>)
+ =
+ results
+ |> Seq.collect (fun (s, h, r) ->
+ r
+ |> Seq.map (fun kv -> NStepAheadStats.Row(s, h, kv.Key.Value, "RMSE", (snd kv.Value).RMSE)))
+
+ /// Save an individual n-step prediction (excluding statistics) for an
+ /// individual subject and hypothesis / model.
+ /// The file directory to save to
+ /// Subject ID
+ /// Model ID or Hypothesis ID
+ /// The unique reference of the estimation result used to generate n-step predictions
+ /// The number of time-steps ahead predicted
+ /// The n-step prediction from Bristlecone.oneStepAhead or similar
+ let save directory subjectId modelId analysisId stepsAhead result =
+ let csv =
+ new NStepAhead(result |> toCsvRows subjectId modelId analysisId stepsAhead)
+
+ let filePath =
+ Config.filePath directory subjectId modelId analysisId (Config.DataType.StepAhead stepsAhead)
+
+ csv.Save(filePath)
+
+ /// Saves a full ensemble overview of the root mean squared error
+ /// of n-step ahead predictions.
+ /// The file directory to save to
+ /// A sequence of subjectID * hypothesisID * n-step result
+ let saveAllRMSE directory results =
+ let csv = new NStepAheadStats(results |> toStatCsvRows)
+ let filePath = Config.filePathEnsemble directory (Config.EnsembleType.RMSE)
+ csv.Save(filePath)
diff --git a/src/Bristlecone/Diagnostics.fs b/src/Bristlecone/Diagnostics.fs
index fa617a4..7a2a174 100644
--- a/src/Bristlecone/Diagnostics.fs
+++ b/src/Bristlecone/Diagnostics.fs
@@ -4,12 +4,12 @@ module Bristlecone.Diagnostics
open Bristlecone.ModelSystem
-/// Convergence diagnostics for monte-carlo markov chain (MCMC) analyses.
+/// Convergence diagnostics for monte-carlo markov chain (MCMC) analyses.
module Convergence =
open Bristlecone.ModelSelection
- /// A per-parameter convergence statistic. The statistic used is given in `StatisticName`.
+ /// A per-parameter convergence statistic. The statistic used is given in `StatisticName`.
type ConvergenceStatistic =
{ Subject: string
HypothesisId: string
@@ -17,19 +17,28 @@ module Convergence =
StatisticName: string
StatisticValue: float }
- /// Calculate the Gelman-Rubin statistic for each parameter in the given
+ /// Calculate the Gelman-Rubin statistic for each parameter in the given
/// `ResultSet`. The statistic tends downwards to one, with one indicating
- /// perfect convergence between all chains.
- let gelmanRubin n (result: ResultSet.ResultSet) =
- let subjectId, _, hi, r = result
- printfn "Calculating Rhat for %s %s" subjectId hi
-
- match r with
- | Some(results, mle) ->
+ /// perfect convergence between all chains.
+ /// How many recent iterations to use from the trace.
+ /// The subject identifier
+ /// The hypothesis identifier
+ /// A result set (of 1 .. many results) for a particular subject and hypothesis
+ /// If more than one replicate, the R-hat convergence statistic across replicates
+ let gelmanRubin nMostRecent subjectId hypothesisId (result: ResultSet.ResultSet<'subject, 'hypothesis>) =
+ printfn "Calculating Rhat for %s (H = %s)" subjectId hypothesisId
+
+ if result.AllResults |> Seq.isEmpty || result.BestResult.IsNone then
+ None
+ else
let chains =
- results
+ result.AllResults
|> Seq.sortBy (fun x -> x.Likelihood)
- |> (fun c -> if c |> Seq.length > n then c |> Seq.take n else c)
+ |> (fun c ->
+ if c |> Seq.length > nMostRecent then
+ c |> Seq.take nMostRecent
+ else
+ c)
|> Seq.map (fun r -> r.Trace |> Seq.map snd)
let minChainLength = chains |> Seq.minBy (fun x -> x |> Seq.length) |> Seq.length
@@ -39,78 +48,95 @@ module Convergence =
else
let chains' = chains |> Seq.map (fun c -> c |> Seq.take minChainLength)
- mle.Parameters
+ result.BestResult.Value.Parameters
|> Parameter.Pool.toList
|> Seq.mapi (fun i (code, p) ->
{ Subject = subjectId
- HypothesisId = hi
+ HypothesisId = hypothesisId
Parameter = code
- StatisticName = "rHat"
+ StatisticName = "R-hat"
StatisticValue =
chains'
|> Seq.map (fun c -> c |> Seq.map (fun x -> x.[i]))
|> Statistics.Convergence.GelmanRubin.rHat })
|> Some
- | None -> None
-
- /// Calculate the Gelman-Rubin statistic for each parameter in all of the
- /// given `ResultSet`. The statistic tends downwards to one, with one indicating
- /// perfect convergence between all chains.
- let gelmanRubinAll n (results: ResultSet.ResultSet list) =
- results |> Seq.choose (gelmanRubin n)
-
-/// Logging functions to output the internal dynamics of model systems.
+ /// Calculate the Gelman-Rubin statistic for each parameter in all of the given
+ /// `ResultSet`. The statistic tends downwards to one, with one indicating
+ /// perfect convergence between all chains.
+ /// How many recent iterations to use from the trace.
+ /// A function to retrieve a subject ID from a subject
+ /// A function to retrieve a hypothesis ID from a hypothesis
+ /// A result set (of 1 .. many results) for a particular subject and hypothesis
+ /// If more than one replicate, the R-hat convergence statistic across replicates
+ let gelmanRubinAll nMostRecent subject hypothesis (results: ResultSet.ResultSet<'subject, 'hypothesis> seq) =
+ results
+ |> Seq.choose (fun r -> gelmanRubin nMostRecent (subject r.Subject) (hypothesis r.Hypothesis) r)
+ |> Seq.concat
+
+
+/// Logging functions to output the internal dynamics of model systems.
module ModelComponents =
open Bristlecone.ModelSelection
open Bristlecone.EstimationEngine
- type IComponentLogger<'data> =
- abstract member StoreValue: string -> float -> 'data -> 'data
- abstract member GetAll: unit -> Map>
-
- /// A component logger stores the value at each time t for each
- /// component specified by a `componentId`.
- type ComponentLogger<'data>() =
- let mutable (data: Map>) = [] |> Map.ofList
-
- interface IComponentLogger<'data> with
-
- member __.StoreValue componentId t v =
- let existing = data |> Map.tryFind componentId
-
- match existing with
- | Some e -> data <- data |> Map.add componentId (e |> Map.add t v)
- | None -> data <- data |> Map.add componentId ([ t, v ] |> Map.ofList)
-
- v
-
- member __.GetAll() = data
-
- /// A component logger that does not store any values.
- type PassThrough<'data>() =
- interface IComponentLogger<'data> with
- member __.StoreValue _ _ v = v
- member __.GetAll() = [] |> Map.ofList
-
- /// Log out components specified in a model by disabling optimisation.
- /// The model will only be computed once.
+ /// Contains types that can collect internal
+ /// components of a model during computation.
+ module Loggers =
+
+ /// A component logger can store value for a
+ /// component's value at a particular time.
+ type IComponentLogger<'data> =
+ abstract member StoreValue: componentName: string -> time: float -> value: 'data -> 'data
+ abstract member GetAll: unit -> Map>
+
+ /// A component logger stores the value at each time t for each
+ /// component specified by a `componentId`.
+ type ComponentLogger<'data>() =
+ let mutable (data: Map>) = [] |> Map.ofList
+
+ interface IComponentLogger<'data> with
+
+ member __.StoreValue componentId t v =
+ let existing = data |> Map.tryFind componentId
+
+ match existing with
+ | Some e -> data <- data |> Map.add componentId (e |> Map.add t v)
+ | None -> data <- data |> Map.add componentId ([ t, v ] |> Map.ofList)
+
+ v
+
+ member __.GetAll() = data
+
+ /// A component logger that does not store any values.
+ type PassThrough<'data>() =
+ interface IComponentLogger<'data> with
+ member __.StoreValue _ _ v = v
+ member __.GetAll() = [] |> Map.ofList
+
+ /// Log out components specified in a model by disabling optimisation. The model will only be computed once.
+ /// The function used
+ ///
+ ///
+ ///
+ ///
+ ///
+ ///
+ ///
let calculateComponents
fitFn
engine
- (result: ResultSet.ResultSet<'subject, IComponentLogger<'data> -> ModelSystem>)
+ (result: ResultSet.ResultSet<'subject, Loggers.IComponentLogger<'data> -> ModelSystem>)
=
- let subject, hypothesis, _, estimate = result
-
- match estimate with
+ match result.BestResult with
| None -> [] |> Map.ofList
- | Some(_, mle) ->
+ | Some mle ->
let eng =
{ engine with
OptimiseWith = Optimisation.None.none }
- let cLog = ComponentLogger<'data>() :> IComponentLogger<'data>
+ let cLog = Loggers.ComponentLogger<'data>() :> Loggers.IComponentLogger<'data>
let p =
mle.Parameters
@@ -123,6 +149,9 @@ module ModelComponents =
|> Option.map (fun v -> k, v))
|> Parameter.Pool.fromList
- let mleHypothesis = { hypothesis cLog with Parameters = p }
- fitFn subject mleHypothesis eng |> ignore
+ let mleHypothesis =
+ { result.Hypothesis cLog with
+ Parameters = p }
+
+ fitFn result.Subject mleHypothesis eng |> ignore
cLog.GetAll()
diff --git a/src/Bristlecone/Equations.fs b/src/Bristlecone/Equations.fs
index 13d95a7..d355c28 100755
--- a/src/Bristlecone/Equations.fs
+++ b/src/Bristlecone/Equations.fs
@@ -4,7 +4,7 @@ namespace Bristlecone.ModelLibrary
///
///
/// Pre-built model parts for use in Bristlecone
-///
+///
module Likelihood =
open Bristlecone
diff --git a/src/Bristlecone/EstimationEngine.fs b/src/Bristlecone/EstimationEngine.fs
index 5e24bd5..9039c05 100755
--- a/src/Bristlecone/EstimationEngine.fs
+++ b/src/Bristlecone/EstimationEngine.fs
@@ -20,8 +20,11 @@ module ModelSystem =
{ 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
+ type ParameterValueAccessor =
+ | ParameterValueAccessor of (string -> float)
+
+ member this.Get name =
+ let (ParameterValueAccessor v) = this in v name
/// A function that computes the likelihood of a set of parameters.
type LikelihoodFn = ParameterValueAccessor -> CodedMap -> float
@@ -34,7 +37,7 @@ module ModelSystem =
{ Parameters: Parameter.Pool
Equations: CodedMap
Measures: CodedMap
- Likelihood: LikelihoodFn }
+ NegLogLikelihood: LikelihoodFn }
type FitValue = { Fit: float; Obs: float }
type FitSeries = TimeSeries
diff --git a/src/Bristlecone/Language.fs b/src/Bristlecone/Language.fs
index b33a140..d3e66ac 100644
--- a/src/Bristlecone/Language.fs
+++ b/src/Bristlecone/Language.fs
@@ -46,8 +46,11 @@ module Language =
| Divide of ModelExpression * ModelExpression
| Arbitrary of (float -> float -> Parameter.Pool -> CodedMap -> float) * list
| Mod of ModelExpression * ModelExpression
- | Exponent of ModelExpression * ModelExpression
+ | Power of ModelExpression * ModelExpression
+ | Logarithm of ModelExpression
+ | Exponential of ModelExpression
| Conditional of ((ModelExpression -> float) -> ModelExpression)
+ | Label of string * ModelExpression
| Invalid
static member (*)(e1, e2) =
@@ -62,7 +65,9 @@ module Language =
static member (-)(e1, e2) = Subtract(e1, e2)
static member (%)(e1, remainder) = Mod(e1, remainder)
- static member Pow(e1, pow) = Exponent(e1, pow)
+ static member Pow(e1, pow) = Power(e1, pow)
+ static member Log(e) = Logarithm(e)
+ static member Exp(e) = Exponential(e)
static member (~-) e =
Multiply[e
@@ -84,7 +89,10 @@ module Language =
| Parameter name ->
match pool |> Parameter.Pool.tryGetRealValue name with
| Some est -> est
- | None -> failwithf "The equation could not be calculated. The parameter '%s' has not been set up (or has yet to be estimated)." name
+ | None ->
+ failwithf
+ "The equation could not be calculated. The parameter '%s' has not been set up (or has yet to be estimated)."
+ name
| Constant n -> n
| Add list ->
match list with
@@ -98,10 +106,14 @@ module Language =
| Divide(l, r) -> compute x t pool environment l / compute x t pool environment r
| Arbitrary(fn, _) -> fn x t pool environment
| Mod(e, m) -> (compute x t pool environment e) % (compute x t pool environment m)
- | Exponent(e, m) -> (compute x t pool environment e) ** (compute x t pool environment m)
+ | Power(e, m) -> (compute x t pool environment e) ** (compute x t pool environment m)
+ | Logarithm e -> log (compute x t pool environment e)
+ | Exponential e -> exp (compute x t pool environment e)
| Conditional m -> m (compute x t pool environment) |> compute x t pool environment
+ | Label(_, m) -> m |> compute x t pool environment
| Invalid -> nan
+
module Writer =
type Writer<'a, 'L> = AWriter of 'a * List<'L>
@@ -145,10 +157,13 @@ module Language =
| _ -> failwith "List was empty"
| Divide(l, r) -> describe l |> flatMap (fun () -> describe r)
| Arbitrary(fn, reqs) -> bind ((), "custom component - TODO may use additional parameters?")
+ | Logarithm(e) -> describe e
+ | Exponential(e) -> describe e
| Mod(e, _) -> describe e
- | Exponent(e, _) -> describe e
+ | Power(e, _) -> describe e
| Invalid -> bind ((), "invalid model")
| Conditional _ -> bind ((), "conditional element") // TODO
+ | Label(_, e) -> describe e
type Requirement =
| ParameterRequirement of string
@@ -163,24 +178,23 @@ module Language =
| Parameter name -> ParameterRequirement name :: reqs
| Constant _ -> reqs
| Add list
- | Multiply list ->
- list
- |> List.collect(fun l -> requirements l reqs)
- |> List.append reqs
+ | Multiply list -> list |> List.collect (fun l -> requirements l reqs) |> List.append reqs
| Divide(l, r)
- | Subtract(l, r) ->
- [ requirements l reqs; requirements r reqs; reqs ] |> List.concat
+ | Subtract(l, r) -> [ requirements l reqs; requirements r reqs; reqs ] |> List.concat
| Arbitrary(fn, r) ->
- r
- |> List.map(fun r ->
+ r
+ |> List.map (fun r ->
match r with
| ArbitraryEnvironment e -> EnvironmentRequirement e
| ArbitraryParameter p -> ParameterRequirement p)
|> List.append reqs
| Mod(e, _) -> requirements e reqs
- | Exponent(e, _) -> requirements e reqs
+ | Power(e, _) -> requirements e reqs
+ | Logarithm e -> requirements e reqs
+ | Exponential e -> requirements e reqs
| Invalid -> reqs
| Conditional _ -> reqs
+ | Label(_, e) -> requirements e reqs
/// Allows common F# functions to use Bristlecone model expressions.
@@ -290,26 +304,28 @@ module Language =
|> Seq.choose id
|> Map.ofSeq
- if Seq.hasDuplicates (Seq.concat [ Map.keys measures; Map.keys equations ])
- then failwith "Duplicate keys were used within equation and measures. These must be unique."
+ if Seq.hasDuplicates (Seq.concat [ Map.keys measures; Map.keys equations ]) then
+ failwith "Duplicate keys were used within equation and measures. These must be unique."
- if equations.IsEmpty then failwith "No equations specified. You must state at least one model equation."
+ if equations.IsEmpty then
+ failwith "No equations specified. You must state at least one model equation."
- equations
+ equations
|> Map.map (fun _ v -> ExpressionParser.requirements v [])
|> Map.toList
|> List.map snd
|> List.collect id
|> List.distinct
- |> List.iter(fun req ->
+ |> List.iter (fun req ->
match req with
| ExpressionParser.ParameterRequirement p ->
match parameters |> Parameter.Pool.hasParameter p with
| Some p -> ()
- | None -> failwithf "The specified model requires the parameter '%s' but this has not been set up." p
+ | None ->
+ failwithf "The specified model requires the parameter '%s' but this has not been set up." p
| ExpressionParser.EnvironmentRequirement _ -> ())
- { Likelihood = likelihoods |> Seq.head
+ { NegLogLikelihood = likelihoods |> Seq.head
Parameters = parameters
Equations = equations |> Map.map (fun _ v -> (fun pool t x env -> compute x t pool env v))
Measures = measures }
@@ -379,6 +395,10 @@ module Language =
Parameters = (snd comp).Parameters |> Map.add c p }
| None -> failwithf "The code '%s' cannot be used as an identifier. See docs." name
+ /// Types to represent a hypothesis, given that a hypothesis
+ /// is a model system that contains some alternate formulations of certain
+ /// components.
+ []
module Hypotheses =
open Writer
@@ -386,17 +406,41 @@ module Language =
let internal nFirstLetters n (s: string) =
(s.Split(' ') |> Seq.map (fun s -> s |> Seq.truncate n) |> string).ToUpper()
+ /// Represents the name of a swappable component in a model
+ /// and the name of its implementation in a specific case.
type ComponentName =
{ Component: string
Implementation: string }
+ /// A short code that may be used to indicate a particular
+ /// model component by its current implementation.
+ /// A string in the format XX_XXX with the first part
+ /// indicating the model component and the second part the implementation.
member this.Reference =
sprintf "%s_%s" (nFirstLetters 2 this.Component) (nFirstLetters 3 this.Implementation)
- /// Implement a component where a model system requires one. A component is
- /// a part of a model that may be varied, for example between competing
- /// hypotheses.
- let createFromComponent comp x =
+ /// A hypothesis consists of a model system and the names of the
+ /// swappable components within it, alongside the name of their current implementation.
+ type Hypothesis =
+ | Hypothesis of ModelSystem.ModelSystem * list
+
+ member private this.unwrap = this |> fun (Hypothesis(h1, h2)) -> (h1, h2)
+
+ /// Compiles a reference code that may be used to identify (although not necessarily uniquely) this hypothesis
+ /// A string in the format XX_XXX_YY_YYY... where XX_XXX is a singe component with XX the component and XXX the implementation.
+ member this.ReferenceCode =
+ this.unwrap |> snd |> List.map (fun c -> c.Reference) |> String.concat "_"
+
+ member this.Components = this.unwrap |> snd
+ member this.Model = this.unwrap |> fst
+
+ /// Implement a component where a model system requires one. A component is a part of a model that may be varied, for example between competing hypotheses.
+ /// A tuple representing a component scaffolded with the `modelComponent` and `subComponent` functions.
+ /// A builder started with `createFromComponent`
+ ///
+ ///
+ /// A builder to add further components or compile with `Hypothesis.compile`
+ let createFromComponent comp builder =
if snd comp |> List.isEmpty then
failwithf "You must specify at least one implementation for the '%s' component" (fst comp)
@@ -406,29 +450,39 @@ module Language =
comp
|> snd
- |> List.map (fun (n, c) -> bind (x c.Expression, (name n, c.Parameters)))
-
- /// Implement a second or further component on a model system where one is
- /// still required.
- let useAnother comp x =
+ |> List.map (fun (n, c) -> bind (builder c.Expression, (name n, c.Parameters)))
+
+ /// Add a second or further component to a model system where one is still required.
+ /// A tuple representing a component scaffolded with the `modelComponent` and `subComponent` functions.
+ /// A builder started with `createFromComponent`
+ ///
+ ///
+ /// A builder to add further components or compile with `Hypothesis.compile`
+ let useAnother comp builder =
let name i =
{ Component = fst comp
Implementation = i }
- List.allPairs (snd comp) x
+ List.allPairs (snd comp) builder
|> List.map (fun ((n, c), model) ->
model |> flatMap (fun m -> bind (m c.Expression, (name n, c.Parameters))))
- let useAnotherWithName comp x =
+ /// Add a second or further component to a model system where one is still required.
+ /// A tuple representing a component scaffolded with the `modelComponent` and `subComponent` functions.
+ /// A builder started with `createFromComponent`
+ ///
+ ///
+ /// A builder to add further components or compile with `Hypothesis.compile`
+ let useAnotherWithName comp builder =
let name i =
{ Component = fst comp
Implementation = i }
- List.allPairs (snd comp) x
+ List.allPairs (snd comp) builder
|> List.map (fun ((n, c), model) ->
model |> flatMap (fun m -> bind (m (n, c.Expression), (name n, c.Parameters))))
- /// Adds parameters from model components into the base model builder.
+ /// Adds parameters from model components into the base model builder.
let internal addParameters
((modelBuilder, newParams): ModelBuilder.ModelBuilder * List>)
=
@@ -438,14 +492,42 @@ module Language =
(fun mb (name, p) -> mb |> ModelBuilder.add name.Value (ModelBuilder.ParameterFragment p))
modelBuilder
- /// Compiles a suite of competing model hypotheses based on the given components.
+ /// Compiles a suite of competing model hypotheses based on the given components.
/// The compilation includes only the required parameters in each model hypothesis,
- /// and combines all labels into a single model identifier.
- let compile (x: Writer> list) =
- if x |> List.isEmpty then
+ /// and combines all labels into a single model identifier.
+ /// A builder started with `createFromComponent`
+ /// A list of compiled hypotheses for this model system and specified components.
+ let compile (builder: Writer> list) =
+ if builder |> List.isEmpty then
failwith "No hypotheses specified"
- //x |> List.map (map Model.compile >> run)
- x
+
+ builder
|> List.map (fun h ->
let names = run h
- names |> addParameters |> Model.compile, names |> snd |> List.map fst)
+
+ (names |> addParameters |> Model.compile, names |> snd |> List.map fst)
+ |> Hypothesis)
+
+
+// type Hypotheses<'subject, 'hypothesis> =
+// | Hypotheses of Hypothesis<'subject, 'hypothesis> seq
+
+// let private unwrap (ResultSetMany m) = m
+
+// let many sets = ResultSetMany sets
+
+// let subject s sets =
+// sets
+// |> unwrap
+// |> Seq.filter(fun s -> s.Subject = s)
+
+// /// Map a function over results on a per-subject and per-hypothesis basis
+// /// A function to map over the subject-hypothesis groups
+// /// A result set (many)
+// /// A transformed results set
+// let map fn (sets:ResultSetMany<'subject, 'hypothesis>) =
+// sets
+// |> unwrap
+// |> Seq.groupBy(fun s -> s.Subject, s.Hypothesis)
+// |> Seq.map(fun (s, h) -> fn s h)
+// |> ResultSetMany
diff --git a/src/Bristlecone/Library.fs b/src/Bristlecone/Library.fs
index 0e802ce..d1756d7 100755
--- a/src/Bristlecone/Library.fs
+++ b/src/Bristlecone/Library.fs
@@ -3,14 +3,21 @@ namespace Bristlecone
open System
open Bristlecone.Logging
+///
+/// The core library of Bristlecone, containing model-fitting functions.
+///
+///
+/// Main functionality of Bristlecone, including functions to scaffold
+/// `ModelSystem`s and for model-fitting (tests and real fits).
[]
module Bristlecone =
open Bristlecone.Time
open Bristlecone.ModelSystem
open Bristlecone.EstimationEngine
+ open Bristlecone.Statistics
- /// A standard estimation engine using a random-walk monte carlo optimiser.
+ /// A basic estimation engine for discrete-time equations, using a Nelder-Mead optimiser.
let mkDiscrete: EstimationEngine =
{ TimeHandling = Discrete
OptimiseWith = Optimisation.Amoeba.single Optimisation.Amoeba.Solver.Default
@@ -18,7 +25,7 @@ module Bristlecone =
Random = MathNet.Numerics.Random.MersenneTwister(true)
Conditioning = Conditioning.NoConditioning }
- /// A standard `EstimationEngine` for ordinary differential equation models.
+ /// A basic estimation engine for ordinary differential equations, using a Nelder-Mead optimiser.
let mkContinuous =
{ TimeHandling = Continuous <| Integration.MathNet.integrate
OptimiseWith = Optimisation.Amoeba.single Optimisation.Amoeba.Solver.Default
@@ -26,12 +33,19 @@ module Bristlecone =
Random = MathNet.Numerics.Random.MersenneTwister(true)
Conditioning = Conditioning.RepeatFirstDataPoint }
- /// Add a writer
+ /// Substitute a specific logger into
+ ///
+ ///
+ ///
+ ///
+ ///
let withOutput out engine = { engine with LogTo = out }
- /// Use a mersenne twister random number generator
+ /// Use a mersenne twister random number generator
/// with a specific seed.
- let withSeed seed engine = { engine with Random = MathNet.Numerics.Random.MersenneTwister(seed, true) }
+ let withSeed seed engine =
+ { engine with
+ Random = MathNet.Numerics.Random.MersenneTwister(seed, true) }
/// Use a custom integration method
let withContinuousTime t engine =
@@ -51,7 +65,7 @@ module Bristlecone =
let withCustomOptimisation optim engine = { engine with OptimiseWith = optim }
- module Fit =
+ module internal Fit =
/// Places a map of `TimeSeries` into a `TimeFrame` that has data
/// that shares a common timeline. If no timeline is shared, returns
@@ -129,24 +143,29 @@ module Bristlecone =
///
///
///
- let fit engine endCondition timeSeriesData (model: ModelSystem) =
+ let tryFit engine endCondition timeSeriesData (model: ModelSystem) =
// A. Setup initial time point values based on conditioning method.
let t0 = Fit.t0 timeSeriesData engine.Conditioning engine.LogTo
// Check there is time-series data actually included and corresponding to correct equations.
let hasRequiredData =
- if timeSeriesData.IsEmpty then Error "No time-series data was specified"
+ if timeSeriesData.IsEmpty then
+ Error "No time-series data was specified"
+ else if Set.isSubset (model.Equations |> Map.keys |> set) (timeSeriesData |> Map.keys |> set) then
+ Ok timeSeriesData
else
- if Set.isSubset (model.Equations |> Map.keys |> set) (timeSeriesData |> Map.keys |> set)
- then Ok timeSeriesData
- else Error (sprintf "Required time-series data were missing. Need: %A" (model.Equations |> Map.keys |> Seq.map (fun k -> k.Value) |> String.concat " + "))
+ Error(
+ sprintf
+ "Required time-series data were missing. Need: %A"
+ (model.Equations |> Map.keys |> Seq.map (fun k -> k.Value) |> String.concat " + ")
+ )
// B. Create a continuous-time that outputs float[]
// containing only the values for the dynamic variable resolution.
let continuousSolver =
result {
-
+
let! timeSeriesData = hasRequiredData
// 1. Set time-series into common timeline
@@ -256,16 +275,29 @@ module Bristlecone =
InternalDynamics = Some estimatedHighRes }
}
+ /// Fit a time-series model to data.
+ /// An estimation engine configured and tested for the given model.
+ /// The condition at which optimisation should cease.
+ /// Time-series dataset that contains a series for each equation in the model system.
+ /// A model system of equations, likelihood function, estimatible parameters, and optional measures.
+ /// The result of the model-fitting procedure. If an error occurs, throws an exception.
+ let fit engine endCondition timeSeriesData (model: ModelSystem) =
+ tryFit engine endCondition timeSeriesData model |> Result.forceOk
open Test
- /// **Description**
- /// Test that the specified estimation engine can correctly estimate known parameters. Random parameter sets are generated from the given model system.
- /// **Parameters**
- /// * `model` - a `ModelSystem` of equations and parameters
- /// * `testSettings` - settings
- /// * `engine` - an `EstimationEngine`
- let testModel engine (settings: Test.TestSettings) (model: ModelSystem) : Result =
+ /// Tests that the specified estimation engine can correctly
+ /// estimate known parameters given specfici test settings.
+ /// Random parameter sets and resultant fake time-series data are generated
+ /// for the model system by using the rules and noise generation settings
+ /// in the stated test settings.
+ ///
+ ///
+ ///
+ /// A test result that indicates the error structure.
+ /// It is wrapped in an F# Result, indicating if the procedure
+ /// was successful or not.
+ let tryTestModel engine (settings: Test.TestSettings) (model: ModelSystem) =
engine.LogTo <| GeneralEvent "Attempting to generate parameter set."
engine.LogTo
@@ -280,7 +312,7 @@ module Bristlecone =
let! settings = Test.isValidSettings model settings
let! trueData, theta = Test.Compute.tryGenerateData engine settings model settings.Attempts
- let! realEstimate =
+ let realEstimate =
fit
{ engine with
OptimiseWith = Optimisation.None.none }
@@ -289,7 +321,7 @@ module Bristlecone =
{ model with
Parameters = Parameter.Pool.fromEstimated theta }
- let! estimated =
+ let estimated =
fit engine settings.EndCondition (Map.merge trueData settings.EnvironmentalData (fun x y -> x)) model
let paramDiffs: Test.ParameterTestResult list =
@@ -328,30 +360,30 @@ module Bristlecone =
EstimatedLikelihood = estimated.Likelihood }
}
- /// **Description**
- /// Repeat a model fit many times, removing a single data point at random each time.
- /// **Parameters**
- /// * `engine` - parameter of type `EstimationEngine`
- /// * `iterations` - parameter of type `int`
- /// * `burnin` - parameter of type `int`
- /// * `bootstrapCount` - parameter of type `int`
- /// * `hypothesis` - parameter of type `ModelSystem`
- /// * `identifier` - parameter of type `ShortCode`
- /// * `series` - parameter of type `CodedMap>`
- ///
- /// **Output Type**
- /// * `EstimationResult list`
- ///
- /// **Exceptions**
- ///
- let bootstrap engine random iterations bootstrapCount hypothesis (identifier: ShortCode.ShortCode) series =
+ /// Test that the specified estimation engine can correctly
+ /// estimate known parameters. Random parameter sets are generated
+ /// from the given model system.
+ /// An estimation engine containing the method used for model-fitting.
+ /// Test settings that define how the test will be conducted.
+ /// The model system to test against the estimation engine.
+ /// A test result that indicates differences between the expected and actual fit.
+ let testModel engine settings model =
+ tryTestModel engine settings model |> Result.forceOk
+
+ /// Repeat a model fit many times, removing a single data point at random each time.
+ /// The estimation engine / fitting method
+ /// The end condition for each model fit
+ /// Number of times to bootstrap data
+ /// A model system / hypothesis to fit
+ /// Time-series to fit with model
+ /// A list of estimation results (one for each bootstrap) for further analysis
+ let bootstrap (engine: EstimationEngine.EstimationEngine) endCondition bootstrapCount model series =
let rec bootstrap s numberOfTimes solutions =
if (numberOfTimes > 0) then
- let subset = Bootstrap.removeSingle random s
- let result = fit engine iterations subset hypothesis
+ let subset = Bootstrap.removeSingle engine.Random s
+ let result = fit engine endCondition subset model
- engine.LogTo
- <| GeneralEvent(sprintf "%s: completed bootstrap %i" identifier.Value numberOfTimes)
+ engine.LogTo <| GeneralEvent(sprintf "Completed bootstrap %i" numberOfTimes)
bootstrap s (numberOfTimes - 1) (solutions |> List.append [ result ])
else
@@ -359,15 +391,24 @@ module Bristlecone =
bootstrap series bootstrapCount []
- /// "How good am I at predicting the next data point"?
- ///
+ /// Addresses the question: "How good am I at predicting the next data point?. Given fitted parameters,
+ /// assesses how well the model predicts the next data point from each point in the time-series data.
+ /// This approach can provide another indication of model performance.
+ /// The exact estimation engine used for the existing model fit
+ /// The exact model system / hypothesis from which parameters have been already estimated
+ /// A function that may transform each shorter time-series before prediction. This may be needed,
+ /// for example, if there are custom start values that need to be configured in a complex way (e.g. for derived mesaurement
+ /// variables).
+ /// The observed data to predict against.
+ /// A parameter pool containing already estimated parameters from model fitting step
+ /// A time-series for each variable containing a step-ahead prediction
let oneStepAhead
engine
hypothesis
- (preTransform: CodedMap> -> CodedMap>)
+ (preTransform: CodedMap> -> CodedMap>)
(timeSeries)
(estimatedTheta: Parameter.Pool)
- =
+ : CodedMap =
let hypothesisMle: ModelSystem =
{ hypothesis with
Parameters = Parameter.Pool.fromEstimated estimatedTheta }
@@ -386,21 +427,52 @@ module Bristlecone =
seq { 1..timeParcelCount }
|> Seq.map (fun i -> pairedDataFrames |> Map.map (fun _ v -> v |> Seq.item (i - 1)) |> preTransform)
- // TODO Remove this hack:
- // First data point is repeated, then skipped when returned
+ // It is predicting with a repeated first point:
+ // The next point estimate is at t1
+ // The next point observation is at t2
data
- |> Seq.map (fun d ->
- fit
- (engine
- |> withCustomOptimisation Optimisation.None.none
- |> withConditioning Conditioning.RepeatFirstDataPoint)
- (Optimisation.EndConditions.afterIteration 0)
+ |> Seq.collect (fun d ->
+ let est =
+ fit
+ (engine
+ |> withCustomOptimisation Optimisation.None.none
+ |> withConditioning Conditioning.RepeatFirstDataPoint)
+ (Optimisation.EndConditions.afterIteration 0)
+ d
+ hypothesisMle
+
+ let nextObservation =
d
- hypothesisMle)
- |> Seq.toList
+ |> Map.map (fun c ts -> ts |> TimeSeries.toObservations |> Seq.skip 1 |> Seq.head)
+ let paired =
+ nextObservation
+ |> Seq.map (fun kv ->
+ let nextEstimate = (est.Series.[kv.Key].Values |> Seq.head).Fit
+ (kv.Key,
+ { Obs = kv.Value |> fst
+ Fit = nextEstimate },
+ kv.Value |> snd))
+
+ paired)
+ |> Seq.groupBy (fun (k, _, _) -> k)
+ |> Seq.map (fun (tsName, values) ->
+ let sos = values |> Seq.averageBy (fun (_, x, _) -> (x.Obs - x.Fit) ** 2.)
+ tsName, (values |> Seq.map (fun (_, v, t) -> (v, t)) |> TimeSeries.fromObservations, { RMSE = sqrt sos }))
+ |> Map.ofSeq
+
+
+ /// Wrappers for fitting functions that use `Array.Parallel` to run many analyses at once.
module Parallel =
- let fit engine endCondition model growth =
- growth |> Array.Parallel.map (fun g -> fit engine endCondition g model)
+ /// A wrapper for `Bristlecone.fit` that uses Array.Parallel to run many analyses at once.
+ /// An estimation engine
+ /// An end condition
+ /// The model / hypothesis to fit
+ /// Time-series data to fit with the model
+ /// A list of estimation results
+ let fit engine endCondition model timeSeries =
+ timeSeries
+ |> Array.Parallel.map (fun g -> fit engine endCondition g model)
+ |> Array.toList
diff --git a/src/Bristlecone/ModelSelection.fs b/src/Bristlecone/ModelSelection.fs
index 8df219c..37a5e5c 100755
--- a/src/Bristlecone/ModelSelection.fs
+++ b/src/Bristlecone/ModelSelection.fs
@@ -1,34 +1,72 @@
+/// Contains tools for conducting Model Selection across
+/// individual subjects and hypotheses.
module Bristlecone.ModelSelection
open Bristlecone
open Bristlecone.ModelSystem
-/// Organises multiple hypotheses and multiple subjects into
-/// distinct analysis groups.
+/// Organises multiple hypotheses and multiple subjects into distinct analysis groups.
[]
module ResultSet =
- /// A representation of all results for a particular subject and hypothesis
+ /// A representation of all results for a particular subject and hypothesis
type ResultSet<'subject, 'hypothesis> =
- ('subject * 'hypothesis * string * (EstimationResult seq * EstimationResult) option)
-
- /// Arrange estimation results into subject and hypothesis groups.
- let arrangeResultSets subjects hypotheses getResults : ResultSet<'a, 'hypothesis> seq =
- Seq.allPairs subjects (hypotheses |> Seq.mapi (fun i v -> (i + 1, v)))
- |> Seq.map (fun (s, (hi, h)) ->
- let r = getResults s h hi
+ { Subject: 'subject
+ Hypothesis: 'hypothesis
+ BestResult: EstimationResult option
+ AllResults: EstimationResult seq }
+
+ /// Arrange estimation results into subject and hypothesis groups.
+ ///
+ ///
+ ///
+ ///
+ ///
+ /// An estimation result
+ ///
+ let arrangeResultSets<'subject, 'hypothesis> (subjects: 'subject seq) (hypotheses: 'hypothesis seq) getResults =
+ Seq.allPairs subjects hypotheses
+ |> Seq.map (fun (s, h) ->
+ let r = getResults s h
if r |> Seq.isEmpty then
- (s, h, hi.ToString(), None)
+ { Subject = s
+ Hypothesis = h
+ BestResult = None
+ AllResults = [] }
else
let r' = r |> Seq.filter (fun x -> not (System.Double.IsNaN(x.Likelihood)))
if Seq.isEmpty r' then
- (s, h, hi.ToString(), None)
+ { Subject = s
+ Hypothesis = h
+ BestResult = None
+ AllResults = [] }
else
- (s, h, hi.ToString(), (r', r' |> Seq.minBy (fun x -> x.Likelihood)) |> Some))
+ { Subject = s
+ Hypothesis = h
+ BestResult = r' |> Seq.minBy (fun x -> x.Likelihood) |> Some
+ AllResults = r' })
+
+/// Runs a model comparison statistic across a sequence of
+/// `ResultSet`s. Uses the best MLE for each subject * hypothesis group
+/// an runs `comparisonFn` across these results.
+/// The subject, hypothesis code, the original result, and the statistic.
+let internal comparisonStatistic comparisonFn getRefCode (results: ResultSet.ResultSet<'subject, 'hypothesis> seq) =
+ results
+ |> Seq.groupBy (fun resultSet -> getRefCode resultSet.Hypothesis)
+ |> Seq.collect (fun (_, r) ->
+ let weights =
+ r |> Seq.choose (fun resultSet -> resultSet.BestResult) |> comparisonFn
-/// Functions for conducting Akaike Information Criterion (AIC).
+ r
+ |> Seq.zip weights
+ |> Seq.map (fun (w, resultSet) -> resultSet.Subject, getRefCode resultSet.Hypothesis, fst w, snd w)
+ |> Seq.toList)
+ |> Seq.toList
+
+
+/// Functions for conducting Akaike Information Criterion (AIC).
module Akaike =
type AkaikeWeight =
@@ -37,26 +75,23 @@ module Akaike =
AICc: float
Weight: float }
- ///**Description**
- /// The Akaike information criterion, a standardised index of model fit quality for models that have different numbers of parameters.
- ///**Parameters**
- /// * `k` - The number of parameters within the model in question.
- /// * `logLikelihood` - a `float` representing the minimum log-likelihood achieved for the model in question.
+ /// The Akaike information criterion, a standardised index of model fit quality for models that have different numbers of parameters.
+ /// The number of parameters within the model in question.
+ /// a `float` representing the minimum log-likelihood achieved for the model in question.
+ ///
let aic (k: int) logLikelihood = 2. * logLikelihood + 2. * (float k)
-
- ///**Description**
- /// The Akaike information criterion, corrected for small sample sizes.
- /// It represents standardised index of model fit quality for models that have different numbers of parameters.
- ///**Assumptions**
- /// Your model must adhere to the following assumptions:
+ /// The Akaike information criterion, corrected for small sample sizes.
+ /// It represents standardised index of model fit quality for models that have different numbers of parameters.
+ /// Your model must adhere to the following assumptions:
/// - Univariate
/// - Linear in parameters
/// - Normally-distributed residuals
- ///**Parameters**
- /// * `n` - The sample size
- /// * `k` - The number of parameters within the model in question.
- /// * `logLikelihood` - a `float` representing the minimum log-likelihood achieved for the model in question.
+ ///
+ /// The sample size
+ /// The number of parameters within the model in question
+ /// A `float` representing the minimum log-likelihood achieved for the model in question.
+ ///
let aicc n k logLikelihood =
let aic = aic k logLikelihood
@@ -81,15 +116,10 @@ module Akaike =
AICc = aicc n p l
Weight = relative / totalLikelihood })
- /// **Description**
- /// Akaike weights for a set of model results.
- /// The weights can be directly interpreted as conditional probabilities for each model.
- ///
- /// **Output Type**
- /// * A `seq` of estimation results paired to their Akaike weights.
- ///
- /// **Exceptions**
- /// * `ArgumentException` - occurs when there are no observations within an estimation result.
+ /// Akaike weights for a sequence of `EstimationResult`s.
+ /// The input model results
+ /// An (EstimationResult * float) sequence of estimation results paired to their Akaike weights.
+ /// Occurs when there are no observations within an estimation result.
let akaikeWeights (models: seq) =
match models |> Seq.tryHead with
| None -> seq []
@@ -101,41 +131,10 @@ module Akaike =
|> akaikeWeights' n
|> Seq.zip models
-
-/// A record of an individual maximum likelihood estimate for
-/// a particular subject and hypothesis.
-type Result<'subject> =
- { AnalysisId: System.Guid
- Subject: 'subject
- ModelId: string
- Estimate: EstimationResult }
-
-/// Given a list of model predictions, find the best MLE for each
-/// model * subject combination, calculate the weights for this set.
-let calculate (results: seq>) =
- results
- |> Seq.groupBy (fun r -> (r.Subject, r.ModelId))
- |> Seq.map (fun (_, r) -> r |> Seq.minBy (fun x -> x.Estimate.Likelihood))
- |> Seq.groupBy (fun r -> r.Subject)
- |> Seq.collect (fun (_, r) ->
- let weights =
- r |> Seq.map (fun r -> r.Estimate) |> Akaike.akaikeWeights |> Seq.map snd
-
- weights |> Seq.zip r)
-
-///
-let weights (results: ResultSet.ResultSet<'a, 'b> seq) =
- results
- |> Seq.groupBy (fun (identifier, _, _, _) -> identifier)
- |> Seq.collect (fun (_, r) ->
- let weights =
- r
- |> Seq.choose (fun (_, _, _, r) -> r)
- |> Seq.map (fun (_, mle) -> mle)
- |> Akaike.akaikeWeights
-
- r
- |> Seq.zip weights
- |> Seq.map (fun (w, (s, _, hi, _)) -> s, hi, fst w, snd w)
- |> Seq.toList)
- |> Seq.toList
+ /// Akaike weights for a result set.
+ /// A function that gets a short reference code from a hypothesis.
+ /// A sequence of `ResultSet`s, within each the 1 .. many results of a particular subject * hypothesis combination.
+ /// An `(EstimationResult * float) seq` of estimation results paired to their Akaike weights.
+ /// Occurs when there are no observations within an estimation result.
+ let akaikeWeightsForSet getRefCode set =
+ comparisonStatistic akaikeWeights getRefCode set
diff --git a/src/Bristlecone/Objective.fs b/src/Bristlecone/Objective.fs
index f6a732e..ea358ed 100644
--- a/src/Bristlecone/Objective.fs
+++ b/src/Bristlecone/Objective.fs
@@ -47,7 +47,9 @@ module Objective =
point
|> predict system integrate solveDiscrete
|> pairObservationsToExpected observed
- |> system.Likelihood(
- let pool = point |> Parameter.Pool.fromPointInTransformedSpace system.Parameters
- ParameterValueAccessor
- <| fun name -> Parameter.Pool.tryGetRealValue name pool |> Option.get)
+ |> system.NegLogLikelihood(
+ let pool = point |> Parameter.Pool.fromPointInTransformedSpace system.Parameters
+
+ ParameterValueAccessor
+ <| fun name -> Parameter.Pool.tryGetRealValue name pool |> Option.get
+ )
diff --git a/src/Bristlecone/Optimisation.fs b/src/Bristlecone/Optimisation.fs
index c1ecdb8..7cceef1 100755
--- a/src/Bristlecone/Optimisation.fs
+++ b/src/Bristlecone/Optimisation.fs
@@ -65,8 +65,8 @@ module EndConditions =
stationarySquaredJumpDistance' 200 5
/// Convergence of results using the Gelman-Rubin Rhat statistic.
- /// * `thin` - Only test for convergence at multiples of the following intervals (when all chains are ready).
- /// * `chainCount` - The number of chains to test for convergence. This makes the agent wait until results for all chains are in.
+ /// `thin` - Only test for convergence at multiples of the following intervals (when all chains are ready).
+ /// `chainCount` - The number of chains to test for convergence. This makes the agent wait until results for all chains are in.
let convergence thin chainCount : EndCondition =
let assess (chains: Map list>) =
@@ -134,14 +134,9 @@ module EndConditions =
else
false
-
-/// A module containing Monte Carlo Markov Chain (MCMC) methods for optimisation.
-/// An introduction to MCMC approaches is provided by
-/// [Reali, Priami, and Marchetti (2017)](https://doi.org/10.3389/fams.2017.00006)
-module MonteCarlo =
+module Initialise =
open Bristlecone.Statistics.Distributions
- open MathNet.Numerics.LinearAlgebra
/// Generate a random point from bounds specified as a `Domain`.
/// A value for each dimension is drawn from a univariate normal distribution, assuming that
@@ -156,13 +151,13 @@ module MonteCarlo =
(Normal.draw rng (max + (min - max) / 2.) ((min - max) / 6.)) () |]
/// Assesses if theta is valid based on the provided
- /// constraints.
+ /// constraints.
let isInvalidTheta theta constraints =
Seq.zip theta constraints
- |> Seq.map(fun (v,c) ->
+ |> Seq.map (fun (v, c) ->
match c with
| Bristlecone.Parameter.Constraint.Unconstrained -> true
- | Bristlecone.Parameter.Constraint.PositiveOnly -> v > 0. )
+ | Bristlecone.Parameter.Constraint.PositiveOnly -> v > 0.)
|> Seq.contains false
/// Attempts to generate random theta based on starting bounds
@@ -173,15 +168,23 @@ module MonteCarlo =
Error "Could not generate a starting point given the domain"
else
let t = initialise domain random
- let isInvalid = isInvalidTheta t (domain |> Seq.map (fun (_,_,c) -> c))
+ let isInvalid = isInvalidTheta t (domain |> Seq.map (fun (_, _, c) -> c))
+
if isInvalid then
tryGenerateTheta f domain random (n - 1)
+ else if System.Double.IsNaN(f t) || System.Double.IsInfinity(f t) then
+ tryGenerateTheta f domain random (n - 1)
else
- if System.Double.IsNaN(f t) || System.Double.IsInfinity(f t)
- then
- tryGenerateTheta f domain random (n - 1)
- else
- Ok t
+ Ok t
+
+
+/// A module containing Monte Carlo Markov Chain (MCMC) methods for optimisation.
+/// An introduction to MCMC approaches is provided by
+/// [Reali, Priami, and Marchetti (2017)](https://doi.org/10.3389/fams.2017.00006)
+module MonteCarlo =
+
+ open Bristlecone.Statistics.Distributions
+ open MathNet.Numerics.LinearAlgebra
/// Jump in parameter space while reflecting constraints.
let constrainJump initial jump (scaleFactor: float) c =
@@ -193,26 +196,22 @@ module MonteCarlo =
else
(initial + (jump * scaleFactor))
-
- /// A recursive metropolis hastings algorithm that ends when `endCondition` returns true.
- ///
- /// **Parameters**
- /// * `random` - `System.Random` to be used for drawing from a uniform distribution.
- /// * `writeOut` - side-effect function for handling `LogEvent` items.
- /// * `endCondition` - `EndCondition` that dictates when the MH algorithm ends.
- /// * `propose` - proposal `'scale -> 'theta -> 'theta` that generates a jump based on the scale value.
- /// * `tune` - parameter of type `int -> (float * 'a) list -> 'b -> 'b`, where `int` is current iteration,
- /// * `f` - an objective function, `'a -> float`, to optimise.
- /// * `theta1` - initial position in parameter space of type `'a`.
- /// * `l1` - initial value of -log likelihood at theta1 in parameter space
- /// * `d` - history of the chain, of type `(float * 'a) list`. Passing a list here allows continuation of a previous analysis.
- /// * `scale` - a scale of type `'b`, which is compatible with the scale tuning function `tune`
- ///
- /// **Output Type**
- /// * `(float * 'a) list * 'b` - A tuple containing a list of results, and the final scale used in
- /// the analysis. The `(float * 'a) list` represents a list of paired -log likelihood values with
- /// the proposed theta.
- ///
+ /// A recursive metropolis hastings algorithm that ends when `endCondition` returns true.
+ /// `System.Random` to be used for drawing from a uniform distribution.
+ /// side-effect function for handling `LogEvent` items.
+ /// `EndCondition` that dictates when the MH algorithm ends.
+ /// proposal `'scale -> 'theta -> 'theta` that generates a jump based on the scale value.
+ /// parameter of type `int -> (float * 'a) list -> 'b -> 'b`, where `int` is current iteration,
+ /// an objective function, `'a -> float`, to optimise.
+ /// initial position in parameter space of type `'a`.
+ /// initial value of -log likelihood at theta1 in parameter space
+ /// history of the chain, of type `(float * 'a) list`. Passing a list here allows continuation of a previous analysis.
+ /// a scale of type `'b`, which is compatible with the scale tuning function `tune`
+ /// the current iteration number
+ ///
+ /// `(float * 'a) list * 'b` - A tuple containing a list of results, and the final scale used in
+ /// the analysis. The `(float * 'a) list` represents a list of paired -log likelihood values with
+ /// the proposed theta.
let rec metropolisHastings' random writeOut endCondition propose tune f theta1 l1 d scale iteration =
if theta1 |> Array.isEmpty then
invalidOp "Not valid theta"
@@ -601,7 +600,7 @@ module MonteCarlo =
fun random writeOut n domain startPoint f ->
let initialCovariance = TuningMode.covarianceFromBounds 10000 domain random
- match tryGenerateTheta f domain random 10000 with
+ match Initialise.tryGenerateTheta f domain random 10000 with
| Ok theta ->
writeOut <| GeneralEvent(sprintf "[Optimisation] Initial theta is %A" theta)
@@ -637,7 +636,7 @@ module MonteCarlo =
let ``Adaptive-Metropolis-within Gibbs``: Optimiser =
InDetachedSpace
<| fun random writeOut endCon domain startPoint f ->
- match tryGenerateTheta f domain random 10000 with
+ match Initialise.tryGenerateTheta f domain random 10000 with
| Ok theta ->
writeOut <| GeneralEvent(sprintf "[Optimisation] Initial theta is %A" theta)
let sigmas = theta |> Array.map (fun _ -> 0.)
@@ -653,7 +652,7 @@ module MonteCarlo =
let ``Metropolis-within Gibbs``: Optimiser =
InDetachedSpace
<| fun random writeOut endCon domain startPoint (f: Objective) ->
- let theta = initialise domain random
+ let theta = Initialise.initialise domain random
let sigmas = theta |> Array.map (fun _ -> 0.)
let _, result, _ =
@@ -669,7 +668,7 @@ module MonteCarlo =
<| fun random writeOut endCon domain startPoint (f: Objective) ->
// Starting condition
- let initialTheta = initialise domain random
+ let initialTheta = Initialise.initialise domain random
let initialSigma = initialTheta |> Array.map (fun _ -> 0.)
let mwg adapt batchSize currentBatch theta sigmas =
@@ -961,7 +960,7 @@ module MonteCarlo =
// 1. Initial conditions
let draw' = jump random
- let theta1 = initialise domain random
+ let theta1 = Initialise.initialise domain random
let l1 = f theta1
// 2. Chain generator
@@ -1120,6 +1119,12 @@ module MonteCarlo =
MinScaleChange: float
BurnLength: EndCondition<'a> }
+ static member Default =
+ { TuneAfterChanges = 50
+ MaxScaleChange = 100.00
+ MinScaleChange = 0.0010
+ BurnLength = EndConditions.afterIteration 2000 }
+
let filzbach'
settings
(theta: float[])
@@ -1275,7 +1280,7 @@ module MonteCarlo =
filzbach' settings theta random writeOut endCon domain f
| None ->
- match tryGenerateTheta f domain random 10000 with
+ match Initialise.tryGenerateTheta f domain random 10000 with
| Ok theta ->
writeOut <| GeneralEvent(sprintf "[Optimisation] Initial theta is %A" theta)
filzbach' settings theta random writeOut endCon domain f
@@ -1361,28 +1366,30 @@ module Amoeba =
else
shrink a f s
- let initialize (d: Domain) (rng: System.Random) =
- [| for (min, max, _) in d -> min + (max - min) * rng.NextDouble() |]
-
- let solve settings rng writeOut (endWhen: EndCondition) domain startPoint f =
+ let solve settings rng writeOut atEnd domain _ f =
let dim = Array.length domain
let start =
- [| for _ in 1 .. settings.Size -> initialize domain rng |]
+ [| for _ in 1 .. settings.Size -> Initialise.tryGenerateTheta f domain rng 1000 |]
+ |> Array.map (fun r ->
+ match r with
+ | Ok r -> r
+ | Error _ -> failwith "Could not generate theta")
|> Array.map (evaluate f)
|> Array.sortBy fst
let amoeba = { Dim = dim; Solutions = start }
- let rec search i (a: Amoeba) =
- if i > 0 then
- // writeOut <| OptimisationEvent { Iteration = i; Likelihood = a; Theta = thetaAccepted }
- search (i - 1) (update a f settings) // TODO unify array vs list
+ let rec search i trace atEnd (a: Amoeba) =
+ if not <| atEnd trace i then
+ printfn "i %i, val %A" i (trace |> List.tryHead)
+ // writeOut <| OptimisationEvent { Iteration = i; Likelihood = trace; Theta = thetaAccepted }
+ search (i + 1) (a.Best :: trace) atEnd (update a f settings)
else
writeOut <| GeneralEvent(sprintf "Solution: -L = %f" (fst a.Solutions.[0]))
a.Solutions |> Array.toList
- search 50000 amoeba
+ search 0 [] atEnd amoeba
let rec swarm
diff --git a/src/Bristlecone/Parameter.fs b/src/Bristlecone/Parameter.fs
index c005fa4..6c77fba 100755
--- a/src/Bristlecone/Parameter.fs
+++ b/src/Bristlecone/Parameter.fs
@@ -22,8 +22,8 @@ module Parameter =
| PositiveOnly
type Estimation =
- | NotEstimated of lowStartingBound:float * highStartingBound:float
- | Estimated of estimate:float
+ | NotEstimated of lowStartingBound: float * highStartingBound: float
+ | Estimated of estimate: float
type Parameter = private Parameter of Constraint * ConstraintMode * Estimation
@@ -130,7 +130,7 @@ module Parameter =
let c, _, estimate = p |> unwrap
match estimate with
- | NotEstimated (x,y) -> (Parameter(c, Detached, NotEstimated (x,y)), c)
+ | NotEstimated(x, y) -> (Parameter(c, Detached, NotEstimated(x, y)), c)
| Estimated v -> (Parameter(c, Detached, Estimated v), c)
/// Contains the `ParameterPool` type, which represents the set of parameters
@@ -145,7 +145,7 @@ module Parameter =
let toList pool = (pool |> unwrap) |> Map.toList
let hasParameter name pool =
- pool |> unwrap |> Map.tryFindBy(fun k -> k.Value = name)
+ pool |> unwrap |> Map.tryFindBy (fun k -> k.Value = name)
/// Returns Some value if a parameter with the `key`
/// exists in the Pool. The value returned is transformed
@@ -203,12 +203,11 @@ 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)
- | Error _ ->
+ | Error _ ->
printfn "Error in (%A, %A)" c r
None)
|> fromList
diff --git a/src/Bristlecone/Statistics.fs b/src/Bristlecone/Statistics.fs
index 8c64c08..9337a8d 100755
--- a/src/Bristlecone/Statistics.fs
+++ b/src/Bristlecone/Statistics.fs
@@ -1,5 +1,8 @@
namespace Bristlecone.Statistics
+/// Statistics returned during n-step ahead analysis
+type NStepStatistics = { RMSE: float }
+
module Distributions =
open MathNet.Numerics.LinearAlgebra
diff --git a/src/Bristlecone/Test.fs b/src/Bristlecone/Test.fs
index 2e4356a..0194a0b 100644
--- a/src/Bristlecone/Test.fs
+++ b/src/Bristlecone/Test.fs
@@ -222,9 +222,9 @@ module Test =
=
let theta = drawParameterSet engine.Random model.Parameters
- tryGenerateData' engine model settings theta
- |> Result.bind (fun series ->
- let brokeTheRules =
+ match tryGenerateData' engine model settings theta with
+ | Ok series ->
+ let rulesPassed =
settings.GenerationRules
|> List.choose (fun (key, ruleFn) ->
series
@@ -233,12 +233,13 @@ module Test =
Map.find k series |> TimeSeries.toObservations |> Seq.map fst |> ruleFn))
if
- (brokeTheRules |> List.contains false)
- || brokeTheRules.Length <> settings.GenerationRules.Length
+ (rulesPassed |> List.contains false)
+ || rulesPassed.Length <> settings.GenerationRules.Length
then
if attempts = 0 then
Error "Could not generate data that complies with the given ruleset"
else
tryGenerateData engine settings model (attempts - 1)
else
- Ok(series, theta))
+ Ok(series, theta)
+ | Error e -> Error e
diff --git a/src/Bristlecone/Types.fs b/src/Bristlecone/Types.fs
index 84d2160..261a0a4 100755
--- a/src/Bristlecone/Types.fs
+++ b/src/Bristlecone/Types.fs
@@ -304,6 +304,10 @@ module Result =
| Ok x -> Some x
| Error _ -> None
+ let forceOk r =
+ match r with
+ | Ok x -> x
+ | Error e -> failwithf "Error occurred: %s" e
[]
module ListExtensions =
diff --git a/src/Bristlecone/templates/ensemble-convergence.csv b/src/Bristlecone/templates/ensemble-convergence.csv
index e0c56c6..ae35910 100755
--- a/src/Bristlecone/templates/ensemble-convergence.csv
+++ b/src/Bristlecone/templates/ensemble-convergence.csv
@@ -1,2 +1,2 @@
-Subject, Hypothesis, Paramater, Statistic, Value (float)
+Subject, Hypothesis, Parameter, Statistic, Value (float)
YUSL29, some-id, lambda, rHat, 1.245
\ No newline at end of file
diff --git a/src/Bristlecone/templates/ensemble-rmse.csv b/src/Bristlecone/templates/ensemble-rmse.csv
new file mode 100755
index 0000000..0a64011
--- /dev/null
+++ b/src/Bristlecone/templates/ensemble-rmse.csv
@@ -0,0 +1,2 @@
+Subject, Hypothesis, Variable, Statistic, Value (float)
+YUSL29, some-id, x, rmse, 1.245
\ No newline at end of file
diff --git a/src/Bristlecone/templates/individual-n-step-ahead.csv b/src/Bristlecone/templates/individual-n-step-ahead.csv
new file mode 100755
index 0000000..f1b076c
--- /dev/null
+++ b/src/Bristlecone/templates/individual-n-step-ahead.csv
@@ -0,0 +1,2 @@
+Subject, ModelId, AnalysisId, Time (date), Observed (float), StepsAhead (int), NStepAhead (float)
+YUSL21A,some-id,91a95ab8-6c43-45f4-9353-d507cd3a9909, 95, 0.0022, 1, 0.0034
\ No newline at end of file
diff --git a/tests/Bristlecone.Benchmark/Program.fs b/tests/Bristlecone.Benchmark/Program.fs
index f3925f6..77e9c2b 100755
--- a/tests/Bristlecone.Benchmark/Program.fs
+++ b/tests/Bristlecone.Benchmark/Program.fs
@@ -241,7 +241,7 @@ module TimeSeriesTests =
List.allPairs optimFunctions timeModels
|> List.map (fun ((optimName: string, optimise), (modelName, modelFn, startValues)) ->
runReplicated Config.startPointCount (fun () ->
- Bristlecone.Bristlecone.testModel (engine optimise) settings modelFn)
+ Bristlecone.Bristlecone.tryTestModel (engine optimise) settings modelFn)
|> summarise modelName optimName)
@@ -260,16 +260,21 @@ let annealSettings =
let optimFunctions =
[ "amoeba single", Amoeba.single Amoeba.Solver.Default
- // "amoeba swarm", Amoeba.swarm 5 20 Amoeba.Solver.Default
- // "anneal classic", MonteCarlo.SimulatedAnnealing.classicalSimulatedAnnealing 0.01 false annealSettings
- // "anneal cauchy", MonteCarlo.SimulatedAnnealing.fastSimulatedAnnealing 0.01 false annealSettings
- // "filzbach", MonteCarlo.Filzbach.filzbach { TuneAfterChanges = 10000; MaxScaleChange = 0.5; MinScaleChange = 0.5; BurnLength = EndConditions.afterIteration 10000 }
- // "automatic MCMC", MonteCarlo.``Automatic (Adaptive Diagnostics)``
- // "metropolis-gibbs", MonteCarlo.``Metropolis-within Gibbs``
- // "adaptive metropolis", MonteCarlo.adaptiveMetropolis 0.250 500
- // "random walk MCMC", MonteCarlo.randomWalk []
- // "random walk w/ tuning", MonteCarlo.randomWalk [ MonteCarlo.TuneMethod.CovarianceWithScale 0.25, 500, EndConditions.afterIteration 10000 ]
- ]
+ "amoeba swarm", Amoeba.swarm 5 20 Amoeba.Solver.Default
+ "anneal classic", MonteCarlo.SimulatedAnnealing.classicalSimulatedAnnealing 0.01 false annealSettings
+ "anneal cauchy", MonteCarlo.SimulatedAnnealing.fastSimulatedAnnealing 0.01 false annealSettings
+ "filzbach",
+ MonteCarlo.Filzbach.filzbach
+ { TuneAfterChanges = 10000
+ MaxScaleChange = 0.5
+ MinScaleChange = 0.5
+ BurnLength = EndConditions.afterIteration 10000 }
+ "automatic MCMC", MonteCarlo.``Automatic (Adaptive Diagnostics)``
+ "metropolis-gibbs", MonteCarlo.``Metropolis-within Gibbs``
+ "adaptive metropolis", MonteCarlo.adaptiveMetropolis 0.250 500
+ "random walk MCMC", MonteCarlo.randomWalk []
+ "random walk w/ tuning",
+ MonteCarlo.randomWalk [ MonteCarlo.TuneMethod.CovarianceWithScale 0.25, 500, EndConditions.afterIteration 10000 ] ]
let timeModels =
[ "predator-prey (with gaussian noise)",
diff --git a/tests/Bristlecone.Tests/Bristlecone.fs b/tests/Bristlecone.Tests/Bristlecone.fs
index 9623156..8925ec4 100644
--- a/tests/Bristlecone.Tests/Bristlecone.fs
+++ b/tests/Bristlecone.Tests/Bristlecone.fs
@@ -10,8 +10,7 @@ let expectSameFloat a b message =
Expect.isTrue (LanguagePrimitives.GenericEqualityER a b) message
let expectSameFloatList a b message =
- Seq.zip a b
- |> Seq.iter(fun (a,b) -> expectSameFloat a b message)
+ Seq.zip a b |> Seq.iter (fun (a, b) -> expectSameFloat a b message)
@@ -51,47 +50,56 @@ module ``Objective creation`` =
"Objective"
[
- // testProperty "Time-series are paired to correct years" <| fun x ->
- // false
-
- // testProperty "Throws when time-series are not the same length" <| fun x ->
- // false
-
- testPropertyWithConfig Config.config "Likelihood functions use 'real' parameter values"
- <| fun shouldTransform (data: float list) (b1: NormalFloat) (b2: NormalFloat) ->
-
- // Returns the parameter value
- let fakeLikelihood : Bristlecone.ModelSystem.LikelihoodFn =
- fun paramAccessor data ->
- paramAccessor.Get "a"
-
- if b1.Get = b2.Get || b1.Get = 0. || b2.Get = 0.
- then ()
- else
- let b1 = if b1.Get < 0. then b1.Get * -1. else b1.Get
- let b2 = if b2.Get < 0. then b2.Get * -1. else b2.Get
-
- let mode =
- if shouldTransform then Language.notNegative
- else Language.noConstraints
-
- let model =
- Language.Model.empty
- |> Language.Model.addEquation "x" (Language.Parameter "a")
- |> Language.Model.estimateParameter "a" mode (min b1 b2) (max b1 b2)
- |> Language.Model.useLikelihoodFunction fakeLikelihood
- |> Language.Model.compile
+ // testProperty "Time-series are paired to correct years" <| fun x ->
+ // false
- let testObjective =
- Objective.create model (fun _ -> Map.empty) (fun _ _ _ -> [| 2.0 |]) ([ (ShortCode.create "x").Value, data |> List.toArray ] |> Map.ofList)
+ // testProperty "Throws when time-series are not the same length" <| fun x ->
+ // false
- Expect.floatClose
- Accuracy.high
- (testObjective [| (if shouldTransform then Parameter.transformOut mode b1 else b1) |])
- b1
- "The likelihood function did not retrieve the 'real' parameter value"
+ testPropertyWithConfig Config.config "Likelihood functions use 'real' parameter values"
+ <| fun shouldTransform (data: float list) (b1: NormalFloat) (b2: NormalFloat) ->
+
+ // Returns the parameter value
+ let fakeLikelihood: Bristlecone.ModelSystem.LikelihoodFn =
+ fun paramAccessor data -> paramAccessor.Get "a"
+
+ if b1.Get = b2.Get || b1.Get = 0. || b2.Get = 0. then
+ ()
+ else
+ let b1 = if b1.Get < 0. then b1.Get * -1. else b1.Get
+ let b2 = if b2.Get < 0. then b2.Get * -1. else b2.Get
+
+ let mode =
+ if shouldTransform then
+ Language.notNegative
+ else
+ Language.noConstraints
+
+ let model =
+ Language.Model.empty
+ |> Language.Model.addEquation "x" (Language.Parameter "a")
+ |> Language.Model.estimateParameter "a" mode (min b1 b2) (max b1 b2)
+ |> Language.Model.useLikelihoodFunction fakeLikelihood
+ |> Language.Model.compile
+
+ let testObjective =
+ Objective.create
+ model
+ (fun _ -> Map.empty)
+ (fun _ _ _ -> [| 2.0 |])
+ ([ (ShortCode.create "x").Value, data |> List.toArray ] |> Map.ofList)
+
+ Expect.floatClose
+ Accuracy.high
+ (testObjective
+ [| (if shouldTransform then
+ Parameter.transformOut mode b1
+ else
+ b1) |])
+ b1
+ "The likelihood function did not retrieve the 'real' parameter value"
- ]
+ ]
module ``Fit`` =
@@ -101,18 +109,22 @@ module ``Fit`` =
"Conditioning"
[
- testPropertyWithConfig Config.config "Repeating first data point sets t0 as t1"
- <| fun time resolution (data: float list) ->
- if data.IsEmpty || data.Length = 1 then ()
- else
- let data =
- [ (Language.code "x").Value,
- Time.TimeSeries.fromSeq time (Time.FixedTemporalResolution.Years resolution) data ]
- |> Map.ofList
+ testPropertyWithConfig Config.config "Repeating first data point sets t0 as t1"
+ <| fun time resolution (data: float list) ->
+ if data.IsEmpty || data.Length = 1 then
+ ()
+ else
+ let data =
+ [ (Language.code "x").Value,
+ Time.TimeSeries.fromSeq time (Time.FixedTemporalResolution.Years resolution) data ]
+ |> Map.ofList
+
+ let result = Bristlecone.Fit.t0 data Conditioning.RepeatFirstDataPoint ignore
- let result = Bristlecone.Fit.t0 data Conditioning.RepeatFirstDataPoint ignore
- expectSameFloatList
- (result) (data |> Map.map (fun k v -> v.Values |> Seq.head)) "t0 did not equal t1"
+ expectSameFloatList
+ (result)
+ (data |> Map.map (fun k v -> v.Values |> Seq.head))
+ "t0 did not equal t1"
// testProperty "t0 is set as a custom point when specified" <| fun () ->
// false
@@ -130,23 +142,29 @@ module ``Fit`` =
testPropertyWithConfig Config.config "Core fitting functions are reproducible"
<| fun b1 b2 seedNumber (obs: float list) startDate months ->
- if System.Double.IsNaN b1 || b1 = infinity || b1 = -infinity ||
- System.Double.IsNaN b2 || b2 = infinity || b2 = -infinity
- then ()
+ if
+ System.Double.IsNaN b1
+ || b1 = infinity
+ || b1 = -infinity
+ || System.Double.IsNaN b2
+ || b2 = infinity
+ || b2 = -infinity
+ then
+ ()
else
let data: CodedMap> =
[ (Language.code "x").Value,
- Time.TimeSeries.fromSeq startDate (Time.FixedTemporalResolution.Months months) obs ]
+ Time.TimeSeries.fromSeq startDate (Time.FixedTemporalResolution.Months months) obs ]
|> Map.ofList
let result =
Expect.wantOk
- (Bristlecone.fit defaultEngine defaultEndCon data (TestModels.constant b1 b2))
+ (Bristlecone.tryFit defaultEngine defaultEndCon data (TestModels.constant b1 b2))
"Fitting did not happen successfully."
let result2 =
Expect.wantOk
- (Bristlecone.fit
+ (Bristlecone.tryFit
{ defaultEngine with
Random = MathNet.Numerics.Random.MersenneTwister(seedNumber, true) }
defaultEndCon
@@ -155,9 +173,21 @@ module ``Fit`` =
""
expectSameFloat result.Likelihood result2.Likelihood "Different likelihoods"
- expectSameFloat result.InternalDynamics result.InternalDynamics "Different internal dynamics"
+
+ expectSameFloat
+ result.InternalDynamics
+ result.InternalDynamics
+ "Different internal dynamics"
+
expectSameFloat result.Parameters result2.Parameters "Different parameters"
- expectSameFloatList (result.Series |> Seq.collect(fun kv -> kv.Value.Values |> Seq.map(fun v -> v.Fit))) (result2.Series |> Seq.collect(fun kv -> kv.Value.Values |> Seq.map(fun v -> v.Fit))) "Different expected series"
+
+ expectSameFloatList
+ (result.Series
+ |> Seq.collect (fun kv -> kv.Value.Values |> Seq.map (fun v -> v.Fit)))
+ (result2.Series
+ |> Seq.collect (fun kv -> kv.Value.Values |> Seq.map (fun v -> v.Fit)))
+ "Different expected series"
+
expectSameFloat result.Trace result2.Trace "Different traces"
// testProperty "Time-series relating to model equations must overlap"
@@ -188,13 +218,18 @@ module ``Fit`` =
"Setting up parameter constraints"
[
- testPropertyWithConfig Config.config "Positive only parameter is transformed when optimising in transformed space"
+ testPropertyWithConfig
+ Config.config
+ "Positive only parameter is transformed when optimising in transformed space"
<| fun (data: float list) startDate months (b1: NormalFloat) (b2: NormalFloat) ->
- let testModel b1 b2 = TestModels.twoEquationConstant Language.notNegative b1 b2
- if b1.Get = b2.Get || b1.Get = 0. || b2.Get = 0.
- then
- Expect.throws (fun () -> testModel b1.Get b2.Get |> ignore) "Model compiled despite having no difference between parameter bounds"
- else
+ let testModel b1 b2 =
+ TestModels.twoEquationConstant Language.notNegative b1 b2
+
+ if b1.Get = b2.Get || b1.Get = 0. || b2.Get = 0. then
+ Expect.throws
+ (fun () -> testModel b1.Get b2.Get |> ignore)
+ "Model compiled despite having no difference between parameter bounds"
+ else
let b1 = if b1.Get < 0. then b1.Get * -1. else b1.Get
let b2 = if b2.Get < 0. then b2.Get * -1. else b2.Get
let mutable inOptimMin = nan
@@ -210,19 +245,21 @@ module ``Fit`` =
{ defaultEngine with
OptimiseWith = optimTest }
- let data =
+ let data =
[ (ShortCode.create "x").Value; (ShortCode.create "y").Value ]
- |> Seq.map(fun c -> c, Time.TimeSeries.fromSeq startDate (Time.FixedTemporalResolution.Months months) data)
+ |> Seq.map (fun c ->
+ c,
+ Time.TimeSeries.fromSeq startDate (Time.FixedTemporalResolution.Months months) data)
|> Map.ofSeq
let result =
Expect.wantOk
- (Bristlecone.fit engine defaultEndCon data (testModel b1 b2))
+ (Bristlecone.tryFit engine defaultEndCon data (testModel b1 b2))
"Errored when should be OK"
Expect.equal
inOptimMin
- (min (log(b1)) (log(b2)))
+ (min (log (b1)) (log (b2)))
"The lower bound was not transformed inside the optimiser" ]
]
diff --git a/tests/Bristlecone.Tests/Config.fs b/tests/Bristlecone.Tests/Config.fs
index abc9306..2073381 100644
--- a/tests/Bristlecone.Tests/Config.fs
+++ b/tests/Bristlecone.Tests/Config.fs
@@ -53,7 +53,7 @@ type BristleconeTypesGen() =
}
|> Arb.fromGen
- static member CodedMap<'snd> () =
+ static member CodedMap<'snd>() =
gen {
let! n = Gen.choose (1, 100)
let! codes = Gen.listOfLength n Arb.generate
diff --git a/tests/Bristlecone.Tests/Language.fs b/tests/Bristlecone.Tests/Language.fs
index d956689..8fda765 100644
--- a/tests/Bristlecone.Tests/Language.fs
+++ b/tests/Bristlecone.Tests/Language.fs
@@ -55,19 +55,22 @@ let modelExpressions =
<| fun (c: NormalFloat) x t pool env -> Constant c.Get |> compute x t pool env = c.Get
testPropertyWithConfig Config.config "Getting parameter value fails when parameter not present"
- <| fun (code:ShortCode.ShortCode) x t pool e ->
- let f () = Parameter code.Value |> compute x t pool e
+ <| fun (code: ShortCode.ShortCode) x t pool e ->
+ let f () =
+ Parameter code.Value |> compute x t pool e
- match pool |> Parameter.Pool.toList |> List.tryFind (fun (k, _) -> k.Value = code.Value) with
+ match
+ pool
+ |> Parameter.Pool.toList
+ |> List.tryFind (fun (k, _) -> k.Value = code.Value)
+ with
| Some p ->
- if p |> snd |> Parameter.isEstimated
- then
- Expect.equal
- (f ())
- (p |> snd |> Parameter.getTransformedValue)
- "Did not fail when parameter was not present"
- | None ->
- Expect.throws (fun () -> f () |> ignore) "Parameter was not present"
+ if p |> snd |> Parameter.isEstimated then
+ Expect.equal
+ (f ())
+ (p |> snd |> Parameter.getTransformedValue)
+ "Did not fail when parameter was not present"
+ | None -> Expect.throws (fun () -> f () |> ignore) "Parameter was not present"
testPropertyWithConfig Config.config "Getting parameter values returns real value when present"
<| fun pool x t e ->
@@ -76,20 +79,27 @@ let modelExpressions =
|> Gen.sample 1 1
|> List.head
- if not (pool |> Parameter.Pool.hasParameter selectedCode.Value |> Option.map Parameter.isEstimated).Value
- then ()
+ if
+ not
+ (pool
+ |> Parameter.Pool.hasParameter selectedCode.Value
+ |> Option.map Parameter.isEstimated)
+ .Value
+ then
+ ()
else
- let existingValue =
- pool |> Parameter.Pool.tryGetRealValue selectedCode.Value |> Option.get
+ let existingValue =
+ pool |> Parameter.Pool.tryGetRealValue selectedCode.Value |> Option.get
- let result = Parameter selectedCode.Value |> compute x t pool e
+ let result = Parameter selectedCode.Value |> compute x t pool e
- Expect.equal result existingValue "The parameter value was not correct"
+ Expect.equal result existingValue "The parameter value was not correct"
testPropertyWithConfig Config.config "Fails when environmental (aka time-varying) data is not present"
- <| fun (code:ShortCode.ShortCode) x t pool e ->
- let f () = Environment code.Value |> compute x t pool e
+ <| fun (code: ShortCode.ShortCode) x t pool e ->
+ let f () =
+ Environment code.Value |> compute x t pool e
match e |> Map.tryFindBy (fun m -> m.Value = code.Value) with
| Some environ -> Expect.equal (f ()) environ "Did not fail when environment data was not present"
@@ -111,7 +121,7 @@ let modelBuilder =
let f () =
likelihoodFns
|> Seq.fold (fun mb l -> mb |> Model.useLikelihoodFunction l) Model.empty
- |> Model.addEquation "x" (Constant 1.)
+ |> Model.addEquation "x" (Constant 1.)
|> Model.compile
if likelihoodFns |> Seq.length <> 1 then
@@ -137,61 +147,73 @@ let modelBuilder =
if eqs |> Seq.length <> 1 then
Expect.throws (fun () -> fn () |> ignore) "Did not throw when no equations specified"
- testPropertyWithConfig Config.config "Compiles with one likelihood function and one or more equations (no duplicate keys)"
+ testPropertyWithConfig
+ Config.config
+ "Compiles with one likelihood function and one or more equations (no duplicate keys)"
<| fun l (eqs: ShortCode.ShortCode list) ->
- if eqs |> Seq.hasDuplicates then ()
+ if eqs |> Seq.hasDuplicates then
+ ()
else
- let mb =
- eqs
- |> List.fold
- (fun mb k -> mb |> Model.addEquation k.Value (Constant 1.))
- (Model.empty |> Model.useLikelihoodFunction l)
-
- if eqs |> Seq.isEmpty then
- Expect.throws (fun () -> mb |> Model.compile |> ignore) "Did not error when no equations specified"
- else
- mb |> Model.compile |> ignore
+ let mb =
+ eqs
+ |> List.fold
+ (fun mb k -> mb |> Model.addEquation k.Value (Constant 1.))
+ (Model.empty |> Model.useLikelihoodFunction l)
+
+ if eqs |> Seq.isEmpty then
+ Expect.throws
+ (fun () -> mb |> Model.compile |> ignore)
+ "Did not error when no equations specified"
+ else
+ mb |> Model.compile |> ignore
testPropertyWithConfig Config.config "Compiles whether measures are present or not"
- <| fun likelihood (measures:CodedMap) ->
- if measures.Keys |> Seq.hasDuplicates then ()
+ <| fun likelihood (measures: CodedMap) ->
+ if measures.Keys |> Seq.hasDuplicates then
+ ()
else
- let model =
- Model.empty
- |> Model.useLikelihoodFunction likelihood
- |> Model.addEquation "eq1" (Constant 1.)
+ let model =
+ Model.empty
+ |> Model.useLikelihoodFunction likelihood
+ |> Model.addEquation "eq1" (Constant 1.)
- measures
- |> Map.fold (fun mb n m -> mb |> Model.includeMeasure n.Value m) model
- |> Model.compile
- |> ignore
+ measures
+ |> Map.fold (fun mb n m -> mb |> Model.includeMeasure n.Value m) model
+ |> Model.compile
+ |> ignore
testPropertyWithConfig Config.config "Doesn't compile if duplicate keys exist"
- <| fun likelihood (eqs: ShortCode.ShortCode list) (measures: (ShortCode.ShortCode * ModelSystem.MeasureEquation) list) ->
- if eqs.IsEmpty then ()
- else
- let compile () =
- eqs
- |> Seq.fold
- (fun mb n -> mb |> Model.addEquation n.Value This)
- (Model.empty |> Model.useLikelihoodFunction likelihood)
- |> fun mb ->
- Seq.fold
- (fun mb (n: ShortCode.ShortCode, eq) -> mb |> Model.includeMeasure n.Value eq)
- mb measures
- |> Model.compile
-
- let keys = [ eqs; (measures |> List.map fst) ] |> List.concat
- if keys |> Seq.hasDuplicates then
- Expect.throws (compile >> ignore) "Duplicate keys existed"
- else
- compile () |> ignore
-
- // testProperty "Only compiles when all required parameters are specified" <| fun (pool:Parameter.Pool) ->
-
- // testProperty "Only compiles when all specified parameters are used" <| fail
-
- // testProperty "Equations in the built model have the correct result" <| fail
+ <| fun
+ likelihood
+ (eqs: ShortCode.ShortCode list)
+ (measures: (ShortCode.ShortCode * ModelSystem.MeasureEquation) list) ->
+ if eqs.IsEmpty then
+ ()
+ else
+ let compile () =
+ eqs
+ |> Seq.fold
+ (fun mb n -> mb |> Model.addEquation n.Value This)
+ (Model.empty |> Model.useLikelihoodFunction likelihood)
+ |> fun mb ->
+ Seq.fold
+ (fun mb (n: ShortCode.ShortCode, eq) -> mb |> Model.includeMeasure n.Value eq)
+ mb
+ measures
+ |> Model.compile
+
+ let keys = [ eqs; (measures |> List.map fst) ] |> List.concat
+
+ if keys |> Seq.hasDuplicates then
+ Expect.throws (compile >> ignore) "Duplicate keys existed"
+ else
+ compile () |> ignore
+
+ // testProperty "Only compiles when all required parameters are specified" <| fun (pool:Parameter.Pool) ->
+
+ // testProperty "Only compiles when all specified parameters are used" <| fail
+
+ // testProperty "Equations in the built model have the correct result" <| fail
]
diff --git a/tests/Bristlecone.Tests/Test.fs b/tests/Bristlecone.Tests/Test.fs
index 23ac1d0..18f2218 100644
--- a/tests/Bristlecone.Tests/Test.fs
+++ b/tests/Bristlecone.Tests/Test.fs
@@ -3,4 +3,4 @@ module TestTests
// module ``When testing model-fitting`` =
// []
-// let
\ No newline at end of file
+// let
diff --git a/tests/Bristlecone.Tests/Workflow.fs b/tests/Bristlecone.Tests/Workflow.fs
index 063d453..0efa2dc 100755
--- a/tests/Bristlecone.Tests/Workflow.fs
+++ b/tests/Bristlecone.Tests/Workflow.fs
@@ -21,4 +21,4 @@ module Orchestration =
// failwith "Not implemented" ]
- ]
\ No newline at end of file
+ ]