examples/machine_learning/geneticalgorithm.cpp

#include <arrayfire.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "ppm_utils.h"

#define OFFSET n/100

using namespace af;

array searchSpace_x_display = 0;
array searchSpace_y_display = 0;
array searchSpace_z_display = 0;
array searchSpace;

array samplex;
array sampley;
array samplez;

int n = 386;
int nsamples = 16;


void update() {
    array result(nsamples, 1, f32);

    gfor(array i, nsamples)
    result(i) = searchSpace(samplex(i).asfloat(), sampley(i).asfloat());
    samplez = result;
}


array selectFittest() {
    //pick top 50% fittest
    array sorted = sort(samplez);
    array indices = where(samplez >= sorted.row(nsamples / 2));

    if (indices.elements() > nsamples / 2) {
        indices = indices.rows(indices.elements() - nsamples / 2, indices.elements() - 1);
    }

    return indices;
}


void reproduce() {
    //Get fittest parents
    array selection = selectFittest();
    array parentsx = samplex(selection);
    array parentsy = sampley(selection);

    //Divide selection in two
    array parentsx1 = parentsx.rows(0, parentsx.elements() / 2 - 1);
    array parentsx2 = parentsx.rows(parentsx.elements() / 2, parentsx.elements() - 1);
    array parentsy1 = parentsy.rows(0, parentsy.elements() / 2 - 1);
    array parentsy2 = parentsy.rows(parentsy.elements() / 2, parentsy.elements() - 1);

    //Get crossover points (at which bit to crossover) and construct bit masks from them
    array crossover = rand(nsamples / 4) % (log((double)n) / log(2.0));
    array lowermask = pow2(crossover.asfloat()).as(u32) - 1;
    array uppermask = INT_MAX - lowermask;

    //Create children as the cross between two parents
    array children_x1 = (parentsx1 & uppermask) + (parentsx2 & lowermask);
    array children_y1 = (parentsy1 & uppermask) + (parentsy2 & lowermask);

    array children_x2 = (parentsx2 & uppermask) + (parentsx1 & lowermask);
    array children_y2 = (parentsy2 & uppermask) + (parentsy1 & lowermask);

    //Join two new sets
    samplex = join(0, children_x1, children_x2);
    sampley = join(0, children_y1, children_y2);

    //Create mutant children
    array mutantx = samplex.copy();
    array mutanty = sampley.copy();

    //Flip a random bit to vary the gene pool
    int bits = (int)(log((double)n) / log(2.0));
    mutantx = mutantx ^ (pow2((rand(nsamples / 2) % bits).asfloat())).as(u32);
    mutanty = mutanty ^ (pow2((rand(nsamples / 2) % bits).asfloat())).as(u32);

    samplex = join(0, samplex, mutantx);
    sampley = join(0, sampley, mutanty);

    //Update the value of each sample with the new coordinates
    update();
}

void init_samples() {
    samplex = rand(nsamples) % n;
    sampley = rand(nsamples) % n;
    update();
}


void init() {
    //initialize space
    searchSpace = tile(array(seq(n / 2)), 1, n / 2) + tile(array(array(seq(n / 2)), 1, n / 2), n / 2);

    searchSpace = join(0, searchSpace, flip(searchSpace,0));
    searchSpace = join(1, searchSpace, flip(searchSpace,1));

    //initialize display data
    searchSpace_x_display = array(seq(n));
    searchSpace_y_display = array(seq(n));

    searchSpace_x_display = searchSpace_x_display.T();

    searchSpace_x_display = tile(searchSpace_x_display, n, 1);
    searchSpace_y_display = tile(searchSpace_y_display, 1, n);

    searchSpace_x_display = flat(searchSpace_x_display);
    searchSpace_y_display = flat(searchSpace_y_display);

    searchSpace_z_display = flat(searchSpace);

    //initalize searchers
    init_samples();
}


void geneticSearch(bool console) {
    init();

    float maximumo = max<float>(samplez);
    array where_mo = where(samplez == maximumo);
    printf("Current max at at (%d,%d): %f\n", samplex(where_mo).scalar<unsigned>(), sampley(where_mo).scalar<unsigned>(), maximumo);

    float truemax = max<float>(searchSpace);
    while (!console || max<float>(samplez) < truemax * 0.99) {
        //Timer
        timer time_start = timer::start();

        //Runs reproduction algorithm
        reproduce();

        float maximum = max<float>(samplez);
        array where_m = where(samplez == maximum);
        printf("Current max at at (%d,%d): %f\n", samplex(where_m).scalar<unsigned>(), sampley(where_m).scalar<unsigned>(), maximum);

        if (!console) {
            fig("sub",1,2,1); plot3(samplex, samplez + OFFSET, sampley);
            fig("sub",1,2,2); plot3(searchSpace_x_display, searchSpace_z_display, searchSpace_y_display);
            fig("draw");
            if (max<float>(samplez) > (truemax) * 0.99) {
                printf("Finished\n");
                float maximum = max<float>(samplez);
                array where_m = where(samplez == maximum);
                printf("Max found at (%d,%d): %f\n", samplex(where_m).scalar<unsigned>(), sampley(where_m).scalar<unsigned>(), maximum);
                printf("Press [ENTER] to continue or [CTRL+C] to quit...\n");
                getchar();

                //Create a new batch of samples
                init_samples();
            }
            //spin if going too fast
            else
                while (timer::stop(time_start) < 0.5);
        }
    }

    if (!console) {
        fig("sub",1,2,1); plot3(samplex, samplez + OFFSET, sampley);
        fig("sub",1,2,2); plot3(searchSpace_x_display, searchSpace_z_display, searchSpace_y_display);

        printf("Press [ENTER] to continue or [CTRL+C] to quit...\n");
        getchar();
    }
}


int main(int argc, char** argv) {

    bool console = false;
    if (argc > 2 || (argc == 2 && strcmp(argv[1], "-"))) {
        printf("usage: %s [-]\n", argv[0]);
        return -1;
    } else if (argc == 2 && argv[1][0] == '-') {
        console = true;
    }

    try {
        printf("** ArrayFire Genetic Algorithm Search Demo **\n\n");
        geneticSearch(console);

    } catch (af::exception& e) {
        fprintf(stderr, "%s\n", e.what());
    }

    if (!console) {
        printf("hit [enter]...");
        getchar();
    }
    return 0;
}