zhukov000

HyperBall

Oct 21st, 2021
921
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1.  
  2. def calcVNum(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 = calcVNum(10**5, 3)
  22.    >>> V
  23.    4.19011
  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 calcVNumLongdouble(x, D):
  39.     """
  40.    Один шаг расчета объема гиперсферы радиусом 1 методом Монте Карло
  41.  
  42.    Функция генерирует заданное количество векторов `x` в пространстве размерности `D`
  43.    Считает количество векторов, чья длина меньше радиуса гиперсферы (меньше 1)
  44.    Возвращает отношение количества посчитанных векторов на общее количество
  45.    
  46.    При расчете используется тип данных np.longdouble
  47.  
  48.    Parameters
  49.    ----------
  50.    x: количество точек для метода Монте Карло
  51.    D: число измерений гиперсферы
  52.  
  53.    Returns
  54.    -------
  55.    V: объем гиперсферы
  56.  
  57.    Examples
  58.    --------
  59.    >>> V = calcVNumLongdouble(10**5, 3)
  60.    >>> V
  61.    4.19011
  62.    """
  63.     import numpy as np
  64.     from random import choice, uniform
  65.  
  66.     if D == 1:  return 2
  67.     genLongDouble = lambda: '0.' + ''.join(choice('0123456789') for _ in range(18))
  68.     RSp = np.longdouble(1.0)   # радиус сферы
  69.     PtsInSphere = 0
  70.  
  71.     Xi = [ [ np.longdouble(uniform(0, 1)) for j in range(D) ] for i in range(x) ]
  72.     for v in Xi:
  73.         s = np.longdouble(0.0)
  74.         for xi in v:
  75.             s += xi ** np.longdouble(2.0)
  76.         Dist = s ** np.longdouble(0.5)
  77.         if Dist < RSp:
  78.             PtsInSphere += 1
  79.  
  80.     return 2**D * PtsInSphere / x
  81.    
  82.  
  83. def showLog(log):
  84.     """
  85.    Отображает информацию о расчете методом Монте Карло в виде таблицы
  86.  
  87.    Для отображения таблицы используется модуль tabulate
  88.  
  89.    Parameters
  90.    ----------
  91.    log: вектор кортежей, содержащих значения в следующем порядке
  92.          - размерность (D)
  93.          - расчетный объем (Vol)
  94.          - количество точек (Pts)
  95.          - время выполнения расчета (Time)
  96.          - абсолютная погрешность от реального значения (Abs)
  97.          - относительная погрешность от реального значения (Rel)
  98.    """
  99.     from tabulate import tabulate
  100.    
  101.     print(tabulate(log, headers=["D", "Vol", "Pts", "Time", "Abs", "Rel"]))
  102.  
  103.  
  104. def showPlot(data):
  105.     """
  106.    Построить графики:
  107.    1) объем гиперсферы в зависимости от числа измерений
  108.    2) количество точек в методе Монте Карло в зависимости от числа измерений
  109.  
  110.    Для отображения графиков используется библиотека pylab
  111.  
  112.    Parameters
  113.    ----------
  114.    data: данные расчета методом Монте Карло - каждая строка содержит
  115.          - размерность пространства (d)
  116.          - расчитанный объем (v)
  117.          - число точек (p)
  118.    """
  119.  
  120.     import numpy as np
  121.     import pylab
  122.  
  123.     D, V, P = np.array(data).T
  124.     # График 1: расчитанный объем шара в зависимости от размерности
  125.     fig, axs = pylab.subplots(2, 1, figsize=(10,20))
  126.     axs[0].set_title('Расчитанный (x) и реальный объем (o)', color='C0')
  127.     axs[0].plot(D, V, 'C0', marker='o')
  128.     axs[0].plot(D, V, 'C2', marker='x')
  129.     # График 2: количество точек в зависимости от размерности пространства
  130.     axs[1].set_title('Количество точек', color='C1')
  131.     axs[1].plot(D, P, marker='x')
  132.  
  133.     pylab.show()
  134.  
  135.  
  136. def calculateHyperballVolume(MaxD, Eps):
  137.     """
  138.    Расчитать объем гиперсферы для нескольких размерностей с заданной точностью
  139.  
  140.    Расчет объема гиперсфер с размерностями от 1 до `MaxD` с точностью `Eps`
  141.    выполняется по методу Монте Карло
  142.    Условие окончания итерационного процесса:
  143.        abs(Vnm - Vn) / s < 1.0
  144.    где Vn - значение объема на предыдущем шаге
  145.        Vnm - усредненное значение между предыдущим и текущим шагами
  146.        s - величина ошибки, по формуле s = abs(Vnm) * Eps + Eps
  147.  
  148.    Parameters
  149.    ----------
  150.    MaxD: максимальная размерность
  151.    Eps: точность
  152.  
  153.    Returns
  154.    -------
  155.    список данных с расчетами для каждой размерности от 1 до `MaxD` в виде
  156.    [размерность, объем, число точек, время, абс.погрешность, относительная погрешность]
  157.    """
  158.     import time
  159.     from random import randint
  160.     from math import pi
  161.  
  162.     MaxPointsCnt = 10**9 # максимальное число точек
  163.  
  164.     # реальные объемы
  165.     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]
  166.     data = []
  167.  
  168.     for d in range(1, MaxD+1):          
  169.         startTime = time.time()          
  170.         n = 10**6
  171.         m = 0
  172.         V1 = Vn = calcVNum(n, d)
  173.         while n <= MaxPointsCnt:
  174.               m = randint(-10000, 10000) + 10**5
  175.               Vm = calcVNum(m, d)
  176.               Vnm = (n * Vn + m * Vm) / (n + m)
  177.  
  178.               V1 = Vnm          
  179.               dV = abs(Vnm - Vn)
  180.               s = abs(Vnm) * Eps + Eps
  181.  
  182.               if dV / s < 1.0:
  183.                 break
  184.               n += m
  185.               Vn = Vnm
  186.         else:
  187.               n -= m
  188.         endTime = time.time()
  189.  
  190.         epsA = abs(V1 - Vreal[d])          
  191.         epsR = epsA / Vreal[d]            
  192.    
  193.         data.append([d, V1, n, endTime-startTime, epsA, epsR])
  194.         print("{:2d} из {:d}".format(d, MaxD))
  195.  
  196.     return data
  197.    
  198.  
  199. data = calculateHyperballVolume(16, 10**-7)
  200. showLog(data)
  201. showPlot([row[:3] for row in data])
  202.  
RAW Paste Data