Merge branch 'master' of https://git.frank-ebner.de/FHWS/Indoor

2018-01-17 16:30:04 +01:00
parent e589a40da6 3591106b38
commit e81e8c112d
89 changed files with 6232 additions and 70 deletions
--- a/math/Distributions.h
+++ b/math/Distributions.h
@@ -10,5 +10,6 @@
 #include "distribution/Triangle.h"
 #include "distribution/NormalN.h"
 #include "distribution/Rectangular.h"
+#include "distribution/NormalCDF.h"

 #endif // DISTRIBUTIONS_H
--- a/math/DrawList.h
+++ b/math/DrawList.h
@@ -3,7 +3,7 @@

 #include <vector>

-#include "Random.h"
+#include "random/RandomGenerator.h"
 #include "../Assertions.h"

 /**
@@ -41,13 +41,13 @@ private:
 	std::vector<Entry> elements;

 	/** the used random number generator */
-	RandomGenerator& gen;
+    Random::RandomGenerator& gen;


 private:

 	/** default random generator. fallback */
-	RandomGenerator defRndGen;
+    Random::RandomGenerator defRndGen;


 public:
@@ -63,7 +63,7 @@ public:
 	}

 	/** ctor with custom RandomNumberGenerator */
-	DrawList(RandomGenerator& gen) : cumProbability(0), gen(gen) {
+    DrawList(Random::RandomGenerator& gen) : cumProbability(0), gen(gen) {
 		;
 	}

--- a/math/Random.h
+++ b/math/Random.h
@@ -1,28 +0,0 @@
-#ifndef RANDOM_H
-#define RANDOM_H
-
-#include <cmath>
-#include <random>
-#include "../misc/Time.h"
-
-#ifdef USE_FIXED_SEED
-	#define RANDOM_SEED 1234
-#else
-	#define RANDOM_SEED ( std::chrono::system_clock::now().time_since_epoch() / std::chrono::milliseconds(1) )
-#endif
-
-//using RandomGenerator = std::minstd_rand;
-
-class RandomGenerator : public std::minstd_rand {
-
-public:
-
-	/** ctor with default seed */
-    RandomGenerator() : std::minstd_rand(RANDOM_SEED) {;}
-
-	/** ctor with custom seed */
-    RandomGenerator(result_type seed) : std::minstd_rand(seed) {;}
-
-};
-
-#endif // RANDOM_H
--- a/math/boxkde/BoxGaus.h
+++ b/math/boxkde/BoxGaus.h
@@ -0,0 +1,130 @@
+#pragma once
+
+#include <cmath>
+#include <vector>
+
+#include "Image2D.h"
+#include "BoxSizes.h"
+
+template <class T>
+struct BoxGaus
+{
+    void boxfilter(std::vector<T>& input, size_t w, size_t h, unsigned filterSize)
+    {
+        assertMsg((filterSize % 2) == 1, "filterSize must be odd");
+
+        unsigned radius = filterSize / 2;
+
+        std::vector<T> buffer(input.size());
+        boxBlur(input, buffer, w, h, radius);
+        boxBlur(buffer, input, w, h, radius);
+    }
+    
+    void approxGaus(Image2D<T>& input, T sigmaX, T sigmaY, unsigned nFilt)
+    {
+        approxGaus(input.data(), input.width, input.height, sigmaX, sigmaY, nFilt);
+    }
+
+	void approxGaus(std::vector<T>& input, size_t w, size_t h, T sigmaX, T sigmaY, unsigned nFilt)
+	{
+        BoxSizes<T> bsX(sigmaX, nFilt);
+        BoxSizes<T> bsY(sigmaY, nFilt);
+        std::vector<T> buffer(input.size());
+
+        assertMsg((2 * bsX.wl + 1 < w) && (2 * bsX.wl + 1 < h), "Box-Filter size in X direction is too big");
+        assertMsg((2 * bsX.wu + 1 < w) && (2 * bsX.wu + 1 < h), "Box-Filter size in X direction is too big");
+        assertMsg((2 * bsY.wl + 1 < w) && (2 * bsY.wl + 1 < h), "Box-Filter size in Y direction is too big");
+        assertMsg((2 * bsY.wu + 1 < w) && (2 * bsY.wu + 1 < h), "Box-Filter size in Y direction is too big");
+
+        // if equal, we can save some cond's inside the loop
+        if (bsX.m == bsY.m)
+        {
+            const size_t m = bsX.m;
+
+            for (size_t i = 0; i < m; i++)
+            {
+                boxBlur(input, buffer, w, h, bsY.wl);
+                boxBlur(buffer, input, w, h, bsX.wl);
+            }
+
+            for (size_t i = 0; i < nFilt - m; i++)
+            {
+                boxBlur(input, buffer, w, h, bsY.wu);
+                boxBlur(buffer, input, w, h, bsX.wu);
+            }
+        }
+        else
+        {
+            for (size_t i = 0; i < nFilt; i++)
+            {
+                boxBlur(input, buffer, w, h, (i < bsY.m ? bsY.wl : bsY.wu) );
+                boxBlur(buffer, input, w, h, (i < bsX.m ? bsX.wl : bsX.wu));
+            }
+        }
+
+	}
+
+private:
+    void boxBlur(const std::vector<T> &src, std::vector<T> &dst, size_t w, size_t h, size_t r)
+    {
+        T iarr = (T)1.0 / (r + r + 1);
+
+        for (size_t i = 0; i < w; i++)
+        {
+            // Init indexes
+            size_t ti = i;
+            size_t li = ti;
+            size_t ri = ti + r*w;
+
+            // Init values
+            T fv = src[ti]; 			// first values
+            T lv = src[ti + w*(h - 1)]; // last values
+            T val = fv * (r + 1); // overhang over image border
+
+            for (size_t j = 0; j < r; j++)
+            {
+                val += src[ti + j*w];  // col sum
+            }
+
+            // <20>berhangbereich links vom Bild
+            for (size_t j = 0; j <= r; j++)
+            {
+                val += src[ri] - fv;
+                dst[i + j*w] = val * iarr;
+
+                ri += w;
+                ti += w;
+            }
+
+            // Bildbereich
+            for (size_t j = r + 1; j < h - r; j++)
+            {
+                val += src[ri] - src[li];
+                dst[i + j*w] = val * iarr;
+
+                li += w;
+                ri += w;
+                ti += w;
+            }
+
+            // <20>berhangbereich rechts vom Bild
+            for (size_t j = h - r; j < h; j++)
+            {
+                val += lv - src[li];
+                dst[i + j*w] = val * iarr;
+
+                li += w;
+                ti += w;
+            }
+        }
+
+        int test = 0;
+    }
+};
+
+
+
+
+
+
+
--- a/math/boxkde/BoxSizes.h
+++ b/math/boxkde/BoxSizes.h
@@ -0,0 +1,85 @@
+#pragma once
+
+#include "DataStructures.h"
+
+template<typename T>
+struct BoxSizes
+{
+    static_assert(std::is_floating_point<T>::value, "This class only works with floats.");
+
+    T sigma, sigmaActual;
+    T wIdeal, mIdeal;
+
+    unsigned n, m, wl, wu;
+
+    BoxSizes(T sigma, unsigned n)
+        : sigma(sigma), n(n)
+    {
+        assertMsg(sigma >= 0.8, "Sigma values below about 0.8 cannot be represented");
+
+        wIdeal = sqrt(12 * sigma*sigma / n + 1); // Ideal averaging filter width
+
+        // wl is first odd valued integer less than wIdeal
+        wl = (unsigned)floor(wIdeal);
+        if (wl % 2 == 0)
+            wl = wl - 1;
+
+        // wu is the next odd value > wl
+        wu = wl + 2;
+
+        // Compute m.Refer to the tech note for derivation of this formula
+        mIdeal = (12 * sigma*sigma - n*wl*wl - 4 * n*wl - 3 * n) / (-4 * wl - 4);
+        m = (unsigned)round(mIdeal);
+
+        assertMsg(!(m > n || m < 0), "calculation of m has failed");
+
+        // Compute actual sigma that will be achieved
+        sigmaActual = sqrt((m*wl*wl + (n - m)*wu*wu - n) / (T)12.0);
+    }
+};
+
+
+template<typename T>
+struct ExBoxSizes
+{
+    static_assert(std::is_floating_point<T>::value, "This class only works with floats.");
+
+    T sigma, sigma_actual;
+    unsigned n, r;
+    T alpha, c, c1, c2;
+
+    T r_f;
+
+    // Special case for h == 1
+    ExBoxSizes(T sigma, unsigned n)
+        : sigma(sigma), n(n)
+    {
+        T var = sigma*sigma;
+        r_f = 0.5*sqrt((12 * var) / n + 1) - T(0.5);
+        r = (unsigned)std::floor(r_f);
+
+        alpha = (2 * r + 1) * ( (r*r + r - 3*var/n) / (6 * (var/n - (r+1)*(r+1))) );
+        
+        c1 = alpha / (2*alpha + 2*r + 1);
+        c2 = (1 - alpha) / (2 * alpha + 2 * r + 1);
+        c = c1 + c2;
+    }
+
+    ExBoxSizes(T sigma, unsigned d, T h)
+        : sigma(sigma), n(d)
+    {
+        T v = sigma*sigma;
+        r_f = sqrt((12 * v) / n + 1)/(2*h) - T(0.5); // (7)
+        r = (unsigned)std::floor(std::max(T(0), r_f));
+
+        alpha = (2 * r + 1) * (((r*r + r) - (v*3*h)/(d*h*h*h)) / (6 * ( v/(d*h*h) - (r+1)*(r+1)) )); // (8) (14)
+        c1 = alpha / (h*(2 * r + 2 * alpha + 1));       // (8) (13)
+        c2 = (1 - alpha) / (h*(2 * r + 2 * alpha + 1)); // (8) (13)
+
+        c = c1 + c2;
+
+        T lambda = h*(2*r + 1 + 2*alpha); // (8)
+        T variance_actual = (d*h*h*h) / (3 * lambda) * (2*r*r*r + 3*r*r + r + 6*alpha*(r+1)*(r+1)); // (14)
+        sigma_actual = sqrt(variance_actual);
+    }
+};
--- a/math/boxkde/DataStructures.h
+++ b/math/boxkde/DataStructures.h
@@ -0,0 +1,115 @@
+#pragma once
+
+#include <cassert>
+#include <limits>
+#include <memory>
+#include <sstream>
+#include <vector>
+
+template <class T>
+struct BoundingBox
+{
+	static_assert(std::is_arithmetic<T>::value, "This class only works with floats or integers.");
+
+	T MinX, MaxX, MinY, MaxY;
+	
+	BoundingBox(T MinX = std::numeric_limits<T>::max(),
+		        T MaxX = std::numeric_limits<T>::lowest(),
+		        T MinY = std::numeric_limits<T>::max(),
+		        T MaxY = std::numeric_limits<T>::lowest())
+		: MinX(MinX), MaxX(MaxX), MinY(MinY), MaxY(MaxY)
+	{ }
+
+	T width () const { return MaxX - MinX; }
+	T heigth() const { return MaxY - MinY; }
+	T area  () const { return width()*heigth(); }
+
+	bool isInside(T x, T y) const { return (x >= MinX && x <= MaxX) && (y >= MinY && y <= MaxY); }
+
+    // Expands the size of the BB if the given values are extreme
+    void expand(T x, T y)
+    {
+        if (x < MinX) MinX = x;
+        if (x > MaxX) MaxX = x;
+        if (y < MinY) MinY = y;
+        if (y > MaxY) MaxY = y;
+    }
+
+    // Enlarges the BB in both direction along an axis.
+    void inflate(T szX, T szY)
+    {
+        MinX -= szX;
+        MinY -= szY;
+
+        MaxX += szX;
+        MaxY += szY;
+    }
+};
+
+
+template <class T>
+struct Point2D
+{
+	static_assert(std::is_arithmetic<T>::value, "This class only works with floats and integers.");
+
+	T X, Y;
+
+	Point2D(T x = 0, T y = 0)
+		: X(x), Y(y)
+	{ }
+};
+
+template <class T>
+struct Size2D
+{
+    static_assert(std::is_arithmetic<T>::value, "This class only works with floats and integers.");
+
+    T sX, sY;
+
+
+    Size2D(T all)
+        : sX(all), sY(all)
+    { }
+
+    Size2D(T x = 0, T y = 0)
+        : sX(x), sY(y)
+    { }
+};
+
+
+#ifdef NDEBUG
+#define assertCond(_EXPR_) (false)
+#define assertThrow(_arg_) ((void)0)
+#define assert_throw(...)  ((void)0)
+#define assertMsg(...)     ((void)0)
+#else
+// Evaluates the expression. Ifdef NDEBUG returns always false
+    #define assertCond(_EXPR_) (!(_EXPR_))
+// Throws a excpetion with the argument as message by prepending the current file name and line number
+    #define assertThrow(_std_string_) assert_throw( (_std_string_), __FILE__, __LINE__)
+    inline void assert_throw(const std::string& message, const char* file, int line)
+    {
+        std::stringstream ss;
+        ss << file << ":" << line << ": " << message;
+        throw std::invalid_argument(ss.str());
+    }
+
+    #define assertMsg(_EXPR_, _MSG_) if (!(_EXPR_)) assertThrow(_MSG_)
+#endif
+
+// ostream overloads
+template<typename T>
+std::ostream& operator<<(std::ostream& os, const Point2D<T>& pt)
+{
+    return os << "(" << pt.X << "; " << pt.Y << ")";
+}
+
+template<typename T>
+std::ostream& operator<<(std::ostream& os, const BoundingBox<T>& bb)
+{
+    return os << "(X: " << bb.MinX << " - " << bb.MaxX << ";" 
+              << " Y: " << bb.MinY << " - " << bb.MaxY << ")";
+}
+
+
+
--- a/math/boxkde/GausLib.cpp
+++ b/math/boxkde/GausLib.cpp
@@ -0,0 +1,83 @@
+#include "GausLib.h"
+
+#include "benchmark.h"
+#include "DataStructures.h"
+#include "Image2D.h"
+#include "BoxGaus.h"
+#include "Grid2D.h"
+#include <cstring>
+
+// X sind 2D Samples [x,x,x,.... y,y,y]
+// sizeX = Anzahl Samples
+// hx, hy Bandbreite KDE (default=100)
+// nBins = Anzahl der Bins in die jeweilige Richtung
+// boundingBox ist Größe des Gebäudes (40 m x 60 m)
+// algorithm = define method  
+
+void gausKde2D_simple(float* X, int sizeX, float* weights, float hx, float hy, int nBinsX, int nBinsY, c_bbox* boundingBox, int algorithm, double* out_maxValue, double* out_maxPosX, double* out_maxPosY, double* out_runTimeInNS, float* out_density)
+{
+    BenchResult bench = benchmarkEx("", 1, [&](){
+        // Create bounding box
+        BoundingBox<float> bb(boundingBox->minX, boundingBox->maxX, boundingBox->minY, boundingBox->maxY);
+
+        // Create histogram
+        Grid2D<float> grid(bb, nBinsX, nBinsY);
+        for (size_t i = 0; i < sizeX; i++)
+        {
+            float x = X[i];
+            float y = X[i + sizeX];
+            float w = weights[i];
+
+            grid.add(x, y, w);
+        }
+
+        int nFilt = 3;
+
+        float sigmaX = hx / grid.binSizeX;
+        float sigmaY = hy / grid.binSizeY;
+
+        // Apply gaus filter
+        switch (algorithm)
+        {
+            default:
+            case ALGO_BOX:
+            {
+                BoxGaus<float> boxGaus;
+                boxGaus.approxGaus(grid.image(), sigmaX, sigmaY, nFilt);
+                break;
+            }
+            case ALGO_EXBOX:
+            {
+                //ExtendedBox<float> boxGaus;
+                //boxGaus.approxGaus(grid.image(), sigmaX, sigmaY, nFilt);
+                break;
+            }
+            case ALGO_BOX_SIMD:
+            {
+                //BoxGausSIMD<float> boxGaus;
+                //boxGaus.approxGaus(grid.image(), sigmaX, sigmaY, nFilt);
+                break;
+            }
+            case ALGO_BOX_CL:
+            {
+                //clBox->writeBuffer(grid.image().data().data());
+                //clBox->execute(sigmaX, nFilt);
+                //clBox->readBuffer(grid.image().data().data());
+                break;
+            }
+        }
+
+
+        Point2D<float> maxPos;
+        *out_maxValue = grid.maximum(maxPos);
+        *out_maxPosX = static_cast<double>(maxPos.X);
+        *out_maxPosY = static_cast<double>(maxPos.Y);
+
+        if (out_density)
+        {
+            memcpy(out_density, grid.image().data().data(), sizeof(float)*grid.image().data().size());
+        }
+    });
+
+    *out_runTimeInNS = bench.min;
+}
--- a/math/boxkde/GausLib.h
+++ b/math/boxkde/GausLib.h
@@ -0,0 +1,34 @@
+// The following ifdef block is the standard way of creating macros which make exporting 
+// from a DLL simpler. All files within this DLL are compiled with the GAUSLIB_EXPORTS
+// symbol defined on the command line. This symbol should not be defined on any project
+// that uses this DLL. This way any other project whose source files include this file see 
+// GAUSLIB_API functions as being imported from a DLL, whereas this DLL sees symbols
+// defined with this macro as being exported.
+
+    typedef struct c_bbox {
+        double minX;
+        double maxX;
+        double minY;
+        double maxY;
+        double minZ;
+        double maxZ;
+    } c_bbox;
+
+    
+#define ALGO_BOX 0
+#define ALGO_EXBOX 1
+#define ALGO_BOX_SIMD 2
+#define ALGO_BOX_CL 3
+    
+    void gausKde2D_simple(float* X, int sizeX,
+                          float* weights,
+                          float hx, float hy,
+                          int nBinsX, int nBinsY,
+                          c_bbox* boundingBox,
+                          int algorithm,
+                          double* out_maxValue,
+                          double* out_maxPosX, double* out_maxPosY,
+                          double* out_runTimeInNS,
+                          float* out_density);
+
+   
--- a/math/boxkde/Grid2D.h
+++ b/math/boxkde/Grid2D.h
@@ -0,0 +1,211 @@
+#pragma once
+
+#include <sstream>
+
+#include "DataStructures.h"
+#include "Image2D.h"
+
+template<class T>
+struct Grid2D
+{
+	static_assert(std::is_floating_point<T>::value, "Grid2D only supports float values");
+
+    BoundingBox<T> bb;
+    size_t numBinsX, numBinsY;
+    T binSizeX, binSizeY;
+private:
+	      Image2D<T> data;
+
+public:
+
+    Grid2D(){   } //TODO: fast hack
+
+    Grid2D(BoundingBox<T> bb, size_t numBins)
+        : Grid2D(bb, numBins, numBins)
+    {
+    }
+
+	Grid2D(BoundingBox<T> bb, size_t numBinsX, size_t numBinsY)
+		: bb(bb),
+          numBinsX(numBinsX),
+          numBinsY(numBinsY),
+          binSizeX((bb.width()) / (numBinsX-1)),
+          binSizeY((bb.heigth()) / (numBinsY-1)),
+          data(numBinsX, numBinsY)
+	{
+	}
+
+          Image2D<T>& image()       { return data; }
+    const Image2D<T>& image() const { return data; }
+
+          T& operator() (size_t x, size_t y)       { return data(x, y); }
+    const T& operator() (size_t x, size_t y) const { return data(x, y); }
+
+    void clear() { data.clear(); }
+
+	void fill(const std::vector<Point2D<T>>& samples)
+	{
+		assertMsg(!samples.empty(), "Samples must be non-empty");
+
+		data.clear();
+
+		const T weight = T(1.0) / samples.size();
+		for (const auto& pt : samples)
+		{
+			add(pt.X, pt.Y, weight);
+		}
+	}
+
+	void fill(const std::vector<Point2D<T>>& samples, const std::vector<T>& weights)
+	{
+        assertMsg(!samples.empty(), "Samples must be non-empty");
+        assertMsg(weights.size() == samples.size(), "Weights must have the same size as samples");
+
+		data.clear();
+
+		for (size_t i = 0; i < samples.size(); i++)
+		{
+			add(samples[i].X, samples[i].Y, weights[i]);
+		}
+	}
+
+	Point2D<size_t> add(T x, T y, T w)
+	{
+        if (assertCond(bb.isInside(x,y)))
+        {
+            std::stringstream ss;
+            ss << "Point " << Point2D<T>(x, y) << " is out of bounds. " << bb;
+            assertThrow(ss.str());
+        }
+
+		return add_simple_bin(x, y, w);
+		//return add_linear_bin(x, y, w);
+	}
+
+    void add_linear(T x, T y, T w)
+    {
+        if (assertCond(bb.isInside(x, y)))
+        {
+            std::stringstream ss;
+            ss << "Point " << Point2D<T>(x, y) << " is out of bounds. " << bb;
+            assertThrow(ss.str());
+        }
+
+        add_linear_bin(x, y, w);
+    }
+
+    T fetch(T x, T y) const
+    {
+        size_t bin_x = (size_t)((x - bb.MinX) / binSizeX);
+        size_t bin_y = (size_t)((y - bb.MinY) / binSizeY);
+
+        //if (bin_x == data.width)  bin_x--;
+        //if (bin_y == data.height) bin_y--;
+
+        return data(bin_x, bin_y);
+    }
+
+    // Returns the summation of all bin values
+    T sum() const
+    {
+        return data.sum();
+    }
+
+    // Takes a point in input space and converts it to grid space coordinates
+    Point2D<size_t> asGridSpace(T x, T y) const
+    {
+        size_t bin_x = (size_t)((x - bb.MinX) / binSizeX);
+        size_t bin_y = (size_t)((y - bb.MinY) / binSizeY);
+
+        if (bin_x == data.width)  bin_x--;
+        if (bin_y == data.height) bin_y--;
+
+        return Point2D<size_t>(bin_x, bin_y);
+    }
+
+    // Takes a Size2D in input space and converts it to grid space coordiantes
+    Size2D<size_t> asGridSpace(Size2D<T> sz) const
+    {
+        return Size2D<size_t>(sz.sX / binSizeX, sz.sY / binSizeY);
+    }
+
+    // Takes a point in grid space and converts it to input space coordinates
+    Point2D<T> asInputSpace(Point2D<size_t> pt) const
+    {
+        return asInputSpace(pt.X, pt.Y);
+    }
+
+    // Takes a point in grid space and converts it to input space coordinates
+    Point2D<T> asInputSpace(size_t x, size_t y) const
+    {
+        return Point2D<T>(x * binSizeX + bb.MinX + T(0.5)*binSizeX,
+                          y * binSizeY + bb.MinY + T(0.5)*binSizeY);
+    }
+
+
+    T maximum(Point2D<T>& pt) const
+    {
+        Point2D<size_t> gridPt;
+
+        T maxValue = image().maximum(gridPt);
+        pt = asInputSpace(gridPt);
+        return maxValue;
+    }
+
+private:
+    Point2D<size_t> add_simple_bin(T x, T y, T w)
+	{
+		size_t bin_x = (size_t)((x - bb.MinX) / binSizeX);
+		size_t bin_y = (size_t)((y - bb.MinY) / binSizeY);
+
+        if (bin_x == data.width)  bin_x--;
+        if (bin_y == data.height) bin_y--;
+
+		data(bin_x, bin_y) += w;
+
+        return Point2D<size_t>(bin_x, bin_y);
+	}
+
+    void add_linear_bin(T x, T y, T w)
+	{
+        T xpos1 = (x - bb.MinX) / binSizeX;
+        T xpos2 = (y - bb.MinY) / binSizeY;
+        size_t ix1 = (size_t)floor(xpos1);
+        size_t ix2 = (size_t)floor(xpos2);
+        T fx1 = xpos1 - ix1;
+        T fx2 = xpos2 - ix2;
+        
+        const size_t ixmin1 = 0;
+        const size_t ixmax1 = numBinsX - 2;
+        const size_t ixmin2 = 0;
+        const size_t ixmax2 = numBinsY - 2;
+
+        if ( ixmin1 <= ix1 && ixmin2 <= ix2 && ix2 <= ixmax2) 
+        {
+            data(ix1, ix2) += w*(1 - fx1)*(1 - fx2);
+            data(ix1+1, ix2) += w*fx1*(1 - fx2);
+            data(ix1, ix2+1) += w*(1 - fx1)*fx2;
+            data(ix1 + 1, ix2 + 1) += w*fx1*fx2;
+        }
+        else if (ix1 == ixmax1 + 1 && ixmin2 <= ix2 && ix2 <= ixmax2) 
+        {
+            // rechts
+            data(ix1, ix2) += w*(1 - fx1)*(1 - fx2);
+            data(ix1, ix2+1) += w*(1 - fx1)*fx2;
+        }
+        else if (ixmin1 <= ix1 && ix1 <= ixmax1 && ix2 == ixmax2 + 1) 
+        {
+            // unten
+            data(ix1, ix2) += w*(1 - fx1)*(1 - fx2);
+            data(ix1+1, ix2) += w*fx1*(1 - fx2);
+        }
+        else if (ix1 == ixmax1 + 1 && ix2 == ixmax2 + 1) 
+        {
+            // rechts-unten
+            data(ix1, ix2) += w*(1 - fx1)*(1 - fx2);
+        }
+	}
+};
+
+template struct Grid2D<float>;
+template struct Grid2D<double>;
--- a/math/boxkde/Image2D.h
+++ b/math/boxkde/Image2D.h
@@ -0,0 +1,187 @@
+#pragma once
+
+#include <algorithm>
+#include <cassert>
+#include <numeric>
+#include <vector>
+
+#include "DataStructures.h"
+
+enum struct LineDirection { X, Y };
+
+template <typename TValue>
+struct ImageView2D
+{
+    static_assert(std::is_arithmetic<TValue>::value, "Image only supports integers or floats");
+
+    // forward declaration
+    //enum struct LineDirection;
+    template<LineDirection D> struct LineView;
+    template<LineDirection D> struct ConstLineView;
+
+    size_t width, height, size;
+protected:
+    TValue* values; // contains image data row-wise 
+
+public:
+    ImageView2D()
+        : width(0), height(0), size(0), values(nullptr)
+    { }
+
+    ImageView2D(size_t width, size_t height)
+        : width(width), height(height), size(width*height), values(nullptr)
+    { }
+
+    ImageView2D(size_t width, size_t height, TValue* data)
+        : width(width), height(height), size(width*height), values(data)
+    { }
+
+    inline       TValue& operator() (size_t x, size_t y)       { return values[indexFromCoord(x, y)]; }
+    inline const TValue& operator() (size_t x, size_t y) const { return values[indexFromCoord(x, y)]; }
+
+    TValue* val_begin() { return values; }
+    TValue* val_end  () { return values + size; }
+
+    const TValue* val_begin() const { return values; }
+    const TValue* val_end  () const { return values + size; }
+
+    inline size_t indexFromCoord(size_t x, size_t y) const
+    {
+        assertMsg(x < width && y < height, "(x,y) out of bounds");
+        return y * width + x;
+    }
+
+    Point2D<size_t> coordFromIndex(size_t index) const
+    {
+        assertMsg(index < size, "Index out of bounds");
+        return Point2D<size_t>(index % width, index / width);
+    }
+
+    Point2D<size_t> maximum() const
+    {
+        size_t maxValueIndex = std::distance(val_begin(), std::max_element(val_begin(), val_end()));
+        return coordFromIndex(maxValueIndex);
+    }
+
+    TValue maximum(Point2D<size_t>& pt) const
+    {
+        size_t maxValueIndex = std::distance(val_begin(), std::max_element(val_begin(), val_end()));
+        pt = coordFromIndex(maxValueIndex);
+        return values[maxValueIndex];
+    }
+
+    void clear()
+    {
+        fill(TValue(0));
+    }
+
+    void fill(TValue value)
+    {
+        std::fill(val_begin(), val_end(), value);
+    }
+    
+    void assign(const ImageView2D<TValue>& other)
+    {
+        assertMsg(size == other.size, "Other must be of the same size as this");
+        std::copy(other.val_begin(), other.val_end(), val_begin());
+    }
+
+    TValue sum() const
+    {
+        return std::accumulate(val_begin(), val_end(), TValue(0));
+    }
+
+    // Returns a transposed view of this image without copying any data.
+    ImageView2D<TValue> transpose() const
+    {
+        return ImageView2D<TValue>(height, width, values);
+    }
+
+    // Returns the value at (x,y). Throws if out of bounds.
+    const TValue& at(size_t x, size_t y) const
+    {
+        return values[indexFromCoord(x, y)];
+    }
+
+    // Returns the value at (x,y) but falls back to default if out of bounds.
+    const TValue& get(size_t x, size_t y, TValue dflt = 0) const
+    {
+        if (x >= width || y >= height)
+            return dflt;
+        else
+            return values[indexFromCoord(x, y)];
+    }
+    
+    std::vector<Point2D<size_t>> findLocalMaxima(int radiusX = 1, int radiusY = 1) const
+    {
+        std::vector<Point2D<size_t>> result;
+        findLocalMaxima(result, radiusX, radiusY);
+        return result;
+    }
+
+    void findLocalMaxima(std::vector<Point2D<size_t>>& result, int radiusX = 1, int radiusY = 1) const
+    {
+        assertMsg(radiusX > 0, "Search radius must be greater than 0.");
+
+        for (size_t y = 0; y < height; y++)
+        {
+            for (size_t x = 0; x < width; x++)
+            {
+                TValue val = at(x, y);
+                bool lclMax = true;
+
+                for (int ry = -radiusY; ry <= radiusY; ry++)
+                {
+                    for (int rx = -radiusX; rx <= radiusX; rx++)
+                    {
+                        TValue other = get(x + rx, y + ry);
+                        if (!(rx == 0 && ry == 0) && val <= other)
+                        {
+                            lclMax = false;
+                        }
+                    }
+                }
+
+                if (lclMax)
+                {
+                    result.emplace_back(x, y);
+                }
+            }
+        }
+    }
+};
+
+template <class TValue>
+struct Image2D : public ImageView2D<TValue>
+{
+	static_assert(std::is_arithmetic<TValue>::value, "Image only supports integers or floats");
+
+	std::vector<TValue> values_vec;
+
+	Image2D()
+	{
+	}
+
+	Image2D(size_t width, size_t height)
+        : ImageView2D<TValue>(width, height), values_vec(width * height)
+	{
+        this->values = values_vec.data();
+    }
+
+    Image2D(size_t width, size_t height, std::vector<TValue>& data)
+        : ImageView2D<TValue>(width, height), values_vec(data)
+    {
+        assertMsg(data.size() == width*height, "Sizes must be the same");
+        this->values = values_vec.data();
+    }
+
+          std::vector<TValue>& data()       { return values_vec; }
+    const std::vector<TValue>& data() const { return values_vec; }
+};
+
+
+template struct ImageView2D<float>;
+template struct ImageView2D<double>;
+
+template struct Image2D<float>;
+template struct Image2D<double>;
--- a/math/boxkde/benchmark.h
+++ b/math/boxkde/benchmark.h
@@ -0,0 +1,108 @@
+#pragma once
+#include <chrono>
+#include <string>
+#include <iostream>
+#include <vector>
+#include <algorithm>
+#include <numeric>
+
+template <typename F>
+static void benchmark(std::string name, size_t count, F&& lambda)
+{
+	auto start = std::chrono::high_resolution_clock::now();
+
+	for (size_t i = 0; i < count; i++)
+	{
+		lambda();
+	}
+
+	auto duration = std::chrono::duration_cast<std::chrono::microseconds>(std::chrono::high_resolution_clock::now() - start);
+	long long durationAVG = duration.count() / count;
+
+	std::cout << "Function " << name << " took avg. " << durationAVG << "us (" << durationAVG / 1000.0f << "ms) over " << count << " runs." << std::endl;
+}
+
+
+struct BenchResult
+{
+    const long long min, max;
+    const double mean, median;
+
+    BenchResult()
+        : min(0), max(0), mean(0), median(0)
+    {}
+
+    BenchResult(long long min, long long max, double mean, double median)
+        : min(min), max(max), mean(mean), median(median)
+    {}
+};
+
+template <typename F>
+static BenchResult benchmarkEx(std::string name, size_t count, F&& lambda)
+{
+    std::vector<long long> durations(count);
+    
+    for (size_t i = 0; i < count; i++)
+    {
+        auto start = std::chrono::high_resolution_clock::now();
+        lambda();
+        auto stop = std::chrono::high_resolution_clock::now();
+        auto duration = std::chrono::duration_cast<std::chrono::nanoseconds>(stop - start);
+        durations[i] = duration.count();
+    }
+
+    long long min = *std::min_element(durations.begin(), durations.end());
+    long long max = *std::max_element(durations.begin(), durations.end());
+    double average = std::accumulate(durations.begin(), durations.end(), 0.0) / durations.size();
+
+    double median;
+    std::sort(durations.begin(), durations.end());
+
+    if (durations.size() % 2 == 0)
+        median = (double) (durations[durations.size() / 2 - 1] + durations[durations.size() / 2]) / 2;
+    else
+        median = (double) durations[durations.size() / 2];
+
+    std::cout << "Function " << name << " took avg. " << average << "ns (" << average / 1000.0f << "us) over " << count << " runs." << std::endl;
+
+    return BenchResult(min, max, average, median);
+}
+
+struct StopWatch
+{
+	typedef std::chrono::high_resolution_clock clock;
+	typedef clock::time_point time_point;
+
+	time_point startTime; // Point in time when start() was called
+	time_point lapTime;   // Point in time when operator() was called
+
+public:
+	StopWatch()
+	{
+		reset();
+	}
+
+	void reset()
+	{
+		startTime = clock::now();
+		lapTime = startTime;
+	}
+
+    std::chrono::microseconds stop()
+	{
+		auto duration = std::chrono::duration_cast<std::chrono::microseconds>(clock::now() - startTime);
+
+		std::cout << "Total time: " << duration.count() << "us (" << duration.count() / 1000.0f << "ms)" << std::endl;
+
+        return duration;
+	}
+
+	void operator()(const std::string& str = "")
+	{
+		auto duration = std::chrono::duration_cast<std::chrono::microseconds>(clock::now() - lapTime);
+
+		std::cout << str << (str.empty() ? "" : " ") << "took " << duration.count() << "us (" << duration.count() / 1000.0f << "ms)" << std::endl;
+
+		lapTime = clock::now();
+	}
+};
--- a/math/distribution/Exponential.h
+++ b/math/distribution/Exponential.h
@@ -3,7 +3,7 @@

 #include <cmath>
 #include <random>
-#include "../Random.h"
+#include "../random/RandomGenerator.h"

 namespace Distribution {

@@ -14,7 +14,7 @@ namespace Distribution {

 		const T lambda;

-		RandomGenerator gen;
+        Random::RandomGenerator gen;
 		std::exponential_distribution<T> dist;

 	public:
--- a/math/distribution/KernelDensity.h
+++ b/math/distribution/KernelDensity.h
@@ -8,7 +8,7 @@
 #include <eigen3/Eigen/Dense>

 #include "../../Assertions.h"
-#include "../Random.h"
+#include "../random/RandomGenerator.h"


 namespace Distribution {
--- a/math/distribution/Normal.h
+++ b/math/distribution/Normal.h
@@ -3,7 +3,7 @@

 #include <cmath>
 #include <random>
-#include "../Random.h"
+#include "../random/RandomGenerator.h"
 #include "../../Assertions.h"

 namespace Distribution {
@@ -17,7 +17,7 @@ namespace Distribution {
 		const T sigma;
 		const T _a;

-		RandomGenerator gen;
+        Random::RandomGenerator gen;
 		std::normal_distribution<T> dist;

 	public:
--- a/math/distribution/NormalCDF.h
+++ b/math/distribution/NormalCDF.h
@@ -0,0 +1,65 @@
+#ifndef NORMALCDF_H
+#define NORMALCDF_H
+
+#include <algorithm>
+#include <random>
+#include "../../Assertions.h"
+
+namespace Distribution {
+
+        /** cumulative density version of the normal distribution */
+        template <typename T> class NormalCDF {
+
+        private:
+
+            const T mu;
+            const T sigma;
+
+            static T RationalApproximation(T t)
+            {
+                // Abramowitz and Stegun formula 26.2.23.
+                // The absolute value of the error should be less than 4.5 e-4.
+                T c[] = {2.515517, 0.802853, 0.010328};
+                T d[] = {1.432788, 0.189269, 0.001308};
+                return t - ((c[2]*t + c[1])*t + c[0]) /
+                           (((d[2]*t + d[1])*t + d[0])*t + 1.0);
+            }
+
+        public:
+
+                /** create a new normally distributed CDF */
+                NormalCDF(const T mu, const T sigma) : mu(mu), sigma(sigma) {
+                        ;
+                }
+
+                /** get the probability for val within the underlying CDF */
+                T getProbability(const T val) const {
+                        return getProbability(mu, sigma, val);
+                }
+
+                /** calculate the probability within the underlying CDF */
+                static T getProbability(const T mu, const T sigma, const T val) {
+                    return (1.0 + std::exp( (val - mu) / (sigma * std::sqrt(2)) ) ) / 2.0;
+                }
+
+                /** get the inverse CDF (https://en.wikipedia.org/wiki/Probit)*/
+                static T getProbit(T p){
+
+                    Assert::isBetween(p, 0.0, 1.0, "value not between");
+
+                    // See: https://www.johndcook.com/blog/normal_cdf_inverse/
+                    if (p < 0.5)
+                    {
+                        // F^-1(p) = - G^-1(p)
+                        return -RationalApproximation( sqrt(-2.0*log(p)) );
+                    }
+                    else
+                    {
+                        // F^-1(p) = G^-1(1-p)
+                        return RationalApproximation( sqrt(-2.0*log(1-p)) );
+                    }
+                }
+    };
+}
+
+#endif // NORMALCDF_H
--- a/math/distribution/NormalN.h
+++ b/math/distribution/NormalN.h
@@ -7,7 +7,7 @@
 #include <eigen3/Eigen/Dense>

 #include "../../Assertions.h"
-#include "../Random.h"
+#include "../random/RandomGenerator.h"

 namespace Distribution {

@@ -24,7 +24,7 @@ namespace Distribution {
        const Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eigenSolver;
        Eigen::MatrixXd transform; //can i make this const?

-        RandomGenerator gen;
+        Random::RandomGenerator gen;
        std::normal_distribution<> dist;

    public:
--- a/math/distribution/Rectangular.h
+++ b/math/distribution/Rectangular.h
@@ -4,7 +4,7 @@

 #include <cmath>
 #include <random>
-#include "../Random.h"
+#include "../random/RandomGenerator.h"
 #include "../../Assertions.h"
 #include "Normal.h"

--- a/math/distribution/Region.h
+++ b/math/distribution/Region.h
@@ -3,7 +3,7 @@

 #include <cmath>
 #include <random>
-#include "../Random.h"
+#include "../random/RandomGenerator.h"
 #include "../../Assertions.h"
 #include "Normal.h"

--- a/math/distribution/Triangle.h
+++ b/math/distribution/Triangle.h
@@ -3,7 +3,7 @@

 #include <cmath>
 #include <random>
-#include "../Random.h"
+#include "../random/RandomGenerator.h"
 #include "../../Assertions.h"
 #include "Normal.h"

--- a/math/distribution/Uniform.h
+++ b/math/distribution/Uniform.h
@@ -3,7 +3,7 @@

 #include <cmath>
 #include <random>
-#include "../Random.h"
+#include "../random/RandomGenerator.h"

 #include <type_traits>

@@ -14,7 +14,7 @@ namespace Distribution {

 	private:

-		RandomGenerator gen;
+        Random::RandomGenerator gen;

 		/** depending on T, Dist is either a uniform_real or uniform_int distribution */
 		typedef typename std::conditional< std::is_floating_point<T>::value, std::uniform_real_distribution<T>, std::uniform_int_distribution<T> >::type Dist;
--- a/math/random/DrawList.h
+++ b/math/random/DrawList.h
@@ -0,0 +1,169 @@
+#ifndef K_MATH_RANDOM_DRAWLIST_H
+#define K_MATH_RANDOM_DRAWLIST_H
+
+#include <vector>
+
+/**
+ * add entries to a list and be able to draw from them depending oh their probability
+ */
+namespace Random {
+
+	/**
+	 * this class represents one entry within a DrawList.
+	 * such an entry consists of userData denoted by the template argument
+	 * and a probability.
+	 */
+	template <typename Entry> struct DrawListEntry {
+
+	template <typename> friend class DrawList;
+	friend class DrawList_Cumulative_Test;
+
+	public:
+
+		/** the user data behind this entry */
+		Entry entry;
+
+		/** this entry's probability */
+		double probability;
+
+	private:
+
+		/** the cumulative probability, tracked among all entries within the DrawList */
+		double cumulativeProbability;
+
+	public:
+
+		/** empty ctor */
+		DrawListEntry() :
+			entry(), probability(0), cumulativeProbability(0) {;}
+
+		/** ctor */
+		DrawListEntry(const Entry& entry, double probability) :
+			entry(entry), probability(probability), cumulativeProbability(0) {;}
+
+		/** compare this entrie's summed probability with the given probability */
+		bool operator < (double probability) const {return cumulativeProbability < probability;}
+
+	private:
+
+		/** ctor */
+		DrawListEntry(Entry& e, double probability, double summedProbability) :
+			entry(entry), probability(probability), cumulativeProbability(summedProbability) {;}
+
+
+	};
+
+
+
+	/**
+	 * a DrawList is a data-structure containing entries that have a
+	 * probability assigned to them.
+	 * using the draw() function one may draw from these entries according
+	 * to their assigned probability in O(log(n))
+	 */
+	template <typename Entry> class DrawList {
+
+	friend class DrawList_Cumulative_Test;
+
+	public:
+
+		/** ctor */
+		DrawList() : sumValid(false) {;}
+
+
+		/** append a new entry to the end of the list */
+		void push_back(const Entry& entry, const double probability) {
+			const DrawListEntry<Entry> dle(entry, probability);
+			entries.push_back(dle);
+			sumValid = false;
+		}
+
+		/** change the entry at the given position. ensure the vector is resized()! */
+		void set(const uint32_t idx, const Entry& entry, const double probability) {
+			entries[idx].entry = entry;
+			entries[idx].probability = probability;
+			sumValid = false;
+		}
+
+		/** resize the underlying vector to hold the given number of entries */
+		void resize(const uint32_t numEntries) {entries.resize(numEntries);}
+
+		/** clear all currently inserted entries */
+		void clear() {entries.clear();}
+
+		/** does the underlying vector contain any entries? */
+		bool empty() const {return entries.empty();}
+
+		/** the number of entries */
+		uint32_t size() const {return entries.size();}
+
+		/** draw a random entry from the draw list */
+		Entry& draw() {
+
+			// random value between [0, 1]
+			double rand01 = double(rand()) / double(RAND_MAX);
+
+			return draw(rand01);
+
+		}
+
+		/** draw an entry according to the given probability [0,1] */
+		Entry& draw(double rand01) {
+
+			// sanity check
+			assert(!entries.empty());
+
+			ensureCumulativeProbability();
+
+			// random value between [0, summedProbability]
+			// (this prevents us from norming the list to [0, 1])
+			double rand = rand01 * entries[entries.size()-1].cumulativeProbability;
+
+			// binary search for the matching entry O(log(n))
+			auto tmp = std::lower_bound(entries.begin(), entries.end(), rand);
+			return (*tmp).entry;
+
+//			// O(n)
+//			for (DrawListEntry<Entry>& dle : entries) {
+//				if (dle.cumulativeProbability > rand) {return dle.entry;}
+//			}
+//			return entries[this->size()-1].entry;
+
+		}
+
+	private:
+
+		/** ensure the cumulative probability is valid. if not -> calculate it */
+		void ensureCumulativeProbability() {
+
+			// already valid?
+			if (sumValid) {return;}
+
+			// calculate the cumulative probability
+			double sum = 0;
+			for (DrawListEntry<Entry>& dle : entries) {
+				sum += dle.probability;
+				dle.cumulativeProbability = sum;
+			}
+
+			// the sum is now valid
+			sumValid = true;
+
+		}
+
+	private:
+
+		/** all entries within the DrawList */
+		std::vector<DrawListEntry<Entry>> entries;
+
+		/** track wether the summedProbability is valid or not */
+		bool sumValid;
+
+	};
+
+
+
+
+}
+
+#endif // K_MATH_RANDOM_DRAWLIST_H
--- a/math/random/DrawWheel.h
+++ b/math/random/DrawWheel.h
@@ -0,0 +1,174 @@
+#ifndef K_MATH_RANDOM_DRAWLIST_H
+#define K_MATH_RANDOM_DRAWLIST_H
+
+#include <vector>
+#include "../distribution/Uniform.h"
+
+/**
+ * add entries to a list and be able to draw from them depending oh their probability
+ *
+ * souces:
+ * https://www.udacity.com/course/viewer#!/c-cs373/l-48704330/e-48748082/m-48740082
+ * https://www.youtube.com/watch?list=PLpUPoM7Rgzi_7YWn14Va2FODh7LzADBSm&feature=player_detailpage&v=kZhOJgooMxI#t=567
+ */
+namespace Random {
+
+	/**
+	 * this class represents one entry within a DrawWheel.
+	 * such an entry consists of userData denoted by the template argument
+	 * and a probability.
+	 */
+	template <typename Entry> struct DrawWheelEntry {
+
+	template <typename> friend class DrawWheel;
+
+	public:
+
+		/** the user data behind this entry */
+		Entry entry;
+
+		/** this entry's probability */
+		double probability;
+
+	public:
+
+		/** empty ctor */
+		DrawWheelEntry() : entry(), probability(0) {;}
+
+		/** ctor */
+		DrawWheelEntry(const Entry& entry, double probability) : entry(entry), probability(probability) {;}
+
+	};
+
+	/**
+	 * a DrawWheel is a data-structure containing entries that have a
+	 * probability assigned to them.
+	 * using the draw() function one may draw from these entries according
+	 * to their assigned probability in ~O(n)
+	 */
+	template <typename Entry> class DrawWheel {
+
+	private:
+
+		/** all entries within the DrawWheel */
+		std::vector<DrawWheelEntry<Entry>> entries;
+
+		/** is the current maximum valid? */
+		bool maxValid = true;
+
+		/** the maximum weight among all entries */
+		double curMax = 0;
+
+		/** the current index within the wheel */
+		int curIdx = 0;
+
+		/** the current offset at the wheel's index */
+		double curOffset = 0;
+
+		/** draw random numbers for the offset */
+		K::UniformDistribution dist;
+
+
+	public:
+
+		/** ctor */
+		DrawWheel() {;}
+
+
+		/** append a new entry to the end of the list */
+		void push_back(const Entry& entry, const double probability) {
+			entries.push_back( DrawWheelEntry<Entry>(entry, probability) );
+			if (curMax < probability) {curMax = probability;}
+		}
+
+		/** change the entry at the given position. ensure the vector is resized()! */
+		void set(const uint32_t idx, const Entry& entry, const double probability) {
+			entries[idx].entry = entry;
+			entries[idx].probability = probability;
+			maxValid = false;
+		}
+
+		/** resize the underlying vector to hold the given number of entries */
+		void resize(const uint32_t numEntries) {entries.resize(numEntries);}
+
+		/** clear all currently inserted entries */
+		void clear() {entries.clear();}
+
+		/** does the underlying vector contain any entries? */
+		bool empty() const {return entries.empty();}
+
+		/** the number of entries */
+		uint32_t size() const {return entries.size();}
+
+
+		/** call this once before drawing anything */
+		void init() {
+
+			// ensure the maximum number is correct
+			ensureMaxProbability();
+
+			// setup the distribution to draw a new offset
+			dist.reset(0, 2 * curMax);
+
+			// draw starting values
+			curIdx = K::UniformDistribution::draw( (int)0, (int)entries.size() - 1);
+			curOffset = dist.draw();
+
+		}
+
+		/** draw a random entry from the wheel */
+		Entry& draw() {
+
+			while(true) {
+
+				// found a suitable particle? use it and draw the next random number
+				if (entries[curIdx].probability >= curOffset) {
+
+					// next offset
+					curOffset += dist.draw();
+
+					// return
+					return entries[curIdx].entry;
+
+				// weight to small, subtract the elements weight and move on to the next element
+				} else {
+
+					// remove the current entries probability
+					curOffset -= entries[curIdx].probability;
+
+					// resume with the next one along the wheel
+					curIdx = (curIdx + 1) % ((int)entries.size());
+
+				}
+
+			}
+
+		}
+
+
+
+	private:
+
+		/** ensure the max probability is valid. if not -> calculate it */
+		void ensureMaxProbability() {
+
+			// valid?
+			if (maxValid) {return;}
+
+			// comparisen
+			const auto lambda = [] (const DrawWheelEntry<Entry>& e1, const DrawWheelEntry<Entry>& e2) {return e1.probability < e2.probability;};
+
+			// find the largest entry
+			const DrawWheelEntry<Entry> max = *std::max_element(entries.begin(), entries.end(), lambda);
+			this->curMax = max.probability;
+
+			// the max is now valid
+			maxValid = true;
+
+		}
+
+	};
+
+}
+
+#endif // K_MATH_RANDOM_DRAWLIST_H
--- a/math/random/Normal.h
+++ b/math/random/Normal.h
@@ -0,0 +1,27 @@
+#ifndef K_MATH_RND_NORMAL_H
+#define K_MATH_RND_NORMAL_H
+
+namespace Random {
+
+
+	/**
+	 * provides some common functions
+	 * for handling normal-distributed random numbers
+	 */
+	class Normal {
+
+	public:
+
+		/** get normal-distributed random number for given mu/sigma */
+		static double get(double mu, double sigma) {
+			static std::random_device rd;
+			static std::mt19937 gen(rd());
+			std::normal_distribution<> d(mu, sigma);
+			return d(gen);
+		}
+
+	};
+
+}
+
+#endif // K_MATH_RND_NORMAL_H
--- a/math/random/RandomGenerator.h
+++ b/math/random/RandomGenerator.h
@@ -0,0 +1,29 @@
+#ifndef RANDOM_Random::RandomGenerator_H
+#define RANDOM_Random::RandomGenerator_H
+
+#include <random>
+#include <chrono>
+
+#ifdef USE_FIXED_SEED
+	#define RANDOM_SEED 1234
+#else
+	#define RANDOM_SEED ( std::chrono::system_clock::now().time_since_epoch() / std::chrono::milliseconds(1) )
+#endif
+
+namespace Random {
+
+    class RandomGenerator : public std::minstd_rand {
+
+    public:
+
+        /** ctor with default seed */
+        RandomGenerator() : std::minstd_rand(RANDOM_SEED) {;}
+
+        /** ctor with custom seed */
+        RandomGenerator(result_type) : std::minstd_rand(RANDOM_SEED) {;}
+
+    };
+
+}
+
+#endif // K_MATH_RANDOM_Random::RandomGenerator_H
--- a/math/random/RandomIterator.h
+++ b/math/random/RandomIterator.h
@@ -0,0 +1,118 @@
+#ifndef RANDOMITERATOR_H
+#define RANDOMITERATOR_H
+
+#include <vector>
+#include <random>
+#include "../../Assertions.h"
+
+namespace Random {
+
+	template <typename Element> class RandomIterator {
+
+	private:
+
+		/** the user's data-vector to randomly iterate */
+		const std::vector<Element>& vec;
+
+		/** the number of random indices */
+		int cnt;
+
+		/** the random number generator */
+		std::minstd_rand gen;
+		bool isRandomized = false;
+
+		/** X random indices */
+		std::vector<int> indices;
+
+
+	public:
+
+		/** ctor */
+		RandomIterator(const std::vector<Element>& vec, const int cnt) : vec(vec), cnt(cnt) {
+
+			//const uint64_t ts = std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count();
+			static int seed = 0; ++seed;
+			gen.seed(seed);
+
+			// sanity check
+			if ((int)vec.size() < cnt)	{throw Exception("number of elements in list is smaller than the requested number to draw");}
+			if (cnt == 0)				{throw Exception("number of elements in list must be at least 1");}
+			if (vec.empty())			{throw Exception("empty input vector given");}
+
+			indices.resize(cnt);
+
+		}
+
+		/** create random samples (vector-indicies) that are hereafter available for iteration */
+		void randomize() {
+
+			// random-number generator between [0:size-1]
+			std::uniform_int_distribution<int> dist(0, (int) vec.size()-1);
+
+			// ensure we use each vector-index only ONCE
+			bool used[vec.size()] = {false};
+
+			// draw X random samples
+			for (int i = 0; i < cnt; ) {
+				const int rnd = dist(gen);
+				if (used[rnd]) {continue;}		// already used? try again!
+				used[rnd] = true;				// mark as used
+				indices[i] = rnd;				// add to the index
+				++i;
+			}
+
+			// the vector is setup correctly
+			isRandomized = true;
+
+		}
+
+		/** the iterator state */
+		struct Iterator {
+
+			/** the current position within "indicies" */
+			int pos;
+
+			/** the vector with the user-data to randomly iterate */
+			const std::vector<Element>& vec;
+
+			/** the vector containing the random indices */
+			const std::vector<int>& indices;
+
+			/** ctor */
+			Iterator(const int pos, const std::vector<Element>& vec, const std::vector<int>& indices) : pos(pos), vec(vec), indices(indices) {;}
+
+			/** end-of-iteration? */
+			bool operator != (const Iterator& o) const {return pos != o.pos;}
+
+			/** next value */
+			Iterator& operator ++ () {++pos; return *this;}
+
+			/** get the user-data */
+			Element operator * () {return vec[indices[pos]];}
+
+		};
+
+		//const Element& operator [] (const int idx) const {return vec[indices[idx]]; }
+
+		/** number of available random entries */
+		size_t size() const {return cnt;}
+
+
+		/** for-each access */
+		Iterator begin() const	{ ensureRandomized(); return Iterator(0, vec, indices); }
+
+		/** for-each access */
+		Iterator end() const	{ ensureRandomized(); return Iterator(cnt, vec, indices); }
+
+private:
+
+		/** ensure the coder called randomize() before using the iterator */
+		void ensureRandomized() const {
+			Assert::isTrue(isRandomized, "call randomize() before using the iterator!");
+		}
+
+	};
+
+}
+
+#endif // RANDOMITERATOR_H
--- a/math/random/Uniform.h
+++ b/math/random/Uniform.h
@@ -0,0 +1,43 @@
+#ifndef K_MATH_RANDOM_UNIFORM_H
+#define K_MATH_RANDOM_UNIFORM_H
+
+#include "RandomGenerator.h"
+
+namespace Random {
+
+	/**
+	 * provides some common functions
+	 * for handling normal-distributed random numbers
+	 */
+	template <typename T> class Uniform {
+
+	private:
+
+		Random::RandomGenerator gen;
+
+		/** depending on T, Dist is either a uniform_real or uniform_int distribution */
+		typedef typename std::conditional< std::is_floating_point<T>::value, std::uniform_real_distribution<T>, std::uniform_int_distribution<T> >::type Dist;
+		Dist dist;
+
+	public:
+
+		/** ctor */
+		Uniform(const T min, const T max) : gen(RANDOM_SEED), dist(min, max) {
+
+		}
+
+		/** get a uniformaly distributed random number */
+		T draw() {
+			return dist(gen);
+		}
+
+		/** set the seed to use */
+		void setSeed(const long seed) {
+			gen.seed(seed);
+		}
+
+	};
+
+}
+
+#endif // K_MATH_RANDOM_UNIFORM_H
--- a/math/random/Unique.h
+++ b/math/random/Unique.h
@@ -0,0 +1,26 @@
+#ifndef K_MATH_RND_UNIQUE_H
+#define K_MATH_RND_UNIQUE_H
+
+namespace Random {
+
+	/**
+	 * provides some common functions
+	 * for handling uniquely distributed random numbers
+	 */
+	class Unique {
+
+	public:
+
+		/** get uniquely distributed random number between min and max */
+		static double getBetween(double min, double max) {
+			double rnd = (double) rand() / (double) RAND_MAX;
+			rnd *= (max-min);
+			rnd += min;
+			return rnd;
+		}
+
+	};
+
+}
+
+#endif // K_MATH_RND_UNIQUE_H
--- a/math/stats/Statistics.h
+++ b/math/stats/Statistics.h
@@ -0,0 +1,169 @@
+#ifndef STATISTICS_H
+#define STATISTICS_H
+
+#include <set>
+#include <cstdint>
+#include <sstream>
+#include <cmath>
+
+namespace Stats {
+
+        /**
+         * store values here and get statistics about their
+         * avg, median, std-dev, etc.
+         */
+        template <typename T> class Statistics {
+
+        public:
+
+                /** ctor */
+                Statistics() : sum(), sumSquared(), cnt() {;}
+
+                /** dtor */
+                ~Statistics() {;}
+
+                /** add a new value */
+                void add(T value) {
+                        ++cnt;
+                        sum += value;
+                        sumSquared += (value*value);
+                        values.insert(value);
+                }
+
+                /** add all values from the given statistics instance */
+                void add(const Statistics<T>& other) {
+                        for (const T val : other.values) {
+                                this->add(val);
+                        }
+                }
+
+                /** get the std-dev of all values */
+                T getStdDev() const {
+                        double E1 = sumSquared / (double) cnt;
+                        double E2 = getAvg();
+                        return std::sqrt(E1 - (E2*E2));
+                }
+
+                /** get average value */
+                T getAvg() const {
+                        return sum / (double) cnt;
+                }
+
+                /** get the given quantile (e.g. 0.5 for median) */
+                T getQuantile(const double q) const {
+                        if (q < 0) {return *values.begin();}
+                        if (q >= 1) {return *(--values.end());}
+                        uint32_t pos = cnt * q;
+                        uint32_t curPos = 0;
+                        for (auto val : values) {
+                                if (curPos == pos) {return val;}
+                                ++curPos;
+                        }
+                        return -1;
+                }
+
+                /** get the median value (= Quantile 0.5) */
+                T getMedian() const {
+                        return getQuantile(0.5f);
+                }
+
+                /** get smallest value */
+                T getMin() const {
+                        if (values.empty()) {return -1;}
+                        return *(values.begin());
+                }
+
+                /** get largest value */
+                T getMax() const {
+                        if (values.empty()) {return -1;}
+                        return *(--values.end());
+                }
+
+                /** get the range between min an max */
+                T getRange() const {
+                        return getMax() - getMin();
+                }
+
+                /** get the squared sum */
+                T getSquaredSum() const {
+                        return sumSquared;
+                }
+
+
+                /** get the sum of all values */
+                T getSum() const {
+                        return sum;
+                }
+
+                /** get number of stored values */
+                uint32_t getCount() const {
+                        return cnt;
+                }
+
+                /** get the number of values that are below "val" */
+                T getNumValuesBelow(uint32_t val) const {
+                        uint32_t numFound = 0;
+                        for (auto curVal : values) {
+                                if (curVal > val) {return numFound;}
+                                ++numFound;
+                        }
+                        return numFound;
+                }
+
+                /** get all contained values in ascending order */
+                const std::multiset<T>& getAll() const {
+                        return values;
+                }
+
+                /** get as string */
+                std::string asString() const {
+                        std::stringstream ss;
+                        appendTo(ss);
+                        return ss.str();
+                }
+
+                /** send to the given stream */
+                void appendTo(std::ostream& out) const {
+                        out.precision(6);
+                        out.setf( std::ios::fixed, std:: ios::floatfield );
+                        out << "cnt(" << getCount() << ")\t";
+                        out << "min(" << getMin() << ")\t";
+                        out << "max(" << getMax() << ")\t";
+                        out << "range(" << getRange() << ")\t";
+                        out << "med(" << getMedian() << ")\t";
+                        out << "avg(" << getAvg() << ")\t";
+                        out << "stdDev(" << getStdDev() << ")\t";
+                }
+
+                /** send to stream */
+                inline std::ostream& operator << (std::ostream& out) const {
+                        appendTo(out); return out;
+                }
+
+                /** reset all statistics */
+                void reset() {
+                        sum = T();
+                        sumSquared = T();
+                        cnt = 0;
+                        values.clear();
+                }
+
+        private:
+
+                /** sum of all added values */
+                T sum;
+
+                /** squared sum of all added values (for std-dev) */
+                T sumSquared;
+
+                /** the number of added values */
+                uint32_t cnt;
+
+                /** multiset to sort all values */
+                std::multiset<T> values;
+
+        };
+
+}
+
+#endif // STATISTICS_H