Package csb :: Package statmech :: Module ensembles

[frames] | no frames]

# Source Code for Module csb.statmech.ensembles

      1  """
      2  Statistical Ensembles
      3  """
      4
      5  from csb.numeric import log, exp
      6  from abc import ABCMeta, abstractmethod
      7
      8
      9 -class StatisticalEnsemble(object):
     10
     11      __metaclass__ = ABCMeta
     12
     13 -    def __call__(self, raw_energies):
     14          return exp(-self.energy(raw_energies))
     15
     16 -    def log_prob(self, raw_energies):
     17          return -self.energy(raw_energies)
     18
     19      @abstractmethod
     20 -    def energy(self, raw_energies):
     21          """
     22          Transforms the raw energies as if they were observed
     23          in this statistical ensemble
     24          """
     25          pass
     26
     27 -    def gradient(self, raw_energies):
     28          raise NotImplementedError()
     29
     30
     31 -class BoltzmannEnsemble(StatisticalEnsemble):
     32
     33 -    def __init__(self, beta=1.):
     34
     35          self._beta = float(beta)
     36
     37      @property
     38 -    def beta(self):
     39          """
     40          Inverse temperature
     41          """
     42          return self._beta
     43      @beta.setter
     44 -    def beta(self, value):
     45          value = float(value)
     46          if value <= 0.:
     47              raise ValueError("Inverse temperature {0} < 0".formate(value))
     48          self._beta = value
     49
     50 -    def energy(self, raw_energies):
     51          return raw_energies * self._beta
     52
     53 -class FermiEnsemble(BoltzmannEnsemble):
     54
     55 -    def __init__(self, beta=1., e_max=0.):
     56
     57          super(FermiEnsemble, self).__init__(beta)
     58          self._e_max = float(e_max)
     59
     60      @property
     61 -    def e_max(self):
     62          """
     63          Maximum energy
     64          """
     65          return self._e_max
     66      @e_max.setter
     67 -    def e_max(self, value):
     68          self._e_max = float(value)
     69
     70 -    def energy(self, raw_energies):
     71
     72          from numpy import isinf
     73
     74          if isinf(self.beta):
     75              m = (raw_energies >= self.e_max).astype('f')
     76              return - m * log(0.)
     77          else:
     78              x = 1 + exp(self.beta * (raw_energies - self.e_max))
     79              return log(x)
     80
     81 -class TsallisEnsemble(StatisticalEnsemble):
     82
     83 -    def __init__(self, q=1., e_min=0.):
     84
     85          self._q = q
     86          self._e_min = e_min
     87
     88      @property
     89 -    def q(self):
     90          """
     91          q-analoge of the temperature
     92          """
     93          return self._q
     94      @q.setter
     95 -    def q(self, value):
     96          if value <= 0.:
     97              raise ValueError("Inverse temperature {0} < 0".formate(value))
     98          self._q = value
     99
    100      @property
    101 -    def e_min(self):
    102          """
    103          lower bound of the energy
    104          """
    105          return self._e_min
    106      @e_min.setter
    107 -    def e_min(self, value):
    108          self._e_min = value
    109
    110 -    def energy(self, raw_energies):
    111          q = self.q
    112          e_min = self.e_min
    113
    114          if (q < 1 + 1e-10):
    115              return raw_energies * q
    116          else:
    117              return log(1 + (raw_energies - e_min) * (q - 1)) * q / (q - 1) + e_min
    118
    119
    120 -class CompositeEnsemble(StatisticalEnsemble):
    121
    122 -    def __init__(self, ensembles=[]):
    123
    124          self._ensembles = ensembles
    125
    126      @property
    127 -    def ensembles(self):
    128          """
    129          Collection of statistical ensembles
    130          """
    131          return self._ensembles
    132      @ensembles.setter
    133 -    def ensembles(self, value):
    134          if not isinstance(value, list):
    135              if len(value) > 0:
    136                  if not isinstance(value[0], StatisticalEnsemble):
    137                      raise  ValueError("Not a list of statistical ensembles")
    138                  else:
    139                      self._enesmbles = value
    140              else:
    141                  self._enesmbles = value
    142
    143 -    def energy(self, raw_energies):
    144          return sum([self._ensembles[i].energy(raw_energies[i])
    145                      for i in range(len(self.ensembles))], 0)
    146
    147 -    def gradient(self, raw_energies):
    148          return sum([self._ensembles[i].gradient(raw_energies[i])
    149                      for i in range(len(self.ensembles))], 0)
    150