diff --git a/model/src/main/scala/fr/iscpif/spacematters/model/container/Container.scala b/model/src/main/scala/fr/iscpif/spacematters/model/container/Container.scala
new file mode 100644
index 0000000..f2334cb
--- /dev/null
+++ b/model/src/main/scala/fr/iscpif/spacematters/model/container/Container.scala
@@ -0,0 +1,51 @@
+/*
+ * Copyright (C) 2015 Romain Reuillon
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Affero General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+package fr.iscpif.spacematters.model.container
+
+import fr.iscpif.spacematters.model.Cell
+
+import scala.util.Random
+
+trait Container {
+
+  def size : Int
+
+  def container(implicit rng: Random): Seq[Seq[Cell]]
+
+
+
+  /**
+   * Stringify just to check validity of generators ; has no sense in general context as a new instance will be randomly generated at each
+   * call of container.
+   *
+   * @return
+   */
+  override def toString: String = {
+    var res = ""
+    container(new Random).foreach(
+      (row : Seq[Cell]) => {
+          row.foreach(
+            (c : Cell) => {res = res +c.capacity+" | "}
+          )
+          res = res + "\n"
+       }
+    )
+    res
+  }
+
+
+}
diff --git a/model/src/main/scala/fr/iscpif/spacematters/model/container/ExpMixtureContainer.scala b/model/src/main/scala/fr/iscpif/spacematters/model/container/ExpMixtureContainer.scala
new file mode 100644
index 0000000..0bb2358
--- /dev/null
+++ b/model/src/main/scala/fr/iscpif/spacematters/model/container/ExpMixtureContainer.scala
@@ -0,0 +1,38 @@
+package fr.iscpif.spacematters.model.container
+
+import fr.iscpif.spacematters.model.Cell
+
+import scala.util.Random
+
+trait ExpMixtureContainer <: Container {
+
+    /** maximal capacity C_m */
+    def maxCapacity : Int
+
+    /** Size of exponential kernels, of the form C_m*exp(-||x-x_0||/r_0) */
+    def kernelRadius : Double
+
+    /** Number of exponential kernels */
+    def centersNumber : Int
+
+    def container(implicit rng: Random) = {
+      val arrayVals = Array.fill[Cell](size, size) {
+        new Cell(0, 0, 0)
+      }
+
+      // generate random center positions
+      val centers = Array.fill[Int](centersNumber, 2) {
+        rng.nextInt(size)
+      }
+
+      for (i <- 0 to size - 1; j <- 0 to size - 1) {
+        for (c <- 0 to centersNumber - 1) {
+          arrayVals(i)(j).capacity = arrayVals(i)(j).capacity + maxCapacity * math.exp(-math.sqrt(math.pow((i - centers(c)(0)), 2) + math.pow((j - centers(c)(1)), 2)) / kernelRadius)
+        }
+      }
+
+      Seq.tabulate(size,size){(i:Int,j:Int)=>Cell(arrayVals(i)(j).capacity,0,0) }
+
+    }
+
+ }
diff --git a/model/src/main/scala/fr/iscpif/spacematters/model/container/PrefAttDiffusionContainer.scala b/model/src/main/scala/fr/iscpif/spacematters/model/container/PrefAttDiffusionContainer.scala
new file mode 100644
index 0000000..b53cf0c
--- /dev/null
+++ b/model/src/main/scala/fr/iscpif/spacematters/model/container/PrefAttDiffusionContainer.scala
@@ -0,0 +1,109 @@
+package fr.iscpif.spacematters.model.container
+
+import fr.iscpif.spacematters.model.Cell
+
+import scala.util.Random
+
+trait PrefAttDiffusionContainer <: Container {
+
+    /** sum of all capacities */
+    def totalCapacity : Double
+
+    /** diffusion parameters */
+    def diffusion : Double
+    def diffusionSteps : Int
+
+    /** Growth rate */
+    def growthRate : Int
+
+    /** Preferential attachment parameter */
+    def alphaAtt : Double
+
+
+
+
+  /**
+   *
+   * @param rng
+   * @return
+   */
+    def container(implicit rng: Random) = {
+      var arrayVals = Array.fill[Cell](size, size){new Cell(0,0,0)}
+      var population :Double = 0
+
+      while(population < totalCapacity) {
+
+        // add new population following pref att rule
+        if(population==0){
+          //choose random patch
+          for(_<- 1 to growthRate){val i = rng.nextInt(size);val j = rng.nextInt(size) ; arrayVals(i)(j).capacity =  arrayVals(i)(j).capacity + 1 }
+        } else{
+          val oldPop = arrayVals.clone()
+          val ptot = oldPop.flatten.map((c:Cell)=>math.pow(c.capacity / population , alphaAtt)).sum
+
+          /*
+          val drawings = Array.fill[Double](growthRate){rng.nextDouble()}
+          Sorting.quickSort(drawings)
+          var s = 0.0; var i =0;var j =0;var k=0;
+          oldPop.flatten.foreach(
+            (c:Cell)=>{
+              s=s+(math.pow(c.capacity/population,alphaAtt)/ptot)
+              while(s<drawings(k)){arrayVals(i)(j).capacity = arrayVals(i)(j).capacity + 1;k=k+1;}
+              j=j+1;if(j==size){j=0;i=i+1};
+            }
+          )
+          */
+
+          for(_<- 1 to growthRate){
+            var s = 0.0 ; val r = rng.nextDouble() ; var i =0;var j =0
+            while(s<r){s=s+(math.pow(oldPop(i)(j).capacity/population,alphaAtt)/ptot);j=j+1;if(j==size){j=0;i=i+1}}
+            if(j==0){j=size-1;i = i - 1}else{j=j-1};
+            arrayVals(i)(j).capacity = arrayVals(i)(j).capacity + 1
+          }
+        }
+
+        // diffuse
+        for (_ <- 1 to diffusionSteps) {
+          arrayVals = diffuse(arrayVals, diffusion)
+        }
+
+        // update total population
+        population = arrayVals.flatten.map(_.capacity).sum
+
+        //println("population : "+population)
+      }
+
+      Seq.tabulate(size,size){(i:Int,j:Int)=>Cell(arrayVals(i)(j).capacity,0,0) }
+
+
+    }
+
+
+
+
+
+
+  /**
+   * Diffuse to neighbors proportion alpha of capacities
+   * @param a
+   */
+  def diffuse(a : Array[Array[Cell]],alpha : Double): Array[Array[Cell]] = {
+    val newVals = a.clone()
+    for (i <- 0 to size - 1; j <- 0 to size - 1) {
+      if(i>=1){newVals(i-1)(j).capacity = newVals(i-1)(j).capacity + (alpha / 8)*a(i)(j).capacity ; newVals(i)(j).capacity = newVals(i)(j).capacity  - (alpha / 8)*a(i)(j).capacity }
+      if(i<size-1){newVals(i+1)(j).capacity = newVals(i+1)(j).capacity + (alpha / 8)*a(i)(j).capacity ; newVals(i)(j).capacity = newVals(i)(j).capacity  - (alpha / 8)*a(i)(j).capacity }
+      if(j>=1){newVals(i)(j-1).capacity = newVals(i)(j-1).capacity + (alpha / 8)*a(i)(j).capacity ; newVals(i)(j).capacity = newVals(i)(j).capacity  - (alpha / 8)*a(i)(j).capacity }
+      if(j<size-1){newVals(i)(j+1).capacity = newVals(i)(j+1).capacity + (alpha / 8)*a(i)(j).capacity ; newVals(i)(j).capacity = newVals(i)(j).capacity  - (alpha / 8)*a(i)(j).capacity }
+      if(i>=1&&j>=1){newVals(i-1)(j-1).capacity = newVals(i-1)(j-1).capacity + (alpha / 8)*a(i)(j).capacity ; newVals(i)(j).capacity = newVals(i)(j).capacity  - (alpha / 8)*a(i)(j).capacity }
+      if(i>=1&&j<size-1){newVals(i-1)(j+1).capacity = newVals(i-1)(j+1).capacity + (alpha / 8)*a(i)(j).capacity ; newVals(i)(j).capacity = newVals(i)(j).capacity  - (alpha / 8)*a(i)(j).capacity }
+      if(i<size-1&&j>=1){newVals(i+1)(j-1).capacity = newVals(i+1)(j-1).capacity + (alpha / 8)*a(i)(j).capacity ; newVals(i)(j).capacity = newVals(i)(j).capacity  - (alpha / 8)*a(i)(j).capacity }
+      if(i<size-1&&j<size-1){newVals(i+1)(j+1).capacity = newVals(i+1)(j+1).capacity + (alpha / 8)*a(i)(j).capacity ; newVals(i)(j).capacity = newVals(i)(j).capacity  - (alpha / 8)*a(i)(j).capacity }
+    }
+    newVals
+  }
+
+
+
+
+
+ }
diff --git a/model/src/main/scala/fr/iscpif/spacematters/model/container/RandomContainer.scala b/model/src/main/scala/fr/iscpif/spacematters/model/container/RandomContainer.scala
new file mode 100644
index 0000000..d6a5e6e
--- /dev/null
+++ b/model/src/main/scala/fr/iscpif/spacematters/model/container/RandomContainer.scala
@@ -0,0 +1,19 @@
+package fr.iscpif.spacematters.model.container
+
+import fr.iscpif.spacematters.model.Cell
+
+import scala.util.Random
+
+
+trait RandomContainer <: Container {
+
+   def maxCapacity : Int
+
+   def container(implicit rng: Random) = {
+     Seq.fill(size, size) {
+       val capacity: Double = rng.nextInt(maxCapacity).toDouble
+       Cell(capa = capacity, green = 0, red = 0)
+     }
+   }
+
+}
diff --git a/model/src/main/scala/fr/iscpif/spacematters/model/test/TestContainer.scala b/model/src/main/scala/fr/iscpif/spacematters/model/test/TestContainer.scala
new file mode 100644
index 0000000..789e821
--- /dev/null
+++ b/model/src/main/scala/fr/iscpif/spacematters/model/test/TestContainer.scala
@@ -0,0 +1,36 @@
+package fr.iscpif.spacematters.model.test
+
+import fr.iscpif.spacematters.model.container._
+
+import scala.util.Random
+
+object TestContainer extends App {
+
+  val t = System.currentTimeMillis()
+
+  implicit val rng = new Random
+
+  val cont = new  PrefAttDiffusionContainer{
+    override def size: Int = 50
+
+    /*// exp Mixture params
+    override def kernelRadius = 0.5
+    override def centersNumber = 2
+    */
+
+    //prefAttdiff
+    override def totalCapacity :Double = 10000
+    override def diffusion : Double = 0.02
+    override def diffusionSteps : Int = 2
+    override def growthRate : Int = 100
+    override def alphaAtt : Double = 1.1
+  }
+
+  // generate instance of container
+  println(cont)
+  // an other
+  println(cont)
+
+  //println((System.currentTimeMillis()-t)/1000.0)
+
+}
diff --git a/scala/model/src/main/scala/fr/iscpif/spacematters/model/Cell.scala b/scala/model/src/main/scala/fr/iscpif/spacematters/model/Cell.scala
index e0a035b..0bd701a 100644
--- a/scala/model/src/main/scala/fr/iscpif/spacematters/model/Cell.scala
+++ b/scala/model/src/main/scala/fr/iscpif/spacematters/model/Cell.scala
@@ -20,7 +20,10 @@ object Cell {
   val empty = Cell(0, 0, 0)
 }
 
