API tools faq
paste
SpaceQuester

Untitled

SpaceQuester
Dec 17th, 2020
299
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 19.70 KB | None | 0 0
raw download clone embed print report
  1. #define _USE_MATH_DEFINES
  2. #include "math.h"
  3. #include <stdlib.h>
  4. #include <stdio.h>
  5. #include <locale.h>
  6. #include <time.h>
  7. #include <stdbool.h>
  8.  
  9. #define Node_count 10*10*2 // ������ ���� �����, ��� �� �� ����� �����, ������������ �� ���, ���������� ������
  10. #define Node_count_half Node_count / 2
  11. #define Node_wire_width (int)sqrt(Node_count_half)
  12.  
  13. #define Equations_per_node 4
  14. #define Equations_count Node_count * Equations_per_node
  15.  
  16. double f[Equations_count];
  17. double f_diff[Equations_count];
  18.  
  19. double C_m = 1; // muF/cm^2
  20. double g_K = 35; // mS/cm^2
  21. double g_Na = 40; // mS/cm^2
  22. double g_L = 0.3; // mS/cm^2
  23. double E_K = -77; // mV
  24. double E_Na = 55; // mV
  25. double E_L = -65; // mV
  26.  
  27. double g_syn; // 0.1 // 0.04 // 0.2
  28. double k_syn = 0.2; // 0.2
  29. double E_syn[Node_count];
  30.  
  31. double I_app[Node_count];
  32.  
  33. /*const double sigma_G = 0.0;
  34. const double sigma_GN = 0.0;
  35. const double sigma_N = 1.0;*/
  36.  
  37. double** A;
  38. double** B;
  39. double* C;
  40. double** sigma;
  41. double tau[Node_count][Node_count];
  42.  
  43. #define tau_min 2 // ms
  44. #define tau_max 12 // ms
  45.  
  46. #define ms_to_step 200 // (0.001 / dt) !!! Don't forget !!!
  47.  
  48. #define Max_delay tau_max * ms_to_step
  49. double V_old_array[Node_count][Max_delay];
  50.  
  51. const double Freq = 500; // Hz
  52. const double Min_magintude = -0.2; // pA
  53. const double Max_magintude = 0.2; // pA
  54. const double Duration = 0.010; // sec
  55.  
  56. double Meander_start_from_zero[Node_count];
  57. double Meander_width[Node_count];
  58. double Meander_height[Node_count];
  59. double Meander_interval[Node_count];
  60. double last_meander_end[Node_count];
  61.  
  62. double I_stim(int i, double t)
  63. {
  64. if (t < Meander_start_from_zero[i])
  65. return 0;
  66.  
  67. t -= Meander_start_from_zero[i];
  68. t = fmod(t, Meander_width[i] + Meander_interval[i]);
  69.  
  70. return t < Meander_width[i] ? Meander_height[i] : 0;
  71. }
  72.  
  73. double V(int i)
  74. {
  75. return f[i * 4];
  76. }
  77.  
  78. void SetV(int i, double value)
  79. {
  80. f[i * 4] = value;
  81. }
  82.  
  83. double m(int i)
  84. {
  85. return f[i * 4 + 1];
  86. }
  87.  
  88. void Setm(int i, double value)
  89. {
  90. f[i * 4 + 1] = value;
  91. }
  92.  
  93. double n(int i)
  94. {
  95. return f[i * 4 + 2];
  96. }
  97.  
  98. void Setn(int i, double value)
  99. {
  100. f[i * 4 + 2] = value;
  101. }
  102.  
  103. double h(int i)
  104. {
  105. return f[i * 4 + 3];
  106. }
  107.  
  108. void Seth(int i, double value)
  109. {
  110. f[i * 4 + 3] = value;
  111. }
  112.  
  113. double V_old(int i, int delay)
  114. {
  115. return V_old_array[i][Max_delay - 1 - delay];
  116. }
  117.  
  118. int RandomI(int min, int max)
  119. {
  120. return ((double)rand() / (RAND_MAX - 1)) * (max - min) + min;
  121. }
  122.  
  123. double RandomD(double min, double max)
  124. {
  125. return ((double)rand() / RAND_MAX) * (max - min) + min;
  126. }
  127.  
  128. double alpha_m(double* f, int i)
  129. {
  130. return 0.182 * (V(i) + 35) / (1 - exp(-(V(i) + 35) / 9));
  131. }
  132.  
  133. double beta_m(double* f, int i)
  134. {
  135. return -0.124 * (V(i) + 35) / (1 - exp((V(i) + 35) / 9));
  136. }
  137.  
  138. double alpha_n(double* f, int i)
  139. {
  140. return 0.02 * (V(i) - 25) / (1 - exp(-(V(i) - 25) / 9));
  141. }
  142.  
  143. double beta_n(double* f, int i)
  144. {
  145. return -0.002 * (V(i) - 25) / (1 - exp((V(i) - 25) / 9));
  146. }
  147.  
  148. double alpha_h(double* f, int i)
  149. {
  150. return 0.25 * exp(-(V(i) + 90) / 12);
  151. }
  152.  
  153. double beta_h(double* f, int i)
  154. {
  155. return 0.25 * exp((V(i) + 62) / 6) / exp((V(i) + 90) / 12);
  156. }
  157.  
  158. double HodgkinHuxley(int i, double* f, double t)
  159. {
  160. int in = i / 4;
  161. int il = i % 4;
  162.  
  163. switch (il)
  164. {
  165. case 0:
  166. {
  167. double sum = 0;
  168.  
  169. /*for (int j = 0; j < Node_count; j++)
  170. {
  171. //sum += A[in][j] * g_syn * (V(in) - V_old(j, tau[in][j]));
  172. //sum += A[in][j] * g_syn * (V(j) - V(in));
  173. //sum += A[in][j] * g_syn * (V(in) - E_syn[in]) / (1 + exp(-V_old(j, tau[in][j]) / k_syn));
  174. //sum += 1 / (0.2 * Node_count_half) * A[in][j] * g_syn * (V(in) - E_syn[in]) / (1 + exp(-V(j) / k_syn)); // i up, j down
  175. sum += 1 / (0.2 * Node_count_half) * A[in][j] * g_syn * (V(j) - E_syn[j]) / (1 + exp(-V(in) / k_syn)); // j up, i down
  176. //printf("in = %d\t Node_count = %d\t A[in][j] = %f\t V(in) = %f\t E_syn[in] = %f\t sum = %f\n", in, j, A[in][j], V(in), E_syn[in], sum);
  177. }*/
  178.  
  179. /*for (int j = 0; j < C[in]; j++)
  180. {
  181. sum += sigma[in][(int)B[in][j]] * (V((int)B[in][j]) - V(in));
  182. }*/
  183.  
  184. for (int j = 0; j < C[in]; j++)
  185. {
  186. //sum += g_syn * (V((int)B[in][j]) - V(in)); // ���������� �����, ����� ��� �������� ���������� �����
  187. //sum += A[in][j] * g_syn * (V(in) - E_syn[in]) / (1 + exp(-V_old((int)B[in][j]) / k_syn)); // ���������� �����
  188. //sum += A[in][(int)B[in][j]] * g_syn * (V(in) - E_syn[in]) / (1 + exp(-V_old((int)B[in][j], tau[in][(int)B[in][j]]) / k_syn));
  189. //sum += A[in][(int)B[in][j]] * g_syn * (V((int)B[in][j]) - E_syn[(int)B[in][j]]) / (1 + exp(-V_old(in, tau[in][(int)B[in][j]]) / k_syn));
  190. //sum += 1 / (0.2 * Node_count_half) * A[in][(int)B[in][j]] * g_syn * (V((int)B[in][j]) - E_syn[(int)B[in][j]]) / (1 + exp(-V(in) / k_syn)); // j up, i down
  191. sum += 1 / (0.2 * Node_count_half) * A[in][(int)B[in][j]] * g_syn * (E_syn[in] - V(in)) / (1 + exp(-V((int)B[in][j]) / k_syn)); // i up, j down
  192.  
  193. //Vold = V_old((int)B[in][j], tau[in][(int)B[in][j]]);
  194. //ee = exp(-V_old((int)B[in][j], tau[in][(int)B[in][j]]) / k_syn);
  195. //printf("i = %d\t V_old = %f\t exp = %f\n", in, Vold, ee);
  196. }
  197. //printf("i = %d\t sum = %f\n", in, sum);
  198.  
  199. return 1000 * ((g_Na * pow(m(in), 3) * h(in) * (E_Na - V(in)) + g_K * n(in) * (E_K - V(in)) + g_L * (E_L - V(in)) + I_app[in] /*+ I_stim(in, t)*/ + sum ) / C_m); // V
  200. }
  201.  
  202. case 1:
  203. {
  204. return 1000 * (alpha_m(f, in) * (1 - m(in)) - beta_m(f, in) * m(in)); // m
  205. }
  206.  
  207. case 2:
  208. {
  209. return 1000 * (alpha_n(f, in) * (1 - n(in)) - beta_n(f, in) * n(in)); // n
  210. }
  211.  
  212. case 3:
  213. {
  214. return 1000 * (alpha_h(f, in) * (1 - h(in)) - beta_h(f, in) * h(in)); // h
  215. }
  216. }
  217.  
  218. return 0;
  219. }
  220.  
  221. void RungeKutta(double t, double dt, double* f, double* f_next)
  222. {
  223. double k[Equations_count][4];
  224.  
  225. // k1
  226. for (int i = 0; i < Equations_count; i++)
  227. k[i][0] = HodgkinHuxley(i, f, t) * dt;
  228.  
  229. double phi_k1[Equations_count];
  230. for (int i = 0; i < Equations_count; i++)
  231. phi_k1[i] = f[i] + k[i][0] / 2;
  232.  
  233. // k2
  234. for (int i = 0; i < Equations_count; i++)
  235. k[i][1] = HodgkinHuxley(i, phi_k1, t) * dt;
  236.  
  237. double phi_k2[Equations_count];
  238. for (int i = 0; i < Equations_count; i++)
  239. phi_k2[i] = f[i] + k[i][1] / 2;
  240.  
  241. // k3
  242. for (int i = 0; i < Equations_count; i++)
  243. k[i][2] = HodgkinHuxley(i, phi_k2, t) * dt;
  244.  
  245. double phi_k3[Equations_count];
  246. for (int i = 0; i < Equations_count; i++)
  247. phi_k3[i] = f[i] + k[i][2] / 2;
  248.  
  249. // k4
  250. for (int i = 0; i < Equations_count; i++)
  251. k[i][3] = HodgkinHuxley(i, phi_k3, t) * dt;
  252.  
  253. for (int i = 0; i < Equations_count; i++)
  254. f_next[i] = f[i] + (k[i][0] + 2 * k[i][1] + 2 * k[i][2] + k[i][3]) / 6;
  255. }
  256.  
  257. void CopyArray(double* source, double* target, int N)
  258. {
  259. for (int i = 0; i < N; i++)
  260. target[i] = source[i];
  261. }
  262.  
  263. bool Approximately(double a, double b)
  264. {
  265. if (a < 0)
  266. a = -a;
  267.  
  268. if (b < 0)
  269. b = -b;
  270.  
  271. return a - b <= 0.000001;
  272. }
  273.  
  274. // http://preshing.com/20111007/how-to-generate-random-timings-for-a-poisson-process/
  275. double nextTime(double rateParameter)
  276. {
  277. return -log(1.0 - (double)rand() / (RAND_MAX)) / rateParameter;
  278. }
  279.  
  280. void GenerateRandomMeander(int i, double min_start_time)
  281. {
  282. double offset = nextTime(Freq);
  283.  
  284. if (offset < 0)
  285. {
  286. int a = 0;
  287. }
  288.  
  289. Meander_start_from_zero[i] = min_start_time + offset;
  290. Meander_width[i] = Duration;
  291. Meander_height[i] = RandomD(Min_magintude, Max_magintude);
  292. }
  293.  
  294. void FillAMatrixZero()
  295. {
  296. for (int i = 0; i < Node_count; i++)
  297. {
  298. for (int j = 0; j < Node_count; j++)
  299. {
  300. A[i][j] = 0;
  301. }
  302. }
  303. }
  304.  
  305. void FillBCMatrix()
  306. {
  307. for (int i = 0; i < Node_count; i++)
  308. {
  309. int bIndex = 0;
  310. C[i] = 0;
  311. for (int j = 0; j < Node_count; j++)
  312. {
  313. if (A[i][j] == 1)
  314. {
  315. B[i][bIndex] = j;
  316. bIndex++;
  317. C[i]++;
  318. }
  319. }
  320. }
  321. }
  322.  
  323. /*void FillSigmaMatrix()
  324. {
  325. for (int i = 0; i < Node_count; i++)
  326. {
  327. for (int j = 0; j < Node_count; j++)
  328. {
  329. if ((i >= 0 && i < Node_count_half) && (j >= 0 || j < Node_count_half))
  330. {
  331. sigma[i][j] = sigma_G;
  332. }
  333. if ((((i >= Node_count_half && i < Node_count) && (j >= 0 && j < Node_count_half)) || ((i >= 0 && i < Node_count_half) && (j >= Node_count_half && j < Node_count))))
  334. {
  335. sigma[i][j] = sigma_GN;
  336. }
  337. if ((i >= Node_count_half && i < Node_count) && (j >= Node_count_half && j < Node_count))
  338. {
  339. sigma[i][j] = sigma_N;
  340. }
  341. }
  342. }
  343. }*/
  344.  
  345. void FillRandomMatrix()
  346. {
  347. double p_links = 0.2; // 0.2
  348. //for (int k = 0; k < 2 * Node_count_half; k++)
  349. for (int k = 0; k < p_links * Node_count_half * (Node_count_half - 1); k++)
  350. {
  351. int i, j;
  352.  
  353. do
  354. {
  355. i = RandomI(Node_count_half, Node_count);
  356. j = RandomI(Node_count_half, Node_count);
  357. } while ((i == j) || (A[i][j] == 1));
  358.  
  359. A[i][j] = 1;
  360. //A[j][i] = 1;
  361. }
  362. }
  363.  
  364. void FillVOldFromCurrent()
  365. {
  366. for (int i = 0; i < Node_count; i++)
  367. for (int j = 0; j < Max_delay; j++)
  368. V_old_array[i][j] = V(i);
  369. }
  370.  
  371. void UpdateVOld()
  372. {
  373. for (int i = 0; i < Node_count; i++)
  374. {
  375. for (int j = 1; j < Max_delay; j++)
  376. V_old_array[i][j - 1] = V_old_array[i][j];
  377.  
  378. V_old_array[i][Max_delay - 1] = V(i);
  379. }
  380. }
  381.  
  382. /*void FillFullTauMatrix()
  383. {
  384. for (int i = 0; i < Node_count; i++)
  385. {
  386. for (int j = 0; j < Node_count; j++)
  387. {
  388. if (i < Node_count_half || j < Node_count_half)
  389. {
  390. tau[i][j] = 0;
  391. continue;
  392. }
  393.  
  394. int i_neuron = i - Node_count_half;
  395. int j_neuron = j - Node_count_half;
  396.  
  397. int i_wire_x = i_neuron / Node_wire_width;
  398. int i_wire_y = i_neuron % Node_wire_width;
  399.  
  400. int j_wire_x = j_neuron / Node_wire_width;
  401. int j_wire_y = j_neuron % Node_wire_width;
  402.  
  403. double distance_max = sqrt(2.) * (Node_wire_width - 1);
  404. double distance = sqrt((i_wire_x - j_wire_x) * (i_wire_x - j_wire_x) + (i_wire_y - j_wire_y) * (i_wire_y - j_wire_y));
  405.  
  406. tau[i][j] = (tau_min + distance / (distance_max) * (tau_max - tau_min)) * ms_to_step;
  407. }
  408. }
  409. }*/
  410.  
  411. void FillTauMatrix()
  412. {
  413. for (int i = 0; i < Node_count; i++)
  414. {
  415. for (int j = 0; j < Node_count; j++)
  416. {
  417. if (i < Node_count_half || j < Node_count_half || i == j || A[i][j] == 0)
  418. {
  419. tau[i][j] = 0;
  420. continue;
  421. }
  422.  
  423. int i_neuron = i - Node_count_half;
  424. int j_neuron = j - Node_count_half;
  425.  
  426. int i_wire_x = i_neuron / Node_wire_width;
  427. int i_wire_y = i_neuron % Node_wire_width;
  428.  
  429. int j_wire_x = j_neuron / Node_wire_width;
  430. int j_wire_y = j_neuron % Node_wire_width;
  431.  
  432. double distance_max = sqrt(2.) * (Node_wire_width - 1);
  433. double distance = sqrt((i_wire_x - j_wire_x) * (i_wire_x - j_wire_x) + (i_wire_y - j_wire_y) * (i_wire_y - j_wire_y));
  434.  
  435. double t = (distance - 1) / (distance_max - 1);
  436. tau[i][j] = (tau_min + t * (tau_max - tau_min)) * ms_to_step;
  437. }
  438. }
  439. }
  440.  
  441. int main(int argc, char *argv[])
  442. {
  443.  
  444. sscanf(argv[1], "%lf", &g_syn);
  445.  
  446. FILE *fp_g_syn;
  447. fp_g_syn = fopen("g_syn.txt", "w");
  448. fprintf(fp_g_syn, "%f\t", g_syn);
  449. fclose(fp_g_syn);
  450.  
  451. FILE *fp_g_syn_omega;
  452. fp_g_syn_omega = fopen("g_syn_omega.txt", "a");
  453. fprintf(fp_g_syn_omega, "%f\t", g_syn);
  454. fclose(fp_g_syn_omega);
  455.  
  456. double g_syn_real;
  457. g_syn_real = 1 / (0.2 * Node_count_half) * g_syn;
  458. printf("g_syn_real = %f\n", g_syn_real);
  459.  
  460. FILE *fp0;
  461. //FILE *fp_Ca;
  462. //FILE *fp_IP3;
  463. //FILE *fp_z;
  464. //FILE *fp_G;
  465. FILE *fp_I_stim;
  466. FILE *fp_V;
  467. FILE *fp_m;
  468. FILE *fp_n;
  469. FILE *fp_h;
  470. FILE *fp_V_spikes;
  471. srand(time(NULL));
  472.  
  473. //for (int i = 0; i < Node_count; i++)
  474. // V_old_length[i] = 0;
  475.  
  476. A = malloc(Node_count * sizeof(double));
  477. for (int i = 0; i < Node_count; i++)
  478. A[i] = malloc(Node_count * sizeof(double));
  479.  
  480. B = malloc(Node_count * sizeof(double));
  481. for (int i = 0; i < Node_count; i++)
  482. B[i] = malloc(Node_count * sizeof(double));
  483.  
  484. C = malloc(Node_count * sizeof(double));
  485.  
  486. sigma = malloc(Node_count * sizeof(double));
  487. for (int i = 0; i < Node_count; i++)
  488. sigma[i] = malloc(Node_count * sizeof(double));
  489.  
  490. FillAMatrixZero();
  491. FillRandomMatrix();
  492. FillBCMatrix();
  493. //FillSigmaMatrix();
  494. FillTauMatrix();
  495.  
  496. fp0 = fopen("A.txt", "w+");
  497. for (int i = 0; i < Node_count; i++)
  498. {
  499. for (int j = 0; j < Node_count; j++)
  500. {
  501. fprintf(fp0, "%d\t", (int)A[i][j]);
  502. }
  503. fprintf(fp0, "\n");
  504. }
  505. fclose(fp0);
  506.  
  507. fp0 = fopen("tau.txt", "w+");
  508. for (int i = 0; i < Node_count; i++)
  509. {
  510. for (int j = 0; j < Node_count; j++)
  511. {
  512. fprintf(fp0, "%f\t", tau[i][j] / ms_to_step);
  513. }
  514. fprintf(fp0, "\n");
  515. }
  516. fclose(fp0);
  517.  
  518. //����� � ���� ����� ������ � ������� ����������� � ��������� ���������
  519. fp0 = fopen("links.txt", "w+");
  520. for (int i = Node_count_half; i < Node_count; i++)
  521. {
  522. int links_count = 0;
  523. for (int j = Node_count_half; j < Node_count; j++)
  524. {
  525. if (A[i][j] == 1)
  526. {
  527. links_count++;
  528. }
  529. }
  530. fprintf(fp0, "%d\n", (int)links_count);
  531.  
  532. }
  533. fclose(fp0);
  534.  
  535. //setlocale(LC_NUMERIC, "French_Canada.1252");
  536. fp0 = fopen("test_Poisson.txt", "w+");
  537. for (int i = 0; i < 1000; i++)
  538. fprintf(fp0, "%f\n", nextTime(Freq));
  539. fclose(fp0);
  540.  
  541. fp0 = fopen("B.txt", "w+");
  542. for (int i = 0; i < Node_count; i++)
  543. {
  544. for (int j = 0; j < C[i]; j++)
  545. {
  546. fprintf(fp0, "%d\t", (int)B[i][j]);
  547. }
  548. fprintf(fp0, "\n");
  549. }
  550. fclose(fp0);
  551.  
  552. fp0 = fopen("C.txt", "w+");
  553. for (int i = 0; i < Node_count; i++)
  554. {
  555. fprintf(fp0, "%d\n", (int)C[i]);
  556. }
  557. fclose(fp0);
  558.  
  559. /*fp0 = fopen("sigma.txt", "w+");
  560. for (int i = 0; i < Node_count; i++)
  561. {
  562. for (int j = 0; j < Node_count; j++)
  563. {
  564. fprintf(fp0, "%f\t", sigma[i][j]);
  565. }
  566. fprintf(fp0, "\n");
  567. }
  568. fclose(fp0);*/
  569.  
  570. // Initial values
  571. /*for (int i = 0; i < Equations_count; i++)
  572. {
  573. f[i] = 0;
  574. }
  575.  
  576. for (int i = 0; i < Equations_count; i++)
  577. {
  578. I_app[i] = RandomD(9, 40);
  579. }*/
  580.  
  581. double percent_stable_state = 0.40; // 0.40
  582. double eps_persent = 0.05;
  583.  
  584. double V0 = -58.7085;
  585. double m0 = 0.0953;
  586. double n0 = 0.000913;
  587. double h0 = 0.3662;
  588.  
  589. double V1 = 14.8409;
  590. double m1 = 0.9174;
  591. double n1 = 0.0140;
  592. double h1 = 0.0539;
  593.  
  594. for (int i = 0; i < Node_count; i++) // init array for all nodes
  595. {
  596. SetV(i, 0); // V
  597. Setm(i, 0); // m
  598. Setn(i, 0); // n
  599. Seth(i, 0); // h
  600. }
  601. for (int i = Node_count_half; i < Node_count; i++) // init only neuron nodes
  602. {
  603. double random = RandomD(0, 1);
  604.  
  605. SetV(i, random < percent_stable_state ? V0 + RandomD(-V0 * eps_persent, V0 * eps_persent) : V1 + RandomD(-V1 * eps_persent, V1 * eps_persent)); // V
  606. Setm(i, random < percent_stable_state ? m0 + RandomD(-m0 * eps_persent, m0 * eps_persent) : m1 + RandomD(-m1 * eps_persent, m1 * eps_persent)); // m
  607. Setn(i, random < percent_stable_state ? n0 + RandomD(-n0 * eps_persent, n0 * eps_persent) : n1 + RandomD(-n1 * eps_persent, n1 * eps_persent)); // n
  608. Seth(i, random < percent_stable_state ? h0 + RandomD(-h0 * eps_persent, h0 * eps_persent) : h1 + RandomD(-h1 * eps_persent, h1 * eps_persent)); // h
  609. }
  610.  
  611. for (int i = 0; i < Node_count_half; i++)
  612. {
  613. I_app[i] = 0; // init for 1st half array
  614. }
  615.  
  616. for (int i = Node_count_half; i < Node_count; i++)
  617. {
  618. I_app[i] = 0.93; // init for neurons; Bifurcation point: I_app = 0.82; I_app_up = 1.04
  619. }
  620.  
  621. double percent_excitable = 0.8; // 0.8
  622. double E_syn0 = 0;
  623. double E_syn1 = -90;
  624.  
  625. for (int i = 0; i < Node_count; i++)
  626. {
  627. E_syn[i] = 0; // init for neurons all array
  628. }
  629.  
  630. for (int i = Node_count_half; i < Node_count; i++)
  631. {
  632. double random = RandomD(0, 1);
  633. E_syn[i] = random < percent_excitable ? E_syn0 : E_syn1;
  634. //printf("i = %d\t E_syn = %f\n", i, E_syn[i]);
  635. }
  636.  
  637. for (int i = 0; i < Node_count; i++)
  638. {
  639. GenerateRandomMeander(i, 0);
  640. last_meander_end[i] = Meander_start_from_zero[i] + Duration;
  641. }
  642.  
  643. const double t_start = 0;
  644. const double t_max = 5.0; // 100 msec = 0.1 sec
  645. const double dt = 0.000005; // 0.01 msec = 0.00001 sec; 1 msec = 0.001 sec
  646.  
  647. double t = t_start;
  648.  
  649. //fp0 = fopen("results.txt", "w+");
  650. //setlocale(LC_NUMERIC, "French_Canada.1252");
  651.  
  652. clock_t start_rk4, end_rk4;
  653. start_rk4 = clock();
  654. int lastPercent = -1;
  655.  
  656. FillVOldFromCurrent();
  657.  
  658. //fp_Ca = fopen("results_Ca.txt", "w+");
  659. //fp_IP3 = fopen("results_IP3.txt", "w+");
  660. //fp_z = fopen("results_z.txt", "w+");
  661. //fp_G = fopen("results_G.txt", "w+");
  662. fp_I_stim = fopen("results_I_stim.txt", "w+");
  663. fp_V = fopen("results_V.txt", "w+");
  664. fp_m = fopen("results_m.txt", "w+");
  665. fp_n = fopen("results_n.txt", "w+");
  666. fp_h = fopen("results_h.txt", "w+");
  667. fp_V_spikes = fopen("results_V_spikes.txt", "w+");
  668.  
  669. while (t < t_max || Approximately(t, t_max))
  670. {
  671. //fprintf(fp_Ca, "%f\t", t);
  672. //fprintf(fp_IP3, "%f\t", t);
  673. //fprintf(fp_z, "%f\t", t);
  674. //fprintf(fp_G, "%f\t", t);
  675. fprintf(fp_I_stim, "%f\t", t);
  676. fprintf(fp_V, "%f\t", t);
  677. fprintf(fp_m, "%f\t", t);
  678. fprintf(fp_n, "%f\t", t);
  679. fprintf(fp_h, "%f\t", t);
  680. fprintf(fp_V_spikes, "%f\t", t);
  681.  
  682. for (int i = 0; i < Node_count; i++)
  683. {
  684. if (t > last_meander_end[i])
  685. {
  686. GenerateRandomMeander(i, t);
  687. last_meander_end[i] = Meander_start_from_zero[i] + Duration;
  688. }
  689.  
  690. fprintf(fp_I_stim, "%f\t", I_stim(i, t));
  691. }
  692. fprintf(fp_I_stim, "\n");
  693.  
  694. // ������ �� �������
  695. //for (int i = 0; i < Equations_count; i += 8)
  696. // fprintf(fp_Ca, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // Ca
  697.  
  698. //for (int i = 1; i < Equations_count; i += 8)
  699. // fprintf(fp_IP3, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // IP3
  700.  
  701. //for (int i = 2; i < Equations_count; i += 8)
  702. // fprintf(fp_z, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // z
  703.  
  704. //for (int i = 3; i < Equations_count; i += 8)
  705. // fprintf(fp_G, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // G
  706.  
  707. for (int i = 0; i < Equations_count; i += 4)
  708. fprintf(fp_V, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // V
  709.  
  710. for (int i = 1; i < Equations_count; i += 4)
  711. fprintf(fp_m, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // m
  712.  
  713. for (int i = 2; i < Equations_count; i += 4)
  714. fprintf(fp_n, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // n
  715.  
  716. for (int i = 3; i < Equations_count; i += 4)
  717. fprintf(fp_h, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // h
  718.  
  719. //fprintf(fp_Ca, "\n");
  720. //fprintf(fp_IP3, "\n");
  721. //fprintf(fp_z, "\n");
  722. //fprintf(fp_G, "\n");
  723. fprintf(fp_V, "\n");
  724. fprintf(fp_m, "\n");
  725. fprintf(fp_n, "\n");
  726. fprintf(fp_h, "\n");
  727.  
  728. double f_next[Equations_count];
  729.  
  730. RungeKutta(t, dt, f, f_next);
  731.  
  732. for (int i = 0; i < Equations_count; i += 4)
  733. {
  734. double diff = f_next[i] - f[i];
  735.  
  736. fprintf(fp_V_spikes, i == Equations_count - 1 ? "%d" : "%d\t", diff < 0 && f_diff[i] > 0 && f[i] > -10 ? 1 : 0);
  737.  
  738. f_diff[i] = diff;
  739. }
  740.  
  741. fprintf(fp_V_spikes, "\n");
  742.  
  743. CopyArray(f_next, f, Equations_count);
  744.  
  745. t += dt;
  746.  
  747. int percent = (int)(100 * (t - t_start) / (t_max - t_start));
  748. if (percent != lastPercent)
  749. {
  750. printf("Progress: %d%%\n", percent);
  751. lastPercent = percent;
  752. }
  753.  
  754. //printf("V(24) = %f\t V_old(24) = %f\n", f[24*4], V_old(24));
  755. UpdateVOld();
  756. }
  757.  
  758. //fclose(fp_Ca);
  759. //fclose(fp_IP3);
  760. //fclose(fp_z);
  761. //fclose(fp_G);
  762. fclose(fp_I_stim);
  763. fclose(fp_V);
  764. fclose(fp_m);
  765. fclose(fp_n);
  766. fclose(fp_h);
  767. fclose(fp_V_spikes);
  768.  
  769. end_rk4 = clock();
  770. double extime_rk4 = (double)(end_rk4 - start_rk4) / CLOCKS_PER_SEC;
  771. int minutes = (int)extime_rk4 / 60;
  772. int seconds = (int)extime_rk4 % 60;
  773. printf("\nExecution time is: %d minutes %d seconds\n ", minutes, seconds);
  774.  
  775. fp0 = fopen("time_exec.txt", "w+");
  776. fprintf(fp0, "%f\n", extime_rk4);
  777. fclose(fp0);
  778.  
  779. for (int i = 0; i < Node_count; i++)
  780. free(A[i]);
  781. free(A);
  782.  
  783. for (int i = 0; i < Node_count; i++)
  784. free(B[i]);
  785. free(B);
  786.  
  787. for (int i = 0; i < Node_count; i++)
  788. free(sigma[i]);
  789. free(sigma);
  790.  
  791. free(C);
  792. }
  793.  
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!