current implementation, but does not seem to reduce allocations

src/integrating.jl

@@ -151,9 +151,10 @@ function Base.show(io::IO, integrand_values::IntegrandValues)
         "\nintegrand:\n", integrand_values.integrand)
 end
 
-mutable struct SavingIntegrandAffect{IntegrandFunc, integrandType}
+mutable struct SavingIntegrandAffect{IntegrandFunc, integrandType, integrandCacheType}
     integrand_func::IntegrandFunc
     integrand_values::IntegrandValues{integrandType}
+    integrand_cache::integrandCacheType
 end
 
 function (affect!::SavingIntegrandAffect)(integrator)
@@ -166,8 +167,8 @@ function (affect!::SavingIntegrandAffect)(integrator)
         if DiffEqBase.isinplace(integrator.sol.prob)
             curu = first(get_tmp_cache(integrator))
             integrator(curu, t_temp)
-            integral .+= gauss_weights[n][i] * 
-                         affect!.integrand_func(curu, t_temp, integrator)
+            affect!.integrand_func(affect!.integrand_cache, curu, t_temp, integrator)
+            integral .+= gauss_weights[n][i] * affect!.integrand_cache
         else
             integral .+= gauss_weights[n][i] * 
                          affect!.integrand_func(integrator(t_temp), t_temp, integrator)
@@ -205,8 +206,8 @@ The outputted values are saved into `integrand_values`. Time points are found vi
     This method is currently limited to ODE solvers of order 10 or lower. Open an issue if other
     solvers are required.
 """
-function IntegratingCallback(integrand_func, integrand_values::IntegrandValues)
-    affect! = SavingIntegrandAffect(integrand_func, integrand_values)
+function IntegratingCallback(integrand_func, integrand_values::IntegrandValues, cache)
+    affect! = SavingIntegrandAffect(integrand_func, integrand_values, cache)
     condition = (u, t, integrator) -> true
     DiscreteCallback(condition, affect!, save_positions=(false,false))
 end

test/interpolating_tests.jl

@@ -55,10 +55,14 @@ function callback_saving(u, t, integrator, sol)
     temp = sol(t)
     return vjp((x) -> lotka_volterra(temp, x, t), integrator.p, u)[1]
 end
+function callback_saving_inplace(du, u, t, integrator, sol)
+    temp = sol(t)
+    du .= vjp((x) -> lotka_volterra(temp, x, t), integrator.p, u)[1]
+end
 cb = IntegratingCallback((u, t, integrator) -> callback_saving(u, t, integrator, sol), 
-    integrand_values)
-cb_inplace = IntegratingCallback((u, t, integrator) -> callback_saving(u, t, integrator, sol), 
-            integrand_values_inplace)
+    integrand_values, 0.0)
+cb_inplace = IntegratingCallback((du, u, t, integrator) -> callback_saving_inplace(du, u, t, integrator, sol), 
+            integrand_values_inplace, zeros(length(p)))
 prob_adjoint = ODEProblem((u, p, t) -> adjoint(u, p, t, sol), 
     [0.0, 0.0], 
     (tspan[end], tspan[1]), 
@@ -100,6 +104,11 @@ function adjoint_linear(u, p, t, sol)
     return -[0 b; -a 0] * u - 2.0 * (sol(t) .- 1.0)
 end
 
+function adjoint_linear_inplace(du, u, p, t, sol)
+    a, b = p
+    du .= -[0 b; -a 0] * u - 2.0 * (sol(t) .- 1.0)
+end
+
 u0 = [1.0, 1.0] #initial condition
 tspan = (0.0, 10.0) #simulation time
 p = [1.0, 2.0] # parameters
@@ -166,18 +175,33 @@ function analytical_derivative(p, t)
 end
 
 integrand_values = IntegrandValues(Vector{Float64})
+integrand_values_inplace = IntegrandValues(Vector{Float64})
 function callback_saving_linear(u, t, integrator, sol)
     return [-sol(t)[2] 0; 0 sol(t)[1]]' * u
 end
+function callback_saving_linear_inplace(du, u, t, integrator, sol)
+    du .= [-sol(t)[2] 0; 0 sol(t)[1]]' * u
+end
 cb = IntegratingCallback((u, t, integrator) -> callback_saving_linear(u, t, integrator, sol),
-    integrand_values)
+    integrand_values, 0.0)
+cb_inplace = IntegratingCallback((du, u, t, integrator) -> callback_saving_linear_inplace(du, u, t, integrator, sol),
+    integrand_values_inplace, zeros(length(p)))
 prob_adjoint = ODEProblem((u, p, t) -> adjoint_linear(u, p, t, sol),
     [0.0, 0.0],
     (tspan[end], tspan[1]),
     p,
     callback = cb)
+prob_adjoint_inplace = ODEProblem((du, u, p, t) -> adjoint_linear_inplace(du, u, p, t, sol),
+    [0.0, 0.0],
+    (tspan[end], tspan[1]),
+    p,
+    callback = cb_inplace)
 sol_adjoint = solve(prob_adjoint, Tsit5(), abstol = 1e-14, reltol = 1e-14)
+sol_adjoint_inplace = solve(prob_adjoint_inplace, Tsit5(), abstol = 1e-14, reltol = 1e-14)
 
 dGdp_new = compute_dGdp(integrand_values)
+dGdp_new_inplace = compute_dGdp(integrand_values_inplace)
 dGdp_analytical = analytical_derivative(p,tspan[end])
+
 @test isapprox(dGdp_analytical, dGdp_new, atol = 1e-11, rtol = 1e-11)
+@test isapprox(dGdp_analytical, dGdp_new_inplace, atol = 1e-11, rtol = 1e-11)