1"""Contains the classes used to generate pseudo-random numbers.
2
3Copyright (C) 2013, Joshua More and Michele Ceriotti
4
5This program is free software: you can redistribute it and/or modify
6it under the terms of the GNU General Public License as published by
7the Free Software Foundation, either version 3 of the License, or
8(at your option) any later version.
9
10This program is distributed in the hope that it will be useful,
11but WITHOUT ANY WARRANTY; without even the implied warranty of
12MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13GNU General Public License for more details.
14
15You should have received a copy of the GNU General Public License
16along with this program. If not, see <http.//www.gnu.org/licenses/>.
17
18
19Allows the user to specify a seed for the random number generator.
20These are used in initialising the velocities and in stochastic thermostats.
21The state of the random number generator is kept track of, so that the if the
22simulation is restarted from a checkpoint, we will see the same dynamics as if
23it had not been stopped.
24
25Classes:
26   Random: An interface between the numpy.random module and the user.
27"""
28
29__all__ = ['Random']
30
31import numpy as np
32import math
33
34class Random(object):
35   """Class to interface with the standard pseudo-random number generator.
36
37   Initializes the standard numpy pseudo-random number generator from a seed
38   at the beginning of the simulation, and keeps track of the state so that
39   it can be output to the checkpoint files throughout the simulation.
40
41   Attributes:
42      rng: The random number generator to be used.
43      seed: The seed number to start the generator.
44      state: A tuple of five objects giving the current state of the random
45         number generator. The first is the type of random number generator,
46         here 'MT19937', the second is an array of 624 integers, the third
47         is the current position in the array that is being read from, the
48         fourth gives whether it has a gaussian random number stored, and
49         the fifth is this stored Gaussian random number, or else the last
50         Gaussian random number returned.
51   """
52
53   def __init__(self, seed=12345, state=None):
54      """Initializes Random.
55
56      Args:
57         seed: An optional seed giving an integer to initialize the state with.
58         state: An optional state tuple to initialize the state with.
59      """
60
61      self.rng = np.random.mtrand.RandomState(seed=seed)
62      self.seed = seed
63      if state is None:
64         self.rng.seed(seed)
65      else:
66         self.state = state
67
68   def get_state(self):
69      """Interface to the standard get_state() function."""
70
71      return self.rng.get_state()
72
73   def set_state(self, value):
74      """Interface to the standard set_state() function.
75
76      Should only be used with states generated from another similar random
77      number generator, such as one from a previous run.
78      """
79
80      return self.rng.set_state(value)
81
82   state=property(get_state, set_state)
83
84   @property
85   def u(self):
86      """Interface to the standard random_sample() function.
87
88      Returns:
89         A pseudo-random number from a uniform distribution from 0-1.
90      """
91
92      return self.rng.random_sample()
93
94   @property
95   def g(self):
96      """Interface to the standard standard_normal() function.
97
98      Returns:
99         A pseudo-random number from a normal Gaussian distribution.
100      """
101
102      return self.rng.standard_normal()
103
104   def gamma(self, k, theta=1.0):
105      """Interface to the standard gamma() function.
106
107      Args:
108         k: Shape parameter for the gamma distribution.
109         theta: Mean of the distribution.
110
111      Returns:
112         A random number from a gamma distribution with a shape k and a
113         mean value theta.
114      """
115
116      return self.rng.gamma(k,theta)
117
118   def gvec(self, shape):
119      """Interface to the standard_normal array function.
120
121      Args:
122         shape: The shape of the array to be returned.
123
124      Returns:
125         An array with the required shape where each element is taken from
126         a normal Gaussian distribution.
127      """
128
129      return self.rng.standard_normal(shape)
130