Explicit Jacobian on DAEs?

Specific Domains Numerics
41

What would you like me to say in the issue?
LocalSensitivityAnalysis does not work with DAe or MassMatix ODE?
for adjoint can you use the final value of one of your variables say u[18] as the cost function?

How would you approach this using Forwarddiff?
Would Sobol or Morris method work/

thanks-

A

42

Sobol and Morris will work. They give qualitatively different answer though since they are global and not local.

http://docs.juliadiffeq.org/latest/analysis/sensitivity.html#Examples-using-ForwardDiff.jl-1

43

Thanks -
Yes I know they will be different - my main goal is to see how changing parameters such as kon, koff, IC50 of drug, etc influences tumor size -
Thatā€™s why I was looking for local sensitivity, but I wonder if this the right approcah - most papers use sobol now a days, but it seems harder to interpret.

Take care.

44

Global sensitivity is probably a better measure for this kind of thing. Itā€™s a little harder to interpret but ā€œa lot more correctā€ in some sense.

45

Thanks-
Iā€™ll go that way then, thanks for the advice!

A

46

Hi -
Sorry to keep bothering you, but A quick question -
Using adjoint_sensitivities - I get a type error, and I am not understanding how it arises. Is again something I shouldnā€™t be doing with MassMatrix ODE?

This is the error:
TypeError: in typeassert, expected Float64, got ForwardDiff.Dual{Nothing,Float64,1}
TypeError: in setindex!, in typeassert, expected Float64, got ForwardDiff.Dual{Nothing,Float64,1}

Here is the code:


# Generate System
X=vcat(U,P)                         # X concatenates U and P to generate equations as function of U and P
Eqs=lambdify(A,X)                   # Generate matrix with all system equations based on X
 
# Generate Jacobian w/r to vars
Jvars=jacobian(A*U,U)               # Symbolic jacobian of system for fixed parameters (p) with respct to variables (U)
Jac_vars=lambdify(Jvars,X)          # Evaluates symbolic jacobian of system at fixed p for given U=u, for example u=IC

# Generate Jacobian w/r to params
Jparams=jacobian(A*U,P)             # Symbolic jacobian of system for fixed variables (u) with respct to parameters (P)
Jac_params=lambdify(Jparams,X)      # Evaluates symbolic jacobian of system at fixed u for given P=p


s=size(A)
varnum=21
parnum=32
algnum=6
dv=varnum - algnum
diff_vars=[if iā‰¤dv true else false end for i in 1:varnum ]
MM = zeros(varnum,varnum)
[MM[i,i] = 1 for i in 1:dv]
c=zeros(varnum)

function ODE_Sym(du,u,p,t)
    E1=convert(Array{Float64,2}, Eqs(vcat(u,p)...))
    for i in 1:s[1]
        if (iā‰¤dv)
            du[i]=(BLAS.gemv('N', E1, u)[i]) - c[i]
        else
            du[i]=(BLAS.gemv('N', E1, u)[i]) - c[i]
        end
    end
#     mul!(du,Eqs(vcat(u,p)...),u) 
end

function Jac_ODE_Sym(J,u,p,t)
    J1=convert(Array{Float64,2}, Jac_vars(vcat(u,p)...))
    for i in 1:s[1]
        for j in 1:s[2]
            J[i,j]=J1[i,j]
        end
    end
end


function paramJac_ODE_Sym(pJ,u,p,t)
    pJ1=convert(Array{Float64,2}, Jac_params(vcat(u,p)...))
    for i in 1:varnum
        for j in 1:parnum
            pJ[i,j]=pJ1[i,j]
        end
    end
end


ODE_test = ODEFunction(ODE_Sym;mass_matrix=MM,jac=Jac_ODE_Sym,paramjac=paramJac_ODE_Sym)
ODE_test_Prob = ODEProblem(ODE_test,IC,tspan,p)
ODE_test_sol = DifferentialEquations.solve(ODE_test_Prob,Rodas5())#, alghint=:stiff)

This solves no problem -
And i get a nice solution fro ODE_test_sol
when I try to use adjoint sensitivity


# objective function is to obtain sensitivity with respect to tumor ratio u[19] I could change to tumor volume (u[16])

g(u,p,t) = (sum(u[19]).^2) ./ 2
function dg(out,u,p,t)
  # [out[i]= 0 for i in vcat(1:18,20:21)]
  # out[19]= u[19]
  [out[i]= u[19] for i in 1:varnum]
end

res = adjoint_sensitivities(ODE_test_sol, Rodas5(),g,nothing,dg)

