# Copyright (C) Kevin R. Coombes, 2007-2012

###
### GENALG.R
###


########################################################################
# Here are the functions to implement the genetic algorithm generically

# The first basic "mating" function, which implements crossover as 
# a single cut in the list of features. The default is not to
# shuffle the feature order, and simply make a single cut.
crossover <- function(a, b, shuffle=FALSE) {
  cutpoint <- sample(1:(length(a)-1), 1)
  if (shuffle) {
    which <- rep(FALSE, length(a))
    which[sample(1:length(a), cutpoint)] <- TRUE
    alpha <- c(a[which], b[!which])
    beta  <- c(b[which], a[!which])
  } else {
    alpha <- c(a[1:cutpoint], b[(cutpoint+1):length(a)])
    beta <- c(b[1:cutpoint], a[(cutpoint+1):length(a)])
  }
  list(a=alpha, b=beta)
}

# The second basic "mating" function. A discrete cumulative probability
# vector is supplied as input, and we convert from a uniform distribution
# to select a pair of individuals from the population.
selectPair <- function(probvec) {
  i <- 1 + sum(probvec < runif(1))
  if (i > length(probvec)) {
    i <- 1
  }
  j <- 1 + sum(probvec < runif(1))
  if (i == j) {
    j <- j + 1
  }
  if (j > length(probvec)) {
    j <- 1
  }
  list(i=i, j=j)
}

# The "culling" routine. The use of "solve" in the computation of the
# mahalanobis distance means that the covariance matrix for the selected
# features has to be nonsingular. Although we should probably test that
# directly, we have so far only run into problems when we manage to select
# the same feature twice in the same set. We work around this by "mutating"
# any individuals that suffer from this 'fatal' gene combination.
removeDuplicates <- function(arow, mf, context) {
  needs <- length(arow) - length(unique(arow))
  while (needs > 0) {
    arow <- c(unique(arow), unlist(lapply(rep(0,needs), mf, context)))
    needs <- length(arow) - length(unique(arow))
  }
  arow
}

# This routine serves two roles. First, it is a constructor function for
# GenAlg objects. You have to supply a data matrix, where we think of each
# row as an individual in the population, and the columns mean whatever you
# want them to mean. In our application, the entries are the indexes into the
# list of features. You also have to supply
#	fitfun = a fuction that will be applied to each row in order to compute
#		the fitness of that individual
#	mutfun = a function that will be applied to each individual produced
#		by mating to cause random mutations in the inherited alleles
#       context = a list containing auxiilary data to be passed throught
#               to the the mutation and fitness functions
#	pm = the per-allele probability of a mutation
#       pc = the per-mating probability of a crossover event
#	gen = an integer represnting which generation this is in the running of
#		the algorithm.
#
# When the GenAlg is initialized, we compute the fitness of all the
# individuals in the population. Thus, the GenAlg function also serves as a
# critical routine when you iterate the generate algorithm to process several
# generations.
setClass("GenAlg",
         slots = c(data="matrix",
                   fitfun="function",
                   mutfun="function",
                   p.mutation="numeric",
                   p.crossover="numeric",
                   generation="numeric",
                   fitness="numeric",
                   best.fit="numeric",
                   best.individual="matrix",
                   context="list"))

GenAlg <- function(data, fitfun, mutfun, context, pm=0.001, pc=0.50, gen=1) {
  fitness <- apply(data, 1, fitfun, context=context)
  best.fit <- max(fitness)
  cheat <- (1:length(fitness))[fitness==best.fit][1]
  b1 <- as.matrix(data[cheat,])
  new("GenAlg",
      data=data, fitfun=fitfun, mutfun=mutfun, p.mutation=pm,
      p.crossover=pc, generation=gen, fitness=fitness,
      best.fit=best.fit, best.individual=b1, context=context)
}

