# How to implement AD for this function?

I am solving a system of equations with the iterative proportional fitting procedure. That means that at each iteration I solve for a subset of the variables `μ0m`, then hold those fixed and solve for the other subset of variables `μ0f`. I am storing the residuals of the system in preallocated arrays `resm` and `resf`. I programmed a function `sys!` that computes the residuals of the system for each of those subsets:

``````const szA = 9
const szE = 2
const szG = 2
function sys!(resm, resf, μ0m, μ0f; Sm=Sm, Sf=Sf, Π=Π)
μ0 = reshape(vcat(vec(μ0m), vec(μ0f)), (szA, szE, szG))
μ0m′ = μ0[:, :, 1]
μ0f′ = μ0[:, :, 2]

μ0m′[1, szE] = 0.0
μ0f′[1, szE] = 0.0

prodsummand = similar(μ0, (szA, szE, szA, szE))
prodsummand!(prodsummand, μ0m′, μ0f′; Sm=Sm, Sf=Sf, Π=Π)

for idx in CartesianIndices(resm)
am, em = Tuple(idx)
resm[am, em] = Sm[am, em] - abs(μ0m′[am,em]) - sum( prodsummand[am,em,af,ef] for af in 1:szA, ef in 1:szE )
end
for idx in CartesianIndices(resf)
af, ef = Tuple(idx)
resf[af, ef] = Sf[af, ef] - abs(μ0f′[af,ef]) - sum( prodsummand[am,em,af,ef] for am in 1:szA, em in 1:szE )
end
return nothing
end

function prodsummand!(prodsum, μ0m′, μ0f′; Sm=Sm, Sf=Sf, Π=Π)
for idx in CartesianIndices(prodsum)
am, em, af, ef = Tuple(idx)
prodsum[am, em, af, ef] = Π[am,em,af,ef]*sqrt(Sm[am,em]*Sf[af,ef])*prod( auxprod((μ0m′[am+k,em]*μ0f′[af+k,ef])/(Sm[am+k,em]*Sf[af+k,ef]))^(0.5*(β*(1-δ))^k) for k in 0:(maxT(am,af; Z=szA)-1) )
end
return prodsum
end

function maxT(am, af; Z=szA)
Z - max(am, af) + 1
end

auxprod(x) = isfinite(abs(x)) ? abs(x) : 1.0
``````

To solve the system, I am using `NLsolve` with finite differences, and it finds a solution. I am now wanting to speed up the computation by using automatic differentiation. I am having trouble specifying the Jacobian with AD. This is the algorithm that solves the system with finite diff:

``````using NLsolve
using LinearAlgebra

function ipfp(init, sys!; Sm=Sm, Sf=Sf, Π=Π, method=:trust_region, autodiff=:finite, outmaxiter=5_000, outtol=1e-8, opts...)
μ0f0 = vec(init)

μ0m1 = similar(init)
μ0f1 = similar(init)
d = Inf
i = 0
while (i <= outmaxiter) && (outtol < d)
solm = nlsolve( (resm, μ0m) -> sys!(resm, resf, μ0m, μ0f0; Sm=Sm, Sf=Sf, Π=Π), init; method=method, autodiff=autodiff, opts... )
if converged(solm)
μ0m1 = abs.(vec(solm.zero))
else
println("solving for μ0m given μ0f did not converge at iteration ", i)
return nothing
end

solf = nlsolve( (resf, μ0f) -> sys!(resm, resf, μ0m1, μ0f; Sm=Sm, Sf=Sf, Π=Π), init; method=method, autodiff=autodiff, opts... )
if converged(solf)
μ0f1 = abs.(vec(solf.zero))
else
println("solving for μ0f given μ0m did not converge at iteration ", i)
return nothing
end
d = norm(μ0f1 - μ0f0)
μ0f0 = copy(μ0f1)
i += 1
end
ret = reshape(vcat(μ0m1, μ0f1), (szA, szE, szG))
ret[1, szE, 1] = 0.0
ret[1, szE, 2] = 0.0
return ret
end
``````

and here are some objects that complete a working example:

``````Sm = rand(szA, szE)
Sf = rand(szA, szE)
Π = rand(szA, szE, szA, szE)
resm = similar(Sm)
resf = similar(Sf)
init = rand(szA, szE)
``````

I tried adding AD-able functions to the algorithm like

