quasirandom.py 7.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179
  1. import torch
  2. from typing import Optional
  3. class SobolEngine(object):
  4. r"""
  5. The :class:`torch.quasirandom.SobolEngine` is an engine for generating
  6. (scrambled) Sobol sequences. Sobol sequences are an example of low
  7. discrepancy quasi-random sequences.
  8. This implementation of an engine for Sobol sequences is capable of
  9. sampling sequences up to a maximum dimension of 21201. It uses direction
  10. numbers from https://web.maths.unsw.edu.au/~fkuo/sobol/ obtained using the
  11. search criterion D(6) up to the dimension 21201. This is the recommended
  12. choice by the authors.
  13. References:
  14. - Art B. Owen. Scrambling Sobol and Niederreiter-Xing points.
  15. Journal of Complexity, 14(4):466-489, December 1998.
  16. - I. M. Sobol. The distribution of points in a cube and the accurate
  17. evaluation of integrals.
  18. Zh. Vychisl. Mat. i Mat. Phys., 7:784-802, 1967.
  19. Args:
  20. dimension (Int): The dimensionality of the sequence to be drawn
  21. scramble (bool, optional): Setting this to ``True`` will produce
  22. scrambled Sobol sequences. Scrambling is
  23. capable of producing better Sobol
  24. sequences. Default: ``False``.
  25. seed (Int, optional): This is the seed for the scrambling. The seed
  26. of the random number generator is set to this,
  27. if specified. Otherwise, it uses a random seed.
  28. Default: ``None``
  29. Examples::
  30. >>> soboleng = torch.quasirandom.SobolEngine(dimension=5)
  31. >>> soboleng.draw(3)
  32. tensor([[0.0000, 0.0000, 0.0000, 0.0000, 0.0000],
  33. [0.5000, 0.5000, 0.5000, 0.5000, 0.5000],
  34. [0.7500, 0.2500, 0.2500, 0.2500, 0.7500]])
  35. """
  36. MAXBIT = 30
  37. MAXDIM = 21201
  38. def __init__(self, dimension, scramble=False, seed=None):
  39. if dimension > self.MAXDIM or dimension < 1:
  40. raise ValueError("Supported range of dimensionality "
  41. f"for SobolEngine is [1, {self.MAXDIM}]")
  42. self.seed = seed
  43. self.scramble = scramble
  44. self.dimension = dimension
  45. cpu = torch.device("cpu")
  46. self.sobolstate = torch.zeros(dimension, self.MAXBIT, device=cpu, dtype=torch.long)
  47. torch._sobol_engine_initialize_state_(self.sobolstate, self.dimension)
  48. if not self.scramble:
  49. self.shift = torch.zeros(self.dimension, device=cpu, dtype=torch.long)
  50. else:
  51. self._scramble()
  52. self.quasi = self.shift.clone(memory_format=torch.contiguous_format)
  53. self._first_point = (self.quasi / 2 ** self.MAXBIT).reshape(1, -1)
  54. self.num_generated = 0
  55. def draw(self, n: int = 1, out: Optional[torch.Tensor] = None,
  56. dtype: torch.dtype = torch.float32) -> torch.Tensor:
  57. r"""
  58. Function to draw a sequence of :attr:`n` points from a Sobol sequence.
  59. Note that the samples are dependent on the previous samples. The size
  60. of the result is :math:`(n, dimension)`.
  61. Args:
  62. n (Int, optional): The length of sequence of points to draw.
  63. Default: 1
  64. out (Tensor, optional): The output tensor
  65. dtype (:class:`torch.dtype`, optional): the desired data type of the
  66. returned tensor.
  67. Default: ``torch.float32``
  68. """
  69. if self.num_generated == 0:
  70. if n == 1:
  71. result = self._first_point.to(dtype)
  72. else:
  73. result, self.quasi = torch._sobol_engine_draw(
  74. self.quasi, n - 1, self.sobolstate, self.dimension, self.num_generated, dtype=dtype,
  75. )
  76. result = torch.cat((self._first_point, result), dim=-2)
  77. else:
  78. result, self.quasi = torch._sobol_engine_draw(
  79. self.quasi, n, self.sobolstate, self.dimension, self.num_generated - 1, dtype=dtype,
  80. )
  81. self.num_generated += n
  82. if out is not None:
  83. out.resize_as_(result).copy_(result)
  84. return out
  85. return result
  86. def draw_base2(self, m: int, out: Optional[torch.Tensor] = None,
  87. dtype: torch.dtype = torch.float32) -> torch.Tensor:
  88. r"""
  89. Function to draw a sequence of :attr:`2**m` points from a Sobol sequence.
  90. Note that the samples are dependent on the previous samples. The size
  91. of the result is :math:`(2**m, dimension)`.
  92. Args:
  93. m (Int): The (base2) exponent of the number of points to draw.
  94. out (Tensor, optional): The output tensor
  95. dtype (:class:`torch.dtype`, optional): the desired data type of the
  96. returned tensor.
  97. Default: ``torch.float32``
  98. """
  99. n = 2 ** m
  100. total_n = self.num_generated + n
  101. if not (total_n & (total_n - 1) == 0):
  102. raise ValueError("The balance properties of Sobol' points require "
  103. "n to be a power of 2. {0} points have been "
  104. "previously generated, then: n={0}+2**{1}={2}. "
  105. "If you still want to do this, please use "
  106. "'SobolEngine.draw()' instead."
  107. .format(self.num_generated, m, total_n))
  108. return self.draw(n=n, out=out, dtype=dtype)
  109. def reset(self):
  110. r"""
  111. Function to reset the ``SobolEngine`` to base state.
  112. """
  113. self.quasi.copy_(self.shift)
  114. self.num_generated = 0
  115. return self
  116. def fast_forward(self, n):
  117. r"""
  118. Function to fast-forward the state of the ``SobolEngine`` by
  119. :attr:`n` steps. This is equivalent to drawing :attr:`n` samples
  120. without using the samples.
  121. Args:
  122. n (Int): The number of steps to fast-forward by.
  123. """
  124. if self.num_generated == 0:
  125. torch._sobol_engine_ff_(self.quasi, n - 1, self.sobolstate, self.dimension, self.num_generated)
  126. else:
  127. torch._sobol_engine_ff_(self.quasi, n, self.sobolstate, self.dimension, self.num_generated - 1)
  128. self.num_generated += n
  129. return self
  130. def _scramble(self):
  131. g: Optional[torch.Generator] = None
  132. if self.seed is not None:
  133. g = torch.Generator()
  134. g.manual_seed(self.seed)
  135. cpu = torch.device("cpu")
  136. # Generate shift vector
  137. shift_ints = torch.randint(2, (self.dimension, self.MAXBIT), device=cpu, generator=g)
  138. self.shift = torch.mv(shift_ints, torch.pow(2, torch.arange(0, self.MAXBIT, device=cpu)))
  139. # Generate lower triangular matrices (stacked across dimensions)
  140. ltm_dims = (self.dimension, self.MAXBIT, self.MAXBIT)
  141. ltm = torch.randint(2, ltm_dims, device=cpu, generator=g).tril()
  142. torch._sobol_engine_scramble_(self.sobolstate, ltm, self.dimension)
  143. def __repr__(self):
  144. fmt_string = [f'dimension={self.dimension}']
  145. if self.scramble:
  146. fmt_string += ['scramble=True']
  147. if self.seed is not None:
  148. fmt_string += [f'seed={self.seed}']
  149. return self.__class__.__name__ + '(' + ', '.join(fmt_string) + ')'