HyperBall

def calcVNumPy(x, D):
    """
    Один шаг расчета объема гиперсферы радиусом 1 методом Монте Карло

    Функция генерирует заданное количество векторов `x` в пространстве размерности `D`
    Считает количество векторов, чья длина меньше радиуса гиперсферы (меньше 1)
    Возвращает отношение количества посчитанных векторов на общее количество

    Parameters
    ----------
    x: количество точек для метода Монте Карло
    D: число измерений гиперсферы

    Returns
    -------
    V: объем гиперсферы

    Examples
    --------
    >>> V = calcVNumPy(10**5, 3)
    >>> V

    """
    import numpy as np

    RSp = 1   # радиус сферы
    Xi = np.random.uniform(0, 1, x * D).reshape(x, D)
    if D > 1:
      Dist = np.sum(Xi**2, axis=1)**0.5
      PtsInSphere = np.sum(Dist < RSp)
    elif D == 1:
      PtsInSphere = x

    return 2**D * PtsInSphere / x


def showLog(log):
    """
    Отображает информацию о расчете методом Монте Карло в виде таблицы

    Для отображения таблицы используется модуль tabulate

    Parameters
    ----------
    log: вектор кортежей, содержащих значения в следующем порядке
          - размерность (D)
          - расчетный объем (Vol)
          - количество точек (Pts)
          - время выполнения расчета (Time)
          - абсолютная погрешность от реального значения (Abs)
          - относительная погрешность от реального значения (Rel)
    """
    from tabulate import tabulate

    print(tabulate(log, headers=["D", "Vol", "Pts", "Time", "Abs", "Rel"]))


def showPlot(data):
    """
    Построить графики:
    1) объем гиперсферы в зависимости от числа измерений
    2) количество точек в методе Монте Карло в зависимости от числа измерений

    Для отображения графиков используется библиотека pylab

    Parameters
    ----------
    data: данные расчета методом Монте Карло - каждая строка содержит
          - размерность пространства (d)
          - расчитанный объем (v)
          - число точек (p)
    """

    import numpy as np
    import pylab

    D, V, P = np.array(data).T
    # График 1: расчитанный объем шара в зависимости от размерности
    fig, axs = pylab.subplots(2, 1, figsize=(10,20))
    axs[0].set_title('Расчитанный (x) и реальный объем (o)', color='C0')
    axs[0].plot(D, V, 'C0', marker='o')
    axs[0].plot(D, V, 'C2', marker='x')
    # График 2: количество точек в зависимости от размерности пространства
    axs[1].set_title('Количество точек', color='C1')
    axs[1].plot(D, P, marker='x')

    pylab.show()


def calculateHyperballVolume(MaxD, Eps):
    """
    Расчитать объем гиперсферы для нескольких размерностей с заданной точностью

    Расчет объема гиперсфер с размерностями от 1 до `MaxD` с точностью `Eps`
    выполняется по методу Монте Карло
    Условие окончания итерационного процесса:
        abs(Vnm - Vn) / s < 1.0
    где Vn - значение объема на предыдущем шаге
        Vnm - усредненное значение между предыдущим и текущим шагами
        s - величина ошибки, по формуле s = abs(Vnm) * Eps + Eps

    Parameters
    ----------
    MaxD: максимальная размерность
    Eps: точность

    Returns
    -------
    список данных с расчетами для каждой размерности от 1 до `MaxD` в виде
    [размерность, объем, число точек, время, абс.погрешность, относительная погрешность]
    """
    import time
    from math import pi

    # реальные объемы
    Vreal = [1, 2, pi, 4*pi/3, pi**2/2, 8*pi**2/15, pi**3/6, 16*pi**3/105, pi**4/24, 32*pi**4/945, pi**5/120, 64*pi**5/10395, pi**6/720, 128*pi**6/135135, pi**7/5040, 15*pi**7/118771, 0.23533063035889312]

    data = []

    for d in range(1, MaxD+1):
        startTime = time.time()
        m = n = 10**5
        Vn = calcVNumPy(n, d)
        while True:
              Vm = calcVNumPy(m, d)
              Vnm = (n * Vn + m * Vm) / (n + m)

              V1 = Vnm
              dV = abs(Vnm - Vn)
              s = abs(Vnm) * Eps + Eps

              if dV / s < 1.0:
                break
              n += m
              Vn = Vnm
        endTime = time.time()

        epsA = abs(V1 - Vreal[d])
        epsR = epsA / Vreal[d]

        data.append([d, V1, n, endTime-startTime, epsA, epsR])
        print("{:2d} из {:d}".format(d, MaxD))

    return data


data = calculateHyperballVolume(15, 10**-6)
showLog(data)
showPlot([row[:3] for row in data])