text: from __future__ import division
  import os
  import glob
  import string
  import random
  import numpy as np
  from scipy.odr import *
  import json
  import matplotlib.pyplot as plt
  import CoolProp
  import CoolProp.CoolProp as CP
  
  LIBRARY = [i / 6.0 for i in range(1, 151)] + [0.35 + i / 2000 for i in range(1, 100)] + [0.05 + 0.001 * i for i in range(1, 100)] + [i + 0.5 for i in range(10)]
  #LIBRARY = [i/1000 for i in range(1,20000)]
  
  
  class Sample(object):
      def __init__(self, v):
          self.v = v
  
  
  class GeneticAncillaryFitter(object):
      def __init__(self,
                 num_samples=600,  # Have this many chromos in the sample group
                 num_selected=60,  # Have this many chromos in the selected group
                 mutation_factor=2,  # Randomly mutate 1/n of the chromosomes
                 num_powers=5,  # How many powers in the fit
                 Ref='R407C',
                 value='rhoV',
                 addTr=True,
                 values=None,
                 Tlims=None
                  ):
          self.num_samples = num_samples
          self.num_selected = num_selected
          self.mutation_factor = mutation_factor
          self.num_powers = num_powers
          self.addTr = addTr
          self.value = value
          self.Ref = Ref
  
          # Thermodynamics
          from CoolProp.CoolProp import PropsSI
  
          if values is None:
              self.Tc = PropsSI(Ref, 'Tcrit')
              self.pc = PropsSI(Ref, 'pcrit')
              self.rhoc = PropsSI(Ref, 'rhomolar_critical')
              self.Tmin = PropsSI(Ref, 'Tmin')
              if Tlims is None:
                  self.T = np.append(np.linspace(self.Tmin + 1e-14, self.Tc - 1, 150), np.logspace(np.log10(self.Tc - 1), np.log10(self.Tc) - 1e-15, 40))
              else:
                  self.T = np.linspace(Tlims[0], Tlims[1])
              self.pL = np.array(PropsSI('P', 'T', self.T, 'Q', [0] * len(self.T), Ref))
              self.pV = np.array(PropsSI('P', 'T', self.T, 'Q', [1] * len(self.T), Ref))
              self.rhoL = PropsSI('Dmolar', 'T', self.T, 'Q', [0] * len(self.T), Ref)
              self.rhoV = PropsSI('Dmolar', 'T', self.T, 'Q', [1] * len(self.T), Ref)
          else:
              self.Tc = values['Tcrit']
              self.pc = values['pcrit']
              self.rhoc = values['rhocrit']
              self.Tmin = values['Tmin']
              self.T = np.array(values['T'])
              self.p = np.array(values['p'])
              self.pL = np.array(values['p'])
              self.pV = np.array(values['p'])
              self.rhoL = np.array(values['rhoL'])
              self.rhoV = np.array(values['rhoV'])
  
          self.logpLpc = (np.log(self.pL) - np.log(self.pc))
          self.logpVpc = (np.log(self.pV) - np.log(self.pc))
          self.rhoLrhoc = np.array(self.rhoL) / self.rhoc
          self.rhoVrhoc = np.array(self.rhoV) / self.rhoc
          self.logrhoLrhoc = np.log(self.rhoL) - np.log(self.rhoc)
          self.logrhoVrhoc = np.log(self.rhoV) - np.log(self.rhoc)
  
          self.x = 1.0 - self.T / self.Tc
  
          MM = PropsSI(Ref, 'molemass')
          self.T_r = self.Tc
  
          if self.value == 'pL':
              self.LHS = self.logpLpc.copy()
              self.EOS_value = self.pL.copy()
              if self.addTr == False:
                  self.description = "p' = pc*exp(sum(n_i*theta^t_i))"
              else:
                  self.description = "p' = pc*exp(Tc/T*sum(n_i*theta^t_i))"
              self.reducing_value = self.pc
          elif self.value == 'pV':
              self.LHS = self.logpVpc.copy()
              self.EOS_value = self.pV.copy()
              if self.addTr == False:
                  self.description = "p'' = pc*exp(sum(n_i*theta^t_i))"
              else:
                  self.description = "p'' = pc*exp(Tc/T*sum(n_i*theta^t_i))"
              self.reducing_value = self.pc
          elif self.value == 'rhoL':
              self.LHS = self.logrhoLrhoc.copy()
              self.EOS_value = self.rhoL
              if self.addTr == False:
                  self.description = "rho' = rhoc*exp(sum(n_i*theta^t_i))"
              else:
                  self.description = "rho' = rhoc*exp(Tc/T*sum(n_i*theta^t_i))"
              self.reducing_value = self.rhoc
          elif self.value == 'rhoV':
              self.LHS = self.logrhoVrhoc.copy()
              self.EOS_value = self.rhoV
              if self.addTr == False:
                  self.description = "rho'' = rhoc*exp(sum(n_i*theta^t_i))"
              else:
                  self.description = "rho'' = rhoc*exp(Tc/T*sum(n_i*theta^t_i))"
              self.reducing_value = self.rhoc
          elif self.value == 'rhoLnoexp':
              self.LHS = (self.rhoLrhoc - 1).copy()
              self.EOS_value = self.rhoL
              self.description = "rho' = rhoc*(1+sum(n_i*theta^t_i))"
              self.reducing_value = self.rhoc
          else:
              raise ValueError
  
          if self.value == 'rhoLnoexp' and self.addTr:
              raise ValueError('Invalid combination')
  
          if self.addTr:
              self.LHS *= self.T / self.Tc
  
      def generate_random_chromosomes(self,):
          '''
          Create a list of random chromosomes to seed our algorithm.
          '''
          chromos = []
          while len(chromos) < self.num_samples:
              chromos.append(Sample(sorted(random.sample(LIBRARY, self.num_powers))))
          return chromos
  
      def fitness(self, chromo):
          '''
          Fitness of a chromo is the sum of the squares of the error of the correlation
          '''
  
          # theta^t where the i-th row is equal to theta^t[i]
          # We need these terms later on to build the A and b matrices
          theta_t = (self.x.reshape(-1, 1)**chromo.v).T
  
          # TODO: more numpy broadcasting should speed this up even more
          # Another few times faster ought to be possible
          I = len(chromo.v)
          A = np.zeros((I, I))
          b = np.zeros((I, 1))
          for i in range(I):
              for j in range(I):
                  A[i, j] = np.sum(theta_t[i] * theta_t[j])
              b[i] = np.sum(theta_t[i] * self.LHS)
  
          # If you end up with a singular matrix, quit this run
          try:
              n = np.linalg.solve(A, b).T
          except np.linalg.linalg.LinAlgError as E:
              chromo.fitness = 1e99
              return
  
          chromo.beta = n
          RHS = np.sum(n * self.x.reshape(-1, 1)**chromo.v, axis=1)
  
          if self.addTr:
              RHS *= self.Tc / self.T
  
          if self.value in ['pL', 'pV']:
              fit_value = np.exp(RHS) * self.pc
          elif self.value in ['rhoL', 'rhoV']:
              fit_value = np.exp(RHS) * self.rhoc
          elif self.value == 'rhoLnoexp':
              fit_value = self.rhoc * (1 + RHS)
          else:
              raise ValueError
  
          max_abserror = np.max(np.abs((fit_value / self.EOS_value) - 1) * 100)
  
          chromo.fitness = max_abserror
          chromo.fit_value = fit_value
          chromo.max_abserror = max_abserror
  
          return chromo.fitness
  
      def tourny_select_chromo(self, samples):
          '''
          Randomly select two chromosomes from the samples, then return the one
          with the best fitness.
          '''
          a = random.choice(samples)
          b = random.choice(samples)
          if a.fitness < b.fitness:
              return a
          else:
              return b
  
      def breed(self, a, b):
          '''
          Breed two chromosomes by splicing them in a random spot and combining
          them together to form two new chromos.
          '''
          splice_pos = random.randrange(len(a.v))
          new_a = a.v[:splice_pos] + b.v[splice_pos:]
          new_b = b.v[:splice_pos] + a.v[splice_pos:]
          return Sample(sorted(new_a)), Sample(sorted(new_b))
  
      def mutate(self, chromo):
          '''
          Mutate a chromosome by changing one of the parameters, but only if it improves the fitness
          '''
          v = chromo.v
          if hasattr(chromo, 'fitness'):
              old_fitness = chromo.fitness
          else:
              old_fitness = self.fitness(chromo)
  
          for i in range(10):
              pos = random.randrange(len(chromo.v))
              chromo.v[pos] = random.choice(LIBRARY)
              new_fitness = self.fitness(chromo)
              if new_fitness < old_fitness:
                  return chromo
              else:
                  return Sample(sorted(v))
  
      def run(self):
          # Create a random sample of chromos
          samples = self.generate_random_chromosomes()
  
          # Calculate the fitness for the initial chromosomes
          for chromo in samples:
              self.fitness(chromo)
  #         print '#'
  
          decorated = sorted([(sample.fitness, sample) for sample in samples])
          samples = [s for sv, s in decorated]
          values = [sv for sv, s in decorated]
          plt.plot(values[0:len(values) // 2])
          plt.close()
  
          # Main loop: each generation select a subset of the sample and breed from
          # them.
          generation = -1
          while generation < 0 or samples[0].fitness > 0.02 or (generation < 3 and generation < 15):
              generation += 1
  
              # Generate the selected group from sample- take the top 10% of samples
              # and tourny select to generate the rest of selected.
              ten_percent = int(len(samples) * .1)
              selected = samples[:ten_percent]
              while len(selected) < self.num_selected:
                  selected.append(self.tourny_select_chromo(samples))
  
              # Generate the solution group by breeding random chromos from selected
              solution = []
              while len(solution) < self.num_samples:
                  solution.extend(self.breed(random.choice(selected),
                                           random.choice(selected)))
  
              # Apply a mutation to a subset of the solution set
              mutate_indices = random.sample(range(len(solution)), len(solution) // self.mutation_factor)
              for i in mutate_indices:
                  solution[i] = self.mutate(solution[i])
  
              for chromo in solution:
                  self.fitness(chromo)
  #             print '#'
  
              decorated = sorted([(sample.fitness, sample) for sample in solution])
              samples = [s for sv, s in decorated]
  
  #             print '------------------  Top 10 values  ---------------'
  #             for sample in samples[0:10]:
  #                 print sample.v, sample.fitness, sample.max_abserror
  
  #             print '// Max error is ',samples[0].max_abserror,'% between',np.min(self.T),'and',np.max(self.T),'K'
  #             print str(samples[0].v), samples[0].beta.tolist()
  
              # Print useful stats about this generation
              (min, median, max) = [samples[0].fitness, samples[len(samples) // 2].fitness, samples[-1].fitness]
  #             print("{0} best value: {1}. fitness: best {2}, median {3}, worst {4}".format(generation, samples[0].v, min, median, max))
  
              # If the string representations of all the chromosomes are the same, stop
              if len(set([str(s.v) for s in samples[0:5]])) == 1:
                  break
  
          self.fitness(samples[0])
  
          print(self.value)
          print('// Max error is ' + str(samples[0].max_abserror) + ' % between ' + str(np.min(self.T)) + ' and ' + str(np.max(self.T)) + ' K')
  
          self.fit_value = samples[0].fit_value
  
          j = dict()
          j['n'] = samples[0].beta.squeeze().tolist()
          j['t'] = samples[0].v
          j['Tmin'] = np.min(self.T)
          j['Tmax'] = np.max(self.T)
          j['type'] = self.value
          j['using_tau_r'] = self.addTr
          j['reducing_value'] = self.reducing_value
          j['T_r'] = self.Tc
  
          # Informational, not used
          j['max_abserror_percentage'] = samples[0].max_abserror
          j['description'] = self.description
  
          return j
  
  
  def build_ancillaries(name, **kwargs):
  
      j = dict()
      j['ANCILLARIES'] = dict()
      gaf = GeneticAncillaryFitter(Ref=name, value='pL', addTr=True, num_powers=6, **kwargs)
      j['ANCILLARIES']['pL'] = gaf.run()
      gaf = GeneticAncillaryFitter(Ref=name, value='pV', addTr=True, num_powers=6, **kwargs)
      j['ANCILLARIES']['pV'] = gaf.run()
      gaf = GeneticAncillaryFitter(Ref=name, value='rhoLnoexp', addTr=False, num_powers=6, **kwargs)
      j['ANCILLARIES']['rhoL'] = gaf.run()
      gaf = GeneticAncillaryFitter(Ref=name, value='rhoV', addTr=True, num_powers=6, **kwargs)
      j['ANCILLARIES']['rhoV'] = gaf.run()
  
      fp = open(os.path.join('ancillaries', name + '_anc.json'), 'w')
      print(json.dumps(j, indent=2), file=fp)
      fp.close()
  
  
  def build_all_ancillaries():
      for fluid in sorted(CoolProp.__fluids__):
          print(fluid)
          if fluid in ['SES36']:
              build_ancillaries(fluid, Tlims=[CP.Props(fluid, 'Ttriple'), CP.Props(fluid, 'Tcrit') - 1])
          elif fluid == 'R507A':
              build_ancillaries(fluid, Tlims=[CP.Props(fluid, 'Ttriple'), CP.Props(fluid, 'Tcrit') - 0.1])
          elif fluid == 'R407F':
              build_ancillaries(fluid, Tlims=[CP.Props(fluid, 'Ttriple'), CP.Props(fluid, 'Tcrit') - 2])
          else:
              build_ancillaries(fluid)
  
  
  if __name__ == "__main__":
  
      fluid = 'Methanol'
      RPfluid = fluid
      build_ancillaries(RPfluid, Tlims=[CP.PropsSI(fluid, 'Ttriple'), CP.PropsSI(fluid, 'Tcrit') - 0.01])
  
      #~ build_all_ancillaries()
  #     inject_ancillaries()
  