``````    lfunm! = (resm, μ0m) -> sys!(resm, resf, μ0m, μ0f0; Sm=Sm, Sf=Sf, Π=Π)
sys_jacm! = ForwardDiff.jacobian!(lfunm!, resm, init)
``````

but I get an error:

``````ERROR: MethodError: no method matching jacobian!(::var"#94#96"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, ::Matrix{Float64}, ::Matrix{Float64})
``````

How can I specify the (partial) Jacobian of `sys!` using AD to pass it to `NLsolve`?

I believe you accidentally swapped the first and second argument of `jacobian!()`.

With this algorithm

``````function ipfp_ad(init, sys!; Sm=Sm, Sf=Sf, Π=Π, method=:trust_region, outmaxiter=5_000, outtol=1e-8, opts...)
μ0f0 = vec(init)

μ0m1 = similar(init)
μ0f1 = similar(init)
d = Inf
i = 0
while (i <= outmaxiter) && (outtol < d)
lfunm! = (resm, μ0m) -> sys!(resm, resf, μ0m, μ0f0; Sm=Sm, Sf=Sf, Π=Π)
sys_jacm! = ForwardDiff.jacobian!(resm, lfunm!, init)
solm = nlsolve( lfunm!, sys_jacm!, init; method=method, opts... )
if converged(solm)
μ0m1 = abs.(vec(solm.zero))
else
println("solving for μ0m given μ0f did not converge at iteration ", i)
return nothing
end

lfunm! = (resf, μ0f) -> sys!(resm, resf, μ0m1, μ0f; Sm=Sm, Sf=Sf, Π=Π)
sys_jacf! = ForwardDiff.jacobian!(resf, lfunf!, μ0f)
solf = nlsolve( lfunf!, sys_jacf!, init; method=method, opts... )
if converged(solf)
μ0f1 = abs.(vec(solf.zero))
else
println("solving for μ0f given μ0m did not converge at iteration ", i)
return nothing
end
d = norm(μ0f1 - μ0f0)
μ0f0 = copy(μ0f1)
i += 1
end
ret = reshape(vcat(μ0m1, μ0f1), (szA, szE, szG))
ret[1, szE, 1] = 0.0
ret[1, szE, 2] = 0.0
return ret
end
``````

I still get an error:

``````ERROR: MethodError: no method matching (::var"#71#73"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)})(::Matrix{ForwardDiff.Dual{ForwardDiff.Tag{var"#71#73"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, Float64}, Float64, 9}})
``````

` sys_jacf!` is not a function, it is an array. However, nlsolve expects a function to calculate the jacobian. Please also provide the full stack trace. From the error it impossible to figure out on which line of your code the problem is

1 Like

Here I post the full script:

``````using NLsolve
using LinearAlgebra
using ForwardDiff

const szA = 9
const szE = 2
const szG = 2

function sys!(resm, resf, μ0m, μ0f; Sm=Sm, Sf=Sf, Π=Π)
μ0 = reshape(vcat(vec(μ0m), vec(μ0f)), (szA, szE, szG))
μ0m′ = μ0[:, :, 1]
μ0f′ = μ0[:, :, 2]

μ0m′[1, szE] = 0.0
μ0f′[1, szE] = 0.0

prodsummand = similar(μ0, (szA, szE, szA, szE))
prodsummand!(prodsummand, μ0m′, μ0f′; Sm=Sm, Sf=Sf, Π=Π)

for idx in CartesianIndices(resm)
am, em = Tuple(idx)
resm[am, em] = Sm[am, em] - abs(μ0m′[am,em]) - sum( prodsummand[am,em,af,ef] for af in 1:szA, ef in 1:szE )
end
for idx in CartesianIndices(resf)
af, ef = Tuple(idx)
resf[af, ef] = Sf[af, ef] - abs(μ0f′[af,ef]) - sum( prodsummand[am,em,af,ef] for am in 1:szA, em in 1:szE )
end
return nothing
end

function prodsummand!(prodsum, μ0m′, μ0f′; Sm=Sm, Sf=Sf, Π=Π)
for idx in CartesianIndices(prodsum)
am, em, af, ef = Tuple(idx)
prodsum[am, em, af, ef] = Π[am,em,af,ef]*sqrt(Sm[am,em]*Sf[af,ef])*prod( auxprod((μ0m′[am+k,em]*μ0f′[af+k,ef])/(Sm[am+k,em]*Sf[af+k,ef]))^(0.5*(β*(1-δ))^k) for k in 0:(maxT(am,af; Z=szA)-1) )
end
return prodsum
end

