Diversity loss simulation - http://valjok.blogspot.com/2014/07/simulating-diversity-loss-mitochondrial.html

ArgParser.sc

class Args (val args: Array[String]) extends Enumeration {
	
    class OptVal extends Val {
      override def toString = "-" + super.toString
    }
	
    val nopar, silent, reportSamples = new OptVal() { // boolean options
      def apply(): Boolean = args.contains(toString)
    }
	
    val sample, maxgen = new OptVal() { // integer options
      def apply(default: Int) = { val i = args.indexOf(toString) ;  if (i == -1) default else args(i+1).toInt}
      def apply(): Int = apply(-1)
    }
	
}

Diversity3.sc

/**
 * Computing the average number of generations to looze the diversity for populations of size 3.
 * Partial sum converges to 3.58
 **/
object Diversity3 {

	// In our model it is 1/9 probability to pass from diversity 3 into diversity 1 and 2/3 into diversity 2
	// It is 1/3 to pass from diversity 2 to diversity 1 and 2/3 is prob to stay at diversity 2.
	// Expectation is prob monopolize at first gen + 2 * prob mono at second gen + 3 * reaching diversity 1 at 3rd gen and so on
	val recursiveExpectation = (1 to 100).foldLeft(0.0, (0.0, 0.0, 1.0)){case ((acc, diversity), i) =>
			val _3 = diversity._3 * 2.0/9
			val _2 = 2.0/3 * (diversity._3 + diversity._2)
			val _1 = /*diversity._1 +*/ diversity._2 * 1.0/3 + diversity._3 * 1.0/9
			//println(diversity)
		(acc + i*_1, (_1, _2, _3))
	}._1                                      //> recursiveExpectation  : Double = 3.8571428571428537
	
	val strightforwardExpectation = (1 to 100).foldLeft(0.0){case (acc, n) =>
		val prob = 3.0/4 * Math.pow(2.0/3, n) - 7.0/4 * Math.pow(2.0/9, n)
		acc + n * prob
	}                                         //> strightforwardExpectation  : Double = 3.8571428571428537


}

DiversityIndex.sc

/**https://en.wikipedia.org/wiki/Diversity_index*/
class DiversityIndex(population: GenSeq[Char]) extends Enumeration {

  type Distribution = GenMap[Char, GenSeq[Char]]
  class Val(val f: Distribution => Number) extends super.Val {	def apply() = f(population groupBy identity)}
  import scala.math._

  def normalize(spectrum: Distribution) = spectrum.values.map( i => i.length.toDouble / population.length )
	
  val richness = new Val (m => m.size)
  val entropy = new Val (m => normalize(m).map( p => -p * log10(p) / log10(2)).sum)
  val simpson = new Val (m =>normalize(m).map(p => Math.pow(p, -2)).sum)


  implicit def conver1t(v: Value) = v.asInstanceOf[Val]

  val valueList: List[Value] = values.toList // next call may produce values in different order
	
  def row: String = {
    //val applied = valueList map (v => /*v.toString + ": "+ */"%1.4e" format v())
    val applied = valueList map (v => {
      val applied = v()
      if (applied.isInstanceOf[Int]) applied + ""
        //else "%1.4e" format applied
        else "%1.2f" format applied
    })
    applied mkString " , "
	
  }
  def header: String = {
    valueList mkString " , "
  }
	
}
/*
val arg = "abc".toList
DiversityIndex.richness() // implicit is not needed

// implicit is needed
DiversityIndex.values map (v => (v, v(arg)))*/

DiversityLoss.sc

import scala.collection._

object LossOfDiversity {

  // We have a mitochondrial Eve mother (and Y-chromosome father), common for all people. This is not
  // one individual. We/ can demonstrate that a population of fixed size devenrates into single genome if
  // we start with a collection of absolutely different genes and, at every step, will produce a new genera
  // tion of genes which consists of randomly chosen individuals of previous generation. Despite all indivi
  // duals have equal change to reproduce, this simulation demonstrates that after the number of steps
  // proportional to the population size, all the diversity is lost and we have a/ population of identical
  // individuals!
  // This simulation proves that you do not need the population bottleneck.
  // Idea http://www.youtube.com/watch?v=CBC2YI4GzNM, similar one http://lex-kravetski.livejournal.com/226581.html (Rus)
	
  val random = new scala.util.Random()
	
  def main(args: Array[String]) {
    if (args.length == 0) {
      println("This program shows how single genome evicts competitor alleles in just a couple of generations from a\n"
      + "population of limited size over generations. The question was raised to confirm the Mitochondric Eve is inevitable\n"
      + "even without the genome bottleneck. Supply the population size and number of generations to simulate.\n"
      + "Usage: <population size> [-nopar] [-reportSamples] [-sample m] [-maxgen n]\n"
      + "Ex: 20 -reportSamples -maxgen 100"
      );
      System.exit(1);
    }

    val opts = new Args(args)
    val v = args(0).toInt
	
          //('a' to 'z').toList
          //(' ' to (' ' + v).toChar).toList.reverse
    val population = Array.tabulate(v){x => (' ' + x).toChar} // arrays are much faster (4 sec vs 17 sec) than lists and can be parallelized
	
    val maxgen = opts.maxgen(10000)
	
    val padding = "%"+maxgen.toString.length+"d"
	
    def row(print: => Unit) = {if (opts.sample() == -1) print}
 
    row{println("#gen => "+"generation diversity = "  + new DiversityIndex(population).header)}
	
    def experiment: Int =  {
      (1 to maxgen).toList.foldLeft(population) { (population, i) =>
	
        val diversity = new DiversityIndex(population)
        row {printf(padding, i); println ( " => " + (population mkString "") + " diversity " + " = " + diversity.row)} // no details if series of experiments
        if (diversity.richness() == 1) {
          //println("population of size " + args(0) + " converges in " + i + " generations")
          return i
        }
        //List.tabulate(population.size)(_ => population(random.nextInt(population.size)))
        population map (_ => population(random.nextInt(population.size)))
      }
      throw new Exception("Failed to converge in " + maxgen + " generations. Increase maxgen?")
    }
	
    val start = System.currentTimeMillis()

    val samples = {
      val s0: List[Int] = List.fill(opts.sample(1)) (0)
      val s1: GenSeq[Int] = if (opts.nopar()) s0 else s0.par
      s1 map (_ => experiment)
    }
	
    if (samples.length > 1) {
      println("population of size " + args(0) + " converges in " + (if (opts.reportSamples()) samples mkString ", " else "") + " generations")
      val N = samples.size
      val avg = samples.sum/N
      val variance = samples.foldLeft(0)((variance, i) => {val diff = avg-i; variance + diff*diff})
      println("average = " + avg + ", unbiased variance=" + variance/(N-1)+ ", std deviation = " + Math.sqrt(variance/N) + " in " + (System.currentTimeMillis() - start)/1000 + " sec")
    }
  }
  //main("2 -maxgen 1000 -sample 2 -silent".split(" +").toArray)
}