Remove inplace functions

Author: briochemc

Project.toml

@@ -1,7 +1,7 @@
 name = "AIBECS"
 uuid = "ace601d6-714c-11e9-04e5-89b7fad23838"
 authors = ["Benoit Pasquier <[email protected]>"]
-version = "0.12.0"
+version = "0.13.0"
 
 [deps]
 Bijectors = "76274a88-744f-5084-9051-94815aaf08c4"

src/multiTracer.jl

@@ -37,42 +37,6 @@ function AIBECSFunction(T::Function, Gs::Tuple, nb::Int)
     AIBECSFunction((T,), Gs, nb, nt, Tidx)
 end
 function AIBECSFunction(Ts::Tuple, Gs::Tuple, nb::Int, nt::Int=length(Ts), Tidx::AbstractVector=1:nt)
-    if all(isinplace(G, nt + 2) for G in Gs) # +2 to account for dx and p
-        iipAIBECSFunction(Ts, Gs, nb, nt, Tidx)
-    else
-        oopAIBECSFunction(Ts, Gs, nb, nt, Tidx)
-    end
-end
-
-function iipAIBECSFunction(Ts, Gs, nb, nt=length(Ts), Tidx=1:nt)
-    tracers(u) = state_to_tracers(u, nb, nt)
-    tracer(u, i) = state_to_tracer(u, nb, nt, i)
-    function G(du, u, p)
-        for j in 1:nt
-            Gs[j](tracer(du, j), tracers(u)..., p)
-        end
-        du
-    end
-    function f(du, u, p, t=0)
-        G(du, u, p)
-        for jT in eachindex(Ts)
-            op = Ts[jT](p)
-            for j in findall(Tidx .== jT)
-                mul!(tracer(du, j), op, tracer(u, j), -1, 1)
-            end
-        end
-        du
-    end
-    # Jacobian
-    ∇ₓG(u, p) = inplace_local_jacobian(Gs, u, p, nt, nb)
-    function T(p)
-        uniqueTs = [Tⱼ(p) for Tⱼ in Ts]
-        blockdiag([uniqueTs[Tidx[j]] for j in 1:nt]...)
-    end
-    jac(u, p, t=0) = ∇ₓG(u, p) - T(p)
-    return ODEFunction{true}(f, jac=jac)
-end
-function oopAIBECSFunction(Ts, Gs, nb, nt=length(Ts), Tidx=1:nt)
     tracers(u) = state_to_tracers(u, nb, nt)
     tracer(u, i) = state_to_tracer(u, nb, nt, i)
     G(u, p) = reduce(vcat, Gⱼ(tracers(u)..., p) for Gⱼ in Gs)
@@ -99,54 +63,16 @@ end
 function AIBECSFunction(Ts, Gs, nb::Int, ::Type{P}) where {P <: APar}
     AIBECSFunction(AIBECSFunction(Ts, Gs, nb), P)
 end
-function AIBECSFunction(fun::ODEFunction{false}, ::Type{P}) where {P <: APar}
+function AIBECSFunction(fun::ODEFunction, ::Type{P}) where {P <: APar}
     jac(u, p::P, t=0) = fun.jac(u, p, t)
     jac(u, λ::Vector, t=0) = fun.jac(u, λ2p(P, λ), t)
     f(u, p::P, t=0) = fun.f(u, p, t)
     f(u, λ::Vector, t=0) = fun.f(u, λ2p(P, λ), t)
     return ODEFunction{false}(f, jac=jac)
 end
-function AIBECSFunction(fun::ODEFunction{true}, ::Type{P}) where {P <: APar}
-    jac(u, p::P, t=0) = fun.jac(u, p, t)
-    jac(u, λ::Vector, t=0) = fun.jac(u, λ2p(P, λ), t)
-    f(du, u, p::P, t=0) = fun.f(du, u, p, t)
-    f(du, u, λ::Vector, t=0) = fun.f(du, u, λ2p(P, λ), t)
-    return ODEFunction{true}(f, jac=jac)
-end
 
 export AIBECSFunction
 
