sciml / multiscalearrays.jl

A framework for developing multi-scale arrays for use in scientific machine learning (SciML) simulations

MultiScaleArrays

MultiScaleArrays.jl allows you to easily build multiple scale models which are fully compatible with native Julia scientific computing packages like DifferentialEquations.jl or Optim.jl. These models utilize a tree structure to describe phenomena of multiple scales, but the interface allows you to describe equations on different levels, using aggregations from lower levels to describe complex systems. Their structure allows for complex and dynamic models to be developed with only a small performance difference. In the end, they present themselves as an AbstractArray to standard solvers, allowing them to be used in place of a Vector in any appropriately made Julia package.

Tutorials and Documentation

For information on using the package, see the stable documentation. Use the in-development documentation for the version of the documentation, which contains the unreleased features.

Example

The usage is best described by an example. Here we build a hierarchy where Embryos contain Tissues which contain Populations which contain Cells, and the cells contain proteins whose concentrations are modeled as simply a vector of numbers (it can be anything linearly indexable).

using MultiScaleArrays
struct Cell{B} <: AbstractMultiScaleArrayLeaf{B}
    values::Vector{B}
end
struct Population{T <: AbstractMultiScaleArray, B <: Number} <: AbstractMultiScaleArray{B}
    nodes::Vector{T}
    values::Vector{B}
    end_idxs::Vector{Int}
end
struct Tissue{T <: AbstractMultiScaleArray, B <: Number} <: AbstractMultiScaleArray{B}
    nodes::Vector{T}
    values::Vector{B}
    end_idxs::Vector{Int}
end
struct Embryo{T <: AbstractMultiScaleArray, B <: Number} <: AbstractMultiScaleArrayHead{B}
    nodes::Vector{T}
    values::Vector{B}
    end_idxs::Vector{Int}
end

This setup defines a type structure which is both a tree and an array. A picture of a possible version is the following:

Let's build a version of this. Using the constructors we can directly construct leaf types: ```julia cell1 = Cell([1.0; 2.0; 3.0]) cell2 = Cell([4.0; 5.0]) ``` and build types higher up in the hierarchy by using the `constuct` method. The method is `construct(T::AbstractMultiScaleArray, nodes, values)`, though if `values` is not given it's taken to be empty. ```julia cell3 = Cell([3.0; 2.0; 5.0]) cell4 = Cell([4.0; 6.0]) population = construct(Population, deepcopy([cell1, cell3, cell4])) population2 = construct(Population, deepcopy([cell1, cell3, cell4])) population3 = construct(Population, deepcopy([cell1, cell3, cell4])) tissue1 = construct(Tissue, deepcopy([population, population2, population3])) # Make a Tissue from Populations tissue2 = construct(Tissue, deepcopy([population2, population, population3])) embryo = construct(Embryo, deepcopy([tissue1, tissue2])) # Make an embryo from Tissues ```

multiscalearrays.jl's Issues

MethodError: no method matching recursivecopy!(::Float64, ::Float64)

I get this error when I use the tuple_nodes.jl test example and change line 38 to

community = construct(Community, (deepcopy(plant1), deepcopy(plant2), ),rand(10))

i.e. adding a vector of values. Same if you do that for next level (scenario)

sidequestion: would it be possible to use namedTuples instead of an vector? I tried but wasn't sure how to specify the type parametrizations. I'll keep trying :-)

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Define y and eltype for multiscale model

See #5

Change Allocation Structure in ODEs/SDEs

ODEs/SDEs only need to update du using its indexing structure. The noise term also only needs to match in indices. After the MMM indexing change, straight usage of MMM's are essentially just as fast as arrays, but there is a difference inside of the diffeq solvers when allocations occur. This can be fixed by making items which are not user facing only match in linear indexing structure.

Thoughts

Hi Chris

Gitter is down atm (and github very slow) so I use an issue to write down stuff, hope thats ok.

1. The idea of multiscalemodels is brilliant. Your approach made it so much more general too.

2. I assume one can use e.g. space as Head node? for linear space its analog:

immutable linSpace{T<:AbstractMultiScaleModel} <: MultiScaleModelHead
  x::Vector{T}
  end_idxs::Vector{Int}
end

