gistfile1.rb

class HPLattice

  attr_reader :residues

  POPULATION_SIZE = 100
  BUMP_PENALTY    = 1000000000
  H_BLANK_SCORE   = 1
  P_BLANK_SCORE   = -1
  HH_SCORE        = -1
  HP_SCORE        = 0
  PP_SCORE        = 0

  REMUTATE_PROBABILITY = 0.9 # probability to mutate again after a mutation

  POPULATION_RECOMBINATION_OLD = 70
  POPULATION_RECOMBINATION_NEW = 30

  GENERATE_EACH_CYCLE = 300
  
  DIRECTIONS = {
    "U" => [0,0,-1],
    "D" => [0,0,1],
    "N" => [0,1,0],
    "S" => [0,-1,0],
    "W" => [-1,0,0],
    "E" => [1,0,0]
  }

  def initialize(residues)
    @residues = residues
    @@instance = self # ghetto-style singleton..
  end
  def self.instance; @@instance; end


  #class Godwin^W^W
  class Darwin
    def initialize
      residues = HPLattice.instance.residues
      @population = []
      @genome_size = residues.size - 1 # we need this many directions in the genome.
      GENERATE_EACH_CYCLE.times do
        @population << Organism.random(@genome_size)
      end
    end

    def evolve
      new_organisms = []
      POPULATION_SIZE.times do
        parent1 = @population[rand POPULATION_SIZE]
        parent2 = @population[rand POPULATION_SIZE]
        new_organisms << parent1.mutate.mate(parent2.mutate)
      end

      best_new = new_organisms.sort_by(&:score)[0...POPULATION_RECOMBINATION_NEW]
      best_old = @population.sort_by(&:score)[0...POPULATION_RECOMBINATION_OLD]

      @population = best_new + best_old

    end

    def best
      @population.sort_by{|e|e.score}.first
    end
    
  end

  
  class Organism
    attr_accessor :genome
    
    def initialize(genome)
      if genome.kind_of? Organism
        genome = genome.genome
      end
      @genome = genome || ""
    end

    # generates a random organism with genome length `size`
    def self.random(size)
      x = Organism.new("")
      size.times do
        x.genome << DIRECTIONS.keys[rand 6]
      end
      x
    end
    
    # Is this how GA's are supposed to work? It seems like a
    # bad idea for this problem, but let's see how it goes...
    def mutate
      mutationpoint = rand(@genome.size)
      new_genome = @genome
      new_genome[mutationpoint] = DIRECTIONS.keys[rand 6]
      return Organism.new(rand < REMUTATE_PROBABILITY ? mutate : new_genome) # 3/4 chance to make another mutation..
    end

    def mate(other)
      mixpoint = rand(@genome.size)
      new_genome = [self.genome.each_char.to_a[0...mixpoint],
                    other.genome.each_char.to_a[mixpoint..-1]].
        flatten
      return Organism.new(new_genome.join(''))
    end

    def score
      @score ||= HPLattice.instance.calculate_score(@genome)
    end
    
  end
  
  
  
  def calculate_score(folds)
    mx = matricize(folds, @residues)

    score = 0

    mx.lattice.each do |i,x|
      x.each do |j,y|
        y.each do |k,z|
          score += BUMP_PENALTY*(z.size-1)

          # TODO -- don't count HH bonds twice. How? Mark LatticeifiedResidue I guess.
          DIRECTIONS.each do |name, delta|
            z1 = z.first
            z2 = mx[*[i,j,k].zip(delta).map{|e|e[0]+e[1]}].first rescue nil
            if z2.nil?
              score += (z1.hydrophobic? ? H_BLANK_SCORE : P_BLANK_SCORE)
            elsif z2.hydrophobic? && z1.hydrophobic?
              unless z1.directions.include?(name) or z2.directions.include?(inverse_direction(name))
                score += HH_SCORE
                z1.directions << name
                z2.directions << inverse_direction(name)
              end
            end
          end
          
        end
      end
    end

    score
  end

  def inverse_direction(name)
    {
      "U" => "D",
      "D" => "U",
      "N" => "S",
      "S" => "N",
      "W" => "E",
      "E" => "W"
    }[name]
  end
  
  def matricize(folds, residues)
    mx = Lattice.new
    coords = [0,0,0]
    if folds.kind_of? String
      folds_arr = folds.each_char.to_a
    else
      folds_arr = folds
    end
    folds_arr << nil
    folds_arr.zip(residues.split('')).each do |z|
      fold, residue = *z
      mx[*coords] = LatticeifiedResidue.new(residue).use(fold)
      coords = mutate_coords(coords, fold)
    end
    return mx
  end

  def mutate_coords(coords, fold)
    m_arr = DIRECTIONS[fold]
    return coords.zip(m_arr).map{|e|e.inject(:+)} rescue coords
  end
  
  class Lattice
    attr_reader :lattice

    def each
      @lattice.each
    end
    
    def initialize
      @lattice = {}
    end
      
    def [](x,y,z)
      return @lattice[x][y][z] rescue nil
    end

    def []=(x,y,z, data)
      @lattice[x]       ||= {}
      @lattice[x][y]    ||= {}
      @lattice[x][y][z] ||= []
      @lattice[x][y][z] << data
    end
  end

  class LatticeifiedResidue
    attr_reader :directions
    
    def initialize(residue)
      @residue = residue
      @directions = []
    end

    # disable a direction for counting. Eg. if we form a HH bond with this
    # residue from a different residue. Also used to disable matching along
    # the backbone.
    def use(direction)
      @directions << direction
      self
    end

    def hydrophobic?
      true if "ACGILMFPTV".each_char.include?(@residue)
    end
    
  end
  
end

if __FILE__ == $0

  seq = ARGV[0] || "AAWWAAQWOJRTBJQSLQJASJFGLWE"
  
  gl = HPLattice.new(seq)

  darwin = HPLattice::Darwin.new

  1000.times do
    darwin.evolve
    b = darwin.best
    puts "#{b.genome} : #{b.score}"
  end
  
end