-"""
-    F, ∇ₓF = F_and_∇ₓF(Ts, Gs, nb)
-
-Returns the state function `F` and its Jacobian, `∇ₓF`.
-
-    F, ∇ₓF = F_and_∇ₓF(T, Gs, nb)
-
-Returns the state function `F` and its Jacobian, `∇ₓF` (with all tracers transported by single `T`).
-
-This function is deprecated. Use
-
-    F = AIBECSFunction(Ts, Gs, nb)
-
-instead.
-You can then call `F(x,p)` for the tendencies, and `F(Val{:jac},x,p)` for the Jacobian.
-"""
-function F_and_∇ₓF(fun::ODEFunction)
-    Base.depwarn("""Deprecation:
-        F, ∇ₓF = F_and_∇ₓF(args...)
-    is deprecated. Use
-        F = AIBECSFunction(args...)
-    instead! (And then you can still call `F(x,p)`.)
-    """, :F_and_∇ₓF, force=true)
-    fun.f, fun.jac
-end
-F_and_∇ₓF(args...) = F_and_∇ₓF(AIBECSFunction(args...))
-export F_and_∇ₓF
-
-
-
-
 """
     localderivative(G, x, p)
     localderivative(Gᵢ, xs, i, p)
@@ -166,9 +92,6 @@ end
 function localderivative(Gᵢ, xs, j, p) # for multiple tracers
     return ForwardDiff.derivative(λ -> Gᵢ(perturb_tracer(xs, j, λ)..., p), 0.0)
 end
-function localderivative(Gᵢ!, dx, xs, j, p) # if Gᵢ are in-place
-    return ForwardDiff.derivative((dx, λ) -> Gᵢ!(dx, perturb_tracer(xs, j, λ)..., p), dx, 0.0)
-end
 perturb_tracer(xs, j, λ) = (xs[1:j - 1]..., xs[j] .+ λ, xs[j + 1:end]...)
 
 
@@ -202,19 +125,11 @@ export split_state_function_and_Jacobian
 function local_jacobian(Gs, x, p, nt, nb)
     return reduce(vcat, local_jacobian_row(Gⱼ, x, p, nt, nb) for Gⱼ in Gs)
 end
-function inplace_local_jacobian(Gs, x, p, nt, nb)
-    return reduce(vcat, inplace_local_jacobian_row(Gⱼ!, x, p, nt, nb) for Gⱼ! in Gs)
-end
 
 function local_jacobian_row(Gᵢ, x, p, nt, nb)
     tracers(x) = state_to_tracers(x, nb, nt)
     return reduce(hcat, sparse(Diagonal(localderivative(Gᵢ, tracers(x), j, p))) for j in 1:nt)
 end
-function inplace_local_jacobian_row(Gᵢ!, x, p, nt, nb)
-    tracers(x) = state_to_tracers(x, nb, nt)
-    dx = Vector{Float64}(undef, nb)
-    return reduce(hcat, sparse(Diagonal(localderivative(Gᵢ!, dx, tracers(x), j, p))) for j in 1:nt)
-end
 
 #= ============================================
 Generate 𝑓 and derivatives from user input

test/bgc_functions.jl

@@ -49,32 +49,10 @@ function sms_POP(DIP, DOP, POP, p)
 end
 sms_all = (sms_DIP, sms_DOP, sms_POP) # bundles all the source-sink functions in a tuple
 
-#===========================================
-In-place sources minus sinks
-===========================================#
-
-function G_DIP!(dx, DIP, DOP, POP, p)
-    @unpack τgeo, xgeo, τDIP, k, τDOP = p
-    @. dx = (xgeo - DIP) / τgeo - (DIP ≥ 0) / τDIP * DIP^2 / (DIP + k) * (ztop == 0) + DOP / τDOP
-end
-function G_DOP!(dx, DIP, DOP, POP, p)
-    @unpack τgeo, xgeo, τDIP, k, τDOP, τPOP, σ = p
-    @. dx = σ * (DIP ≥ 0) / τDIP * DIP^2 / (DIP + k) * (ztop == 0) - DOP / τDOP + POP / τPOP
-end
-function G_POP!(dx, DIP, DOP, POP, p)
-    @unpack τgeo, xgeo, τDIP, k, τDOP, τPOP, σ = p
-    @. dx = (1 - σ) * (DIP ≥ 0) / τDIP * DIP^2 / (DIP + k) * (ztop == 0) - POP / τPOP
-end
-Gs = (G_DIP!, G_DOP!, G_POP!)
-
-
 #===========================================
 AIBECS F and ∇ₓF
 ===========================================#
 fun = AIBECSFunction(T_all, sms_all, nb)
-fun! = AIBECSFunction((T_D, T_POP), Gs, nb, 3, [1, 1, 2])
-F, ∇ₓF = F_and_∇ₓF(T_all, sms_all, nb) # deprecated
-F!, ∇ₓF! = F_and_∇ₓF((T_D, T_POP), Gs, nb, 3, [1, 1, 2]) # deprecated
 
 #===========================================
 AIBECS split operators for CNLF
@@ -126,19 +104,6 @@ Tests
     end
 end
 
-
-@testset "In-place functions" begin
-    nt = length(T_all)
-    n = nt * nb
-    @unpack xgeo = p
-    x = xgeo * ones(n) .* exp.(cos.(collect(1:n)) / 10)
-    dx = similar(x)
-    @testset "in-place F ≈ out-of-place F" begin
-        @test fun(x, p) ≈ fun!(dx, x, p) rtol = 1e-10
-        @test fun.jac(x, p) ≈ fun!.jac(x, p) rtol = 1e-10
-    end
-end
-
 @testset "Linear/nonlinear split functions" begin
     nt = length(T_all)
     n = nt * nb

test/parameters.jl

@@ -27,7 +27,7 @@ function prior(::Type{T}, s::Symbol) where {T<:AbstractParameters}
             σ = 10.0 # Assumes that a sensible unit is chosen (i.e., that within 10.0 * U)
             return Normal(μ, σ)
         else
-            return LocationScale(lb, ub-lb, LogitNormal()) # <- The LogitNormal works well for Optim?
+            return lb + (ub-lb) * LogitNormal() # <- The LogitNormal works well for Optim?
         end
     else
         return nothing

test/solvers.jl

@@ -16,13 +16,6 @@ algs = [CTKAlg]
             @test s isa SteadyStateSolution
             @test oopprob isa SteadyStateProblem
             @test norm(s.u) / norm(fun.f(s.u, testp, 0)) > ustrip(upreferred(1e5u"Myr"))
-
-            iipprob = SteadyStateProblem(fun!, nx * x, testp)
-            s = solve(iipprob, alg())
-            @test s isa SteadyStateSolution
-            @test iipprob isa SteadyStateProblem
-            ds = copy(s.u)
-            @test norm(s.u) / norm(fun!.f(ds, s.u, testp, 0)) > ustrip(upreferred(1e5u"Myr"))
         end
     end
 end