I guess the best (efficient) for 2D would be to just add 2 vectors holding the positions so you can calculate e.g. distances. Thats even general to handle both regular grids and random positions in 2d. I guess you could even make positions dynamic for mobile organisms...(makes me think if one could combine e.g. an ODE for e.g. cell growth/mobility with a PDE for diffusion of chemicals and have them interact....i.e. cells react to and affect chemical gradients)

3. can you mix leaf types? for example have predator and prey types that dispatch the growth dynamics function based on type?

4. parameters. I used parametrised functions to give parameters (e.g. traits) to the ode function. Will these AbstractMultiScaleModel types now hold these parameters?

5. when you apply the function to the leaf vector (Proteins in your example) would parameters in the parent node be available? more curious than any direct implementation in mind...Actually, if space is head node you could have environmental parameters such as temperature or nutrients there that determine growth rates of organisms living in that space-node. Actually, nutrients would also have to be dynamic but the value is not in a leaf node, unless you can have a nutrient leaf node mixed with others such as organisms...is this to confusing?

6. assume I have an species (population of) and I track abundance of individuals but I also want to use quantitative genetics to have a trait, e.g. mean size, evolve. I guess thats related to 3) as you'd want one abundance type and one size type as leaf for parent species and then have one function describing change in abundance and one describing change in mean size

This stuff has big potential :-)

Jon

Define similar

This is required for the solvers

Package compatibility caps

Ref: https://discourse.julialang.org/t/package-compatibility-caps/15301

Fast path for broadcasting with an array

When an array is in the broadcast, it falls back to a full indexing form. A better thing to do would be recursive views via @generated functions, but to get v0.7 working this is being skipped over for now.

no stable docs

add_node! not outputting expected result with integrator input

First off, this framework for multi-scale modeling is fantastic. A lot of the software I've tried does not have an adequate ODE solver suite for complex simulations.

I was writing a callback event that adds a cell to a population using the add_nodes! function and it added a cell with different values than what I gave it. My best guess for the cause is using the similar function in the diffeq.jl file. Here is a minimal working example of what I mean:

using DifferentialEquations, MultiScaleArrays

# Create Structures that define a Cell and Population
struct Cell{B} <: AbstractMultiScaleArrayLeaf{B}
    values::Vector{B}
end

struct Population{T<:AbstractMultiScaleArray,B<:Number} <: AbstractMultiScaleArrayHead{B}
    nodes::Vector{T}
    values::Vector{B}
    end_idxs::Vector{Int}
end


nCells=2
cellConstructor = [Cell([1.0, 2.0]) for i=1:nCells]
cellPop = construct(Population, deepcopy(cellConstructor))


function StateModel(dpop,pop,p,t) #This is arbitrary
    for (cell,dcell) in LevelIter(1,pop,dpop)
        dcell = 1.0
    end
end

prob = ODEProblem(StateModel,cellPop,(0.0,1.0))


integrator = init(prob,Tsit5())

add_node!(integrator,Cell([5.0, 20.0]))

integrator.u.nodes[end].values .== [5.0, 20.0] #Both false

Addat deleteat on a vector

It should fine the same indices and perform the addat/deleteat. That would greatly simplify the resizing code.

Segfault on 1.10 due to invalid getrs! inputs

It looks like MultiScaleArrays.jl is invoking getrs! (from ldiv!(::LU)) with inputs that cause Julia to crash: JuliaLang/julia#49489. Although Julia shouldn't crash, it seems like the inputs are invalid, and would cause issues with other BLAS implementations too.

This also involves LinearSolve and OrdinaryDiffEq, so I'll leave it up to the SciML devs to further triage this issue. Details can be found in the Julia issue linked above.

Stiff solvers cannot resize

This is broken on the latest version. The design of ForwardDiff causes a problem here since the configs are immutable and we need to modify the duals in there, so it needs to propagate down a mutation call to each dual cache.

Define easier "level looping"

You should be able to choose a level to loop at, like at the population level or the tissue level.

Solving a MultiScaleArrays which leafs contain additional information (e.g. celltype::Symbol)

Hi,
First thanks a lot for developing this! This is just amazing.
I am trying to develop a model of symbiosis similar to the "embryo" one you developed.
I wrote a "MultiScaleArrays" describing my model, and a function computing the derivative.
When I try to solve it, I get this error message:
ERROR: MethodError: Cannot convert an object of type Array{Float64,1} to an object of type CellGenotype CellGenotype being my leaftype
This may have arisen from a call to the constructor CellGenotype(...), since type constructors fall back to convert methods.

