// Author: John McCullock
// Date: 04-21-06
// Description: Genetic Algorithm Example 2

#include <iostream>
#include <iomanip>
#include <cmath>
#include <string>
#include <ctime>
#include <cstdlib>

using namespace std;

const double target = 100.0;    // The number for the algorithm to find.
const int maxInputs = 50;       // Number of chomosomes in population.
const int maxEpochs = 10000;    // Arbitrary number of cycles to stop after.
const double mRate = 0.001;     // Mutation Rate.
const double mRate2 = 0.03;     // Stagnation avoidance mutation rate.

const bool showVerboseResults = false;

/* 0000 0000 0000 0000 0000 0000 0000 0000 0000

Four bits are required to represent the range of characters used:
 Digits:                Operators:
    0: 0000                 +: 1010
    1: 0001                 +: 1011
    2: 0010                 -: 1100
    3: 0011                 -: 1101
    4: 0100                 *: 1110
    5: 0101                 /: 1111
    6: 0110
    7: 0111
    8: 1000
    9: 1001
*/

const int maxSize = 36;

int digits[10][4] = {{0, 0, 0, 0},
                    {0, 0, 0, 1},
                    {0, 0, 1, 0},
                    {0, 0, 1, 1},
                    {0, 1, 0, 0},
                    {0, 1, 0, 1},
                    {0, 1, 1, 0},
                    {0, 1, 1, 1},
                    {1, 0, 0, 0},
                    {1, 0, 0, 1}};

int operatrs[6][4] = {{1, 0, 1, 0},
                      {1, 0, 1, 1},
                      {1, 1, 0, 0},
                      {1, 1, 0, 1},
                      {1, 1, 1, 0},
                      {1, 1, 1, 1}};

int inputs[maxInputs][maxSize];		// Current population.
int nextGen[maxInputs][maxSize];	// Next population.
int totals[maxInputs];				// Decoded values.
int fitness[maxInputs];				// Fitness as percentage.
bool selected[maxInputs];			// Eligible parents.
int childCount = 0;
int nextMutation;                   // For scheduling mutations.
int mutations = 0;
int NaNAvoidanceCalls = 0;
bool stopTest = false;

void geneticAlgorithm();
void prepNextEpoch();
void selection();
void mating();
int chooseParent();
int chooseParent(int parentA);
void crossover(int chromA, int chromB, int newChrom[]);
int decodeInput(int inputIndex);
void getFitness();
void breakNaN();
void mutation(int childIndex);
void getDecodedFunction();
int getRandomNumber(int low, int high);
int minimum(int intArray[]);
int maximum(int intArray[]);
void initializeChromosomes();
inline bool IsNaN(double x);

int main()
{
    srand((unsigned)time(0));
    
    geneticAlgorithm();

	return 0;
}

void geneticAlgorithm()
{
    int epoch;
    bool done = false;

    initializeChromosomes();

    epoch = 0;

    do
    {
        for(int i = 0; i <= maxInputs - 1; i++)
        {
            totals[i] = decodeInput(i);

            if(abs(totals[i] - target) <= 0.001 || epoch == maxEpochs)
            {
                done = true;
            }
        } // i

        getFitness();

        if((epoch == 0) || (showVerboseResults == true) || (done == true))
        {
            for(int i = 0; i <= maxInputs - 1; i++)
            {
                for(int j = 0; j <= maxSize - 1; j++)
                {
                    cout << inputs[i][j];
                } // j
                cout << " = " << totals[i];
                cout << "\t" << fitness[i] << "%" << endl;
            } // i
            cout << "\n";
        }

        selection();

        mating();

        prepNextEpoch();

        epoch += 1;
        //This is here simply to show the runtime status.
        cout << "Epoch: " << epoch << endl;

        if(stopTest == true)
        {
            done = true;
        }

    }while(!done);

    cout << "Done." << endl;
    if(epoch != maxEpochs)
    {
        getDecodedFunction();
    }
    cout << "Completed " << epoch << " epochs." << endl;

    cout << "Encountered " << mutations << " mutations in " 
                    << childCount << " offspring." << endl;

    cout << "NaN avoidance routine called " << NaNAvoidanceCalls 
                                        << " times." << endl;
    
    return;
}

