Writing Mat Files to Disk

New to Julia
1

I’m trying to save the outputs for my ODE simulations into a file I can open in MATLAB for analysis. I’m using the MAT library. I’m obtaining the following error (with stacktrace). The code that produce the error follows it. I’m not sure how to interpret this error. Ideas?

ERROR: LoadError: This is the write function for CompositeKind, but the input doesn't fit
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] m_write(mfile::MAT.MAT_HDF5.MatlabHDF5File, parent::HDF5.Group, name::String, s::Nothing)
    @ MAT.MAT_HDF5 ~/.julia/packages/MAT/wBZA8/src/MAT_HDF5.jl:532
  [3] m_write(mfile::MAT.MAT_HDF5.MatlabHDF5File, parent::HDF5.Group, name::String, k::Vector{String}, v::Vector{Any})
    @ MAT.MAT_HDF5 ~/.julia/packages/MAT/wBZA8/src/MAT_HDF5.jl:520
  [4] m_write(mfile::MAT.MAT_HDF5.MatlabHDF5File, parent::HDF5.Group, name::String, s::ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float32}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float32, Float32}, true, Vector{Float64}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.FullSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float32}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float32}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats})
    @ MAT.MAT_HDF5 ~/.julia/packages/MAT/wBZA8/src/MAT_HDF5.jl:535
  [5] m_write(mfile::MAT.MAT_HDF5.MatlabHDF5File, parent::HDF5.Group, name::String, data::Vector{ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float32}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float32, Float32}, true, Vector{Float64}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.FullSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float32}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float32}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}})
    @ MAT.MAT_HDF5 ~/.julia/packages/MAT/wBZA8/src/MAT_HDF5.jl:481
  [6] m_write(mfile::MAT.MAT_HDF5.MatlabHDF5File, parent::HDF5.File, name::String, k::Vector{String}, v::Vector{Any})
    @ MAT.MAT_HDF5 ~/.julia/packages/MAT/wBZA8/src/MAT_HDF5.jl:520
  [7] m_write(mfile::MAT.MAT_HDF5.MatlabHDF5File, parent::HDF5.File, name::String, s::EnsembleSolution{Float64, 3, Vector{ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float32}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float32, Float32}, true, Vector{Float64}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.FullSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float32}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float32}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}}})
    @ MAT.MAT_HDF5 ~/.julia/packages/MAT/wBZA8/src/MAT_HDF5.jl:535
  [8] write(parent::MAT.MAT_HDF5.MatlabHDF5File, name::String, thing::EnsembleSolution{Float64, 3, Vector{ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float32}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float32, Float32}, true, Vector{Float64}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.FullSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float32}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float32}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}}})
    @ MAT.MAT_HDF5 ~/.julia/packages/MAT/wBZA8/src/MAT_HDF5.jl:548
  [9] top-level scope
    @ ~/Dropbox/work/marder/test2.jl:77
 [10] include(fname::String)
    @ Base.MainInclude ./client.jl:476
 [11] top-level scope
    @ REPL[1]:1

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 = -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];

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

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

using MAT
outFile = matopen("sims.mat","w");
write(outFile,"sims",sim);
close(outFile);
2

You might want to just use HDF5 directly via HDF5.jl. HDF5 support in MATLAB is decent these days. As you may notice here, v7.3 MAT files are just a particular form of HDF5 file.

3

I think the problem is the output type of the simulation, based on the error. This is the type of variable that the output simulation are:

EnsembleSolution Solution of length 1000 with uType:
ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float32}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float32, Float32}, true, Vector{Float64}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.FullSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float32}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float32}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}

Are there functions that extract the fields that are available in this object? I’m having trouble navigating and indexing into the object.

1 Like
4

you are right. this is the problem. what data do you want to extract?

5

well i just want time and the 6 dimensions of u. but im curious now about what else is being returned.

6

part of the reason it’s not a 6d array is that every simulation in the ensemble will have different time steps. as such, they may have different lengths. this basically stores for each simulation, the parameters of the problem, the time steps, the u values, and probably a bunch more that I’m forgetting about.

7

The output from DifferentialEquation.jl simulations is a structure of interpolation functions. Does MATLAB support such structures in *.mat files? And can MATLAB understand the resulting interpolation functions? [I haven’t used MATLAB in ages…]

If MATLAB does not support the interpolation structures of DifferentialEquation.jl, you will have to convert the structures to matrices of data before saving them into *.mat files. This is pretty straighttforward.

8

Ya I forced solve to create outputs at particular times using the saveat parameter. The output is still in a strange format, I had to loop through all the data and manually load it into an Array which I then saved.

9

The saveat option does not make the result into a matrix. The solution will still be a structure of interpolation functions.

There is no need to “loop through all the data” and “manually load it into an Array”. Given that your solution is named sim, you will find that sim.t provides a vector of the time instances for the solution, while Array(sim) gives a matrix of the solutions in the time instances of sim.t.