function maxT(am, af; Z=szA)
Z - max(am, af) + 1
end

auxprod(x) = isfinite(abs(x)) ? abs(x) : 1.0

function ipfp_ad(init, sys!; Sm=Sm, Sf=Sf, Π=Π, method=:trust_region, autodiff=:finite, outmaxiter=5_000, outtol=1e-8, opts...)
μ0f0 = vec(init)

μ0m1 = similar(init)
μ0f1 = similar(init)
d = Inf
i = 0
while (i <= outmaxiter) && (outtol < d)
lfunm! = (resm, μ0m) -> sys!(resm, resf, μ0m, μ0f0; Sm=Sm, Sf=Sf, Π=Π)
sys_jacm! = ForwardDiff.jacobian!(resm, lfunm!, init)
solm = nlsolve( (resm, μ0m) -> sys!(resm, resf, μ0m, μ0f0; Sm=Sm, Sf=Sf, Π=Π), init; method=method, autodiff=autodiff, opts... )
if converged(solm)
μ0m1 = abs.(vec(solm.zero))
else
println("solving for μ0m given μ0f did not converge at iteration ", i)
return nothing
end

lfunf! = (resf, μ0f) -> sys!(resm, resf, μ0m0, μ0f; Sm=Sm, Sf=Sf, Π=Π)
sys_jacf! = ForwardDiff.jacobian!(resf, lfunf!, μ0m1)
solf = nlsolve( (resf, μ0f) -> sys!(resm, resf, μ0m1, μ0f; Sm=Sm, Sf=Sf, Π=Π), init; method=method, autodiff=autodiff, opts... )
if converged(solf)
μ0f1 = abs.(vec(solf.zero))
else
println("solving for μ0f given μ0m did not converge at iteration ", i)
return nothing
end
d = norm(μ0f1 - μ0f0)
μ0f0 = copy(μ0f1)
i += 1
end
ret = reshape(vcat(μ0m1, μ0f1), (szA, szE, szG))
ret[1, szE, 1] = 0.0
ret[1, szE, 2] = 0.0
return ret
end

Sm = rand(szA, szE)
Sf = rand(szA, szE)
Π = rand(szA, szE, szA, szE)
resm = similar(Sm)
resf = similar(Sf)
init = rand(szA, szE)
const β = 0.92
const δ = 0.10
``````

Here is the error

``````julia> ipfp_sol = ipfp_ad(init, sys!; Sm=Sm, Sf=Sf, store_trace=true, extended_trace=true, method=:trust_region)
ERROR: MethodError: no method matching (::var"#15#19"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)})(::Matrix{ForwardDiff.Dual{ForwardDiff.Tag{var"#15#19"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, Float64}, Float64, 9}})
Closest candidates are:
(::var"#15#19")(::Any, ::Any) at /Users/amrods/discourseQ4.jl:53
Stacktrace:
[1] chunk_mode_jacobian!(result::Matrix{Float64}, f::var"#15#19"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, x::Matrix{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{var"#15#19"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, Float64}, Float64, 9, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{var"#15#19"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, Float64}, Float64, 9}}})
@ ForwardDiff ~/.julia/packages/ForwardDiff/QOqCN/src/jacobian.jl:223
[2] jacobian!(result::Matrix{Float64}, f::Function, x::Matrix{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{var"#15#19"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, Float64}, Float64, 9, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{var"#15#19"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, Float64}, Float64, 9}}}, ::Val{true})
@ ForwardDiff ~/.julia/packages/ForwardDiff/QOqCN/src/jacobian.jl:60
[3] jacobian!(result::Matrix{Float64}, f::Function, x::Matrix{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{var"#15#19"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, Float64}, Float64, 9, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{var"#15#19"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, Float64}, Float64, 9}}}) (repeats 2 times)
@ ForwardDiff ~/.julia/packages/ForwardDiff/QOqCN/src/jacobian.jl:56
[4] ipfp_ad(init::Matrix{Float64}, sys!::typeof(sys!); Sm::Matrix{Float64}, Sf::Matrix{Float64}, Π::Array{Float64, 4}, method::Symbol, autodiff::Symbol, outmaxiter::Int64, outtol::Float64, opts::Base.Iterators.Pairs{Symbol, Bool, Tuple{Symbol, Symbol}, NamedTuple{(:store_trace, :extended_trace), Tuple{Bool, Bool}}})
@ Main ~/discourseQ4.jl:54
[5] top-level scope
@ REPL[1]:1
``````

