#ifndef lint
static char *rcsid = "$Header: /tmp_mnt/vida/disks/disk5/Users/terry/r/gassy/RCS/mutate.c,v 1.4 1992/10/09 06:43:50 terry Exp terry $";
#endif


#include <assert.h>
#include <math.h>
#include "types.h"
#include "mutate.h"

VOID *
no_mutation(nselections, context)
INT nselections;
CONTEXT *context;
{
    return 0;
}

#if defined(APP_INDIVIDUALS_ARE_STRINGS) && defined(APP_INDIVIDUALS_HAVE_CONSTANT_SIZE)

#if defined(APP_INDIVIDUALS_ARE_BINARY)

#if defined(APP_FORCE_MUTATION)
#define mutate_locus(genome, locus, alleles) genome[locus] = (genome[locus] == '1') ? '0' : '1'
#else
#if defined(APP_INDIVIDUALS_HAVE_EQUAL_PROB_ALLELES) || (defined(APP_ALLELE_PROBS_CHANGE) && defined(APP_ALLELE_PROBS_CHANGE_TO_EQUAL))
#define mutate_locus(genome, locus, alleles) genome[locus] = uniform(2) ? '0' : '1'
#else
#define mutate_locus(genome, locus, alleles) genome[locus] = (knuth_random() < alleles[locus].probabilities[0]) ? '0' : '1'
#endif
#endif

#else /* Non Binary. Use a function. */
extern VOID mutate_locus();
#endif

#if ! defined(APP_ONE_ALLELE_LOCI_ALLOWED)

/*
 * Mutate the population. There are no single valued loci.
 * We are passed the numebr of loci to skip until the next
 * mutation. Return this (adjusted) number.
 */
VOID *
simple_mutation(nselections, context, time_to_next_mutation)
INT nselections;
CONTEXT *context;
INT time_to_next_mutation;
{
    register INT population_index;
    register INT genome_index;
    register INT loci_in_population = context->num_mutatable_loci_in_population;
    register INT genome_size;
    
    if (time_to_next_mutation >= loci_in_population){
	return (VOID *)(time_to_next_mutation - loci_in_population);
    }
    
    population_index = (INT) 0;
    genome_index = (INT) 0;
    genome_size = context->genome_size;
    
    do {
	char original_allele;
	population_index += time_to_next_mutation / genome_size;
	genome_index += time_to_next_mutation % genome_size;
	
	if (genome_index >= genome_size){
	    genome_index -= genome_size;
	    population_index++;
	}
	
	loci_in_population -= time_to_next_mutation;
	time_to_next_mutation = when_is_next_mutation(context->mutation_prob);
	original_allele = context->population[population_index].genome[genome_index];
	mutate_locus(context->population[population_index].genome, genome_index, context->alleles);
	
	if (original_allele != context->population[population_index].genome[genome_index]){
	    context->population[population_index].evaluated = FALSE;
	}
    } while (time_to_next_mutation < loci_in_population);

    return (VOID *) (time_to_next_mutation - loci_in_population);
}
#else

#if defined(APP_INDIVIDUALS_HAVE_EQUAL_PROB_ALLELES) || (defined(APP_ALLELE_PROBS_CHANGE) && defined(APP_ALLELE_PROBS_CHANGE_TO_EQUAL))
#define locus_has_one_allele(allele) (allele.nalleles == 1)
#else
#define locus_has_one_allele(allele) (allele.possible_alleles[0] != '\0' && allele.possible_alleles[1] == '\0')
#endif

#if defined(APP_HAS_CONSTANT_NUM_OF_ONE_ALLELE_LOCI)

/*
 * Mutate the population. There are a constant number of single valued loci.
 * We are passed the numebr of loci to skip until the next
 * mutation. Return this (adjusted) number.
 */
