Euler's method for neuroscience in C

Question

I need to estimate for my internship the firing rate of neurons that follows an ODE. I code at first in python and it goes pretty well but for better performance, my supervisor told to write the same code in C. However, I never coded in C so I am a very beginner and the file in which i want to have the values for the firing rate is full of zero... Can someone help me ?

Thank you very much

So here is my python code :

import numpy as np
import matplotlib.pyplot as plt
from math import cos, sin, sqrt, pi, exp as cos, sin, sqrt, pi, exp
#parameters
eps = 0.05
f = 0.215
mu = 1.1
D = 0.001
DeltaT = 0.01
timewindow = 40
num_points = int(timewindow/DeltaT)
T = np.linspace(0, timewindow, num_points)
#signal
s = [sin(2*3.14*f*t) for t in T]
N=30000
compteur=np.zeros(num_points)
v = np.zeros((num_points,N))
samples = np.random.normal(0, 1, (num_points,N))
for i in range(1,num_points):
     for j in range(N):
        v[i,j] = v[i-1,j] + DeltaT *(-v[i-1,j]+ mu + eps*s[i-1]) + \
                sqrt(2*D*DeltaT)*samples[i,j]
        if v[i,j]>1:
            v[i,j]=0
            compteur[i]+=1/DeltaT/N


plt.plot(T,compteur)
plt.show()

and here is my "translation" in C :

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define PI 3.1415926536
float s(float x, float);
double AWGN_generator();
FILE* fopen(const char* nomDuFichier, const char* modeOuverture);
int fclose(FILE* pointeurSurFichier);

int main(int argc, char *argv[])
{
   //parameters
   double eps = 0.05;
   float f = 0.215 ;
   double mu = 1.1 ;
   double D = 0.001 ;
   int time_window = 90;
   int num_points = 1000;
   long num_neurons = 1000;
   double deltaT = time_window/num_points;
   int i ;

  //time 
   double Time[num_points] ;

   Time[0]= 0.0;

   for (i = 1 ;  i < num_points ; i++ )

      {   Time[i] = Time[i-1] + deltaT;
      }

   //opening file for saving data
   FILE* fichier = NULL;

   fichier = fopen("challala.txt", "w");

    if (fichier != NULL)
   {
       double v[num_points][num_neurons] ;
       memset(v, 0, num_points*num_neurons*sizeof(long) );
       long  compteur[num_points];
       memset( compteur, 0, num_points*sizeof(long) );

       int pos_1 ;
       int  pos_2 ;
     //Euler's method
       for (pos_1 = 1 ; pos_1 < num_points ; pos_1 ++)

      {
           for (pos_2 = 0 ; pos_2<num_neurons ; pos_2 ++)
          {
               float t = Time[pos_1-1] ;
               v[pos_1][pos_2] = v[pos_1-1][pos_2] + deltaT *(-v[pos_1-1] 
            [pos_2]+ mu + eps*s(t, f))+ sqrt(2*D*deltaT)*AWGN_generator();

                if (v[pos_1][pos_2]>1)
                   {
                      v[pos_1][pos_2]=0 ;
                      compteur[pos_1]+=1/deltaT/num_neurons ;
                   }
           }
           fprintf(fichier, "%ld",compteur[pos_1]);

      }

   fclose(fichier);
   printf("ca a marche test.txt");

  }
  else
    {
        // On affiche un message d'erreur si on veut
        printf("Impossible d'ouvrir le fichier test.txt");
    }

    return 0;
 }

  float s(float x, float f)
 {
  return sin(2*M_PI*f*x);
 }



 double AWGN_generator()
 {/* Generates additive white Gaussian Noise samples with zero mean and a 
  standard deviation of 1. */

  double temp1;
  double temp2;
  double result;
  int p;

  p = 1;

  while( p > 0 )
  {
    temp2 = ( rand() / ( (double)RAND_MAX ) ); /*  rand() function generates an
                                                   integer between 0 and  
                                                                 RAND_MAX,
                                                   which is defined in 
                                                     stdlib.h.
                                               */

   if ( temp2 == 0 )
    {// temp2 is >= (RAND_MAX / 2)
     p = 1;
    }// end if
   else
   {// temp2 is < (RAND_MAX / 2)
      p = -1;
    }// end else

  }// end while()

  temp1 = cos( ( 2.0 * (double)PI ) * rand() / ( (double)RAND_MAX ) );
  result = sqrt( -2.0 * log( temp2 ) ) * temp1;

  return result;    // return the generated random sample to the caller

 }

Answer 1

Code repaired as able with old code commented out. See comments for change details.

Need to see challala.txt for a definitive test.

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

// Why use a coarse approximation?
//#define PI 3.1415926536
#define PI 3.1415926535897932384626433832795

// Let us stick to double only.
//float s(float x, float);
double s(double x, double);

// Add void, else declaration does not check the parameters.
//double AWGN_generator();
double AWGN_generator(void);

// These should have already been declared in <stdio.h>
// FILE* fopen(const char* nomDuFichier, const char* modeOuverture);
// int fclose(FILE* pointeurSurFichier);

