Izhikevich_test

/*
Izhikevich in integer math
Zach Fredin
zach@neurotinker.com
1/18/2016

Copyright (c) 2016 by Zach Fredin

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.


This program attemts to convert Izhikevich's neuron model to purely
integer math by iterating through test values and improving the
model "fit" versus his original equations. I'm not sure it will work.

[test.c is the baseline using floating point math!]
*/

/*  Izhikevich model

    dv/dt = 0.04v^2 + 5v +140 - u + I
    du/dt = a(bv - u)
    if v >= 30 (neuron fires)
    v <- c; u <- u + d

    Typical parameters for a regular spiking neuron:
    a   0.02
    b   0.2
    c   -65
    d   8
    I   6

    Model start at t=0
    v0 = -65
    u0 = 0

    Arbitrary names for the coefficients:
    dv/dt = E*v^2 + F*v + G - u + I
    du/dt = a(bv - u)
    if v >= H (neuron fires)
    v <- c; u <- u + d
*/

/*  Measuring model "fit"
    Various checks based on typical parameters listed above
    Check conditions
        0 <= t < 500
        slope threshold = 1.0*t

    Number of firings       6
    Fire ratio          0.132
    Linear region avg slope     0.252
    Linear region range pct     0.180

    Optimize for:
    - fire ratio at a given slope
    - number of firings
    - linear region avg slope
    - linear region range pct
    - not overflowing 16-bit signed integers
*/

#include <stdio.h>
#include <time.h>
#include <tgmath.h>
#include <stdlib.h>

/*  Izhikevich model inputs */
const float inp_I = 6;
const float inp_v = -65; //starting point for v
const float inp_u = 0; //starting point for u

/*  Izhikevich model variable array */
volatile float var_v[500];
volatile float var_u[500];

/*  Izhikevich model parameters (constants for now) */
const float para_a = 0.02;
const float para_b = 0.2;
const float para_c = -65;
const float para_d = 8;

/*  Izhikevich model coefficients (iteratively optimized) */
volatile float coef_E; //v^2
volatile float coef_Ebest;
volatile float coef_F; //v
volatile float coef_Fbest;
volatile float coef_G; //offset
volatile float coef_Gbest;
volatile float coef_H; //threshold to reset v and u
volatile float coef_Hbest;

/*  Izhikevich model coefficient limts (all minimums are 0) */
const float max_E = 0.20;
const float max_F = 20;
const float max_G = 500;
const float max_H = 200;

/*  Slope threshold */
volatile float model_slopeThresh = 1.0;

/*  Calculated outputs */
volatile int model_numFirings;
volatile float model_fireRatio;
volatile float model_linSlopeAvg;
volatile float model_linRangePct;

/*  Ideal values */
const int ideal_numFirings = 6;
const float ideal_fireRatio = 0.136;
const float ideal_linSlopeAvg = 0.252;
const float ideal_linRangePct = 0.180;

/*  Optimization parameters */
const int iterations = 1000000;
const float wgt_numFirings = 0.25;
const float wgt_fireRatio = 0.25;
const float wgt_linSlopeAvg = 0.25;
const float wgt_linRangePct = 0.25;

/*  Model fit variable (current and best) */
volatile float fit = 0;
volatile float fitbest = 0;

/*  Execution time variables */
clock_t start, end;
double cpu_time_used;

double randomDouble(double scale) {
    srand((double)clock());
    return ((double)rand()/(double)(RAND_MAX)) * scale;
}

void initializeOutputs(void) {
    int i;
    model_numFirings = 0;
    model_fireRatio = 0;
    model_linSlopeAvg = 0;
    model_linRangePct = 0;
}

void calcFit(void) {
    float fit_numFirings = wgt_numFirings * fabs((float)model_numFirings - (float)ideal_numFirings) / (float)ideal_numFirings;
    float fit_fireRatio = wgt_fireRatio * fabs(model_fireRatio - ideal_fireRatio) / ideal_fireRatio;
    float fit_linSlopeAvg = wgt_linSlopeAvg * fabs(model_linSlopeAvg - ideal_linSlopeAvg) / ideal_linSlopeAvg;
    float fit_linRangePct = wgt_linRangePct * fabs(model_linRangePct - ideal_linRangePct) / ideal_linRangePct;
    fit = 1 - (fit_numFirings + fit_fireRatio + fit_linSlopeAvg + fit_linRangePct);
}

void calcOutputs(void) {
    initializeOutputs();
    int i;
    int counter_fireRatio = 0;
    float acc_linSlopeAvg = 0;
    float lin_min = 0;
    float lin_max = 0;
    for(i=1;i<500;i++) {
        // calculate fire ratio
        if(((var_v[i] - var_v[i-1]) > model_slopeThresh) | ((var_v[i-1] - var_v[i]) > model_slopeThresh)) {
            counter_fireRatio++;
        }

        else {
            //calculate linear region average slope
            acc_linSlopeAvg += (var_v[i] - var_v[i-1]);

            //record min/max for linear region
            if(lin_min == 0 | (var_v[i] < lin_min)) {
                lin_min = var_v[i];
            }
            if(lin_max == 0 | (var_v[i] > lin_max)) {
                lin_max = var_v[i];
            }
        }

        // calculate number of firings
        if(var_v[i] > coef_H) {
            model_numFirings++;
        }
    }

    model_fireRatio = (float)counter_fireRatio / 500;
    model_linSlopeAvg = acc_linSlopeAvg / (500 - (float)counter_fireRatio);
    model_linRangePct = (lin_max - lin_min) / (coef_H - lin_min);
}

void runModel(void) {
    int i;
    var_v[0] = inp_v;
    var_u[0] = inp_u;
    for(i=1;i<500;i++) {
        if(var_v[i-1] > coef_H) {
            var_v[i] = para_c;
            var_u[i] = var_u[i-1] + para_d;
        }
        else {
            var_v[i] = var_v[i-1] + coef_E * var_v[i-1] * var_v[i-1] + coef_F * var_v[i-1] + coef_G - var_u[i-1] + inp_I;
            var_u[i] = var_u[i-1] + para_a * (para_b * var_v[i-1] - var_u[i-1]);
        }
    }
}

int main(void) {
    int i;
    start = clock();
    printf("Izhikevich model optimizer \n");

    for (i=0; i<iterations; i++) {
        coef_E = randomDouble(max_E);
        coef_F = randomDouble(max_F);
        coef_G = randomDouble(max_G);
        coef_H = randomDouble(max_H);
        runModel();
        calcOutputs();
        calcFit();
        if(fit > fitbest) {
            fitbest = fit;
            coef_Ebest = coef_E;
            coef_Fbest = coef_F;
            coef_Gbest = coef_G;
            coef_Hbest = coef_H;
        }
    }

    printf("Model iterations: %i\n", iterations);
    printf("Best fit: %f\n", fitbest);
    printf("    E %f\n", coef_Ebest);
    printf("    F %f\n", coef_Fbest);
    printf("    G %f\n", coef_Gbest);
    printf("    H %f\n", coef_Hbest);

    end = clock();
    printf("execution time (s): %f\n", ((double)(end-start))/CLOCKS_PER_SEC);
    return 0;
}