change definition of ∂

Author: rveltz

src/Utils.jl

@@ -22,7 +22,7 @@ end
 ####################################################################################################
 # iterated derivatives
 ∂(f) = x -> ForwardDiff.derivative(f, x)
-∂(f, n::Int) = n == 0 ? f : ∂(∂(f), n-1)
+∂(f, ::Val{n}) where {n} = n == 0 ? f : ∂(∂(f), Val(n-1))
 ####################################################################################################
 function print_nonlinear_step(step, residual, itlinear = 0, lastRow = false)
     if lastRow

src/periodicorbit/PeriodicOrbitCollocation.jl

@@ -1156,7 +1156,7 @@ function compute_error!(coll::PeriodicOrbitOCollProblem, x::AbstractVector{𝒯}
     # sol is the piecewise polynomial approximation of y.
     # However, sol is of degree m, hence ∂(sol, m+1) = 0
     # we thus estimate yᵐ⁺¹ using ∂(sol, m)
-    dmsol = ∂(sol, m)
+    dmsol = ∂(sol, Val(m))
     # we find the values of vm := ∂m(x) at the mid points
     τsT = getmesh(coll) .* period
     vm = [ dmsol( (τsT[i] + τsT[i+1]) / 2 ) for i = 1:Ntst ]

src/periodicorbit/PeriodicOrbitTrapeze.jl

@@ -304,6 +304,7 @@ function potrap_functional_jac!(pb::AbstractPOFDProblem, out, u, par, du)
 end
 
 (pb::PeriodicOrbitTrapProblem)(u::AbstractVector, par) = potrap_functional!(pb, similar(u), u, par)
+residual(pb::PeriodicOrbitTrapProblem, u::AbstractVector, par) = pb(u, par)
 (pb::PeriodicOrbitTrapProblem)(u::AbstractVector, par, du) = potrap_functional_jac!(pb, similar(du), u, par, du)
 
 ####################################################################################################
@@ -436,7 +437,7 @@ function (pb::PeriodicOrbitTrapProblem)(::Val{:JacFullSparse}, u0::AbstractVecto
     AγBlock = jacobian_potrap_block(pb, u0, par; γ = γ)
 
     # we now set up the last line / column
-    @views ∂TGpo = (pb(vcat(u0[begin:end-1], T + δ), par) .- pb(u0, par)) ./ δ
+    @views ∂TGpo = (residual(pb, vcat(u0[begin:end-1], T + δ), par) .- residual(pb, u0, par)) ./ δ
 
     # this is "bad" for performance. Get converted to SparseMatrix at the next line
     Aγ = block_to_sparse(AγBlock) # most of the computing time is here!!
@@ -493,7 +494,7 @@ This method returns the jacobian of the functional G encoded in PeriodicOrbitTra
             # J0[(M-1)*N+1:(M)*N, (M-1)*N+1:(M)*N] .= In
 
         # we now set up the last line / column
-        ∂TGpo = (pb(vcat(u0[1:end-1], T + δ), par) .- pb(u0, par)) ./ δ
+        ∂TGpo = (residual(pb,vcat(u0[1:end-1], T + δ), par) .- residual(pb,u0, par)) ./ δ
         J0[:, end] .=  ∂TGpo
 
         # this following does not depend on u0, so it does not change. However we update it in case the caller updated the section somewhere else
@@ -547,7 +548,7 @@ end
 
     if updateborder
         # we now set up the last line / column
-        ∂TGpo = (pb(vcat(u0[1:end-1], T + δ), par) .- pb(u0, par)) ./ δ
+        ∂TGpo = (residual(pb, vcat(u0[1:end-1], T + δ), par) .- residual(pb, u0, par)) ./ δ
         J0[:, end] .=  ∂TGpo
 
         # this following does not depend on u0, so it does not change. However we update it in case the caller updated the section somewhere else
@@ -776,7 +777,7 @@ function (J::POTrapJacobianBordered)(u0::AbstractVector, par; δ = convert(eltyp
     T = extract_period_fdtrap(J.Aγ.prob, u0)
     # we compute the derivative of the problem w.r.t. the period TODO: remove this or improve!!
     # TODO REMOVE CE vcat!
-    @views J.∂TGpo .= (J.Aγ.prob(vcat(u0[begin:end-1], T + δ), par) .- J.Aγ.prob(u0, par)) ./ δ
+    @views J.∂TGpo .= (residual(J.Aγ.prob, vcat(u0[begin:end-1], T + δ), par) .- residual(J.Aγ.prob, u0, par)) ./ δ
 
     J.Aγ(u0, par) # update Aγ
 
@@ -848,10 +849,10 @@ function _newton_trap(probPO::PeriodicOrbitTrapProblem,
             _J =  probPO(Val(:JacFullSparse), orbitguess, getparams(probPO.prob_vf)) |> Array
             jac = (x, p) -> probPO(Val(:JacFullSparseInplace), _J, x, p)
         elseif jacobianPO == :DenseAD
-            jac = (x, p) -> ForwardDiff.jacobian(z -> probPO(z, p), x)
+            jac = (x, p) -> ForwardDiff.jacobian(z -> residual(probPO, z, p), x)
         elseif jacobianPO == :FullMatrixFreeAD
-            jac = (x, p) -> dx -> ForwardDiff.derivative(t -> probPO(x .+ t .* dx, p), 0)
-        else
+            jac = (x, p) -> dx -> ForwardDiff.derivative(t -> residual(probPO, x .+ t .* dx, p), 0)
+        else 
              jac = (x, p) -> ( dx -> probPO(x, p, dx))
         end
 
@@ -984,9 +985,9 @@ function continuation_potrap(prob::PeriodicOrbitTrapProblem,
             _J =  prob(Val(:JacFullSparse), orbitguess, getparams(prob.prob_vf)) |> Array
             jac = (x, p) -> (prob(Val(:JacFullSparseInplace), _J, x, p); FloquetWrapper(prob, _J, x, p));
         elseif jacobianPO == :DenseAD
-            jac = (x, p) -> FloquetWrapper(prob, ForwardDiff.jacobian(z -> prob(z, p), x), x, p)
+            jac = (x, p) -> FloquetWrapper(prob, ForwardDiff.jacobian(z -> residual(prob, z, p), x), x, p)
         elseif jacobianPO == :FullMatrixFreeAD
-            jac = (x, p) -> FloquetWrapper(prob, dx -> ForwardDiff.derivative(t->prob(x .+ t .* dx, p), 0), x, p)
+            jac = (x, p) -> FloquetWrapper(prob, dx -> ForwardDiff.derivative(t->residual(prob, x .+ t .* dx, p), 0), x, p)
         else
              jac = (x, p) -> FloquetWrapper(prob, x, p)
         end

src/periodicorbit/PeriodicOrbits.jl

@@ -596,7 +596,7 @@ function continuation(br::AbstractResult{PeriodicOrbitCont, Tprob},
     # perform continuation
     pbnew = set_params_po(pbnew, setparam(br, newp))
 
-    pbnew(orbitguess, setparam(br, newp))[end] |> abs > 1 && @warn "PO constraint not satisfied"
+    residual(pbnew, orbitguess, setparam(br, newp))[end] |> abs > 1 && @warn "PO constraint not satisfied"
 
     branch = continuation( pbnew, orbitguess, alg, _contParams;
         kwargs..., # put this first to be overwritten just below!

test/stuartLandauCollocation.jl

@@ -51,7 +51,7 @@ _orbit(t) = [cos(2pi * t), 0, 0] * sqrt(par_sl.r / par_sl.c3)
 _ci = BK.generate_solution(prob_col, _orbit, 1.)
 BK.get_periodic_orbit(prob_col, _ci, par_sl)
 BK.getmaximum(prob_col, _ci, par_sl)
-@test BK.∂(sin, 2)(0.) == 0
+@test BK.∂(sin, Val(2))(0.) == 0
 prob_col(_ci, par_sl) #|> scatter
 BK.get_time_slices(prob_col, _ci)