Fix for new AIBECS version (#10)

* Fix for new AIBECS version

* Fix email

* Fix ReadMe

* Deploy docs in PRs

Author: briochemc

.github/workflows/CompatHelper.yml

@@ -0,0 +1,16 @@
+name: CompatHelper
+on:
+  schedule:
+    - cron: '00 00 * * *'
+  workflow_dispatch:
+jobs:
+  CompatHelper:
+    runs-on: ubuntu-latest
+    steps:
+      - name: Pkg.add("CompatHelper")
+        run: julia -e 'using Pkg; Pkg.add("CompatHelper")'
+      - name: CompatHelper.main()
+        env:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+          COMPATHELPER_PRIV: ${{ secrets.DOCUMENTER_KEY }}
+        run: julia -e 'using CompatHelper; CompatHelper.main()'

Project.toml

@@ -1,7 +1,7 @@
 name = "F1Method"
 uuid = "d5605bae-6d76-11e9-2fd7-71bcf42edbaa"
-authors = ["Benoit Pasquier <pasquieb@uci.edu>"]
-version = "0.3.2"
+authors = ["Benoit Pasquier <briochemc@gmail.com>"]
+version = "0.4.0"
 
 [deps]
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
@@ -11,7 +11,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 [compat]
 DiffEqBase = "6"
 ForwardDiff = "0.10"
-julia = "1.0"
+julia = "1"
 
 [extras]
 AIBECS = "ace601d6-714c-11e9-04e5-89b7fad23838"

README.md

@@ -5,7 +5,7 @@
 
 <p>
   <a href="https://briochemc.github.io/F1Method.jl/stable/">
-    <img src=https://img.shields.io/badge/docs-stable-important.svg?style=flat-square&label=Documentation&logo=Read%20the%20Docs>
+    <img src="https://img.shields.io/github/workflow/status/briochemc/F1Method.jl/Documentation?style=for-the-badge&label=Documentation&logo=Read%20the%20Docs&logoColor=white">
   </a>
 </p>
 
@@ -19,8 +19,14 @@
 </p>
 
 <p>
-  <a href="https://travis-ci.com/briochemc/F1Method.jl">
-    <img alt="Build Status" src="https://img.shields.io/travis/com/briochemc/F1Method.jl/master?label=OSX/Linux/Windows&logo=travis&logocolor=white&style=flat-square">
+  <a href="https://github.com/briochemc/F1Method.jl/actions">
+    <img src="https://img.shields.io/github/workflow/status/briochemc/F1Method.jl/Mac%20OS%20X?label=OSX&logo=Apple&logoColor=white&style=flat-square">
+  </a>
+  <a href="https://github.com/briochemc/F1Method.jl/actions">
+    <img src="https://img.shields.io/github/workflow/status/briochemc/F1Method.jl/Linux?label=Linux&logo=Linux&logoColor=white&style=flat-square">
+  </a>
+  <a href="https://github.com/briochemc/F1Method.jl/actions">
+    <img src="https://img.shields.io/github/workflow/status/briochemc/F1Method.jl/Windows?label=Windows&logo=Windows&logoColor=white&style=flat-square">
   </a>
   <a href="https://codecov.io/gh/briochemc/F1Method.jl">
     <img src="https://img.shields.io/codecov/c/github/briochemc/F1Method.jl/master?label=Codecov&logo=codecov&logoColor=white&style=flat-square">
@@ -75,17 +81,16 @@ Once initial values for the state, `x`, and parameters, `p`, are chosen, simply
 
 ```julia
 # Initialize the cache for storing reusable objects
-mem = F1Method.initialize_mem(x, p)
+mem = initialize_mem(F, ∇ₓf, ∇ₓF, x, p, alg; options...)
 ```
 
 wrap the functions into functions of `p` only via
 
 ```julia
-
 # Wrap the objective, gradient, and Hessian functions
-objective(p) = F1Method.objective(f, F, ∇ₓF, mem, p, myAlg(); my_options...)
-gradient(p) = F1Method.gradient(f, F, ∇ₓf, ∇ₓF, mem, p, myAlg(); my_options...)
-hessian(p) = F1Method.hessian(f, F, ∇ₓf, ∇ₓF, mem, p, myAlg(); my_options...)
+objective(p) = F1Method.objective(f, F, ∇ₓF, mem, p, alg; options...)
+gradient(p) = F1Method.gradient(f, F, ∇ₓf, ∇ₓF, mem, p, alg; options...)
+hessian(p) = F1Method.hessian(f, F, ∇ₓf, ∇ₓF, mem, p, alg; options...)
 ```
 
 and compute the objective, gradient, or Hessian via either of
@@ -101,10 +106,9 @@ hessian(p)
 That's it.
 You were told it was simple, weren't you?
 Now you can test how fast and accurate it is!
-(Or trust our published work, *[Pasquier and Primeau](https://www.bpasquier.com/publication/pasquier_primeau_sisc_2019/)* (in preparation).)
 
 ## Citing the software
 
 If you use this package, or implement your own package based on the F-1 algorithm please cite us.
-If you use the F-1 algorithm, please cite *[Pasquier and Primeau](https://www.bpasquier.com/publication/pasquier_primeau_sisc_2019/)* (in review for publication in the *SIAM Journal on Scientific Computing*).
-If you also use this package directly, please cite us using the [CITATION.bib](./CITATION.bib), which contains a bibtex entry for the software.
+If you use the F-1 algorithm, please cite *[Pasquier and Primeau](https://www.bpasquier.com/publication/pasquier_primeau_sisc_2019/)* (in prep.).
+If you also use this package directly, please cite it! (Use [the Zenodo link](https://doi.org/10.5281/zenodo.2667835) or the [CITATION.bib file](./CITATION.bib), which contains a bibtex entry.)

docs/Project.toml

@@ -4,4 +4,4 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 
 [compat]
-Documenter = "~0.24"
+Documenter = "~0.25"

docs/make.jl

@@ -9,4 +9,5 @@ makedocs(
 
 deploydocs(
     repo = "github.com/briochemc/F1Method.jl.git",
+    push_preview = true
 )

docs/src/index.md

@@ -123,7 +123,7 @@ Finally, we wrap the objective, gradient, and Hessian functions defined by the F
 
 ```jldoctest usage
 # Initialize the cache for storing reusable objects
-mem = F1Method.initialize_mem(x₀, p₀)
+mem = F1Method.initialize_mem(F, ∇ₓf, ∇ₓF, x₀, p₀, MyAlg())
 # Define the functions via the F1 method
 F1_objective(p) = F1Method.objective(f, F, ∇ₓF, mem, p, MyAlg())
 F1_gradient(p) = F1Method.gradient(f, F, ∇ₓf, ∇ₓF, mem, p, MyAlg())

src/F1Method.jl

@@ -21,32 +21,32 @@ Contains
 - `p`   the parameters 𝒑
 The `Mem`-type object should be initialized with `initialize_mem`.
 """
-mutable struct Mem
-    s     # 𝒔(𝒑)
-    A     # factors of 𝐀 = ∇ₓ𝑭(𝒔,𝒑)
-    ∇s    # ∇𝒔(𝒑)
-    ∇ₓf   # ∇ₓ𝑓(𝒔,𝒑)
-    p     # 𝒑
+mutable struct Mem{Ts, TA, T∇s, T∇ₓf, Tp}
+    s::Ts     # 𝒔(𝒑)
+    A::TA     # factors of 𝐀 = ∇ₓ𝑭(𝒔,𝒑)
+    ∇s::T∇s   # ∇𝒔(𝒑)
+    ∇ₓf::T∇ₓf # ∇ₓ𝑓(𝒔,𝒑)
+    p::Tp     # 𝒑 that matches all but 𝒔(𝒑)
+    psol::Tp  # 𝒑 that matches 𝒔(𝒑)
 end
 
 
 function update_mem!(f, F, ∇ₓf, ∇ₓF, mem, p, alg; options...)
     if p ≠ mem.p                      # only update mem if 𝒑 has changed
         update_solution!(F, ∇ₓF, mem, p, alg; options...)
-        s, m = mem.s.u, length(p)
-        ∇ₚF = ForwardDiff.jacobian(λ -> F(s,p+λ), zeros(m))
-        mem.A = factorize(∇ₓF(s,p))   # update factors of ∇ₓ𝑭(𝒔,𝒑)
-        mem.∇s .= mem.A \ -∇ₚF        # update ∇𝒔
-        mem.∇ₓf .= ∇ₓf(s,p)           # update ∇ₓ𝑓(𝒔,𝒑)
-        mem.p = p                     # update 𝒑
+        ∇ₚF = ForwardDiff.jacobian(p -> F(mem.s, p), p)
+        mem.A = factorize(∇ₓF(mem.s, p)) # update factors of ∇ₓ𝑭(𝒔,𝒑)
+        mem.∇s .= mem.A \ -∇ₚF           # update ∇𝒔
+        mem.∇ₓf .= ∇ₓf(mem.s, p)         # update ∇ₓ𝑓(𝒔,𝒑)
+        mem.p .= p                  # update 𝒑 for the variables above
     end
 end
 
 function update_solution!(F, ∇ₓF, mem, p, alg; options...)
-    if ~(mem.s isa SteadyStateSolution) || p ≠ mem.s.prob.p
-        mem.s isa SteadyStateSolution ? x = mem.s.u : x = mem.s
-        prob = SteadyStateProblem(F, ∇ₓF, x, p)       # define problem
-        mem.s = solve(prob, alg; options...)          # update 𝒔
+    if p ≠ mem.psol
+        prob = SteadyStateProblem(F, ∇ₓF, mem.s, p) # define problem
+        mem.s .= solve(prob, alg; options...).u      # update 𝒔
+        mem.psol .= p                          # update 𝒑 for 𝒔
     end
 end
 
@@ -63,7 +63,7 @@ The Jacobian, `∇ₓF`, and the memory cache `mem` must be supplied.
 """
 function objective(f, F, ∇ₓF, mem, p, alg; options...)
     update_solution!(F, ∇ₓF, mem, p, alg; options...)
-    return f(mem.s,p)
+    return f(mem.s, p)
 end
 
 """
@@ -74,7 +74,7 @@ Returns the gradient of the `objective` function using the F-1 method.
 function gradient(f, F, ∇ₓf, ∇ₓF, mem, p, alg; options...)
     update_mem!(f, F, ∇ₓf, ∇ₓF, mem, p, alg; options...)
     s, ∇s, m = mem.s, mem.∇s, length(p)
-    ∇ₚf = ForwardDiff.jacobian(λ -> [f(s,p+λ)], zeros(m))
+    ∇ₚf = ForwardDiff.jacobian(p -> [f(s,p)], p)
     return mem.∇ₓf * ∇s + ∇ₚf
 end
 
@@ -96,9 +96,15 @@ end
 
 Initializes the memory cache for the F-1 method.
 """
-function initialize_mem(x, p)
-    n, m = length(x), length(p)
-    return Mem(copy(x), nothing, zeros(n,m), zeros(1,n), nothing)
+function initialize_mem(F, ∇ₓf, ∇ₓF, x, p, alg; options...)
+    x = copy(x)
+    p = copy(p)
+    psol = copy(p)
+    prob = SteadyStateProblem(F, ∇ₓF, x, p)
+    s = solve(prob, alg; options...).u
+    A = factorize(∇ₓF(s, p))
+    ∇ₚF = ForwardDiff.jacobian(p -> F(s, p), p)
+    return Mem(s, A, A \ -∇ₚF, ∇ₓf(s, p), p, psol)
 end
 
 end

test/AIBECS_test.jl

@@ -1,17 +1,30 @@
 
 # This AIBECS test is derived from the P-model tutorial in AIBECS.jl
 # But it uses one of the samll circulations to test derivatives of the cost function
+module SubModuleAIBECS
 
+using Test
+using F1Method
+using AIBECS
+import AIBECS: @units, units
+import AIBECS: @limits, limits
+import AIBECS: @initial_value, initial_value
+import AIBECS: @flattenable, flattenable, flatten
+import AIBECS: @prior, prior
+using Unitful: m, d, s, yr, Myr, mol, mmol, μmol, μM
+using Distributions
+using WorldOceanAtlasTools
+using FiniteDiff
 
 # AIBECS model
 grd, T_Circ = Primeau_2x2x2.load()
 T_DIP(p) = T_Circ
 T_POP(p) = transportoperator(grd, z -> w(z,p))
 function w(z,p)
     @unpack w₀, w′ = p
-    return @. w₀ + w′ * z
+    w₀ + w′ * z
 end
-const z_top = topdepthvec(grd) # uptake only in top layer
+z_top = topdepthvec(grd) # uptake only in top layer
 function U(x,p)
     @unpack σ, τ_DIP, k = p
     return @. σ * x/τ_DIP * x/(x+k) * (z_top==0) * (x≥0)
@@ -28,11 +41,8 @@ function G_POP(DIP, POP, p)
     @unpack τ_geo = p
     return @. $U(DIP,p) - $R(POP,p) - POP / τ_geo
 end
-import AIBECS: @units, units
-import AIBECS: @limits, limits
-import AIBECS: @initial_value, initial_value
-import AIBECS: @flattenable, flattenable, flatten
-const ∞ = Inf
+
+∞ = Inf
 @initial_value @units @flattenable @limits struct PmodelParameters{U} <: AbstractParameters{U}
     w₀::U       |  0.64 | m/d      | true  | (0,∞)
     w′::U       |  0.13 | m/d/m    | true  | (0,∞)
@@ -43,7 +53,6 @@ const ∞ = Inf
     DIP_geo::U  |  2.12 | mmol/m^3 | true  | (-∞,∞)
     σ::U        |  0.3  | NoUnits  | true  | (0,1)
 end
-import AIBECS: @prior, prior
 function prior(::Type{T}, s::Symbol) where {T<:AbstractParameters}
     if flattenable(T, s)
         lb, ub = limits(T, s)
@@ -65,12 +74,13 @@ prior(::T, s::Symbol) where {T<:AbstractParameters} = prior(T,s)
 prior(::Type{T}) where {T<:AbstractParameters} = Tuple(prior(T,s) for s in AIBECS.symbols(T))
 prior(::T) where {T<:AbstractParameters} = prior(T)
 p = PmodelParameters()
+λ = p2λ(p)
 nb = sum(iswet(grd))
-𝐹, ∇ₓ𝐹 = state_function_and_Jacobian((T_DIP, T_POP), (G_DIP, G_POP), nb)
+F, ∇ₓF = F_and_∇ₓF((T_DIP, T_POP), (G_DIP, G_POP), nb, PmodelParameters)
 @unpack DIP_geo = p
 x = DIP_geo * ones(2nb) # initial guess
-prob = SteadyStateProblem(𝐹, ∇ₓ𝐹, x, p)
-sol = solve(prob, CTKAlg()).u
+prob = SteadyStateProblem(F, ∇ₓF, x, p)
+sol = solve(prob, CTKAlg())
 
 # Observations from World Ocean Atlas used to evaluate the
 # AIBECS model mismatch with observations
@@ -81,42 +91,23 @@ modify(DIP, POP) = (DIP,)
 ωs = (1.0,) # the weight for the mismatch (weight of POP = 0)
 ωp = 1.0       # the weight for the parameters prior estimates
 obs = (DIPobs,)
-𝑓, ∇ₓ𝑓, ∇ₚ𝑓 = generate_objective_and_derivatives(ωs, ωp, grd, modify, obs)
+f, ∇ₓf = f_and_∇ₓf(ωs, ωp, grd, modify, obs, PmodelParameters)
 
 # Now we apply the F1 method
-mem = F1Method.initialize_mem(x, p)
-objective(p) = F1Method.objective(𝑓, 𝐹, ∇ₓ𝐹, mem, p, CTKAlg(), τstop=ustrip(u"s", 1e3u"Myr"))
-gradient(p) = F1Method.gradient(𝑓, 𝐹, ∇ₓ𝑓, ∇ₓ𝐹, mem, p, CTKAlg(), τstop=ustrip(u"s", 1e3u"Myr"))
-hessian(p) = F1Method.hessian(𝑓, 𝐹, ∇ₓ𝑓, ∇ₓ𝐹, mem, p, CTKAlg(), τstop=ustrip(u"s", 1e3u"Myr"))
-
-# and convert p::PmodelParameters to λ::Vector according to AIBECS change of variables
-λ2p = subfun(typeof(p))
-∇λ2p = ∇subfun(typeof(p))
-∇²λ2p = ∇²subfun(typeof(p))
-p2λ = invsubfun(typeof(p))
-λ = p2λ(p)
-function obj(λ)
-    return objective(λ2p(λ))
-end
-function grad(λ)
-    return gradient(λ2p(λ)) * Diagonal(vec(∇λ2p(λ)))
-end
-function hess(λ)
-    p = λ2p(λ)
-    ∇p = Diagonal(vec(∇λ2p(λ)))
-    ∇²p = ∇²λ2p(λ)
-    G = gradient(p)
-    H = hessian(p)
-    return ∇p * H * ∇p + Diagonal(vec(G)) * ∇²p
-end
+τ = ustrip(u"s", 1e3u"Myr")
+mem = F1Method.initialize_mem(F, ∇ₓf, ∇ₓF, x, λ, CTKAlg(), τstop=τ)
+objective(λ) = F1Method.objective(f, F, ∇ₓF, mem, λ, CTKAlg(), τstop=τ)
+gradient(λ) = F1Method.gradient(f, F, ∇ₓf, ∇ₓF, mem, λ, CTKAlg(), τstop=τ)
+hessian(λ) = F1Method.hessian(f, F, ∇ₓf, ∇ₓF, mem, λ, CTKAlg(), τstop=τ)
 
 # Finally we test the result with the "reliable" FiniteDiff :)
-λ = p2λ(p)
-@test FiniteDiff.finite_difference_gradient(obj, λ) ≈ grad(λ)' rtol=1e-5
-@test FiniteDiff.finite_difference_hessian(obj, λ) ≈ hess(λ) rtol=1e-3
-
+@test FiniteDiff.finite_difference_gradient(objective, 2λ) ≈ gradient(2λ)' rtol=1e-3
+@test FiniteDiff.finite_difference_hessian(objective, 2λ) ≈ hessian(2λ) rtol=1e-3
+@test FiniteDiff.finite_difference_gradient(objective, λ) ≈ gradient(λ)' rtol=1e-3
+@test FiniteDiff.finite_difference_hessian(objective, λ) ≈ hessian(λ) rtol=1e-3
 
 
 
 
+end #submodule