int main(int argc, char *argv[]) {
  //parameters
  double eps = 0.05;
  // float f = 0.215;
  double f = 0.215;
  double mu = 1.1;
  double D = 0.001;
  int time_window = 90;

  // Unclear why `int/long` used here.  size_t would be idiomatic for array sizing.
  int num_points = 1000;
  long num_neurons = 1000;

  // Avoid integer division when a FP quotinet is desired
  // double deltaT = time_window / num_points;
  double deltaT = 1.0*time_window / num_points;
  int i;

  //time
  double Time[num_points];

  Time[0] = 0.0;

  for (i = 1; i < num_points; i++)     
  {
    Time[i] = Time[i - 1] + deltaT;
  }

  //opening file for saving data
  FILE* fichier = NULL;
  fichier = fopen("challala.txt", "w");
  if (fichier != NULL) {
    double v[num_points][num_neurons];
    // zero fill use wrong type in sizeof.  
    // Avoid type in sizeof, better to use sizeof object
    // memset(v, 0, num_points * num_neurons * sizeof(long));
    memset(v, 0, sizeof v);

    // Let us use FP here.
    //long compteur[num_points];
    double compteur[num_points];
    memset(compteur, 0, sizeof compteur);

    int pos_1;
    int pos_2;
    //Euler's method
    for (pos_1 = 1; pos_1 < num_points; pos_1++)
    {
      for (pos_2 = 0; pos_2 < num_neurons; pos_2++) {
        // float t = Time[pos_1 - 1];
        double t = Time[pos_1 - 1];
        v[pos_1][pos_2] = v[pos_1 - 1][pos_2]
            + deltaT * (-v[pos_1 - 1][pos_2] + mu + eps * s(t, f))
            + sqrt(2 * D * deltaT) * AWGN_generator();

        if (v[pos_1][pos_2] > 1) {
          v[pos_1][pos_2] = 0;
          compteur[pos_1] += 1 / deltaT / num_neurons;
        }
      }
      // Change of type
      // fprintf(fichier, "%ld", compteur[pos_1]);
      fprintf(fichier, " %g", compteur[pos_1]);
    }

    fclose(fichier);
    printf("ca a marche test.txt");

  } else {

    // Was not the file another name?
    // printf("Impossible d'ouvrir le fichier test.txt");
    printf("Impossible d'ouvrir le fichier \"%s\"\n", challala.txt);
  }

  return 0;
}

//float s(float x, float f) {
double s(double x, double f) {
  // M_PI is not defined in the standard C library, although common in extensions.
  //return sin(2 * M_PI * f * x);
  return sin(2 * PI * f * x);
}

double AWGN_generator() {    
  double temp1;
  double temp2;
  double result;

  int p;
  p = 1;
  while (p > 0) {
    temp2 = (rand() / ((double) RAND_MAX));
    if (temp2 == 0) {  // temp2 is >= (RAND_MAX / 2)
      p = 1;
    }  // end if
    else {  // temp2 is < (RAND_MAX / 2)
      p = -1;
    }  // end else
  }  // end while()

  temp1 = cos((2.0 * (double) PI) * rand() / ((double) RAND_MAX));
  result = sqrt(-2.0 * log(temp2)) * temp1;
  return result;    // return the generated random sample to the caller
}

---

Minor and advanced numeric issue:

Realize that quotient 1.0*time_window / num_points maybe be a little different than mathematical expected due to finite precision of double. This is, at worst, expected to be a very small amount, maybe about 0.5 parts in 2⁵³.

Yet the repetitive additions accumulate the error.

  double deltaT = 1.0*time_window / num_points;
  int i;
  double Time[num_points];
  Time[0] = 0.0;
  for (i = 1; i < num_points; i++) {
    Time[i] = Time[i - 1] + deltaT;
  }

To avoid that accumulated error, code can re-calculate Time[i] anew on each iteration.

  double deltaT = 1.0*time_window / num_points;
  double Time[num_points];
  for (int i = 0; i < num_points; i++) {
    Time[i] = deltaT*i;
  }

Of course, such small errors are often ignorable, but mayne not when num_points is large enough. This happens when your good code is applied to ever larger tasks.

Answer 2

regarding the following 3 statements

int time_window = 90;
int num_points = 1000;

double deltaT = time_window/num_points;

Since time_window and num_points are integers, the division is performed as an integer divide.

In an integer divide, all fractions are truncated.

the expression: time_window/num_points is actually:

90 / 1000

the resulting fraction has everything right of the decimal point truncated, so the result is 0

so: Time[0] + 0 results in 0.0.

The same (calculated) value: 0.0 is then propagated thorough the whole array

suggest changing:

int time_window = 90;
int num_points = 1000;

to

double time_window = 90.0;
double num_points  = 1000.0;

regarding:

memset(v, 0, num_points*num_neurons*sizeof(long) );

this statement may (or may not) perform the desired functionality. It depends on if the size of double is the same as the size of long

Suggest using:

memset( v, 0, sizeof( v ) );