prepping the documenation and cleaning up the repo

docs/make.jl

@@ -1,22 +1,27 @@
-push!(LOAD_PATH, joinpath(@__DIR__, "../src/"))
 using Documenter, WENDy
-# CI_FLG = get(ENV, "CI", nothing) == "true"
+
+# DocMeta.setdocmeta(WENDy, :DocTestSetup, :(using WENDy); recursive=true)
 
 makedocs(
     modules = [WENDy],
+    checkdocs=:exports,
+    # doctest=true, 
+    # linktest=true, 
+    # authors="Nicholas Rummel <[email protected]>",
+    # repo="https://github.com/nrummel/WENDy.jl/blob/{commit}{path}#{line}",
+    sitename="WENDy.jl",
     format = Documenter.HTML(
-        # prettyurls = CI_FLG,
-        canonical = "http://weavejl.mpastell.com/stable/",
+        prettyurls =  get(ENV, "CI", nothing) == "true",
+        canonical = "https://nrummel.github.io/WENDy.jl/",
     ),
-    sitename = "WENDy.jl",
     pages = [
-        "index.md",
-        "getting_started.md",
-        "usage.md",
+        "Home"=>"index.md",
+        "Examples"=>"examples.md",
+        "Reference"=>"reference.md",
     ],
 )
 
 deploydocs(
-    repo = "github.com/JunoLab/Weave.jl.git",
+    repo = "github.com/nrummel/WENDy.jl",
     push_preview = true,
 )

docs/src/index.md

@@ -1,36 +1,42 @@
 # WENDy.jl - Weak Form Estimation of Nonlinear Dynamics
 
