added von mises distributionb
quick&dirty: added activity to the grid-walkers
This commit is contained in:
57
math/distribution/Bessel.h
Normal file
57
math/distribution/Bessel.h
Normal file
@@ -0,0 +1,57 @@
|
||||
#ifndef BESSEL_H
|
||||
#define BESSEL_H
|
||||
|
||||
namespace Distribution {
|
||||
|
||||
/** helper class */
|
||||
template <typename T> class Bessel {
|
||||
|
||||
public:
|
||||
|
||||
/**
|
||||
* calculation for I_v(z)
|
||||
* https://en.wikipedia.org/wiki/Bessel_function
|
||||
* http://www.boost.org/doc/libs/1_35_0/libs/math/doc/sf_and_dist/html/math_toolkit/special/bessel/mbessel.html
|
||||
*/
|
||||
T getModified(const int v, const T z) const {
|
||||
|
||||
// running factorials (updated within for-loops)
|
||||
uint64_t runFac1 = 1;
|
||||
uint64_t runFac2 = factorial(v); // Delta(k+v+1) = (k+v+1-1)! = (k+v)! [k starts at 0] -> (v)!
|
||||
|
||||
// running potential
|
||||
T runPot = 1;
|
||||
const T pot = T(0.25) * z * z;
|
||||
|
||||
// start for-loop at k=1 instead of k=0. allows for running factorial without using if (k==0)
|
||||
T sum = 0;
|
||||
for (int k = 0; k < 16; ) {
|
||||
sum += runPot / runFac1 / runFac2;
|
||||
++k;
|
||||
runFac1 *= k;
|
||||
runFac2 *= k;
|
||||
runPot *= pot;
|
||||
}
|
||||
|
||||
// done
|
||||
return std::pow( (T(0.5) * z), v) * sum;
|
||||
|
||||
}
|
||||
|
||||
private:
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
static uint64_t factorial(const int n) {
|
||||
uint64_t res = 1;
|
||||
for (int i = 2; i <= n; ++i) {res *= i;}
|
||||
return res;
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif // BESSEL_H
|
||||
58
math/distribution/LUT.h
Normal file
58
math/distribution/LUT.h
Normal file
@@ -0,0 +1,58 @@
|
||||
#ifndef LUT_H
|
||||
#define LUT_H
|
||||
|
||||
#include <vector>
|
||||
#include "../Math.h"
|
||||
|
||||
namespace Distribution {
|
||||
|
||||
template <typename Scalar> class LUT {
|
||||
|
||||
private:
|
||||
|
||||
Scalar min;
|
||||
Scalar max;
|
||||
Scalar diff;
|
||||
|
||||
/** number of samples between min and max */
|
||||
int numSamples;
|
||||
|
||||
/** the look-up-table */
|
||||
std::vector<Scalar> samples;
|
||||
|
||||
|
||||
public:
|
||||
|
||||
/** ctor */
|
||||
template <typename Distribution> LUT(const Distribution dist, const Scalar min, const Scalar max, const int numSamples) :
|
||||
min(min), max(max), diff(max-min), numSamples(numSamples) {
|
||||
|
||||
buildLUT(dist);
|
||||
|
||||
}
|
||||
|
||||
/** get the probability of val from the LUT */
|
||||
Scalar getProbability(const Scalar val) const {
|
||||
int idx = ((val - min) * numSamples / diff);
|
||||
idx = Math::clamp(0, numSamples-1, idx);
|
||||
return samples[idx];
|
||||
}
|
||||
|
||||
private:
|
||||
|
||||
/** build the look-up-table */
|
||||
template <typename Distribution> void buildLUT(const Distribution dist) {
|
||||
|
||||
samples.resize(numSamples);
|
||||
for (int idx = 0; idx < numSamples; ++idx) {
|
||||
const Scalar val = min + (idx * diff / numSamples);
|
||||
samples[idx] = dist.getProbability(val);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif // LUT_H
|
||||
56
math/distribution/VonMises.h
Normal file
56
math/distribution/VonMises.h
Normal file
@@ -0,0 +1,56 @@
|
||||
#ifndef K_MATH_DISTRIBUTION_VONMISES_H
|
||||
#define K_MATH_DISTRIBUTION_VONMISES_H
|
||||
|
||||
#include <cmath>
|
||||
#include <random>
|
||||
|
||||
#include "../Math.h"
|
||||
#include "Bessel.h"
|
||||
#include "LUT.h"
|
||||
|
||||
namespace Distribution {
|
||||
|
||||
/** von-mises distribution */
|
||||
template <typename T> class VonMises {
|
||||
|
||||
private:
|
||||
|
||||
/** center of the distribution */
|
||||
const T mu;
|
||||
|
||||
/** like 1.0/variance of the distribution */
|
||||
const T kappa;
|
||||
|
||||
/** pre-calcuated look-up-table */
|
||||
std::vector<T> lut;
|
||||
|
||||
/** helper for modified bessel I_0(kappa) */
|
||||
Bessel<T> bessel;
|
||||
|
||||
public:
|
||||
|
||||
/** ctor */
|
||||
VonMises(const T mu, const T kappa) : mu(mu), kappa(kappa) {
|
||||
|
||||
}
|
||||
|
||||
LUT<T> getLUT() const {
|
||||
return LUT<T>(*this, -M_PI, +M_PI, 10000);
|
||||
}
|
||||
|
||||
/** get probability for the given value */
|
||||
T getProbability(const T _val) const {
|
||||
|
||||
const T val = Math::clamp((T)-M_PI, (T)+M_PI, _val);
|
||||
|
||||
return
|
||||
std::exp(kappa * std::cos(val - mu)) /
|
||||
(2 * M_PI * bessel.getModified(0, kappa));
|
||||
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif // K_MATH_DISTRIBUTION_VONMISES_H
|
||||
Reference in New Issue
Block a user