void initializeChromosomes()
{
    int getRand, l;

    for(int i = 0; i <= maxInputs - 1; i++)
    {
        l = 0;
        for(int j = 0; j <= 8; j++)
        {
            if((j % 2) == 0)
            {
                // j is even; this is an operand
                getRand = getRandomNumber(0, 9);
                for(int k = 0; k <= 3; k++)
                {
                    inputs[i][l] = digits[getRand][k];
                    l += 1;
                } // k
            }else{
                // j is odd; this is an operator
                getRand = getRandomNumber(0, 5);
                for(int k = 0; k <= 3; k++)
                {
                    inputs[i][l] = operatrs[getRand][k];
                    l += 1;
                } // k
            }
        } // j
    } // i

    return;
}

void getFitness()
{
/*
    Renders errors into percentage of correctness.
    Lowest errors = 100%, Highest errors = 0%
*/
    int bestScore, worstScore;
    int NaNCount = 0; // Monitors fitness scores for non-real numbers

    // The worst score would be the one furthest from the Target.
    worstScore = abs(static_cast<int>(ceil(target)) - totals[maximum(totals)]);

    // The best would be the closest.
    bestScore = abs(static_cast<int>(ceil(target)) - totals[minimum(totals)]);

    // Convert to a weighted percentage.
    bestScore = worstScore - 
                abs(static_cast<int>(ceil(target)) - totals[minimum(totals)]);
    
    for(int i = 0; i <= maxInputs - 1; i++)     
    {
        if(bestScore != 0)
		{
			fitness[i] = (worstScore - (abs(target - totals[i]))) * 100 / bestScore;
		}

    } // i

    // Prepare to shake up population if all fitness scores become NaN.
    if(abs(maximum(totals) - minimum(totals)) < 1)
    {
        breakNaN();
    }

    return;

}

void breakNaN()
{
// If the entire population achieves a fitness score of NaN,
// shake up the monotany with several mutations.
    int getRand;

    for(int i = 0; i <= maxInputs - 1; i++)
    {
        getRand = getRandomNumber(0, 100);
        if(getRand <= (1 / mRate2))
        {
            mutation(i);
        }
    } // i

    NaNAvoidanceCalls += 1;
}

void selection()
{
/*
'We start out with n individuals, and will stay with n.
'To do this, pick out the most fit, mate them, then replace the least fit
'with the new offspring.  The parents will remain for the next
'population. Basically, the least fit are always being weeded out.
*/
    int rouletteSpin;

    for(int i = 0; i <= maxInputs - 1; i++)
    {
        rouletteSpin = getRandomNumber(0, 100);
        if(fitness[i] >= rouletteSpin){
            selected[i] = true;
        }else{
            selected[i] = false;
        }
    } // i

    return;

}

void mating()
{
    int pointer1, pointer2;
    int maxChild = 0;
    int canMate = 0;
    int cannotMate[maxInputs];
    int parentA, parentB;
    int newChild[maxSize];
    
    for(int i = 0; i <= maxInputs - 1; i++)
    {
        cannotMate[i] = -1;
    }

    //Determine total who can mate.
    pointer1 = 0;
    pointer2 = 0;
    
    for(int i = 0; i <= maxInputs - 1; i++)
    {
        if(selected[i] == true){
            canMate += 1;
            // Copy selected individuals to next generation.
            for(int j = 0; j <= maxSize - 1; j++)
            {
                nextGen[pointer1][j] = inputs[i][j];
            } // j
            pointer1 += 1;
        }else{
            // Cannot mate.
            cannotMate[pointer2] = i;
            pointer2 += 1;
        }
    } // i

    maxChild = maxInputs - canMate; // Total number of offspring to be created.

    if(canMate > 1 && pointer2 > 0)
    {
        for(int i = 0; i <= maxChild - 1; i++)
        {
            parentA = chooseParent();
            parentB = chooseParent(parentA);
            crossover(parentA, parentB, newChild);
            for(int j = 0; j <= maxSize - 1; j++)
            {
                nextGen[pointer1][j] = newChild[j];
            } // j

            if(childCount == nextMutation)
            {
                mutation(pointer1);
            }

            pointer1 += 1;
            childCount += 1;

            // Schedule next mutation.
            if(childCount % static_cast<int>(1 / mRate) == 0)
            {
                nextMutation = childCount + 
                               getRandomNumber(0, static_cast<int>(1.0 / mRate));
            }

        } // i
    }

    return;

}