Line 53 is where I define `lfunm!`.

so this should do it?

``````sys_jacm! = (resm, x) -> ForwardDiff.jacobian!(resm, lfunm!, x)
``````

I changed the lines defining the Jacobian:

``````using NLsolve
using LinearAlgebra
using ForwardDiff

const szA = 9
const szE = 2
const szG = 2
const β = 0.92
const δ = 0.10

function sys!(resm, resf, μ0m, μ0f; Sm=Sm, Sf=Sf, Π=Π)
μ0 = reshape(vcat(vec(μ0m), vec(μ0f)), (szA, szE, szG))
μ0m′ = μ0[:, :, 1]
μ0f′ = μ0[:, :, 2]

μ0m′[1, szE] = 0.0
μ0f′[1, szE] = 0.0

prodsummand = similar(μ0, (szA, szE, szA, szE))
prodsummand!(prodsummand, μ0m′, μ0f′; Sm=Sm, Sf=Sf, Π=Π)

for idx in CartesianIndices(resm)
am, em = Tuple(idx)
resm[am, em] = Sm[am, em] - abs(μ0m′[am,em]) - sum( prodsummand[am,em,af,ef] for af in 1:szA, ef in 1:szE )
end
for idx in CartesianIndices(resf)
af, ef = Tuple(idx)
resf[af, ef] = Sf[af, ef] - abs(μ0f′[af,ef]) - sum( prodsummand[am,em,af,ef] for am in 1:szA, em in 1:szE )
end
return nothing
end

function prodsummand!(prodsum, μ0m′, μ0f′; Sm=Sm, Sf=Sf, Π=Π)
for idx in CartesianIndices(prodsum)
am, em, af, ef = Tuple(idx)
prodsum[am, em, af, ef] = Π[am,em,af,ef]*sqrt(Sm[am,em]*Sf[af,ef])*prod( auxprod((μ0m′[am+k,em]*μ0f′[af+k,ef])/(Sm[am+k,em]*Sf[af+k,ef]))^(0.5*(β*(1-δ))^k) for k in 0:(maxT(am,af; Z=szA)-1) )
end
return prodsum
end

function maxT(am, af; Z=szA)
Z - max(am, af) + 1
end

auxprod(x) = isfinite(abs(x)) ? abs(x) : 1.0

function ipfp_ad(init, sys!; Sm=Sm, Sf=Sf, Π=Π, method=:trust_region, outmaxiter=5_000, outtol=1e-8, opts...)
μ0f0 = vec(init)

μ0m1 = similar(init)
μ0f1 = similar(init)
d = Inf
i = 0
while (i <= outmaxiter) && (outtol < d)
lfunm!(resx, x) = sys!(resx, resf, x, μ0f0; Sm=Sm, Sf=Sf, Π=Π)
sys_jacm!(resx, x) = ForwardDiff.jacobian!(resx, lfunm!, x)
solm = nlsolve( lfunm!, sys_jacm!, init; method=method, opts... )
if converged(solm)
μ0m1 = abs.(vec(solm.zero))
else
println("solving for μ0m given μ0f did not converge at iteration ", i)
return nothing
end

lfunf!(resx, x) =  sys!(resm, resx, μ0m1, x; Sm=Sm, Sf=Sf, Π=Π)
sys_jacf!(resx, x) = ForwardDiff.jacobian!(resx, lfunf!, x)
solf = nlsolve( lfunf!, sys_jacf!, init; method=method, opts... )
if converged(solf)
μ0f1 = abs.(vec(solf.zero))
else
println("solving for μ0f given μ0m did not converge at iteration ", i)
return nothing
end
d = norm(μ0f1 - μ0f0)
μ0f0 = copy(μ0f1)
i += 1
end
ret = reshape(vcat(μ0m1, μ0f1), (szA, szE, szG))
ret[1, szE, 1] = 0.0
ret[1, szE, 2] = 0.0
return ret
end

Sm = rand(szA, szE)
Sf = rand(szA, szE)
Π = rand(szA, szE, szA, szE)
resm = similar(Sm)
resf = similar(Sf)
init = rand(szA, szE)

sol = ipfp_ad(init, sys!; Sm=Sm, Sf=Sf, Π=Π, method=:trust_region, outmaxiter=5_000, outtol=1e-8)
``````

I get this error:

