欧拉在 C 语言中的神经科学方法 [关闭]

Question

我需要为我的实习估算 ODE 之后神经元的放电率。我一开始是用 python 编写代码，它运行得很好，但是为了获得更好的性能，我的主管告诉我用 C 编写相同的代码。但是，我从来没有用 C 编写过代码，所以我是一个非常初学者和我想要的文件射击率的值全为零...有人可以帮助我吗？

非常感谢

这是我的python代码：

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()

这是我在 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

代码已修复，旧代码已被注释掉。有关更改的详细信息，请参阅 cmets。

需要查看challala.txt 进行最终测试。

#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
}

---

次要和高级数字问题：

由于double 的有限精度，请意识到商1.0*time_window / num_points 可能与数学预期略有不同。这在最坏的情况下预计会是一个非常小的数量，在 2⁵³ 中可能约为 0.5 份。

然而，重复的添加会累积错误。

  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;
  }

为避免累积错误，代码可以在每次迭代时重新计算 Time[i]。

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

当然，这样的小错误通常是可以忽略的，但当num_points 足够大时可能不会。当您的优秀代码应用于更大的任务时，就会发生这种情况。

Answer 2

关于以下3个陈述

int time_window = 90;
int num_points = 1000;

double deltaT = time_window/num_points;

由于time_window 和num_points 是整数，因此除法是作为整数除法执行的。

在整数除法中，所有分数都被截断。

表达式：time_window/num_points 实际上是：

90 / 1000

得到的分数的小数点右边的所有内容都被截断，所以结果为 0

so: Time[0] + 0 结果为 0.0。

相同的（计算的）值：0.0 然后传播到整个数组

建议更改：

int time_window = 90;
int num_points = 1000;

到

double time_window = 90.0;
double num_points  = 1000.0;

关于：

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

此语句可能（或可能不会）执行所需的功能。这取决于double的大小是否与long的大小相同

建议使用：

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