# Domain issue when working with nlsolve

**URL:** <https://discourse.julialang.org/t/domain-issue-when-working-with-nlsolve/98193>\
**Category:** New to Julia\
**Tags:** question, optimization, modelling\
**Created:** [May 2, 2023, 12:23pm UTC](https://discourse.julialang.org/t/domain-issue-when-working-with-nlsolve/98193 "2023-05-02T12:23:28Z")\
**Posts on this page:** 17\
**Page:** 1

<div class="post-metadata">

**Author:** ![cdestael](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cdestael/32/46700_2.png) [@cdestael](https://discourse.julialang.org/u/cdestael)\
**Post date:** [May 2, 2023, 12:23pm UTC](https://discourse.julialang.org/t/domain-issue-when-working-with-nlsolve/98193/1 "2023-05-02T12:23:28Z")

</div>

I am quite stumped on this problem so appreciate in advance any suggestions. I am currently working with a trade model and trying to solve a system of 4 equations using “nlsolve.” I am running into issues with the domain and have tried a few workarounds, but new errors pop up. Below is a MWE of the larger problem. Thanks again.

Function used in code:

```julia
function sys_tHE_AB!(F, Z_AB, k_e, fe, Fe, param, j) # k_e: (1000, 10, 11) fe: (11, 11) Fe: (11, 11)
    μ = param.μ
    σ = param.σ
    μj = μ[j]
    σj = σ[j]
    ρ = param.ρ
    δ = param.δ
    β = sqrt(2*(ρ + δ)/(σj^2) + (μj/(σj^2) - 1/2)^2) - (μj/(σj^2) - 1/2)
    α = -sqrt(2*(ρ + δ)/(σj^2) + (μj/(σj^2) - 1/2)^2) - (μj/(σj^2) - 1/2)
    disc = ρ + δ - μj
    F[1] = Z_AB[3]*Z_AB[1]^α + k_e*Z_AB[1]/disc - fe/disc - Fe - Z_AB[4]*Z_AB[1]^β
    F[2] = Z_AB[3]*Z_AB[2]^α + k_e*Z_AB[2]/disc - fe/disc - Z_AB[4]*Z_AB[2]^β
    F[3] = α*Z_AB[3]*Z_AB[1]^(α-1) + k_e/disc - β*Z_AB[4]*Z_AB[1]^(β-1)
    F[4] = α*Z_AB[3]*Z_AB[2]^(α-1) + k_e/disc - β*Z_AB[4]*Z_AB[2]^(β-1)
end

```

And the main code that calls the function:

```julia
N = 10 
USA = N + 1
ρ = 0.05 # rate of time preference
δ = 0.05 # exogenous death rate (here also endogenous exit)
η = 5.0 # elasticity of substitution
H = η^(-η)*(η-1)^(η-1)
P = ones(1,N+1) # price indexes
n = 1000
T = 31 # number of periods
s = 100
b = 1 # lower bound of the (Pareto) productivity distribution
θ = 3 # shape parameter of the (Pareto) productivity distribution 
z=10*ones(1000,1)
Y_simu=ones(1,11)
Q=10*ones(1,11)
fh=10*ones(1,10)
Fh=10*ones(1,10)
τ=ones(11,11)
w=ones(1,11)
fe=ones(11,11)
Fe=ones(11,11)
μ_USj=10*ones(11,1)
σ_USj=10*ones(11,1)
struct Parameter 
    μ::Matrix{Float64}
    σ::Matrix{Float64}
    ρ::Float64
    δ::Float64
    η::Float64
end
β = sqrt.(2*(ρ + δ)./(σ_USj.^2) .+ (μ_USj./(σ_USj.^2) .- 1/2).^2) .- (μ_USj./(σ_USj.^2) .- 1/2)
α = -sqrt.(2*(ρ + δ)./(σ_USj.^2) .+ (μ_USj./(σ_USj.^2) .- 1/2).^2) .- (μ_USj./(σ_USj.^2) .- 1/2)
disc = ρ + δ .- μ_USj
disc = disc'

# Preallocation (n firms, N host countries + the US)
k_e = zeros(n,N,N+1)
Bo_0 = zeros(n,N,N+1)
Z_HE_0 = 10*zeros(n,N,N+1)
Z_EH_0 = 10*zeros(n,N,N+1)
init_HE = zeros(n,N,N+1,4)
thresholds_HE_sol = zeros(n,N,N+1,4)
init_OH = zeros(n,N,4)
thresholds_OH_sol = zeros(n,N,T,4)
status_exp = zeros(n,N,N+1,T)
status = zeros(n,N,T)
k_h = H*(1 ./ repeat(z, 1, N)) .^(1-η) .* repeat((P[:,1:10] .^η) .* Q[:,1:10],n)
Z_OH_0 = repeat(β[1:N]' ./ (β[1:N].-1)', n, 1) .* (k_h.^-1) .* (repeat(fh,n,1) + (ρ.+δ.-repeat(μ_USj[1:10]',n,1)) .* repeat(Fh,n,1))
Z_HO_0 = (repeat(β[1:N]' ./ (β[1:N].-1)', n, 1).* (k_h.^-1) .* repeat(fh,n,1)) 
Bo_0_H= 100*ones(1000,10) # k_h ./ repeat(disc[:,1:N],n,1) ./ repeat(β[1:N]',n,1) ./ Z_OH_0.^(repeat(β[1:N]',n,1).-1)
for c in 1:N    
    for d in 1:N+1
        if d != c
            k_e[:,c,d] = H*(τ[c,d]*w[c]./z./w[d]).^(1-η).*P[d]^η*Q[d] 
            Z_HE_0[:,c,d] = (β[d]/(β[d].-1).*k_e[:,c,d].^-1 .* (fe[c,d] + (ρ+δ-μ_USj[d])*Fe[c,d]))   
            Z_EH_0[:,c,d] = (β[d]/(β[d].-1).*k_e[:,c,d].^-1 .* fe[c,d])        
            Bo_0[:,c,d] = k_e[:,c,d]./disc[d]./β[d]./(Z_HE_0[:,c,d].^(β[d].-1))
        end
    end
end

# SOLVING FOR THE THRESHOLDS AND THE PARAMETERS OF THE VALUE FUNCTIONS
for j in 1:2
    for c in 1:N    
        for d in 1:N+1
            if d != c
                k_e[:,c,d] = H*(τ[c,d]*w[c]./z./w[d]).^(1-η).*P[d]^η*Q[d] # 1000 x 1 vector
                Z_HE_0[:,c,d] = (β[d]/(β[d].-1).*k_e[:,c,d].^-1 .* (fe[c,d] + (ρ+δ-μ_USj[d])*Fe[c,d])) # 1000 x 1 vector   
                Z_EH_0[:,c,d] = (β[d]/(β[d].-1).*k_e[:,c,d].^-1 .* fe[c,d]) # 1000 x 1 vector        
                Bo_0[:,c,d] = k_e[:,c,d]./disc[d]./β[d]./(Z_HE_0[:,c,d].^(β[d].-1)) # 1000 x 1 vector
            end
        end
    end
    # initialization
    k = N+1
    thresholds_HE_sol = zeros(n,N,N+1,4)
    for i in eachindex(z)
        if k != j
            # creating initial Hessian
            init_HE[i,j,k,:] = [Z_HE_0[i,j,k] Z_EH_0[i,j,k] 0 Bo_0[i,j,k]] # 1 x 4 vector
            if i > 20 && j != 3 && isnan(sum(thresholds_HE_sol[i-1,j,k,:])) == 0
                init_HE[i,j,k,:] = thresholds_HE_sol[i-1,j,k,:]
            end
            if init_HE[i,j,k,3] < 0
                init_HE[i,j,k,3] = 0.05
            end

            # initial matrices and function to be minimized
            initial_x = init_HE[i,j,k,:]    
            f_thresholds_HE!(F, Z_AB) = sys_tHE_AB!(F, Z_AB, k_e[i,j,k], fe[j,k], Fe[j,k], param, j)           
            result = nlsolve(f_thresholds_HE!, initial_x, method = :trust_region, xtol=10e-5, ftol=10e-8, iterations=10000, store_trace=true, extended_trace=true).zero
        end
    end
end

```

The analogous code in Matlab (below) runs without error. Thank you in advance for helping me understand the difference.

```julia
f_thresholds_HE = @(Z_thresholds_exp)sys_tHE_AB(Z_thresholds_exp, k_e(i,j,k), fe(j,k), Fe(j,k),param,k);
[thresholds_HE_sol(i,j,k,:),fval] = fsolve(f_thresholds_HE,init_HE(i,j,k,:),options);

```

---

<div class="post-metadata">

**Author:** ![Krastanov](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/krastanov/32/6817_2.png) [@Krastanov](https://discourse.julialang.org/u/Krastanov)\
**Post date:** [May 2, 2023, 2:45pm UTC](https://discourse.julialang.org/t/domain-issue-when-working-with-nlsolve/98193/2 "2023-05-02T14:45:17Z")

</div>

It might take a bit for someone to attempt to help, as the example you are sharing is very heavy. In the meantime, something that can help the helpers:

- Any chance you can simplify this further? Could you give just one single `nlsolve` call that fails in an obvious way. Currently you have a very large amount of code dedicated to just setup and the nlsolve is called multiple times inside of a deep loop.

- Could you elaborate on how it fails? Is there an error being raised (in which case, could you share the exception stack trace)? Is there an answer, but the answer is obviously wrong (in which case, could you explain why it is obviously wrong)?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 2, 2023, 3:47pm UTC](https://discourse.julialang.org/t/domain-issue-when-working-with-nlsolve/98193/3 "2023-05-02T15:47:02Z")

</div>

It’s hard to see where exactly it is because the error message is missing from your post, but probably the simplest way to solve it could be to use NaNMath.jl on whichever of the `^` calls is failing so that it throws a NaN instead of a DomainError, which would help the system recover and re-adapt its steps.

---

<div class="post-metadata">

**Author:** ![cdestael](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cdestael/32/46700_2.png) [@cdestael](https://discourse.julialang.org/u/cdestael)\
**Post date:** [May 2, 2023, 8:00pm UTC](https://discourse.julialang.org/t/domain-issue-when-working-with-nlsolve/98193/4 "2023-05-02T20:00:08Z")

</div>

Thank you both @Krastanov and @ChrisRackauckas for your feedback. @ChrisRackauckas, thanks so much for the NaNMath.jl suggestion!! I implemented it, but now am getting a new error, which I think reflects the fact that nlsolve is not ignoring the “NaN”? Any ideas on how to fix this new problem? I am elaborating below on the error code and including the stack trace. Thanks a lot again.

Code:

```julia
# Description: this file calibrates trade costs from each country to US
using MAT, Statistics, DataFrames, LinearAlgebra
using NLsolve
using NaNMath; nm=NaNMath

function sys_tHE_AB!(F, Z_AB, k_e, fe, Fe, param, j) # k_e: (1000, 10, 11) fe: (11, 11) Fe: (11, 11)
    μ = param.μ
    σ = param.σ
    μj = μ[j]
    σj = σ[j]
    ρ = param.ρ
    δ = param.δ
    β = sqrt(2*(ρ + δ)/(σj^2) + (μj/(σj^2) - 1/2)^2) - (μj/(σj^2) - 1/2)
    α = -sqrt(2*(ρ + δ)/(σj^2) + (μj/(σj^2) - 1/2)^2) - (μj/(σj^2) - 1/2)
    disc = ρ + δ - μj
    
    F[1] = Z_AB[3]*nm.Z_AB[1]^α + k_e*Z_AB[1]/disc - fe/disc - Fe - Z_AB[4]*nm.Z_AB[1]^β
    F[2] = Z_AB[3]*nm.Z_AB[2]^α + k_e*Z_AB[2]/disc - fe/disc - Z_AB[4]*nm.Z_AB[2]^β
    F[3] = α*Z_AB[3]*nm.Z_AB[1]^(α-1) + k_e/disc - β*Z_AB[4]*nm.Z_AB[1]^(β-1)
    F[4] = α*Z_AB[3]*nm.Z_AB[2]^(α-1) + k_e/disc - β*Z_AB[4]*nm.Z_AB[2]^(β-1)
end

N = 10 
USA = N + 1
ρ = 0.05 # rate of time preference
δ = 0.05 # exogenous death rate (here also endogenous exit)
η = 5.0 # elasticity of substitution
H = η^(-η)*(η-1)^(η-1)
P = ones(1,N+1) # price indexes
n = 1000
T = 31 # number of periods
s = 100
b = 1 # lower bound of the (Pareto) productivity distribution
θ = 3 # shape parameter of the (Pareto) productivity distribution 
z=10*ones(1000,1)
Y_simu=ones(1,11)
Q=10*ones(1,11)
fh=10*ones(1,10)
Fh=10*ones(1,10)
τ=ones(11,11)
w=ones(1,11)
fe=ones(11,11)
Fe=ones(11,11)
μ_USj=10*ones(11,1)
σ_USj=10*ones(11,1)

# creating a structure of the parameters
struct Parameter 
    μ::Matrix{Float64}
    σ::Matrix{Float64}
    ρ::Float64
    δ::Float64
    η::Float64
end
param=Parameter(matread("brownian.mat")["mu_USj"], matread("brownian.mat")["sigma_USj"], ρ, δ, η)
β = sqrt.(2*(ρ + δ)./(σ_USj.^2) .+ (μ_USj./(σ_USj.^2) .- 1/2).^2) .- (μ_USj./(σ_USj.^2) .- 1/2)
α = -sqrt.(2*(ρ + δ)./(σ_USj.^2) .+ (μ_USj./(σ_USj.^2) .- 1/2).^2) .- (μ_USj./(σ_USj.^2) .- 1/2)
disc = ρ + δ .- μ_USj
disc = disc'

# Preallocation (n firms, N host countries + the US)
k_e = zeros(n,N,N+1)
Bo_0 = zeros(n,N,N+1)
Z_HE_0 = 10*zeros(n,N,N+1)
Z_EH_0 = 10*zeros(n,N,N+1)
init_HE = zeros(n,N,N+1,4)
thresholds_HE_sol = zeros(n,N,N+1,4)
init_OH = zeros(n,N,4)
thresholds_OH_sol = zeros(n,N,T,4)
status_exp = zeros(n,N,N+1,T)
status = zeros(n,N,T)
k_h = H*(1 ./ repeat(z, 1, N)) .^(1-η) .* repeat((P[:,1:10] .^η) .* Q[:,1:10],n)
Z_OH_0 = repeat(β[1:N]' ./ (β[1:N].-1)', n, 1) .* (k_h.^-1) .* (repeat(fh,n,1) + (ρ.+δ.-repeat(μ_USj[1:10]',n,1)) .* repeat(Fh,n,1))
Z_HO_0 = (repeat(β[1:N]' ./ (β[1:N].-1)', n, 1).* (k_h.^-1) .* repeat(fh,n,1)) 
Bo_0_H= 100*ones(1000,10) # k_h ./ repeat(disc[:,1:N],n,1) ./ repeat(β[1:N]',n,1) ./ Z_OH_0.^(repeat(β[1:N]',n,1).-1)
for c in 1:N    
    for d in 1:N+1
        if d != c
            k_e[:,c,d] = H*(τ[c,d]*w[c]./z./w[d]).^(1-η).*P[d]^η*Q[d] 
            Z_HE_0[:,c,d] = (β[d]/(β[d].-1).*k_e[:,c,d].^-1 .* (fe[c,d] + (ρ+δ-μ_USj[d])*Fe[c,d]))   
            Z_EH_0[:,c,d] = (β[d]/(β[d].-1).*k_e[:,c,d].^-1 .* fe[c,d])        
            Bo_0[:,c,d] = k_e[:,c,d]./disc[d]./β[d]./(Z_HE_0[:,c,d].^(β[d].-1))
        end
    end
end

# SOLVING FOR THE THRESHOLDS AND THE PARAMETERS OF THE VALUE FUNCTIONS
for j in 1:2
    for c in 1:N    
        for d in 1:N+1
            if d != c
                k_e[:,c,d] = H*(τ[c,d]*w[c]./z./w[d]).^(1-η).*P[d]^η*Q[d] # 1000 x 1 vector
                Z_HE_0[:,c,d] = (β[d]/(β[d].-1).*k_e[:,c,d].^-1 .* (fe[c,d] + (ρ+δ-μ_USj[d])*Fe[c,d])) # 1000 x 1 vector   
                Z_EH_0[:,c,d] = (β[d]/(β[d].-1).*k_e[:,c,d].^-1 .* fe[c,d]) # 1000 x 1 vector        
                Bo_0[:,c,d] = k_e[:,c,d]./disc[d]./β[d]./(Z_HE_0[:,c,d].^(β[d].-1)) # 1000 x 1 vector
            end
        end
    end
    # initialization
    k = N+1
    thresholds_HE_sol = zeros(n,N,N+1,4)
    for i in eachindex(z)
        if k != j
            # creating initial Hessian
            init_HE[i,j,k,:] = [Z_HE_0[i,j,k] Z_EH_0[i,j,k] 0 Bo_0[i,j,k]] # 1 x 4 vector
            if i > 20 && j != 3 && isnan(sum(thresholds_HE_sol[i-1,j,k,:])) == 0
                init_HE[i,j,k,:] = thresholds_HE_sol[i-1,j,k,:]
            end
            if init_HE[i,j,k,3] < 0
                init_HE[i,j,k,3] = 0.05
            end

            # initial matrices and function to be minimized
            initial_x = init_HE[i,j,k,:]    
            Z_AB=ones(4)
            f_thresholds_HE!(F, Z_AB) = sys_tHE_AB!(F, Z_AB, k_e[i,j,k], fe[j,k], Fe[j,k], param, j)           
            result = nlsolve(f_thresholds_HE!, initial_x, method = :trust_region, xtol=10e-5, ftol=10e-8, iterations=10000, store_trace=true, extended_trace=true).zero
        end
    end
end

```

Stack trace:

```julia
ERROR: UndefVarError: Z_AB not defined
Stacktrace:
  [1] getproperty(x::Module, f::Symbol)
    @ Base ./Base.jl:31
  [2] sys_tHE_AB!(F::Vector{Float64}, Z_AB::Vector{Float64}, k_e::Float64, fe::Float64, Fe::Float64, param::Parameter, j::Int64)
    @ Main ~/Dropbox/CORINNE_BU_RA/Julia/calibration/test.jl:17
  [3] (::var"#f_thresholds_HE!#15"{Int64, Int64, Int64})(F::Vector{Float64}, Z_AB::Vector{Float64})
    @ Main ~/Dropbox/CORINNE_BU_RA/Julia/calibration/test.jl:116
  [4] (::NLSolversBase.var"#fj_finitediff!#21"{var"#f_thresholds_HE!#15"{Int64, Int64, Int64}, FiniteDiff.JacobianCache{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, UnitRange{Int64}, Nothing, Val{:central}(), Float64}})(F::Vector{Float64}, J::Matrix{Float64}, x::Vector{Float64})
    @ NLSolversBase ~/.julia/packages/NLSolversBase/kavn7/src/objective_types/oncedifferentiable.jl:138
  [5] value_jacobian!!(obj::OnceDifferentiable{Vector{Float64}, Matrix{Float64}, Vector{Float64}}, F::Vector{Float64}, J::Matrix{Float64}, x::Vector{Float64})
    @ NLSolversBase ~/.julia/packages/NLSolversBase/kavn7/src/interface.jl:124
  [6] value_jacobian!!
    @ ~/.julia/packages/NLSolversBase/kavn7/src/interface.jl:122 [inlined]
  [7] trust_region_(df::OnceDifferentiable{Vector{Float64}, Matrix{Float64}, Vector{Float64}}, initial_x::Vector{Float64}, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, factor::Float64, autoscale::Bool, cache::NLsolve.NewtonTrustRegionCache{Vector{Float64}})
    @ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/solvers/trust_region.jl:119
  [8] trust_region (repeats 2 times)
    @ ~/.julia/packages/NLsolve/gJL1I/src/solvers/trust_region.jl:235 [inlined]
  [9] nlsolve(df::OnceDifferentiable{Vector{Float64}, Matrix{Float64}, Vector{Float64}}, initial_x::Vector{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
 [10] nlsolve(f::Function, initial_x::Vector{Float64}; method::Symbol, autodiff::Symbol, inplace::Bool, kwargs::Base.Pairs{Symbol, Real, NTuple{5, Symbol}, NamedTuple{(:xtol, :ftol, :iterations, :store_trace, :extended_trace), Tuple{Float64, Float64, Int64, Bool, Bool}}})
    @ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/nlsolve/nlsolve.jl:52
 [11] top-level scope
    @ ~/Dropbox/CORINNE_BU_RA/Julia/calibration/test.jl:117

```

---

<div class="post-metadata">

**Author:** ![Krastanov](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/krastanov/32/6817_2.png) [@Krastanov](https://discourse.julialang.org/u/Krastanov)\
**Post date:** [May 2, 2023, 8:29pm UTC](https://discourse.julialang.org/t/domain-issue-when-working-with-nlsolve/98193/5 "2023-05-02T20:29:14Z")

</div>

Check which line of your code causes that error to appear (and generally share that when reporting errors). On that line you are probably trying to access a variable called `Z_AB`. The error message is saying that you have not defined such a variable in that scope. Maybe due to a typo or a mixup between local and global variables?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 2, 2023, 8:40pm UTC](https://discourse.julialang.org/t/domain-issue-when-working-with-nlsolve/98193/6 "2023-05-02T20:40:34Z")

</div>

> [@cdestael](#):
>
> `nm.Z_AB`

That doesn’t make sense. It would be `nm.^`, `Z_AB` should be unchanged.

---

<div class="post-metadata">

**Author:** ![cdestael](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cdestael/32/46700_2.png) [@cdestael](https://discourse.julialang.org/u/cdestael)\
**Post date:** [May 3, 2023, 1:04pm UTC](https://discourse.julialang.org/t/domain-issue-when-working-with-nlsolve/98193/7 "2023-05-03T13:04:14Z")

</div>

Thanks very much @ChrisRackauckas!! The (correct) use of `nm.` solved the domain error.

Now when I run the code, it gives a method matching error. I’ve been reading online and experimenting with changes, but it is not clear to me how to resolve it. Any suggestions on what to do?  
Thank you so much again.

Updated function:

```julia
function sys_tHE_AB!(F, Z_AB, k_e, fe, Fe, param, j) 
    μ = param.μ
    σ = param.σ
    μj = μ[j]
    σj = σ[j]
    ρ = param.ρ
    δ = param.δ
    β = sqrt(2*(ρ + δ)/(σj^2) + (μj/(σj^2) - 1/2)^2) - (μj/(σj^2) - 1/2)
    α = -sqrt(2*(ρ + δ)/(σj^2) + (μj/(σj^2) - 1/2)^2) - (μj/(σj^2) - 1/2)
    disc = ρ + δ - μj

    F[1] = Z_AB[3]*Z_AB[1]nm.^α + k_e*Z_AB[1]/disc - fe/disc - Fe - Z_AB[4]*Z_AB[1]nm.^β
    F[2] = Z_AB[3]*Z_AB[2]nm.^α + k_e*Z_AB[2]/disc - fe/disc - Z_AB[4]*Z_AB[2]nm.^β
    F[3] = α*Z_AB[3]*Z_AB[1]nm.^(α-1) + k_e/disc - β*Z_AB[4]*Z_AB[1]nm.^(β-1)
    F[4] = α*Z_AB[3]*Z_AB[2]nm.^(α-1) + k_e/disc - β*Z_AB[4]*Z_AB[2]nm.^(β-1)

    F = [F[1] F[2] F[3] F[4]] 
end

```

Stack trace:

```julia
ERROR: MethodError: no method matching length(::Module)
Closest candidates are:
  length(::Union{Base.KeySet, Base.ValueIterator}) at abstractdict.jl:58
  length(::Union{HDF5.Attribute, HDF5.Dataset}) at ~/.julia/packages/HDF5/HtnQZ/src/dataspaces.jl:289
  length(::Union{HDF5.File, HDF5.Group}) at ~/.julia/packages/HDF5/HtnQZ/src/groups.jl:70
  ...
Stacktrace:
  [1] _similar_shape(itr::Module, #unused#::Base.HasLength)
    @ Base ./array.jl:663
  [2] _collect(cont::UnitRange{Int64}, itr::Module, #unused#::Base.HasEltype, isz::Base.HasLength)
    @ Base ./array.jl:718
  [3] collect(itr::Module)
    @ Base ./array.jl:712
  [4] broadcastable(x::Module)
    @ Base.Broadcast ./broadcast.jl:704
  [5] broadcasted(::Function, ::Module, ::Float64)
    @ Base.Broadcast ./broadcast.jl:1301
  [6] sys_tHE_AB!(F::Vector{Float64}, Z_AB::Vector{Float64}, k_e::Float64, fe::Float64, Fe::Float64, param::Parameter, j::Int64)
    @ Main ~/Dropbox/CORINNE_BU_RA/Julia/calibration/sys_tHE_AB.jl:15
  [7] (::var"#f_thresholds_HE!#11"{Int64, Int64, Int64})(F::Vector{Float64}, Z_AB::Vector{Float64})
    @ Main ~/Dropbox/CORINNE_BU_RA/Julia/calibration/1_calibrate_fe_US.jl:114
  [8] (::NLSolversBase.var"#fj_finitediff!#21"{var"#f_thresholds_HE!#11"{Int64, Int64, Int64}, FiniteDiff.JacobianCache{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, UnitRange{Int64}, Nothing, Val{:central}(), Float64}})(F::Vector{Float64}, J::Matrix{Float64}, x::Vector{Float64})
    @ NLSolversBase ~/.julia/packages/NLSolversBase/kavn7/src/objective_types/oncedifferentiable.jl:138
  [9] value_jacobian!!(obj::OnceDifferentiable{Vector{Float64}, Matrix{Float64}, Vector{Float64}}, F::Vector{Float64}, J::Matrix{Float64}, x::Vector{Float64})
    @ NLSolversBase ~/.julia/packages/NLSolversBase/kavn7/src/interface.jl:124
 [10] value_jacobian!!
    @ ~/.julia/packages/NLSolversBase/kavn7/src/interface.jl:122 [inlined]
 [11] trust_region_(df::OnceDifferentiable{Vector{Float64}, Matrix{Float64}, Vector{Float64}}, initial_x::Vector{Float64}, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, factor::Float64, autoscale::Bool, cache::NLsolve.NewtonTrustRegionCache{Vector{Float64}})
    @ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/solvers/trust_region.jl:119
 [12] trust_region (repeats 2 times)
    @ ~/.julia/packages/NLsolve/gJL1I/src/solvers/trust_region.jl:235 [inlined]
 [13] nlsolve(df::OnceDifferentiable{Vector{Float64}, Matrix{Float64}, Vector{Float64}}, initial_x::Vector{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
 [14] nlsolve(f::Function, initial_x::Vector{Float64}; method::Symbol, autodiff::Symbol, inplace::Bool, kwargs::Base.Pairs{Symbol, Real, NTuple{5, Symbol}, NamedTuple{(:xtol, :ftol, :iterations, :store_trace, :extended_trace), Tuple{Float64, Float64, Int64, Bool, Bool}}})
    @ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/nlsolve/nlsolve.jl:52
 [15] top-level scope
    @ ~/Dropbox/CORINNE_BU_RA/Julia/calibration/1_calibrate_fe_US.jl:115

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 3, 2023, 1:26pm UTC](https://discourse.julialang.org/t/domain-issue-when-working-with-nlsolve/98193/8 "2023-05-03T13:26:56Z")

</div>

Just do `const ^ = NanMath.^` before the function definition.

---

<div class="post-metadata">

**Author:** ![cdestael](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cdestael/32/46700_2.png) [@cdestael](https://discourse.julialang.org/u/cdestael)\
**Post date:** [May 3, 2023, 2:10pm UTC](https://discourse.julialang.org/t/domain-issue-when-working-with-nlsolve/98193/9 "2023-05-03T14:10:38Z")

</div>

Sorry for being dense about this @ChrisRackauckas, but when I put `const ^ = NanMath .^` 1) before the function definition, or 2) in the function, before `F[1]=...`, or 3) in the main file beginning, I get the following errors. Do I need to define it as a function?

1. 

```julia
julia> const ^ = NanMath .^
ERROR: syntax: incomplete: premature end of input
Stacktrace:
 [1] top-level scope
   @ ~/Dropbox/CORINNE_BU_RA/Julia/calibration/1_calibrate_fe_US.jl:5

```

1. 

```julia
ERROR: syntax: invalid function name ".^" around /Users/corinnes/Dropbox/CORINNE_BU_RA/Julia/calibration/sys_tHE_AB.jl:15
Stacktrace:
 [1] top-level scope
   @ ~/Dropbox/CORINNE_BU_RA/Julia/calibration/sys_tHE_AB.jl:2

```

1. 

```julia
julia> # Description: this file calibrates trade costs from each country to US
       using MAT, Statistics, DataFrames, LinearAlgebra
       using NLsolve
       using NaNMath
       nm = NaNMath
       const ^ = NanMath .^
       include("/Users/corinnes/Dropbox/CORINNE_BU_RA/Julia/calibration/sys_tHE_AB.jl")
ERROR: cannot assign a value to variable Base.^ from module Main
Stacktrace:
 [1] top-level scope
   @ ~/Dropbox/CORINNE_BU_RA/Julia/calibration/1_calibrate_fe_US.jl:6

```

---

<div class="post-metadata">

**Author:** ![stillyslalom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stillyslalom/32/45687_2.png) [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)\
**Post date:** [May 3, 2023, 2:51pm UTC](https://discourse.julialang.org/t/domain-issue-when-working-with-nlsolve/98193/10 "2023-05-03T14:51:27Z")

</div>

You have a space between `NaNMath` and `.^`. As Chris said, it should be `NanMath.^`.

---

<div class="post-metadata">

**Author:** ![cdestael](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cdestael/32/46700_2.png) [@cdestael](https://discourse.julialang.org/u/cdestael)\
**Post date:** [May 3, 2023, 7:20pm UTC](https://discourse.julialang.org/t/domain-issue-when-working-with-nlsolve/98193/11 "2023-05-03T19:20:51Z")

</div>

Apologies, that was a typo when I was writing about it here. However even without a space yields same errors (depending on placement):

```julia
julia> const ^ = NanMath.^
ERROR: syntax: incomplete: premature end of input
Stacktrace:
 [1] top-level scope
   @ ~/Dropbox/CORINNE_BU_RA/Julia/calibration/1_calibrate_fe_US.jl:5

```

or (placed in file before function definition)

```julia
julia> using NaNMath
       include("/Users/corinnes/Dropbox/CORINNE_BU_RA/Julia/calibration/sys_tHE_AB.jl")
ERROR: LoadError: UndefVarError: ^ not defined
Stacktrace:
 [1] top-level scope
   @ ~/Dropbox/CORINNE_BU_RA/Julia/calibration/sys_tHE_AB.jl:1
in expression starting at /Users/corinnes/Dropbox/CORINNE_BU_RA/Julia/calibration/sys_tHE_AB.jl:1

```

Any thoughts? Thanks!

---

<div class="post-metadata">

**Author:** ![DanDeepPhase](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dandeepphase/32/24076_2.png) [@DanDeepPhase](https://discourse.julialang.org/u/DanDeepPhase)\
**Post date:** [May 3, 2023, 8:03pm UTC](https://discourse.julialang.org/t/domain-issue-when-working-with-nlsolve/98193/12 "2023-05-03T20:03:19Z")

</div>

You have

`using NaNMath`  
`const ^ = NanMath.^`

where the second has a lower case n. The package name should be the same in both lines.

---

<div class="post-metadata">

**Author:** ![cdestael](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cdestael/32/46700_2.png) [@cdestael](https://discourse.julialang.org/u/cdestael)\
**Post date:** [May 3, 2023, 8:17pm UTC](https://discourse.julialang.org/t/domain-issue-when-working-with-nlsolve/98193/13 "2023-05-03T20:17:03Z")

</div>

Thanks @DanDeepPhase! The error persists with the correct spelling (and no space):  
`const ^ = NaNMath.^`.

```julia
ERROR: UndefVarError: ^ not defined
Stacktrace:
 [1] top-level scope
   @ ~/Dropbox/CORINNE_BU_RA/Julia/calibration/1_calibrate_fe_US.jl:5

```

Or:

```julia
ERROR: syntax: invalid function name ".^" around /Users/corinnes/Dropbox/CORINNE_BU_RA/Julia/calibration/sys_tHE_AB.jl:3
Stacktrace:
 [1] top-level scope
   @ ~/Dropbox/CORINNE_BU_RA/Julia/calibration/sys_tHE_AB.jl:2

```

---

<div class="post-metadata">

**Author:** ![DanDeepPhase](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dandeepphase/32/24076_2.png) [@DanDeepPhase](https://discourse.julialang.org/u/DanDeepPhase)\
**Post date:** [May 3, 2023, 9:52pm UTC](https://discourse.julialang.org/t/domain-issue-when-working-with-nlsolve/98193/15 "2023-05-03T21:52:24Z")

</div>

NaNMath doesn’t have a `^` function, but it does have a `pow` function.

Maybe try: `const ^ = NaNMath.pow`

debugging the debugging…Maybe that gets you to the fail state Chris was proposing.

---

<div class="post-metadata">

**Author:** ![cdestael](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cdestael/32/46700_2.png) [@cdestael](https://discourse.julialang.org/u/cdestael)\
**Post date:** [May 3, 2023, 10:24pm UTC](https://discourse.julialang.org/t/domain-issue-when-working-with-nlsolve/98193/16 "2023-05-03T22:24:46Z")

</div>

Thanks a lot for the followup, @DanDeepPhase!! I’m not having much success with `const ^ = NaNMath.pow` either (see error below). However now I’ve edited the code to use `max(var,0)` as a workaround (based on this post: [DifferentialEquations.jl, Error: Exponentiation yielding a complex result requires a complex argument - #6 by Gregers](https://discourse.julialang.org/t/differentialequations-jl-error-exponentiation-yielding-a-complex-result-requires-a-complex-argument/44935/6)).

```julia
ERROR: cannot assign a value to variable Base.^ from module Main
Stacktrace:
 [1] top-level scope
   @ ~/Dropbox/CORINNE_BU_RA/Julia/calibration/1_calibrate_fe_US.jl:5

```

Updated function:

```julia
function sys_tHE_AB!(F, Z_AB, k_e, fe, Fe, ρ, δ, μ, σ, μj, σj, j) # k_e: (1000, 10, 11) fe: (11, 11) Fe: (11, 11)
    
    ρ = 0.05 # rate of time preference
    δ = 0.05 # exogenous death rate (here also endogenous exit)
    μ=matread("parameters.mat")["mu"]
    σ=matread("parameters.mat")["sigma"]
    μj = μ[j]
    σj = σ[j]
    β = sqrt(2*(ρ + δ)/(σj^2) + (μj/(σj^2) - 1/2)^2) - (μj/(σj^2) - 1/2)
    α = -sqrt(2*(ρ + δ)/(σj^2) + (μj/(σj^2) - 1/2)^2) - (μj/(σj^2) - 1/2)
    disc = ρ + δ - μj

    F[1] = Z_AB[3]*max(Z_AB[1],0)^α + k_e*Z_AB[1]/disc - fe/disc - Fe - Z_AB[4]*max(Z_AB[1],0)^β
    F[2] = Z_AB[3]*max(Z_AB[2],0)^α + k_e*Z_AB[2]/disc - fe/disc - Z_AB[4]*max(Z_AB[2],0)^β
    F[3] = α*Z_AB[3]*max(Z_AB[1],0)^(α-1) + k_e/disc - β*Z_AB[4]*max(Z_AB[1],0)^(β-1)
    F[4] = α*Z_AB[3]*max(Z_AB[2],0)^(α-1) + k_e/disc - β*Z_AB[4]*max(Z_AB[2],0)^(β-1)

    N = [F[1] F[2] F[3] F[4]] 

end

```

---

<div class="post-metadata">

**Author:** ![sijo](https://avatars.discourse-cdn.com/v4/letter/s/da6949/32.png) [@sijo](https://discourse.julialang.org/u/sijo)\
**Post date:** [May 16, 2023, 11:57am UTC](https://discourse.julialang.org/t/domain-issue-when-working-with-nlsolve/98193/17 "2023-05-16T11:57:36Z")

</div>

I see you found another solution but for future reference:

The error `cannot assign a value to variable Base.^ from module Main` occurs when you try to define `^` after you have already used the standard `^` in Main. To make it work, restart Julia and make sure you define ` const ^ = NaNMath.pow` before using `^`.

---

<div class="post-metadata">

**Author:** ![cdestael](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cdestael/32/46700_2.png) [@cdestael](https://discourse.julialang.org/u/cdestael)\
**Post date:** [May 16, 2023, 1:58pm UTC](https://discourse.julialang.org/t/domain-issue-when-working-with-nlsolve/98193/18 "2023-05-16T13:58:27Z")

</div>

That is incredibly helpful @sijo – thank you!! This works, and likely more “correct” that using the max operator. Thanks again.