Birth Death processes...?

[zahachtah]

Assume I simulate individual organisms. Sometimes they die and their mass goes to zero. In that case I don't want the effect of them to be calculated. I now do that by

function f(i,j,d2,d)
    d[i]+= u[i]>0.0 ? exp(-d2*dc)*u[j]*(u[j]/u[i])^ac : 0.0
    return d
end

in the function called by CellListMap.map_pairwise! which speeds it up somewhat. But can I also keep a list of organisms which are dead and make CellListMap.map_pairwise! just skip the call to the function altogether? When I have a birth of an organism, I take one dead organism and revive it with the new genetic makeup (traits) by setting its mass>0.

So the short question is, can I somehow keep a list of which organisms are alive without having to redo the entire SArray of positions for every death/birth event?

[zahachtah]

Just had a thought, could one do it as a view of a SArray?

xpositions = rand(SVector{2,Float64},10)

alive=[2,5,6,7]

system = PeriodicSystem(
           xpositions = view(xpositions,alive), 
           unitcell=[1.0,1.0], 
           cutoff = 0.1, 
           output = 0.0,
           output_name = :organism
       );

Now one could update the alive vector and this should make CellListMap.map_pairwise! skip pairwise calculations of any non-existing organisms....maybe.... (late now, will test tomorrow if you haven't commented yet :-) )

[lmiq]

That won't really bring you any benefit, because the coordinates are internally copied anyway to make them local to each computation.

The simplest approach is to simply put the live individuals at the beginning of the array and resize it. Partition the array in live/death is a relatively cheap operation (actually there is an internal function in CellListMap that can almost do that, here is a variation of it:

julia> function keep_alive!(alive, x::AbstractVector)
           iswap = 1
           @inbounds for i in eachindex(alive,x)
               if alive[i]
                   if iswap != i
                       x[iswap], x[i] = x[i], x[iswap]
                       alive[iswap], alive[i] = alive[i], alive[iswap]
                   end
                   iswap += 1
               end
           end
           return resize!(x,iswap - 1), resize!(alive,iswap - 1)
       end
keep_alive! (generic function with 1 method)

julia> x = rand(SVector{2,Float64}, 100); alive = rand(Bool, 100);

julia> keep_alive!(alive, x)
(SVector{2, Float64}[[0.29119093606094626, 0.14799284331418305], [0.4972730380281174, 0.09269265678798855], [0.9159907295858599, 0.023425334184397073], [0.6070998442833461, 0.5991916701778053], [0.3160812299339907, 0.3164821161774142], [0.7439114803685659, 0.11559108939285456], [0.6670382934847978, 0.9839740944475307], [0.657358664500978, 0.32044636378412406], [0.6082878186504177, 0.50068291214486], [0.08662133450538845, 0.4684000312786505]  …  [0.3151220417919539, 0.4605952237829716], [0.2875938461476413, 0.5834502287377114], [0.030323370160795338, 0.6514632971028673], [0.11688254221953154, 0.6127507999882679], [0.44709649959706943, 0.41760496956657467], [0.7726898156341978, 0.6743612022838364], [0.189673060868509, 0.2515943452830405], [0.6597804541954473, 0.028298360650823362], [0.002315181581732384, 0.30165618551226736], [0.9019060110665619, 0.9487908898003208]], Bool[1, 1, 1, 1, 1, 1, 1, 1, 1, 1  …  1, 1, 1, 1, 1, 1, 1, 1, 1, 1])

Then, for example, let us say that you initialized the system with the original coordinates:

julia> x = rand(SVector{2,Float64}, 100); alive = rand(Bool, 100);

julia> sys = PeriodicSystem(xpositions=x,unitcell=[1,1],cutoff=0.05,output=0.0)
PeriodicSystem1{output} of dimension 2, composed of:
    Box{OrthorhombicCell, 2}
      unit cell matrix = [ 1.0, 0.0; 0.0, 1.0 ]
      cutoff = 0.05
      number of computing cells on each dimension = [22, 22]
      computing cell sizes = [0.05, 0.05] (lcell: 1)
      Total number of cells = 484
    CellList{2, Float64}
      100 real particles.
      87 cells with real particles.
      123 particles in computing box, including images.
    Parallelization auxiliary data set for: 
      Number of batches for cell list construction: 1
      Number of batches for function mapping: 1
    Type of output variable (output): Float64

Note that the system has 100 real particles.

If you now "kill" the individuals:

julia> keep_alive!(alive, x);

The next time map_pairwise! is called the system will be updated for the new set of coordinates:

julia> PeriodicSystems.map_pairwise!((x,y,i,j,d2,out) -> out += d2, sys)
0.00416479370378711

julia> sys
PeriodicSystem1{output} of dimension 2, composed of:
    Box{OrthorhombicCell, 2}
      unit cell matrix = [ 1.0, 0.0; 0.0, 1.0 ]
      cutoff = 0.05
      number of computing cells on each dimension = [22, 22]
      computing cell sizes = [0.05, 0.05] (lcell: 1)
      Total number of cells = 484
    CellList{2, Float64}
      45 real particles.
      41 cells with real particles.
      55 particles in computing box, including images.
    Parallelization auxiliary data set for: 
      Number of batches for cell list construction: 1
      Number of batches for function mapping: 1
    Type of output variable (output): Float64

(note now that there are 45 real particles, the ones that were alive).

Note that that partition and resizing is cheap, and won't affect the overall performance of almost any simulation. Here with 10_000 particles:

julia> x = rand(SVector{2,Float64}, 10^4); alive = rand(Bool, 10^4);

julia> sys = PeriodicSystem(xpositions=x,unitcell=[1,1],cutoff=0.05,output=0.0);

julia> @btime PeriodicSystems.map_pairwise!((x,y,i,j,d2,out) -> out += d2, $sys)
  1.049 ms (144 allocations: 14.78 KiB)
489.7633383677577

julia> @btime keep_alive!(alive, x) setup=(x = rand(SVector{2,Float64}, 10^4); alive = rand(Bool, 10^4)) evals=1;
  46.064 μs (0 allocations: 0 bytes)

julia> keep_alive!(alive, x);

julia> @btime PeriodicSystems.map_pairwise!((x,y,i,j,d2,out) -> out += d2, $sys)
  499.010 μs (144 allocations: 14.78 KiB)
119.80104258632134

[zahachtah]

Great explanation,thank you very much. I think I can adapt this to birth events too. Enjoying the package, very efficient!