# Domain Troubles with Differential Equations Solver (with non-auto-differentiable function inside the ODE def)

**URL:** <https://discourse.julialang.org/t/domain-troubles-with-differential-equations-solver-with-non-auto-differentiable-function-inside-the-ode-def/13650>\
**Category:** Numerics\
**Tags:** question, diffeq\
**Created:** [August 18, 2018, 5:06am UTC](https://discourse.julialang.org/t/domain-troubles-with-differential-equations-solver-with-non-auto-differentiable-function-inside-the-ode-def/13650 "2018-08-18T05:06:31Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![shoshievass](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shoshievass/32/3260_2.png) [@shoshievass](https://discourse.julialang.org/u/shoshievass)\
**Post date:** [August 18, 2018, 5:06am UTC](https://discourse.julialang.org/t/domain-troubles-with-differential-equations-solver-with-non-auto-differentiable-function-inside-the-ode-def/13650/1 "2018-08-18T05:06:31Z")

</div>

I’m trying to solve an ODE that uses a non-auto-differentiable function inside its definition (in particular, I solve a constrained quadratic program with JuMP to define `du/dt` in each iteration). In principle, the solution to the ODE should be bounded between the 0 and the boundary condition (the first element of t is the largest t in the domain, and the last is the smallest)

If I pretend that the quadratic program inside is not constrained, I can write down an analytical solution for `du/dt` and solving the ode with a stiff solver works almost perfectly, without any explicit domain constraints.

With the constraint (and the JuMP model), however, using the `alg_hints=[:stiff]` option gives me a ForwardDiff error, which I take to mean that the stiff solver requires some automatic-derivatives that it can’t obtain.

Without the `stiff` option, I noticed that the solver kept taking `u` into negative numbers, where the ode didn’t make sense. I tried to restrict `u` to a positive domain in two ways:

(1) I added the callback `cb = PositiveDomain([1.0])`  
(2) I added the out\_of\_domain option: `isoutofdomain=(u,p,t) -> any(x -> (x < 0), u)`

For some reason that I don’t understand, using only the callback, didn’t actually keep `u` from going negative, but using only the `isoutofdomain` function produced other errors. Using both options together made it work sometimes - and it looked correct when it worked.

In the majority of cases, though, the solver still kept breaking, and I noticed that the `u` values were getting shot up way above the boundary value just before the solver broke.

I tried to solve that by expanding the `isoutofdomain` function to: `isoutofdomain=(u,p,t) -> any(x -> ((x < 0) | (x > s_max)), u)`.

This didn’t seem to help - `u` kept getting shot up and then turning to NaN:

```julia
u: 1156.621761774482t: 1.113545608146277
u: 1155.2373082699469t: 1.1123344073491215
u: 1157.9898666009408t: 1.1081535756818317
u: 1225.539274307506t: 1.1075696780978779
u: 6.158591645740778e19t: 1.1074239366474008
u: NaNt: 1.1074239366474008
...some warnings...
Solution interpolation cannot extrapolate past the final timepoint. Either solve on a longer timespan or use the local extrapolation from the integrator interface.

```

It’s a bit hard to reproduce a complete example that breaks (I tried to make a baby example but it worked fine…) but this is the main idea [I can produce a fuller example if need be]

```julia
function du_dt( t, u, args)

    b_opt = getBidVector( args ) # This calls a JuMP model that solves a quadratic program and outputs a vector b_opt
    db_du = getDbDsVec(args, b_opt) # This calls a simple (analytically-defined) function using args and b_opt

    println("u: ", u, "t: ", t)

    denominator = db_du' * ( f(t, u, args, b_opt ) ); #f(.) is also analytically-defined
    numerator = g(t, u, args, b_opt) # some function 

    du_dt = numerator / denominator; # Note: this should be a scalar

    return du_dt

end

and then:

ode(u, p, t) = du_dt(t, u, args)

prob = ODEProblem(ode, u_max, [t_max : 20 : t_min])

cb = PositiveDomain([1.0])

sol = DifferentialEquations.solve(prob, callback=cb, isoutofdomain=(u,p,t) -> any(x -> ((x < 0) | (x > u_max)), u))

```

Might anyone have any idea of what might be going wrong and/or how to fix it? It would be very much appreciated! Perhaps you might know off-hand @ChrisRackauckas ?

---

<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:** [August 18, 2018, 11:47am UTC](https://discourse.julialang.org/t/domain-troubles-with-differential-equations-solver-with-non-auto-differentiable-function-inside-the-ode-def/13650/2 "2018-08-18T11:47:01Z")

</div>

Is it non-differentiable or non-auto-differentiable? If it’s non-differentiable you need to be careful with what algorithms you chose (and about existence of a solution).

What happened with only `isoutofdomain`? I would try only `isoutofdomain` since I think there may be an issue with the `PositiveDomain` callback but I’m not certain.

Is the JuMP model converging to a unique solution to high enough tolerance every time?

What happens inside of the step which suddenly jumps up? What value of the derivative is large or changing? I cannot check this without a runnable example, but you can print stuff out.

Also,

> [@shoshievass](#):
>
> ode(u, p, t) = du\_dt(t, u, args)

that’s a little odd: why not `ode(u,args,t)`?

Hopefully these narrow down the problem.

---

<div class="post-metadata">

**Author:** ![shoshievass](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shoshievass/32/3260_2.png) [@shoshievass](https://discourse.julialang.org/u/shoshievass)\
**Post date:** [August 18, 2018, 12:47pm UTC](https://discourse.julialang.org/t/domain-troubles-with-differential-equations-solver-with-non-auto-differentiable-function-inside-the-ode-def/13650/3 "2018-08-18T12:47:42Z")

</div>

Thank you very much for the response!

It is non-_auto_-differentiable. It is analytically differentiable inside the (positive) domain and I compute that derivative in `db_du` inside of `du_dt`, but I’m not sure how to pass that information on further to the integrator…

The JuMP model converges most of the time - there are occasionally blips but the integrator seems to do ok with those, and the jumps to NaN don’t seem to coincide…

I tried taking out the PositiveDomain() callback (but keeping the isoutofdomain option) and for some reason, now, certain cases that worked before (w/ the cb) had the same problem as the other ones…

With some additional print statements:

```julia
u: 1126.407942124624; t: 1.0849745558374158
numerator: 409.7772562472494
denominator: 0.6585804546718615
du/dt: 622.2129025244476

u: 1126.5099110176955; t: 1.0849745558374158
numerator: 439.4244401845725
denominator: 0.6584945168617354
du/dt: 667.3167793086404

u: 1125.7547334090814; t: 1.083842892986983
numerator: 584.4985936643423
denominator: 0.6581408156313718
du/dt: 888.1056755363476

u: 1124.4554650417726; t: 1.0826760853275306
numerator: 550.4443067543003
denominator: 0.6582149155978451
du/dt: 836.26835811574

u: 1126.9772934735458; t: 1.0786484902138787
numerator: 46014.01725027083
denominator: 0.6442312490613255
du/dt: 71424.68999651192

u: 1177.2474896547437; t: 1.078085993743522
numerator: 4.0819649439384765e18
denominator: 0.6012901278348439
du/dt: 6.788677802905269e18

u: 1.348924369022055e15; t: 1.0779455940334857
numerator: NaN
denominator: NaN
du/dt: NaN

u: NaN; t: 1.0779455940334857
numerator: 6.209605055194862e120
denominator: -2.52642256119077e-10
du/dt: -2.457864788963929e130

```

So `du/dt` is jumping around a lot, mostly because the `numerator` function jumps around - this is because it includes an `exp(.)` function, which sadly makes the ODE harder to stably solve. I guess what I’d like to tell the solver to do is to roll back if it sees crazy values of `u` or `du/dt` like it does before it crashes here. That’s what I was hoping to do with `isoutofdomain=(u,p,t) -> any(x -> ((x < 0) | (x > s_max)), u)`, but maybe that gets overridden?

I don’t have a deep understanding, but maybe this is why stiff solvers do better on the unconstrained version…  
Is there anyway to use a stiff solver here? This is the error I get when I try:

```julia
MethodError: Cannot `convert` an object of type ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.TimeDerivativeWrapper{#ode_multi_item_w_risk_aversion_subs#215{Auction,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Float64,Void},Float64},Float64,1} to an object of type Float64
This may have arisen from a call to the constructor Float64(...),
since type constructors fall back to convert methods.

```

Re `ode(u, p, t) = du_dt(t, u, args)` - I was just using this as a shortcut (alternate function definition) since the args don’t change between calls and putting in a generic value of 1 for `p`; not idiomatic 😅 but I’ve used similar things with other ODEs that seem to work fine.

Thank you very much once again!

---

<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:** [August 18, 2018, 12:58pm UTC](https://discourse.julialang.org/t/domain-troubles-with-differential-equations-solver-with-non-auto-differentiable-function-inside-the-ode-def/13650/4 "2018-08-18T12:58:06Z")

</div>

> [@shoshievass](#):
>
> I don’t have a deep understanding, but maybe this is why stiff solvers do better on the unconstrained version…  
> Is there anyway to use a stiff solver here?

Yes, just tell it to not autodiff. `Rosenbrock23(autodiff=false)` or `Rodas5(autodiff=false)`.

> [@shoshievass](#):
>
> The JuMP model converges most of the time - there are occasionally blips but the integrator seems to do ok with those, and the jumps to NaN don’t seem to coincide…

This could be the issue. Is it used in the numerator? You need to find out what suddenly gives `numerator: 46014.01725027083` and whether that is a fluke. A method for stiff equations can handle it if it’s not a fluke, but if it is a fluke then that looks like the start of some odd divergence.

> [@shoshievass](#):
>
> So `du/dt` is jumping around a lot, mostly because the `numerator` function jumps around - this is because it includes an `exp(.)` function, which sadly makes the ODE harder to stably solve. I guess what I’d like to tell the solver to do is to roll back if it sees crazy values of `u` or `du/dt` like it does before it crashes here. That’s what I was hoping to do with `isoutofdomain=(u,p,t) -> any(x -> ((x < 0) | (x > s_max)), u)` , but maybe that gets overridden?

It would only be overridden if your derivative function wasn’t deterministic. The fact that the JuMP model can diverge could be making it non-deterministic and then the derivative is 🤷‍♂️.

---

<div class="post-metadata">

**Author:** ![shoshievass](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shoshievass/32/3260_2.png) [@shoshievass](https://discourse.julialang.org/u/shoshievass)\
**Post date:** [August 18, 2018, 7:27pm UTC](https://discourse.julialang.org/t/domain-troubles-with-differential-equations-solver-with-non-auto-differentiable-function-inside-the-ode-def/13650/5 "2018-08-18T19:27:24Z")

</div>

> [@ChrisRackauckas](#):
>
> Yes, just tell it to not autodiff. `Rosenbrock23(autodiff=false)` or `Rodas5(autodiff=false)` .

Is this enough? I tried changing my ODE call to:

```julia
sol = DifferentialEquations.solve(prob, Rosenbrock23(autodiff=false)) 

```

for example, and still got the following error (I’ll copy the whole thing this time):

```julia
u: 1251.8518364213967; t: 1.2102828512236605
numerator: -1.176735245048757e-12
denominator: 0.654046236490126
du/dt: -1.7991621683561538e-12

u: 1251.8518364213967; t: 1.2102828512236605
numerator: -1.176735245048757e-12
denominator: 0.654046236490126
du/dt: -1.7991621683561538e-12

u: 1251.8518364213967; t: 1.2102818512236606
numerator: 0.01367538009225448
denominator: 0.6540454086369006
du/dt: 0.020908915362245886

MethodError: Cannot `convert` an object of type ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.TimeDerivativeWrapper{#ode_multi_item_w_risk_aversion_subs#313{Auction,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Float64,Void},Float64},Float64,1} to an object of type Float64
This may have arisen from a call to the constructor Float64(...),
since type constructors fall back to convert methods.

Stacktrace:
 [1] push!(::Array{Float64,1}, ::ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.TimeDerivativeWrapper{#ode_multi_item_w_risk_aversion_subs#313{Auction,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Float64,Void},Float64},Float64,1}) at ./array.jl:646
 [2] parseNLExpr_runtime(::JuMP.Model, ::ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.TimeDerivativeWrapper{#ode_multi_item_w_risk_aversion_subs#313{Auction,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Float64,Void},Float64},Float64,1}, ::Array{ReverseDiffSparse.NodeData,1}, ::Int64, ::Array{Float64,1}) at /Users/shoshievass/.julia/v0.6/JuMP/src/parsenlp.jl:196
 [3] macro expansion at /Users/shoshievass/.julia/v0.6/JuMP/src/parseExpr_staged.jl:489 [inlined]
 [4] macro expansion at /Users/shoshievass/.julia/v0.6/JuMP/src/parsenlp.jl:226 [inlined]
 [5] macro expansion at /Users/shoshievass/.julia/v0.6/JuMP/src/macros.jl:1157 [inlined]
 [6] getBidVector(::Float64, ::ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.TimeDerivativeWrapper{#ode_multi_item_w_risk_aversion_subs#313{Auction,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Float64,Void},Float64},Float64,1}, ::Auction, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}) at /Users/shoshievass/Git/Procurement_Risk_Aversion/julia_cost_estimation/counterfactual/aggregatedDataCF/quadprogCounterfactualFunctions.jl:110
 [7] ode_multi_item_w_risk_aversion(::Auction, ::ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.TimeDerivativeWrapper{#ode_multi_item_w_risk_aversion_subs#313{Auction,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Float64,Void},Float64},Float64,1}, ::Float64, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}) at /Users/shoshievass/Git/Procurement_Risk_Aversion/julia_cost_estimation/counterfactual/aggregatedDataCF/quadprogCounterfactualFunctions.jl:167
 [8] perform_step!(::OrdinaryDiffEq.ODEIntegrator{OrdinaryDiffEq.Rosenbrock23{0,false,DiffEqBase.LinSolveFactorize{Base.LinAlg.#lufact!},DataType},Float64,Float64,Void,Float64,Float64,Float64,Array{Float64,1},DiffEqBase.ODESolution{Float64,1,Array{Float64,1},Void,Void,Array{Float64,1},Array{Array{Float64,1},1},DiffEqBase.ODEProblem{Float64,Float64,false,Void,#ode_multi_item_w_risk_aversion_subs#313{Auction,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Void,Void,UniformScaling{Int64},DiffEqBase.StandardODEProblem},OrdinaryDiffEq.Rosenbrock23{0,false,DiffEqBase.LinSolveFactorize{Base.LinAlg.#lufact!},DataType},OrdinaryDiffEq.InterpolationData{#ode_multi_item_w_risk_aversion_subs#313{Auction,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Array{Float64,1},Array{Float64,1},Array{Array{Float64,1},1},OrdinaryDiffEq.Rosenbrock23ConstantCache{Float64,DiffEqDiffTools.TimeDerivativeWrapper{#ode_multi_item_w_risk_aversion_subs#313{Auction,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Float64,Void},DiffEqDiffTools.UDerivativeWrapper{#ode_multi_item_w_risk_aversion_subs#313{Auction,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Float64,Void}}}},#ode_multi_item_w_risk_aversion_subs#313{Auction,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Void,OrdinaryDiffEq.Rosenbrock23ConstantCache{Float64,DiffEqDiffTools.TimeDerivativeWrapper{#ode_multi_item_w_risk_aversion_subs#313{Auction,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Float64,Void},DiffEqDiffTools.UDerivativeWrapper{#ode_multi_item_w_risk_aversion_subs#313{Auction,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Float64,Void}},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,DiffEqBase.#ODE_DEFAULT_NORM,DiffEqBase.CallbackSet{Tuple{},Tuple{}},DiffEqBase.#ODE_DEFAULT_ISOUTOFDOMAIN,DiffEqBase.#ODE_DEFAULT_PROG_MESSAGE,DiffEqBase.#ODE_DEFAULT_UNSTABLE_CHECK,DataStructures.BinaryHeap{Float64,DataStructures.GreaterThan},DataStructures.BinaryHeap{Float64,DataStructures.GreaterThan},Void,Void,Int64,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Float64}, ::OrdinaryDiffEq.Rosenbrock23ConstantCache{Float64,DiffEqDiffTools.TimeDerivativeWrapper{#ode_multi_item_w_risk_aversion_subs#313{Auction,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Float64,Void},DiffEqDiffTools.UDerivativeWrapper{#ode_multi_item_w_risk_aversion_subs#313{Auction,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Float64,Void}}, ::Bool) at /Users/shoshievass/.julia/v0.6/OrdinaryDiffEq/src/perform_step/rosenbrock_perform_step.jl:175
 [9] solve!(::OrdinaryDiffEq.ODEIntegrator{OrdinaryDiffEq.Rosenbrock23{0,false,DiffEqBase.LinSolveFactorize{Base.LinAlg.#lufact!},DataType},Float64,Float64,Void,Float64,Float64,Float64,Array{Float64,1},DiffEqBase.ODESolution{Float64,1,Array{Float64,1},Void,Void,Array{Float64,1},Array{Array{Float64,1},1},DiffEqBase.ODEProblem{Float64,Float64,false,Void,#ode_multi_item_w_risk_aversion_subs#313{Auction,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Void,Void,UniformScaling{Int64},DiffEqBase.StandardODEProblem},OrdinaryDiffEq.Rosenbrock23{0,false,DiffEqBase.LinSolveFactorize{Base.LinAlg.#lufact!},DataType},OrdinaryDiffEq.InterpolationData{#ode_multi_item_w_risk_aversion_subs#313{Auction,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Array{Float64,1},Array{Float64,1},Array{Array{Float64,1},1},OrdinaryDiffEq.Rosenbrock23ConstantCache{Float64,DiffEqDiffTools.TimeDerivativeWrapper{#ode_multi_item_w_risk_aversion_subs#313{Auction,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Float64,Void},DiffEqDiffTools.UDerivativeWrapper{#ode_multi_item_w_risk_aversion_subs#313{Auction,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Float64,Void}}}},#ode_multi_item_w_risk_aversion_subs#313{Auction,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Void,OrdinaryDiffEq.Rosenbrock23ConstantCache{Float64,DiffEqDiffTools.TimeDerivativeWrapper{#ode_multi_item_w_risk_aversion_subs#313{Auction,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Float64,Void},DiffEqDiffTools.UDerivativeWrapper{#ode_multi_item_w_risk_aversion_subs#313{Auction,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Float64,Void}},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,DiffEqBase.#ODE_DEFAULT_NORM,DiffEqBase.CallbackSet{Tuple{},Tuple{}},DiffEqBase.#ODE_DEFAULT_ISOUTOFDOMAIN,DiffEqBase.#ODE_DEFAULT_PROG_MESSAGE,DiffEqBase.#ODE_DEFAULT_UNSTABLE_CHECK,DataStructures.BinaryHeap{Float64,DataStructures.GreaterThan},DataStructures.BinaryHeap{Float64,DataStructures.GreaterThan},Void,Void,Int64,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Float64}) at /Users/shoshievass/.julia/v0.6/OrdinaryDiffEq/src/solve.jl:341
 [10] #solve#1469(::Array{Any,1}, ::Function, ::DiffEqBase.ODEProblem{Float64,Float64,false,Void,#ode_multi_item_w_risk_aversion_subs#313{Auction,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Void,Void,UniformScaling{Int64},DiffEqBase.StandardODEProblem}, ::OrdinaryDiffEq.Rosenbrock23{0,false,DiffEqBase.LinSolveFactorize{Base.LinAlg.#lufact!},DataType}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /Users/shoshievass/.julia/v0.6/OrdinaryDiffEq/src/solve.jl:7
 [11] solve(::DiffEqBase.ODEProblem{Float64,Float64,false,Void,#ode_multi_item_w_risk_aversion_subs#313{Auction,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Void,Void,UniformScaling{Int64},DiffEqBase.StandardODEProblem}, ::OrdinaryDiffEq.Rosenbrock23{0,false,DiffEqBase.LinSolveFactorize{Base.LinAlg.#lufact!},DataType}) at /Users/shoshievass/.julia/v0.6/OrdinaryDiffEq/src/solve.jl:6
 [12] getScalingAuctionEquilibrium(::Int64, ::Auction, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}) at /Users/shoshievass/Git/Procurement_Risk_Aversion/julia_cost_estimation/counterfactual/aggregatedDataCF/quadprogCounterfactualFunctions.jl:217
 [13] getMinimalAuctionCFcomparison(::Auction, ::Int64, ::Int64) at /Users/shoshievass/Git/Procurement_Risk_Aversion/julia_cost_estimation/counterfactual/aggregatedDataCF/quadprogCounterfactualFunctions.jl:240
 [14] include_string(::String, ::String) at ./loading.jl:522

```

I’ll try to figure out the JuMP thing - I think the problem should always have a unique solution so this shouldn’t be the bottleneck. Being able to use a stiff solver seems like it might help a lot though…

---

<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:** [August 18, 2018, 7:35pm UTC](https://discourse.julialang.org/t/domain-troubles-with-differential-equations-solver-with-non-auto-differentiable-function-inside-the-ode-def/13650/6 "2018-08-18T19:35:07Z")

</div>

Can you post your `Pkg.status()`? That shouldn’t have ForwardDiff involved at all.

---

<div class="post-metadata">

**Author:** ![shoshievass](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shoshievass/32/3260_2.png) [@shoshievass](https://discourse.julialang.org/u/shoshievass)\
**Post date:** [August 18, 2018, 7:36pm UTC](https://discourse.julialang.org/t/domain-troubles-with-differential-equations-solver-with-non-auto-differentiable-function-inside-the-ode-def/13650/7 "2018-08-18T19:36:41Z")

</div>

> [@ChrisRackauckas](#):
>
> Pkg.status()

```julia
Pkg.status(): 

11 required packages:
 - CSV 0.2.5
 - DataFrames 0.11.7
 - DifferentialEquations 4.5.0
 - Distributions 0.15.0
 - Hiccup 0.1.1
 - IJulia 1.9.3
 - Ipopt 0.4.0
 - IterableTables 0.7.3
 - JuMP 0.18.2
 - NLsolve 1.0.1
 - Plots 0.17.4
109 additional packages:
 - BandedMatrices 0.4.0
 - BinDeps 0.8.10
 - BinaryProvider 0.3.3
 - BoundaryValueDiffEq 1.0.0
 - Calculus 0.4.1
 - CategoricalArrays 0.3.13
 - ChunkedArrays 0.1.1
 - CodecZlib 0.4.4
 - ColorTypes 0.6.7
 - Colors 0.8.2
 - CommonSubexpressions 0.1.0
 - Compat 1.0.1
 - Conda 1.0.1
 - Contour 0.4.0
 - DataStreams 0.3.6
 - DataStructures 0.8.4
 - DataValues 0.3.3
 - DelayDiffEq 3.6.0
 - DiffBase 0.3.2
 - DiffEqBase 3.13.3
 - DiffEqBiological 2.3.2
 - DiffEqCallbacks 1.1.1
 - DiffEqDevTools 1.0.2
 - DiffEqDiffTools 0.4.1
 - DiffEqFinancial 1.0.1
 - DiffEqJump 4.5.1
 - DiffEqMonteCarlo 0.10.1
 - DiffEqNoiseProcess 1.0.2
 - DiffEqOperators 1.3.0
 - DiffEqPDEBase 0.4.0
 - DiffEqParamEstim 1.1.2
 - DiffEqPhysics 1.0.0
 - DiffEqSensitivity 1.2.0
 - DiffEqUncertainty 0.1.0
 - DiffResults 0.0.3
 - DiffRules 0.0.7
 - DimensionalPlotRecipes 0.0.2
 - Distances 0.6.0
 - DistributedArrays 0.4.0
 - EllipsisNotation 0.3.0
 - FillArrays 0.1.0
 - FixedPointNumbers 0.4.6
 - ForwardDiff 0.7.5
 - FunctionWrappers 1.0.0
 - GR 0.32.3
 - GenericSVD 0.1.0
 - InternedStrings 0.6.2
 - IteratorInterfaceExtensions 0.0.2
 - JSON 0.17.2
 - Juno 0.4.1
 - LearnBase 0.1.6
 - LineSearches 3.2.5
 - LinearMaps 1.0.4
 - LossFunctions 0.2.0
 - LsqFit 0.3.0
 - MacroTools 0.4.4
 - MathOptInterface 0.4.1
 - MathProgBase 0.7.2
 - MbedTLS 0.5.12
 - Measures 0.2.0
 - Media 0.3.0
 - Missings 0.2.10
 - MuladdMacro 0.0.2
 - MultiScaleArrays 0.7.1
 - NLSolversBase 4.4.1
 - NaNMath 0.3.2
 - NamedTuples 4.0.2
 - Nullables 0.0.7
 - Optim 0.14.1
 - OptimBase 1.0.0
 - OrdinaryDiffEq 3.21.0
 - PDMats 0.8.0
 - ParameterizedFunctions 3.1.0
 - Parameters 0.9.2
 - PenaltyFunctions 0.0.2
 - PlotThemes 0.2.0
 - PlotUtils 0.4.4
 - PoissonRandom 0.0.1
 - PositiveFactorizations 0.1.0
 - Primes 0.3.0
 - QuadGK 0.3.0
 - RandomNumbers 0.1.1
 - RecipesBase 0.3.1
 - RecursiveArrayTools 0.15.0
 - Reexport 0.1.0
 - Requires 0.4.4
 - ResettableStacks 0.3.1
 - ReverseDiffSparse 0.8.2
 - Rmath 0.4.0
 - Roots 0.7.1
 - SHA 0.5.7
 - Showoff 0.2.1
 - SimpleTraits 0.6.0
 - SortingAlgorithms 0.2.1
 - SpecialFunctions 0.6.0
 - StaticArrays 0.7.2
 - StatsBase 0.23.1
 - StatsFuns 0.6.1
 - SteadyStateDiffEq 0.4.0
 - StochasticDiffEq 4.4.5
 - Sundials 1.6.0
 - SymEngine 0.4.2
 - TableTraits 0.2.0
 - TranscodingStreams 0.5.4
 - URIParser 0.3.1
 - VectorizedRoutines 0.0.2
 - VersionParsing 1.1.2
 - WeakRefStrings 0.4.7
 - ZMQ 0.6.4

```

---

<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:** [August 18, 2018, 7:40pm UTC](https://discourse.julialang.org/t/domain-troubles-with-differential-equations-solver-with-non-auto-differentiable-function-inside-the-ode-def/13650/8 "2018-08-18T19:40:12Z")

</div>

Oh, you’re running into [https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/259](https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/259) . You need to use an inplace function to use finite differencing for now. It’s one of the limitations we are working on.

```julia
ode(du, u, p, t) = (du .= du_dt(t, u, args))

```

is the trivial way to do it (though of course you can dig in and optimize it more)

---

<div class="post-metadata">

**Author:** ![shoshievass](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shoshievass/32/3260_2.png) [@shoshievass](https://discourse.julialang.org/u/shoshievass)\
**Post date:** [August 18, 2018, 7:51pm UTC](https://discourse.julialang.org/t/domain-troubles-with-differential-equations-solver-with-non-auto-differentiable-function-inside-the-ode-def/13650/9 "2018-08-18T19:51:25Z")

</div>

Ah! Cool 😃

Now I’m getting an `InexactError` but maybe that’s some issue on my end? 😕

```julia
InexactError()

Stacktrace:
 [1] convert(::Type{Int64}, ::Float64) at ./float.jl:679
 [2] zeros at ./array.jl:258 [inlined]
 [3] zeros at ./array.jl:259 [inlined]
 [4] zeros at ./array.jl:261 [inlined]
 [5] alg_cache(::OrdinaryDiffEq.Rodas5{0,false,DiffEqBase.LinSolveFactorize{Base.LinAlg.#lufact!},DataType}, ::Float64, ::Float64, ::Type{T} where T, ::Type{T} where T, ::Type{T} where T, ::Float64, ::Float64, ::#ode_multi_item_w_risk_aversion_subs#347{Auction,Array{Float64,1},Array{Float64,1},Array{Float64,1}}, ::Float64, ::Float64, ::Float64, ::Void, ::Bool, ::Type{Val{true}}) at /Users/shoshievass/.julia/v0.6/OrdinaryDiffEq/src/caches/rosenbrock_caches.jl:769
 [6] #init#1470(::Int64, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Void, ::Bool, ::Void, ::Bool, ::Bool, ::Void, ::Bool, ::Bool, ::Float64, ::Bool, ::Rational{Int64}, ::Void, ::Void, ::Int64, ::Rational{Int64}, ::Int64, ::Rational{Int64}, ::Rational{Int64}, ::Bool, ::Int64, ::Void, ::Void, ::Int64, ::Float64, ::Float64, ::DiffEqBase.#ODE_DEFAULT_NORM, ::##344#348, ::DiffEqBase.#ODE_DEFAULT_UNSTABLE_CHECK, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Int64, ::String, ::DiffEqBase.#ODE_DEFAULT_PROG_MESSAGE, ::Void, ::Bool, ::Bool, ::Array{Any,1}, ::DiffEqBase.#init, ::DiffEqBase.ODEProblem{Float64,Float64,true,Void,#ode_multi_item_w_risk_aversion_subs#347{Auction,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Void,Void,UniformScaling{Int64},DiffEqBase.StandardODEProblem}, ::OrdinaryDiffEq.Rodas5{0,false,DiffEqBase.LinSolveFactorize{Base.LinAlg.#lufact!},DataType}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /Users/shoshievass/.julia/v0.6/OrdinaryDiffEq/src/solve.jl:226
 [7] (::DiffEqBase.#kw##init)(::Array{Any,1}, ::DiffEqBase.#init, ::DiffEqBase.ODEProblem{Float64,Float64,true,Void,#ode_multi_item_w_risk_aversion_subs#347{Auction,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Void,Void,UniformScaling{Int64},DiffEqBase.StandardODEProblem}, ::OrdinaryDiffEq.Rodas5{0,false,DiffEqBase.LinSolveFactorize{Base.LinAlg.#lufact!},DataType}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at ./<missing>:0
 [8] #solve#1469(::Array{Any,1}, ::Function, ::DiffEqBase.ODEProblem{Float64,Float64,true,Void,#ode_multi_item_w_risk_aversion_subs#347{Auction,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Void,Void,UniformScaling{Int64},DiffEqBase.StandardODEProblem}, ::OrdinaryDiffEq.Rodas5{0,false,DiffEqBase.LinSolveFactorize{Base.LinAlg.#lufact!},DataType}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /Users/shoshievass/.julia/v0.6/OrdinaryDiffEq/src/solve.jl:6
 [9] (::DiffEqBase.#kw##solve)(::Array{Any,1}, ::DiffEqBase.#solve, ::DiffEqBase.ODEProblem{Float64,Float64,true,Void,#ode_multi_item_w_risk_aversion_subs#347{Auction,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Void,Void,UniformScaling{Int64},DiffEqBase.StandardODEProblem}, ::OrdinaryDiffEq.Rodas5{0,false,DiffEqBase.LinSolveFactorize{Base.LinAlg.#lufact!},DataType}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at ./<missing>:0 (repeats 2 times)
 [10] getScalingAuctionEquilibrium(::Int64, ::Auction, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}) at /Users/shoshievass/Git/Procurement_Risk_Aversion/julia_cost_estimation/counterfactual/aggregatedDataCF/quadprogCounterfactualFunctions.jl:219
 [11] getMinimalAuctionCFcomparison(::Auction, ::Int64, ::Int64) at /Users/shoshievass/Git/Procurement_Risk_Aversion/julia_cost_estimation/counterfactual/aggregatedDataCF/quadprogCounterfactualFunctions.jl:243
 [12] include_string(::String, ::String) at ./loading.jl:522

```

---

<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:** [August 18, 2018, 7:52pm UTC](https://discourse.julialang.org/t/domain-troubles-with-differential-equations-solver-with-non-auto-differentiable-function-inside-the-ode-def/13650/10 "2018-08-18T19:52:30Z")

</div>

Is your uType or tType an integer?

---

<div class="post-metadata">

**Author:** ![shoshievass](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shoshievass/32/3260_2.png) [@shoshievass](https://discourse.julialang.org/u/shoshievass)\
**Post date:** [August 18, 2018, 7:53pm UTC](https://discourse.julialang.org/t/domain-troubles-with-differential-equations-solver-with-non-auto-differentiable-function-inside-the-ode-def/13650/11 "2018-08-18T19:53:02Z")

</div>

No they should both be floats (e.g. `u: 1251.8518364213967; t: 1.2102828512236605`) 😕

---

<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:** [August 18, 2018, 7:54pm UTC](https://discourse.julialang.org/t/domain-troubles-with-differential-equations-solver-with-non-auto-differentiable-function-inside-the-ode-def/13650/12 "2018-08-18T19:54:47Z")

</div>

Can you show me what line is erroring? Without an example this is hard to debug, but if you post where in the code things are erroring I can try to see if I see anything.

---

<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:** [August 18, 2018, 7:56pm UTC](https://discourse.julialang.org/t/domain-troubles-with-differential-equations-solver-with-non-auto-differentiable-function-inside-the-ode-def/13650/13 "2018-08-18T19:56:32Z")

</div>

[https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/v3.21.0/src/caches/rosenbrock\_caches.jl#L769](https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/v3.21.0/src/caches/rosenbrock_caches.jl#L769)

It’s showing that `rate_prototype` is likely integers. What is your initial condition?

---

<div class="post-metadata">

**Author:** ![shoshievass](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shoshievass/32/3260_2.png) [@shoshievass](https://discourse.julialang.org/u/shoshievass)\
**Post date:** [August 18, 2018, 7:57pm UTC](https://discourse.julialang.org/t/domain-troubles-with-differential-equations-solver-with-non-auto-differentiable-function-inside-the-ode-def/13650/14 "2018-08-18T19:57:13Z")

</div>

Can you tell from the error where it’s trying to do the conversion? The reference to `getScalingAuctionEquilibrium` is just the outer function of the ODE call.

I’d be happy to share all the code and an example if it were helpful - I just don’t want to make you dig through my code unnecessarily.

---

<div class="post-metadata">

**Author:** ![shoshievass](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shoshievass/32/3260_2.png) [@shoshievass](https://discourse.julialang.org/u/shoshievass)\
**Post date:** [August 18, 2018, 7:58pm UTC](https://discourse.julialang.org/t/domain-troubles-with-differential-equations-solver-with-non-auto-differentiable-function-inside-the-ode-def/13650/15 "2018-08-18T19:58:38Z")

</div>

Initial conditions:

```julia
init u: 1251.8518364213967
t span: (1.2102828512236605, 0.6313105134942845)

```

---

<div class="post-metadata">

**Author:** ![shoshievass](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shoshievass/32/3260_2.png) [@shoshievass](https://discourse.julialang.org/u/shoshievass)\
**Post date:** [August 18, 2018, 8:29pm UTC](https://discourse.julialang.org/t/domain-troubles-with-differential-equations-solver-with-non-auto-differentiable-function-inside-the-ode-def/13650/16 "2018-08-18T20:29:17Z")

</div>

Okay, here is an example that worked with the callback and domain restriction but doesn’t work now. No worries if you’re too busy to look - I really, really appreciate the help!

```julia
using JuMP
using Ipopt
using DataFrames, CSV
using DifferentialEquations
using Distributions
using NLsolve

type Bidder
    bidder_id::String
    b::Vector{Float64}
    score::Float64
end

type Auction
    contract_no::Int64
    num_bidders::Int64
    T::Int64
    bidders::Array{Bidder}
    q_a::Vector{Float64}
    q_a_model::Vector{Float64}
    q_o::Vector{Float64}
    c::Vector{Float64}
    sigma_t_sq::Vector{Float64}
    gamma::Float64
    item_ids::Vector{Int64}
    optb_denom::Array{Float64}
    optb_num_coeff::Array{Float64}
    optb_num_main::Array{Float64}
    extra_work_payment::Float64
    winning_alpha::Float64
    highest_alpha::Float64
    lowest_alpha::Float64
    alpha_lgmean::Float64
    sigma_lgalpha::Float64
    dollar_scale::Int64
end

type AuctionEquilibrium
    auction_seqid::Int64
    contract_no::Int64
    alpha_types::Vector{Float64}
    score_maps::Vector{Float64}
    expost_cost::Vector{Float64}
end

type AugmentedAuctionEquilibrium
    auction_seqid::Int64
    contract_no::Int64
    alpha_types::Vector{Float64}
    score_maps::Vector{Float64}
    expost_cost::Vector{Float64}
    winning_eq_score::Float64
    winning_eq_cost::Float64
end

auc = Auction(30163,
	 		  6,
			  80, 
			  Bidder[
			  	Bidder("69", [3.0, 50.0, 7.5, 7.5, 1.8, 1.3, 0.25, 0.6, 1.8, 1.4, 5.0, 0.6, 5.0, 5.0, 2.55, 2.55, 6.15, 0.9, 1.6, 0.85, 1.2, 0.3, 0.36, 0.38, 0.65, 0.335, 6.0, 3.5, 3.75, 0.07, 0.8, 2.25, 4.5, 5.0, 0.04, 5.5, 0.95, 0.013, 0.3, 0.3, 0.66, 0.8, 0.18, 0.185, 0.17, 0.07, 0.66, 0.8, 0.49, 0.045, 0.17, 1.25, 1.25, 0.07, 3.1, 1.2, 3.0, 0.11, 0.9, 18.0, 4.5, 0.6, 0.55, 0.9, 0.3, 0.075, 1.45, 0.125, 7.5, 3.0, 0.5, 0.35, 0.14, 0.145, 0.15, 14.5, 3.5, 55.0, 145.0, 336.0], 951.971), 
				Bidder("108", [5.3, 70.0, 15.0, 7.0, 1.2, 1.075, 0.15, 0.95, 1.7, 3.0, 4.0, 1.0, 1.0, 2.0, 0.0001, 0.0001, 0.0001, 0.55, 1.1, 0.85, 0.95, 0.6, 0.35, 0.375, 0.75, 0.4, 0.7, 2.8, 3.5, 0.1, 0.825, 1.85, 4.8, 3.81, 0.06, 3.81, 0.381, 0.0045, 0.07, 0.35, 0.55, 0.71, 0.1, 0.21, 0.175, 0.085, 0.656, 0.786, 0.488, 0.044, 0.166, 1.233, 1.233, 0.07, 3.2, 1.0, 1.875, 0.07, 0.7, 19.0, 2.6, 0.55, 0.47, 0.6, 0.25, 0.1, 1.52, 0.15, 2.65, 4.0, 1.0, 0.35, 0.2, 0.2, 0.25, 17.0, 3.1, 90.0, 200.0, 330.0], 992.433), 
				Bidder("15", [6.0, 100.0, 30.0, 12.0, 6.5, 2.0, 0.2, 0.655, 1.6, 2.0, 1.5, 0.0001, 0.001, 0.001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 1.0, 1.0, 0.575, 0.3, 0.3, 0.985, 0.3, 3.5, 6.25, 3.3, 0.15, 0.5, 3.5, 4.0, 6.0, 0.1, 6.0, 0.6, 0.025, 0.1, 0.25, 0.51, 0.69, 0.26, 0.22, 0.16, 0.12, 0.67, 0.8, 0.5, 0.05, 0.17, 1.25, 1.25, 0.08, 3.2, 0.8, 2.9, 0.082, 0.82, 15.0, 2.1, 1.3, 0.5, 0.85, 0.45, 0.062, 1.3, 0.125, 4.5, 2.0, 1.25, 0.45, 0.2, 0.2, 0.125, 17.0, 3.5, 80.0, 160.0, 291.6], 1010.56), 
				Bidder("113", [5.0, 100.0, 5.0, 5.0, 2.5, 2.5, 0.25, 1.0, 2.5, 2.0, 5.0, 0.75, 4.5, 4.0, 2.3, 2.75, 7.5, 1.0, 1.5, 1.2, 1.2, 2.0, 0.4, 0.4, 1.0, 0.6, 7.5, 3.0, 3.5, 0.05, 0.3, 2.5, 3.0, 6.5, 0.04, 6.5, 0.65, 0.1, 0.1, 0.5, 0.85, 0.85, 0.5, 0.5, 0.5, 0.1, 0.8, 0.8, 0.7, 0.1, 0.15, 1.5, 1.5, 0.1, 2.5, 1.0, 3.5, 0.2, 0.75, 15.0, 3.5, 2.0, 0.9, 0.9, 0.5, 0.175, 1.75, 0.15, 6.0, 2.2, 1.0, 0.4, 0.2, 0.25, 0.3, 27.5, 4.0, 60.0, 150.0, 250.0], 1020.65), 
				Bidder("27", [10.0, 90.0, 20.0, 10.0, 1.0, 1.0, 0.2, 1.0, 1.5, 1.5, 5.0, 0.8, 2.0, 3.0, 2.5, 2.8, 6.0, 1.0, 2.0, 1.0, 1.2, 0.7, 0.25, 0.3, 1.0, 0.15, 5.0, 3.0, 3.5, 0.1, 0.5, 2.5, 4.0, 5.0, 0.1, 5.0, 1.2, 0.05, 0.2, 0.4, 0.82, 1.0, 0.2, 0.25, 0.2, 0.08, 0.6, 1.0, 0.6, 0.05, 0.2, 1.3, 1.3, 0.08, 4.0, 1.5, 3.0, 0.2, 1.1, 15.0, 3.0, 1.0, 1.0, 0.5, 0.3, 0.1, 1.3, 0.1, 6.5, 3.0, 1.5, 0.5, 0.3, 0.3, 0.1, 18.0, 2.5, 60.0, 165.0, 350.0], 1052.92), 
				Bidder("30", [3.5, 63.25, 15.0, 12.25, 1.2, 1.5, 0.24, 0.8, 2.5, 0.6, 1.5, 0.55, 1.0, 1.95, 0.5515, 1.65, 2.3, 2.3, 2.3, 1.2, 1.2, 0.5, 0.35, 0.25, 1.0, 0.45, 3.1385, 1.87, 3.1, 0.05, 3.0, 2.25, 3.45, 5.4, 0.076, 4.8, 0.48, 0.025, 0.03, 0.65, 0.85, 1.0, 0.5, 0.25, 0.23, 0.1, 0.6562, 0.7875, 0.488, 0.044, 0.166, 1.233, 1.233, 0.0657, 3.1, 0.65, 2.4, 0.135, 0.88, 12.0, 2.625, 2.0, 0.41, 0.65, 0.5, 0.15, 0.66, 0.13, 7.7, 2.0, 1.5, 0.45, 0.16, 0.16, 0.644, 20.91, 3.325, 201.3, 186.0, 358.5], 1169.53)
				],
				[1.0, 1.0, 1.0, 1.0, 4.305, 5.64, 0.66, 6.3, 4.876, 1.313, 1.0, 1.7, 0.0, 1.936, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 3.5, 0.0, 3.0, 2.0, 2.56, 2.0, 2.026, 7.758, 2.0, 0.0, 2.0, 2.359, 0.8361, 1.919, 7.0, 2.794, 1.6, 7.57, 2.7, 0.5, 5.52, 5.14, 1.0, 2.0, 3.0, 7.3, 0.76, 1.14, 2.0, 2.0, 4.0, 2.0, 2.0, 4.3, 0.493, 1.0, 1.191, 0.0, 1.1, 1.45, 1.144, 0.954, 0.0, 1.33, 2.0, 2.0, 1.37, 6.0, 1.342, 2.0, 1.0968, 1.4, 2.78, 2.83, 0.0, 4.2301, 2.464, 1.0, 1.0, 1.0], 
				[0.966373, 1.03436, 1.03854, 1.03706, 0.899599, 10.2544, 2.6365, 0.791545, 7.95488, 1.5645, 1.05316, 1.98792, 0.97056, 1.51926, 7.57679, 7.51483, 1.00723, 0.132271, 1.00289, 2.05508, 1.17663, 0.896151, 2.10946, 2.15079, 1.92008, 0.329052, 2.08897, 5.54745, 0.325804, 2.04633, 2.04073, 2.0651, 1.27009, 1.57682, 0.723703, 1.96701, 1.73725, 2.73545, 1.96947, 1.29554, 5.99257, 4.8658, 1.05778, 2.09554, 2.09466, 6.9172, 1.03985, 0.372803, 2.29, 2.16601, 8.72816, 2.18974, 2.17816, 5.25118, 1.12995, 0.754065, 1.13315, 2.08308, 2.28083, 1.54813, 0.12893, 0.126055, 0.133338, 0.579984, 1.00542, 1.01635, 1.31789, 6.07841, 1.25241, 1.87727, 0.588432, 2.41505, 2.52315, 3.44419, 1.02171, 5.31657, 2.50114, 1.02333, 1.03668, 0.99921], 
				[1.0, 1.0, 1.0, 1.0, 5.5, 9.95, 2.5, 4.6, 7.7, 1.5, 1.0, 2.0, 1.0, 1.5, 7.5, 7.5, 1.0, 1.0, 1.0, 2.0, 2.0, 1.0, 2.0, 2.0, 2.3, 2.0, 2.0, 5.5, 2.0, 2.0, 2.0, 2.0, 1.15, 1.45, 8.5, 1.83, 1.6, 3.0, 2.0, 1.3, 6.0, 5.0, 1.0, 2.0, 2.0, 7.0, 1.0, 0.4, 2.0, 2.0, 8.0, 2.0, 2.0, 5.5, 1.2, 0.9, 1.1, 2.0, 2.2, 1.2, 0.8, 1.0, 1.0, 1.1, 1.0, 1.0, 1.4, 6.0, 1.35, 2.0, 0.575, 2.4, 2.4, 3.3, 1.0, 5.25, 2.43, 1.0, 1.0, 1.0], [5.0, 67.481, 40.0, 37.45, 1.4, 2.2, 0.2, 0.55, 1.2, 1.25, 5.411, 0.25, 2.0, 3.5, 0.5, 0.75, 0.65, 0.65, 0.65, 0.9, 0.9, 0.45, 0.25, 0.3, 0.75, 0.35, 15.0, 2.0, 5.0, 0.05, 0.75, 3.0, 2.7, 3.5, 0.044, 3.2, 0.66, 0.004, 0.09, 0.35, 0.7, 0.85, 0.025, 0.22, 0.15, 0.1, 0.4, 0.5, 0.525, 0.06, 0.1, 1.0, 1.0, 0.05, 3.5, 2.0, 2.5, 0.07, 0.5, 12.0, 1.8, 1.5, 0.36, 0.6, 0.35, 0.09, 1.0, 0.25, 7.0, 7.0, 0.6, 0.6, 0.125, 0.125, 0.075, 16.4, 3.3, 25.0, 95.15, 220.335], 
				# [6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5, 6.10217e-5],
				[0.266744, 0.0273225, 0.0421256, 0.0369319, 14.6476, 22.9854, 1.60495, 11.302, 15.6259, 0.615968, 0.0455594, 0.886161, 0.184613, 0.474107, 12.6213, 11.5017, 0.21842, 0.399528, 0.195885, 0.947867, 1.32736, 0.278993, 0.974505, 0.972442, 1.60564, 1.62458, 0.819462, 6.38442, 1.75707, 0.96529, 0.883053, 0.878317, 0.291857, 0.580641, 27.6958, 1.03356, 0.698904, 3.12308, 1.43728, 0.539892, 11.1259, 7.61012, 0.233439, 1.05164, 1.05777, 13.2366, 0.225638, 0.0384456, 0.574783, 0.655254, 11.0342, 0.462189, 0.446732, 7.68193, 0.450646, 0.312049, 0.307968, 1.05974, 1.27164, 0.367653, 0.34818, 0.433924, 0.482024, 0.375989, 0.201127, 0.325242, 0.667156, 9.42226, 0.654951, 0.85768, 0.104165, 1.50268, 1.75021, 3.37543, 0.227249, 5.03697, 1.43975, 0.0304967, 0.0331136, 0.161771], 
				0.6826130247891279, 
				[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80], 
				[-1.5131, 3.89486, -1.24173, -9.59145, 0.747652, 1.44794, -2.50622, -0.157109, 0.343891, -2.98074, -32.7223, -3.67103, -7.41064, -1.99665, -0.532382, -0.38287, -7.30409, -3.69846, -8.21912, -2.76577, -1.71772, -5.77714, -3.31557, -3.27313, -1.73864, -1.78881, 10.7598, 0.503338, 3.02247, -3.54961, -3.18483, -0.956048, -4.14559, -0.838532, -0.489197, 0.123921, -3.31727, -1.66486, -2.32753, -3.8333, -0.236908, -0.291461, -7.41735, -3.08405, -3.13491, -0.818766, -7.29965, -17.5757, -5.52017, -5.24277, -1.15959, -6.51784, -6.77796, -1.19387, -1.12625, -3.01074, -3.70539, -3.20879, -2.50567, 6.32945, -2.1918, -2.50377, -3.24424, -4.48277, -8.28796, -5.25166, -2.64572, -0.856315, 3.41898, 2.94877, -8.9902, -2.17478, -2.25734, -1.57351, -7.57005, 14.5892, 0.367736, -31.9679, 42.6842, 209.596], 
				[0.00804566, 0.0785483, 0.0509461, 0.0581106, 0.00080585, 0.000929025, 0.003343, 0.000873496, 0.00105756, 0.00522625, 0.0471063, 0.00484367, 0.011625, 0.00679004, 0.00127531, 0.00139944, 0.00982573, 0.00537167, 0.0109561, 0.00452834, 0.00323368, 0.00769242, 0.00440456, 0.00441391, 0.00307423, 0.00264208, 0.00523791, 0.00184883, 0.00244285, 0.00444661, 0.00486072, 0.00488693, 0.00845639, 0.00535941, 0.000658661, 0.0037999, 0.00491315, 0.00206156, 0.00298639, 0.00516765, 0.00115737, 0.00141005, 0.00919358, 0.00408151, 0.00405786, 0.00113496, 0.00951142, 0.022329, 0.00746763, 0.00655054, 0.00155598, 0.00928683, 0.00960816, 0.00153656, 0.00571483, 0.00618979, 0.00766555, 0.00405031, 0.00371292, 0.00700486, 0.0049311, 0.00494588, 0.00445234, 0.00627877, 0.0106705, 0.00659858, 0.00450358, 0.00136664, 0.00442366, 0.00500451, 0.0118468, 0.00342769, 0.00294291, 0.00209817, 0.00944397, 0.0022369, 0.00362224, 0.0703728, 0.0648113, 0.0132665], 
				[-0.0705317, 1.79982, 1.26654, 1.40457, -0.30851, 0.0196543, 0.108125, -0.330944, 0.0236111, 0.140789, 1.49203, 0.0197956, -0.0792282, 0.0875003, 0.0148866, 0.0109489, 0.100937, -2.13481, 0.090367, 0.0893706, -0.597983, -0.31913, 0.142725, 0.185531, -0.215395, -1.01031, 0.144729, 0.0201933, -0.93597, 0.0786908, 0.0796696, 0.107854, 0.469833, 0.2554, -0.276229, 0.158789, 0.230287, -0.0704781, -0.000628079, 0.0274074, 0.00732072, -0.00790154, 0.310958, 0.119024, 0.117503, 0.0015787, 0.242251, -0.553278, 0.556083, 0.298562, 0.0767308, 0.474636, 0.465135, -0.0217838, -0.116005, -0.424942, 0.160538, 0.106354, 0.0891904, 0.995242, -1.89333, -1.97991, -1.76723, -1.33972, 0.100579, 0.0958125, -0.0919954, 0.0177554, -0.118465, -0.108551, 0.21072, 0.033676, 0.0906773, 0.0571992, 0.160741, 0.0286561, 0.0744131, 1.25088, 1.55519, 0.0866853], 
				32.887551, 
				1.513837763112417, 
				1.5896136462082893, 
				1.1563027292545096, 
				0.18945827107692492, 
				0.25238806156330307, 
				1000)

function getBoundarySMax(alpha_max, auc, q_b, q_e, sigma_sq)
    function ce_fun(svec)
        s_max_var = svec[1]
        return getCE_riskaverse_scale(alpha_max, s_max_var, auc, q_b, q_e,sigma_sq)
    end
    risky_cost_term = alpha_max*(auc.c' * (q_b - (auc.gamma * sigma_sq * 0.5)))

    s_max_sol = nlsolve(ce_fun, [risky_cost_term]; inplace = false)
    return s_max_sol.zero[1]
end

function getCE_riskaverse_scale(alpha, s, auc, q_b, q_e , sigma_sq)

    b_opt = getBidVector( s, alpha, auc, q_b, q_e, sigma_sq)
    b_min_c = b_opt - (alpha*auc.c);

    profit_ce = (q_b' * b_min_c) - ((auc.gamma) * 0.5 * ((sigma_sq .* b_min_c)' * b_min_c));

    return profit_ce
end

function getBidVector( s, alpha, auc, q_b, q_e, sigma_sq)
    # Compute the optimal unit bid vector for each score by solving the quadratic optimization problem

    m = Model(solver=IpoptSolver(print_level=0))

    	# m = Model(solver=GurobiSolver(OutputFlag=0))
    @variables m begin
    b[i=1:auc.T] >= 0
   	end

    @NLobjective( m,
    			Min,
   			 	(sum( ( (auc.gamma^2) * 0.5 * sigma_sq[i] * (b[i]^2) - auc.gamma * (q_b[i] + auc.gamma * sigma_sq[i] * alpha * auc.c[i]) * b[i] ) for i=1:auc.T))
   				)
				 
   	@constraint(m, sum( b[i] * q_e[i] for i=1:auc.T ) == s)
			  
	status = JuMP.solve(m)
	
	b_min = getvalue(b)
	
    return b_min

end

function getDbDsVec(b, alpha, auc, q_b, q_e, sigma_sq)
# %Computes anayltical derivative of the optimal bid function of item t
# % w.r.t. score s given alpha

    pos_b = [ifelse(b[t] > 0, 1, 0) for t=1:auc.T]
	
	if (sum(sigma_sq) > 0)
    	denom = 1.0/(sum((q_e[t]^2 * pos_b[t]) / sigma_sq[t] for t in 1:auc.T))
		db_ds = [((q_e[t] * pos_b[t] / sigma_sq[t]) * denom) for t in 1:auc.T]
	else
		db_ds = [ifelse(b[t] > 0, (q_b[t]/q_e[t]), 0) for t=1:auc.T]
	end

    return(db_ds)

end

function ode_multi_item_w_risk_aversion( auc, alpha, s, qb, qo, sigma_sq)

    alpha_dist = LogNormal(auc.alpha_lgmean, auc.sigma_lgalpha)

    f(a) = pdf(alpha_dist, a);
    F(a) = cdf(alpha_dist, a);

    b_opt = getBidVector( s, alpha, auc, qb, qo, sigma_sq)
    db_ds = getDbDsVec(b_opt, alpha, auc, qb, qo, sigma_sq)

    println("u: ", s, "; t: ", alpha)

    b_min_c = b_opt - (alpha*auc.c);

    profit_ce = (qb' * b_min_c) - ((auc.gamma) * 0.5 * ((sigma_sq .* b_min_c)' * b_min_c));

    profit_term = (exp(auc.gamma * (profit_ce)) - 1)

    numerator = (f(alpha) ./ (1 - F(alpha) - 0.000000001)) * (auc.num_bidders - 1);

    denominator = db_ds' * ((auc.gamma .* qb) - (auc.gamma^2 .* sigma_sq .* b_min_c) );
    numerator = numerator * profit_term;
	println("numerator: ", numerator)
	println("denominator: ", denominator)

    ds_dc = numerator / denominator; # Note: this should be a scalar
	println("du/dt: ", ds_dc)
	println("\n")

    return ds_dc

end

function getExPostCost(auc, alpha, s, qb, qo, sigma_sq)
    b_opt = getBidVector( s, alpha, auc, qb, qo, sigma_sq)

    total_cost = b_opt' * auc.q_a
    return(total_cost)
end

function getScalingAuctionEquilibrium(auc_seqid, auc, qb_input, qe_input, sigma_sq_input, alpha_grid)

    alpha_min = alpha_grid[1]
    alpha_max = alpha_grid[length(alpha_grid)]
    alpha_span = (alpha_max,alpha_min);

    s_max = getBoundarySMax(alpha_max, auc, qb_input, qe_input, sigma_sq_input)

	ode_multi_item_w_risk_aversion_subs(ds, s, p, alpha) = (ds .= ode_multi_item_w_risk_aversion(auc, alpha, s, qb_input, qe_input, sigma_sq_input))

	println("init u: ", s_max)
	println("t span: ", alpha_span)
    eq_prob = ODEProblem(ode_multi_item_w_risk_aversion_subs, s_max, alpha_span)

    # eq_sol = DifferentialEquations.solve(eq_prob, isoutofdomain=(u,p,t) -> any(x -> ((x < 0) | (x > s_max)), u), alg_hints=[:stiff])
	eq_sol = DifferentialEquations.solve(eq_prob, Rodas5(autodiff=false), isoutofdomain=(u,p,t) -> any(x -> (x < 0)))	
	
	eq_sol_save = [eq_sol(a) for a in alpha_grid]
    eq_cost_save = [getExPostCost(auc, a, eq_sol(a) , qb_input, qe_input, sigma_sq_input) for a in alpha_grid]

    score_at_winning_alpha = eq_sol(auc.winning_alpha)
    cost_at_winning_alpha = getExPostCost(auc, auc.winning_alpha, score_at_winning_alpha, qb_input, qe_input, sigma_sq_input)

    return(AugmentedAuctionEquilibrium(auc_seqid, auc.contract_no, alpha_grid, eq_sol_save, eq_cost_save, score_at_winning_alpha, cost_at_winning_alpha))
end

function getMinimalAuctionCFcomparison(auc, gridlength=50)
    alpha_max = auc.highest_alpha*1.2;
    alpha_min = auc.lowest_alpha*0.8;

    alpha_grad = (alpha_max-alpha_min)/gridlength;

    alpha_grid = alpha_min:alpha_grad:(alpha_max); #Using a fairly coarse grid for this demo though :(

    data_eq = getScalingAuctionEquilibrium(1, auc, auc.q_a_model, auc.q_o, auc.sigma_t_sq, alpha_grid)

    results = DataFrame(data_alpha = data_eq.alpha_types,
                        data_score = data_eq.score_maps,
                        data_cost = data_eq.expost_cost
						)

    return(results, data_eq)
end

getMinimalAuctionCFcomparison(auc)

```

---

<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:** [August 18, 2018, 9:49pm UTC](https://discourse.julialang.org/t/domain-troubles-with-differential-equations-solver-with-non-auto-differentiable-function-inside-the-ode-def/13650/17 "2018-08-18T21:49:02Z")

</div>

oh, the issue is that it’s a scalar equation so it can’t do inplace. The fix for the autodiff stuff is on v0.7/v1.0 only and would be tough to backport. I can just make a branch that should work though. Do you need to be on v0.7?

---

<div class="post-metadata">

**Author:** ![shoshievass](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shoshievass/32/3260_2.png) [@shoshievass](https://discourse.julialang.org/u/shoshievass)\
**Post date:** [August 18, 2018, 9:50pm UTC](https://discourse.julialang.org/t/domain-troubles-with-differential-equations-solver-with-non-auto-differentiable-function-inside-the-ode-def/13650/18 "2018-08-18T21:50:17Z")

</div>

Ideally I’d be able to stick with 0.6.4 since JuMP doesn’t work with 0.7 yet. Would that be doable?

---

<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:** [August 18, 2018, 9:59pm UTC](https://discourse.julialang.org/t/domain-troubles-with-differential-equations-solver-with-non-auto-differentiable-function-inside-the-ode-def/13650/19 "2018-08-18T21:59:37Z")

</div>

This branch should do what you want:

[https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/tree/06finite](https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/tree/06finite)

---

<div class="post-metadata">

**Author:** ![shoshievass](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shoshievass/32/3260_2.png) [@shoshievass](https://discourse.julialang.org/u/shoshievass)\
**Post date:** [August 18, 2018, 10:01pm UTC](https://discourse.julialang.org/t/domain-troubles-with-differential-equations-solver-with-non-auto-differentiable-function-inside-the-ode-def/13650/20 "2018-08-18T22:01:39Z")

</div>

Awesome, thank you very much! Do you by chance know a reference for how to integrate branches of packages like this vs. using Pk.add(“”)?

[Next page](https://discourse.julialang.org/t/domain-troubles-with-differential-equations-solver-with-non-auto-differentiable-function-inside-the-ode-def/13650.md?page=2)