class Config(object):
   """
      A configuration object for an optimization algorithm. Each optimizer is created with a configuration object that is used to parameterize the optimizer instance. A configuration object may have arbitrary properties and methods. Default versions of methods used by several optimizers are provided.
      
      The initializer allows arbitrary keywords so that config objects
      can be quickly created, e.g.
      
      from pyec.config import Config as _
      opt = SomeOptimizer[_(p=.5,q=.2)]
   """
   def __init__(self, **kwargs):
      self.__dict__['__properties__'] = {}
      for k,v in kwargs.iteritems():
         setattr(self, k, v)
   
   def __getattr__(self, key):
      if not self.__dict__.has_key(key):
         return None
      return self.__dict__[key]

   def __setattr__(self, key, val):
      self.__dict__[key] = val
      self.__dict__['__properties__'][key] = val
      
   def merge(self, cfg):
      """
         Given a second config object, produce a new config object
         that contains the properties of this config and the argument,
         with the argument given precedence.
         
         The copy is shallow, so all objects will share any internal
         state.
         
         :param cfg: A second configuration
         :type cfg: :class:`Config`
         :returns: A third :class:`Config` object with the desired properties
      """
      ret = self.__class__(**self.__properties__)
      for k,v in cfg.__properties__.iteritems():
          setattr(ret, k, v)
      return ret class Convolution(PopulationDistribution, HistoryMapper):
   """Build a convolved optimizer.
   
   :param subs: The optimizer instances to convolve
   :type subs: A list of :class:`PopulationDistribution` objects
   
   """
   config = Config()
   
   def __init__(self, subs, **kwargs):
      if not len(subs):
          raise ValueError("Convolution requires at least one optimizer")
          
      kwargs["history"] = self.makeHistory(subs)
      super(Convolution, self).__init__(**kwargs)
      self.subs = subs
      self.history = None
      
      # force any self-convolutions to self checkpoint
      # since they are now inside of a convolution
      for sub in subs:
          if hasattr(sub, 'checkpoint') and not sub.checkpoint:
              sub.checkpoint = True
      
   def makeHistory(self, subs):
      def generator(config):
          hs = [opt.config.history for opt in subs]
          return CheckpointedMultipleHistory(config, *hs)
      return generator

   def compatible(self, history):
      """A convolution is compatible with a CheckpointedHistory
      only if it contains a history that is compatible for each
      suboptimizer.
      
      """
      return (isinstance(history, CheckpointedMultipleHistory) and
          np.array([self.mapHistory(s, history) is not None 
                    for s in self.subs]).all()) 

   def batch(self, popSize):
       pop = None
       self.history.checkpoint()
       for sub in self.subs:
           if pop is not None:
               fitness = sub.needsScores() and self.fitness or None
               self.history.update(pop, fitness, sub.config.space, sub)
           sub.update(self.mapHistory(sub), self.fitness)
           pop = sub()
             
       self.history.rollback()
       return pop
   
   def update(self, history, fitness):
      super(Convolution, self).update(history, fitness)
      for sub in self.subs:
          sub.update(self.mapHistory(sub), fitness)  
      
   def needsScores(self):
      return self.subs[0].needsScores() class StructureProposal(StructureSearch,
                        ProposalDistribution,
                        PopulationDistribution):
   """A proposal distribution for structure searching with a heuristic
   algorithm such as simulated annealing.
   
   Config params:
   
   * branchFactor The maximum number of parents per variable
   * remove_edge_prob A callable object that takes a network and returns a
                      probability of removing a random edge.
   * add_edge_prob A callable object that takes a network and returns a
                   probability of adding a random edge
   * reverse_edge_prob A callable object that takes a network and
                       returns a probability of reversing a random edges
   
   """
   config = Config(
      branchFactor = 5,
      remove_edge_prob = ProbDecayInvEdgeRatio(), 
      add_edge_prob = ProbDecayEdgeRatio(), 
      reverse_edge_prob = ProbDecayEdgeRatio(),
      history = MarkovHistory,
   )
   
   def __init__(self, **kwargs): 
      PopulationDistribution.__init__(self,**kwargs)
      self.prem = self.config.remove_edge_prob
      self.padd = self.config.add_edge_prob
      self.prev = self.config.reverse_edge_prob
      self.network = None
      self.data = None
      self.history = None

   def compatible(self, history):
      return hasattr(history, 'lastPopulation')

   """ # obsolete? Doesn't make sense...
   def __getstate__(self):
      self.data = None
      self.network = None
      return self.__dict__
      
   def __setstate__(self,state):
      state['network'] = None
      state['data'] = None
      self.__dict__.update(state)
   """

   def __getstate__(self):
      return {'prem':self.prem,'padd':self.padd,'prev':self.prev}
      
   def __setstate__(self, state):
      self.data = None
      self.network = None
      self.config = None
      self.prem = state['prem']
      self.prev = state['prev']
      self.padd = state['padd']
         
   def canReverse(self, newChild, newParent):
      check = super(StructureProposal, self).canReverse(newChild, newParent)
      if len(newChild.parents) >= self.config.branchFactor:
         return False
      return check
   
   def maybeChange(self, p, changer, existing=True):
      changed = False
      self.network.computeEdgeStatistics()
      try:
         if random.random_sample() < p(self.network):
            # pick a random edge
            if existing:
               if len(self.network.edges) == 0:
                  return False
               index = random.randint(0, len(self.network.edges))
               parent, child = self.network.edges[index]
               changer(parent.index, child, self.data)
            else:
               indexFrom = indexTo = 0
               exists = False
               attempts = 0
               maxAttempts = self.network.numVariables ** 2
               while (indexFrom == indexTo or exists) and (attempts < maxAttempts):
                  indexFrom = random.randint(0, len(self.network.variables))
                  indexTo = random.randint(0, len(self.network.variables))
                  exists = False
                  
                  for edge in self.network.edges:
                     if edge[0] == self.network.variables[indexFrom] and edge[1] == self.network.variables[indexTo]:
                        exists = True
                        break   
                
                  #check branch factor
                  if len(self.network.variables[indexTo].parents) > self.config.branchFactor:
                     exists = True
                     
                  if self.config.cmpIdx is not None:
                     if self.network.variables[indexFrom].index != self.config.cmpIdx and self.network.variables[indexTo].index != self.config.cmpIdx:
                        exists = True
                
                  attempts += 1
                  
               if attempts >= maxAttempts or exists: 
                  return False
               changer(self.network.variables[indexTo], self.network.variables[indexFrom], self.data)
            changed = True
            self.network.computeEdgeStatistics()
      except CyclicException, msg:
         pass
      except DuplicateEdgeException, msg:
         pass
      except IrreversibleEdgeException, msg:
         pass
      except Exception, msg:
         print Exception, msg
         traceback.print_exc(file=sys.stdout)
      return changed

   def batch(self, size, networks = None, data = None, **kwargs):
      self.network = None
      if networks is None:
         if (self.history is not None and
             self.history.lastPopulation() is not None):
            networks = [self.config.space.copy(x)
                        for x,s in self.history.lastPopulation()]
         else:
            networks = [None for i in xrange(size)]
      
      ret = [self.search(net, data, **kwargs) for net in networks]
      self.network = None
      return ret

   def sample(self):
      return self.search()

   def search(self, network=None, data=None, **kwargs):
      if network:
         self.network = network
      else:
         self.network = BayesNet(**self.config.space.config.__properties__)
         self.network.updateVariables(self.config.data)
      if data:
         self.data = data
      else:
         self.data = self.config.data
      self.network.computeEdgeStatistics()
      
      total = 0
      changes = 1
      while total == 0: # changes > 0:
         changes = 0
         
         # remove an edge
         if self.maybeChange(self.prem, self.removeEdge):
            changes += 1 
         
         # reverse an edge
         if self.maybeChange(self.prev, self.reverseEdge):
            changes += 1
         
         # add an edge
         if self.maybeChange(self.padd, self.addEdge, False):
            changes += 1
         total += changes
      self.network.sort()
      network = self.network
      self.network = None
      return network
      
                     
   def adjust(self, acceptanceRatio):
      """Called by simulated annealing to update generation statistics"""

   def densityRatio(self, x, y, i):
      """Called by simulated annealing to adjust for assymetric densities
         The density governing this proposal should be symmetric.
         Still need to check this claim
      """
      return 1.0 class Selection(PopulationDistribution):
   """A selection method (Abstract class)"""
   config = Config(history=SortedMarkovHistory)
   
   def needsScores(self):
      return True
   
   def compatible(self, history):
      return hasattr(history, 'lastPopulation') class Mutation(PopulationDistribution):
   """
      Base class for mutation operators. Provides a method ``mutate`` that performs mutation on a single member of a population.
      
   """
   config = Config(history=MarkovHistory)
   
   def needsScores(self):
      return False
   
   def mutate(self, x):
      """
         Perform mutation on a single organism.
         
         :param x: The organism to mutate.
         :type x: varied
         :returns: The mutated organism.
      """
      return x

   def batch(self, popSize):
      pop = []
      for i, val in enumerate(self.history.lastPopulation()):
         x,s = val
         self.idx = i
         z = self.mutate(x)
         cnt = 0
         while not self.config.space.in_bounds(z):
            cnt += 1
            if cnt > 10000:
               raise RuntimeError("Rejection sampling in mutation failed"
                                  " to generate a point in the space after"
                                  " 10,000 attempts.")
            z = self.mutate(x)
         pop.append(z)
      return pop

   def compatible(self, historyClass):
      return hasattr(historyClass, 'lastPopulation') class Bernoulli(Mutation):
   """
      Bernoulli mutation on a :class:`TernaryString` object. Randomly flips bits on a fully specified binary string.
      
      Config parameters
      * p -- The probability of flipping each bit.   
 
   """
   config = Config(p=0.01)

   def __init__(self, **kwargs):
      super(Bernoulli, self).__init__(**kwargs)
      if (not inspect.isclass(self.config.space.type) or
          not issubclass(self.config.space.type, TernaryString)):
         raise ValueError("Space must produce points descending "
                          "from TernaryString")
      p = self.config.p
      if ((isinstance(p, np.ndarray) and (p > 1.0).any()) or
          (not isinstance(p, np.ndarray) and p > 1.0)):
         raise ValueError("bit flip probability > 1.0: " + str(p))

   def p(self):
      return self.config.p

   def mutate(self, x):
      """This method uses an iterative algorithm to flip bits in
      a TernaryString with a given probability. The algorithm essentially
      uses the binary representation of the bit-flipping probability in order
      to convert random byte sampling (``p=.5``) to account for any probability,
      with resolution of ``1e-16``.
      
      """
      flipped = TernaryString.bernoulli(self.p(), self.config.space.dim)
      base = x.base & ~flipped.base | ~x.base & flipped.base
      known = x.known & flipped.known
      return TernaryString(base, known, self.config.space.dim) class Crossover(PopulationDistribution):
   """
      Performs recombination using a crossover policy given by
      a :class:`Crosser` instance.
   
      Config parameters
      
      * crossoverProb -- A probability that is sampled; if a bernouli test
                         on this value fails, then no crossover is performed,
                         and the first selected organism is chosen
      * order -- The number of organisms to crossover at each step
      * crosser -- A :class:`Crosser` to perform the crossover
      
   """
   config = Config(crossoverProb=1.0, # probability of performing crossover
                   order=2, # number of solutions to crossover
                   crosser=UniformCrosser, # descendant of Crosser
                   history=MultiStepMarkovHistory) # num steps must match order
   
   def __init__(self, **kwargs):
      super(Crossover, self).__init__(**kwargs)
      self.dual = (hasattr(self.config.crosser, 'dual')
                   and self.config.crosser.dual)
   
   def compatible(self, history):
      return (hasattr(history, 'populations') and hasattr(history, 'order')
              and history.order() == self.config.order)
   
   def batch(self, popSize):
      if self.dual:
         psize = popSize / 4
      else:
         psize = popSize
      crossoverProb = self.config.crossoverProb
      if crossoverProb < 1e-16:
         return [x for x,s in self.history.populations[0]]
      self.crosser = self.config.crosser(self.config)
      pops = self.history.populations
      newpop = [self.crosser([org for org, s in orgs], crossoverProb)
                for orgs in zip(*pops)]
      if self.dual:
         pop = []
         for x,y in newpop:
            pop.append(x)
            pop.append(y)
         newpop = pop
      return newpop class Gaussian(Mutation):
   """
      Gaussian mutation. Randomly mutates a random selection of components
      of a real vector.
      
      Config parameters
      * sd -- The standard deviation to use. 
      * p -- The probability of mutating each component; defaults to 1.0
      
  """
   config = Config(sd=0.01,
                   p=1.0)
   
   def __init__(self, **kwargs):
      super(Gaussian, self).__init__(**kwargs)
      if self.config.space.type != np.ndarray:
         raise ValueError("Space for Gaussian must have type numpy.ndarray")
         
   def sd(self):
      return self.config.sd
         
   def mutate(self, x):
      p = np.random.binomial(1, self.config.p, len(x))
      ret = x + p * np.random.randn(len(x)) * self.sd()
      if not self.config.space.in_bounds(ret):
         scale = self.config.space.scale
         center = self.config.space.center
         ret = np.maximum(np.minimum(ret, scale+center),center-scale)
      return ret
      
   def density(self, center, point):
      sd = self.sd()
      return np.exp(-(1. / (2. * sd * sd)) *
                    ((point - center) ** 2).sum()) / sd / np.sqrt(2*np.pi) class AreaSensitiveGaussian(Gaussian):
   """
      Gaussian mutation with a variance that depends on the size of a
      partition of the space (see :class:`pyec.util.partitions.Partition`).
      
      Config parameters:
      * jogo2012 -- Whether to use the methods from Lockett & Miikkulainen, 2012
      * decay -- A function ``decay(n,config)`` to compute the multiplier
                 that controls the rate of decrease in standard deviation.
                 Faster decay causes faster convergence, but may miss the
                 optimum. Default is ``((1/generations))``
      
      Standard deviation is determined by::
      
         sd = .5 * ((upper-lower) * decay(n)
         
      where ``dim`` is the dimension of the space (``config.dim``) and 
      ``area`` is the volume of the partition region for the object being mutated.
      
   """
   config = Config(jogo2012=False,
                   decay=lambda n,cfg: (1./n)) # **(1./cfg.space.dim))
   
   def __init__(self, **kwargs):
      super(AreaSensitiveGaussian, self).__init__(**kwargs)
      if self.config.jogo2012:
         self.config.decay = lambda n,cfg: (1/n)**.5
   
   def compatible(self, history):
      # N.B. we use a Markov history, but assume we are passed the partition nodes
      return hasattr(history, 'lastPopulation')
   
   def mutate(self, x):
      """Apply area-sensitive gaussian mutation.
         Assume that the passed value is a :class:`Point` object
         
         :param x: The point to mutate
         :type x: :class:`Point`
         :returns: A numpy.ndarray mutation the ``point`` property of ``x``.
      
      """
      center = self.config.space.center
      scale = self.config.space.scale
      
      y, node = x
      y = y.point
      if self.config.jogo2012:
         area = node.area
         # sd = 2 * scale / (-log(area))
         if area < 0.0: 
            area = 1e-100
         sd = self.config.sd * scale * (area ** (1./len(y)))
         sd *= self.config.decay(self.history.updates, self.config)
      else:
         area = node
         lower, upper = area.bounds.extent()
         sd = .5 * (upper - lower)
         sd = np.minimum(sd, scale)
         sd *= self.config.decay(self.history.updates, self.config)
      
      ret = y + np.random.randn(len(y)) * sd
      if not self.config.space.in_bounds(ret):
         ret = np.maximum(np.minimum(ret, scale+center),center-scale)
      
      return ret

   def density(self, info, point):
      scale = self.config.space.scale
      center, area = info
      sd = self.config.sd * scale * (area ** (1./len(point)))
      
      return exp(-(1./ (2 * sd * sd)) *
                 ((point - center) ** 2).sum()) / sd / np.sqrt(2*np.pi)  class AreaSensitiveBernoulli(Bernoulli):
   """
      Like :class:`AreaSensitiveGaussian` but with a Bernoulli.
      
      Because a bit space is finite, there is no need to decay, but
      we still want to make exiting the area less likely to promote
      exploration within the selected area.
      
      Config parameters
      * p -- The bit-flipping prob
      
   """
   _p = None
   
   def p(self):
      return self._p
   
   def density(self, x, z):
      y = x[0].point
      area = x[1]
      bf = self._p
      prod = 1.0
      for i in xrange(self.config.dim):
         if z[i] != y[i]:
            prod *= bf
         else:
            prod *= (1. - bf)
      return prod
      
   def mutate(self, x):
      y = x[0].point
      area = x[1]
      #lower, upper = area.bounds.extent()
      #diff = TernaryString(lower.base ^ upper.base,-1L,lower.length).toArray()
      #logArea = self.config.space.dim - diff.sum()
      #self._p = (diff * self.config.p  +
      #           (1.0 - diff) * self.config.p ** logArea)
      self._p = self.config.p
      return super(AreaSensitiveBernoulli, self).mutate(y) class AddChainLayerMutation(Mutation):
    """Add a layer to an RNN copied from another layer,
    with a link going from the copied layer to the new layer
    whose weights are the identity matrix, and a
    link going to a random layer with random weights.
   
    Config parameters:
   
      * layer_creation_prob - The probability of adding a layer
      * sd - The std dev for the Gaussian to generate random weights
      
    """
    config = Config(layer_creation_prob=.01,
                    sd=1.0)
   
    def mutate(self, net):
        if net.changed:
            return net
      
        if np.random.random_sample() > self.config.layer_creation_prob:
            return net
      
        net = self.config.space.copy(net)
        satisfactory = lambda lay: (not isinstance(lay, RnnLayerInput) and
                                    not isinstance(lay, RnnLayerBias))
        targets = [layer for layer in net.layers if satisfactory(layer)]
        target = targets[np.random.randint(0,len(targets))]
        source = net.layers[np.random.randint(0,len(targets))]
        middle = RnnLayer(source.size,activator=self.config.space.activator)
        net.addLayer(middle)
        net.connect(source, middle, np.identity(source.size))
        #if (source, target) in net.links:
        #    net.connect(middle, target, net.links[(source, target)].copy())
        #else:
        w = self.config.sd * np.random.randn(target.size, middle.size)
        net.connect(middle, target, w)
        net.changed = True
        #net.checkLinks()
        return net class AddRadialLayerMutation(Mutation):
    """Add a layer to an RNN copied from another layer,
    with a link going from the copied layer to the new layer
    whose weights are the identity matrix, and a
    link going to a random layer with random weights.
   
    Config parameters:
   
      * layer_creation_prob - The probability of adding a layer
      * sd - The std dev for the Gaussian to generate random weights
      
    """
    config = Config(layer_creation_prob=.01,
                    sd=1.0)
   
    def mutate(self, net):
        if net.changed:
            return net
      
        if np.random.random_sample() > self.config.layer_creation_prob:
            return net
      
        net = self.config.space.copy(net)
        satisfactory = lambda lay: (not isinstance(lay, RnnLayerInput) and
                                    not isinstance(lay, RnnLayerBias))
        targets = [layer for layer in net.layers if satisfactory(layer)]
        target = targets[np.random.randint(0,len(targets))]
        source = net.layers[np.random.randint(0,len(targets))]
        middle = RnnLayer(source.size,activator=RADIAL)
        net.addLayer(middle)
        net.connect(source, middle, np.identity(source.size))
        #if (source, target) in net.links:
        #    net.connect(middle, target, net.links[(source, target)].copy())
        #else:
        w = self.config.sd * np.random.randn(target.size, middle.size)
        net.connect(middle, target, w)
        net.changed = True
        #net.checkLinks()
        return net class RemoveLayerMutation(Mutation):
    """Remove a random hidden layer.
   
    Config parameters:
   
      * layer_deletion_prob - The probability of removing a layer
   
    """
    config = Config(layer_deletion_prob=0.01)
   
    def mutate(self, net):
        if net.changed:
            return net
        if np.random.random_sample() > self.config.layer_deletion_prob:
            return net
        net = self.config.space.copy(net)
        hiddenLayers = net.hiddenLayers()
        if len(hiddenLayers) == 0:
            return net
        toRemove = hiddenLayers[np.random.randint(0,len(hiddenLayers))]
        net.removeLayer(toRemove)
        net.changed = True
        #net.checkLinks()
        return net class AddNodeMutation(Mutation):
   """Add a node to a neural network.
   
   Config parameters:
      
      * node_creation_prob - The probability of adding a node
      * sd - The standard deviation for the Gaussian used to make weights
                  for any new nodes
   
   """
   config = Config(node_creation_prob=.01,
                   sd=1.0)
   
   def mutate(self, net):
      if net.changed:
         return net
      p = self.config.node_creation_prob
      if np.random.random_sample() > p:
         return net
      
      sd = self.config.sd
      hiddenLayers = net.hiddenLayers()
      if len(hiddenLayers) == 0:
         return net
      
      
      net = self.config.space.copy(net)
      layer = net.hiddenLayers()[np.random.randint(0,len(hiddenLayers))]
      layer.size += 1
      for inLayer in layer.inLinks:
          w0 = net.links[(inLayer, layer)]
          w = np.zeros((layer.size, inLayer.size))
          if inLayer == layer:
              w[:(layer.size-1),:(layer.size-1)] = w0
              w[:,layer.size-1] = sd * np.random.randn(layer.size)
          else:
              w[:(layer.size-1),:] = w0
          w[layer.size-1,:] = sd * np.random.randn(inLayer.size)
          net.connect(inLayer, layer, w)
      
      for outLayer in layer.outLinks:
          if outLayer == layer:
              continue # we already did this!
          w0 = net.links[(layer, outLayer)]
          w = np.zeros((outLayer.size, layer.size))
          w[:,:(layer.size-1)] = w0
          w[:,layer.size-1] = sd * np.random.randn(outLayer.size)
          net.connect(layer, outLayer, w)
      
      net.changed = True
      #net.checkLinks()
      return net class RemoveNodeMutation(Mutation):
   """Delete a node from a neural network.
   
   Config parameters:
      
      * node_deletion_prob - The probability of deleting a node
   
   """
   config = Config(node_deletion_prob=.01)
   
   def mutate(self, net):
      if net.changed:
         return net
      p = self.config.node_deletion_prob
      if np.random.random_sample() > p:
         return net
      
      hiddenLayers = net.hiddenLayers()
      if len(hiddenLayers) == 0:
         return net
      
      net = self.config.space.copy(net)
      
      layer = net.hiddenLayers()[np.random.randint(0,len(hiddenLayers))]
      if layer.size <= 1:
         return net
      
      indexToRem = np.random.randint(0, layer.size)
      layer.size -= 1
      
      for inLayer in layer.inLinks:
          w0 = net.links[(inLayer, layer)]
          w = np.delete(w0, indexToRem, axis=0)
          if inLayer == layer:
              w = np.delete(w, indexToRem, axis=1)
          net.connect(inLayer, layer, w)
      
      for outLayer in layer.outLinks:
          if outLayer == layer:
              continue
          w0 = net.links[(layer, outLayer)]
          w = np.delete(w0, indexToRem, axis=1)
          net.connect(layer, outLayer, w)
      
      net.changed = True
      #net.checkLinks()
      return net class AddLinkMutation(Mutation):
   """Add a link to a network.
   
   Config parameters:
   
      * link_creation_prob - The probability of adding a random link
      * sd - The std dev for the Gaussian used to create the weights
   
   """
   config = Config(link_creation_prob=0.025,
                   sd=1.0)
   
   def mutate(self, net):
      if net.changed:
         return net
      
      p = self.config.link_creation_prob
      if np.random.random_sample() > p:
         return net
      
      # randomly pick two layers
      satisfactory = lambda lay: (not isinstance(lay, RnnLayerInput) and
                                  not isinstance(lay, RnnLayerBias))
      targets = [layer for layer in net.layers if satisfactory(layer)]
      target = targets[np.random.randint(0,len(targets))]
      source = net.layers[np.random.randint(0,len(targets))]
      if not (source, target) in net.links:
         net = self.config.space.copy(net)
         source = [layer for layer in net.layers if layer.id == source.id][0]
         target = [layer for layer in net.layers if layer.id == target.id][0]
         w = self.config.sd * np.random.randn(target.size, source.size)
         net.connect(source, target, w)
         net.changed = True
         #net.checkLinks()
      
      return net class RemoveLinkMutation(Mutation):
   """Remove a random link from a network.
   
   Config parameters:
   
      * link_deletion_prob - The probability of removing a random link
   
   """
   config = Config(link_deletion_prob=0.025)
   
   def mutate(self, net):
      if net.changed:
         return net
      
      p = self.config.link_deletion_prob
      if np.random.random_sample() > p:
         return net
      
      # randomly pick two layers
      satisfactory = lambda lay: (not isinstance(lay, RnnLayerInput) and
                                  not isinstance(lay, RnnLayerBias))
      targets = [layer for layer in net.layers if satisfactory(layer)]
      target = targets[np.random.randint(0,len(targets))]
      source = net.layers[np.random.randint(0,len(targets))]
      if (source, target) in net.links:
         net = self.config.space.copy(net)
         source = [layer for layer in net.layers if layer.id == source.id][0]
         target = [layer for layer in net.layers if layer.id == target.id][0]
         net.disconnect(source, target)
         net.changed = True
         #net.checkLinks()
      return net class NetAreaStripper(Mutation):
    def mutate(self, x):
        point, node = x
        net = self.config.space.copy(point.point)
        net.changed = False
        return net class AreaSensitiveGaussianWeightMutation(GaussianWeightMutation):
    """Weight mutation, sensitive to area partition.
    
    Config parameters:
      * decay -- A function ``decay(n,config)`` to compute the multiplier
                 that controls the rate of decrease in standard deviation.
                 Faster decay causes faster convergence, but may miss the
                 optimum. Default is ``((1/generations))``
    
    """
    config = Config(decay=lambda n,cfg: (1./(2*np.log(n))))
    
    def gaussInt(self, z):
        # x is std normal from zero to abs(z)
        x = .5 * erf(np.abs(z)/np.sqrt(2))
        return .5 + np.sign(z) * x
    
    def sd(self, net):
        try:
            area, path = self.config.segment.partitionTree.traverse(net)
            lower, upper = area.bounds.extent()
            sd = .5 * (upper - lower)
            scale = self.config.space.scale
            sd = self.gaussInt(upper/scale) - self.gaussInt(lower/scale)
            sd *= self.config.decay(self.history.updates, self.config)
        except Exception:
            n = self.history.updates
            return self.config.sd * self.config.decay(n, self.config) class UniformOrIntermediateRnnCrosser(UniformRnnCrosser):
    def __init__(self, *args, **kwargs):
        super(UniformOrIntermediateRnnCrosser, self).__init__(*args, **kwargs)
        self._intermediate = np.random.random_sample() < .4
        
    def intermediate(self, w1, w2):
        return .5 * (w1+w2)
    
    def crossLink(self, w1, w2):
        if self._intermediate:
            return self.intermediate(w1, w2)
        else:
            return super(UniformOrIntermediateRnnCrosser, self).crossLink(w1, w2) class UniformRnnLinkCrosser(UniformRnnCrosser):
    """Crossover for RNNs with entire matrices chosen from one RNN or the other.
    
    """
    def crossLink(self, w1, w2):
        """Just picks one of the two links; no crossover within matrices."""
        if np.random.random_sample() < .5:
            return w1.copy()
        return w2.copy() class History(object):
    """A History used to track the progress of an optimization algorithm.
    
    Different algorithms should extend this class in order to define
    the minimal amount of history that needs to be stored in order for the
    algorthm to operate.
    
    """
    useCache = True
    attrs = set()
    sorted = False
    root = True # whether the current history is the top level root for
                # all histories of this algorithm
    
    def __init__(self, config):
        super(History, self).__init__()
        self.config = config
        self.evals = 0
        self.minSolution = None
        self.minScore = np.inf
        self.maxSolution = None
        self.maxScore = -np.inf 
        self._empty = True
        if not hasattr(self, 'cache'):
            self.cache = LRUCache() # 10,000 items by default
        self.updates = 0
        #how often to print generation report
        self.printEvery = config.printEvery or 1000000000000L 
        self.attrs = set(["evals","minSolution","minScore","maxScore","attrs",
                          "minSolution", "maxSolution","_empty","cache",
                          "updates","printEvery", "useCache"])
    
    def __getstate__(self):
        """Used by :class:`CheckpointedHistory`
        and :class:`CheckpointedMultipleHistory` to checkpoint a history
        so it can be rolled back after updates.
        
        Should return all objects in the history in a dictionary, sensitive
        to the fact that object references may need to be copied.
        
        :returns: A dictionary with the state of the history
        
        """
        state = {}
      
        for attr in self.attrs:
            val = getattr(self, attr)
            if isinstance(val, list):
                val = [x for x in val]
            state[attr] = val
         
        state['cfg'] = self.config.__properties__
        return state

    def __setstate__(self, state):
        for attr in self.attrs:
            val = state[attr]
            setattr(self, attr, val)
         
        import pyec.config
        self.config = pyec.config.Config(**state['cfg'])  
    
    def empty(self):
        """Whether the history has been used or not."""
        return self._empty
    
    def better(self, score1, score2):
        """Return whether one score is better than another.
        
        Uses ``config.minimize`` to decide whether lesser or greater
        numbers are better.
        
        :param score1: the score in question (floating point number)
        :type score1: ``float``
        :param score2: the score being compared
        :type score2: ``float``
        :returns: whether ``score1`` is better than ``score2``
        """
        if self.config.minimize:
            return score1 < score2
        else:
            return score1 > score2
    
    def best(self):
        """Get the best solution, whether minimizing or maximizing.
        
        Same as ``optimal``
        
        """
        if self.config.minimize:
            return self.minimal()
        else:
            return self.maximal()
    
    optimal = best
    
    def minimal(self):
        """Get the minimal solution and its score.
        
        :returns: A tuple of two item with the solution object first
                  and the score of that item second 
                  
        """
        return self.minSolution, self.minScore
        
    def maximal(self):
        """Get the maximal solution and its score.
        
        :returns: A tuple of two item with the solution object first
                  and the score of that item second 
                  
        """
        return self.maxSolution, self.maxScore
        
    def num_evaluations(self):
        """Get the number of function evaluations performed in this 
        history object.
        
        :returns: An integer representing the number of times the fitness
                  function or objective has been evaluated
                  
        """
        return self.evals
        
    def update(self, population, fitness, space, opt):
        """
         Update the state of the :class:`History` with the latest population 
         and its fitness scores. Subclasses should probably override
         ``internalUpdate`` rather than ``update``, unless they want to
         change how the min/max are tracked.
         
         Returns the history for use in continuations.
         
         If ``population`` is ``None``, then this method does nothing; this
         is so that you can set up a loop like to run an optimizer like::
         
         p = None
         f = lambda x: 
         t = History()
         o = some_optimizer
         s = o.config.space
         for i in xrange(generations):
             p = some_optimizer[t.update(p,f,s,0), f]()
         t.update(p,f,s,o)
          
         :params population: The previous population.
         :type population: list of points in the search domain
         :params fitness: The fitness / cost / objective function
         :type fitness: Any callable object
         :params space: The search domain
         :type space: :class:`Space`
         :params opt: The optimizer reporting this population
         :type opt: :class:`PopulationDistribution`
         :returns: The history (``self``), for continuations
         
        """
        if population is None:
            return
        
        #self.config.stats.start(repr(self) + "history.update.all")
        self._empty = False
        self.evals += len(population)
        self.updates += 1
        
        #self.config.stats.start(repr(self) + "history.update.scoreall")
        # score the sample
        pop  = population
        scored = [(x, self.score(x, fitness, space)) for x in pop]
        #self.config.stats.stop(repr(self) + "history.update.scoreall")
        #self.config.stats.start(repr(self) + "history.update.findbest")
        
        if self.root and self.config.observer is not None:
            self.config.observer.report(opt, scored)
        
        for x,s in scored:
            if s > self.maxScore:
                 self.maxScore = s
                 self.maxSolution = x
               
            if s < self.minScore:
                 self.minScore = s
                 self.minSolution = x
        #self.config.stats.stop(repr(self) + "history.update.findbest")
        
        if not (self.updates % self.printEvery):
            genmin = min([s for x,s in scored])
            genmax = max([s for x,s in scored])
            genavg = np.average([s for x,s in scored])
            print self.updates, ": min", self.minScore, " max", self.maxScore,
            print " this generation (min, avg, max): ", genmin, genavg, genmax
        
        #self.config.stats.start(repr(self) + "history.update.internal")       
        self.internalUpdate(scored)
        #self.config.stats.stop(repr(self) + "history.update.internal")
        #self.config.stats.stop(repr(self) + "history.update.all")
        return self
        
    def internalUpdate(self, population):
        """
         Update the state of the :class:`History` with the latest population 
         and its fitness scores. This is an internal call intended for
         overridden by subclasses. One of the important functions is to
         delete points no longer needed by the history.
          
         :params population: The previous population with its fitness scores. 
         :type population: list of (point, score) tuples
         
        """    
        pass

    def score(self, point, fitness, space):
        """Get the fitness score, caching where possible.
      
         :param point: A valid point in the space
         :type point: Must match ``space.type``
         :param fitness: The fitness function
         :type fitness: Any callable
         :param space: The space to which the point belongs
         :type space: :class:`Space`
         :returns: The fitness value, cached if possible
         
        """
        if fitness is None: 
            return None
        
        #self.config.stats.start("history.score")
            
        if self.useCache: 
            try:
                hashed = space.hash(point)
                if self.cache.has_key(hashed):
                    ret = self.cache[hashed]
                    #self.config.stats.stop("history.score")
                    return ret
            except Exception:
                pass
        
        if not space.in_bounds(point):
            # use NaN so that the result is less than nor greater than
            # any other score, and therefore NEVER optimal
            s = np.inf - np.inf
        else:
            try:
                s = fitness(space.convert(point))
            except ValueError:
                s = np.inf - np.inf
        
        if self.useCache:
            try:
                hashed = space.hash(point)
                self.cache[hashed] = s
            except Exception:
                pass
        #self.config.stats.stop("history.score")
        return s
    
    def setCache(self, cache):
        self.cache = cache class Euclidean(Space):
    """A Euclidean space of fixed finite dimension.
    
    Uses a numpy array with dtype 64-bit floats.
    
    :param dim: The dimension of the space
    :type dim: ``int``
    :param scale: The scale for the space, integer or array
    :type scale: ``int`` for spherical space, or ``numpy.ndarray``
    
    """
    def __init__(self, dim=1, center=0.0, scale=1.0):
        super(Euclidean, self).__init__(np.ndarray)
        if not isinstance(center, np.ndarray):
            center = center * np.ones(dim)
        if not isinstance(scale, np.ndarray):
            scale = scale * np.ones(dim)
        try:
            bottom = center - scale
        except Exception:
            raise ValueError("Mismatched center or scale in Euclidean")
        self.center = center
        self.scale = scale
        self.dim = dim
        
        if (dim,) != np.shape(center):
            raise ValueError("Dimension of center doesn't match dim")
        
        if (dim,) != np.shape(scale):
            raise ValueError("Dimension of scale array doesn't match dim")
    
    def gaussInt(self, z):
        # x is std normal from zero to abs(z)
        x = .5 * erf(np.abs(z)/np.sqrt(2))
        return .5 + np.sign(z) * x
    
    def proportion(self, smaller, larger, index):
        """Assume ``smaller`` is a hyperrectangle, and ``larger`` is either
        a euclidean space or a hyperrectangle containing ``smaller``.
        
        To handle something more general, we would need to integrate somehow,
        either monte carlo or decomposing into hyperrectangles.
        
        :param smaller: A Hyperrectangle in this space
        :type smaller: :class:`Hyperrectangle`
        :param larger: A Hyperrectangle or this space, either way containing
                       ``smaller``
        :type larger: :class:`Hyperrectangle` or :class:`Euclidean`
        :returns: The ratio of ``smaller``'s volume over ``larger``'s.
        """
        center = self.center[index]
        scale = self.scale[index]
        slower = smaller.lower[index]
        supper = smaller.upper[index]
        if isinstance(larger, Hyperrectangle):
            llower = larger.lower[index]
            lupper = larger.upper[index]
        else: # Euclidean
            llower = -np.inf
            lupper = np.inf
        slow = self.gaussInt((slower - center) / scale) 
        shigh = self.gaussInt((supper - center) / scale)
        llow = self.gaussInt((llower - center) / scale) 
        lhigh = self.gaussInt((lupper - center) / scale)
        return (shigh - slow) / (lhigh - llow)
      
    def extent(self):
        upper = np.zeros(self.dim)
        lower = np.zeros(self.dim)
        upper.fill(np.inf)
        lower.fill(-np.inf)
        return lower, upper
      
    def random(self):
        """Get a random point in Euclidean space. Use the constraint to 
        generate a random point in the space if possible, otherwise
        use a zero-centered elliptical gaussian scaled by ``self.scale``.
        
        """
        test = self.center + self.scale * np.random.randn(self.dim)
        return test
                
    def hash(self, point):
        parts = [((point+i)**2).sum() for i in np.arange(10)]
        return ",".join([str(pt) for pt in parts])
    
    def in_bounds(self, point, **kwargs):
        return isinstance(point, self.type) and np.shape(point) == (self.dim,)

    def copy(self, point):
        return np.copy(point) class Binary(Space):
    """A binary space of fixed finite dimension.
    
    Uses a :class:`TernaryString` as a representation.
    
    :param dim: The dimension of the space
    :type dim: ``int``
    
    """
    _area = None
    
    def __init__(self, dim=1):
        super(Binary, self).__init__(TernaryString)
        self.dim = dim
        
    def area(self):
        return 1.0
            
    def extent(self):
        return TernaryString(0L, 0L, self.dim), TernaryString(-1L, 0L, self.dim)
            
    def random(self):
        """Get a random point in binary space. Use the constraint to 
        generate a random point in the space if possible, otherwise
        use a random byte string.
        
        """
        return TernaryString.random(self.dim)
    
    def in_bounds(self, x, **kwargs):
        return isinstance(x, TernaryString) and x.length == self.dim class fact(object):
   def __init__(self):
   def __call__(self, n):
   def choose(self, n,m): class Segment(object):

   def __init__(self, config):
      self.config = config
      self.taylorCenter = config.taylorCenter or 1.0
      self.taylorDepth = config.taylorDepth or 0
      
      self.points = []
      self.partitionTree = Partition(self,config)
      self.scoreTree = ScoreTree(self,config)
      
   def clearSegment(self):
      del self.points
      del self.partitionTree
      del self.scoreTree
      self.points = []
      self.partitionTree = Partition(self,self.config)
      self.scoreTree = ScoreTree(self,self.config) class Partition(object):
   areaTree = None
   
   def __init__(self, segment, config):
      self.config = config
      self.segment = segment
               
      self.root = PartitionTreeNode(
         segment = self.segment,
         index = 0,
         parent = None,
         upper = None,
         lower = None,
         area = config.space.area(),
         bounds = config.space,
         point = None
      )
      
      self.areaTree = AreaTree()
      
   def save(self):
      pass

   def traverse(self, point):
      """
         Traverse the partition tree to find the partition within which "point" is located.
         
         :param tree: The partition tree
         :type tree: :class:`Partition`
         :param point: The point for searching.
         :type point: ``config.space.type``
         :returns: A tuple with (partitionNode, path); with the
                   :class:`PartitionTreeNode` for the matching partition, and
                   the list of tree nodes traversed

      """
      if isinstance(point, Point):
         point = point.point
      
      # get the top parent
      current = self.root
      region = self.root.bounds
      path = [current]
      
      last = None
      children = current.children
      while len(children) > 0:
         if current == last:
            print "partition.traverse looped, exception"
            print current.bounds.extent()
            print point.point
            raise Exception("Loop in partition.traverse!")
         last = current
         
         for child in children:
            if child.bounds.in_bounds(point, index=child.index):
               current = child
               path.append(current)
               children = current.children
               break
            
      node = current
      return node, path

   def propagateScoreSum(self, node, config):
      """
         Propagate the score up the tree
         
         node - the node to start at (bottom-up traversal)
         config - an evo.config.Config
         
      """
      current = node
      while current is not None:
         current.score_sum = sum([child.score_sum for child in current.children])
         current = current.parent
      
   def computeDepth(self, point):
      """
         compute the depth of the leaf that points to the db_point with id pointId
         
      """
      depth = 0
      current = point.partition_node
      while current is not None:
         current = current.parent
         if current is not None:
            depth += 1
      
      return depth
      
   def printTree(self, segment, node=None, indent=""):
      if node is None:
         node = self.root
      print indent, node.id, ": ", node.layer, node.index, node.bounds, node.area
      children = node.children
      for child in children:
         self.printTree(segment, child, indent + "\t")      

   def largestArea(self):
      try:
         return self.areaTree.largest()
      except AreaException:
         return 1.0 class ScoreTree(object):
   
   def __init__(self, segment, config):
      self.config = config
      self.segment = segment
      self.root = ScoreTreeNode(
         segment = self.segment,
         parent = None,
         point = None
      )
      
   def save(self):
      pass

   def printTree(self, segment, node=None, indent=''):
      """
         Recursively print the score tree.
         
         segment - the segment for the tree to print
         node - the node to start at, or None for the root
         indent - a prefix for printing; tabs added for each new level
      """
      if node is None: 
         node = self.root
      
      print indent, node.id, ',', node.balance, ",", node.area, ",", node.height,
      print ",", node.child_count, ",", node.min_score, " - ", node.max_score
      
      children = node.children
      if len(children) > 0:
         for child in sorted(children, key=lambda c: not c.left):
            self.printTree(segment, child, indent + '\t')

   def checkBalance(self, node=None):
      if node is None:
         node = self.root
      
      try:
         children = node.children
         if len(children) > 0:
            children = sorted(children, key=lambda c: not c.left)
            assert (node.balance == (children[0].height - children[1].height))
            for child in children:
               self.checkBalance(child)
      except AssertionError:
         print "NODE: ", node.id, node.balance, children[0].height, children[1].height
         self.printTree(self.segment)
         raise RuntimeError("Boo")
      except:
         raise

   def rotateLeft(self, node, config):
      """
         Perform a left rotation at a node in the tree. For balancing.
         
         node - the node to rotate left.
         config - a config object with property "shiftToDb"
         
         Returns -1, the change in tree height. This return will only be correct
         if used in the context of an AVL tree balancing algorithm - otherwise
         the height change would need to be computed differently (it is possible).
      """
      center = node
      children = sorted(node.children, key=lambda c: not c.left)
      left = children[0]
      right = children[1]
      
      right.left = center.left
      center.left = True
      left.left = True
      
      children = sorted(right.children, key=lambda c: not c.left)
      rightleft = children[0]
      rightright = children[1]   
      
      rightleft.parent = center
      rightleft.left = False
      right.parent = center.parent
      center.parent = right

      center.area = left.area + rightleft.area
      right.area = center.area + rightright.area
      
      center.taylor = left.taylor + rightleft.taylor
      right.taylor = center.taylor + rightright.taylor
      
      center.min_score = rightleft.min_score
      center.max_score = left.max_score
      right.min_score = rightright.min_score
      right.max_score = center.max_score
      
      
      center.height = max([left.height, rightleft.height]) + 1
      right.height = max([center.height, rightright.height]) + 1

      center.child_count = left.child_count + rightleft.child_count + 1
      right.child_count = center.child_count + rightright.child_count + 1

      center.balance = center.balance + 1
      if right.balance < 0:
         center.balance = center.balance - right.balance
      
      right.balance = right.balance + 1
      if center.balance > 0:
         right.balance = right.balance + center.balance
      
        
      center.children = [left,rightleft]
      right.children = [center,rightright]
      if right.parent is None:
         right.segment.scoreTree.root = right
      else:
         newChildren = [right]
         for child in right.parent.children:
            if child != center:
               newChildren.append(child)
         right.parent.children = newChildren 
         
      return -1

   def rotateRight(self, node, config):
      """
         Perform a right rotation at a node in the tree. For balancing.
         
         node - the node to rotate right.
         config - a config object with property "shiftToDb"
         
         Returns -1, the change in tree height. This return will only be correct
         if used in the context of an AVL tree balancing algorithm - otherwise
         the height change would need to be computed differently (it is possible).
      """
      center = node
      children = sorted(node.children, key=lambda c: not c.left)
      
      left = children[0]
      right = children[1]
      
      left.left = center.left
      center.left = False
      right.left = False
      
      children = sorted(left.children, key=lambda c: not c.left)
      leftleft = children[0]
      leftright = children[1]

      leftright.parent = center
      leftright.left = True
      left.parent = center.parent
      center.parent = left

      center.area = right.area + leftright.area
      left.area = center.area + leftleft.area
      
      center.taylor = right.taylor + leftright.taylor
      left.taylor = center.taylor + leftleft.taylor
   
      center.min_score = right.min_score
      center.max_score = leftright.max_score
      left.min_score = center.min_score
      left.max_score = leftleft.max_score
      

      center.height = max([right.height, leftright.height]) + 1
      left.height = max([center.height, leftleft.height]) + 1

      center.child_count = leftright.child_count + right.child_count + 1
      left.child_count = center.child_count + leftleft.child_count + 1

      center.balance = center.balance - 1
      if left.balance > 0:
         center.balance = center.balance - left.balance

      left.balance = left.balance - 1
      if center.balance < 0:
         left.balance = left.balance + center.balance

      center.children = [leftright,right]      
      left.children = [leftleft,center]
      if left.parent is None:
         left.segment.scoreTree.root = left
      else:
         newChildren = [left]
         for child in left.parent.children:
            if child != center:
               newChildren.append(child)
         left.parent.children = newChildren 
      
      return -1
      
      
   def traverse(self, point, config):
      """
         Traverse the score tree.
         
         point - the point to compare
         config - an evo.config.Config instance with property "shiftToDb"
         
         Returns node in the tree that must be split to insert the new point into the tree.
      """
      current = self.root
      while True:
         children = sorted(current.children, key=lambda c: not c.left)
         if len(children) == 0:
            break
         if children[1].max_score < point.score:
            current = children[0]
         else:
            current = children[1]
      return current

   def removeLeaf(self, node, config):
      if node.parent is None:
         if node.point is not None:
            node.point.score_node = None
         node.point = None
         return
          
      parent = node.parent
      other = [child for child in node.parent.children if child != node][0]
      
      left = True
      for i,child in enumerate(parent.parent.children):
         if child == parent:
            left = i == 0 
            break
      
      parent.area = other.area
      parent.point = other.point
      parent.min_score = other.min_score
      parent.max_score = other.max_score
      parent.taylor = other.taylor
      parent.balance = 0
      parent.children = []
      del node
      del other
      parent.height -= 1
      if parent.parent is not None:
         #print "before, ", parent.parent_id
         #self.printTree(parent.segment)
      
      
         superparent = parent.parent
         if left:
            superparent.balance -= 1
         else:
            superparent.balance += 1
         if superparent.balance == 0:
            superparent.height -= 1
         
         self.propagateArea(superparent, config, -1)
         
         #print "after, ", superparent.id
         #self.printTree(parent.segment)
      

   def propagateAreaOnly(self, node, config, diff):
      """
         Propagate the change in area given by "diff" 
      """
      parent = node
      while parent is not None:
         parent.area -= float(diff)
         parent = parent.parent
      

   def propagateArea(self, node, config, inserted=1):
      """
         Propagate the area and other features up the score tree, balancing as needed using an AVL algorithm.
         
         For balancing, positive numbers are heavy on the left, negative numbers are heavy on the right.
         
         node - the node to start at (bottom-up traversal)
         config - an evo.config.Config object with property "shiftToDb"
         inserted - how many new levels have been added to the tree, 1 or 0. If 0, no balancing is needed.
      """
      
      if inserted > 0:
         heightChange = inserted
      elif inserted < 0:
         if node.balance == 0:
            heightChange = -1
         else:
            heightChange = 0
      else:
         heightChange = 0
      current = node
      
      children = sorted(current.children, key=lambda c: not c.left)
      while True:
      
         left = children[0]
         right = children[1]
         
         if left.min_score < right.max_score:
            self.printTree(node.segment)
            raise Exception, "Violated ordered assumption!"
            
         
         # is the current the right or the left child?
         onLeft = current.left
         
         
         rotated = False
         if inserted != 0:
            if current.balance > 1:
               if left.balance != 0:
                  heightChange -= 1
               if left.balance >= 0:
                  rotated = True
                  self.rotateRight(current, config)
               else:# left.balance < 0:
                  rotated = True
                  self.rotateLeft(left, config)
                  self.rotateRight(current, config)
            elif current.balance < -1:
               if right.balance != 0:
                  heightChange -= 1
               if right.balance <= 0:
                  rotated = True
                  self.rotateLeft(current, config)
               else:# right.balance > 0:
                  rotated = True
                  self.rotateRight(right, config)
                  self.rotateLeft(current, config)
               
         
         if not rotated:
            current.area = left.area + right.area
            current.taylor = left.taylor + right.taylor
            current.child_count = left.child_count + right.child_count + 1
            current.min_score = right.min_score
            current.max_score = left.max_score
         
         next = None
         if current.parent is not None:
            next = current.parent
            children = sorted(next.children, key=lambda c: not c.left)
         
         if inserted != 0 and next is not None:
            oldbal = next.balance
            if onLeft:
               newbal = next.balance + heightChange
            else:
               newbal = next.balance - heightChange
            next.balance = newbal
            if inserted > 0:
               if newbal == 0:
                  heightChange = 0
            else:
               if newbal != 0:
                  heightChange = 0
            next.height += heightChange
         
         if next is None:
            break
         current = next
         
      #if inserted != 0:
      #   print "after, ", node.id, inserted
      #   self.printTree(node.segment)
      
      
   def insert(self, point, config, stats):
      """
         Insert a new point into the score tree, then propagate area, taylor scores, and rebalance.
         Assumes the point has already been inserted into the partition tree
         
         point - the point to insert
         config - an evo.config.Config object with properties "shiftToDb", "taylorDepth", "taylorCenter", "selection"
         stats - an evo.trainer.RunStats object
      """
      
      stats.start("insert.traverse")
      node = self.traverse(point, config)
      stats.stop("insert.traverse")
      stats.start("insert.main")
      lr = config.learningRate
      
      ns = [i+0. for i in xrange(self.segment.taylorDepth)]
      fs = ns * 1
      if config.taylorDepth > 0:
         fs[0] = 1.0
         for i in xrange(self.segment.taylorDepth - 1):
            fs[i+1] *= fs[i]
      ns = array(ns)
      fs = array(fs)
      center = self.segment.taylorCenter      
      
      if node.point is None:
         node.point = point
         node.area = point.partition_node.area
         node.min_score = point.score
         node.max_score = point.score
         score = (config.minimize and -1 or 1) * point.score * lr 
         taylor = nan_to_num(score ** ns) / fs
         taylor *= node.area
         taylor *= nan_to_num(exp(score/center)) 
         node.taylor = nan_to_num(taylor)
         return
      
      
      other = node.point
      node.point = None
      
      if other.score > point.score:
         upPoint = other
         downPoint = point
         node.max_score = other.score
         node.min_score = point.score
      else:
         upPoint = point
         downPoint = other
         node.max_score = point.score
         node.min_score = other.score
      
      node.height = 1
      
      upArea = upPoint.partition_node.area
      downArea = downPoint.partition_node.area
      
      score1 = (config.minimize and -1 or 1) * upPoint.score * lr    
      score2 = (config.minimize and -1 or 1) * downPoint.score * lr    
      taylor1 = (score1 ** ns) / fs
      taylor1 *= upArea
      taylor1 *= (exp(score1/center)) 
      taylor2 = (score2 ** ns) / fs
      taylor2 *= downArea
      taylor2 *= (exp(score2/center))
      
      if self.segment.taylorDepth > 0:
         if (abs(taylor1) < 1e-300).all() and (abs(taylor2) < 1e-300).all():
            node.point = other
            node.height = 0
            node.max_score = other.score
            node.min_score = other.score
            # node has not been saved to the DB, so we don't need to rollback
            #node.save()
            #transaction.rollback()
            raise SeparationException, "Zero taylor coeffs"
      
      n1 = ScoreTreeNode(
         segment=node.segment,
         point=upPoint,
         area = upArea,
         parent = node,
         max_score = upPoint.score,
         min_score = upPoint.score,
         left = upPoint.score >= downPoint.score,
         taylor = taylor1
      )
      
      n2 = ScoreTreeNode(
         segment=node.segment,
         point=downPoint,
         area = downArea,
         parent = node,
         max_score = downPoint.score,
         min_score = downPoint.score,
         left = upPoint.score < downPoint.score,
         taylor = taylor2
      )
      node.children = [n1,n2]
      upPoint.score_node = n1
      downPoint.score_node = n2
      
      stats.stop("insert.main")
      stats.start("insert.propagate")
      self.propagateArea(node, config)
      stats.stop("insert.propagate")

   def resetTaylor(self, segment, temp, config):
      """
         Recompute the taylor coefficients on the Score tree.
         
         This should be done whenever 1/temp gets more than 0.5 away
         from the inverse of the taylor center for accuracy.
         
         Requires a loop through all points, but only needs to be done
         at a logarithmically often pace.
         
         segment - the segment to recompute
         temp - the new temperature center
         config - an evo.config.Config object with "shiftToDb", "taylorCenter", "taylorDepth"
         
         This method sets config.taylorCenter to temp when complete.
      """
      logg.info("resetTaylor: segment %s, new center %s, old center %s" % (segment, temp, config.taylorCenter))
      
      ns = [i+0. for i in xrange(segment.taylorDepth)]
      fs = ns * 1
      if len(fs) > 0:
         fs[0] = 1.0
         for i in xrange(segment.taylorDepth - 1):
            fs[i+1] *= fs[i]
      ns = array(ns)
      fs = array(fs)
      lr = config.learningRate
      center = temp
      
      next = []
      for point in segment.points:
         if point.alive and point.score_node is not None:
            next.append(point.score_node)
      
      heights = {0: next}
      height = 0
      mult = config.minimize and -1 or 1
      while len(heights) > 0:
         nodes = set(heights[height])
         for node in nodes:
            if node.point is not None:
               score = mult * lr * node.point.score
               taylor = nan_to_num(score ** ns) / fs
               taylor *= node.area
               taylor *= nan_to_num(exp(score/center)) 
               node.taylor = nan_to_num(taylor)
            else:
               node.taylor = zeros(segment.taylorDepth)
               for child in node.children:
                  node.taylor += child.taylor
            if node.parent is not None:
               if not heights.has_key(node.parent.height):
                  heights[node.parent.height] = [node.parent]
               else:
                  if node.parent not in heights[node.parent.height]:
                     heights[node.parent.height].append(node.parent)
         del heights[height]
         height += 1
       
      segment.taylorCenter = temp class AreaTree(object):
      
    root = None
    map = {}
    
    def traverse(self, area):
        if self.root is None:
           raise AreaException("Cannot traverse empty tree")
        #elif self.root.low > area or self.root.high < area: # if we want this check, it needs to get the high from the space
        #   err = "Area {0} out of bounds ({1},{2})"
        #   err = err.format(area, self.root.low, self.root.high)
        #   raise AreaException(err)
           
        next = self.root
        while len(next.children):
            if area >= next.children[0].low:
                next = next.children[0]
            else:
                next = next.children[1]
        
        return next
    
    def insert(self, partitionNode):
        if self.root is None:
            self.root = AreaTreeNode()
            self.root.val = partitionNode
            assert self.root.val is not None
            return
            
        parent = self.traverse(partitionNode.area)
        
        mid = partitionNode.area
        
        left = AreaTreeNode()
        left.left = True
        left.parent = parent
        left.high = parent.high
        left.low = mid
        
        right = AreaTreeNode()
        right.left = False
        right.parent = parent
        right.high = mid
        right.low = parent.low
        
        if parent.val.area >= partitionNode.area:
            right.val = partitionNode
            left.val = parent.val
        else:
            left.val = partitionNode
            right.val = parent.val
        
        assert left.val is not None
        assert right.val is not None
        assert left.low >= right.high
        assert left.val.area >= right.val.area    
                    
        parent.val = None
        parent.children = [left,right]
        
        self.map[left.val.id] = left
        self.map[right.val.id] = right
        
        #print "AFTER INSERT: "
        #self.printTree()
        
    def remove(self, partitionNode):        
        if not self.map.has_key(partitionNode.id):
            err = "Area tree has no node {0}".format(partitionNode.id)
            raise AreaException(err)
            
        target = self.map[partitionNode.id]
        #print "Removing ", target.id
        if target == self.root:
            raise AreaException("Cannot remove root")
        
        if target.left:
            other = target.parent.children[1]
        else:
            other = target.parent.children[0]
        
        #print "Other: ", other.id
        #print "Parent: ", target.parent.id
        
        #assert len(other.children) or other.val is not None
        
        target.parent.val = other.val
        target.parent.children = other.children
        for child in other.children:
            child.parent = target.parent
        other.val = None
        target.val = None
        if target.parent.val is not None:
            self.map[target.parent.val.id] = target.parent

        #assert len(target.parent.children) or target.parent.val is not None 
        
        del other
        del target
        del self.map[partitionNode.id]
              
        #print "AFTER REMOVE: "
        #self.printTree()
      
    def largest(self):
        if self.root is None:
            raise AreaException("Area tree has no root")
        
        next = self.root
        while len(next.children):
            #print "id: ", (next.low, next.high)
            next = next.children[0]
        
        return next.val.area
    
    def printTree(self, node=None, indent=""):
        if node is None:
            node = self.root
        print indent, node.id, node.low, node.high,
        print node.val is None and "None" or node.val.area
        children = node.children
        for child in children:
            self.printTree(child, indent + "\t")  class Point(object):
   
   def __init__(self, segment, point=None, statistics = None,
                score=0.0, alive=True, count=1):
      self.segment = segment
      self.point = point
      self.statistics = statistics
      self.score = score
      self.alive = alive
      self.score_node = None
      self.count = count
      self.partition_node = None
    
   @property
   def id(self):
      return self   
      
   def separable(self, config):
      """
         Get a representation of the point suitable for insertion in the Partition tree.
         
         Called by PartitionTree.separate()
         
         DEPRECATED, will be removed.
      """
      if self.point is not None: 
         return self.point
      elif self.bayes is not None:
         rep = self.bayes.edgeBinary()
         return rep
      elif self.binary:
         if self.segment.config.binaryPartition:
            return self.binary
         else:
            return self.binary.toArray(self.segment.config.dim)
      elif self.other:
         if config.behavioralPartition:
            return array(self.statistics)
         topology, angles = str(other).split("::")
         rnn = cpp.createFromTopology(topology)
         rnn = cpp.setBinaryAnglesAtDepth(rnn, RnnAngles.parse(angles).toCpp(),16)
         #rnn = cpp.convert(str(self.other))
         weights = cpp.getAngles(rnn)
         cpp.discard(rnn)
         del rnn
         return weights
      else:
         return None

   @classmethod
   def bulkSave(cls, segment, points, stats):
      segment.points.extend(points)
      
      sep = segment.config.separator(segment.config)
      for gp in points:
         if gp.score != gp.score:
            # don't insert NaN
            continue
         try:
            stats.start("separate")
            sep.separate(segment.partitionTree, gp)
            stats.stop("separate")
            stats.start("insert")
            segment.scoreTree.insert(gp, segment.config, stats)
            #segment.scoreTree.checkBalance()
            stats.stop("insert")
         except SeparationException:
            logg.debug("Exception when separating or inserting points")
            gp.alive = False
         except IndexError:
            segment.scoreTree.printTree(segment)
            
      #segment.scoreTree.printTree(segment)

   @classmethod
   def sampleTournament(cls, segment, temp, config):
      """
         Sample evolutionary annealing in tournament mode by walking the (balanced) score tree
      
         segment - the segment to sample 
         temp - the temperature at which to sample
         config - a evo.config.Config object containing parameters "pressure" and "learningRate"
         
         if config.shiftToDb is True, this will attempt to call sampleTournamentInDb.
      """
      current = segment.scoreTree.root
      min = int(config.minimize)
      
      while True:
         children = sorted(current.children, key=lambda c: not c.left)
      
         if len(children) == 0:
            break
         
         try:
            if children[min].max_score == children[1-min].min_score:
               p = .5
            else:
               p = 1. / (1. + ((1. - config.pressure) ** (config.learningRate / temp * (2 ** children[0].height))))
         except:
            p = 1.0
         
         if p < 1.0:
            p = (p / (1. - p))
            p *= children[min].area
            div = p + children[1-min].area
            if div > 0.0:
               p /= div
            else:
               p = 0.5
         
         rnd = random.random_sample()
         if rnd < p:
            current = children[min]
         else:
            current = children[1-min]
         
      return current.point, current.point.partition_node

   @classmethod
   def sampleProportional(cls, segment, temp, config):
      """
         Sample evolutionary annealing in proportional mode by walking the (balanced) score tree.
      
         segment - the segment to sample 
         temp - the temperature at which to sample
         config - a evo.config.Config object containing parameters "pressure" and "learningRate"
         
         if config.shiftToDb is True, this will attempt to use a stored procedure for postgres.
      """
      # choose the proper center
      # this function has better approximates to the right of the center
      center = segment.taylorCenter
      offset = 0
      depth = segment.taylorDepth
      
      current = segment.scoreTree.root
      
      children = sorted(current.children, key=lambda c: not c.left)
      ns = array([i+0. for i in xrange(depth)])
      diffs = (1./temp - 1./center) ** ns
      while len(children) > 0:
         # build choice
         p1 = (children[0].taylor[offset:offset+depth] * diffs).sum()
         p2 = (children[1].taylor[offset:offset+depth] * diffs).sum()
         if p1 > 0.0 or p2 > 0.0:
            p1 = p1 / (p1 + p2)
         else:
            p1 = .5
         idx = 1
         if random.random_sample() < p1:
            idx = 0
         current = children[idx]
         children = sorted(current.children, key=lambda c: not c.left)
      
      return current.point, current.point.partition_node

   @classmethod
   def sampleRegionalTournament(cls, segment, temp, config):
      """
         Sample regions from partition tree using a tournament at each node
      """
      current = segment.partitionTree.root
      currentArea = current.area
      children = current.children
      q = .5 + .5 * (config.pressure ** (1. / (temp*config.learningRate)))
      while len(children) > 0:
         try:
            if children[1].best_score == children[0].best_score:
               p = children[0].area / (children[0].area + children[1].area)
            else:
               if children[0].best_score > children[1].best_score:
                  p0 = q * children[0].area
                  p1 = (1. - q) * children[1].area
               else: #  children[1].best_score > children[0].best_score:
                  p0 = (1. - q) * children[0].area
                  p1 = q * children[1].area
               p = p0 / (p0 + p1)
         except:
            p = .5
         if random.random_sample() < p:
            current = children[0]
            currentArea = children[0].area
         else:
            current = children[1]
            currentArea = children[1].area
         children = current.children
      
      return current.point, current  class SeparationAlgorithm(object):
   """A algorithm that enables traversal and separation of the partition tree.
   This class abstracts the space-dependent portions of the partitioning
   logic.
   
   :param config: Configuration parameters
   :type config: :class:`Config`
   
   """
   def __init__(self, config):
      if not self.compatible(config.space):
         err = "{0} is not compatible with space {1}"
         raise ValueError(err.format(self.__class__.__name__, config.space))
      self.config = config
      
   def compatible(self, space):
      """Check whether this partitioning method is compatible with a given space
      
      :param space: The space to check
      :type space: class:`Space`
      
      """
      return True
   
   def firstPoint(self, tree, node, point):
      node.point = point
      point.partition_node = node
      tree.areaTree.insert(node)
   
   def separate(self, tree, point):
      """Insert the point into the partition tree,
         separating the current partition that contains it.
         
         :param tree: The partition tree
         :type tree: :class:`Partition`
         :param point: The point to insert
         :type point: ``config.space.type``
         :param config: A :class:`Config` object
         :type config: :class:`Config`
         :param stats: - an :class:`RunStats` object
         :type stats: :class:`RunStats`
         
      """
      self.config.stats.start("separate.traverse")
      lr = self.config.learningRate
      node, path = tree.traverse(point)
      self.config.stats.stop("separate.traverse")
      if node.point is None and node.bounds is self.config.space:
         self.firstPoint(tree, node, point)
         return
      
      self.config.stats.start("separate.main")
      self.config.stats.start("separate.prepare")
      
      other = node.point
      n1, n2 = self.split(point, node)
         
      # update the best score
      if path is not None:
         if self.config.minimize and (other is None or
                                      point.score < other.score):
            for n in path:
               n.best_score = min([point.score,n.best_score])
         elif other is None or point.score > other.score:
            for n in path:
               n.best_score = max([point.score,n.best_score])
      
      self.config.stats.stop("separate.compute")
      self.config.stats.stop("separate.main")
      
      self.config.stats.start("separate.areaTree")
      
      try:
         if n1.point is not None:
             tree.areaTree.insert(n1)
         if n2.point is not None:
             tree.areaTree.insert(n2)
         if not isinstance(node.bounds, Complement):
             tree.areaTree.remove(node)
      except AreaException:
         traceback.print_exc()
         raise
      
      self.config.stats.stop("separate.areaTree")
      self.config.stats.start("separate.propagate")
      
      # correct area in score tree
      if other is not None and other.score_node is not None:
         sn = other.score_node
         diff = sn.area - other.partition_node.area
         sn.segment.scoreTree.propagateAreaOnly(sn, self.config, diff)
      
      #print "end separate"
      self.config.stats.stop("separate.propagate")
      
   def split(self, point, node):
      """Do the actual work to partition a region.
      
      :param point: The point that resulted in the split
      :type point: ``config.space.type``
      :param node: A node from the partition tree corresponding to a region
                   of the space
      :type node: :class:`PartitionTreeNode`
      :throws: :class:`SeparationException` when separation of the node fails
      
      """
      raise NotImplementedError("Abstract Separation Algorithm; use a subclass") class VectorSeparationAlgorithm(SeparationAlgorithm):
   def compatible(self, space):
      return space.type == ndarray
   
   def cast(self, point):
      return point.point
   
   def rectify(self, pointrep, otherrep):
      try:
         if len(pointrep) == len(otherrep):
            return pointrep, otherrep
         
         pointrep = list(pointrep)
         otherrep = list(otherrep)
      
         while len(pointrep) < len(otherrep):
            pointrep.append(0)
         while len(pointrep) > len(otherrep):
            otherrep.append(0)
            
      except:
         pass
      
      return pointrep, otherrep
   
   def chooseIndex(self, pointrep, otherrep, lower, upper):
      newDiff = 0
      newIndex = 0
      for i in xrange(len(pointrep)):
         if lower[i] == upper[i]:
            continue
         diff = abs(pointrep[i] - otherrep[i])
         if diff > newDiff:
            newDiff = diff
            newIndex = i
      
      return newIndex, newDiff
   
   def makeParts(self, point, pointrep, other, otherrep, newIndex, lower, upper):
      if pointrep[newIndex] > otherrep[newIndex]:
         upPoint = point
         downPoint = other
         midPoint = (otherrep[newIndex] 
                     + (pointrep[newIndex] - otherrep[newIndex])/2.)
      else:
         upPoint = other
         downPoint = point
         midPoint = (pointrep[newIndex] 
                     + (otherrep[newIndex] - pointrep[newIndex])/2.)
      
      downUpper = upper.copy()
      downLower = lower.copy()
      downUpper[newIndex] = midPoint
      down = Hyperrectangle(downLower, downUpper)
      down.parent = other.partition_node.bounds
      down.owner = self.config.space
      
      upUpper = upper.copy()
      upLower = lower.copy()
      upLower[newIndex] = midPoint
      up = Hyperrectangle(upLower, upUpper)
      up.parent = other.partition_node.bounds
      up.owner = self.config.space
      
      return down, up, downPoint, upPoint
   
   def split(self, point, node):
      """A vector-based separation algorithm that generates
      Hyperrectangular partitions parallel to the axes.
      
      """
      other = node.point
      node.point = None
      try:
         newIndex = 0
         newDiff = 0
         midPoint = 0
         upPoint = other
         downPoint = point
         self.config.stats.stop("separate.prepare")
         self.config.stats.start("separate.separable")
         pointrep = self.cast(point)
         otherrep = self.cast(other)
         pointrep, otherrep = self.rectify(pointrep, otherrep)
         lower, upper = node.bounds.extent()
         newIndex, newDiff = self.chooseIndex(pointrep, otherrep, lower, upper)
         self.config.stats.stop("separate.separable")
         if newDiff == 0.0:
            raise SeparationException("No difference in points")
      
         self.config.stats.start("separate.compute")
         down, up, downPoint, upPoint = self.makeParts(point, pointrep,
                                                       other, otherrep,
                                                       newIndex,
                                                       lower, upper)
      
         downArea = float(down.area(index=newIndex))      
         upArea = float(up.area(index=newIndex))
      
         if upArea != upArea or downArea != downArea:
            raise SeparationException("NaN area!")
      
      except:
         node.point = other
         raise
      
      n1 = PartitionTreeNode(
         segment=node.segment,
         point=upPoint,
         area = upArea,
         bounds = up,
         index = newIndex,
         upper=None,
         lower=None,
         best_score = upPoint.score,
         score_sum = upPoint.score,
         parent = node,
         layer = 0,
         layerrep = ''
      )
    
      n2 = PartitionTreeNode(
         segment=node.segment,
         point=downPoint,
         area = downArea,
         bounds = down,
         index = newIndex,
         upper=None,
         lower=None,
         best_score = downPoint.score,
         score_sum = downPoint.score,
         parent = node,
         layer = 0,
         layerrep = ''
      )

      node.children = [n1,n2]
      upPoint.partition_node = n1
      downPoint.partition_node = n2
      
      return n1, n2 class LongestSideVectorSeparationAlgorithm(VectorSeparationAlgorithm):
   """Like :class:`VectorSeparationAlgorithm`, except that it always
   partitions the longest side. This can be helpful to guarantee that the
   longest side length of any partition hyperrectangle is reduced
   evenly.
   
   """
   def chooseIndex(self, pointrep, otherrep, lower, upper):
      newIndex = None
      newDiff = 0.0
      infinity = False
      for i in xrange(len(pointrep)):
         if pointrep[i] == otherrep[i]:
            continue
         if upper[i] == inf and lower[i] == -inf:
            return i, inf
         elif upper[i] == inf:
            if infinity:
               if -lower[i] > newDiff:
                  newDiff = -lower[i]
                  newIndex = i
            else:
               infinity = True
               newDiff = -lower[i]
               newIndex = i
         elif lower[i] == -inf:
            if infinity:
               if upper[i] > newDiff:
                  newDiff = upper[i]
                  newIndex = i
            else:
               infinity = True
               newDiff = upper[i]
               newIndex = i
         elif not infinity:
            diff = upper[i] - lower[i]
            if diff > newDiff:
               newDiff = diff
               newIndex = i
            
      if infinity:
         return newIndex, inf
      
      if newIndex is None:
         raise SeparationException("No difference between points")
      return newIndex, newDiff class BinarySeparationAlgorithm(VectorSeparationAlgorithm):
   def compatible(self, space):
      return space.type == TernaryString
   
   def rectify(self, pointrep, otherrep):
      if pointrep.length < otherrep.length:
         pointrep.known &= (1L << (pointrep.length)) - 1L
         pointrep.length = otherrep.length
      
      if otherrep.length < pointrep.length:
         otherrep.known &= (1L << (otherrep.length)) - 1L
         otherrep.length = pointrep.length
         
      return pointrep, otherrep

   def chooseIndex(self, pointrep, otherrep, lower, upper):
      if (pointrep.known & pointrep.base) == (otherrep.known & otherrep.base):
         raise SeparationException("No difference between points")
      
      diff = pointrep ^ otherrep
      
      # randomly wrap the bits, and then take the first difference
      wrap = random.randint(0,diff.length)
      maskAll = (1L << diff.length) - 1L
      maskLow = 1L << wrap
      wrapHigh = (diff.base & maskAll) // maskLow
      wrapKnownH = (diff.known & maskAll) // maskLow
      maskLow -= 1L
      wrapLow = diff.base & maskLow
      wrapKnownL = diff.known & maskLow
      rebuilt = TernaryString(wrapHigh | (wrapLow << wrap),
                              wrapKnownH | (wrapKnownL << wrap),
                              diff.length)
      mask = 1L
      for i in xrange(rebuilt.length):
         if rebuilt[i]:
            return ((i + (diff.length - wrap)) % diff.length), 1.0
      
      raise SeparationException("No difference in points")
   
   def makeParts(self, point, pointrep, other, otherrep, newIndex, lower, upper):
      if pointrep[newIndex]:
         upPoint = point
         downPoint = other
      else:
         upPoint = other
         downPoint = point
      
      mask = 1L << newIndex
      
      downSpec = TernaryString((lower.base & lower.known) & ~mask,
                               lower.known | mask,
                               lower.length)
      down = BinaryRectangle(downSpec)
      down.parent = other.partition_node.bounds
      
      upSpec = TernaryString((upper.base & upper.known) | mask,
                             upper.known | mask,
                             upper.length)
      up = BinaryRectangle(upSpec)
      up.parent = other.partition_node.bounds
      
      return down, up, downPoint, upPoint class BayesSeparationAlgorithm(VectorSeparationAlgorithm):
   """Separation algorithm for Bayesian networks based on structure."""
   def compatible(self, space):
      from pyec.distribution.bayes.net import BayesNet
      return space.type == BayesNet
   
   def rectify(self, pointrep, otherrep):
      return pointrep, otherrep
   
   def chooseIndex(self, pointrep, otherrep, lower, upper):
      pointrep.computeEdgeStatistics()
      otherrep.computeEdgeStatistics()
      for edge in pointrep.edges:
         if edge not in otherrep.edges:
            if edge not in lower and edge not in upper:
               return edge, 1.0
         
      for edge in otherrep.edges:
         if edge not in pointrep.edges:
            if edge not in lower and edge not in upper:
               return edge, 1.0
      
      raise SeparationException("No difference in points")

   def makeParts(self, point, pointrep, other, otherrep, newIndex, lower, upper):
      from pyec.distribution.bayes.space import BayesNetFixedEdges
      pointrep.computeEdgeStatistics()
      otherrep.computeEdgeStatistics()
      
      if newIndex in pointrep.edges:
         upPoint = point
         downPoint = other
      else:
         upPoint = other
         downPoint = point
      
      downEdges = scopy.copy(upper)
      downNedges = scopy.copy(lower) + [newIndex]
      down = BayesNetFixedEdges(self.config.space.space, downEdges, downNedges)
      down.parent = other.partition_node.bounds
      
      upEdges = scopy.copy(upper) + [newIndex]
      upNedges = scopy.copy(lower)
      up = BayesNetFixedEdges(self.config.space.space, upEdges, upNedges)
      up.parent = other.partition_node.bounds
      
      return down, up, downPoint, upPoint class LayeredSeparationAlgorithm(SeparationAlgorithm):
   """Used to separate spaces defined by a layer hierarchy; written
   for neural networks initially. Assumes that the space is divided
   into layered strata; see :class:`LayeredSpace`. The bottom layer
   is separated differently from the higher layers, and so the config
   should contain a property named "secondary_separator"
   that can be used to separate points in the lowest rung of the
   hierarchy.
   
   """
   def compatible(self, space):
      return isinstance(space, LayeredSpace)
   
   def firstPoint(self, tree, node, point):
      layers = self.config.space.extractLayers(point.point)
      other = None
      segment = node.segment
      parent = node.parent
      area = 1.0
      for layer in layers:
         area *= hasattr(layer, 'layerFactor') and layer.layerFactor() or 1.0
         node.bounds = layer
         node.best_score = point.score
         node.score_sum = point.score
         if other is not None:
            other.bounds = Complement(node.parent.bounds, layer)
         if parent is not None:
            parent.children = [node, other]
         next = PartitionTreeNode(segment, 0, node, None, None, 1.0, None)
         other = PartitionTreeNode(segment, 0, node, None, None, 0.0, None)
         parent = node
         node = next
      node = next.parent
      node.bounds.owner = node.bounds
      node.point = point
      node.area = area
      point.partition_node = node
      tree.areaTree.insert(node)
      while node.parent is not None:
         node = node.parent
         node.area = area
   
   def split(self, point, node):
      if node.point is None and isinstance(node.bounds, Complement):
         # we hit a complement node
         layers = self.config.space.extractLayers(point.point)
         area = 1.0
         while node.bounds.subtrahend.__class__ != layers[0].__class__:
            layer = layers.pop(0)
            area *= hasattr(layer, 'layerFactor') and layer.layerFactor() or 1.0
         segment = node.segment
         parent = node
         node = PartitionTreeNode(segment, 0, node, None, None, 0.0, None)
         other = PartitionTreeNode(segment, 0, node, None, None, 0.0, None)
         node.best_score = point.score
         node.score_sum = point.score
         area = 1.0
         for layer in layers:
            area *= hasattr(layer, 'layerFactor') and layer.layerFactor() or 1.0
            node.bounds = layer
            node.best_score = point.score
            node.score_sum = point.score
            other.bounds = Complement(node.parent.bounds, layer)
            parent.children = [node, other]
            next = PartitionTreeNode(segment, 0, node, None, None, 0.0, None)
            other = PartitionTreeNode(segment, 0, node, None, None, 0.0, None)
            parent = node
            node = next
         node = parent
         node.bounds.owner = node.bounds
         node.point = point
         node.area = area
         point.partition_node = node
         current = node
         while current.parent is not None:
            current = current.parent
            current.area += area
         return node, other
      else:
         cfg = self.config.merge(Config(space=node.bounds.owner))
         sep = self.config.secondary_separator(cfg)
         # problem: point.point is LayeredRnnGenotype;
         # VectorSeparationAlgorithm expects np.ndarray for Euclidean
         # need to remap point, here and in traverse
         # or create a space wrapper
         pt = point.point
         other = node.point
         otherPt = other.point
         point.point = self.config.space.layers(pt)[-1]
         other.point = self.config.space.layers(otherPt)[-1]
         bounds = node.bounds
         node.bounds = bounds.wrapped
         try:
            n1,n2 = sep.split(point, node)
            n1.bounds = self.config.space.wrapLayer(n1.bounds)
            n2.bounds = self.config.space.wrapLayer(n2.bounds)
         except:
            node.point = other
            raise
         finally:
            node.bounds = bounds
            point.point = pt
            other.point = otherPt
         return n1, n2 class RunStats(object):
   """
      A simple stats recording tool that can be used to aggregate the
      amount of time spent in various methods of an optimizer. Keeps track 
      of variables by key names and outputs the time spent between `start` 
      and `stop` for each key. For a recorded key, the average time spent 
      between `start` and `stop` can be retrieved by [], like so::
      
         def checkStats():     
            stats = RunStats()
            # start the timer
            stats.start("test")
            .. .
            # stop the timer
            stats.stop("test")
            # print just the key "test"
            print stats["test"]
            # print all
            print stats


   """
   totals = {}
   times = {}
   counts = {}
   recording = True
   
   def start(self, key):
      """
         Start recording time for `key`.
         
         :param key: A name for the key.
         :type key: str
      """
      if not self.recording: return
      if not self.totals.has_key(key):
         self.totals[key] = 0.0
         self.counts[key] = 0
      self.times[key] = time()

   def stop(self, key):
      """
         Stop recording time for `key`.
         
         :param key: The previously started key that is to be stopped.
         :type key: str
      """
      if not self.recording: return
      now = time()
      self.totals[key] += now - self.times[key]
      del self.times[key]
      self.counts[key] += 1
      
   def __getitem__(self, key):
      return self.totals[key] / self.counts[key]
      
   def __str__(self):
      ret = ""
      for key,val in sorted(self.totals.items(), key=lambda x: x[0]):
         ret += "%s: %.9f\n" % (key, self[key])
      return ret