int chooseParent()
{
    //Overloaded function, see also "chooseParent(int parentA)".
    int parent;
    bool done = false;

    do
    {
        // Randomly choose an eligible parent.
        parent = getRandomNumber(0, maxInputs - 1);
        if(selected[parent] == true)
        {
            done = true;
        }
    }while(!done);

    return parent;

}

int chooseParent(int parentA)
{
    //Overloaded function, see also "chooseParent()".
    int parent;
    bool done = false;

    do
    {
        //Randomly choose an eligible parent.
        parent = getRandomNumber(0, maxInputs - 1);
        if(parent != parentA)
        {
            if(selected[parent] == true)
            {
                done = true;
            }
        }
    }while(!done);

    return parent;

}

void crossover(int chromA, int chromB, int newChrom[])
{
// select a random gene along the length of the chromosomes and
// swap all genes after that point.
    int randomPoint;

    // We want the point to be at a logical place, so that 
    // valid operands and operators are kept intact.
    randomPoint = getRandomNumber(0, 8);
    randomPoint *= 4;

    for(int i = 0; i <= maxSize - 1; i++)
    {
        if(i < randomPoint){
            newChrom[i] = inputs[chromA][i];
        }else if(i >= randomPoint){
            newChrom[i] = inputs[chromB][i];
        }
    } // i

    return;

}

void mutation(int childIndex)
{
    int j;
    int randomPoint;
    int randomItem;

    // We want the point to be at a logical place, so that valid operands and
    // operators are kept intact.
    randomPoint = getRandomNumber(0, 8);

    if(randomPoint % 2 == 0)
    {
        // randomPoint is even; this is an operand.
        randomItem = getRandomNumber(0, 9);
    }else{
        // randomPoint is odd; this is an operator.
        randomItem = getRandomNumber(0, 5);
    }

    j = 0;
    for(int i = 0; i <= maxSize - 1; i++)
    {
        if((i > (randomPoint * 4) - 1) && (i < (randomPoint * 4) + 4))
        {
            if(randomPoint % 2 == 0)
            {
                nextGen[childIndex][i] = digits[randomItem][j];
            }else{
                nextGen[childIndex][i] = operatrs[randomItem][j];
            }
            j += 1;
        }else{
            nextGen[childIndex][i] = inputs[childIndex][i];
        }
    } // i

    mutations += 1;

}

int decodeInput(int inputIndex)
{
/* 
    Take a chromosome, decode it, evaluate it mathematically,
    and return the answer.
    Ignore the usual rules of algebraic evaluation, and
    simple go from left to right.
*/
    int pointer;
    bool done = false;
    int operatr;
    int operand;
    double retMod;
    double total = 0.0;
/*
    '0000 0000 0000 0000 0000 0000 0000 0000 0000

    '0-3:   operand
    '4-7:   operator
    '8-11:  operand
    '12-15: operator
    '16-19: operand
    '20-23: operator
    '24-27: operand
    '28-31: operator
    '32-35: operand
*/
    pointer = 0;
    // Get first operand...
    for(int i = 3; i >= 0; i--)
    {
        //Convert from binary to decimal.
        total += inputs[inputIndex][pointer] * pow(2, i);
        pointer += 1;
    } // i

    done = false;
    do
    {
        //Get next operator...
        operatr = 0;
        for(int i = 3; i >= 0; i--)
        {
            // Convert from binary to decimal.
            operatr += inputs[inputIndex][pointer] * static_cast<int>(pow(2, i));
            pointer += 1;
        } // i

        // Get next operand...
        operand = 0;
        for(int i = 3; i >= 0; i--)
        {
            //Convert from binary to decimal.
            operand += inputs[inputIndex][pointer] * static_cast<int>(pow(2, i));
            pointer += 1;
        } // i

        if(operatr == 10 || operatr == 11){
            // Addition
            total += operand;
        }else if(operatr == 12 || operatr == 13){
            // Subtraction
            total -= operand;
        }else if(operatr == 14){
            // Multiplication
            total *= operand;
        }else{
            // Division
            if(operand != 0){ // Avoid divide-by-zero errors.
                total /= operand;
            }
        }

        if(pointer >= 35){
            done = true;
        }

    }while(!done);

    // Round total to the nearest whole number.
    if(modf(total, &retMod) < 0.5){
        total = floor(total);
    }else if(modf(total, &retMod) >= 0.5){
        total = ceil(total);
    }
    
    return static_cast<int>(total);

}

