cesaraustralia / DynamicGrids.jl

Grid-based simulations in Julia
Other
225 stars 7 forks source link

Compatibility with complex object such as YAXArrays #254

Open remylegoff opened 1 month ago

remylegoff commented 1 month ago

I want use dispersion in a complex macroevolution simulation. To do so, I use YAXArrays to handle the high dimensionnality of the data in a spatial context. Then, there is a version compatibility issue between DynamicGrids and YAXArrays. DynamicGrids only work on JULIA version 1.9 whereas YAXArrays need JULIA version 1.10. Then DynamicGrids only take simple object as input which reduce the application possibility for complex scenari. Is there any development plan for this kind of issue ?

rafaqz commented 1 month ago

Yeah, check out the dev version of DynamicGrids.jl for now.

But I'm not sure what you mean by simple and complex objects ? Like mathematical complex? That should totally work.

As for YAX I pretty much always use some kind of AbstractDimArray like YAX gives you in DynamicGrids.jl. And usually with some complicated object in the grid rather than just floating point values. But you may need to move the data to memory as YAX is usually disk based.

Maybe be a bit more specific about what simple and complex means to you and what exactly is breaking.

remylegoff commented 1 month ago

By simple object, I mean that tuple and matrix are common and from base JULIA. I name complex object such as YAXArrays or raster that are build from the simple object and add features useful in specific task.

rafaqz commented 1 month ago

I basically always use a Raster! Rasters.jl was originally written for this use-case.