-This is the documentation of [WENDy].jl](https://github.com/nrummel/WENDy.jl).
+This is the documentation of [WENDy.jl](https://github.com/nrummel/WENDy.jl). The work comes from the the MathBio Group at University of Colorado Boulder. For further reading find our paper at [arxiv link](https://arxiv.org/).
+
+## Problem Statement
 WENDy is an algorithm that can estimate unknown parameters for ordinary differential equations given noisey data.
 
 The set up for this algorithm is to assume that a physical system with state variable, $\boldsymbol{u} \in \mathbb{R}^D$, is governed by a system of ordinary differential equation with true parameters, $\mathbf{p}^* \in \mathbb{R}^J$:
 ```math
     \dot{\boldsymbol{u}}(t) = f(\boldsymbol{u}(t), t, \mathbf{p}^*)
 ```
-The user has observed data of this system. The data has been corrupted by that measurement error:
-
+The user has observed data of this system on a uniform grid, $$\{t_m \mathbf{u}_m\}_{m=0}^M$$. The data has been corrupted by noise:
+- **Additive Gaussian Case:**
 ```math
     \begin{align*}
-    \text{\textbf{Additive Gaussian Case: } }\\ 
-    \{\mathbf{u}_m &= \boldsymbol{u}(t_m, p^*) + \epsilon_m \}_{m=0}^M \\
-    \epsilon &\sim \mathcal{N}(\mathbf{0}, \mathbb{I}_D)\\
-
-    \text{\textbf{Multiplicative LogNormal Case: } }\\ 
-    \{\mathbf{u}_m &= \boldsymbol{u}(t_m, p^*) \eta_m \}_{m=0}^M \\
-    \log(\eta) &\sim \mathcal{N}(\mathbf{0}, \mathbb{I}_D)\\
+        \{\mathbf{u}_m &= \boldsymbol{u}(t_m, p^*) + \epsilon_m \}_{m=0}^M \\
+        \epsilon_m &\stackrel{iid}{\sim} \mathcal{N}(\mathbf{0}, \mathbb{I}_D)\\
     \end{align*}
 ```
-The goal of the algorithm is that recover parameters $\hat{\mathbf{p}}$. In other words, we hope that if the one were to solve the system of differential equations with the estimated parameters then it would match the true state, then 
+- **Multiplicative LogNormal Case:**
 ```math
-    \frac{\| \boldsymbol{u}(t; \hat{\mathbf{p}}) - \boldsymbol{u}(t; \mathbf{p}^*)\|}{\|\boldsymbol{u}(t; \mathbf{p}^*)\|} \ll 1 \\
+    \begin{align*}
+        \{\mathbf{u}_m &= \boldsymbol{u}(t_m, p^*) \circ \eta_m \}_{m=0}^M \\
+        \log(\eta) &\stackrel{iid}{\sim} \mathcal{N}(\mathbf{0}, \mathbb{I}_D)\\
+    \end{align*}
+```
+*Note*: The Hammard product $$\circ$$ is the element-wise multiplication on the two vectors. 
+
+The goal of the algorithm is that recover unknown parameters $\mathbf{p}$. In other words, we hope that if the one were to solve the system of differential equations with the estimated parameters then it would match the true state, then 
+```math
+    \frac{\| \boldsymbol{u}(t; \mathbf{p}) - \boldsymbol{u}(t; \mathbf{p}^*)\|}{\|\boldsymbol{u}(t; \mathbf{p}^*)\|} \ll 1 \\
 ```
 This is done by leveraging a estimated distribution of the weak form residual, $\mathbf{r}$, and then approximating a maximum likelihood estimate: 
 ```math
     \mathbf{S}(\mathbf{p})^{-\tfrac{1}{2}} \mathbf{r}(\mathbf{p}) \sim \mathcal{N}(0, \mathbb{I})
 ```
-The work comes from the the MathBio Group at University of Colorado Boulder. For further reading find our paper at [arxiv link](https://arxiv.org/).
 
-**Current features**
+
+## Current features
 
 - Estimation of parameters for ordinary differential equations
 - Supports: 
@@ -40,12 +46,4 @@ The work comes from the the MathBio Group at University of Colorado Boulder. For
     - Box constraints for parameter spaces
 - Provides acceleration for problems that are linear in parameters 
 - Directly calls robust optimization algorithms that are well suited to an non-convex problems.
-- Creates efficient Julia functions for the likelihood function and its derivatives with minimal inputs from the end user.
-
-## Index
-
-```@contents
-Pages = [
-    "index.md",
-    "getting_started.md"
-]
+- Creates efficient Julia functions for the likelihood function and its derivatives with minimal inputs from the end user.

examples/exampleGoodwin.jl

@@ -4,10 +4,10 @@ using WENDy
 using OrdinaryDiffEq: ODEProblem
 using OrdinaryDiffEq: solve as solve_ode
 ## Define rhs, ic, time domain, and length of parameters
-function f!(du, u, w, t)
-    du[1] = w[1] / (w[3] + u[3]^w[4]) - w[2] * u[1]
-    du[2] = w[5]*u[1]- w[6]*u[2]
-    du[3] = w[7]*u[2]-w[8]*u[3]
+function f!(du, u, p, t)
+    du[1] = p[1] / (p[3] + u[3]^p[4]) - p[2] * u[1]
+    du[2] = p[5]*u[1]- p[6]*u[2]
+    du[3] = p[7]*u[2]-p[8]*u[3]
     nothing
 end
 tRng  = (0.0, 80.0)
@@ -26,29 +26,32 @@ ode = ODEProblem(
 tt    = tRng[1]:dt:tRng[end]
 Mp1   = length(tt)
 Ustar = reduce(vcat, um' for um in solve_ode(ode, saveat=dt).u)
-snr   = 0.01
+snr   = 0.1
 U     = Ustar.* exp.(snr*randn(size(Ustar)))
 ## Create wendy problem struct
-wendyProb = WENDyProblem(tt, U, f!, J, Val(false)#=Nonlinear in parameters=#, Val(LogNormal)#=LogNormalNoise=#; ll=Logging.Info);
+wendyProb = WENDyProblem(tt, U, f!, J, Val(false)#=Nonlinear in parameters=#, Val(LogNormal)#=LogNormalNoise=#, 
+    WENDyParameters(radiusMinTime=dt, radiusMaxTime=tRng[end]/5); # be sure to set the min and max testFunction radii to something reasonable
+    ll=Logging.Info # turn on logging information
+);
 ## Solve the wendy problm given an intial guess for the parameters 
 J = length(wstar)
-w0 = wstar + 0.5*randn(J).*abs.(wstar)
-# w0 = [0.7336635446929285, 0.34924148955718615, 13.976533673829369, 6.282647403425017, 0.9149290036628314, 0.19280811624587546, 0.7701803478856664, 0.7805192636751863]
-relErr = norm(w0 - wstar) / norm(wstar)
+p₀ = wstar + 0.5*randn(J).*abs.(wstar)
+relErr = norm(p₀ - wstar) / norm(wstar)
 @info "Initializing with Relative Coefficient Error = $(relErr)"
 @info "Solving wendy problem ..."
-@time what = WENDy.solve(wendyProb, w0)
+@time what = WENDy.solve(wendyProb, p₀)
 relErr = norm(what - wstar) / norm(wstar)
 @info "Relative Coefficient Error = $(relErr)"
 ## plot the resutls 
 odeprob = ODEProblem(f!, u₀, tRng, what)
 sol = solve_ode(odeprob; saveat=dt)
 Uhat = reduce(vcat, um' for um in sol.u)
+colors = ["red", "blue", "green"]
 plot(
     reduce(vcat, [
-        scatter(x=tt, y=Uhat[:,d], name="estimate"),
-        scatter(x=tt, y=U[:,d], name="data", mode="markers"),
-        scatter(x=tt, y=Ustar[:,d], name="truth")
+        scatter(x=tt, y=Uhat[:,d],line_color=colors[d], line_dash="dash", name="estimate", legendgroup=d, legendgrouptitle_text="u[$d]"),
+        scatter(x=tt, y=U[:,d], marker_color=colors[d], name="data", mode="markers", legendgroup=d ),
+        scatter(x=tt, y=Ustar[:,d],line_color=colors[d], name="truth", legendgroup=d )
     ] for d in 1:D),
-    Layout(title="Results",xaxis_title="time(s)", yaxis_title="State")
+    Layout(title="Goodwin Example",xaxis_title="time(s)", yaxis_title="State")
 )

examples/exampleLogistic.jl

@@ -4,8 +4,8 @@ using WENDy
 using OrdinaryDiffEq: ODEProblem
 using OrdinaryDiffEq: solve as solve_ode
 ## Define rhs, ic, time domain, and length of parameters
-function f!(du, u, w, t)
-    du[1] = w[1] * u[1] - w[2] * u[1]^2
+function f!(du, u, p, t)
+    du[1] = p[1] * u[1] - p[2] * u[1]^2
     nothing
 end
 tRng = (0.0, 10.0)
@@ -26,10 +26,10 @@ U = Ustar + 0.01*randn(size(Ustar))
 ## Create wendy problem struct
 wendyProb = WENDyProblem(tt, U, f!, J, Val(true), Val(Normal); ll=Logging.Info)
 ## Solve the wendy problm given an intial guess for the parameters 
-w0 = [0.5, 0.5]
-relErr = norm(w0 - wstar) / norm(wstar)
+p₀ = [0.5, 0.5]
+relErr = norm(p₀ - wstar) / norm(wstar)
 @info "Init Relative Coefficient Error = $(relErr)"
-what = WENDy.solve(wendyProb, w0)
+what = WENDy.solve(wendyProb, p₀)
 relErr = norm(what - wstar) / norm(wstar)
 @info "Relative Coefficient Error = $(relErr)"
 ## plot the resutls 

examples/exampleSIRTDI.jl

@@ -4,11 +4,11 @@ using WENDy
 using OrdinaryDiffEq: ODEProblem
 using OrdinaryDiffEq: solve as solve_ode
 ## Define rhs, ic, time domain, and length of parameters
-function f!(du, u, w, t)
-    β = (w[1] * exp(-w[1] * w[2])) / (1 - exp(-w[1] * w[2]))
-    du[1] = -w[1] * u[1] + w[3] * u[2] + β * u[3] 
-    du[2] = w[1] * u[1] - w[3] * u[2] - w[4] * (1 - exp(-w[5]  * t^2)) * u[2] 
-    du[3] = w[4] * (1 - exp(-w[5]  * t^2)) * u[2] - β * u[3]
+function f!(du, u, p, t)
+    β = (p[1] * exp(-p[1] * p[2])) / (1 - exp(-p[1] * p[2]))
+    du[1] = -p[1] * u[1] + p[3] * u[2] + β * u[3] 
+    du[2] = p[1] * u[1] - p[3] * u[2] - p[4] * (1 - exp(-p[5]  * t^2)) * u[2] 
+    du[3] = p[4] * (1 - exp(-p[5]  * t^2)) * u[2] - β * u[3]
 end
 tRng  = (0.0, 50.0)
 dt    = 0.1
@@ -45,11 +45,11 @@ params = WENDyParameters(
 wendyProb = WENDyProblem(tt, U, f!, J, Val(false)#=Nonlinear in parameters=#, Val(LogNormal)#=LogNormalNoise=#, params; ll=Logging.Info, constraints=constraints);
 ## Solve the wendy problm given an intial guess for the parameters 
 J = length(wstar)
-w0 = wstar + 0.5*randn(J).*abs.(wstar)
-relErr = norm(w0 - wstar) / norm(wstar)
+p₀ = wstar + 0.5*randn(J).*abs.(wstar)
+relErr = norm(p₀ - wstar) / norm(wstar)
 @info "Initializing with Relative Coefficient Error = $(relErr)"
 @info "Solving wendy problem ..."
-@time what = WENDy.solve(wendyProb, w0)
+@time what = WENDy.solve(wendyProb, p₀)
 relErr = norm(what - wstar) / norm(wstar)
 @info "Relative Coefficient Error = $(relErr)"
 ## plot the resutls 

src/wendyDataTypes.jl

@@ -1,8 +1,31 @@
 """
-    WENDyParameters
-
 Hyper-parameters for the WENDy Algorithm
 
+# Constructor 
+    WENDyParameters(;   
+        diagReg::Real=1.0e-10,
+        radiusMinTime::Real=0.01,
+        radiusMaxTime::Real=5.0,
+        numRadii::Int=100,
+        radiiParams::AbstractVector{<:Real}=2 .^(0:3),
+        testFunSubRate::Real=2.0,
+        maxTestFunCondNum::Real=1e4,
+        minTestFunInfoNum::Real=0.95,
+        Kmax::Int=200,
+        Kᵣ::Union{Nothing,Int}=100,
+        fsAbstol::Real=1e-8,
+        fsReltol::Real=1e-8,
+        nlsAbstol::Real=1e-8,
+        nlsReltol::Real=1e-8,
+        nlsMaxiters::Int=1000,
+        optimAbstol::Real=1e-8,
+        optimReltol::Real=1e-8,
+        optimMaxiters::Int=500,
+        optimTimelimit::Real=200.0,
+        fsAlg::OrdinaryDiffEqAlgorithm=Rodas4P(),
+        fsU0Free::Bool=true
+    )
+
 # Fields
 - diagReg::Real = 1.0e-10 : Regularization constant for the covariance computations
 - radiusMinTime::Real = 0.01 : Minimum test function radius (in time units matching _tt)
@@ -16,13 +39,13 @@ Hyper-parameters for the WENDy Algorithm
 - Kᵣ::Union{Nothing,Int} = 100 : how many test function to spread through the time domain in the min radii detection sub-algorithm
 - fsAbstol::Real = 1e-8 : forward solve absolute tolerance for solving ordinary differential equation
 - fsReltol::Real = 1e-8 : forward solve relative tolerance for solving ordinary differntial equation
-- nlsAbstol::Real = 1e-8 : nonlinear least squares absolute tolerance (only used in the IRWLS WENDy algorithm)
-- nlsReltol::Real = 1e-8 : nonlinear least squares relative tolerance (only used in the IRWLS WENDy algorithm)
-- nlsMaxiters::Int = 1000 : nonlinear least squares maximum iterations (only used in the IRWLS WENDy algorithm)
+- nlsAbstol::Real = 1e-8 : nonlinear least squares absolute tolerance (only used in the IRLS WENDy algorithm)
+- nlsReltol::Real = 1e-8 : nonlinear least squares relative tolerance (only used in the IRLS WENDy algorithm)
+- nlsMaxiters::Int = 1000 : nonlinear least squares maximum iterations (only used in the IRLS WENDy algorithm)
 - optimAbstol::Real = 1e-8 : absolute tolerance (used by all other optimization algorithms)
 - optimReltol::Real = 1e-8 : relative tolerance (used by all other optimization algorithms)
 - optimMaxiters::Int = 200 : maximum iterations (used by all other optimization algorithms)
-- optimTimelimit::Real = 200.0 : maximum time in seconds (used by all other optimization algorithms)
+- optimTimelimit::Real = 500.0 : maximum time in seconds (used by all other optimization algorithms)
 - fsAlg::OrdinaryDiffEqAlgorithm = Rodas4P() :  forward solve algorithm used by the forward solve nonlinear least squares algorithm
 - fsU0Free::Bool = true : Specifies if the forward solve algorithm should also optimize over the initial condition
 """
@@ -68,16 +91,16 @@ struct WENDyInternals{lip, DistType}
     Hₚf!::Function
     Hₚ∇ₓf!::Function
 end
-## IRWLS 
-abstract type IRWLS_Iter end 
+## IRLS 
+abstract type IRLSIter end 
 ##
-struct Linear_IRWLS_Iter <: IRWLS_Iter
+struct LinearIRLSIter <: IRLSIter
     b₀::AbstractVector{<:AbstractFloat}
     G0::AbstractMatrix{<:AbstractFloat}
     Rᵀ!::Function 
 end 
 ##
-struct NLS_iter <: IRWLS_Iter
+struct NonLinearIRLSIter <: IRLSIter
     b₀::AbstractVector
     Rᵀ!::Function 
     r!::Function
@@ -106,10 +129,33 @@ struct LeastSquaresCostFunction <: CostFunction
     KD::Int
 end 
 """
-    WENDyProblem{lip, DistType}
+    WENDyProblem{lip, DistType}(...)
 
 A WENDyProblem struct pre-computes and allocates data structures for efficient solving of the parameter inverse problem
 
+# Constructor 
+    WENDyProblem(
+        _tt::AbstractVector{<:Real}, 
+        U::AbstractVecOrMat{<:Real}, 
+        _f!::Function, 
+        J::Int, 
+        ::Val{lip}=Val(false), 
+        ::Val{DistType}=Val(Normal), 
+        params::WENDyParameters=WENDyParameters(); 
+        constraints::Union{Nothing,AbstractVector{Tuple{<:Real,<:Real}}}=nothing, 
+        ll::LogLevel=Warn
+    )
+## Arguments 
+- _tt::AbstractVector{<:Real} : vector of times (equispaced)
+- U::AbstractVecOrMat{<:Real} : Corrupted state variable data 
+- _f!::Function : Right hand-side of the differential equation
+- J::Int : number of parameters (to be estimated)
+- ::Val{lip}=Val(false) : (optional) specify whether the right hand side is `linear in parameters' for improved computational efficiency
+- ::Val{DistType}=Val(Normal) : (optional) specify the distribution of the measurement noise. Choose either Val(Normal) for additive Gaussian noise of Val(LogNormal) for multiplicative LogNormal noise.
+- params::WENDyParameters : (optional) struct of hyper-parameters for the WENDy Algorithm (see the doc for WENDyParameters)
+- constraints=nothing : (optional) Linear box constraints for each parameter, ∀j ∈ [1, ⋯,J], ℓⱼ ≤ pⱼ ≤ uⱼ. Accepts constraints as a list of tuples, [(ℓ₁,u₁), ⋯]. Note: this only is compatible with the TrustRegion solver.
+- ll::LogLevel=Warn : (optional) see additional algorithm information by setting ll=Info
+
 # Fields
 - D::Int : number of state variables
 - J::Int : number of parameters (to be estimated)
@@ -118,23 +164,10 @@ A WENDyProblem struct pre-computes and allocates data structures for efficient s
 - u₀::AbstractVector{<:Real} : Initial Condition of the ODE (Necessary for the forward solver)
 - constraints : vector of tuples containing linear constraints for each parameter
 - data : Internal data structure 
-- fslsq::LeastSquaresCostFunction : Cost function for the comparison method
+- oels::LeastSquaresCostFunction : Cost function for the comparison method
 - wlsq::LeastSquaresCostFunction : Cost function for the weak form least squares problem
 - wnll::SecondOrderCostFunction : Cost function for the weak form negative log-likelihood 
 
-# Constructor 
-    WENDyProblem(_tt::AbstractVector{<:Real}, U::AbstractVecOrMat{<:Real}, _f!::Function, J::Int, ::Val{lip}=Val(false), ::Val{DistType}=Val(Normal), params::WENDyParameters=WENDyParameters(); constraints::Union{Nothing,AbstractVector{Tuple{T1,T2}}}=nothing, ll::LogLevel=Warn)
-It is recommend that you use the built-in constructor to build your WENDyProblem
-## Arguments 
-- _tt::AbstractVector{<:Real} : vector of times (equispaced)
-- U::AbstractVecOrMat{<:Real} : Corrupted state variable data 
-- _f!::Function : Right hand-side of the differential equation
-- J::Int : number of parameters (to be estimated)
-- ::Val{lip}=Val(false) : (optional) specify whether the right hand side is `linear in parameters' for improved computational efficiency
-- ::Val{DistType}=Val(Normal) : (optional) specify the distribution of the measurement noise
-- params::WENDyParameters : (optional) struct of hyper-parameters for the WENDy Algorithm (see the doc for WENDyParameters)
-- constraints=nothing : Linear box constraints for each parameter
-- ll::LogLevel=Warn : (optional) see additional algorithm information by setting ll=Info
 """
 struct WENDyProblem{lip, DistType}
     D::Int # number of state variables
@@ -145,7 +178,7 @@ struct WENDyProblem{lip, DistType}
     constraints::Union{Nothing,AbstractVector{Tuple{<:Real,<:Real}}}
     data::WENDyInternals{lip, DistType}
     # Cost functions 
-    fslsq::LeastSquaresCostFunction
+    oels::LeastSquaresCostFunction
     wlsq::LeastSquaresCostFunction
     wnll::SecondOrderCostFunction 
 end