when I tried with Vern9 or even DP8 on a server no result has shown up after 2 hours
My system has 21 equation (6 algebraic), 21 variables and 32 parameters
i appreciate any advice-

Thanks,

A

Thanks

PS- Once I am done with adjoint i am moving to sobol, but i was asked for sensitivity against u[19] explicitly

47

Adjoints will have the same issue with DDEs. We addressed this in our paper on this actually: [1812.01892] A Comparison of Automatic Differentiation and Continuous Sensitivity Analysis for Derivatives of Differential Equation Solutions . So we should just error on any mass matrix until we implement the extended sensitivity equation.

48

Thanks and thanks for the reference -
Iā€™ll read it today :slight_smile:

In the implemntation of Morris sensitivity do you have recomended values for len_trajectory ,num_trajectory ,total_num_trajectory (i see from the source code that defualts are 10 10 and 5*10 -
but I was wondering if those were just placeholders (in your example the numbers are much larger)

also if i indurtand correctly for everyparameter you need to provide an interval and number of samples would it be acceptable if I want to fix a parameter to provide something like this for param_range and steps

[[-1,1],[1,1],[5,5]],[10,10,1]

Thanks,

take care-

49

Hi, so you would want to set the len_trajectory and num_trajectory such that you have a decent spread over the parameter space but also taking into consideration the amount of computation you want to be doing as the number of simulations depends on the size of the design matrix generated which is what these arguments control. I hope the argument names are sufficiently clear to what these represent. I think the providing a range like [5,5] to keep the parameters should just work. Do let us know by opening an issue if it doesnā€™t, it should be a quick fix in that case. Iā€™ll also leave a link to the docs for Morris Method if you havenā€™t looked at them, maybe that will help in clearing more doubts http://docs.juliadiffeq.org/latest/analysis/global_sensitivity.html#Morris-Method-1. Please ping me if you have further doubts with GSA.

50

Thanks for the reply and the docs. Yes the anmes are pretty clear :slight_smile: just wanted to make sure I was thinking correctly - I checked the source code to figure things out, but never found the detailed docs, so thanks!

I already tried the [5,5] and [1] in number of steps it gave error -
I tried [5,5] and 1 for steps , also [5,5] and 0 for steps.
where should I raise the issue?

This is the error:

BoundsError: attempt to access 1-element StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}} at index [2]

Take care

51

Another question - about the range and steps -
Lets say I want to sample a parameter over a log scale -
Letā€™s say

p=1
p_range = [p/100, p*100]
psteps =11

# whihc I assume is equivalent to
collect(range(1.0/100, length=11, stop=1.0*100))

samples through

11-element Array{Float64,1}:
0.01
10.009
20.008
30.007
40.006
50.005
60.004
70.003
80.002
90.001
100.0

Is this correct?

If this is true consider that a logscale would be a better range to sample trhough the different
order of magnitud for the parameter in question-
for example-
something like this

10 .^range(log10(p/100), length=11, stop=log10(p*100))

samples through
11-element Array{Float64,1}:
0.01
0.025118864315095794
0.06309573444801933
0.15848931924611132
0.3981071705534972
1.0
2.51188643150958
6.309573444801933
15.848931924611133
39.810717055349734
100.0

Am I right that the parameter sapce is divided linearly and not logarithmically?

52

This is a completely new question. Youā€™ll likely find more help if you make a new topic to ask this question in. No DiffEq experience is needed to answer this so there may be a lot of bites. I am on a road trip and cannot give much attention.

1 Like
53

Iā€™ll send it around, I placed here because it directly deals with the DiffEqsensitivity.Morris inputs-
Have fun : ) and enjoy your trip!

A

54

Hi Chirs,
super quick question -
Is saving as dataframe supported for problems with mass matrices?
When I copy your code form the ODE page to save the linear function example you provide - it works perfectly. But, when I try to apply the same steps to my DAE (ODE with mass matrix) solution the df I get is transposed (iterations are columns and variables are rows and timestamp is nowhere to be found)

Any advice how can i include the time stamp in my df?

Thanks,

a

55

Yes

Example?

57

Feel free to share it privately as well. I wonā€™t get to it today though.

58

I actually donā€™t know how -
If you send me an address Iā€™ll be happy to send it there

Thanks,

a

59

I think I privately message you - Ill remove the post now -
Please let me know if you canā€™t find my message

A

60

Hi -
sorry to bother you again, but I still canā€™t figure out why DataFrame(sol) of ODEs will not include the time span- Any ideas how to fix that? any particular version of DataFrame or DifferentialEqautions I should use? My julia is 1.0 - I can always upgrade if you think that help.

Thanks,

A

61

Iā€™ll need an example