-
-
Notifications
You must be signed in to change notification settings - Fork 16
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Thoughts #5
Comments
As soon as this works I'll try it on a couple of classical multi scale models in eco-evo field, e.g. evolutionary food webs, spatial evolutionary food webs, spatial vegetation models. |
Hey, Sorry this turned into a novel. Now that I know it works (tried it in DifferentialEquations.jl, yay it's working on my The Structure - why it's fast, and the constraintsLet me first explain the constraint that this has. If you take a look at an arbitrary part of a loop in DifferentialEquations.jl, you will see it's littered with these: for i in uidx
tmp[i] = u[i]+Δt*(a71*k1[i]+a72*k2[i]+a73*k3[i]+a74*k4[i]+a75*k5[i]+a76*k6[i])
end So anything that is
Will work, and will actually work without allocating (i.e. it will be a very fast loop). Here are the two main caveats:
So that's the main constraint: the ODE/SDE/etc. solvers are only doing the continuous updates on the terms top.x.x.x.x.x.... everything has to be in the Now that's the only real speed constraint. By design, since Lastly, talking about performance, if you have your type as immutable, then carnaval's PR (Github is down for me while writing this, mention to be this and I'll replace this with the link to the PR) will turn this into essentially a 0 overhead loop. So this is designed with all of the enhancements of v0.6 in mind so that way it boils down to something with no overhead in the most constrained case (using immutable types, always strict typing). The beauty of Julia is that, at any spot, we can "break" these stricter rules at the cost of some performance for some generality (and the further discussion will give some examples). How to add generality, and the costsSo there's speed in every place you look, with the generality of the hierarchy, but with the
Yes, you can use any
You can. immutable Organism{T} <: MultiScaleModelLeaf
x::Vector{T}
predator::Bool
end Now your leaf type will automatically update both predators and prey in the ODE, and you can distinguish between them inside of So those are two tools to use this for more generality.
You can already do this with immutable Population{T<:AbstractMultiScaleModel} <: AbstractMultiScaleModel
x::Vector{T}
y::Vector{BottomType}
end_idxs::Vector{Int}
end where
Would size need to be a type, or just a value in SpaceAnd for the last one, there's a lot to unfold.
Depends on the model. Remember that, other than this
The first way has one caveat: the changes are only applied after a step. So for example if you're taking really large steps, if you use the first way you may have particles/organisms/etc. to move every step because they should actually "be moving in-between the steps". You can constrain the stepsize or use a fixed-timestep ODE method (you can change any method to fixed timesteps with Events are nice because they will handle the discontinuity well. For example, your event could be "every second update the positions". So you could do what you wanted to do with 1, but instead check if the time has gone another full minute, and put an event there if it has. Side-note: ODE Callback/Event ExtensionUsing the notation of the events, here's how I am extending events to handle this. Instead of having default_callback = @ode_callback begin
@ode_savevalues
end cache=MyCache where default_callback = @ode_callback begin
@ode_savevalues
end just defaults to having the cache be some empty immutable cache (0 overhead). Then extend the event interface to be: function event_f(t,u,f,user_cache)
u[1]
end and function apply_event!(u,f,user_cache,solver_cache)
u[2] = -u[2]
end ( This will take less than an hour to do, and it's a breaking change so I'm doing this Back to the MultiScaleModelsNow here's how you can use this in a multi-scale model to do discrete timestep updates function event_f(t,u,f,user_cache)
60*user_cache.last_minute-t
end Notice that this is positive until So then the [Note that this is general What about continuous updates to space? Actually, diffusion of discrete particles However, if events happen too often, they constrain the solver's timesteps to be
Beautiful idea that I hadn't thought of, but the If there was more than one background protein, you can make f = function(t,u,du)
### Diffuse
for i in eachindex(u.y)
du.y[i] = A[i]*u.y[i]
end
### Now do stuff on `u.x`
end where now I have Defining the ODE, the internal looping structure
Yes, I think this is a good lead-in to describing the looping structure for At the highest level and lowest level this is easy. To loop over everything in the lowest for i in eachindex(u)
u[i] # Will get each thing that at a leaf, or in a `y` (when that's implemented)
end At the highest level, you can grab each child by for i in eachindex(u.x)
u.x[i] # So if the top is population, this is a cell
end I want to extend the indexing so that way, if you have Tissue -> Population -> Cells -> Protein concentrations (which are the numbers at the bottom, so the leaf type is Cell), then I want to add: for i in eachindex(u,2)
u[i] # u is a tissue, so level 2 is cells, so this loops through each cells
# i = (i1,i2), so i[1] would tell you which population, i[2] is which cell
# together, i just ends up giving the cell for convenience.
# Using this, you can change `u.y[i]` from the cells, so have the cells
# influence the PDE!
end This is actually a simple fix from what I have, and would handle these constraints. for i in eachindex(u)
u[i] # Will get each thing that at a leaf, or in a `y` (when that's implemented)
end so this generality is simply added as convenience for users (i.e. you can already Overarching Philosophy: Mixing Discrete and the Continuous BlurHere's the main idea: this let's you add structures, but also "blur" whenever necessary. You want to pick a middle, say cells, and then have the protein concentrations of the That's what this framework should One final point: Linear indexabilityOne point to note is that |
Update on the Is this still type inferrable? I wonder if the king of Julia arrays, @timholy, would Here's the indexing code. It's actually really succinct (other than my 4 AM bisection where I gave up and made it a vectorized calll... I'll fix that soon hahaha). https://github.com/JuliaDiffEq/MultiScaleModels.jl/blob/master/src/indexing.jl#L42 |
Great read and very promising! Looking forward to testing |
So yeah, I'm closing this to get more specific in issues. It works in ODEs and SDEs. The changing size is possible, but currently you just have to know which cache variables to change. I'll make the cache iterators soon, it's just that making them will take quite a bit more typing than just using them directly so behind the scenes I've been avoiding this! The big change from this idea, adding There's still some performance improvements that could occur (especially in SDEs), but the full functionality should be available very soon (it really needed the new integrator + callback interfaces!) |
Hi Chris
Gitter is down atm (and github very slow) so I use an issue to write down stuff, hope thats ok.
The idea of multiscalemodels is brilliant. Your approach made it so much more general too.
I assume one can use e.g. space as Head node? for linear space its analog:
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)
can you mix leaf types? for example have predator and prey types that dispatch the growth dynamics function based on type?
parameters. I used parametrised functions to give parameters (e.g. traits) to the ode function. Will these AbstractMultiScaleModel types now hold these parameters?
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?
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
The text was updated successfully, but these errors were encountered: