#include "Population.h"
#include "Parameters.h"
#include "Random.h"

#include <cstdio>
#include <algorithm>

const FloatType Population::limitProbabilityMin = static_cast<FloatType>( 0.0 );
const FloatType Population::limitProbabilityMax = static_cast<FloatType>( 1.0 );

Population::Population( unsigned numberOfGenes_, AlleleType *numberOfAlleles_ ) :
    numberOfGenes( numberOfGenes_ ),
    numberOfAlleles( numberOfAlleles_ ),
    maxAlleles( *std::max_element( numberOfAlleles_, numberOfAlleles_ + numberOfGenes_ ) ),
    data( new MatrixType( boost::extents[numberOfGenes][maxAlleles] ) ),
    convergedAlleles( numberOfGenes_ ) {
  for ( unsigned locus = 0; locus < numberOfGenes; locus++ ) {
    FloatType probability = static_cast<FloatType>( 1.0 ) / static_cast<FloatType>( numberOfAlleles[locus] );
    for ( AlleleType allele = 0; allele < numberOfAlleles[locus]; allele++ ) {
      ( *data )[locus][allele] = probability;
    }
    convergedAlleles.set( locus, probability > Parameters::getInstance()->getConvergenceThreshold() );
  }
}

Population::~Population() {
  delete data;
}

AlleleType Population::getRandomAllele( unsigned locus, Random *random ) const {
  return random->getInt( 0, numberOfAlleles[locus] - 1 );
}

AlleleType Population::getProportionateAllele( unsigned locus, Random *random ) const {
  const unsigned myNumAlleles = numberOfAlleles[locus];
  unsigned allele;
  FloatType randVal;
  FloatType total = static_cast<FloatType>( 0.0 );
  FloatType cumulative[myNumAlleles];

  // Build an array of cumulative probabilities for the alleles at this locus
  // TODO: Maybe cache these arrays...test the app with a profiler to see if that's necessary.
  for ( allele = 0; allele < myNumAlleles; allele++ ) {
    total += ( *data )[locus][allele];
    cumulative[allele] = total;
  }

  // Pick a random number and compare it to the cumulative probabilities
	randVal = random->get01() * total; // Random value between 0 and total
  for ( allele = 0; allele < myNumAlleles; allele++ )
    if ( randVal < cumulative[allele] )
      return static_cast<AlleleType>( allele );
  // It can happen that randVal is infinitesimally >= total. In that case, just return the max. allele.
  return static_cast<AlleleType>( myNumAlleles - 1 );
}

void Population::selectIndividual( AlleleType * destination ) const {
  const Parameters * const parameters = Parameters::getInstance();
	Random *random = Random::getInstance();
  for ( unsigned locus = 0; locus < numberOfGenes; locus++ ) {
		if ( random->get01() < parameters->getMutationRate() ) {
      // Do mutation
      destination[locus] = getRandomAllele( locus, random );
    } else {
      // Use probability
      destination[locus] = getProportionateAllele( locus, random );
    }
  }
}

void Population::rewardGenes( const AlleleType * phenoType ) {
  adjustProbabilities( phenoType, Parameters::getInstance()->getReward() );
}

void Population::penalizeGenes( const AlleleType * phenoType ) {
  adjustProbabilities( phenoType, -Parameters::getInstance()->getPenalty() );
}

void Population::adjustProbabilities( const AlleleType * phenoType, FloatType epsilon ) {
  const FloatType cT = Parameters::getInstance()->getConvergenceThreshold();
  for ( unsigned locus = 0; locus < numberOfGenes; locus++ ) {
    FloatType *cell = &(( *data )[ locus ][ phenoType[locus] ] );
    FloatType oldValue = *cell, newValue = oldValue + epsilon;
    if ( newValue < limitProbabilityMax && newValue > limitProbabilityMin ) {
      *cell = newValue; // Don't update probability if new value is out of range
      if ( oldValue < cT && newValue > cT )
        convergedAlleles.set( locus );
      else if ( oldValue > cT && newValue < cT )
        convergedAlleles.reset( locus );
    }
  }
}

unsigned Population::getNumberOfConvergedAlleles() const {
	return convergedAlleles.count();
}

void Population::print() const {
  std::printf( "Virtual Population :: %ux%u matrix", numberOfGenes, maxAlleles );
  for ( unsigned locus = 0; locus < numberOfGenes; locus++ ) {
    std::printf( "\n%5u |", locus );
    for ( AlleleType allele = 0; allele < numberOfAlleles[locus]; allele++ ) {
      std::printf( " %.3f ", ( *data )[locus][allele] );
    }
  }
  std::printf( "\n" );
}