# This is the iterative routine that controls the overall genetic algorithm.
# Having first initialized the genetic algorithm to produce a GenAlg object,
# you feed it to this routine and get back anotheer GenAlg object that
# represents the next generation.
newGeneration <- function(ga) {
  x <- ga@data
  fit <- ga@fitness
  probs <- cumsum(fit)/sum(fit) # prob of selection given by relative fitness
  temp <- matrix(unlist(lapply(1:(nrow(x)/2), function(i, p, data) {
    ij <- selectPair(p)
    if (runif(1) <=  ga@p.crossover) {
      ab <- crossover(data[ij$i,], data[ij$j,])
    } else {
      ab <- list(a=data[ij$i,], b=data[ij$j,])
    }
  }, probs, x)), nrow=nrow(x), ncol=ncol(x), byrow=T)
  targets <- sample(1:(nrow(x)*ncol(x)), ceiling(ga@p.mutation*nrow(x)*ncol(x)))
  start.size <- dim(temp)[2]
  temp[targets] <- ga@mutfun(temp[targets], ga@context)
  temp <- as.matrix(t(apply(temp, 1, removeDuplicates, ga@mutfun, ga@context)))
  GenAlg(temp, ga@fitfun, ga@mutfun, ga@context,
         ga@p.mutation, ga@p.crossover, ga@generation+1)
}

setMethod("summary", signature(object="GenAlg"),
          function(object, ...) {
  cat(paste("An object representing generation", object@generation,
            "in a genetic algorithm.\n"))
  cat(paste("Population size:", nrow(object@data), "\n"))
  cat(paste("Mutation probability:", object@p.mutation, "\n"))
  cat(paste("Crossover probability:", object@p.crossover, "\n"))
  cat("Fitness distribution:\n")
  print(summary(object@fitness))
})

setMethod('as.data.frame', signature(x='GenAlg'),
          function(x, row.names=NULL, optional=FALSE, ...) {
  val <- data.frame(x@fitness, x@data)
  size <- dim(x@data)[2]
  colnames(val) <- c('Fitness', paste('Feature', 1:size, sep=''))
  val
})

setMethod("as.matrix", signature(x="GenAlg"),
          function(x, ...) {
	as.matrix(as.data.frame(x))
})

## old version. Idea is to count, for each pair of feature selectors,
## the number of "new" features in the second selector. Then this value
## averaged. The double loop (well, sapply in a for loop) has been
## identified as a bottleneck when trying to compute diveristy in a loop.
## tests show it can be speeded up be a factor of about 6 by vectorizing
## (matrifying?) the computation.
popDiv <- function(x) {
  N <- nrow(x)
  ndiff <- 0
  for (i in 1:(N-1)) {
    more <- sapply((i+1):N, function(j) {
      length(unique(c(x[i,], x[j,]))) - ncol(x)
    })
    ndiff <- ndiff+sum(more)
  }
  ndiff/(N*(N-1)/2)  
}

## Vectorized version.
pd <- function(x) {
  N <- nrow(x)
  Z <- matrix(0, ncol = max(x), nrow = N)
  for (I in 1:N) Z[I, x[I,]] <- 1
  W <- ncol(x) - (Z %*% t(Z))
  sum(W[upper.tri(W)])/(N * (N-1) / 2)
}

popDiversity <- function(ga) {
  pd(ga@data)
}

########################################################################
# this is a residue of a previous demonstration version of this code.
# it only used alleles taking values 0 or 1 in a binary representation
# of an approximate real number solution to a problem. The fitness function
# that acompanied it has long since vanished.
simpleMutate <- function(allele, context) { 1-allele }

########################################################################
# Now we add additional functions for a specific genetic algorithm
# that we are often interested in. Note that the fitness and mutation
# functions assume that they have access to the following variables inside
# the 'context' list:
#	dataset = the matrix of feature values for individual samples. The
#		entries in the genetic algorithm data matrix (ie, the alleles)
#		are indices into the rows of the dataset matrix.
#	gps = a logical vector classifying the columns of the dataset into
#		two distinct groups, such as cancer or normal.
#
# In order to create an instance of the class of genetic algorithms that
# uses the fitness and mutation functions defined in this section, you must
# create the "dataset" and "gps" objects, and must then also initialize
# a starting population for the gentic algorithm.


# This function uses Mahalanobis distance to compute the fitness of
# an individual.
selectionFitness <- function(arow, context) {
  maha(t(context$dataset[arow,]), context$gps, method='var')
}

# This is the mutation routine. It assumes that every feature is equally
# likely to mutate into every other feature.
selectionMutate <- function(allele, context) {
	sample(1:nrow(context$dataset),1)
}