Typing in Julia

I ran the following code to create a vector of initial conditions:

IC_V = collect(-75.0:15.0:75.0);
IC_n = collect(0:.2:1);
IC_gCa = collect(0:p[6]/5:p[6]);
IC_gK = collect(0:p[7]/5:p[7]);
IC_Ca = collect(0.0:50.0:300.0);
IC_nHalf = collect(0:p[8]/5:p[8]);

AllICs = collect(Iterators.product(IC_V,IC_n,IC_gCa,IC_gK,IC_Ca,IC_nHalf));

I want to create an Ensemble Problem by constructing an ODEproblem with different initial conditions. I think there’s a type problem though. AllICs[1] returns a NTuple{6, Float64} but I need the initial conditions as Matrix{Float64}. Naively, I imagine the following should work, but it does not: convert(Matrix{Float64},AllICs[1]). How do I convert my tuple into what I need?

Sorry, I’m not used to dealing with typing this stuff – Matlab hides this all :smiley: . More generically, are there any quick tutorials on managing these numerical issues (Julia or otherwise)?

Can you provide the definition of p? I’m just using a placeholder: p = rand(0:9, 8), for now.

Firstly, make sure to not use collect where you don’t need them. There is not reason to collect your initial ranges, so do:

IC_V = -75.0:15.0:75.0;
IC_n = 0:.2:1;
IC_gCa = 0:p[6]/5:p[6];
IC_gK = 0:p[7]/5:p[7];
IC_Ca = 0.0:50.0:300.0;
IC_nHalf = 0:p[8]/5:p[8];

In general, don’t use collect at all, unless you really need it, which I guess you do in the next step.

The product iterator returns tuples, but tuples are turned into vectors by collect itself, so you can just collect each element, like this (notice the dot):

AllICs = collect.(Iterators.product(IC_V, IC_n, IC_gCa, IC_gK, IC_Ca, IC_nHalf));

But, why do you need them as vectors instead of tuples, and do you really need to collect this into AllICs? Most of the time it’s better to not collect, and instead just iterate and do your calculations one-by-one.

At least, ask yourself this question. Can you just avoid the collection, and can you use tuples instead of vectors?

2 Likes

Ah sorry, p is just a vector of parameters. 'p[6] = 3, 'p[7] = 6, and p[8] = 20.

To answer the question why I’m collecting all my initial condition:

The documentation for how to use the ensemble problem is here. In particular one needs to define an EnsembleProblem which is created using a constructor. prob_func is part of this constructor – it takes the originally defined ODE and reformulates it. I was going to define all of my intial conditions beforehand, then index into the array of intial conditions using i in prob_func. Full code below:

using DifferentialEquations

#################################################
##### ODE params

# initialize parameter vector
p = fill(NaN,8);

# Time Constants for Homeostasis -- (g_Ca, g_K, nHalf)
p[1] = 500; p[2] = 500; p[3] = 500;

# C_T & Delta
p[4] = 20; p[5] = 5; 

# Sigmoid Ranges -- (g_Ca, g_K, nHalf)
p[6] = 3; p[7] = 6; p[8] = 20;

# Calcium Handling
A = -1; k = -1/.1;

I_app = 0;
C = 1; g_L = .5;
E_L = -50; E_Ca = 100; E_K = -70;

#################################################
## Initial Condition Testing -- Tests for Multistability

IC_V = collect(-75.0:15.0:75.0);
IC_n = collect(0:.2:1);
IC_gCa = collect(0:p[6]/5:p[6]);
IC_gK = collect(0:p[7]/5:p[7]);
IC_Ca = collect(0.0:50.0:300.0);
IC_nHalf = collect(0:p[8]/5:p[8]);

AllICs = collect(Iterators.product(IC_V,IC_n,IC_gCa,IC_gK,IC_Ca,IC_nHalf));

#################################################

# Define functions
sigmoid(x) = 1/(1+exp(-x));
tau(V, halfN) = 3/cosh(((V-halfN)-10)/29);
Ninf(V, halfN) = sigmoid(((V-halfN)-10)/7.25);

# Define Currents
Ica(g_Ca,V) = g_Ca*(sigmoid((V+1)/7.5)+.1)*(V-E_Ca);
Ik(g_K,n,V) = g_K*n*(V-E_K);
Il(V) = g_L*(V-E_L);

# u[1] - V cell 1, u[2] - slow var cell, u[3] - g_Ca, u[4] - g_K, u[5] - [Ca], u[6] - halfN
# p[1] - tau_L.gCa , p[2] - tau_L.gK, p[3] - tau_L.nHalf, p[4] - C_T, p[5] - Delta, p[6] - g_Ca, p[7] - g_K, p[8] - nHalf
function f(du,u,p,t)
    du[1] = (-1*Il(u[1]) + -1*Ica(u[3],u[1]) + -1*Ik(u[4],u[2],u[1]) + I_app)/C
    du[2] = (Ninf(u[1],u[6])-u[2])/tau(u[1],u[6])
    du[3] = (p[6]*sigmoid((p[4]-u[5])/p[5]) - u[3])/p[1]
    du[4] = (p[7]*sigmoid((u[5]-p[4])/p[5]) - u[4])/p[2]
    du[5] = -1*k*(A*Ica(u[3],u[1]) - u[5])
    du[6] = (p[8]*(sigmoid((p[4]-u[5])/p[5])) - u[6])/p[3]
end

tspan = (0.0f0,8000.0f0);
prob = ODEProblem(f,ICs,tspan,p);
function prob_func(prob,i,repeat)
  remake(prob,u0=AllICs[i])
end
ensemble_prob = EnsembleProblem(prob, prob_func = prob_func, safetycopy=false)
@time sim = solve(ensemble_prob,Tsit5(),EnsembleThreads(),trajectories=1000)

Maybe there’s a better way to do this? I don’t know.

If u0 is required to be a vector here (which seems to be the case based on your description), you can change this to

remake(prob,u0=collect(AllICs[i]))

to change the tuples to vectors on the fly. This will avoid having to collect all the elements in one go (leading to a lot of memory usage), but give similar results.

1 Like

But I agree with @DNF that all the other collects in these lines:

seem unnecessary. Just IC_V = -75.0:15.0:75.0 will work for this usage, and similarly for the other lines. Generally, as a good practice, you should avoid collects unless you get an error about needing a Vector or Array (as you presumably did in this case). Generally it shouldn’t be necessary.

1 Like