Applying GSA to Data Model

[lhdp0110]

I am trying to run a PRCC and an eFAST GSA on a data model that is written using the InformationGeometry.jl package. I am a beginner in both data modeling and in Julia, so I may just be misunderstanding how to implement the gsa methods, but I can't figure out what the proper syntax would be. The code below runs up until it reaches the gsa command, where it throws an error saying that there is no method matching what I have written. I would greatly appreciate any help with this!

using InformationGeometry, OrdinaryDiffEq, Plots, CSV, DataFrames
using DiffEqCallbacks, GlobalSensitivity, Distributions

function pathwayEq!(du,u,p,t)
    DP,cN,cNstar,phiR,deltaA = p

    V = 1
    VN = 0.1
    m = 2
    k = 2

    cP = 0.1                                 
    cI = 0.02                                
    cB = 0.044    
    cH = 0.05
    cA = 0.041
    cM = 0.00001              

    phiP = 15                     
    phiI = 50                      
    phiM = 10                            

    deltaP = 0.01 
    deltaPstar = 0.01
    deltaR = 0.01
    deltaN = 0.01
    deltaH = 0.05
    
    rhoIstar = 0.01
    KR = 2
    KM = 10    
    IT = 1                                            
    QM = 1

    QP = deltaP
    QR = deltaR/V
    QH = deltaH/V
    QA = deltaA/V
    
    du[1] = QP - m * cP * exp(-phiP*u[8]) * u[1]^k * u[10] - deltaP * u[1]    # dP/dt
    du[2] = cP * exp(-phiP*u[8]) * u[1]^k * u[10] - deltaPstar * u[2]         # dPstar/dt
    du[3] = cI * exp(-phiI*u[8]) * (IT - u[3]) * u[2] - rhoIstar * u[3]       # dIstar/dt
    du[4] = QR / V + (phiR*u[7])/(KR+u[7]) - cB * u[3] * u[4] - deltaR * u[4] # dRB/dt
    du[5] = cB * u[3] * u[4] - DP * u[5] - deltaR * u[5]                      # dRP/dt
    du[6] = DP * u[5] / VN - cN * u[6] + cNstar * u[7] - deltaN * u[6]        # dRN/dt
    du[7] = cN * u[6] - cNstar * u[7]                                         # dRNstar/dt
    du[8] = QH / V + cH * u[7] - deltaH * u[8]                                # dH/dt
    du[9] = QA / V + cA * u[7] - deltaA * u[9]                                # dA/dt
    du[10] = QM * u[10] + (phiM *u[9]/(KM+u[9])) * u[10] - cM * u[10]^2       # DG/dt
end

P0 = 1
Pstar0 = 0.1
Istar0 = 0.1
RB0 = 1
RP0 = 0.1
RN0 = 0.1
RNstar0 = 0.1
H0 = 1
A0 = 0.48
M0 = 1
u0 = [P0;Pstar0;Istar0;RB0;RP0;RN0;RNstar0;H0;A0;M0]

AMPs = [0.475,4.85,3.32,4.66,4.47,2.37,1.14]
confInts = [1.42,2.25,2.80,1.86,2.20,2.32,1.95]

hours = [0, 120, 240, 480, 720, 1440, 2880]
IMDdataSet = DataSet(hours, AMPs, confInts; xnames=["Time (seconds)"], ynames=["AMPs"])    
params = [5.61,9.60,1.96,0.0707,0.0826] 

@time IMDdataModel = DataModel(IMDdataSet, ODEFunction(pathwayEq!), u0, [9], params; tol=0.1)   

plot(IMDdataModel)

f1 = function(params)
    prob = remake(IMDdataModel;params=params)
    sol = solve(prob, Tsit5())
  end;

bounds = [[4,22],[8,14],[1,3],[0,1],[0,1]]

m_regression = gsa(f1, RegressionGSA(rank=true), bounds)
m_regression.partial_correlation
m_regression.partial_rank_correlation

m_efast = gsa(f1, eFAST(), bounds)
m_efast.S1
m_efast.ST ```

[Vaibhavdixit02]

I tried running your model -

DP,cN,cNstar,phiR,deltaA = p

expects p to have 5 elements but params = [5.61,9.60,0.0707,0.0826] only has 4, can you clarify this?

[lhdp0110]

Yes this is a typo sorry, I just fixed it

[Vaibhavdixit02]

I am still running into an error while trying to run the model, can you please confirm whether it works

julia> @time IMDdataModel = DataModel(IMDdataSet, ODEFunction(pathwayEq!), u0, [9], params; tol=0.1)
ERROR: MethodError: no method matching lu!(::StaticArraysCore.MArray{Tuple{10, 10}, Float64, 2, 100}, ::LinearAlgebra.RowMaximum)
Closest candidates are:
  lu!(::StridedMatrix{T}, ::LinearAlgebra.RowMaximum; check) where T<:Union{Float32, Float64, ComplexF32, ComplexF64} at /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/lu.jl:80
  lu!(::Union{LinearAlgebra.Hermitian{T, S}, LinearAlgebra.Symmetric{T, S}} where {T, S}, ::Union{LinearAlgebra.NoPivot, LinearAlgebra.RowMaximum}; check) at /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/lu.jl:88
  lu!(::StridedMatrix{T} where T, ::Union{LinearAlgebra.NoPivot, LinearAlgebra.RowMaximum}; check) at /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/lu.jl:134
  ...
Stacktrace:
  [1] do_factorization(alg::LinearSolve.LUFactorization{LinearAlgebra.RowMaximum}, A::StaticArraysCore.MArray{Tuple{10, 10}, Float64, 2, 100}, b::StaticArraysCore.MArray{Tuple{10}, Float64, 1, 10}, u::StaticArraysCore.MArray{Tuple{10}, Float64, 1, 10})
    @ LinearSolve ~/.julia/packages/LinearSolve/WhI1M/src/factorization.jl:

[lhdp0110]

Yes, that is the same error I am getting; to clarify, the model does work (I just checked) and produces a plot when you run it only up to line 73 plot(IMDdataModel)
The blocks that follow the line above are the ones I am having an issue with

[Vaibhavdixit02]

I got it at the line before it, on running @time IMDdataModel = DataModel(IMDdataSet, ODEFunction(pathwayEq!), u0, [9], params; tol=0.1)