The versions I am running are
Julia : Version 0.6.0 (2017-06-19 13:05 UTC) on x86_64-pc-linux-gnu
DifferentialEquations : v2.2.1
MultiScaleArrays : v0.4.0

Below there is a very simplified version of my script allowing to repeat the error:

using DifferentialEquations
using MultiScaleArrays

struct CellGenotype{B<:Float64} <: AbstractMultiScaleArrayLeaf{B}
    Iam::Symbol
    values::Vector{B} 
end

struct Host{T<:AbstractMultiScaleArray,B<:Float64} <: AbstractMultiScaleArrayHead{B}
    nodes::Vector{T}
    values::Vector{B}
    end_idxs::Vector{Int}
end

host = construct(Host,[
        CellGenotype(# Genotype present at start for symbiont 1
        :aSymbiontGenotype, # I am
        [Float64(1)],
        ), 
        CellGenotype(# Genotype present at start for symbiont 1
        :aSymbiontGenotype, # I am
        [Float64(1)],
        )]
    )

evolve = function (t,u::Host,du)
  for i in 1:length(u)
    du[i] = 2*u[i]
  end
end

## Solve
tspan = (0.0,50)
prob = ODEProblem(evolve, host, tspan)
sol = solve(prob)

The problem disappears if I remove the field Iam::Symbol of the leaf CellGenotype. But my model being much more complex than what is shown here, I really need to have these kind of additional fields in the leafs.
The use of these additional fields is described here:

struct Cell{B} <: AbstractMultiScaleArrayLeaf{B}
    values::Vector{B}
    celltype::Symbol
end

But this tutorial (here) cruelly lack the solving step.
((Some additional minor comments on the tutorial:

• I think all the level_iter(embryo,2) should be replaced by level_iter(embryo,3)
• The variable d_embryo is not explained
• Close to d_embryo, function f is not explained
• I think the sentence

Since em will be the "vector" for the differential equation or otimization problem, it will be the value passed to the event handling.
should be
Since emBRYO will be the "vector" for the differential equation or oPtimization problem, it will be the value passed to the event handling.

• I think the code
  tissue3 = construct(Tissue, deepcopy([population; population2]))
  should be

                                               \|/
tissue3 = construct(Tissue, deepcopy([population, population2]))

))

Resizing Macros for DifferentialEquations.jl

Resizing macros should be made to resize the caches in events with DifferentialEquations.jl since the interface is slightly different. This just amounts to making new macros where the name of the function that's called in the loop is changed.

MultiScaleModels Profiling

I optimized the ODE solvers around these guys a bit. As a test case, I took a linear system on the embryo from the tests. The embryo is 4 layers and a total length of 13, and so it uses very small leaf nodes. This means that, with a cheap problem and a high hierarchy, it should accentuate any differences and give an upper bound on the cost of using a MultiScaleModel.

I then benchmarked the same problem solved with MultiScaleModels and then just having it be a regular length 13 contiguous array:

prob = ODEProblem(f,em,(0.0,1500.0))
@benchmark sol1 = solve(prob,Tsit5(),save_timeseries=false)
prob = ODEProblem(f,em[:],(0.0,1500.0))
@benchmark sol2 = solve(prob,Tsit5(),save_timeseries=false)

to get the results:

BenchmarkTools.Trial: 
  memory estimate:  230.28 kb
  allocs estimate:  10920
  --------------
  minimum time:     16.469 ms (0.00% GC)
  median time:      16.942 ms (0.00% GC)
  mean time:        17.033 ms (0.46% GC)
  maximum time:     26.385 ms (33.77% GC)
  --------------
  samples:          294
  evals/sample:     1
  time tolerance:   5.00%
  memory tolerance: 1.00%

BenchmarkTools.Trial: 
  memory estimate:  196.37 kb
  allocs estimate:  10409
  --------------
  minimum time:     7.500 ms (0.00% GC)
  median time:      7.846 ms (0.00% GC)
  mean time:        7.907 ms (0.51% GC)
  maximum time:     13.730 ms (42.01% GC)
  --------------
  samples:          632
  evals/sample:     1
  time tolerance:   5.00%
  memory tolerance: 1.00%