``````ERROR: LoadError: MethodError: no method matching (::var"#lfunm!#15"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)})(::Matrix{ForwardDiff.Dual{ForwardDiff.Tag{var"#lfunm!#15"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, Float64}, Float64, 9}})
Closest candidates are:
(::var"#lfunm!#15")(::Any, ::Any) at /Users/amrods/discourseQ4.jl:55
Stacktrace:
[1] chunk_mode_jacobian!(result::Matrix{Float64}, f::var"#lfunm!#15"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, x::Matrix{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{var"#lfunm!#15"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, Float64}, Float64, 9, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{var"#lfunm!#15"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, Float64}, Float64, 9}}})
@ ForwardDiff ~/.julia/packages/ForwardDiff/QOqCN/src/jacobian.jl:223
[2] jacobian!(result::Matrix{Float64}, f::Function, x::Matrix{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{var"#lfunm!#15"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, Float64}, Float64, 9, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{var"#lfunm!#15"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, Float64}, Float64, 9}}}, ::Val{true})
@ ForwardDiff ~/.julia/packages/ForwardDiff/QOqCN/src/jacobian.jl:60
[3] jacobian!(result::Matrix{Float64}, f::Function, x::Matrix{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{var"#lfunm!#15"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, Float64}, Float64, 9, Matrix{ForwardDiff.Dual{ForwardDiff.Tag{var"#lfunm!#15"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, Float64}, Float64, 9}}}) (repeats 2 times)
@ ForwardDiff ~/.julia/packages/ForwardDiff/QOqCN/src/jacobian.jl:56
[4] (::var"#sys_jacm!#16"{var"#lfunm!#15"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}})(resx::Matrix{Float64}, x::Matrix{Float64})
@ Main ~/discourseQ4.jl:56
[5] (::NLSolversBase.var"#fj!#7"{var"#lfunm!#15"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}, var"#sys_jacm!#16"{var"#lfunm!#15"{Matrix{Float64}, Matrix{Float64}, Array{Float64, 4}, typeof(sys!)}}})(fx::Matrix{Float64}, jx::Matrix{Float64}, x::Matrix{Float64})
@ NLSolversBase ~/.julia/packages/NLSolversBase/geyh3/src/objective_types/abstract.jl:7
[6] value_jacobian!!(obj::OnceDifferentiable{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}, F::Matrix{Float64}, J::Matrix{Float64}, x::Matrix{Float64})
@ NLSolversBase ~/.julia/packages/NLSolversBase/geyh3/src/interface.jl:124
[7] value_jacobian!!
@ ~/.julia/packages/NLSolversBase/geyh3/src/interface.jl:122 [inlined]
[8] trust_region_(df::OnceDifferentiable{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}, initial_x::Matrix{Float64}, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, factor::Float64, autoscale::Bool, cache::NLsolve.NewtonTrustRegionCache{Matrix{Float64}})
@ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/solvers/trust_region.jl:119
[9] trust_region (repeats 2 times)
@ ~/.julia/packages/NLsolve/gJL1I/src/solvers/trust_region.jl:235 [inlined]
[10] nlsolve(df::OnceDifferentiable{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}, initial_x::Matrix{Float64}; method::Symbol, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, linesearch::Static, linsolve::NLsolve.var"#27#29", factor::Float64, autoscale::Bool, m::Int64, beta::Int64, aa_start::Int64, droptol::Float64)
@ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/nlsolve/nlsolve.jl:26
[11] nlsolve(f::Function, j::Function, initial_x::Matrix{Float64}; inplace::Bool, kwargs::Base.Iterators.Pairs{Symbol, Symbol, Tuple{Symbol}, NamedTuple{(:method,), Tuple{Symbol}}})
@ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/nlsolve/nlsolve.jl:67
[12] ipfp_ad(init::Matrix{Float64}, sys!::typeof(sys!); Sm::Matrix{Float64}, Sf::Matrix{Float64}, Π::Array{Float64, 4}, method::Symbol, outmaxiter::Int64, outtol::Float64, opts::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ Main ~/discourseQ4.jl:57
[13] top-level scope
@ ~/discourseQ4.jl:91
in expression starting at /Users/amrods/discourseQ4.jl:91
``````

Line 55 is `lfunm!(resx, x) = sys!(resx, resf, x, μ0f0; Sm=Sm, Sf=Sf, Π=Π)`

yes that looks good

You’re still calling `jacobian!()` incorrectly. For your case you need to provide 4 arguments. Please read the documentation—it is well explained there.

2 Likes

yeah, i missed that bang at the end, see A bug in ForwardDiff? - #4 by amrods