# Simplifying simulation by clusters

**URL:** <https://discourse.julialang.org/t/simplifying-simulation-by-clusters/35870>\
**Category:** Modelling & Simulations\
**Tags:** question\
**Created:** [March 12, 2020, 7:26am UTC](https://discourse.julialang.org/t/simplifying-simulation-by-clusters/35870 "2020-03-12T07:26:21Z")\
**Posts on this page:** 19\
**Page:** 1

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [March 12, 2020, 7:26am UTC](https://discourse.julialang.org/t/simplifying-simulation-by-clusters/35870/1 "2020-03-12T07:26:21Z")

</div>

This is partly a question asking for recommendations on an algorithm, ideally something that is implemented some Julia package, or can be done simply using one.

I am simulating a process in \mathbb{R}^3 (I can constrain them to a box, so think [0,1]^3 if that helps), with a number of agents. They _split_ following different paths, and thus their number can quickly explode. I only care about their distribution (their position captures all their state), so I would like to replace “close” points with their means and continue the calculation with a lower number of points.

To make thinks concrete, consider the MWE

```julia
using StaticArrays
l(x) = 1 / (1 + exp(-x)) # logistic, to keep samples in box
n = 1000
xs = SVector.(l.(randn(n)), l.(randn(n) .+ 0.1), l.(randn(n)) ./ 2)
ys = sort(vcat([x => 0.2 for x in xs], # weighted, pretend to split
               [x .* 0.9 .+ 0.5 * 0.1 => 0.8 for x in xs]),
          by = first)

```

I would like to find clusters of “close” points in `ys`, remove them, and replace them with their mean and added weight, so that I end up with, say, 100–500 points that I can cope with more easily.

“Close” can be any reasonable metric that makes the exercise simple.

My fallback solution is to repeatedly slice up the box along coordinates by using weighted medians (not unlike a decision tree). Eg

```julia
function median_split(by, x)
    # quick and dirty example, please don't focus on making this faster
    y = sort(x, by = by ∘ first)
    cum_weight = cumsum(last.(y))
    half_weight = cum_weight[end] / 2
    i = findfirst(x -> x > half_weight, cum_weight)
    y[1:i], y[(i+1):end]
end

median_split(x -> x[2], ys) # split on second coordinate

```

Any advice (including vague ideas, references to papers) is appreciated.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [March 12, 2020, 10:21am UTC](https://discourse.julialang.org/t/simplifying-simulation-by-clusters/35870/2 "2020-03-12T10:21:46Z")

</div>

In fluid simulation there are usually two approaches, particle based or grid based. There are also hybrids, where one switches from grid based inside a medium, to particle based near the interface to another medium. Maybe some inspiration to be drawn from there?  
Related video on the topic

[![](https://global.discourse-cdn.com/julialang/original/3X/8/d/8db1a58685d7f402bb62c85a407083e041b61f06.jpeg "Simulating Grains of Sand, Now 6 Times Faster") ](https://www.youtube.com/watch?v=e_9f5Z0sMYE)

and

[![](https://global.discourse-cdn.com/julialang/original/3X/5/e/5e22a9e4c1f1e815f7e784707a79b23281eddb71.jpeg "Narrow Band Liquid Simulations | Two Minute Papers #61") ](https://www.youtube.com/watch?v=nfPBT71xYVQ)

---

<div class="post-metadata">

**Author:** ![Sukera](https://avatars.discourse-cdn.com/v4/letter/s/ce7236/32.png) [@Sukera](https://discourse.julialang.org/u/Sukera)\
**Post date:** [March 12, 2020, 11:15am UTC](https://discourse.julialang.org/t/simplifying-simulation-by-clusters/35870/3 "2020-03-12T11:15:26Z")

</div>

Do you care about the total distribution or only some subset of the agents that satisfy some condition?

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [March 12, 2020, 12:00pm UTC](https://discourse.julialang.org/t/simplifying-simulation-by-clusters/35870/4 "2020-03-12T12:00:23Z")

</div>

> [@Sukera](#):
>
> Do you care about the total distribution or only some subset of the agents that satisfy some condition?

The whole problem works like this, conceptually.

Individuals move in a space (a, i)\_t \in \mathbb{R}^3 \times {1, \dots, N}. I have an initial discrete distribution over (a, i)\_0.

The law of motion for the a part is

a' = f(a, i)

As for the i, it is stochastic, with a distribution

i' \sim g(a, i, \theta)

where \theta are parameters to the problem I am estimating. And I care about linear functions of

h(a, i, t, \theta)

(eg a sum/integral over everything, but also split by i and/or t).

It is really important to _keep everything differentiable_ in \theta. Naturally it is, but there are about 6 million combinations by the time I am done simulating. I am estimating \theta, so I prefer sub-minute runtimes for the whole thing.

My initial idea was to pretend an equi-distributed g over i (or fix \theta), explore the space for a using this to see where I would bin, then approximate the f part using bins.

We have a running version in Julia (will be public when we submit the paper), it is just that the approximation part is very inaccurate (but differentiably inaccurate, which is a big win 😉).

> [@baggepinnen](#):
>
> In fluid simulation there are usually two approaches, particle based or grid based.

This could be yet another instance of engineers/physicists being light years ahead of economists so I want to read up on this. What would be a good crash course in the grid based methods? Is there a tutorial paper or a textbook? I can deal with the math and (classical) physics.

---

<div class="post-metadata">

**Author:** ![ericphanson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ericphanson/32/215186_2.png) [@ericphanson](https://discourse.julialang.org/u/ericphanson)\
**Post date:** [March 12, 2020, 12:04pm UTC](https://discourse.julialang.org/t/simplifying-simulation-by-clusters/35870/5 "2020-03-12T12:04:38Z")

</div>

> [@Tamas\_Papp](#):
>
> I would like to find clusters of “close” points in `ys` , remove them, and replace them with their mean and added weight, so that I end up with, say, 100–500 points that I can cope with more easily.

It’s not my area, but I’ve heard of this kind of idea by the name “sequential importance resampling”, so maybe those keywords could help find some literature on it.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [March 13, 2020, 5:53am UTC](https://discourse.julialang.org/t/simplifying-simulation-by-clusters/35870/6 "2020-03-13T05:53:57Z")

</div>

Unfortunately, my only exposure to this is the two-minute papers channel on YouTube so I can’t be much help regarding references. My first thought was actually to look at heat transfer, and discretization of PDEs in general.

Like Eric said though, it does have some similarities with sequential Monte-Carlo as well.

---

<div class="post-metadata">

**Author:** ![Vasily\_Pisarev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vasily_pisarev/32/7929_2.png) [@Vasily\_Pisarev](https://discourse.julialang.org/u/Vasily_Pisarev)\
**Post date:** [March 13, 2020, 10:59am UTC](https://discourse.julialang.org/t/simplifying-simulation-by-clusters/35870/7 "2020-03-13T10:59:24Z")

</div>

Will a [cell list](http://cacs.usc.edu/education/cs596/01-1LinkedListCell.pdf) approach be helpful for the task?

That’s a standard technique to facilitate neighbor search in molecular simulations, then the neighbors lists can be used to search clusters of points (using Stillinger’s criterion?)

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [March 13, 2020, 3:59pm UTC](https://discourse.julialang.org/t/simplifying-simulation-by-clusters/35870/8 "2020-03-13T15:59:44Z")

</div>

The objects I am keeping track of do not interact, so I am not sure.

---

<div class="post-metadata">

**Author:** ![Vasily\_Pisarev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vasily_pisarev/32/7929_2.png) [@Vasily\_Pisarev](https://discourse.julialang.org/u/Vasily_Pisarev)\
**Post date:** [March 13, 2020, 5:47pm UTC](https://discourse.julialang.org/t/simplifying-simulation-by-clusters/35870/9 "2020-03-13T17:47:32Z")

</div>

The interaction isn’t implied, the cell lists are only used to spped up finding all pairs (i, j) of particles in [0; L]^3 such that ||\textbf{r}\_i - \textbf{r}\_j|| \< r\_0.

And I remembered [this article](https://drops.dagstuhl.de/opus/volltexte/2019/10879/pdf/OASIcs-SLATE-2019-12.pdf) with an implementation, among others, of gravitational simulation with an octree space partitioning, which should be useful for finding clusters.

---

<div class="post-metadata">

**Author:** ![Balinus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/balinus/32/243_2.png) [@Balinus](https://discourse.julialang.org/u/Balinus)\
**Post date:** [March 13, 2020, 5:48pm UTC](https://discourse.julialang.org/t/simplifying-simulation-by-clusters/35870/10 "2020-03-13T17:48:29Z")

</div>

As mentioned by @ericphanson and @baggepinnen, the problem of distances of your agents does sounds like a Particle filter problem.

Maybe take a look at this paper, specifically section 3 about reducing the variance in the weight. i find it very easy to approach. It is targeted at geophysical applications but it might be relevent.

van Leeuwen, P.J., 2009: [Particle Filtering in Geophysical Systems.](https://journals.ametsoc.org/doi/abs/10.1175/2009MWR2835.1)_Mon. Wea. Rev.,_ **137** , 4089–4114, [https://doi.org/10.1175/2009MWR2835.1](https://doi.org/10.1175/2009MWR2835.1)

---

<div class="post-metadata">

**Author:** ![francesco.alemanno](https://avatars.discourse-cdn.com/v4/letter/f/e8c25b/32.png) [@francesco.alemanno](https://discourse.julialang.org/u/francesco.alemanno)\
**Post date:** [March 13, 2020, 7:23pm UTC](https://discourse.julialang.org/t/simplifying-simulation-by-clusters/35870/11 "2020-03-13T19:23:37Z")

</div>

would a markov model be of any use?  
you could think about discretizing your space, and modeling the probabilities of an agent moving from a particular cell into one of the neighbouring cells.  
i think this can work nicely since you have an equation that directly gives a state transition

(a',i')\sim T(a,i;\theta)

the memory cost for this approach would scale cubically w.r.t the lattice size, and the whole algorithm will surely be differentiable by \theta.

---

<div class="post-metadata">

**Author:** ![francesco.alemanno](https://avatars.discourse-cdn.com/v4/letter/f/e8c25b/32.png) [@francesco.alemanno](https://discourse.julialang.org/u/francesco.alemanno)\
**Post date:** [March 14, 2020, 9:13am UTC](https://discourse.julialang.org/t/simplifying-simulation-by-clusters/35870/12 "2020-03-14T09:13:27Z")

</div>

made a very toyish working example of what i mean:

[https://github.com/francescoalemanno/TestMarkov.jl/](https://github.com/francescoalemanno/TestMarkov.jl/)

core function, keep in mind i coded it this morning, and the code is not polished in any way

```julia
module TestMarkov 

#........

 function simulate(θ)
        X=0:0.1:10
        Y=0:0.1:10
        L=build_lattice(X,Y)
        normL=0.0
        iS=(1.0,2.5)
        S=iS
        for i in 1:100000
            P=approx_fpos((X,Y),S) # find approximate position of the state inside the lattice
            L[P]+=1 #accumulate probabilistic knowledge
            normL+=1
            S=T(S,θ) # EVOLVE STATE
            if i%1000 == 0
                # after some iterations interrupt this trajectory and start a new one according to the density in L
                # resampling from the knowledge of L
                S=resample_state(L,(X,Y))
                #or just reset state
                #S=iS
            end
        end
        X,Y,L
    end

#........
end

```

this are two runs of the simulation with different theta, particles which get too near are accumulated in the same cell, so the clustering is automatic, and the nice thing is that where the clustering is low, you can look at individual trajectories, without any need to keep track of each individual trajectory, indeed in the code at any point in time the only thing in memory regarding particles is one single state.

 ![paths](https://global.discourse-cdn.com/julialang/original/3X/1/7/17c4813187e389296a0a2a52c91f72ee3303de01.png)

---

<div class="post-metadata">

**Author:** ![Ahmed\_Salih](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmed_salih/32/206579_2.png) [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)\
**Post date:** [April 26, 2020, 1:00am UTC](https://discourse.julialang.org/t/simplifying-simulation-by-clusters/35870/13 "2020-04-26T01:00:54Z")

</div>

Could you try to explain this figure in your link?

 ![image](https://global.discourse-cdn.com/julialang/original/3X/c/e/ce8cc6019b8897f294b9464b078907960e8b9a8a.png)

I am struggling quite a bit to get it, but I understand it as having a “header” which tells me there is four cells, and for example for cell 1, the entry is “E” i.e. empty, since n o particles are stored in this cell. But for cell 0, it says “4”, and then lscl[4] = 2, which makes sense since cell 0 holds particle 2. But how do I then find particle 0 and 4 then too?

Kind regards

---

<div class="post-metadata">

**Author:** ![Vasily\_Pisarev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vasily_pisarev/32/7929_2.png) [@Vasily\_Pisarev](https://discourse.julialang.org/u/Vasily_Pisarev)\
**Post date:** [April 26, 2020, 3:02am UTC](https://discourse.julialang.org/t/simplifying-simulation-by-clusters/35870/14 "2020-04-26T03:02:15Z")

</div>

`head[0]` is 4, meaning that the particle number 4 is in cell 0.  
Then, we go to `lscl[4]` and get 2, then to `lscl[2]` and get 0, then to `lscl[0]` and get -1, meaning that there are no more particles in the cell.  
The chains for other cells are  
-1  
6 → 5 → -1  
and  
7 → 3 → 1 → -1

---

<div class="post-metadata">

**Author:** ![Ahmed\_Salih](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmed_salih/32/206579_2.png) [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)\
**Post date:** [April 26, 2020, 10:21am UTC](https://discourse.julialang.org/t/simplifying-simulation-by-clusters/35870/15 "2020-04-26T10:21:54Z")

</div>

So now comes my follow up question, “how do I know before hand that I need to go to lscl[2] = 0?”

For example taking the chain for cell 0, I see clearly that I can find particle 4 and 2, due to the header and the index lscl[4], but how do I find particle 0, at lscl[2]? How do I know that I have to move to lscl[2]?

This is primarily what bothers me, and it seems simple but for some reason I cannot grasp it - maybe I am confusing construction of the list with using it.

Kind regards

---

<div class="post-metadata">

**Author:** ![Vasily\_Pisarev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vasily_pisarev/32/7929_2.png) [@Vasily\_Pisarev](https://discourse.julialang.org/u/Vasily_Pisarev)\
**Post date:** [April 26, 2020, 11:02am UTC](https://discourse.julialang.org/t/simplifying-simulation-by-clusters/35870/16 "2020-04-26T11:02:18Z")

</div>

`head[0]` is 4, meaning that particle number 4 is in cell 0 _and_ the number of the next particle is in `lscl[4]`. From there, you proceed with the same logic: `lscl[4]` is 2, so particle 2 is in cell and the next particle number is in `lscl[2]` etc. until you hit `EMPTY`.

---

<div class="post-metadata">

**Author:** ![Ahmed\_Salih](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmed_salih/32/206579_2.png) [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)\
**Post date:** [April 26, 2020, 11:07am UTC](https://discourse.julialang.org/t/simplifying-simulation-by-clusters/35870/17 "2020-04-26T11:07:04Z")

</div>

Ah incredible! It finally makes sense now. I can see now why this is a very smart approach since I have to store very little in memory and if used properly it reduces complexity to O(N), without understanding too much what that actually means.

I will try to implement the list constructor for this example now, I might ask more if some issue pops up.

Kind regards

---

<div class="post-metadata">

**Author:** ![Ahmed\_Salih](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmed_salih/32/206579_2.png) [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)\
**Post date:** [April 26, 2020, 3:08pm UTC](https://discourse.julialang.org/t/simplifying-simulation-by-clusters/35870/18 "2020-04-26T15:08:31Z")

</div>

So I managed to get the cell list algorithm to work and here is a showcase:

![image](https://global.discourse-cdn.com/julialang/original/3X/7/d/7decb6fd940705faa626576548e0797c7c4601f3.png)

Where altering colors denote connection a specific cell. My question is now how to go about selecting the correct cells?

Say I want to look at cell zero, then I care about these cells (red circle):

 ![image](https://global.discourse-cdn.com/julialang/original/3X/9/e/9e7b34dca6898a21f40d595aa72ede7b7161f0a4.png)

I can’t figure out a robust way to ensure that no matter where I am in my “grid” that I will only interact with adjacent and diagonal elements - and in the cases where no elements exist (shown by red x), I just want to not include it.

I’ve been looking through the documnt and they mention something about constructing a vector cell index, but I can’t really see how to go about this.

Kind regards

---

<div class="post-metadata">

**Author:** ![Ahmed\_Salih](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmed_salih/32/206579_2.png) [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)\
**Post date:** [April 26, 2020, 3:10pm UTC](https://discourse.julialang.org/t/simplifying-simulation-by-clusters/35870/19 "2020-04-26T15:10:01Z")

</div>

And for anyone interested this is my implementation:

```julia
using Plots
using LoopVectorization
# Data from example
#data =
#[
#0.234 0.318
#0.734 0.811
#0.393 0.396
#0.649 0.604
#0.359 0.147
# 0.68 0.258
#0.814 0.106
#0.882 0.691
#]

N = 2^13

pg = fill(tuple(0.0,0.0,0.0),N)

for i = 1:N
    pg[i] = map(rand, (Float64,Float64,Float64))
end
#pg = map(x-> x .* -1,pg)

# CONSTRUCTION OF CELLS
rc = [0.25;0.25]
# The total length of the simulation domain
La = [1;1]

function CellList2D(pg,rc,La)
    # Calculate the amount of cells in each direction
    Lca = [div(La[1],rc[1]); div(La[2],rc[2])]
    # Convert it to Int64
    Lca = convert(Array{Int64,1},Lca)
    # Number of cells in total
    nCell = Lca[1] * Lca[2]

    np = length(pg)
    # Preallocate header
    head = zeros(Int64,nCell)
    lscl = zeros(Int64,np)
    for i = 1:np
        # Calculate serial index
        c = div(pg[i][1],rc[1])*Lca[2] + div(pg[i][3],rc[2])
        # Convert to Julia indexing from C
        c = convert(Int64,c) + 1
        # Input serial index to header
        lscl[i] = head[c];
        # In header insert the last value, to keep track correctly
        head[c] = i;
    end
    #println(head)
    #println(lscl)
    return head,lscl
end

head,lscl = @time CellList2D(pg,rc,La)

scatter()
a = zeros(Int64,0)
for i in head
    append!(a,i)
    while i != 0
        i = lscl[i]
        append!(a,i)
    end
    pop!(a)
    Plots.display(scatter!(first.(pg[a]),last.(pg[a]),aspect_ratio=0.75,xlims=(0,1),ylims=(0,1),xticks=[0,0.25,0.5,0.75,1], yticks=[0,0.25,0.5,0.75,1],show=true,legend=:outertopright))
    #println(a)
    global a = zeros(Int64,0)
end

```

With performance of:

```julia
@benchmark CellList2D($pg,$rc,$La)
BenchmarkTools.Trial: 
  memory estimate: 64.50 KiB
  allocs estimate: 6
  --------------
  minimum time: 236.100 μs (0.00% GC)
  median time: 246.201 μs (0.00% GC)
  mean time: 254.215 μs (1.24% GC)
  maximum time: 5.767 ms (94.38% GC)
  --------------
  samples: 10000
  evals/sample: 1

```

For N = 2^13