VOID *
simple_mutation(nselections, context, time_to_next_mutation)
INT nselections;
CONTEXT *context;
INT time_to_next_mutation;
{
    register INT population_index;
    register INT genome_index;
    register INT num_mutatable_loci_in_population = context->num_mutatable_loci_in_population;
    register INT num_mutatable_loci_on_genome;

    if (time_to_next_mutation >= num_mutatable_loci_in_population){
	return (VOID *) (time_to_next_mutation - num_mutatable_loci_in_population);
    }

    population_index = (INT) 0;
    genome_index = (INT) 0;
    num_mutatable_loci_on_genome = context->num_mutatable_loci_on_genome;
    
    do {
	char original_allele;
	register INT i = (INT) 0;

	population_index += time_to_next_mutation / num_mutatable_loci_on_genome;
	genome_index += time_to_next_mutation / num_mutatable_loci_on_genome;
	
	if (genome_index >= num_mutatable_loci_on_genome){
	    genome_index -= num_mutatable_loci_on_genome;
	    population_index++;
	}
	
	while (genome_index > 0){
	    if (!locus_has_one_allele(context->alleles[i])){
		genome_index--;
	    }
	    i++;
	}

	mutatable_loci_in_population -= time_to_next_mutation;
	time_to_next_mutation = when_is_next_mutation(context->mutation_prob);
	original_allele = context->population[population_index].genome[genome_index];
	mutate_locus(context->population[population_index].genome, genome_index, context->alleles);
	
	if (original_allele != context->population[population_index].genome[genome_index]){
	    context->population[population_index].evaluated = FALSE;
	}
    } while (time_to_next_mutation < num_mutatable_loci_in_population);

    return (VOID *) (time_to_next_mutation - num_mutatable_loci_in_population);
}
#else

/*
 * Mutate the population. The number of single valued loci is not known.
 * We are passed the numebr of loci to skip until the next
 * mutation. Return this (adjusted) number.
 */
VOID *
simple_mutation(nselections, context, time_to_next_mutation)
INT nselections;
CONTEXT *context;
INT time_to_next_mutation;
{
    register INT population_index = (INT) 0;
    register INT genome_index = (INT) 0;
    register INT genome_size = context->genome_size;
    register char *genome = context->population[0].genome;
    register LOCUS *alleles = context->alleles;
    register INT population_size = context->population_size;
    
    while (1){

	while (locus_has_one_allele(alleles[genome_index])){
	    genome_index++;
	    
	    if (genome_index == genome_size){
		genome_index = (INT) 0;
		population_index++;

		if (population_index == population_size){
		    return (VOID *)time_to_next_mutation;
		}
	    }

	    if (time_to_next_mutation == 0){
		char original_allele;
		
		time_to_next_mutation = when_is_next_mutation(context->mutation_prob);
		original_allele = context->population[population_index].genome[genome_index];
		mutate_locus(context->population[population_index].genome, genome_index, context->alleles);
		
		if (original_allele != context->population[population_index].genome[genome_index]){
		    context->population[population_index].evaluated = FALSE;
		}
	    }
	    else {
		time_to_next_mutation--;
	    }
	}
    }
}
#endif
#endif
#endif

#if ! defined(APP_INDIVIDUALS_ARE_BINARY)
VOID
mutate_locus(genome, locus, alleles)
char *genome;
INT locus;
LOCUS *alleles;
{
    /* Non-binary. */
    
    register INT new_allele;
    register char present_allele = genome[locus];
    register char *possible_alleles = alleles[locus].possible_alleles;

#if defined(APP_INDIVIDUALS_HAVE_EQUAL_PROB_ALLELES) || (defined(APP_ALLELE_PROBS_CHANGE) && defined(APP_ALLELE_PROBS_CHANGE_TO_EQUAL))
    /* Allele probabilities are equal. */

#if defined(APP_FORCE_MUTATION)
    do {
	new_allele = uniform(alleles[locus].nalleles);
    } while(present_allele == possible_alleles[new_allele]);

    genome[locus] = possible_alleles[new_allele];
    assert(genome[locus] != '\0');
#else
    genome[locus] = possible_alleles[uniform(alleles[locus].nalleles)];
    assert(genome[locus] != '\0');
#endif

#else
    /* Allele probabilities are not equal. */
    DOUBLE d;
    PROBABILITY *probabilities = alleles[locus].probabilities;

#if defined(APP_FORCE_MUTATION)
    do {
	new_allele = (INT) 0;
	d = knuth_random();
	
	while(1){
	    if ((d -= probabilities[new_allele]) < 0.0){
		break;
	    }
	    
	    new_allele++;
	}
    } while (present_allele == possible_alleles[new_allele]);
    
    genome[locus] = possible_alleles[new_allele];
    assert(genome[locus] != '\0');
#else
    d = knuth_random();
    new_allele = (INT) 0;
    
    while(1){
	if ((d -= probabilities[new_allele]) < 0.0){
	    genome[locus] = possible_alleles[new_allele];
	    assert(genome[locus] != '\0');
	    break;
	}
	
	new_allele++;
    }
#endif
#endif

    return;
}
#endif

INT
when_is_next_mutation(prob)
PROBABILITY prob;
{
    DOUBLE d;

    if (prob == 0.0){
	return -1;
    }
    
    if (prob == 1.0){
	return 0;
    }

    do {
	d = knuth_random();
    } while (d == 0.0);
    
    return (INT) ceil(log(d) / log(1.0 - prob));
}