# 如果下面出现代码的话, 将会以这个函数代指随机整数的生成方法
# 调用一次将会返回一个随机整数
def random_integer() -> int: ...

def random() -> float:
	return float(random_integer() - MIN) / (MAX - MIN)

import numpy as np
from numpy.random import random     # uniform random number in range [0, 1)
from matplotlib import pyplot as plt


def erfinv(x: float) -> float:
    # Mike Giles - Approximating the erfinv function
    w = -np.log(1. - x*x)
    if w < 5.:
        w -= 2.5
        p = (((((((2.81022636e-08   * w + \
                   3.43273939e-07)  * w + \
                   -3.5233877e-06)  * w + \
                   -4.39150654e-06) * w + \
                   0.00021858087)   * w + \
                   -0.00125372503)  * w + \
                   -0.00417768164)  * w + \
                   0.246640727)     * w + \
                   1.50140941
    else:
        w = np.sqrt(w) - 3.
        p = (((((((-0.000200214257 * w + \
                   0.000100950558) * w + \
                   0.00134934322)  * w + \
                   -0.00367342844) * w + \
                   0.00573950773)  * w + \
                   -0.0076224613)  * w + \
                   0.00943887047)  * w + \
                   1.00167406)     * w + \
                   2.83297682
    return p * x


class RandomDistribution:       # basic interface

    def pdf(self, x: float) -> float:   # probability density function
        raise NotImplementedError

    def __call__(self) -> float:        # ganerate a random number with a distribution
        raise NotImplementedError


class Triangular(RandomDistribution):
    def __init__(self, _min: float, _max: float, _mode: float) -> None:
        assert _min < _max and _min <= _mode <= _max
        self.min = _min
        self.max = _max
        self.mode = _mode
        # precompute some values
        self._call_coff1 = (_mode - _min) / (_max - _min)
        self._call_coff2 = (_max - _min) * (_mode - _min)
        self._call_coff3 = (_max - _min) * (_max - _mode)
        self._pdf_coff1 = 2. / self._call_coff2
        self._pdf_coff2 = 2. / self._call_coff3

    def pdf(self, x: float) -> float:
        if self.min <= x < self.mode:
            return self._pdf_coff1 * (x - self.min)
        elif self.mode <= x < self.max:
            return self._pdf_coff2 * (self.max - x)
        else:
            return 0

    def __call__(self) -> float:
        x = random()
        if x < self._call_coff1:
            return self.min + np.sqrt(self._call_coff2 * x)
        else:
            return self.max - np.sqrt(self._call_coff3 * (1. - x))


class Normal(RandomDistribution):
    def __init__(self, sigma: float, mu: float) -> None:
        assert sigma > 0.
        self.sigma = sigma
        self.mu = mu
        # precompute some values
        self._call_coff1 = sigma * np.sqrt(2.)
        self._pdf_coff1 = 1. / (self._call_coff1 * np.sqrt(np.pi))
        self._pdf_coff2 = 1. / (2 * sigma*sigma)

    def pdf(self, x: float) -> float:
        return self._pdf_coff1 * np.exp(-np.square(x - self.mu) * self._pdf_coff2)

    def __call__(self) -> float:
        return self._call_coff1 * erfinv(2. * random() - 1.) + self.mu


def test_random_distribution(generator: RandomDistribution, _min: float, _max: float, interval: float, n_samples: int) -> None:
    n_counting = int((_max - _min) / interval)
    _coff1 = n_counting / (_max - _min)
    counting = np.zeros([n_counting], float)
    for _ in np.arange(n_samples):
        y = generator()
        index = int(_coff1 * (y - _min))
        if 0 <= index < n_counting:
            counting[index] += 1.
    counting *= 1. / (n_samples * interval)
    x_counting = np.linspace(_min, _max, n_counting)
    pdf = np.array([generator.pdf(x) for x in x_counting])
    plt.bar(x_counting, counting, width=interval)
    plt.plot(x_counting, pdf, color="#FBAF22")
    plt.show()


#####################################################################################


test_random_distribution(Triangular(-2, 1, -1), -2.5, 2.5, 0.01, 1000000)

test_random_distribution(Normal(1., 0.), -4, 4, 0.01, 1000000)