-case class Cell(capacity: Int, green: Int, red: Int) {
+case class Cell(capa: Double, green: Int, red: Int) {
+
+  var capacity = capa
+
   def isFull = population >= capacity
   def isEmpty = population <= 0
   def population = red + green
diff --git a/scala/model/src/main/scala/fr/iscpif/spacematters/model/Schelling.scala b/scala/model/src/main/scala/fr/iscpif/spacematters/model/Schelling.scala
index 4f5320a..d821370 100644
--- a/scala/model/src/main/scala/fr/iscpif/spacematters/model/Schelling.scala
+++ b/scala/model/src/main/scala/fr/iscpif/spacematters/model/Schelling.scala
@@ -18,8 +18,6 @@ package fr.iscpif.spacematters.model
 
 import fr.iscpif.spacematters.model.initial._
 import fr.iscpif.spacematters.model.move._
-import fr.iscpif.spacematters.model.initial.InitialState
-import fr.iscpif.spacematters.model.move.Moves
 
 import scala.util.Random
 
diff --git a/scala/model/src/main/scala/fr/iscpif/spacematters/model/Simulation.scala b/scala/model/src/main/scala/fr/iscpif/spacematters/model/Simulation.scala
index b7a41cc..e822aad 100644
--- a/scala/model/src/main/scala/fr/iscpif/spacematters/model/Simulation.scala
+++ b/scala/model/src/main/scala/fr/iscpif/spacematters/model/Simulation.scala
@@ -6,7 +6,8 @@ import initial._
 import move._
 import metric._
 import stop._
-import Moran._
+import container._
+import metric.Moran
 
 import scalax.io.Resource
 
@@ -16,12 +17,21 @@ object Simulation extends App {
 
   implicit val rng = new Random
 
-  val simulation = new Schelling with RandomState with RandomMoves with SpeilmanStop {
+  val simulation = new Schelling with RandomState with RandomContainer with RandomMoves with SpeilmanStop {
     override def size: Int = 50
     override def greenRatio: Double = 0.5
     override def redRatio: Double = 0.35
     override def maxCapacity: Int = 50
     override def similarWanted: Double = 0.4
+
+    /*
+    override def totalCapacity :Double = 50000
+    override def diffusion : Double = 0.02
+    override def diffusionSteps : Int = 2
+    override def growthRate : Int = 100
+    override def alphaAtt : Double = 1.1
+*/
+
   }
 
   simulation.run
@@ -40,7 +50,7 @@ object Simulation extends App {
     (state, step) <- simulation.states.take(100).zipWithIndex
   } {
     def unsatisfied = simulation.unsatisfieds(state).map(_.number).sum
-    println(s"Step $step: # of unsatisfied: $unsatisfied, Dissimilarity D: ${"%.3f".format(dissimilarity(state, Green, Red))}, Moran I Red: ${"%.3f".format(colorRatioMoran(state, Red))}, Entropy H: ${"%.3f".format(segregationEntropy(state, Green, Red))}, Exposure Reds to Greens :${"%.3f".format(exposureOfColor1ToColor2(state, Red, Green))}, Isolation Reds :${"%.3f".format(isolation(state, Red, Green))}, Concentration Greens : ${"%.3f".format(delta(state, Green, Red))}")
+    println(s"Step $step: # of unsatisfied: $unsatisfied, Dissimilarity D: ${"%.3f".format(dissimilarity(state, Green, Red))}, Moran I Red: ${"%.3f".format(Moran.colorRatioMoran(state, Red))}, Entropy H: ${"%.3f".format(segregationEntropy(state, Green, Red))}, Exposure Reds to Greens :${"%.3f".format(exposureOfColor1ToColor2(state, Red, Green))}, Isolation Reds :${"%.3f".format(isolation(state, Red, Green))}, Concentration Greens : ${"%.3f".format(delta(state, Green, Red))}")
 
     for { (position @ (i, j), c) <- state.cells } {
       def agents = Color.all.map(_.cellColor.get(c)).mkString(",")
@@ -56,7 +66,7 @@ object Simulation extends App {
     val similarWanted = simulation.similarWanted
 
     output2.append(
-      s"""$step, $unsatisfied,${dissimilarity(state, Green, Red)}, ${colorRatioMoran(state, Red)}, ${segregationEntropy(state, Green, Red)}, ${exposureOfColor1ToColor2(state, Red, Green)},${exposureOfColor1ToColor2(state, Green, Red)}, ${isolation(state, Red, Green)}, ${isolation(state, Green, Red)},${delta(state, Red, Green)},${delta(state, Green, Red)}, $size, $greenRatio,$redRatio, $maxCapacity, $similarWanted\n""".stripMargin)
+      s"""$step, $unsatisfied,${dissimilarity(state, Green, Red)}, ${Moran.colorRatioMoran(state, Red)}, ${segregationEntropy(state, Green, Red)}, ${exposureOfColor1ToColor2(state, Red, Green)},${exposureOfColor1ToColor2(state, Green, Red)}, ${isolation(state, Red, Green)}, ${isolation(state, Green, Red)},${delta(state, Red, Green)},${delta(state, Green, Red)}, $size, $greenRatio,$redRatio, $maxCapacity, $similarWanted\n""".stripMargin)
 
   }
 
diff --git a/scala/model/src/main/scala/fr/iscpif/spacematters/model/initial/RandomState.scala b/scala/model/src/main/scala/fr/iscpif/spacematters/model/initial/RandomState.scala
index 08d2722..8e70426 100644
--- a/scala/model/src/main/scala/fr/iscpif/spacematters/model/initial/RandomState.scala
+++ b/scala/model/src/main/scala/fr/iscpif/spacematters/model/initial/RandomState.scala
@@ -17,18 +17,17 @@
 package fr.iscpif.spacematters.model.initial
 
 import fr.iscpif.spacematters.model._
+import fr.iscpif.spacematters.model.container._
 
 import scala.util.Random
 
-trait RandomState <: InitialState { self: Schelling =>
+trait RandomState <: InitialState with Container { self: Schelling =>
 
   def maxCapacity: Int
 
   def initialState(implicit rng: Random) = {
-    val cells = Seq.fill(size, size) {
-      val capacity = rng.nextInt(maxCapacity)
-      Cell(capacity = capacity, green = 0, red = 0)
-    }
+
+    val cells = container(rng)
 
     val totalCapacity = cells.flatten.map(_.capacity).sum
     val greens = (totalCapacity * greenRatio).toInt
diff --git a/scala/model/src/main/scala/fr/iscpif/spacematters/model/metric/Convolution.scala b/scala/model/src/main/scala/fr/iscpif/spacematters/model/metric/Convolution.scala
new file mode 100644
index 0000000..db6f6d3
--- /dev/null
+++ b/scala/model/src/main/scala/fr/iscpif/spacematters/model/metric/Convolution.scala
@@ -0,0 +1,59 @@
+package fr.iscpif.spacematters.model.metric
+
+import org.apache.commons.math3.complex.Complex
+import org.apache.commons.math3.transform.{ TransformType, DftNormalization, FastFourierTransformer }
+import org.apache.commons.math3.util.MathArrays
+
+import scala.math._
+
+object Convolution {
+
+  /**
+   * Generic convol for double Arrays (in O(nlog(n)), using FFT)
+   *
+   * @param x
+   * @param k centered kernel
+   * @return y = x*k with y_i = \sum_{j=1}{|K|}{x_{i-j-|K|/2}*k_j}
+   */
+  def convolution(x: Array[Double], k: Array[Double]): Array[Double] = {
+    val xl = pow(2.0, ceil(log(x.length) / log(2.0)) + 1)
+    val xp = x.padTo(x.length + (xl.toInt - x.length) / 2, 0.0).reverse.padTo(xl.toInt, 0.0).reverse
+    val kp = k.padTo(k.length + (xl.toInt - k.length) / 2, 0.0).reverse.padTo(xl.toInt, 0.0).reverse
+    val tr = new FastFourierTransformer(DftNormalization.STANDARD)
+    val ftx = tr.transform(xp, TransformType.FORWARD)
+    val ftk = tr.transform(kp, TransformType.FORWARD)
+    val real = MathArrays.ebeSubtract(MathArrays.ebeMultiply(ftx.map { z => z.getReal }, ftk.map { z => z.getReal }), MathArrays.ebeMultiply(ftx.map { z => z.getImaginary }, ftk.map { z => z.getImaginary }))
+    val im = MathArrays.ebeAdd(MathArrays.ebeMultiply(ftx.map { z => z.getReal }, ftk.map { z => z.getImaginary }), MathArrays.ebeMultiply(ftx.map { z => z.getImaginary }, ftk.map { z => z.getReal }))
+    val trinv = tr.transform(Array.tabulate(real.length) { i => new Complex(real(i), im(i)) }, TransformType.INVERSE).map { z => z.getReal }
+    trinv.splitAt(trinv.length - x.length / 2)._2 ++ trinv.splitAt(x.length - x.length / 2)._1
+  }
+
+  /**
+   * Square convol (for tests)
+   *
+   * @param x
+   * @param k
+   * @return
+   */
+  def directConvol(x: Array[Double], k: Array[Double]): Array[Double] = {
+    val kl = k.length
+    val xpadded = x.padTo(x.length + kl, 0.0).reverse.padTo(x.length + 2 * kl, 0.0).reverse
+    Array.tabulate(x.length + k.length) { i => MathArrays.ebeMultiply(k.reverse, xpadded.splitAt(i + 1)._2.splitAt(k.length)._1).sum }
+  }
+
+  /**
+   *  2D convolution
+   *  Using bijection [|1,N|]2 ~ [|1,N|] by flattening, after having good paddling
+   *
+   * @param x
+   * @param k
+   */
+  def convolution2D(x: Array[Array[Double]], k: Array[Array[Double]]): Array[Array[Double]] = {
+    val xpad = x.map { row => row.padTo(2 * (row.length / 2) + 1 + k(0).length, 0.0).reverse.padTo(2 * (row.length / 2) + 1 + 2 * k(0).length, 0.0).reverse }.padTo(2 * (x.length / 2) + 1 + k.length, Array.fill(2 * (x(0).length / 2) + 1 + 2 * k(0).length) { 0.0 }).reverse.padTo(2 * (x.length / 2) + 1 + 2 * k.length, Array.fill(2 * (x(0).length / 2) + 1 + 2 * k(0).length) { 0.0 })
+    val xpos = Array.fill(x.length, x(0).length) { 1.0 }.map { row => row.padTo(2 * (row.length / 2) + 1 + k(0).length, 0.0).reverse.padTo(2 * (row.length / 2) + 1 + 2 * k(0).length, 0.0).reverse }.padTo(2 * (x.length / 2) + 1 + k.length, Array.fill(2 * (x(0).length / 2) + 1 + 2 * k(0).length) { 0.0 }).reverse.padTo(2 * (x.length / 2) + 1 + 2 * k.length, Array.fill(2 * (x(0).length / 2) + 1 + 2 * k(0).length) { 0.0 }).flatten
+    val kpad = k.map { row => row.padTo(row.length + (xpad(0).length - row.length) / 2, 0.0).reverse.padTo(xpad(0).length, 0.0).reverse }.padTo(k.length + (xpad.length - k.length) / 2, Array.fill(xpad(0).length) { 0.0 }).reverse.padTo(xpad.length, Array.fill(xpad(0).length) { 0.0 })
+    val flatconv = convolution(xpad.flatten, kpad.flatten)
+    flatconv.zipWithIndex.filter { case (_, j) => xpos(j) > 0 }.map { case (d, _) => d }.sliding(x(0).length, x.length).toArray.reverse
+  }
+
+}
diff --git a/scala/model/src/main/scala/fr/iscpif/spacematters/model/metric/package.scala b/scala/model/src/main/scala/fr/iscpif/spacematters/model/metric/package.scala
index 9ad1ff5..a03da37 100644
--- a/scala/model/src/main/scala/fr/iscpif/spacematters/model/metric/package.scala
+++ b/scala/model/src/main/scala/fr/iscpif/spacematters/model/metric/package.scala
@@ -190,4 +190,50 @@ package object metric {
       if i != oi || j != oj
     } yield state(i + oi)(j + oj)
 
+  /**
+   * Moran index using fast convolution.
+   *
+   * @param matrix
+   * @return
+   */
+  def moran_convol(matrix: Matrix[T], quantity: Quantity[T]): Double = {
+    val conf = matrix.map { row => row.map(quantity).toArray }.toArray
+    val n = conf.length
+    val flatConf = conf.flatten
+    val popMean = flatConf.sum / flatConf.length
+    val centeredConf = conf.map { r => r.map { d => d - popMean } }
+    val variance = MathArrays.ebeMultiply(centeredConf.flatten, centeredConf.flatten).sum
+    val weights = spatialWeights(2 * n - 1)
+    val totWeight = Convolution.convolution2D(Array.fill(n, n) { 1.0 }, weights).flatten.sum
+    flatConf.length / (totWeight * variance) * MathArrays.ebeMultiply(centeredConf.flatten, Convolution.convolution2D(centeredConf, weights).flatten).sum
+  }
+
+  def spatialWeights(n: Int): Array[Array[Double]] = {
+    Array.tabulate(n, n) { (i, j) => if (i == n / 2 && j == n / 2) 0.0 else 1 / Math.sqrt((i - n / 2) * (i - n / 2) + (j - n / 2) * (j - n / 2)) }
+  }
+
+  /**
+   * Mean distance using fast convolution.
+   *
+   * @param matrix
+   * @return
+   */
+  def distance_convol(matrix: Matrix[T], quantity: Quantity[T]): Double = {
+    val conf = matrix.map { row => row.map(quantity).toArray }.toArray
+    val totPop = conf.flatten.sum
+    val dmat = distanceMatrix(2 * conf.length - 1)
+    val conv = Convolution.convolution2D(conf, dmat)
+    math.sqrt(math.Pi) / (conf.length * totPop * totPop) * MathArrays.ebeMultiply(conv.flatten, conf.flatten).sum
+  }
+
+  /**
+   * Distance kernel
+   *
+   * @param n
+   * @return
+   */
+  def distanceMatrix(n: Int): Array[Array[Double]] = {
+    Array.tabulate(n, n) { (i, j) => Math.sqrt((i - n / 2) * (i - n / 2) + (j - n / 2) * (j - n / 2)) }
+  }
+
 }
diff --git a/scala/model/src/main/scala/fr/iscpif/spacematters/model/stop/SpeilmanStop.scala b/scala/model/src/main/scala/fr/iscpif/spacematters/model/stop/SpeilmanStop.scala
index 414a917..45c2c1b 100644
--- a/scala/model/src/main/scala/fr/iscpif/spacematters/model/stop/SpeilmanStop.scala
+++ b/scala/model/src/main/scala/fr/iscpif/spacematters/model/stop/SpeilmanStop.scala
@@ -17,7 +17,6 @@
 package fr.iscpif.spacematters.model.stop
 
 import fr.iscpif.spacematters.model._
-
 import scala.util.{ Random, Try }
 
 trait SpeilmanStop <: Schelling {