Does this package overlap with HierarchicalMatrices.jl? about multiscalearrays.jl

[ChrisRackauckas]

Yes, it's basically a hierarchical vector with a few things:

1. The .y nodes are outside of the hierarchy at each level and are used to store things like cell-specific morphogen gradients (current PDE solutions).
2. Broadcasted operations are quite efficient through recursion: https://github.com/JuliaDiffEq/MultiScaleArrays.jl/blob/master/src/math.jl#L19-L42
3. Resizing is straightforward and cheap https://github.com/JuliaDiffEq/MultiScaleArrays.jl/blob/master/src/diffeq.jl
4. The nodes at the bottom can all be different sizes (for different cell types)
5. You can add new fields to the nodes to be able to signify cell types and track cell lineages.

It's basically for building beefed up simulations of http://docs.juliadiffeq.org/latest/features/callback_functions.html#Example-3:-Growing-Cell-Population-1 where you have entire gene regulatory networks, different cell types, and morphogen gradients across space. I built it for modeling single cell RNA-seq data, so it needs to satisfy those 5 points. I'd be surprised if anything else really hits this dead on because it's a very specific list of requirements for a specific modeling application.

[ChrisRackauckas]

It would be fun to try it with an n-body simulation (which I think is an application of fast multipole? And HierarchicalMatrix solvers == Fast multipole)

Let's take that over to DiffEqPhysics. There might be a GSoC for n-body problems so this would be great to handle.

[ChrisRackauckas]

Does it work inside the ODE and SDE solvers, and does have easy resizing?

[dlfivefifty]

What do you mean by easy resizing?

[ChrisRackauckas]

tissue3 = construct(Tissue, deepcopy([population, population2]))
add_node!(embryo, tissue3) # Adds a new tissue to the embryo
remove_node!(embryo, 2, 1) # Removes population 1 from tissue 2 of the embryo

The main purpose of this package is to be able to add and remove portions of the array with minimal cost.

[dlfivefifty]

On the matrix level what does add_node! do? Adds additional rows/columns?

[ChrisRackauckas]

This is an abstract vector. If it was a matrix it would be a ragged array in general, but it would be like adding rows and columns

[MikaelSlevinsky]

HierarchicalMatrices.jl is meant to create matrices with heterogeneous leaf types (cell types in MultiScaleArrays.jl). Also, in HierarchicalMatrices.jl, the nodes can support further hierarchical nodes and leaves themselves. The main application (so far) is a fast BLAS-overridden A_mul_B! in the 0.6 syntax, but fast linear solves should be no problem either.

It can modify the block types on-the-fly (provided the new type was included in the parse-time list), but it can't modify the number of blocks on-the-fly as it stores a concrete matrix of each type as a placeholder.

[dlfivefifty]

Ahhh.... so it's more of a HierarchicalVector?

https://github.com/JuliaApproximation/SingularIntegralEquations.jl/blob/aa6426f512fb826530be862f12ef6ac9aeb91b14/src/LinearAlgebra/HierarchicalVector.jl

[MikaelSlevinsky]

I mean, it could "remove" a block by replacing it with a type of zero array (or just change the referenced type of block number in the hash table assigned.)

[dlfivefifty]

Does it work inside the ODE and SDE solvers

It would be fun to try it with an n-body simulation (which I think is an application of fast multipole? And HierarchicalMatrix solvers == Fast multipole)

[MikaelSlevinsky]

Yes, that would be fantastic. The example shows usage for a couple different potentials.

[MikaelSlevinsky]

so explicit methods work immediately

[dlfivefifty]

Does this package overlap with HierarchicalMatrices.jl?
it can't modify the number of blocks on-the-fly

So it sounds like the answer is "yes it overlaps, but not 100%".

So the question is whether there is enough overlap to try to make a common interface.

BlockArrays.jl is another overlapping package. I'd suggest at the very least both HierarchicalMatrices.jl and MultiScaleArrays.jl support Block from BlockArrays.jl, so that you can dig into the hierarchy via (MultiScaleArrays.jl) A[Block(1)][Block(3)][Block(2)] and (HierarchicalMatrices.jl) A[Block(1,1)][Block(2,3)][Block(3,4)].

[MikaelSlevinsky]

P.S. How many levels deep do you usually go? I've found that indexing usually slows down past about 10-12 levels, curious what you find.

[ChrisRackauckas]

Usually just 3.