Izhikevich_model

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

The MIT License (MIT)
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.
*/

/*  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

    For speed, I'm going to force E to be a bit shift rather than a divide.
    That means it's limited (in original Izhikevich terms) to 1, 0.25, 0.125, etc:

    dv/dt = ((v * v) >> E) + F*v + G - u + I

    Since a and b are also less than one in all neuron models, these will also be bit shifts:
    a >> 5 (equal to *0.03125)
    b >> 2 (equal to *0.25)
*/

/*  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 signed int inp_I = 30; //I-value for fast spiking
const signed int inp_v = -65; //starting point for v
const signed int inp_u = 0; //starting point for u

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

/*  Izhikevich model parameters (constants for now)
    Arbitrarily scaled to low integer values */
volatile signed int para_a; //remember, this is now a bit shift: 1/a^2
volatile signed int para_abest;
volatile signed int para_b; //remember, this is now a bit shift: 1/b^2
volatile signed int para_bbest;
volatile signed int para_c;
volatile signed int para_cbest;
volatile signed int para_d;
volatile signed int para_dbest;

volatile signed int max_a = 10;
volatile signed int max_b = 5;
volatile signed int max_c = 100; //this gets flipped negative after randomization
volatile signed int max_d = 20;

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

/*  Izhikevich model coefficient limts (all minimums are 0) */
const signed int max_E = 7; //remember, this is now a bit shift: 1/E^2
const signed int max_F = 10;
const signed int max_G = 200;
const signed int max_H = 200;

/*  Slope threshold */
volatile signed int model_slopeThresh = 5;

/*  Calculated outputs */
volatile float model_numFirings;
volatile float model_numFiringsBest;
volatile float model_fireRatio;
volatile float model_fireRatioBest;
volatile float model_linSlopeAvg;
volatile float model_linSlopeAvgBest;
volatile float model_linRangePct;
volatile float model_linRangePctBest;

/*  Ideal values */
const float ideal_numFirings = 12;
const float ideal_fireRatio = 0.136;
const float ideal_linSlopeAvg = 0.252;
const float ideal_linRangePct = 0.180;
const float min_numFirings = 5;
const float max_numFirings = 20;
const float max_fireRatio = 0.5;

/*  Optimization parameters */
const int iterations = 1000000;
const float wgt_numFirings = 0.4;
const float wgt_fireRatio = 0.2;
const float wgt_linSlopeAvg = 0.3;
const float wgt_linRangePct = 0.1;

/* Model error variables */
volatile float err_numFirings;
volatile float err_numFiringsBest;
volatile float err_fireRatio;
volatile float err_fireRatioBest;
volatile float err_linSlopeAvg;
volatile float err_linSlopeAvgBest;
volatile float err_linRangePct;
volatile float err_linRangePctBest;
volatile float totalError;
volatile float totalErrorBest = 10;

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

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

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 calcError(void) {
    err_numFirings = fabs(model_numFirings - ideal_numFirings) / ideal_numFirings;
    if (model_numFirings < min_numFirings) {
        err_numFirings += 100; //sick of lazy neurons, blow this up!
    }
    if (model_numFirings > max_numFirings) {
        err_numFirings += 100;
    }
    err_fireRatio = fabs(model_fireRatio - ideal_fireRatio) / ideal_fireRatio;
    if (model_fireRatio > max_fireRatio) {
        err_fireRatio += 10;
    }
    err_linSlopeAvg = fabs(model_linSlopeAvg - ideal_linSlopeAvg) / ideal_linSlopeAvg;
    if (model_linSlopeAvg == 0) {
        err_linSlopeAvg += 10;
    }
    err_linRangePct = fabs(model_linRangePct - ideal_linRangePct) / ideal_linRangePct;
    if (model_linRangePct == 0 ) {
        err_linRangePct += 10;
    }
    totalError = wgt_numFirings * err_numFirings + wgt_fireRatio * err_fireRatio + wgt_linSlopeAvg * err_linSlopeAvg + wgt_linRangePct * err_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 += ((float)var_v[i] - (float)var_v[i-1]);

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

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

    model_fireRatio = (float)counter_fireRatio / 500;
    model_linSlopeAvg = acc_linSlopeAvg / (500 - (float)counter_fireRatio);
    model_linRangePct = (lin_max - lin_min) / ((float)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] + ((var_v[i-1] * var_v[i-1]) >> coef_E) + coef_F * var_v[i-1] + coef_G - var_u[i-1] + inp_I;
            var_u[i] = var_u[i-1] + (((var_v[i-1] >> para_b) - var_u[i-1]) >> para_a);
            }
    }
}

void printBest(void) {
    printf("        best model  ideal       error       weight\n");
    printf("numFirings  %f", model_numFiringsBest);
    printf("    %f", ideal_numFirings);
    printf("    %f", err_numFiringsBest);
    printf("    %f\n", wgt_numFirings);
    printf("fireRatio   %f", model_fireRatioBest);
    printf("    %f", ideal_fireRatio);
    printf("    %f", err_fireRatioBest);
    printf("    %f\n", wgt_fireRatio);
    printf("linSlopeAvg %f", model_linSlopeAvgBest);
    printf("    %f", ideal_linSlopeAvg);
    printf("    %f", err_linSlopeAvgBest);
    printf("    %f\n", wgt_linSlopeAvg);
    printf("linRangePct %f", model_linRangePctBest);
    printf("    %f", ideal_linRangePct);
    printf("    %f", err_linRangePctBest);
    printf("    %f\n", wgt_linRangePct);
    printf("model fit error: %f\n", totalErrorBest);
}

void updateBest(void) {
    coef_Ebest = coef_E;
    coef_Fbest = coef_F;
    coef_Gbest = coef_G;
    coef_Hbest = coef_H;
    model_numFiringsBest = model_numFirings;
    err_numFiringsBest = err_numFirings;
    model_fireRatioBest = model_fireRatio;
    err_fireRatioBest = err_fireRatio;
    model_linSlopeAvgBest = model_linSlopeAvg;
    err_linSlopeAvgBest = err_linSlopeAvg;
    model_linRangePctBest = model_linRangePct;
    err_linRangePctBest = err_linRangePct;
    totalErrorBest = totalError;
}

void printVU(void) {
    int i;
    printf("t       v(t)        u(t)\n");
    for(i=0;i<500;i++) {
        printf("%i", i);
        printf("        %i", var_v[i]);
        printf("        %i\n", var_u[i]);
    }
}

void printOutputs(void) {
    printf("a = %i\n",para_a);
    printf("b = %i\n",para_b);
    printf("c = %i\n",para_c);
    printf("d = %i\n",para_d);
    printf("I = %i\n",inp_I);
    printf("E = %i\n",coef_Ebest);
    printf("F = %i\n",coef_Fbest);
    printf("G = %i\n",coef_Gbest);
    printf("H = %i\n",coef_Hbest);
}

int main(void) {
    int i;
    start = clock();
    printf("Izhikevich model optimizer \n");
    for (i=0;i<iterations;i++) {
        para_a = 5;
        para_b = 2;
        para_c = -65;
        para_d = 8;
        coef_E = randomInt(max_E);
        coef_F = randomInt(max_F);
        coef_G = randomInt(max_G);
        coef_H = randomInt(max_H);
        runModel();
        calcOutputs();
        calcError();
        if (totalError < totalErrorBest) {
            updateBest();
        }
    }
    printOutputs();
    printBest();
    end = clock();
    printf("iterations: %i\n", iterations);
    printf("execution time (s): %f\n", ((double)(end-start))/CLOCKS_PER_SEC);
    return 0;
}