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

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:

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]

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 ?

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,

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

Hopefully these narrow down the problem.

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:

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:

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 :sweat_smile: but I’ve used similar things with other ODEs that seem to work fine.

Thank you very much once again!

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

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.

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 :man_shrugging:.

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

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

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

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…

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

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

Oh, you’re running into 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.

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)

Ah! Cool :smiley:

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

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

Is your uType or tType an integer?

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

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.

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

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.

Initial conditions:

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

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!

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)

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?

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

This branch should do what you want:

https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/tree/06finite

1 Like

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("")?