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