# Parameter estimation with algebraic output variable and discrete callbacks

**URL:** https://discourse.julialang.org/t/parameter-estimation-with-algebraic-output-variable-and-discrete-callbacks/34996
**Category:** New to Julia
**Tags:** question, diffeq
**Created:** [February 21, 2020, 8:55pm UTC](https://discourse.julialang.org/t/parameter-estimation-with-algebraic-output-variable-and-discrete-callbacks/34996 "2020-02-21T20:55:01Z")
**Posts on this page:** 5
**Page:** 1

<div class="post-metadata">

### Author: ![rotton](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rotton/32/13746_2.png) [@rotton](https://discourse.julialang.org/u/rotton)
#### Post date: [February 21, 2020, 8:55pm UTC](https://discourse.julialang.org/t/parameter-estimation-with-algebraic-output-variable-and-discrete-callbacks/34996/1 "2020-02-21T20:55:01Z")

</div>

Trying to setup a Parameter Estimation in Julia with `DifferentialEquations`, I started reading the documentation and some topics here on Discourse. But it seems I am actually facing _three_ issues at once, which is why I am opening this topic.  
I have a simple ODE with two states. Input parameters change at certain times, which I am handling through `PresetTimeCallback`, with additional suggestions from [here](https://stackoverflow.com/questions/59908626/how-to-read-corresponding-value-at-specified-time-in-discretecallback). (Using Interpolations.jl instead didn’t work, since [`previous` interpolation is yet to be implemented](https://github.com/JuliaMath/Interpolations.jl/pull/338).)  
Consider this MWE:

```julia
using DifferentialEquations
using DelimitedFiles
using Plots

function odm2prod(dx, x, adjParams, t, fixedParams)
	V_liq, X_in, Y_in = fixedParams
	k_1, f_1, q_in = adjParams

	rho_1 = k_1*x[1]
#	q_prod = 0.52*f_1*x[1]
        # Differential Equations
	dx[1] = q_in/V_liq*(X_in - x[1]) - rho_1
	dx[2] = q_in/V_liq*(Y_in - x[2])
end

x0 = [11.5, 11.55*0.2]
tspan = (0.0, 7.0)
fixedParams = [159, 249, 58]
adjParams = [0.19, 29.0, 0]
prob = ODEProblem((dx, x, adjParams, t) -> odm2prod(dx, x, adjParams,
    t, fixedParams), x0, tspan, adjParams)

inputArray = readdlm("rawDataFiles/Inputarray_Odm2prod_2h_7d.csv", ',',
    Float64)
dosetimes = inputArray[:,1]/86400

function affect!(integrator)
	ind_t = findall(t -> t==integrator.t, dosetimes)
	integrator.p[3] = inputArray[ind_t[1], 2]
end
cb = PresetTimeCallback(dosetimes, affect!)
sol = solve(prob, Tsit5(), callback=cb)

```

(Actually, `X_in` and `Y_in` should also change at every callback, but I did not implement this due to problem #3, _see below_)

> **Content of input CSV file**
>
> 3300,430,32.64,32.64  
> 3521.0232558,0,32.64,32.64  
> 3660,60,249.83539555,58.001892497  
> 3828.0000005,0,249.83539555,58.001892497  
> 10860,60,249.83539555,58.001892497  
> 11028,0,249.83539555,58.001892497  
> 18060.005,60,249.83539555,58.001892497  
> 18228.005,0,249.83539555,58.001892497  
> 25260.01,60,249.83539555,58.001892497  
> 25428.01,0,249.83539555,58.001892497  
> 32460.015,60,249.83539555,58.001892497  
> 32628.015,0,249.83539555,58.001892497  
> 39660.02,60,249.83539555,58.001892497  
> 39828.02,0,249.83539555,58.001892497  
> 46860.025,60,249.83539555,58.001892497  
> 47028.025,0,249.83539555,58.001892497  
> 54060.03,60,249.83539555,58.001892497  
> 54228.03,0,249.83539555,58.001892497  
> 61260.035,60,249.83539555,58.001892497  
> 61428.035,0,249.83539555,58.001892497  
> 68460.04,60,249.83539555,58.001892497  
> 68628.04,0,249.83539555,58.001892497  
> 75660.045,60,249.83539555,58.001892497  
> 75828.045,0,249.83539555,58.001892497  
> 82860.05,60,249.83539555,58.001892497  
> 83028.05,0,249.83539555,58.001892497  
> 89700,430,32.64,32.64  
> 89921.023256,0,32.64,32.64  
> 90060,60,249.83539555,58.001892497  
> 90228,0,249.83539555,58.001892497  
> 97260,60,249.83539555,58.001892497  
> 97428,0,249.83539555,58.001892497  
> 104460.005,60,249.83539555,58.001892497  
> 104628.005,0,249.83539555,58.001892497  
> 111660.01,60,249.83539555,58.001892497  
> 111828.01,0,249.83539555,58.001892497  
> 118860.015,60,249.83539555,58.001892497  
> 119028.015,0,249.83539555,58.001892497  
> 126060.02,60,249.83539555,58.001892497  
> 126228.02,0,249.83539555,58.001892497  
> 133260.025,60,249.83539555,58.001892497  
> 133428.025,0,249.83539555,58.001892497  
> 140460.03,60,249.83539555,58.001892497  
> 140628.03,0,249.83539555,58.001892497  
> 147660.035,60,249.83539555,58.001892497  
> 147828.035,0,249.83539555,58.001892497  
> 154860.04,60,249.83539555,58.001892497  
> 155028.04,0,249.83539555,58.001892497  
> 162060.045,60,249.83539555,58.001892497  
> 162228.045,0,249.83539555,58.001892497  
> 169260.05,60,249.83539555,58.001892497  
> 169428.05,0,249.83539555,58.001892497  
> 176100,430,32.64,32.64  
> 176321.02326,0,32.64,32.64  
> 176460,60,249.83539555,58.001892497  
> 176628,0,249.83539555,58.001892497  
> 183660,60,249.83539555,58.001892497  
> 183828,0,249.83539555,58.001892497  
> 190860.005,60,249.83539555,58.001892497  
> 191028.005,0,249.83539555,58.001892497  
> 198060.01,60,249.83539555,58.001892497  
> 198228.01,0,249.83539555,58.001892497  
> 205260.015,60,249.83539555,58.001892497  
> 205428.015,0,249.83539555,58.001892497  
> 212460.02,60,249.83539555,58.001892497  
> 212628.02,0,249.83539555,58.001892497  
> 219660.025,60,249.83539555,58.001892497  
> 219828.025,0,249.83539555,58.001892497  
> 226860.03,60,249.83539555,58.001892497  
> 227028.03,0,249.83539555,58.001892497  
> 234060.035,60,249.83539555,58.001892497  
> 234228.035,0,249.83539555,58.001892497  
> 241260.04,60,249.83539555,58.001892497  
> 241428.04,0,249.83539555,58.001892497  
> 248460.045,60,249.83539555,58.001892497  
> 248628.045,0,249.83539555,58.001892497  
> 255660.05,60,249.83539555,58.001892497  
> 255828.05,0,249.83539555,58.001892497  
> 262500,430,32.64,32.64  
> 262721.02326,0,32.64,32.64  
> 262860,60,249.83539555,58.001892497  
> 263028,0,249.83539555,58.001892497  
> 270060,60,249.83539555,58.001892497  
> 270228,0,249.83539555,58.001892497  
> 277260.005,60,249.83539555,58.001892497  
> 277428.005,0,249.83539555,58.001892497  
> 284460.01,60,249.83539555,58.001892497  
> 284628.01,0,249.83539555,58.001892497  
> 291660.015,60,249.83539555,58.001892497  
> 291828.015,0,249.83539555,58.001892497  
> 298860.02,60,249.83539555,58.001892497  
> 299028.02,0,249.83539555,58.001892497  
> 306060.025,60,249.83539555,58.001892497  
> 306228.025,0,249.83539555,58.001892497  
> 313260.03,60,249.83539555,58.001892497  
> 313428.03,0,249.83539555,58.001892497  
> 320460.035,60,249.83539555,58.001892497  
> 320628.035,0,249.83539555,58.001892497  
> 327660.04,60,249.83539555,58.001892497  
> 327828.04,0,249.83539555,58.001892497  
> 334860.045,60,249.83539555,58.001892497  
> 335028.045,0,249.83539555,58.001892497  
> 342060.05,60,249.83539555,58.001892497  
> 342228.05,0,249.83539555,58.001892497  
> 348900,430,32.64,32.64  
> 349121.02326,0,32.64,32.64  
> 349260,60,249.83539555,58.001892497  
> 349428,0,249.83539555,58.001892497  
> 356460,60,249.83539555,58.001892497  
> 356628,0,249.83539555,58.001892497  
> 363660.005,60,249.83539555,58.001892497  
> 363828.005,0,249.83539555,58.001892497  
> 370860.01,60,249.83539555,58.001892497  
> 371028.01,0,249.83539555,58.001892497  
> 378060.015,60,249.83539555,58.001892497  
> 378228.015,0,249.83539555,58.001892497  
> 385260.02,60,249.83539555,58.001892497  
> 385428.02,0,249.83539555,58.001892497  
> 392460.025,60,249.83539555,58.001892497  
> 392628.025,0,249.83539555,58.001892497  
> 399660.03,60,249.83539555,58.001892497  
> 399828.03,0,249.83539555,58.001892497  
> 406860.035,60,249.83539555,58.001892497  
> 407028.035,0,249.83539555,58.001892497  
> 414060.04,60,249.83539555,58.001892497  
> 414228.04,0,249.83539555,58.001892497  
> 421260.045,60,249.83539555,58.001892497  
> 421428.045,0,249.83539555,58.001892497  
> 428460.05,60,249.83539555,58.001892497  
> 428628.05,0,249.83539555,58.001892497  
> 435300,430,32.64,32.64  
> 435521.02326,0,32.64,32.64  
> 435660,60,249.83539555,58.001892497  
> 435828,0,249.83539555,58.001892497  
> 442860,60,249.83539555,58.001892497  
> 443028,0,249.83539555,58.001892497  
> 450060.005,60,249.83539555,58.001892497  
> 450228.005,0,249.83539555,58.001892497  
> 457260.01,60,249.83539555,58.001892497  
> 457428.01,0,249.83539555,58.001892497  
> 464460.015,60,249.83539555,58.001892497  
> 464628.015,0,249.83539555,58.001892497  
> 471660.02,60,249.83539555,58.001892497  
> 471828.02,0,249.83539555,58.001892497  
> 478860.025,60,249.83539555,58.001892497  
> 479028.025,0,249.83539555,58.001892497  
> 486060.03,60,249.83539555,58.001892497  
> 486228.03,0,249.83539555,58.001892497  
> 493260.035,60,249.83539555,58.001892497  
> 493428.035,0,249.83539555,58.001892497  
> 500460.04,60,249.83539555,58.001892497  
> 500628.04,0,249.83539555,58.001892497  
> 507660.045,60,249.83539555,58.001892497  
> 507828.045,0,249.83539555,58.001892497  
> 514860.05,60,249.83539555,58.001892497  
> 515028.05,0,249.83539555,58.001892497  
> 521700,430,32.64,32.64  
> 521921.02326,0,32.64,32.64  
> 522060,60,249.83539555,58.001892497  
> 522228,0,249.83539555,58.001892497  
> 529260,60,249.83539555,58.001892497  
> 529428,0,249.83539555,58.001892497  
> 536460.005,60,249.83539555,58.001892497  
> 536628.005,0,249.83539555,58.001892497  
> 543660.01,60,249.83539555,58.001892497  
> 543828.01,0,249.83539555,58.001892497  
> 550860.015,60,249.83539555,58.001892497  
> 551028.015,0,249.83539555,58.001892497  
> 558060.02,60,249.83539555,58.001892497  
> 558228.02,0,249.83539555,58.001892497  
> 565260.025,60,249.83539555,58.001892497  
> 565428.025,0,249.83539555,58.001892497  
> 572460.03,60,249.83539555,58.001892497  
> 572628.03,0,249.83539555,58.001892497  
> 579660.035,60,249.83539555,58.001892497  
> 579828.035,0,249.83539555,58.001892497  
> 586860.04,60,249.83539555,58.001892497  
> 587028.04,0,249.83539555,58.001892497  
> 594060.045,60,249.83539555,58.001892497  
> 594228.045,0,249.83539555,58.001892497  
> 601260.05,60,249.83539555,58.001892497  
> 601428.05,0,249.83539555,58.001892497

For completeness, here are the observations, with the columns `t` (in seconds; needs to be divided by 86400, in order to match unit of `tspan`), `x[1]`, `q_prod` (needs to be multiplied by 24, in order to match unit of `q_prod` in code):

> **Content of measurements CSV file**
>
> 0,11.3526,14.02611  
> 3600,11.41204,5.990548  
> 7200,11.48661,14.13363  
> 10800,11.39558,13.81184  
> 14400,11.47054,14.30515  
> 18000,11.37988,13.85038  
> 21600,11.45519,14.30656  
> 25200,11.36489,13.83604  
> 28800,11.44058,14.28658  
> 32400,11.35062,13.81046  
> 36000,11.42665,14.25907  
> 39600,11.33701,13.77983  
> 43200,11.41337,14.22781  
> 46800,11.32404,13.74663  
> 50400,11.40072,14.19458  
> 54000,11.31169,13.71217  
> 57600,11.38867,14.16042  
> 61200,11.29993,13.67726  
> 64800,11.37718,14.12593  
> 68400,11.28871,13.64237  
> 72000,11.36625,14.09166  
> 75600,11.27804,13.6079  
> 79200,11.35585,14.05792  
> 82800,11.26789,13.57412  
> 86400,11.34595,14.02487  
> 90000,11.40542,5.988436  
> 93600,11.48,14.13261  
> 97200,11.38898,13.81054  
> 100800,11.46395,14.30391  
> 104400,11.3733,13.84899  
> 108000,11.44865,14.30534  
> 111600,11.35835,13.8347  
> 115200,11.43404,14.28531  
> 118800,11.34408,13.8091  
> 122400,11.42012,14.25784  
> 126000,11.33049,13.77849  
> 129600,11.40687,14.22666  
> 133200,11.31755,13.74535  
> 136800,11.39422,14.19333  
> 140400,11.30519,13.71084  
> 144000,11.38219,14.15921  
> 147600,11.29345,13.67594  
> 151200,11.37071,14.12471  
> 154800,11.28225,13.64105  
> 158400,11.3598,14.09051  
> 162000,11.2716,13.60663  
> 165600,11.34943,14.05678  
> 169200,11.26147,13.57286  
> 172800,11.33954,14.02368  
> 176400,11.39906,5.985696  
> 180000,11.47366,14.13173  
> 183600,11.38264,13.80937  
> 187200,11.45763,14.3028  
> 190800,11.36698,13.84773  
> 194400,11.44234,14.30423  
> 198000,11.35205,13.83347  
> 201600,11.42775,14.28423  
> 205200,11.3378,13.80791  
> 208800,11.41384,14.25672  
> 212400,11.32421,13.77728  
> 216000,11.40058,14.22547  
> 219600,11.31127,13.74408  
> 223200,11.38796,14.19223  
> 226800,11.29894,13.70962  
> 230400,11.37593,14.15805  
> 234000,11.28719,13.67469  
> 237600,11.36448,14.12362  
> 241200,11.27602,13.63984  
> 244800,11.35358,14.08939  
> 248400,11.26539,13.60541  
> 252000,11.34322,14.05565  
> 255600,11.25527,13.57164  
> 259200,11.33335,14.02261  
> 262800,11.39295,5.979839  
> 266400,11.46755,14.13088  
> 270000,11.37654,13.8083  
> 273600,11.45153,14.30173  
> 277200,11.36089,13.84653  
> 280800,11.43624,14.30311  
> 284400,11.34595,13.83227  
> 288000,11.42167,14.28316  
> 291600,11.33172,13.80673  
> 295200,11.40778,14.25573  
> 298800,11.31816,13.77616  
> 302400,11.39455,14.2245  
> 306000,11.30523,13.74299  
> 309600,11.38193,14.19126  
> 313200,11.29292,13.70855  
> 316800,11.36992,14.1571  
> 320400,11.28119,13.67363  
> 324000,11.35848,14.12262  
> 327600,11.27002,13.63875  
> 331200,11.34759,14.08839  
> 334800,11.2594,13.60431  
> 338400,11.33724,14.05465  
> 342000,11.24929,13.57054  
> 345600,11.32737,14.02153  
> 349200,11.38695,5.98036  
> 352800,11.46156,14.12997  
> 356400,11.37055,13.80708  
> 360000,11.44554,14.30051  
> 363600,11.3549,13.84516  
> 367200,11.43028,14.3019  
> 370800,11.34,13.83092  
> 374400,11.41572,14.28194  
> 378000,11.32578,13.8054  
> 381600,11.40186,14.25456  
> 385200,11.31224,13.77487  
> 388800,11.38864,14.22333  
> 392400,11.29933,13.74171  
> 396000,11.37604,14.19009  
> 399600,11.28703,13.70726  
> 403200,11.36404,14.15591  
> 406800,11.27531,13.67234  
> 410400,11.35261,14.12148  
> 414000,11.26416,13.63749  
> 417600,11.34174,14.08725  
> 421200,11.25356,13.60306  
> 424800,11.33139,14.05345  
> 428400,11.24345,13.56925  
> 432000,11.32154,14.02036  
> 435600,11.38121,5.97516  
> 439200,11.45584,14.12915  
> 442800,11.36483,13.80602  
> 446400,11.43984,14.2995  
> 450000,11.3492,13.84401  
> 453600,11.42458,14.30084  
> 457200,11.3343,13.82977  
> 460800,11.41004,14.28096  
> 464400,11.3201,13.80431  
> 468000,11.39617,14.25352  
> 471600,11.30656,13.77375  
> 475200,11.38297,14.22233  
> 478800,11.29366,13.7406  
> 482400,11.37037,14.18904  
> 486000,11.28136,13.70612  
> 489600,11.35837,14.15484  
> 493200,11.26965,13.67118  
> 496800,11.34695,14.12039  
> 500400,11.2585,13.63631  
> 504000,11.33609,14.08616  
> 507600,11.24791,13.60188  
> 511200,11.32576,14.05241  
> 514800,11.23783,13.5681  
> 518400,11.31593,14.01936  
> 522000,11.37565,5.971109  
> 525600,11.45027,14.12832  
> 529200,11.35927,13.80494  
> 532800,11.43428,14.29842  
> 536400,11.34365,13.84281  
> 540000,11.41904,14.29982  
> 543600,11.32877,13.82862  
> 547200,11.4045,14.27986  
> 550800,11.31457,13.80312  
> 554400,11.39064,14.25243  
> 558000,11.30103,13.77256  
> 561600,11.37744,14.22124  
> 565200,11.28814,13.73942  
> 568800,11.36487,14.18803  
> 572400,11.27587,13.70499  
> 576000,11.3529,14.15386  
> 579600,11.26418,13.67009  
> 583200,11.3415,14.11943  
> 586800,11.25306,13.63524  
> 590400,11.33065,14.08518  
> 594000,11.24247,13.6008  
> 597600,11.32032,14.05139  
> 601200,11.23239,13.567  
> 604800,11.3105,14.01833

Now, since my observations are for `x[1]` and `q_prod`, I need `q_prod` as output variable from the `ODEProblem`. I managed to write this function, which I am not sure will be of any help for this endeavor:

```julia
# Out-function for odm2prod to get algebraic states
function odm2prod_out(sol, adjParams, fixedParams)
	V_liq, X_in, Y_in = fixedParams
	k_1, f_1, q_in = adjParams

	x = Array(sol)
	# Algebraic equation
	q_prod = 0.52*f_1*x[1,:]
	return q_prod
end

y = odm2prod_out(sol, adjParams, fixedParams)

plot(sol.t, y/24)

```

So these are my problems:

1. How to get an algebraic output variable, which is dependent on the states and parameters, as part of `sol`?
2. While the callbacks may occur at any time, the observations are equally spaced, at fixed interval, so the number and location of points I get in `sol` are not the same as my measured times, even with the `saveat=` option enabled (I checked). How would I reconcile both timescales?
3. If I understood correctly, [_all_ parameters](https://discourse.julialang.org/t/parameter-estimation-with-differentialequations-jl-for-a-subset-of-parameters/33851) in `adjParams` would be adjusted by Parameter Estimation routine, but some **must be** adjusted by the callback. How would I separate them?

Thanks in advance for any help, it is very much 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: [February 21, 2020, 9:56pm UTC](https://discourse.julialang.org/t/parameter-estimation-with-algebraic-output-variable-and-discrete-callbacks/34996/2 "2020-02-21T21:56:41Z")

</div>

> [@rotton](#):
>
> How to get an algebraic output variable, which is dependent on the states and parameters, as part of `sol` ?

do the algebraic calculation on the output. We may make a better interface for this in the future, but it’s not a big priority.

> [@rotton](#):
>
> While the callbacks may occur at any time, the observations are equally spaced, at fixed interval, so the number and location of points I get in `sol` are not the same as my measured times, even with the `saveat=` option enabled (I checked). How would I reconcile both timescales?

Use the `save_positions` keyword argument in the callback. Specifically, `save_positions=(false,false)` so the callback doesn’t add any time points.

> [@rotton](#):
>
> If I understood correctly, [_all_ parameters](https://discourse.julialang.org/t/parameter-estimation-with-differentialequations-jl-for-a-subset-of-parameters/33851) in `adjParams` would be adjusted by Parameter Estimation routine, but some **must be** adjusted by the callback. How would I separate them?

Use a constant global ref or something of the sort. This is something I hope to have a better solution for soon.

---

<div class="post-metadata">

### Author: ![rotton](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rotton/32/13746_2.png) [@rotton](https://discourse.julialang.org/u/rotton)
#### Post date: [February 22, 2020, 9:02am UTC](https://discourse.julialang.org/t/parameter-estimation-with-algebraic-output-variable-and-discrete-callbacks/34996/3 "2020-02-22T09:02:18Z")

</div>

> [@ChrisRackauckas](#):
>
> Use the `save_positions` keyword argument in the callback. Specifically, `save_positions=(false,false)` so the callback doesn’t add any time points.

Ok, this worked. 👍

> [@ChrisRackauckas](#):
>
> do the algebraic calculation on the output.

Does this mean I need to write a custom loss function? Or where would I perform this calculation? Inside `odm2prod`? 🤔

> [@ChrisRackauckas](#):
>
> Use a constant global ref or something of the sort.

You mean define the input parameters which are changed by the callback as global variables, and leave the parameters to be optimized inside `adjParams`? I tried, but it doesn’t work. Changes are not having any affect:

```julia
using DifferentialEquations
using DelimitedFiles
using Plots

function odm2prod(dx, x, adjParams, t, fixedParams)
	V_liq = fixedParams
	k_1, f_1 = adjParams

	rho_1 = k_1*x[1]
#	q_prod = 0.52*f_1*x[1]
    # Differential Equations
	dx[1] = q_in/V_liq*(X_in - x[1]) - rho_1
	dx[2] = q_in/V_liq*(Y_in - x[2])
end

global q_in = 0
global X_in = 32.64
global Y_in = 32.64
x0 = [11.5, 11.55*0.2]
tspan = (0.0, 7.0)
const fixedParams = 159
adjParams = [0.19, 29.0, 0]
prob = ODEProblem((dx, x, adjParams, t) -> odm2prod(dx, x, adjParams,
	t, fixedParams), x0, tspan, adjParams)

inputArray = readdlm("rawDataFiles/Inputarray_Odm2prod_2h_7d.csv",
	',', Float64)
dosetimes = inputArray[:,1]/86400

function affect!(integrator)
	ind_t = findall(t -> t==integrator.t, dosetimes)
	q_in = inputArray[ind_t[1], 2]
	X_in = inputArray[ind_t[1], 3]
	Y_in = inputArray[ind_t[1], 4]
end
cb = PresetTimeCallback(dosetimes, affect!, save_positions=(false, false))
sol = solve(prob, Tsit5(), callback=cb, saveat=1/24)

plot(sol, vars=[1, 2])

```

---

<div class="post-metadata">

### Author: ![rotton](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rotton/32/13746_2.png) [@rotton](https://discourse.julialang.org/u/rotton)
#### Post date: [March 31, 2020, 8:15am UTC](https://discourse.julialang.org/t/parameter-estimation-with-algebraic-output-variable-and-discrete-callbacks/34996/4 "2020-03-31T08:15:31Z")

</div>

> [@ChrisRackauckas](#):
>
> Use a constant global ref or something of the sort. This is something I hope to have a better solution for soon.

Ok, I have figured out that I needed to put the `global` keyword also inside the `affect!` function. Then it does what it is supposed to.

```julia
using DifferentialEquations
using DelimitedFiles
using Plots

function odm2prod(dx, x, adjParams, t, fixedParams)
	V_liq = fixedParams
	k_1, f_1 = adjParams

	rho_1 = k_1*x[1]
#	q_prod = 0.52*f_1*x[1]
    # Differential Equations
	dx[1] = q_in/V_liq*(X_in - x[1]) - rho_1
	dx[2] = q_in/V_liq*(Y_in - x[2])
end

global q_in = 0
global X_in = 32.64
global Y_in = 32.64
x0 = [11.5, 11.55*0.2]
tspan = (0.0, 7.0)
fixedParams = 159
adjParams = [0.19, 29.0, 0]
prob = ODEProblem((dx, x, adjParams, t) -> odm2prod(dx, x, adjParams,
	t, fixedParams), x0, tspan, adjParams)

inputArray = readdlm("rawDataFiles/Inputarray_Odm2prod_2h_7d.csv",
	',', Float64)
dosetimes = inputArray[:,1]/86400

function affect!(integrator)
	ind_t = findall(t -> t==integrator.t, dosetimes)
	global q_in = inputArray[ind_t[1], 2]
	global X_in = inputArray[ind_t[1], 3]
	global Y_in = inputArray[ind_t[1], 4]
end
cb = PresetTimeCallback(dosetimes, affect!, save_positions=(false, false))
sol = solve(prob, Tsit5(), callback=cb, saveat=1/24)

plot(sol, vars=[1, 2])

```

But I still don’t know how to

> [@ChrisRackauckas](#):
>
> do the algebraic calculation on the output.

and how to generate the loss function with it.  
Could you give me a hint here, please?

---

<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: [March 31, 2020, 3:12pm UTC](https://discourse.julialang.org/t/parameter-estimation-with-algebraic-output-variable-and-discrete-callbacks/34996/5 "2020-03-31T15:12:25Z")

</div>

> [@rotton](#):
>
> and how to generate the loss function with it.  
> Could you give me a hint here, please?

I just mean you do `xsquared = sol[1,:].^2` to get the timeseries of the first one squared, and stuff like that.