(so maybe post some code that doesn't work and we can workshop it)

remylegoff commented 2 weeks ago

The issue I got is that I need to interact between multiple and growing number of arrays. I'm developing a mechanistic spatial evolutionary model with a growing number of species overtime and interaction between these species. And I can't see how to make it work by using YAXArrays to keep all species in the same object to be able to pick either all cell of a species or all species in a cell

But for now I have an error UndefVarError: `positions` not defined

rafaqz commented 2 weeks ago

By using YAX I think you just mean DimensionalData.jl ? Rather than YAX disk-based tools? YAX dimensions are just DD dimensions extended for Zarr and similar. They are supported here, but not how you think.

You cant have a dimension of an array in a simulation, it breaks the abstraction. It will just try to run in 3d! The init array defines the dimensionality of the space. I guess we could actually allow those to not be the same thing and pass in a view over the extra dimensions?

But you wont be easily able to grow the array to add new species anyway (julia multidimensional arrays don't change size)... so the actual efficient (and reccomended) way to do this is use some kind of StaticArray. So StaticArrays.jl, and LabelledArrays.jl has a labelled version.

Then you allocated the maximum length vector up-front and track how many things are in it as a run-time variable. So make a struct with a StaticArray in it, and a number saying how many things it holds.

Like this:

struct Species{A}
    species::A
    nspecies::Int
end
Base.zero(::Species{A,N}) = Species(zero(species.species), 0)

I use a StaticArray like this for about 50 species in some simulations. If you need 1000, yes you will need something else, likely not DynamicGrids.jl in its current form.

(Additionally if you want to post errors please post MWEs and the full stack trace, then I can actually help)

remylegoff commented 2 weeks ago

Sorry for the error here is MWE :

using Pkg
Pkg.activate("chap3")
Pkg.instantiate()
using DynamicGrids, ColorSchemes, Colors, BenchmarkTools

const DEAD, ALIVE, BURNING = 1, 2, 3
init = fill(ALIVE, 400, 400)
bench_output = ResultOutput(init; tspan=1:200)
setneighbors_rule = let prob_combustion=0.0001, prob_regrowth=0.01
    SetNeighbors(Moore(1)) do data, neighborhood, cell, I
        if cell == DEAD
            if rand() <= prob_regrowth
                data[I...] = ALIVE
            end
        elseif cell == BURNING
            for pos in positions(neighborhood, I)
                if data[pos...] == ALIVE
                    data[pos...] = BURNING
                end
            end
            data[I...] = DEAD
        elseif cell == ALIVE
            if rand() <= prob_combustion 
                data[I...] = BURNING
            end
        end
    end
end
sim!(bench_output, setneighbors_rule)

and the traceback:

ERROR: UndefVarError: `positions` not defined
Stacktrace:
  [1] (::var"#13#14"{…})(data::DynamicGrids.RuleData{…}, neighborhood::Moore{…}, cell::Int64, I::Tuple{…})
    @ Main ~/Documents/PhD/Chapitre_3/Scripts/Gene3sis/Julia-implementation/genesis_julia/dispersion/mwe_position_issue.jl:16
  [2] applyrule!
    @ ~/.julia/packages/DynamicGrids/IIwvg/src/rules.jl:624 [inlined]
  [3] cell_kernel!
    @ ~/.julia/packages/DynamicGrids/IIwvg/src/maprules.jl:240 [inlined]
  [4] maprule!
    @ ~/.julia/packages/DynamicGrids/IIwvg/src/maprules.jl:89 [inlined]
  [5] maprule!
    @ ~/.julia/packages/DynamicGrids/IIwvg/src/maprules.jl:68 [inlined]
  [6] maprule!(simdata::DynamicGrids.SimData{…}, ruletype::Val{…}, rule::SetNeighbors{…})
    @ DynamicGrids ~/.julia/packages/DynamicGrids/IIwvg/src/maprules.jl:46
  [7] maprule!
    @ ~/.julia/packages/DynamicGrids/IIwvg/src/maprules.jl:10 [inlined]
  [8] sequencerules!
    @ ~/.julia/packages/DynamicGrids/IIwvg/src/sequencerules.jl:16 [inlined]
  [9] sequencerules!
    @ ~/.julia/packages/DynamicGrids/IIwvg/src/sequencerules.jl:5 [inlined]
 [10] |>
    @ ./operators.jl:917 [inlined]
 [11] _step!(sd::DynamicGrids.SimData{…}, rules::Tuple{…})
    @ DynamicGrids ~/.julia/packages/DynamicGrids/IIwvg/src/framework.jl:208
 [12] simloop!(output::ResultOutput{…}, simdata::DynamicGrids.SimData{…}, ruleset::Ruleset{…}, fspan::UnitRange{…}; printframe::Bool)
    @ DynamicGrids ~/.julia/packages/DynamicGrids/IIwvg/src/framework.jl:187
 [13] simloop!
    @ ~/.julia/packages/DynamicGrids/IIwvg/src/framework.jl:171 [inlined]
 [14] #runsim!#144
    @ ~/.julia/packages/DynamicGrids/IIwvg/src/framework.jl:160 [inlined]
 [15] runsim!(output::ResultOutput{…}, simdata::DynamicGrids.SimData{…}, ruleset::Ruleset{…}, fspan::UnitRange{…})
    @ DynamicGrids ~/.julia/packages/DynamicGrids/IIwvg/src/framework.jl:156
 [16] sim!(output::ResultOutput{…}, ruleset::Ruleset{…}; init::Matrix{…}, mask::Nothing, tspan::UnitRange{…}, aux::Nothing, fps::Nothing, boundary::Remove{…}, proc::SingleCPU, opt::NoOpt, cellsize::Int64, timestep::Nothing, simdata::Nothing, kw::@Kwargs{})
    @ DynamicGrids ~/.julia/packages/DynamicGrids/IIwvg/src/framework.jl:97
 [17] sim!(output::ResultOutput{Matrix{…}, Vector{…}, Extent{…}}, ruleset::Ruleset{DynamicGrids.SimSettings{…}})
    @ DynamicGrids ~/.julia/packages/DynamicGrids/IIwvg/src/framework.jl:43
 [18] sim!(output::ResultOutput{…}, rules::SetNeighbors{…})
    @ DynamicGrids ~/.julia/packages/DynamicGrids/IIwvg/src/framework.jl:100
 [19] top-level scope
    @ ~/Documents/PhD/Chapitre_3/Scripts/Gene3sis/Julia-implementation/genesis_julia/dispersion/mwe_position_issue.jl:29
Some type information was truncated. Use `show(err)` to see complete types.

DynamicGrids is version 0.21.4

rafaqz commented 2 weeks ago

Yeah, use the dev branch for now...

remylegoff commented 2 weeks ago

I use the dev branch

rafaqz commented 2 weeks ago

Ok sorry I will check this later on

rafaqz commented 2 weeks ago

Oops this is just a discrepancy between Stencils.jl which uses indices and positions from the registered version of DynamicGrids.

Use indices for now. We should add a deprecated method const positions = indices in DynamicGrids.jl so code works accross versions.

using DynamicGrids, ColorSchemes, Colors, BenchmarkTools

const DEAD, ALIVE, BURNING = 1, 2, 3
init = fill(ALIVE, 400, 400)
bench_output = ResultOutput(init; tspan=1:200)
setneighbors_rule = let prob_combustion=0.0001, prob_regrowth=0.01
    SetNeighbors(Moore(1)) do data, neighborhood, cell, I
        if cell == DEAD
            if rand() <= prob_regrowth
                data[I...] = ALIVE
            end
        elseif cell == BURNING
            for pos in indices(neighborhood, I)
                if data[pos...] == ALIVE
                    data[pos...] = BURNING
                end
            end
            data[I...] = DEAD
        elseif cell == ALIVE
            if rand() <= prob_combustion 
                data[I...] = BURNING
            end
        end
    end
end
sim!(bench_output, setneighbors_rule)