zhukov000

HyperBall

Oct 21st, 2021
765
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1.  
  2. def calcVNumPy(x, D):
  3.     """
  4.    Один шаг расчета объема гиперсферы радиусом 1 методом Монте Карло
  5.  
  6.    Функция генерирует заданное количество векторов `x` в пространстве размерности `D`
  7.    Считает количество векторов, чья длина меньше радиуса гиперсферы (меньше 1)
  8.    Возвращает отношение количества посчитанных векторов на общее количество
  9.  
  10.    Parameters
  11.    ----------
  12.    x: количество точек для метода Монте Карло
  13.    D: число измерений гиперсферы
  14.  
  15.    Returns
  16.    -------
  17.    V: объем гиперсферы
  18.  
  19.    Examples
  20.    --------
  21.    >>> V = calcVNumPy(10**5, 3)
  22.    >>> V
  23.  
  24.    """
  25.     import numpy as np
  26.    
  27.     RSp = 1   # радиус сферы
  28.     Xi = np.random.uniform(0, 1, x * D).reshape(x, D)    
  29.     if D > 1:
  30.       Dist = np.sum(Xi**2, axis=1)**0.5                    
  31.       PtsInSphere = np.sum(Dist < RSp)                    
  32.     elif D == 1:
  33.       PtsInSphere = x
  34.  
  35.     return 2**D * PtsInSphere / x                        
  36.  
  37.  
  38. def showLog(log):
  39.     """
  40.    Отображает информацию о расчете методом Монте Карло в виде таблицы
  41.  
  42.    Для отображения таблицы используется модуль tabulate
  43.  
  44.    Parameters
  45.    ----------
  46.    log: вектор кортежей, содержащих значения в следующем порядке
  47.          - размерность (D)
  48.          - расчетный объем (Vol)
  49.          - количество точек (Pts)
  50.          - время выполнения расчета (Time)
  51.          - абсолютная погрешность от реального значения (Abs)
  52.          - относительная погрешность от реального значения (Rel)
  53.    """
  54.     from tabulate import tabulate
  55.    
  56.     print(tabulate(log, headers=["D", "Vol", "Pts", "Time", "Abs", "Rel"]))
  57.  
  58.  
  59. def showPlot(data):
  60.     """
  61.    Построить графики:
  62.    1) объем гиперсферы в зависимости от числа измерений
  63.    2) количество точек в методе Монте Карло в зависимости от числа измерений
  64.  
  65.    Для отображения графиков используется библиотека pylab
  66.  
  67.    Parameters
  68.    ----------
  69.    data: данные расчета методом Монте Карло - каждая строка содержит
  70.          - размерность пространства (d)
  71.          - расчитанный объем (v)
  72.          - число точек (p)
  73.    """
  74.  
  75.     import numpy as np
  76.     import pylab
  77.  
  78.     D, V, P = np.array(data).T
  79.     # График 1: расчитанный объем шара в зависимости от размерности
  80.     fig, axs = pylab.subplots(2, 1, figsize=(10,20))
  81.     axs[0].set_title('Расчитанный (x) и реальный объем (o)', color='C0')
  82.     axs[0].plot(D, V, 'C0', marker='o')
  83.     axs[0].plot(D, V, 'C2', marker='x')
  84.     # График 2: количество точек в зависимости от размерности пространства
  85.     axs[1].set_title('Количество точек', color='C1')
  86.     axs[1].plot(D, P, marker='x')
  87.  
  88.     pylab.show()
  89.  
  90.  
  91. def calculateHyperballVolume(MaxD, Eps):
  92.     """
  93.    Расчитать объем гиперсферы для нескольких размерностей с заданной точностью
  94.  
  95.    Расчет объема гиперсфер с размерностями от 1 до `MaxD` с точностью `Eps`
  96.    выполняется по методу Монте Карло
  97.    Условие окончания итерационного процесса:
  98.        abs(Vnm - Vn) / s < 1.0
  99.    где Vn - значение объема на предыдущем шаге
  100.        Vnm - усредненное значение между предыдущим и текущим шагами
  101.        s - величина ошибки, по формуле s = abs(Vnm) * Eps + Eps
  102.  
  103.    Parameters
  104.    ----------
  105.    MaxD: максимальная размерность
  106.    Eps: точность
  107.  
  108.    Returns
  109.    -------
  110.    список данных с расчетами для каждой размерности от 1 до `MaxD` в виде
  111.    [размерность, объем, число точек, время, абс.погрешность, относительная погрешность]
  112.    """
  113.     import time
  114.     from math import pi
  115.  
  116.     # реальные объемы
  117.     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]
  118.    
  119.     data = []
  120.  
  121.     for d in range(1, MaxD+1):          
  122.         startTime = time.time()          
  123.         m = n = 10**5
  124.         Vn = calcVNumPy(n, d)
  125.         while True:
  126.               Vm = calcVNumPy(m, d)
  127.               Vnm = (n * Vn + m * Vm) / (n + m)
  128.  
  129.               V1 = Vnm          
  130.               dV = abs(Vnm - Vn)
  131.               s = abs(Vnm) * Eps + Eps
  132.  
  133.               if dV / s < 1.0:
  134.                 break
  135.               n += m
  136.               Vn = Vnm
  137.         endTime = time.time()
  138.  
  139.         epsA = abs(V1 - Vreal[d])          
  140.         epsR = epsA / Vreal[d]            
  141.    
  142.         data.append([d, V1, n, endTime-startTime, epsA, epsR])
  143.         print("{:2d} из {:d}".format(d, MaxD))
  144.  
  145.     return data
  146.    
  147.  
  148. data = calculateHyperballVolume(15, 10**-6)
  149. showLog(data)
  150. showPlot([row[:3] for row in data])
RAW Paste Data