# Wrong interpolated value in \`SavingCallback\` output

**URL:** <https://discourse.julialang.org/t/wrong-interpolated-value-in-savingcallback-output/91405>\
**Category:** Modelling & Simulations\
**Created:** [December 8, 2022, 3:59am UTC](https://discourse.julialang.org/t/wrong-interpolated-value-in-savingcallback-output/91405 "2022-12-08T03:59:31Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![cuihantao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cuihantao/32/8783_2.png) [@cuihantao](https://discourse.julialang.org/u/cuihantao)\
**Post date:** [December 8, 2022, 3:59am UTC](https://discourse.julialang.org/t/wrong-interpolated-value-in-savingcallback-output/91405/1 "2022-12-08T03:59:31Z")

</div>

This is a follow-up question on [DiffEq skips over event in the first `solve` but works after `reinit!` - #16 by ChrisRackauckas](https://discourse.julialang.org/t/diffeq-skips-over-event-in-the-first-solve-but-works-after-reinit/90935/16). The code simulates a circuit in which an AC power source is supplying an RL load through a diode. I used a `SavingCallback` to save `u[1]` at fixed intervals, but the result contains data point(s) that are extremely large.

The code below is a bit long, but the suspect issue is related to the `SavingCallback`.

```julia
using DifferentialEquations
using Plots
using BenchmarkTools

function I_load!(du, u, p, t)
    du[1] = (sqrt(2)*230*sin(2*pi*50.0*t) - u[1] * p[1]) *2
end

dt = [];

# test if diode needs to be turned off
function condition_off(u, t, integrator)
    # println("test affect off at $(integrator.t), u=$(integrator.u)")
    # println("In condition_off")
    # @show integrator.t integrator.dt

    # push!(dt, integrator.dt)
    if integrator.p[2] == 1
        u[1]
    else
        return 1
    end
end

function affect_off!(integrator)

    integrator.p[1] = 1e10 # total resistance (load + diode)
    integrator.p[2] = 0
    # println("In affect_off!")

end

# test if diode needs to be turned on
function condition_on(u, t, integrator)
    # println("In condition_on")

    if integrator.p[2] == 0
        # voltage across the diode minus the minimum turn-on voltage (1)
        ret = (sqrt(2)*230*sin(2*pi*50.0*t) - (0.5 * get_du(integrator)[1] + u[1] * 100)) - 1
    else
        ret = -1
    end

end

function affect_on!(integrator)
    # println("applied affect on at t=$(integrator.t)")
    integrator.p[1] = 100 # total resistance (load + diode)
    integrator.p[2] = 1
    
end

cb_off = ContinuousCallback(condition_off, nothing, affect_off!,
             repeat_nudge=1//10) # only when condition goes negative
cb_on = ContinuousCallback(condition_on, affect_on!, nothing,
             repeat_nudge=1//10) # only when condition goes positive

saved_values = SavedValues(Float64, Tuple{Float64})

function save_u(u, t, integrator)
    @show integrator.t u[1] # <- outputs come from here
    (u[1], )
end

cb_save = SavingCallback(save_u, saved_values,
                         saveat=0:1/2000:0.06,
)

cb_set = CallbackSet(cb_off, cb_on, cb_save);

u0 = [0.0]
p = [100., 1.] # resistance, diode flag
tspan = (0.0, 0.06)

prob = ODEProblem(I_load!, u0, tspan, p)

solver = Rodas5P()

integrator = init(prob, solver,
    callback=cb_set,
    abstol=1e-5,
    reltol = 1e-6,
    dtmax=1/60,
    saveat=0:1/2000:0.06;
);

sol = solve!(integrator)
plot(sol, label="I Rodas5P")

```

Output:

```julia
integrator.t = 9.999999999999999e-5
u[1] = 0.0
integrator.t = 0.0010355149302289026
u[1] = 0.02466443779613447
integrator.t = 0.0010355149302289026
u[1] = 0.09489608425988974
integrator.t = 0.0016705724725223623
u[1] = 0.20457147552261043
integrator.t = 0.002348954342345277
u[1] = 0.3470141746597498
integrator.t = 0.0030822014110460615
u[1] = 0.5151192104287117
integrator.t = 0.0030822014110460615

...truncated...

integrator.t = 0.052876422347891805
u[1] = 0.49957763366843355
integrator.t = 0.053380394153775476
u[1] = -1.0304261248426776e40 # <- this does not seem correct
integrator.t = 0.05444431272273977
u[1] = -3.528390922450684e-8
integrator.t = 0.05444431272273977
u[1] = -2.0560975872248562e-8
integrator.t = 0.06
u[1] = -3.213201469371499e-8
integrator.t = 0.06
u[1] = -3.254155700933845e-8
integrator.t = 0.06
u[1] = -3.2109302361073413e-8
integrator.t = 0.06
u[1] = -3.0880836404674164e-8
integrator.t = 0.06
u[1] = -2.890823031148805e-8
integrator.t = 0.06
u[1] = -2.625004076845552e-8
integrator.t = 0.06
u[1] = -2.2971309978110092e-8
integrator.t = 0.06
u[1] = -1.9143565658578353e-8
integrator.t = 0.06
u[1] = -1.4844821043579996e-8
integrator.t = 0.06
u[1] = -1.0159574882427828e-8
integrator.t = 0.06
u[1] = -5.1788114400275675e-9
integrator.t = 0.06
u[1] = -4.968781781758675e-16

```

The issue is that it silently fails. If the `abstol` and `reltol` are set to other values, the extremely large values will appear at other times (but all close to discontinuities).

---

<div class="post-metadata">

**Author:** ![cuihantao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cuihantao/32/8783_2.png) [@cuihantao](https://discourse.julialang.org/u/cuihantao)\
**Post date:** [December 8, 2022, 4:22am UTC](https://discourse.julialang.org/t/wrong-interpolated-value-in-savingcallback-output/91405/2 "2022-12-08T04:22:25Z")

</div>

While debugging this issue, I found another issue. `SavingCallback` saves timestamps at the given `saveat`, and the saved values match that in `sol.u`. The issue is that the reported `integrator.t` from within the saving callback differs from the `saveat` times.

```julia
julia> saved_values.t
121-element Vector{Float64}:
 0.0
 0.0005
 0.001
 0.0015
 0.002
 0.0025
 0.003

... truncated ...

```

```julia
julia> saved_values.saveval
121-element Vector{Tuple{Float64}}:
 (0.0,)
 (0.02466443779613447,)
 (0.09489608425988974,)
 (0.20457147552261043,)
 (0.3470141746597498,)
 (0.5151192104287117,)
 (0.7014919141095118,)

... truncated ...

```

Output from the saving callback:

```julia
integrator.t = 9.999999999999999e-5
u[1] = 0.0
integrator.t = 0.0010355149302289026
u[1] = 0.02466443779613447
integrator.t = 0.0010355149302289026
u[1] = 0.09489608425988974
integrator.t = 0.0016705724725223623
u[1] = 0.20457147552261043
integrator.t = 0.002348954342345277
u[1] = 0.3470141746597498

```

What did I miss here? Is it because `integrator.t` is not modified while interpolating? Thanks.

(edit: compared the saved values with `sol` and found it to be correct, except for the extremely large value).

---

<div class="post-metadata">

**Author:** ![cuihantao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cuihantao/32/8783_2.png) [@cuihantao](https://discourse.julialang.org/u/cuihantao)\
**Post date:** [December 10, 2022, 6:18pm UTC](https://discourse.julialang.org/t/wrong-interpolated-value-in-savingcallback-output/91405/3 "2022-12-10T18:18:07Z")

</div>

> [@cuihantao](#):
>
> ```julia
> u[1] = -1.0304261248426776e40 # <- this does not seem correct
> 
> ```

Could this be a bug in the interpolation routine? I lowered the tolerances and lowered the `repeat_nudge`, but the issue persists.

---

<div class="post-metadata">

**Author:** ![cuihantao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cuihantao/32/8783_2.png) [@cuihantao](https://discourse.julialang.org/u/cuihantao)\
**Post date:** [December 10, 2022, 11:14pm UTC](https://discourse.julialang.org/t/wrong-interpolated-value-in-savingcallback-output/91405/4 "2022-12-10T23:14:37Z")

</div>

I tracked the issue down to the use of `first(get_tmp_cache(integrator))` in `SavingAffect`.

The first value saved in the cache seen by the very last `affect!` function is close to zero. But the retrieved first value by the saving callback from the cache is mutated.

The print-out of `get_tmp_cache(integrator)` from the `ContinuousCallBack` looks like

```julia
curt=0.013, curu=[0.3303721934837517] before giving to integrator
curt=0.013, curu=[0.21664209586008248] inplace update currentu
  In Cond off: ([-1.0658179181507342e-7], [0.0009123193020741382])
  In Cond off: ([-1.0658179181507342e-7], [0.0009123193020741382])
  In Cond off: ([-1.0658179181507342e-7], [0.0009123193020741382])
  In Cond off: ([-1.0658179181507342e-7], [0.0009123193020741382])
  In Cond off: ([-1.0658179181507342e-7], [0.0009123193020741382])
  In Cond off: ([-0.0006075553301955638], [0.0009123193020741382])
  In Cond off: ([-1.2564339208389116e-6], [0.0009123193020741382])
  In Cond off: ([-2.593250008765789e-9], [0.0009123193020741382])
  In Cond off: ([-5.34254518576491e-12], [0.0009123193020741382])
  In Cond off: ([-1.865278874804446e-15], [0.0009123193020741382])
  In Cond off: ([8.97371217660693e-15], [0.0009123193020741382])
  In Cond off: ([3.0533229421397907e-15], [0.0009123193020741382])
  In Cond off: ([1.0908325831142319e-15], [0.0009123193020741382])
  In Cond off: ([-8.701559074838522e-16], [0.0009123193020741382])
  In Cond off: ([9.721148422111273e-17], [0.0009123193020741382])
  In Cond off: ([9.721148422111273e-17], [0.0009123193020741382])
  In Affect off: ([-1.3433720342587264e-9], [2.674489451237605e-5]) # !!<- first here does not match the one below 
curt=0.0132, curu=[-2.4091928323155922e10] before giving to integrator
curt=0.0132, curu=[2.2505962311451964e9] inplace update currentu
  In Cond off: ([7.978805348180991e-8], [-0.7438025523740635])

```

At other normal steps, the value in the cache seen by a `ContinuousCallBack` is what the next `SavingCallBack` retrieves, such as

```julia
  In Cond off: ([-0.0002850830064895942], [3.039545393851428e-8])
  In Cond off: ([-0.00028514796297349465], [3.039545393851428e-8])
  In Cond off: ([-0.0002852128703656998], [3.039545393851428e-8])
  In Cond off: ([-0.0002852128703656998], [3.039545393851428e-8]) # !!<- first here matches below
curt=0.0134, curu=[-0.0002852128703656998] before giving to integrator
curt=0.0134, curu=[-0.0002850108807624704] inplace update currentu

```

I’m stuck here because I’m not sure where the cache gets mutated. Any help is appreciated!

---

<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:** [December 11, 2022, 2:40am UTC](https://discourse.julialang.org/t/wrong-interpolated-value-in-savingcallback-output/91405/5 "2022-12-11T02:40:54Z")

</div>

This is the same issue as [SavedValues not stored at events · Issue #557 · SciML/DifferentialEquations.jl · GitHub](https://github.com/SciML/DifferentialEquations.jl/issues/557) and [FunctionCallingCallback not correctly synchronized with events · Issue #126 · SciML/DiffEqCallbacks.jl · GitHub](https://github.com/SciML/DiffEqCallbacks.jl/issues/126) . I actually plan to fix it next week because I thought of an elegant solution. The issue is that discrete callbacks are applied prior to the continuous callback, so you get values from the “pre pullback” phase in the SavingCallback/FunctionCallingCallback. What actually needs to happen is that they need to happen after the continuous callback has changed the time interval, or have a choice in where to be applied.

---

<div class="post-metadata">

**Author:** ![cuihantao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cuihantao/32/8783_2.png) [@cuihantao](https://discourse.julialang.org/u/cuihantao)\
**Post date:** [December 11, 2022, 2:58am UTC](https://discourse.julialang.org/t/wrong-interpolated-value-in-savingcallback-output/91405/6 "2022-12-11T02:58:29Z")

</div>

Chris, thanks a lot for your reply. I started debugging into OrdinaryDiffEq.jl but will now wait for a fix.

By the way, do you know if the current Julia Debugger in VSCode can be used to step into DiffEq? I got an error in `FunctionWrappers.jl` regarding `llvmcall` not being supported in an interpreter.

---

<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:** [December 11, 2022, 8:39am UTC](https://discourse.julialang.org/t/wrong-interpolated-value-in-savingcallback-output/91405/7 "2022-12-11T08:39:43Z")

</div>

> [@cuihantao](#):
>
> By the way, do you know if the current Julia Debugger in VSCode can be used to step into DiffEq? I got an error in `FunctionWrappers.jl` regarding `llvmcall` not being supported in an interpreter.

Just use `SciMLBase.FullSpecialize` in the construction to avoid the funciton wrappers.

---

<div class="post-metadata">

**Author:** ![cuihantao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cuihantao/32/8783_2.png) [@cuihantao](https://discourse.julialang.org/u/cuihantao)\
**Post date:** [December 11, 2022, 4:34pm UTC](https://discourse.julialang.org/t/wrong-interpolated-value-in-savingcallback-output/91405/8 "2022-12-11T16:34:03Z")

</div>

Thanks! Now I recall reading that in the FAQ.

---

<div class="post-metadata">

**Author:** ![cuihantao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cuihantao/32/8783_2.png) [@cuihantao](https://discourse.julialang.org/u/cuihantao)\
**Post date:** [December 11, 2022, 8:19pm UTC](https://discourse.julialang.org/t/wrong-interpolated-value-in-savingcallback-output/91405/9 "2022-12-11T20:19:29Z")

</div>

> [@ChrisRackauckas](#):
>
> The issue is that discrete callbacks are applied prior to the continuous callback, so you get values from the “pre pullback” phase in the SavingCallback/FunctionCallingCallback.

I’m not sure if that’s the cause of my case. In `handle_callbacks!`, continuous callbacks are applied prior to discrete callbacks. This may be a separate issue.

The `u` value in the cache suddenly changes to a very large value after calling `reeval_internals_due_to_modification!` of the continuous callback. More specifically, `_ode_addsteps!` updated the cached to a large value.

I stepped into `_ode_addsteps!` invoked by `reeval_internals_due_to_modification!`. In the former, `cache.tmp` keeps diverging (as the norm of k1 through k8 keeps increasing). It might be a numerical issue, but I could be wrong.

---

<div class="post-metadata">

**Author:** ![cuihantao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cuihantao/32/8783_2.png) [@cuihantao](https://discourse.julialang.org/u/cuihantao)\
**Post date:** [December 11, 2022, 8:24pm UTC](https://discourse.julialang.org/t/wrong-interpolated-value-in-savingcallback-output/91405/10 "2022-12-11T20:24:12Z")

</div>

How can I force a Jacobian update in `_ode_addsteps!`? I saw a comment in that function saying

```julia
        ### Jacobian does not need to be re-evaluated after an event
        ### Since it's unchanged
        jacobian2W!(W, mass_matrix, dtgamma, J, true)

```

The event in my code changed a parameter and thus the Jacobian. I’m poking around…

---

<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:** [December 11, 2022, 11:57pm UTC](https://discourse.julialang.org/t/wrong-interpolated-value-in-savingcallback-output/91405/11 "2022-12-11T23:57:00Z")

</div>

> [@cuihantao](#):
>
> The event in my code changed a parameter and thus the Jacobian.

No that happens after to it. What’s going on is that you start with an interval `(t, t+dt)`. Then you find an event at `t + theta * dt` for `theta in [0,1]`. You then shrink the interval to `(t, t + theta * dt)` and re-evaluate the stepper to get a valid interpolation over that new interval. In this interval you’d then want to run the callbacks. But they should know that the real `dt` is `theta * dt`.

---

<div class="post-metadata">

**Author:** ![cuihantao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cuihantao/32/8783_2.png) [@cuihantao](https://discourse.julialang.org/u/cuihantao)\
**Post date:** [December 12, 2022, 12:03am UTC](https://discourse.julialang.org/t/wrong-interpolated-value-in-savingcallback-output/91405/12 "2022-12-12T00:03:42Z")

</div>

I figured out how your described event timing and time modification work. I suspect that due to a silent inexact error, `t + theta * dt` ends up on the other side of the zero crossing. In this case, shouldn’t the Jacobian be reevaluated? My problem has large parameters `1e10` so the exact error can be possible.

---

<div class="post-metadata">

**Author:** ![cuihantao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cuihantao/32/8783_2.png) [@cuihantao](https://discourse.julialang.org/u/cuihantao)\
**Post date:** [December 12, 2022, 12:27am UTC](https://discourse.julialang.org/t/wrong-interpolated-value-in-savingcallback-output/91405/13 "2022-12-12T00:27:08Z")

</div>

If I reduce that huge parameter from `1e10` to `1e3`, the result is correct. The solution is different due to the parameter change.

 ![image](https://global.discourse-cdn.com/julialang/original/3X/6/0/60bba0e392668cfc2efe1d02f4b1132c7fce7a53.png)

If the parameter is changed to `1e4`, the issue appears:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/2/c/2c12c5fb9d99e6590d51b85e4d8cbf4beec448c2.png)

If the parameter is increased to `1e5`, the issue appears at different times as shown below. All these times are in the proximity of events.

 ![image](https://global.discourse-cdn.com/julialang/original/3X/c/e/ce60abd4335b4072d56a8e4721f8c56c48fb9a8d.png)

---

<div class="post-metadata">

**Author:** ![cuihantao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cuihantao/32/8783_2.png) [@cuihantao](https://discourse.julialang.org/u/cuihantao)\
**Post date:** [December 12, 2022, 12:31am UTC](https://discourse.julialang.org/t/wrong-interpolated-value-in-savingcallback-output/91405/14 "2022-12-12T00:31:54Z")

</div>

I’ll stop here and wait for your fix. Thank you for your patience. 😄