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