import numpy as np
import logging
from pyec.config import Config
from pyec.distribution.convolution import Convolution
from pyec.distribution.bayes.mutators import *
from pyec.distribution.bayes.structure.proposal import StructureProposal
from pyec.distribution.ec.selectors import Selection
from pyec.distribution.ec.mutators import (Mutation,
                                           Bernoulli,
                                           Crossover,
                                           Gaussian,
                                           AreaSensitiveGaussian,
                                           AreaSensitiveBernoulli)
from pyec.distribution.nn.mutators import (AddChainLayerMutation,
                                           AddRadialLayerMutation,
                                           RemoveLayerMutation,
                                           AddNodeMutation,
                                           RemoveNodeMutation,
                                           AddLinkMutation,
                                           RemoveLinkMutation,
                                           NetAreaStripper,
                                           AreaSensitiveGaussianWeightMutation,
                                           UniformOrIntermediateRnnCrosser,
                                           UniformRnnLinkCrosser)
from pyec.history import History
from pyec.space import Euclidean, Binary
from pyec.util.combinatorics import factorial
from pyec.util.partitions import (Segment,
                                  Partition,
                                  ScoreTree,
                                  AreaTree,
                                  Point,
                                  SeparationAlgorithm,
                                  VectorSeparationAlgorithm,
                                  LongestSideVectorSeparationAlgorithm,
                                  BinarySeparationAlgorithm,
                                  BayesSeparationAlgorithm,
                                  LayeredSeparationAlgorithm)
from pyec.util.RunStats import RunStats
from scipy.special import erf
         self.config.separator = VectorSeparationAlgorithm
   
   def compatible(self, history):
      return isinstance(history, PartitionHistory)
   
   def sample(self):
      """
         Child classes should override this method in order to select a point
         from the active segment in :class:`pyec.util.partitions.Point`.
         
         The actual return should be a partition node.
         
      """
      pass
      
   def batch(self, m):
      return [self.sample() for i in xrange(m)]


class ProportionalAnnealing(Annealing):
   """
      Proportional Annealing, as described in 
      
      Lockett, A. And Miikkulainen, R. Real-space Evolutionary Annealing (2011).
      
      See :class:`Annealing` for more details.
      
      Config parameters:
      
      * taylorCenter -- Center for Taylor approximation to annealing densities.
      * taylorDepth -- The number of coefficients to use for Taylor approximation of the annealing density.
      
   """
   config = Config(taylorCenter = 1.0,
                   taylorDepth = 10,
                   history = TaylorPartitionHistory)
   
   def compatible(self, history):
      return isinstance(history, TaylorPartitionHistory)
         
   def sample(self):
      return Point.sampleProportional(self.history.segment,
                                      self.history.temperature(),
                                      self.config)
  
  
class TournamentAnnealing(Annealing):
   """
      Tournament Annealing, as described in Chapter 11 of
      
      <http://nn.cs.utexas.edu/downloads/papers/lockett.thesis.pdf>
      
      See :class:`Annealing` for more details.
      
   """
   config = Config(pressure=0.025)
   
   def sample(self):
      return Point.sampleTournament(self.history.segment,
                                    self.history.temperature(),
                                    self.config)


