whole lotta changes...

.travis.yml

@@ -4,9 +4,17 @@ python: "3.6"
 install:
   - pip install -U -r requirements.txt
   - pip install -U -r requirements-dev.txt
+  - pip install -U -r requirements-examples.txt
 
 script:
   - python -m pytest --cov-report term --cov=SDM
+  - |
+      python -m ipykernel install --user
+      for i in examples/*.ipynb; do
+         travis_wait 30 jupyter nbconvert --ExecutePreprocessor.timeout=1800 --to markdown --execute --stdout $i > $i.md.travis;
+        # TODO!
+        #diff $i.md.repo $i.md.travis
+      done;
 
 after_success:
   - codecov

SDM/backends/__init__.py

@@ -0,0 +1,6 @@
+"""
+Created at 24.07.2019
+
+@author: Piotr Bartman
+@author: Sylwester Arabas
+"""

SDM/backends/numpy.py

@@ -0,0 +1,137 @@
+"""
+Created at 24.07.2019
+
+@author: Piotr Bartman
+@author: Sylwester Arabas
+"""
+
+import numpy as np
+
+
+# TODO backend.array overrides __getitem__
+
+class Numpy:
+    storage = np.ndarray
+
+    @staticmethod
+    def array(shape, type):
+        if type is float:
+            data = np.full(shape, np.nan, dtype=np.float)
+        elif type is int:
+            data = np.full(shape, -1, dtype=np.int)
+        else:
+            raise NotImplementedError
+        return data
+
+    @staticmethod
+    def from_ndarray(array):
+        if array.ndim > 1:
+            result = array.copy()
+        else:
+            result = np.reshape(array.copy(), (1, -1))
+        return result
+
+    @staticmethod
+    def shuffle(data, length, axis):
+        idx = np.random.permutation(length)
+        Numpy.reindex(data, idx, length, axis=axis)
+
+    @staticmethod
+    def reindex(data, idx, length, axis):
+        if axis == 1:
+            data[:, 0:length] = data[:, idx]
+        else:
+            raise NotImplementedError
+
+    @staticmethod
+    def argsort(idx, data, length):
+        idx[0:length] = data[0:length].argsort()
+
+    @staticmethod
+    def stable_argsort(idx: np.ndarray, data: np.ndarray, length: int):
+        idx[0:length] = data[0:length].argsort(kind='stable')
+
+    @staticmethod
+    def amin(data):
+        result = np.amin(data)
+        return result
+
+    @staticmethod
+    def amax(data):
+        result = np.amax(data)
+        return result
+
+    @staticmethod
+    def transform(data, func, length):
+        data[:length] = np.fromfunction(
+            np.vectorize(func, otypes=(data.dtype,)),
+            (length,),
+            dtype=np.int
+        )
+
+    @staticmethod
+    def foreach(data, func):
+        for i in range(len(data)):
+            func(i)
+
+    @staticmethod
+    def shape(data):
+        return data.shape
+
+    @staticmethod
+    def urand(data, min=0, max=1):
+        data[:] = np.random.uniform(min, max, data.shape)
+
+    # TODO do not create array
+    @staticmethod
+    def remove_zeros(data, idx, length) -> int:
+        for i in range(length):
+            if data[0][idx[0][i]] == 0:
+                idx[0][i] = idx.shape[1]
+        idx.sort()
+        return np.count_nonzero(data)
+
+    @staticmethod
+    def extensive_attr_coalescence(n, idx, length, data, gamma):
+        # TODO in segments
+        for i in range(length // 2):
+            j = 2 * i
+            k = j + 1
+
+            j = idx[j]
+            k = idx[k]
+
+            if n[j] < n[k]:
+                j, k = k, j
+            g = min(gamma[i], n[j] // n[k])
+
+            new_n = n[j] - g * n[k]
+            if new_n > 0:
+                data[:, k] += g * data[:, j]
+            else:  # new_n == 0
+                data[:, j] = g * data[:, j] + data[:, k]
+                data[:, k] = data[:, j]
+
+    @staticmethod
+    def n_coalescence(n, idx, length, gamma):
+        # TODO in segments
+        for i in range(length // 2):
+            j = 2 * i
+            k = j + 1
+
+            j = idx[j]
+            k = idx[k]
+
+            if n[j] < n[k]:
+                j, k = k, j
+            g = min(gamma[i], n[j] // n[k])
+
+            new_n = n[j] - g * n[k]
+            if new_n > 0:
+                n[j] = new_n
+            else:  # new_n == 0
+                n[j] = n[k] // 2
+                n[k] = n[k] - n[j]
+
+
+

SDM/colliders.py

@@ -5,41 +5,52 @@
 @author: Sylwester Arabas
 """
 
-
-import numpy as np
-# import  numba
+from SDM.backends.numpy import Numpy as backend
 
 
 class SDM:
-    def __init__(self, kernel, dt, dv):
-        M = 0  # TODO dependency to state[]
-        N = 1
-        self.probability = lambda sd1, sd2, n_sd: \
-            max(sd1[N], sd2[N]) * kernel(sd1[M], sd2[M]) * dt / dv * n_sd * (n_sd - 1) / 2 / (n_sd//2)
+    def __init__(self, kernel, dt, dv, n_sd):
+        self.probability = lambda sd1n, sd2n, sd1x, sd2x, n_sd: \
+            max(sd1n, sd2n) * kernel(sd1x, sd2x) * dt / dv * n_sd * (n_sd - 1) / 2 / (n_sd // 2)
+        self.rand = backend.array((n_sd // 2,), type=float)
+        self.prob = backend.array((n_sd // 2,), float)
+        self.gamma = backend.array((n_sd // 2,), float)
 
-    # @numba.jit() #TODO
     def __call__(self, state):
-        n_sd = len(state)
-
-        assert np.amin(state.n) > 0
-        if n_sd < 2:
-            return
+        assert state.is_healthy()
 
         # toss pairs
         state.unsort()
 
+        # TODO (segments)
+        # state.sort_by('z', stable=True)  # state.stable_sort_by_segment()
+
         # collide iterating over pairs
-        rand = np.random.uniform(0, 1, n_sd // 2)
+        backend.urand(self.rand)
+
+        backend.transform(self.prob, lambda j: self.probability(state._n[state._idx[2 * j]],
+                                                                state._n[state._idx[2 * j + 1]],
+                                                                state._x[state._idx[2 * j]],
+                                                                state._x[state._idx[2 * j + 1]],
+                                                                state.SD_num),
+                          state.SD_num // 2)
+
+        backend.transform(self.gamma, lambda j: self.prob[j] // 1 + (self.rand[j] < self.prob[j] - self.prob[j] // 1),
+                          state.SD_num // 2)
+
+        # TODO (potential optimisation... some doubts...)
+        # state.sort_by_pairs('n')
 
-        prob_func = np.vectorize(lambda j: self.probability(state[2*int(j)], state[2*int(j)+1], n_sd))
-        prob = np.fromfunction(prob_func, rand.shape, dtype=float)
+        # TODO (when an example with intensive param will be available)
+        # backend.intesive_attr_coalescence(data=state.get_intensive(), gamma=self.gamma)
 
-        # prob = np.empty_like(rand)
-        # for i in range(1, n_sd, 2):
-        #     prob[i // 2] = self.probability(state[idx[i]], state[idx[i - 1]], n_sd)
+        for attrs in state.get_extensive_attrs().values():
+            backend.extensive_attr_coalescence(n=state._n,
+                                               idx=state._idx,
+                                               length=state.SD_num,
+                                               data=attrs,
+                                               gamma=self.gamma)
 
-        gamma = np.floor(prob) + np.where(rand < prob - np.floor(prob), 1, 0)
+        backend.n_coalescence(n=state._n, idx=state._idx, length=state.SD_num, gamma=self.gamma)
 
-        # TODO ! no loops
-        for i in range(1, n_sd, 2):
-            state.collide(i, i - 1, gamma[i//2])
+        state.housekeeping()

SDM/runner.py

@@ -7,6 +7,7 @@
 
 from SDM.stats import Stats
 
+
 class Runner:
 	def __init__(self, state, dynamics):
 		self.state = state