void prepNextEpoch()
{
    //Copy next generation into current generation input.
    for(int i = 0; i <= maxInputs - 1; i++)
    {
        for(int j = 0; j <= maxSize - 1; j++)
        {
            inputs[i][j] = nextGen[i][j];
        } // j
    } // i

    //Reset flags for selected individuals.
    for(int i = 0; i <= maxInputs - 1; i++)
    {
        selected[i] = false;
    } // i

}

void getDecodedFunction()
{
    int pointer, item;
    
    for(int i = 0; i <= maxInputs - 1; i++)
    {
        totals[i] = decodeInput(i);
        if(abs(totals[i] - target) <= 0.001)
        {
            // Print the chromosome.
            for(int j = 0; j <= maxSize - 1; j++)
            {
                cout << inputs[i][j];
            } // j
            cout << "\n";

            //Print the decoded function.
            pointer = 0;
            for(int j = 0; j <= 8; j++)
            {
                item = 0;
                for(int k = 3; k >= 0; k--)
                {
                    item += inputs[i][pointer] * static_cast<int>(pow(2, k)); //B2D
                    pointer += 1;
                } // k
                if(item < 10)
                {
                    cout << item << " ";
                }else{
                    if(item == 10 || item == 11){
                        // Addition
                        cout << "+ ";
                    }else if(item == 12 || item == 13){
                        // Subtraction
                        cout << "- ";
                    }else if(item == 14){
                        // Multiplication
                        cout << "* ";
                    }else{
                        // Division
                        cout << "/ ";
                    }
                }
            } // j
            cout << "= " << totals[i] << endl;
        }
    } // i

    return;

}

int getRandomNumber(int low, int high)
{
    //Generate random value between low and high (inclusive).
    return (rand() % ((high + 1) - low)) + low;
}

int minimum(int intArray[])
{
// Returns an array index.
	int winner = 0;
	bool foundNewWinner;
	bool done = false;

	do
	{
        foundNewWinner = false;
        for(int i = 0; i <= maxInputs - 1; i++)
		{
			if(i != winner)				// 'Avoid self-comparison.
			{
                // The minimum has to be in relation to the Target.
                if(abs(target - intArray[i]) < abs(target - intArray[winner]))
                {
                    winner = i;
                    foundNewWinner = true;
                }
			}
		} // i

        if(foundNewWinner == false)
        {
            done = true;
        }

	}while(!done);

    return winner;

}

int maximum(int intArray[])
{
// Returns an array index.
	int winner = 0;
	bool foundNewWinner;
	bool done = false;

	do
	{
        foundNewWinner = false;
        for(int i = 0; i <= maxInputs - 1; i++)
		{
			if(i != winner)				// 'Avoid self-comparison.
			{
                // The minimum has to be in relation to the Target.
                if(abs(target - intArray[i]) > abs(target - intArray[winner]))
                {
                    winner = i;
                    foundNewWinner = true;
                }
			}
		} // i

        if(foundNewWinner == false)
        {
            done = true;
        }

	}while(!done);

    return winner;

}