class AreaStripper(Mutation):
   """:class:`Annealing` and its descendents pass along a tuple; this class
   extracts the point for use with other methods.
   
   """
   def mutate(self, x):
      return x[0].point

AnnealingCrossover = (
   ((TournamentAnnealing << AreaStripper) <<
    (((TournamentAnnealing >> 1) << AreaStripper) <<
     Crossover))
)

RealEvolutionaryAnnealing = (
   TournamentAnnealing << AreaSensitiveGaussian
)[Config(separator=LongestSideVectorSeparationAlgorithm)]

RealEvolutionaryAnnealingJogo2012 = RealEvolutionaryAnnealing[Config(
   jogo2012=True,
   schedule="log",
   separator=VectorSeparationAlgorithm
)]

CrossedRealEvolutionaryAnnealing = (
   AnnealingCrossover << Gaussian
)[Config(separator=LongestSideVectorSeparationAlgorithm)]

BinaryEvolutionaryAnnealing = (
   TournamentAnnealing << AreaStripper << Bernoulli
)[Config(separator=BinarySeparationAlgorithm)]

CrossedBinaryEvolutionaryAnnealing = (
   AnnealingCrossover << Bernoulli
)[Config(separator=BinarySeparationAlgorithm)]

BayesEvolutionaryAnnealing = (
   AnnealingCrossover << StructureProposal
)[Config(separator=BayesSeparationAlgorithm,
         crosser=Merger, #UniformBayesCrosser,
         schedule="log",
         learningRate=0.1,
         crossoverProb=.25)]

NeuroannealingCrossover = (
   ((TournamentAnnealing << NetAreaStripper) <<
    (((TournamentAnnealing >> 1) << NetAreaStripper) <<
      Crossover[Config(crosser=UniformOrIntermediateRnnCrosser, crossoverProb=.5)]))
)[Config(learningRate=.1, schedule="log")]

Neuroannealing = (
   #TournamentAnnealing << NetAreaStripper << 
   NeuroannealingCrossover <<
   AddChainLayerMutation << RemoveLayerMutation <<
   AddNodeMutation << RemoveNodeMutation <<
   AddLinkMutation << RemoveLinkMutation <<
   AreaSensitiveGaussianWeightMutation

)[Config(separator=LayeredSeparationAlgorithm,