This puts an upper bound for the cost of using MultiScaleModels at just over 2x. But this is without saving: saving the timeseries is more costly with a MMM. Measuring the cost with saving every step:

prob = ODEProblem(f,em,(0.0,1500.0))
@benchmark sol1 = solve(prob,Tsit5())
prob = ODEProblem(f,em[:],(0.0,1500.0))
@benchmark sol2 = solve(prob,Tsit5())

BenchmarkTools.Trial: 
  memory estimate:  3.20 mb
  allocs estimate:  43891
  --------------
  minimum time:     26.454 ms (0.00% GC)
  median time:      27.738 ms (0.00% GC)
  mean time:        28.962 ms (4.21% GC)
  maximum time:     39.190 ms (27.02% GC)
  --------------
  samples:          173
  evals/sample:     1
  time tolerance:   5.00%
  memory tolerance: 1.00%

BenchmarkTools.Trial: 
  memory estimate:  1.52 mb
  allocs estimate:  21060
  --------------
  minimum time:     9.414 ms (0.00% GC)
  median time:      9.910 ms (0.00% GC)
  mean time:        10.252 ms (3.11% GC)
  maximum time:     18.033 ms (35.72% GC)
  --------------
  samples:          488
  evals/sample:     1
  time tolerance:   5.00%
  memory tolerance: 1.00%

This puts an upper bound around 3x. So the maximal cost of the abstraction is about 2x-3x for ODEs. It's actually lower for SDEs and DDEs because more of the calculation in those domains are spent on the interpolation and noise generation, which are actually able to mostly avoid the costs of MMMs. So the final product is something where the abstract cost is less than the performance difference between OrdinaryDiffEq and other packages, meaning using MMMs in OrdinaryDiffEq should still be slightly faster than using other packages with contiguous arrays. I think this counts as well optimized, and so my goal will be to now make this all compatible with the solvers for stiff equations since will make a large impact.

@zahachtah I think you might be interested in these results.

Efficient broadcast iteration

Multidimensional indexing

em[1,3,2]

as in the README needs to be implemented. That should be pretty trivial using recursion, but good tests are needed which will probably take more time.

All resize testing is disabled and broken

Repeat allocation

function similar(m::AbstractMultiScaleModel)
  m_new = construct(typeof(m),deepcopy(m.x),copy(m.y)) # Copy because y is a vector!
  for i in eachindex(m.x) # Is this necessary?
    m_new.x[i] = similar(m.x[i])
  end
  m_new
end

Pictures

Indexing Issues with MultiScaleArrays.jl on v0.7

I'm trying to follow along with the tissue hierarchy of an embryo example in the README and I'm running into some issues when I try to call:

embryo[2,3,1] # Gives the 1st cell in the 3rd population of the second tissue

I get the following error:

BoundsError: attempt to access 2-element Array{Population{Cell{Float64},Float64},1} at index [3]

What is the correct way to reference an individual element? And are there any special syntax considerations when using the MultiScaleArrays with DifferentialEquations/Optim?

The rest of the code seems to work (i.e. adding a node to embryo, etc). Attached is a copy of the code I used to implement the example.

MultiScaleArrays_Test_2018-10-01.pdf

Package causing "unsatisfiable package requirements" error

Getting the following error (I think because I'm using DifferentialEquations.jl).

Don't know how to resolve it 😬

---

ERROR: unsatisfiable package requirements detected: no feasible version could be found for package: MultiScaleArrays
  (you may try increasing the value of the
   JULIA_PKGRESOLVE_ACCURACY environment variable)

Does this package overlap with HierarchicalMatrices.jl?

This package seems to be very similar to HierarchicalMatrices.jl. Maybe it makes sense to unify them?

@MikaelSlevinsky

ERROR: UndefVarError: dims not defined using custom fields

when using some custom fields in the structs I get an error using similar() and therefore when trying to use DifferentialEquations with f(du,u....).

e.g. using the example from readme but changing:

struct Embryo{T<:AbstractMultiScaleArray,B<:Number} <: AbstractMultiScaleArrayHead{B}
           nodes::Vector{T}
           values::Vector{B}
           end_idxs::Vector{Int}
           x::Float64 # <- change
       end

and constructing works fine but similar(embryo) gives

ERROR: UndefVarError: dims not defined

Register

@JuliaRegistrator register()

Improve docstrings

Such that makedocs warnonly = [:missing_docs] can be removed.