API tools faq
paste
SpaceQuester

Untitled

SpaceQuester
Sep 25th, 2020
302
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 49.53 KB | None | 0 0
raw download clone embed print report
  1. #define _CRT_SECURE_NO_WARNINGS
  2.  
  3. #define _USE_MATH_DEFINES
  4. #include "math.h"
  5. #include <stdlib.h>
  6. #include <stdio.h>
  7. #include <locale.h>
  8. #include <time.h>
  9. #include <stdbool.h>
  10. #include <list>
  11. #include <omp.h>
  12. #include <vector>
  13.  
  14. using namespace std;
  15.  
  16. #define Node_count 200
  17.  
  18. #define Equations_per_node 12 // !!! 12 - Don't change !!!
  19. #define Equations_count Node_count * Equations_per_node
  20.  
  21. int MaxDeep = 50;
  22.  
  23. double* f;
  24. double* f_diff;
  25.  
  26. double** k;
  27. double* phi_k1;
  28. double* phi_k2;
  29. double* phi_k3;
  30.  
  31. bool enable_I_syn_out = false;
  32.  
  33. double c_0 = 2; // uM
  34. double c_1 = 0.185;
  35. double v_1 = 6; // s^-1
  36. double v_2 = 0.11; // s^-1
  37. double v_3 = 2.2; // uM/s
  38. double* v_4; // uM/s - Controling parameter // 0.5 //double v_4[Node_count]; // uM/s - Controling parameter //0.495
  39. double v_5 = 0.025; // uM/s
  40. double v_6 = 0.2;  // uM/s
  41. double k_1 = 0.5; // s^-1
  42. double k_2 = 1; // uM
  43. double k_3 = 0.1;
  44. double k_4 = 1.1; // uM/s
  45. double a_2 = 0.14; // uM/s
  46. double d_1 = 0.13; // uM
  47. double d_2 = 1.049; // uM
  48. double d_3 = 0.9434; // uM
  49. double d_5 = 0.082; // uM
  50. double alpha = 0.8;
  51. double tau_IP3 = 7.143; // s
  52. double IP3_star = 0.16; // uM
  53. double d_Ca = 0.001; // 0.001
  54. double d_IP3 = 0.2; // 0.12
  55. double alpha_Glu = 2; // 2
  56. double g_astro; // 3
  57.  
  58. // https://neuronaldynamics.epfl.ch/online/Ch2.S2.html
  59. double C_m = 1; // muF/cm^2
  60. double g_K = 35; // mS/cm^2
  61. double g_Na = 40; // mS/cm^2
  62. double g_L = 0.3; // mS/cm^2
  63. double E_K = -77; // mV
  64. double E_Na = 55; // mV
  65. double E_L = -65; // mV
  66.  
  67. double C_m_P = 1; // muF/cm^2
  68. double g_K_P = 35; // mS/cm^2
  69. double g_Na_P = 40; // mS/cm^2
  70. double g_L_P = 0.3; // mS/cm^2
  71. double E_K_P = -77; // mV
  72. double E_Na_P = 55; // mV
  73. double E_L_P = -65; // mV
  74.  
  75. double p_rewir;
  76. double p_inhib;
  77.  
  78. double I_app_min;
  79. double I_app_max;
  80.  
  81. double g_syn;
  82. double k_syn = 0.2;
  83. double* E_syn;
  84.  
  85. double g_syn_P;
  86. double k_syn_P = 0.2;
  87. double* E_syn_P;
  88.  
  89. double alpha_G_P = 25; //s^-1
  90. double beta_G_P = 500; //s^-1
  91.  
  92. double* I_app;
  93. double* I_app_P;
  94.  
  95. double** A_A;
  96. double** B_A;
  97. double* C_A;
  98.  
  99. double** A_N;
  100. double** B_N;
  101. double* C_N;
  102.  
  103. double** A_N_P;
  104. double** B_N_P;
  105. double* C_N_P;
  106.  
  107. list<double>* V_spikes;
  108. list<double>* V_spikes_Freq;
  109.  
  110. FILE* fp_I_syn;
  111.  
  112. double** tau;
  113.  
  114. #define tau_min  2 // ms
  115. #define tau_max 12 // ms
  116.  
  117. #define ms_to_step 40 // (0.001 / dt) !!! Don't forget !!!
  118.  
  119. #define Max_delay tau_max * ms_to_step
  120. double** V_old_array;
  121.  
  122. double Poisson_Freq; // Hz
  123. //const double Min_magintude = -0.20; // muA/cm^2
  124. double Max_magnitude; // muA/cm^2 // 0.20
  125. const double Duration = 0.002; // sec
  126.  
  127. double* Meander_start_from_zero;
  128. double* Meander_width;
  129. double* Meander_height;
  130. double* Meander_interval;
  131. double* last_meander_end;
  132.  
  133. struct Interval
  134. {
  135.     double T_start;
  136.     double T_end;
  137. };
  138.  
  139. enum Ca_state
  140. {
  141.     Wait_for_up,
  142.     Wait_for_down
  143. };
  144.  
  145.  
  146.  
  147. Ca_state ca_all_state[Node_count];
  148.  
  149. vector<Interval> Ca_all_intervals[Node_count];
  150. vector<Interval> Ca_intervals;
  151.  
  152. //bool thread_count_printed = false;
  153.  
  154. double I_stim(int i, double t)
  155. {
  156.     if (t < Meander_start_from_zero[i])
  157.         return 0;
  158.  
  159.     t -= Meander_start_from_zero[i];
  160.     t = fmod(t, Meander_width[i] + Meander_interval[i]);
  161.  
  162.     return t < Meander_width[i] ? Meander_height[i] : 0;
  163. }
  164.  
  165. double Ca(int i)
  166. {
  167.     return f[i * Equations_per_node];
  168. }
  169.  
  170. void SetCa(int i, double value)
  171. {
  172.     f[i * Equations_per_node] = value;
  173. }
  174.  
  175. double IP3(int i)
  176. {
  177.     return f[i * Equations_per_node + 1];
  178. }
  179.  
  180. void SetIP3(int i, double value)
  181. {
  182.     f[i * Equations_per_node + 1] = value;
  183. }
  184.  
  185. double z(int i)
  186. {
  187.     return f[i * Equations_per_node + 2];
  188. }
  189.  
  190. void Setz(int i, double value)
  191. {
  192.     f[i * Equations_per_node + 2] = value;
  193. }
  194.  
  195. double G_P(int i)
  196. {
  197.     return f[i * Equations_per_node + 3];
  198. }
  199.  
  200. void SetG_P(int i, double value)
  201. {
  202.     f[i * Equations_per_node + 3] = value;
  203. }
  204.  
  205. double V(int i)
  206. {
  207.     return f[i * Equations_per_node + 4];
  208. }
  209.  
  210. void SetV(int i, double value)
  211. {
  212.     f[i * Equations_per_node + 4] = value;
  213. }
  214.  
  215. double m(int i)
  216. {
  217.     return f[i * Equations_per_node + 5];
  218. }
  219.  
  220. void Setm(int i, double value)
  221. {
  222.     f[i * Equations_per_node + 5] = value;
  223. }
  224.  
  225. double n(int i)
  226. {
  227.     return f[i * Equations_per_node + 6];
  228. }
  229.  
  230. void Setn(int i, double value)
  231. {
  232.     f[i * Equations_per_node + 6] = value;
  233. }
  234.  
  235. double h(int i)
  236. {
  237.     return f[i * Equations_per_node + 7];
  238. }
  239.  
  240. void Seth(int i, double value)
  241. {
  242.     f[i * Equations_per_node + 7] = value;
  243. }
  244.  
  245. ///
  246.  
  247. double V_P(int i)
  248. {
  249.     return f[i * Equations_per_node + 8];
  250. }
  251.  
  252. void SetV_P(int i, double value)
  253. {
  254.     f[i * Equations_per_node + 8] = value;
  255. }
  256.  
  257. double m_P(int i)
  258. {
  259.     return f[i * Equations_per_node + 9];
  260. }
  261.  
  262. void Setm_P(int i, double value)
  263. {
  264.     f[i * Equations_per_node + 9] = value;
  265. }
  266.  
  267. double n_P(int i)
  268. {
  269.     return f[i * Equations_per_node + 10];
  270. }
  271.  
  272. void Setn_P(int i, double value)
  273. {
  274.     f[i * Equations_per_node + 10] = value;
  275. }
  276.  
  277. double h_P(int i)
  278. {
  279.     return f[i * Equations_per_node + 11];
  280. }
  281.  
  282. void Seth_P(int i, double value)
  283. {
  284.     f[i * Equations_per_node + 11] = value;
  285. }
  286.  
  287. double V_old(int i, int delay)
  288. {
  289.     return V_old_array[i][Max_delay - 1 - delay];
  290. }
  291.  
  292. int RandomI(int min, int max)
  293. {
  294.     return ((double)rand() / (RAND_MAX - 1)) * (max - min) + min;
  295. }
  296.  
  297. double RandomD(double min, double max)
  298. {
  299.     return ((double)rand() / RAND_MAX) * (max - min) + min;
  300. }
  301.  
  302. double J_channel(double* f, int i)
  303. {
  304.     return c_1 * v_1 * pow(IP3(i), 3) * pow(Ca(i), 3) * pow(z(i), 3) * (c_0 / c_1 - (1 + 1 / c_1) * Ca(i)) / pow((IP3(i) + d_1) * (Ca(i) + d_5), 3);
  305. }
  306.  
  307. double J_PLC(double* f, int i)
  308. {
  309.     return v_4[i] * (Ca(i) + (1 - alpha) * k_4) / (Ca(i) + k_4);
  310. }
  311.  
  312. double J_leak(double* f, int i)
  313. {
  314.     return c_1 * v_2 * (c_0 / c_1 - (1 + 1 / c_1) * Ca(i));
  315. }
  316.  
  317. double J_pump(double* f, int i)
  318. {
  319.     return v_3 * pow(Ca(i), 2) / (pow(k_3, 2) + pow(Ca(i), 2));
  320. }
  321.  
  322. double J_in(double* f, int i)
  323. {
  324.     return v_5 + v_6 * pow(IP3(i), 2) / (pow(k_2, 2) + pow(IP3(i), 2));
  325. }
  326.  
  327. double J_out(double* f, int i)
  328. {
  329.     return k_1 * Ca(i);
  330. }
  331.  
  332. double J_Glu(double* f, int i)
  333. {
  334.     /*double J = 0;
  335.     if (E_syn_P[i] == 0)
  336.     {
  337.         J += alpha_Glu / (1 + exp(-(G_P(i) - 0.25) / 0.01));
  338.     }
  339.     return J;*/
  340.  
  341.     return alpha_Glu / (1 + exp(-(G_P(i) - 0.25) / 0.01));
  342. }
  343.  
  344. double alpha_m(double* f, int i)
  345. {
  346.     return 0.182 * (V(i) + 35) / (1 - exp(-(V(i) + 35) / 9));
  347. }
  348.  
  349. double beta_m(double* f, int i)
  350. {
  351.     return -0.124 * (V(i) + 35) / (1 - exp((V(i) + 35) / 9));
  352. }
  353.  
  354. double alpha_n(double* f, int i)
  355. {
  356.     return 0.02 * (V(i) - 25) / (1 - exp(-(V(i) - 25) / 9));
  357. }
  358.  
  359. double beta_n(double* f, int i)
  360. {
  361.     return -0.002 * (V(i) - 25) / (1 - exp((V(i) - 25) / 9));
  362. }
  363.  
  364. double alpha_h(double* f, int i)
  365. {
  366.     return 0.25 * exp(-(V(i) + 90) / 12);
  367. }
  368.  
  369. double beta_h(double* f, int i)
  370. {
  371.     return 0.25 * exp((V(i) + 62) / 6) / exp((V(i) + 90) / 12);
  372. }
  373.  
  374. //
  375.  
  376. double alpha_m_P(double* f, int i)
  377. {
  378.     return 0.182 * (V_P(i) + 35) / (1 - exp(-(V_P(i) + 35) / 9));
  379. }
  380.  
  381. double beta_m_P(double* f, int i)
  382. {
  383.     return -0.124 * (V_P(i) + 35) / (1 - exp((V_P(i) + 35) / 9));
  384. }
  385.  
  386. double alpha_n_P(double* f, int i)
  387. {
  388.     return 0.02 * (V_P(i) - 25) / (1 - exp(-(V_P(i) - 25) / 9));
  389. }
  390.  
  391. double beta_n_P(double* f, int i)
  392. {
  393.     return -0.002 * (V_P(i) - 25) / (1 - exp((V_P(i) - 25) / 9));
  394. }
  395.  
  396. double alpha_h_P(double* f, int i)
  397. {
  398.     return 0.25 * exp(-(V_P(i) + 90) / 12);
  399. }
  400.  
  401. double beta_h_P(double* f, int i)
  402. {
  403.     return 0.25 * exp((V_P(i) + 62) / 6) / exp((V_P(i) + 90) / 12);
  404. }
  405.  
  406. double UllahJung_HodgkinHuxley(int i, double* f, double t)
  407. {
  408.     int in = i / Equations_per_node;
  409.     int il = i % Equations_per_node;
  410.  
  411.     switch (il)
  412.     {
  413.     case 0: // Ca
  414.     {
  415.         double sum_1 = 0;
  416.  
  417.         /*for (int j = 0; j < Node_count; j++)
  418.         {
  419.           sum_1 += d_Ca * (Ca(j) - Ca(in));
  420.         }*/
  421.  
  422.         for (int j = 0; j < C_A[in]; j++)
  423.         {
  424.             sum_1 += d_Ca * (Ca((int)B_A[in][j]) - Ca(in));
  425.         }
  426.  
  427.         return J_channel(f, in) - J_pump(f, in) + J_leak(f, in) + J_in(f, in) - J_out(f, in) + sum_1;
  428.     }
  429.  
  430.     case 1: // IP3
  431.     {
  432.         double sum_2 = 0;
  433.  
  434.         /*for (int j = 0; j < Node_count; j++)
  435.         {
  436.         sum_2 += d_IP3 * (IP3(j) - IP3(in));
  437.         }*/
  438.  
  439.         for (int j = 0; j < C_A[in]; j++)
  440.         {
  441.             sum_2 += d_IP3 * (IP3((int)B_A[in][j]) - IP3(in));
  442.         }
  443.  
  444.         return (IP3_star - IP3(in)) / tau_IP3 + J_PLC(f, in) + sum_2 + J_Glu(f, in);
  445.     }
  446.  
  447.     case 2: // z
  448.     {
  449.         return a_2 * (d_2 * (IP3(in) + d_1) / (IP3(in) + d_3) * (1 - z(in)) - Ca(in) * z(in));
  450.     }
  451.  
  452.     case 3: // G_P
  453.     {
  454.         return -alpha_G_P * G_P(in) + beta_G_P * (1 / (1 + exp(-V_P(in) / 0.5)));
  455.     }
  456.  
  457.     case 4: // V
  458.     {
  459.         double I_syn = 0;
  460.         double I_syn_P = 0;
  461.  
  462.         /*for (int j = 0; j < Node_count; j++)
  463.         {
  464.         //sum += A[in][j] * g_syn * (V(in) - V_old(j, tau[in][j]));
  465.         //sum += A[in][j] * g_syn * (V(j) - V(in));
  466.         //sum += A[in][j] * g_syn * (V(in) - E_syn[in]) / (1 + exp(-V_old(j, tau[in][j]) / k_syn));
  467.         //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
  468.         //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
  469.           I_syn += A_N[in][j] * g_syn * (E_syn[in] - V(in)) / (1 + exp(-(V(j) / k_syn)));
  470.         //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);
  471.       }*/
  472.  
  473.       /*for (int j = 0; j < C[in]; j++)
  474.       {
  475.       sum += sigma[in][(int)B[in][j]] * (V((int)B[in][j]) - V(in));
  476.       }*/
  477.  
  478.         for (int j = 0; j < C_N[in]; j++)
  479.         {
  480.             if ((1 + g_astro * Ca(in)) > 0)
  481.             {
  482.                 if (Ca(in) >= 0.3)
  483.                 {
  484.                     I_syn += g_syn * (1 + g_astro * Ca(in)) * (V(in) - E_syn[(int)B_N[in][j]]) / (1 + exp(-(V((int)B_N[in][j]) / k_syn))); // версия с V (без V_old)
  485.                 }
  486.                 else
  487.                 {
  488.                     I_syn += g_syn * (V(in) - E_syn[(int)B_N[in][j]]) / (1 + exp(-(V((int)B_N[in][j]) / k_syn))); // версия с V (без V_old)
  489.                 }
  490.             }
  491.             //I_syn += g_syn * (1 + g_astro * Ca(in)) * (E_syn[in] - V(in)) / (1 + exp(-(V_old((int)B_N[in][j], tau[in][(int)B_N[in][j]]) / k_syn)); // версия с V_old
  492.             //I_syn += g_syn * (1 + g_astro * Ca(in)) * (E_syn[in] - V(in)) / (1 + exp(-(V(j) / k_syn))); // версия без с V (без V_old)
  493.             //I_syn += g_syn * (E_syn[in] - V(in)) / (1 + exp(-(V(j) / k_syn))); // версия без с V (без V_old), упрощенная версия
  494.             // sum_3 += g_syn * (1 + g_astro * Ca(in)) * (E_syn[i] - V(i)) / (1 + exp(-(V(j) / k_syn))); // образец из старой версии
  495.             //I_syn += g_syn * (E_syn[in] - V(in)) / (1 + exp(-(V_old((int)B_N[in][j], tau[in][(int)B_N[in][j]])) / k_syn));*/
  496.           //I_syn += g_syn * (E_syn[(int)B_N[in][j]] - V(in)) / (1 + exp(-(V(j) / k_syn))); // версия с V (без V_old)
  497.           //I_syn += g_syn * (E_syn[(int)B_N[in][j]] - V(in)) / (1 + exp(-(V((int)B_N[in][j]) / k_syn))); // версия с V (без V_old), упрощенная версия !!!
  498.           //I_syn += g_syn * (E_syn[(int)B_N[in][j]] - V(in)) / (1 + exp(-(V_old((int)B_N[in][j], tau[in][(int)B_N[in][j]]) / k_syn))); // версия с V (без V_old), упрощенная версия !!!
  499.         //sum += g_syn * (V((int)B[in][j]) - V(in)); // устаревшая часть, нужна для проверки разностной схемы
  500.         //sum += A[in][j] * g_syn * (V(in) - E_syn[in]) / (1 + exp(-V_old((int)B[in][j]) / k_syn)); // устаревшая часть
  501.         //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));
  502.         //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));
  503.         //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
  504.         //sum += 1 / (0.2 * Node_count) * /*(int)A_N[in][(int)B_N[in][j]] * */ g_syn * (1 + g_astro * Ca(in)) * (V(in) - E_syn[in]) / (1 + exp(-V((int)B_N[in][j]) / k_syn)); // i up, j down
  505.           // i up, j down
  506.           /*printf("i = %d\t j = %d\t A[i, j] = %d\n", in, (int)B_N[in][j], (int)A_N[in][(int)B_N[in][j]]);*/
  507.         /*printf("i = %d\t V_old = %f\t exp = %f\n", in, Vold, ee);*/
  508.         }
  509.         /*printf("i = %d\t sum = %f\n", in, sum);*/
  510.  
  511.         /*if ((1 + g_astro * Ca(in)) > 0)
  512.         {
  513.             if (Ca(in) >= 0.3)
  514.             {
  515.                 I_syn_P += g_syn_P * (1 + g_astro * Ca(in)) * (V_P(in) - E_syn_P[in]) / (1 + exp(-(V_P(in) / k_syn_P))); // версия с V (без V_old)
  516.             }
  517.             else
  518.             {
  519.                 I_syn_P += g_syn_P * (V_P(in) - E_syn_P[in]) / (1 + exp(-(V_P(in) / k_syn_P))); // версия с V (без V_old)
  520.             }
  521.         }*/
  522.         I_syn_P += g_syn_P * (V_P(in) - E_syn_P[in]) / (1 + exp(-(V_P(in) / k_syn_P))); // версия с V (без V_old)
  523.  
  524.         if (enable_I_syn_out)
  525.             fprintf(fp_I_syn, i == Equations_count - 1 ? "%f" : "%f\t", I_syn);
  526.  
  527.         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_syn + I_syn_P) / C_m); // V
  528.     }
  529.  
  530.     case 5: // m
  531.     {
  532.         return 1000 * (alpha_m(f, in) * (1 - m(in)) - beta_m(f, in) * m(in)); // m
  533.     }
  534.  
  535.     case 6: // n
  536.     {
  537.         return 1000 * (alpha_n(f, in) * (1 - n(in)) - beta_n(f, in) * n(in)); // n
  538.     }
  539.  
  540.     case 7: // h
  541.     {
  542.         return 1000 * (alpha_h(f, in) * (1 - h(in)) - beta_h(f, in) * h(in)); // h
  543.     }
  544.  
  545.     case 8: // V_P
  546.     {
  547.         return 1000 * ((g_Na_P * pow(m_P(in), 3) * h_P(in) * (E_Na_P - V_P(in)) + g_K_P * n_P(in) * (E_K_P - V_P(in)) + g_L_P * (E_L_P - V_P(in)) + I_app_P[in] + I_stim(in, t)) / C_m_P); // V_P
  548.     }
  549.  
  550.     case 9: // m_P
  551.     {
  552.         return 1000 * (alpha_m_P(f, in) * (1 - m_P(in)) - beta_m_P(f, in) * m_P(in)); // m_P
  553.     }
  554.  
  555.     case 10: // n_P
  556.     {
  557.         return 1000 * (alpha_n_P(f, in) * (1 - n_P(in)) - beta_n_P(f, in) * n_P(in)); // n_P
  558.     }
  559.  
  560.     case 11: // h_P
  561.     {
  562.         return 1000 * (alpha_h_P(f, in) * (1 - h_P(in)) - beta_h_P(f, in) * h_P(in)); // h_P
  563.     }
  564.     }
  565.  
  566.     return 0;
  567. }
  568.  
  569. void RungeKutta(double t, double dt, double* f, double* f_next)
  570. {
  571.     // k1
  572. #pragma omp parallel for
  573.     for (int i = 0; i < Equations_count; i++)
  574.     {
  575.         //if (!thread_count_printed)
  576.         //{
  577.         //  thread_count_printed = true;
  578.         //  printf("Threads = %d\n", omp_get_num_threads());
  579.         //}
  580.  
  581.         k[i][0] = UllahJung_HodgkinHuxley(i, f, t) * dt;
  582.         phi_k1[i] = f[i] + k[i][0] / 2;
  583.         k[i][1] = UllahJung_HodgkinHuxley(i, phi_k1, t) * dt;
  584.         phi_k2[i] = f[i] + k[i][1] / 2;
  585.         k[i][2] = UllahJung_HodgkinHuxley(i, phi_k2, t) * dt;
  586.         phi_k3[i] = f[i] + k[i][2] / 2;
  587.         k[i][3] = UllahJung_HodgkinHuxley(i, phi_k3, t) * dt;
  588.         f_next[i] = f[i] + (k[i][0] + 2 * k[i][1] + 2 * k[i][2] + k[i][3]) / 6;
  589.     }
  590.  
  591.     //for (int i = 0; i < Equations_count; i++)
  592.     //  phi_k1[i] = f[i] + k[i][0] / 2;
  593.  
  594.     // k2
  595.     //for (int i = 0; i < Equations_count; i++)
  596.     //  k[i][1] = UllahJung_HodgkinHuxley(i, phi_k1, t) * dt;
  597.  
  598.  
  599.     //for (int i = 0; i < Equations_count; i++)
  600.     //  phi_k2[i] = f[i] + k[i][1] / 2;
  601.  
  602.     // k3
  603.     //for (int i = 0; i < Equations_count; i++)
  604.     //  k[i][2] = UllahJung_HodgkinHuxley(i, phi_k2, t) * dt;
  605.  
  606.  
  607.     //for (int i = 0; i < Equations_count; i++)
  608.     //  phi_k3[i] = f[i] + k[i][2] / 2;
  609.  
  610.     //enable_I_syn_out = true;
  611.  
  612.     // k4
  613.     //for (int i = 0; i < Equations_count; i++)
  614.     //  k[i][3] = UllahJung_HodgkinHuxley(i, phi_k3, t) * dt;
  615.  
  616.     //enable_I_syn_out = false;
  617.  
  618.     //for (int i = 0; i < Equations_count; i++)
  619.     //  f_next[i] = f[i] + (k[i][0] + 2 * k[i][1] + 2 * k[i][2] + k[i][3]) / 6;
  620. }
  621.  
  622. void CopyArray(double* source, double* target, int N)
  623. {
  624.     for (int i = 0; i < N; i++)
  625.         target[i] = source[i];
  626. }
  627.  
  628. bool Approximately(double a, double b)
  629. {
  630.     if (a < 0)
  631.         a = -a;
  632.  
  633.     if (b < 0)
  634.         b = -b;
  635.  
  636.     return a - b <= 0.000001;
  637. }
  638.  
  639. //bool CheckSameLine(int i, int j)
  640. //{
  641. //  return i / Node_wire_width == j / Node_wire_width;
  642. //}
  643. //
  644. //bool IsWireNeighbors(int i, int j)
  645. //{
  646. //  if (CheckSameLine(i, j) && (i == j - 1 || i == j + 1))
  647. //      return true;
  648. //
  649. //  if (i == j - Node_wire_width || i == j + Node_wire_width)
  650. //      return true;
  651. //
  652. //  return false;
  653. //}
  654.  
  655. // http://preshing.com/20111007/how-to-generate-random-timings-for-a-poisson-process/
  656. double nextTime(double rateParameter)
  657. {
  658.     return -log(1.0 - (double)rand() / (RAND_MAX)) / rateParameter;
  659. }
  660.  
  661. void GenerateRandomMeander(int i, double min_start_time)
  662. {
  663.     double offset = nextTime(Poisson_Freq);
  664.  
  665.     if (offset < 0)
  666.     {
  667.         int a = 0;
  668.     }
  669.  
  670.     Meander_start_from_zero[i] = min_start_time + offset;
  671.     Meander_width[i] = Duration;
  672.     //Meander_height[i] = RandomD(-Max_magnitude, Max_magnitude);
  673.     Meander_height[i] = RandomD(0, Max_magnitude);
  674.     //Meander_height[i] = Max_magnitude;
  675. }
  676.  
  677. void FillAMatrixZero()
  678. {
  679.     for (int i = 0; i < Node_count; i++)
  680.     {
  681.         for (int j = 0; j < Node_count; j++)
  682.         {
  683.             A_A[i][j] = 0;
  684.             A_N[i][j] = 0;
  685.             A_N_P[i][j] = 0;
  686.         }
  687.     }
  688. }
  689.  
  690. void FillAstrociteMatrix()
  691. {
  692.     for (int i = 0; i < Node_count; i++)
  693.     {
  694.         for (int j = 0; j < Node_count; j++)
  695.         {
  696.             if (i == j)
  697.             {
  698.                 A_A[i][j] = 0;
  699.                 continue;
  700.             }
  701.  
  702.             if (i > j)
  703.             {
  704.                 A_A[i][j] = A_A[j][i];
  705.                 continue;
  706.             }
  707.  
  708.             if (i == 0 && j == Node_count - 1)
  709.             {
  710.                 A_A[i][j] = 1;
  711.                 continue;
  712.             }
  713.  
  714.             if (i == Node_count - 1 && j == 0)
  715.             {
  716.                 A_A[i][j] = 1;
  717.                 continue;
  718.             }
  719.  
  720.             if (i == j - 1 || i == j + 1)
  721.             {
  722.                 A_A[i][j] = 1;
  723.                 continue;
  724.             }
  725.         }
  726.     }
  727.     //A_A[0][1] = 0; // only for debug. diffusion Ca test
  728.     //A_A[1][0] = 0; // only for debug. diffusion Ca test
  729.     //A_A[1][3] = 0;
  730.     //A_A[3][1] = 0;
  731.     //A_A[0][2] = 0;
  732.     //A_A[2][0] = 0;
  733. }
  734.  
  735. void FillBCMatrix_A()
  736. {
  737.     for (int i = 0; i < Node_count; i++)
  738.     {
  739.         int bIndex = 0;
  740.         C_A[i] = 0;
  741.         for (int j = 0; j < Node_count; j++)
  742.         {
  743.             if (A_A[i][j] == 1)
  744.             {
  745.                 B_A[i][bIndex] = j;
  746.                 bIndex++;
  747.                 C_A[i]++;
  748.             }
  749.         }
  750.     }
  751. }
  752.  
  753. void FillBCMatrix_N()
  754. {
  755.     for (int i = 0; i < Node_count; i++)
  756.     {
  757.         int bIndex = 0;
  758.         C_N[i] = 0;
  759.         for (int j = 0; j < Node_count; j++)
  760.         {
  761.             if (A_N[i][j] == 1)
  762.             {
  763.                 B_N[i][bIndex] = j;
  764.                 bIndex++;
  765.                 C_N[i]++;
  766.             }
  767.         }
  768.     }
  769. }
  770.  
  771. void FillBCMatrix_N_P()
  772. {
  773.     for (int i = 0; i < Node_count; i++)
  774.     {
  775.         int bIndex = 0;
  776.         C_N_P[i] = 0;
  777.         for (int j = 0; j < Node_count; j++)
  778.         {
  779.             if (A_N_P[i][j] == 1)
  780.             {
  781.                 B_N_P[i][bIndex] = j;
  782.                 bIndex++;
  783.                 C_N_P[i]++;
  784.             }
  785.         }
  786.     }
  787. }
  788.  
  789. bool IsWireNeighbors(int i, int j, int deep)
  790. {
  791.     int j_border_left = j - deep < 0 ? j + Node_count : j;
  792.     int j_border_right = j + deep >= Node_count ? j - Node_count : j;
  793.  
  794.     if (i == j_border_left - deep || i == j_border_right + deep)
  795.     {
  796.         return true;
  797.     }
  798.  
  799.     return false;
  800. }
  801.  
  802. void FillNeuronMatrix()
  803. {
  804.     for (int i = 0; i < Node_count; i++)
  805.     {
  806.         for (int j = 0; j < Node_count; j++)
  807.         {
  808.             if (i == j)
  809.             {
  810.                 A_N[i][j] = 0;
  811.                 continue;
  812.             }
  813.  
  814.             if (i > j)
  815.             {
  816.                 A_N[i][j] = A_N[j][i];
  817.                 continue;
  818.             }
  819.  
  820.             for (int deep = 1; deep <= MaxDeep; deep++)
  821.             {
  822.                 if (IsWireNeighbors(i, j, deep))
  823.                     A_N[i][j] = 1;
  824.             }
  825.         }
  826.     }
  827. }
  828.  
  829. void RandomizeNeuronMatrix()
  830. {
  831.     srand(time(NULL));
  832.  
  833.     for (int i = 0; i < Node_count; i++)
  834.     {
  835.         //if (i == Node_count / 2)
  836.         //  srand(time(NULL));
  837.  
  838.         for (int link = 0; link < MaxDeep * 2; link++)
  839.         {
  840.             double x = RandomD(0, 1);
  841.  
  842.             if (x > p_rewir)
  843.                 continue;
  844.  
  845.             int rndJ;
  846.  
  847.             do
  848.             {
  849.                 rndJ = RandomI(0, Node_count);
  850.             } while (i == rndJ || A_N[i][rndJ] == 1);
  851.  
  852.             int rndJ_last;
  853.  
  854.             do
  855.             {
  856.                 rndJ_last = RandomI(i - MaxDeep - 1, i + MaxDeep + 1);
  857.  
  858.                 if (rndJ_last < 0)
  859.                     rndJ_last += Node_count;
  860.                 else if (rndJ_last >= Node_count)
  861.                     rndJ_last -= Node_count;
  862.  
  863.             } while (i == rndJ_last || A_N[i][rndJ_last] == 0);
  864.  
  865.             A_N[i][rndJ_last] = 0;
  866.  
  867.             A_N[i][rndJ] = 1;
  868.         }
  869.     }
  870. }
  871.  
  872. void FillNeuronPoissonMatrix()
  873. {
  874.     for (int i = 0; i < Node_count; i++)
  875.     {
  876.         for (int j = 0; j < Node_count; j++)
  877.         {
  878.             A_N_P[i][j] = 0;
  879.         }
  880.     }
  881. }
  882.  
  883. void FillVOldFromCurrent()
  884. {
  885.     for (int i = 0; i < Node_count; i++)
  886.         for (int j = 0; j < Max_delay; j++)
  887.             V_old_array[i][j] = V(i);
  888. }
  889.  
  890. void UpdateVOld()
  891. {
  892.     for (int i = 0; i < Node_count; i++)
  893.     {
  894.         for (int j = 1; j < Max_delay; j++)
  895.             V_old_array[i][j - 1] = V_old_array[i][j];
  896.  
  897.         V_old_array[i][Max_delay - 1] = V(i);
  898.     }
  899. }
  900.  
  901. //void FillFullTauMatrix()
  902. //{
  903. //  for (int i = 0; i < Node_count; i++)
  904. //  {
  905. //      for (int j = 0; j < Node_count; j++)
  906. //      {
  907. //          if (i < Node_count || j < Node_count)
  908. //          {
  909. //              tau[i][j] = 0;
  910. //              continue;
  911. //          }
  912. //
  913. //          int i_neuron = i - Node_count;
  914. //          int j_neuron = j - Node_count;
  915. //
  916. //          int i_wire_x = i_neuron / Node_wire_width;
  917. //          int i_wire_y = i_neuron % Node_wire_width;
  918. //
  919. //          int j_wire_x = j_neuron / Node_wire_width;
  920. //          int j_wire_y = j_neuron % Node_wire_width;
  921. //
  922. //          double distance_max = sqrt(2.) * (Node_wire_width - 1);
  923. //          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));
  924. //
  925. //          tau[i][j] = (tau_min + distance / (distance_max) * (tau_max - tau_min)) * ms_to_step;
  926. //      }
  927. //  }
  928. //}
  929. //
  930. //void FillTauMatrix()
  931. //{
  932. //  for (int i = 0; i < Node_count; i++)
  933. //  {
  934. //      for (int j = 0; j < Node_count; j++)
  935. //      {
  936. //          if (i == j || A_N[i][j] == 0)
  937. //          {
  938. //              tau[i][j] = 0;
  939. //              continue;
  940. //          }
  941. //
  942. //          int i_neuron = i;
  943. //          int j_neuron = j;
  944. //
  945. //          int i_wire_x = i_neuron / Node_wire_width;
  946. //          int i_wire_y = i_neuron % Node_wire_width;
  947. //
  948. //          int j_wire_x = j_neuron / Node_wire_width;
  949. //          int j_wire_y = j_neuron % Node_wire_width;
  950. //
  951. //          double distance_max = sqrt(2.) * (Node_wire_width - 1);
  952. //          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));
  953. //
  954. //          double t = (distance - 1) / (distance_max - 1);
  955. //          tau[i][j] = (tau_min + t * (tau_max - tau_min)) * ms_to_step;
  956. //      }
  957. //  }
  958. //}
  959.  
  960. int main(int argc, char* argv[])
  961. {
  962.     // run like: UJ_HH_Ring_acc.out 250 0.2 0.4 0.05 3.0 1.05 1.50 // (1) p_rewir (2) g_syn (3) g_astro (4) Max_magnitude (5) g_syn_P (6) Poisson_Freq
  963.     /*sscanf(argv[1], "%d", &Node_count);
  964.     FILE* fp_Node_count;
  965.     fp_Node_count = fopen("Node_count.txt", "w");
  966.     fprintf(fp_Node_count, "%d\t", Node_count);
  967.     fclose(fp_Node_count);
  968.     printf("Node_count = %d\n", Node_count);*/
  969.  
  970.     sscanf(argv[1], "%lf", &p_rewir);
  971.     FILE* fp_p_rewir;
  972.     fp_p_rewir = fopen("p_rewir.txt", "w");
  973.     fprintf(fp_p_rewir, "%f\t", p_rewir);
  974.     fclose(fp_p_rewir);
  975.     printf("p_rewir = %f\n", p_rewir);
  976.  
  977.     /*sscanf(argv[3], "%lf", &p_inhib);
  978.     FILE* fp_p_inhib;
  979.     fp_p_inhib = fopen("p_inhib.txt", "w");
  980.     fprintf(fp_p_inhib, "%f\t", p_inhib);
  981.     fclose(fp_p_inhib);
  982.     printf("p_inhib = %f\n", p_inhib);*/
  983.  
  984.     sscanf(argv[2], "%lf", &g_syn);
  985.     FILE* fp_g_syn;
  986.     fp_g_syn = fopen("g_syn.txt", "w");
  987.     fprintf(fp_g_syn, "%f\t", g_syn);
  988.     fclose(fp_g_syn);
  989.     printf("g_syn = %f\n", g_syn);
  990.  
  991.     sscanf(argv[3], "%lf", &g_astro);
  992.     FILE* fp_g_astro;
  993.     fp_g_astro = fopen("g_astro.txt", "w");
  994.     fprintf(fp_g_astro, "%f\t", g_astro);
  995.     fclose(fp_g_astro);
  996.     printf("g_astro = %f\n", g_astro);
  997.  
  998.     sscanf(argv[4], "%lf", &Max_magnitude);
  999.     FILE* fp_Max_magnitude;
  1000.     fp_Max_magnitude = fopen("Max_magnitude.txt", "w");
  1001.     fprintf(fp_Max_magnitude, "%f\t", Max_magnitude);
  1002.     fclose(fp_Max_magnitude);
  1003.     printf("I_stim magnitude = %f\n", Max_magnitude);
  1004.  
  1005.     sscanf(argv[5], "%lf", &g_syn_P);
  1006.     FILE* fp_g_syn_P;
  1007.     fp_g_syn_P = fopen("g_syn_P.txt", "w");
  1008.     fprintf(fp_g_syn_P, "%f\t", g_syn_P);
  1009.     fclose(fp_g_syn_P);
  1010.     printf("g_syn_P = %f\n", g_syn_P);
  1011.  
  1012.     sscanf(argv[6], "%lf", &Poisson_Freq);
  1013.     FILE* fp_Poisson_Freq;
  1014.     fp_Poisson_Freq = fopen("Poisson_Freq.txt", "w");
  1015.     fprintf(fp_Poisson_Freq, "%f\t", Poisson_Freq);
  1016.     fclose(fp_Poisson_Freq);
  1017.     printf("Poisson_Freq = %f\n", Poisson_Freq);
  1018.  
  1019.     //sscanf(argv[6], "%lf", &I_app_min);
  1020.     //sscanf(argv[7], "%lf", &I_app_max);
  1021.  
  1022.     for (int i = 0; i < Node_count; i++)
  1023.         ca_all_state[i] = Ca_state::Wait_for_up;
  1024.  
  1025.     f = new double[Equations_count];
  1026.     f_diff = new double[Equations_count];
  1027.     v_4 = new double[Node_count];
  1028.     E_syn = new double[Node_count];
  1029.     I_app = new double[Node_count];
  1030.     E_syn_P = new double[Node_count];
  1031.     I_app_P = new double[Node_count];
  1032.     V_spikes = new list<double>[Node_count];
  1033.     V_spikes_Freq = new list<double>[Node_count];
  1034.  
  1035.     Meander_start_from_zero = new double[Node_count];
  1036.     Meander_width = new double[Node_count];
  1037.     Meander_height = new double[Node_count];
  1038.     Meander_interval = new double[Node_count];
  1039.     last_meander_end = new double[Node_count];
  1040.  
  1041.     tau = new double* [Node_count];
  1042.     for (int i = 0; i < Node_count; i++)
  1043.         tau[i] = new double[Node_count];
  1044.  
  1045.     V_old_array = new double* [Node_count];
  1046.     for (int i = 0; i < Node_count; i++)
  1047.         V_old_array[i] = new double[Max_delay];
  1048.  
  1049.     FILE* fp0;
  1050.     FILE* fp_I_stim;
  1051.     FILE* fp_Ca;
  1052.     FILE* fp_IP3;
  1053.     //FILE *fp_z;
  1054.     //FILE* fp_G_P;
  1055.     FILE* fp_V;
  1056.     FILE* fp_V_P;
  1057.     //FILE *fp_m;
  1058.     //FILE *fp_n;
  1059.     //FILE *fp_h;
  1060.     FILE* fp_V_spikes;
  1061.     FILE* fp_Esyn;
  1062.     FILE* fp_Esyn_P;
  1063.  
  1064.     //FILE* fp_res;
  1065.     srand(time(NULL));
  1066.  
  1067.     //for (int i = 0; i < Node_count; i++)
  1068.     //  V_old_length[i] = 0;
  1069.  
  1070.     A_A = new double* [Node_count];
  1071.     for (int i = 0; i < Node_count; i++)
  1072.         A_A[i] = new double[Node_count];
  1073.  
  1074.     B_A = new double* [Node_count];
  1075.     for (int i = 0; i < Node_count; i++)
  1076.         B_A[i] = new double[Node_count];
  1077.  
  1078.     C_A = new double[Node_count];
  1079.  
  1080.     A_N = new double* [Node_count];
  1081.     for (int i = 0; i < Node_count; i++)
  1082.         A_N[i] = new double[Node_count];
  1083.  
  1084.     B_N = new double* [Node_count];
  1085.     for (int i = 0; i < Node_count; i++)
  1086.         B_N[i] = new double[Node_count];
  1087.  
  1088.     C_N = new double[Node_count];
  1089.  
  1090.     A_N_P = new double* [Node_count];
  1091.     for (int i = 0; i < Node_count; i++)
  1092.         A_N_P[i] = new double[Node_count];
  1093.  
  1094.     B_N_P = new double* [Node_count];
  1095.     for (int i = 0; i < Node_count; i++)
  1096.         B_N_P[i] = new double[Node_count];
  1097.  
  1098.     C_N_P = new double[Node_count];
  1099.  
  1100.     FillAMatrixZero();
  1101.     FillAstrociteMatrix();
  1102.     FillNeuronMatrix();
  1103.     RandomizeNeuronMatrix();
  1104.     FillNeuronPoissonMatrix();
  1105.     FillBCMatrix_A();
  1106.     FillBCMatrix_N();
  1107.     FillBCMatrix_N_P();
  1108.     //FillTauMatrix();
  1109.  
  1110.     fp0 = fopen("A_A.txt", "w+");
  1111.     for (int i = 0; i < Node_count; i++)
  1112.     {
  1113.         for (int j = 0; j < Node_count; j++)
  1114.         {
  1115.             fprintf(fp0, "%d\t", (int)A_A[i][j]);
  1116.         }
  1117.         fprintf(fp0, "\n");
  1118.     }
  1119.     fclose(fp0);
  1120.  
  1121.     fp0 = fopen("A_N.txt", "w+");
  1122.     for (int i = 0; i < Node_count; i++)
  1123.     {
  1124.         for (int j = 0; j < Node_count; j++)
  1125.         {
  1126.             fprintf(fp0, "%d\t", (int)A_N[i][j]);
  1127.         }
  1128.         fprintf(fp0, "\n");
  1129.     }
  1130.     fclose(fp0);
  1131.  
  1132.     fp0 = fopen("tau.txt", "w+");
  1133.     for (int i = 0; i < Node_count; i++)
  1134.     {
  1135.         for (int j = 0; j < Node_count; j++)
  1136.         {
  1137.             fprintf(fp0, "%f\t", tau[i][j] / ms_to_step);
  1138.         }
  1139.         fprintf(fp0, "\n");
  1140.     }
  1141.     fclose(fp0);
  1142.  
  1143.     // Write to file number of links for each neuron
  1144.     fp0 = fopen("links.txt", "w+");
  1145.     for (int i = 0; i < Node_count; i++)
  1146.     {
  1147.         int links_count = 0;
  1148.         for (int j = 0; j < Node_count; j++)
  1149.         {
  1150.             if (A_N[i][j] == 1)
  1151.             {
  1152.                 links_count++;
  1153.             }
  1154.         }
  1155.         fprintf(fp0, "%d\n", (int)links_count);
  1156.     }
  1157.     fclose(fp0);
  1158.  
  1159.     //setlocale(LC_NUMERIC, "French_Canada.1252");
  1160.     fp0 = fopen("test_Poisson.txt", "w+");
  1161.     for (int i = 0; i < 1000; i++)
  1162.         fprintf(fp0, "%f\n", nextTime(Poisson_Freq));
  1163.     fclose(fp0);
  1164.  
  1165.     fp0 = fopen("B_A.txt", "w+");
  1166.     for (int i = 0; i < Node_count; i++)
  1167.     {
  1168.         for (int j = 0; j < C_A[i]; j++)
  1169.         {
  1170.             fprintf(fp0, "%d\t", (int)B_A[i][j]);
  1171.         }
  1172.         fprintf(fp0, "\n");
  1173.     }
  1174.     fclose(fp0);
  1175.  
  1176.     fp0 = fopen("B_N.txt", "w+");
  1177.     for (int i = 0; i < Node_count; i++)
  1178.     {
  1179.         for (int j = 0; j < C_N[i]; j++)
  1180.         {
  1181.             fprintf(fp0, "%d\t", (int)B_N[i][j]);
  1182.         }
  1183.         fprintf(fp0, "\n");
  1184.     }
  1185.     fclose(fp0);
  1186.  
  1187.     fp0 = fopen("C_A.txt", "w+");
  1188.     for (int i = 0; i < Node_count; i++)
  1189.     {
  1190.         fprintf(fp0, "%d\n", (int)C_A[i]);
  1191.     }
  1192.     fclose(fp0);
  1193.  
  1194.     fp0 = fopen("C_N.txt", "w+");
  1195.     for (int i = 0; i < Node_count; i++)
  1196.     {
  1197.         fprintf(fp0, "%d\n", (int)C_N[i]);
  1198.     }
  1199.     fclose(fp0);
  1200.  
  1201.     /*for (int i = 0; i < 6; i++)
  1202.     {
  1203.         v_4[i] = 0.6;
  1204.     }*/
  1205.     for (int i = 0; i < Node_count; i++)
  1206.     {
  1207.         v_4[i] = 0.4; // 0.4
  1208.     }
  1209.  
  1210.     // Initial values
  1211.     /*for (int i = 0; i < Equations_count; i++)
  1212.     {
  1213.       f[i] = 0;
  1214.   }*/
  1215.  
  1216.   //I_app_min = 1.1;
  1217.   //I_app_max = 1.5;
  1218.  
  1219.     FILE* fp_I_app;
  1220.     fp_I_app = fopen("I_app.txt", "w");
  1221.     for (int i = 0; i < Node_count; i++)
  1222.     {
  1223.         I_app[i] = 0.7; //RandomD(I_app_min, I_app_max);
  1224.         fprintf(fp_I_app, "%f\n", I_app[i]);
  1225.     }
  1226.     fclose(fp_I_app);
  1227.  
  1228.     FILE* fp_I_app_P;
  1229.     fp_I_app_P = fopen("I_app_P.txt", "w");
  1230.     for (int i = 0; i < Node_count; i++)
  1231.     {
  1232.         I_app_P[i] = 0.7;
  1233.         fprintf(fp_I_app_P, "%f\n", I_app_P[i]);
  1234.     }
  1235.     fclose(fp_I_app_P);
  1236.  
  1237.     for (int i = 0; i < Node_count; i++) // init array for all nodes
  1238.     {
  1239.         SetG_P(i, 0); // G_P
  1240.     }
  1241.  
  1242.     double percent_stable_state = 0.50; // 0.40
  1243.     double eps_persent = 0.05; //0.05
  1244.  
  1245.     double Ca0 = 0.07;
  1246.     double IP30 = 0.16;
  1247.     double z0 = 0.67;
  1248.  
  1249.     for (int i = 0; i < Node_count; i++)
  1250.     {
  1251.         /*SetCa(i, Ca0 + RandomD(-Ca0 * eps_persent, Ca0 * eps_persent)); // Ca
  1252.         SetIP3(i, IP30 + RandomD(-IP30 * eps_persent, IP30 * eps_persent)); // IP3
  1253.         Setz(i, z0 + RandomD(-z0 * eps_persent, z0 * eps_persent)); // z */
  1254.         SetCa(i, Ca0); // Ca
  1255.         SetIP3(i, IP30); // IP3
  1256.         Setz(i, z0); // z
  1257.     }
  1258.  
  1259.     double V0 = -58.7085;
  1260.     double m0 = 0.0953;
  1261.     double n0 = 0.000913;
  1262.     double h0 = 0.3662;
  1263.  
  1264.     double V1 = 14.8409;
  1265.     double m1 = 0.9174;
  1266.     double n1 = 0.0140;
  1267.     double h1 = 0.0539;
  1268.  
  1269.     /*for (int i = 0; i < Node_count; i++) // init only for neurons
  1270.     {
  1271.       double random = RandomD(0, 1);
  1272.  
  1273.       SetV(i, random < percent_stable_state ? V0 + RandomD(-V0 * eps_persent, V0 * eps_persent) : V1 + RandomD(-V1 * eps_persent, V1 * eps_persent)); // V
  1274.       Setm(i, random < percent_stable_state ? m0 + RandomD(-m0 * eps_persent, m0 * eps_persent) : m1 + RandomD(-m1 * eps_persent, m1 * eps_persent)); // m
  1275.       Setn(i, random < percent_stable_state ? n0 + RandomD(-n0 * eps_persent, n0 * eps_persent) : n1 + RandomD(-n1 * eps_persent, n1 * eps_persent)); // n
  1276.       Seth(i, random < percent_stable_state ? h0 + RandomD(-h0 * eps_persent, h0 * eps_persent) : h1 + RandomD(-h1 * eps_persent, h1 * eps_persent)); // h
  1277.     }*/
  1278.  
  1279.     for (int i = 0; i < Node_count; i++) // init only for neurons
  1280.     {
  1281.         double random = RandomD(0, 1);
  1282.  
  1283.         /*SetV(i, random < percent_stable_state ? V0 : V1); // V
  1284.         Setm(i, random < percent_stable_state ? m0 : m1); // m
  1285.         Setn(i, random < percent_stable_state ? n0 : n1); // n
  1286.         Seth(i, random < percent_stable_state ? h0 : h1); // h*/
  1287.  
  1288.         SetV(i, V0); // V
  1289.         Setm(i, m0); // m
  1290.         Setn(i, n0); // n
  1291.         Seth(i, h0); // h
  1292.  
  1293.         /*SetV_P(i, random < percent_stable_state ? V0 : V1); // V
  1294.         Setm_P(i, random < percent_stable_state ? m0 : m1); // m
  1295.         Setn_P(i, random < percent_stable_state ? n0 : n1); // n
  1296.         Seth_P(i, random < percent_stable_state ? h0 : h1); // h*/
  1297.  
  1298.         SetV_P(i, V0); // V
  1299.         Setm_P(i, m0); // m
  1300.         Setn_P(i, n0); // n
  1301.         Seth_P(i, h0); // h
  1302.     }
  1303.  
  1304.     /*for (int i = 0; i < Node_count; i++) // init only for neurons
  1305.     {*/
  1306.     /*SetV(i, RandomD(-80, 20)); // V
  1307.     Setm(i, RandomD(0, 1)); // m
  1308.     Setn(i, RandomD(0, 1)); // n
  1309.     Seth(i, RandomD(0, 1)); // h*/
  1310.     /*SetV(i, V1); // V
  1311.     Setm(i, m1); // m
  1312.     Setn(i, n1); // n
  1313.     Seth(i, h1); // h
  1314. }*/
  1315.  
  1316.     double E_syn0 = 0; // Excitatory neuron
  1317.     double E_syn1 = -90; // Inhibitory neuron
  1318.  
  1319.     fp_Esyn = fopen("results_E_syn.txt", "w+");
  1320.     for (int i = 0; i < Node_count; i++)
  1321.     {
  1322.         /*E_syn[i] = E_syn0;
  1323.  
  1324.         double x = RandomD(0, 1);
  1325.  
  1326.         if (x > p_inhib)
  1327.         {
  1328.             fprintf(fp_Esyn, "%f\n", E_syn[i]);
  1329.             continue;
  1330.         }*/
  1331.  
  1332.         E_syn[i] = E_syn1;
  1333.  
  1334.         fprintf(fp_Esyn, "%f\n", E_syn[i]);
  1335.  
  1336.         E_syn_P[i] = E_syn0;
  1337.     }
  1338.     fclose(fp_Esyn);
  1339.  
  1340.     for (int i = 0; i < Node_count; i++)
  1341.     {
  1342.         GenerateRandomMeander(i, 0);
  1343.         last_meander_end[i] = Meander_start_from_zero[i] + Duration;
  1344.     }
  1345.  
  1346.     const double t_start = 0;
  1347.     const double t_max = 60; // 100 msec = 0.1 sec // 240 // 270
  1348.     const double dt = 0.00002; // 0.01 msec = 0.00001 sec; 0.1 msec = 0.0001 sec; 1 msec = 0.001 sec // 0.000025
  1349.  
  1350.     double t = t_start;
  1351.  
  1352.     //fp_res = fopen("matlab_res.txt", "a+");
  1353.     //fprintf(fp_res, "%f\t%f\t", Max_magintude, g_syn);
  1354.  
  1355.     //fp0 = fopen("results.txt", "w+");
  1356.     //setlocale(LC_NUMERIC, "French_Canada.1252");
  1357.  
  1358.     double start_rk4, end_rk4;
  1359.     //clock_t start_rk4, end_rk4;
  1360.     //start_rk4 = omp_get_wtime();
  1361.     start_rk4 = clock();
  1362.     int lastPercent = -1;
  1363.  
  1364.     //FillVOldFromCurrent();
  1365.  
  1366.     k = new double* [Equations_count];
  1367.     for (int i = 0; i < Equations_count; i++)
  1368.         k[i] = new double[4];
  1369.  
  1370.     phi_k1 = new double[Equations_count];
  1371.     phi_k2 = new double[Equations_count];
  1372.     phi_k3 = new double[Equations_count];
  1373.  
  1374.     fp_I_stim = fopen("results_I_stim.txt", "w+");
  1375.     //fp_I_syn = fopen("results_I_syn.txt", "w+");
  1376.     fp_Ca = fopen("results_Ca.txt", "w+");
  1377.     //fp_IP3 = fopen("results_IP3.txt", "w+");
  1378.     //fp_z    = fopen("results_z.txt", "w+");
  1379.     //fp_G_P = fopen("results_G_P.txt", "w+");
  1380.     fp_V = fopen("results_V.txt", "w+");
  1381.     fp_V_P = fopen("results_V_P.txt", "w+");
  1382.     //fp_m = fopen("results_m.txt", "w+");
  1383.     //fp_n = fopen("results_n.txt", "w+");
  1384.     //fp_h = fopen("results_h.txt", "w+");
  1385.     fp_V_spikes = fopen("results_V_spikes.txt", "w+");
  1386.     //
  1387.  
  1388.     double* f_next = new double[Equations_count];
  1389.  
  1390.     while (t < t_max || Approximately(t, t_max))
  1391.     {
  1392.         fprintf(fp_I_stim, "%f\t", t);
  1393.         fprintf(fp_Ca, "%f\t", t);
  1394.         //fprintf(fp_IP3, "%f\t", t);
  1395.         //fprintf(fp_z, "%f\t", t);
  1396.         //fprintf(fp_G, "%f\t", t);
  1397.         fprintf(fp_V, "%f\t", t);
  1398.         fprintf(fp_V_P, "%f\t", t);
  1399.         //fprintf(fp_m, "%f\t", t);
  1400.         //fprintf(fp_n, "%f\t", t);
  1401.         //fprintf(fp_h, "%f\t", t);
  1402.         fprintf(fp_V_spikes, "%f\t", t);
  1403.  
  1404.         for (int i = 0; i < Node_count; i++)
  1405.         {
  1406.             if (t > last_meander_end[i])
  1407.             {
  1408.                 GenerateRandomMeander(i, t);
  1409.                 last_meander_end[i] = Meander_start_from_zero[i] + Duration;
  1410.             }
  1411.  
  1412.             fprintf(fp_I_stim, "%f\t", I_stim(i, t));
  1413.         }
  1414.         fprintf(fp_I_stim, "\n");
  1415.  
  1416.         for (int i = 0; i < Equations_count; i += Equations_per_node)
  1417.             fprintf(fp_Ca, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // Ca
  1418.  
  1419.         //for (int i = 1; i < Equations_count; i += Equations_per_node)
  1420.         //  fprintf(fp_IP3, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // IP3
  1421.  
  1422.         //for (int i = 2; i < Equations_count; i += Equations_per_node)
  1423.         //  fprintf(fp_z, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // z
  1424.  
  1425.         //for (int i = 3; i < Equations_count; i += Equations_per_node)
  1426.         //  fprintf(fp_G_P, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // G
  1427.  
  1428.         for (int i = 4; i < Equations_count; i += Equations_per_node)
  1429.         {
  1430.             if (isnan(f[i]))
  1431.                 return 1;
  1432.  
  1433.             fprintf(fp_V, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // V
  1434.         }
  1435.  
  1436.         //for (int i = 5; i < Equations_count; i += Equations_per_node)
  1437.         //  fprintf(fp_m, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // m
  1438.  
  1439.         //for (int i = 6; i < Equations_count; i += Equations_per_node)
  1440.         //  fprintf(fp_n, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // n
  1441.  
  1442.         //for (int i = 7; i < Equations_count; i += Equations_per_node)
  1443.         //  fprintf(fp_h, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // h
  1444.  
  1445.         for (int i = 8; i < Equations_count; i += Equations_per_node)
  1446.             fprintf(fp_V_P, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // V_P
  1447.  
  1448.         //for (int i = 9; i < Equations_count; i += Equations_per_node)
  1449.         //  fprintf(fp_m_P, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // m_P
  1450.  
  1451.         //for (int i = 10; i < Equations_count; i += Equations_per_node)
  1452.         //  fprintf(fp_n_P, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // n_P
  1453.  
  1454.         //for (int i = 11; i < Equations_count; i += Equations_per_node)
  1455.         //  fprintf(fp_h_P, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // h_P
  1456.  
  1457.         fprintf(fp_Ca, "\n");
  1458.         //fprintf(fp_IP3, "\n");
  1459.         //fprintf(fp_z, "\n");
  1460.         //fprintf(fp_G, "\n");
  1461.         fprintf(fp_V, "\n");
  1462.         //fprintf(fp_m, "\n");
  1463.         //fprintf(fp_n, "\n");
  1464.         //fprintf(fp_h, "\n");
  1465.  
  1466.         //fprintf(fp_G_P, "\n");
  1467.         fprintf(fp_V_P, "\n");
  1468.         //fprintf(fp_m_P, "\n");
  1469.         //fprintf(fp_n_P, "\n");
  1470.         //fprintf(fp_h_P, "\n");
  1471.  
  1472.         RungeKutta(t, dt, f, f_next);
  1473.  
  1474. #pragma omp parallel for
  1475.         for (int i = 0; i < Node_count; i++)
  1476.         {
  1477.             int index = Equations_per_node * i + 4;
  1478.             double diff = f_next[index] - f[index];
  1479.  
  1480.             fprintf(fp_V_spikes, i == Equations_count - 1 ? "%d" : "%d\t", diff < 0 && f_diff[index] > 0 && f[index] > -10 && (V_spikes[i].size() == 0 || t - V_spikes[i].back() > 0.001) ? 1 : 0);
  1481.  
  1482.             if (diff < 0 && f_diff[index] > 0 && f[index] > -10 && (V_spikes[i].size() == 0 || t - V_spikes[i].back() > 0.001))
  1483.             {
  1484.                 V_spikes[i].push_back(t);
  1485.             }
  1486.  
  1487.             f_diff[index] = diff;
  1488.         }
  1489.  
  1490.         fprintf(fp_V_spikes, "\n");
  1491.  
  1492.         CopyArray(f_next, f, Equations_count);
  1493.  
  1494.         //printf("V(24) = %f\t V_old(24) = %f\n", f[24*4], V_old(24));
  1495.         //UpdateVOld();
  1496.  
  1497.         //fprintf(fp_I_syn, "\n");
  1498.  
  1499.         for (int i = 0; i < Node_count; i++)
  1500.         {
  1501.             if (ca_all_state[i] == Ca_state::Wait_for_up)
  1502.             {
  1503.                 if (Ca(i) > 0.3)
  1504.                 {
  1505.                     ca_all_state[i] = Ca_state::Wait_for_down;
  1506.                     Interval interval;
  1507.                     interval.T_start = t;
  1508.                     interval.T_end = -1;
  1509.                     Ca_all_intervals[i].push_back(interval);
  1510.                 }
  1511.             }
  1512.             else if (ca_all_state[i] == Ca_state::Wait_for_down)
  1513.             {
  1514.                 if (Ca(i) < 0.3)
  1515.                 {
  1516.                     Ca_all_intervals[i].back().T_end = t;
  1517.                     ca_all_state[i] = Ca_state::Wait_for_up;
  1518.                 }
  1519.             }
  1520.         }
  1521.  
  1522.         t += dt;
  1523.  
  1524.         int percent = (int)(100 * (t - t_start) / (t_max - t_start));
  1525.         if (percent != lastPercent)
  1526.         {
  1527.             printf("Progress: %d%%\n", percent);
  1528.             lastPercent = percent;
  1529.         }
  1530.     }
  1531.  
  1532.     delete[] f_next;
  1533.  
  1534.     double* V_mean_freqs = new double[Node_count];
  1535.     double* V_STD_freqs = new double[Node_count];
  1536.     //list<double> V_mean_freqs;
  1537.  
  1538.     int max_interval_count = 0;
  1539.  
  1540.     for (int i = 0; i < Node_count; i++)
  1541.     {
  1542.         if (Ca_all_intervals[i].size() > max_interval_count)
  1543.         {
  1544.             max_interval_count = Ca_all_intervals[i].size();
  1545.         }
  1546.     }
  1547.  
  1548.     for (int i = 0; i < max_interval_count; i++)
  1549.     {
  1550.         Interval mean_interval;
  1551.         mean_interval.T_start = 0;
  1552.         mean_interval.T_end = 0;
  1553.  
  1554.         int mean_count = 0;
  1555.  
  1556.         for (int j = 0; j < Node_count; j++)
  1557.         {
  1558.             if (Ca_all_intervals[j].size() > i && Ca_all_intervals[j][i].T_end > 0)
  1559.             {
  1560.                 mean_interval.T_start += Ca_all_intervals[j][i].T_start;
  1561.                 mean_interval.T_end += Ca_all_intervals[j][i].T_end;
  1562.  
  1563.                 mean_count++;
  1564.             }
  1565.         }
  1566.  
  1567.         mean_interval.T_start /= mean_count;
  1568.         mean_interval.T_end /= mean_count;
  1569.  
  1570.         Ca_intervals.push_back(mean_interval);
  1571.     }
  1572.  
  1573. #pragma omp parallel for
  1574.     for (int i = 0; i < Node_count; i++)
  1575.     {
  1576.         list<double>::iterator it_V_spikes = V_spikes[i].begin();
  1577.  
  1578.         while (it_V_spikes != V_spikes[i].end() && *it_V_spikes < 0.25 * t_max)
  1579.         {
  1580.             V_spikes[i].pop_front();
  1581.             it_V_spikes = V_spikes[i].begin();
  1582.         }
  1583.  
  1584.         list<double> V_freqs;
  1585.  
  1586.         it_V_spikes = V_spikes[i].begin();
  1587.  
  1588.         V_mean_freqs[i] = 0;
  1589.         V_STD_freqs[i] = 0;
  1590.  
  1591.         if (V_spikes[i].size() <= 1)
  1592.         {
  1593.             continue;
  1594.         }
  1595.         else
  1596.         {
  1597.             for (int j = 1; j < V_spikes[i].size(); j++)
  1598.             {
  1599.                 double first = *it_V_spikes;
  1600.                 advance(it_V_spikes, 1);
  1601.                 double next = *it_V_spikes;
  1602.  
  1603.                 double T = next - first;
  1604.                 V_freqs.push_back(1 / T);
  1605.             }
  1606.         }
  1607.  
  1608.         list<double>::iterator it_Freq = V_freqs.begin();
  1609.  
  1610.         for (int j = 0; j < V_freqs.size(); j++)
  1611.         {
  1612.             V_mean_freqs[i] += *it_Freq;
  1613.             advance(it_Freq, 1);
  1614.         }
  1615.  
  1616.         V_mean_freqs[i] /= V_freqs.size();
  1617.  
  1618.         it_Freq = V_freqs.begin();
  1619.  
  1620.         for (int j = 0; j < V_freqs.size(); j++)
  1621.         {
  1622.             V_STD_freqs[i] += pow(*it_Freq - V_mean_freqs[i], 2);
  1623.             advance(it_Freq, 1);
  1624.         }
  1625.  
  1626.         V_STD_freqs[i] /= V_freqs.size();
  1627.         V_STD_freqs[i] = sqrt(V_STD_freqs[i]);
  1628.     }
  1629.  
  1630.     fp0 = fopen("V_STD_freqs.txt", "w+");
  1631.     for (int i = 0; i < Node_count; i++)
  1632.     {
  1633.         fprintf(fp0, "%f\t", V_STD_freqs[i]);
  1634.     }
  1635.     fclose(fp0);
  1636.  
  1637.     double V_STD_mean_Freq = 0;
  1638.     double V_STD_mean_Freq_not_zero = 0;
  1639.     double V_STD_mean_Freq_not_zero_count = 0;
  1640.  
  1641.     for (int j = 0; j < Node_count; j++)
  1642.     {
  1643.         if (V_STD_freqs[j] != 0)
  1644.         {
  1645.             V_STD_mean_Freq_not_zero_count++;
  1646.             V_STD_mean_Freq += V_STD_freqs[j];
  1647.         }
  1648.     }
  1649.  
  1650.     V_STD_mean_Freq_not_zero = V_STD_mean_Freq / V_STD_mean_Freq_not_zero_count;
  1651.     V_STD_mean_Freq /= Node_count;
  1652.  
  1653.     fp0 = fopen("V_STD_mean_Freq.txt", "w+");
  1654.     fprintf(fp0, "%f\n", V_STD_mean_Freq);
  1655.     //fprintf(fp_res, "%f\t", V_STD_mean_Freq);
  1656.     fclose(fp0);
  1657.  
  1658.     fp0 = fopen("V_STD_mean_Freq_not_zero.txt", "w+");
  1659.     fprintf(fp0, "%f\n", V_STD_mean_Freq_not_zero);
  1660.     //fprintf(fp_res, "%f\t", V_STD_mean_Freq_not_zero);
  1661.     fclose(fp0);
  1662.  
  1663.     double V_mean_mean_Freq = 0;
  1664.     double V_mean_mean_Freq_not_zero = 0;
  1665.     double V_mean_mean_Freq_not_zero_count = 0;
  1666.  
  1667.     for (int j = 0; j < Node_count; j++)
  1668.     {
  1669.         if (V_mean_freqs[j] != 0)
  1670.         {
  1671.             V_mean_mean_Freq_not_zero_count++;
  1672.             V_mean_mean_Freq += V_mean_freqs[j];
  1673.         }
  1674.     }
  1675.  
  1676.     delete[] V_mean_freqs;
  1677.     delete[] V_STD_freqs;
  1678.  
  1679.     V_mean_mean_Freq_not_zero = V_mean_mean_Freq / V_mean_mean_Freq_not_zero_count;
  1680.     V_mean_mean_Freq /= Node_count;
  1681.  
  1682.     fp0 = fopen("V_mean_mean_Freq.txt", "w+");
  1683.     fprintf(fp0, "%f\n", V_mean_mean_Freq);
  1684.     fclose(fp0);
  1685.  
  1686.     fp0 = fopen("V_mean_mean_Freq_not_zero.txt", "w+");
  1687.     fprintf(fp0, "%f\n", V_mean_mean_Freq_not_zero);
  1688.     fclose(fp0);
  1689.  
  1690.     // CHUNKS
  1691.     double chunk_t_start = 0.25 * t_max;
  1692.     double chunk_t_step = 0.5;
  1693.     int chunk_step_count = (0.75 * t_max / chunk_t_step);
  1694.  
  1695.     double dt_chunk = 0.1 / V_mean_mean_Freq_not_zero; // 0.25
  1696.     int chunks_count = chunk_t_step / dt_chunk;
  1697.  
  1698.     double corr_aver_mean = 0, corr_aver_mean_not_zero = 0;
  1699.  
  1700.     FILE* fp_corr_not_zero = fopen("corr_not_zero.txt", "w+");
  1701.     FILE* fp_V_spikes_chunks = fopen("V_spikes_chunks.txt", "w+");
  1702.  
  1703.     vector<double> k_syn;
  1704.     vector<double> k_syn_time;
  1705.  
  1706.     for (int s = 0; s < chunk_step_count; s++)
  1707.     {
  1708.         int** V_spikes_chunks = new int* [Node_count];
  1709.  
  1710. #pragma omp parallel for
  1711.         for (int i = 0; i < Node_count; i++)
  1712.         {
  1713.             V_spikes_chunks[i] = new int[chunks_count];
  1714.  
  1715.             if (V_spikes[i].size() == 0)
  1716.                 for (int ch = 0; ch < chunks_count; ch++)
  1717.                     V_spikes_chunks[i][ch] = 0;
  1718.  
  1719.             double ch_start = chunk_t_start + chunk_t_step * s;
  1720.             double ch_end = ch_start + dt_chunk;
  1721.             list<double>::iterator currentSpike = V_spikes[i].begin();
  1722.  
  1723.             for (int ch = 0; ch < chunks_count; ch++)
  1724.             {
  1725.                 while (currentSpike != V_spikes[i].end() && *currentSpike < ch_start)
  1726.                     currentSpike++;
  1727.  
  1728.                 if (currentSpike == V_spikes[i].end())
  1729.                     V_spikes_chunks[i][ch] = 0;
  1730.                 else
  1731.                     V_spikes_chunks[i][ch] = *currentSpike >= ch_start && *currentSpike <= ch_end;
  1732.  
  1733.                 ch_start += dt_chunk;
  1734.                 ch_end += dt_chunk;
  1735.             }
  1736.         }
  1737.  
  1738.         char buffer[50];
  1739.         //sprintf(buffer, "V_spikes_chunks_%d.txt", s);
  1740.         //fp0 = fopen(buffer, "w+");
  1741.  
  1742.         for (int ch = 0; ch < chunks_count; ch++)
  1743.         {
  1744.             for (int i = 0; i < Node_count; i++)
  1745.                 fprintf(fp_V_spikes_chunks, "%d\t", V_spikes_chunks[i][ch]);
  1746.  
  1747.             fprintf(fp_V_spikes_chunks, "\n");
  1748.         }
  1749.  
  1750.         //fclose(fp0);
  1751.  
  1752.         double corr = 0;
  1753.         double corr_not_zero = 0;
  1754.         double counter = 0;
  1755.         double counter_not_zero = 0;
  1756.  
  1757.         for (int i = 0; i < Node_count; i++)
  1758.         {
  1759.             for (int j = 0; j < Node_count; j++)
  1760.             {
  1761.                 if (i == j)
  1762.                     continue;
  1763.  
  1764.                 int k1 = 0;
  1765.                 int k2 = 0;
  1766.                 int k3 = 0;
  1767.  
  1768.                 for (int l = 0; l < chunks_count; l++)
  1769.                 {
  1770.                     if (V_spikes_chunks[i][l] == 1 && V_spikes_chunks[j][l] == 1)
  1771.                         k1++;
  1772.  
  1773.                     k2 += V_spikes_chunks[i][l];
  1774.                     k3 += V_spikes_chunks[j][l];
  1775.                 }
  1776.  
  1777.                 if (k2 != 0 && k3 != 0)
  1778.                 {
  1779.                     corr += (double)k1 / sqrt((double)k2 * (double)k3);
  1780.                     counter_not_zero++;
  1781.                 }
  1782.  
  1783.                 counter++;
  1784.             }
  1785.         }
  1786.  
  1787.         double corr_aver = corr / counter;
  1788.         double corr_aver_not_zero = corr / counter_not_zero;
  1789.  
  1790.         corr_aver_mean += corr_aver;
  1791.         corr_aver_mean_not_zero += corr_aver_not_zero;
  1792.  
  1793.         //sprintf(buffer, "corr_aver_%d.txt", s);
  1794.  
  1795.         //fp0 = fopen(buffer, "w+");
  1796.         //fprintf(fp0, "%f\n", corr_aver);
  1797.         //fclose(fp0);
  1798.  
  1799.         double ch_start = chunk_t_start + chunk_t_step * s;
  1800.         fprintf(fp_corr_not_zero, "%f\t%f\n", ch_start, corr_aver_not_zero);
  1801.  
  1802.         for (int i = 0; i < Node_count; i++)
  1803.             delete[] V_spikes_chunks[i];
  1804.  
  1805.         delete[] V_spikes_chunks;
  1806.  
  1807.         k_syn.push_back(corr_aver_not_zero);
  1808.         k_syn_time.push_back(ch_start);
  1809.     }
  1810.  
  1811.     fclose(fp_corr_not_zero);
  1812.     fclose(fp_V_spikes_chunks);
  1813.  
  1814.     corr_aver_mean /= chunk_step_count;
  1815.     corr_aver_mean_not_zero /= chunk_step_count;
  1816.  
  1817.     double k_syn_local_max_mean = 0;
  1818.     double k_syn_local_max_mean_count = 0;
  1819.  
  1820.     double k_syn_local_min_mean = 0;
  1821.     double k_syn_local_min_mean_count = 0;
  1822.  
  1823.     fp0 = fopen("k_syn_local_max_or_min.txt", "w+");
  1824.  
  1825.     for (int i = 0; i < Ca_intervals.size(); i++)
  1826.     {
  1827.         if (Ca_intervals[i].T_end < 0)
  1828.             break;
  1829.  
  1830.         if (g_astro > 0)
  1831.         {
  1832.             double k_syn_local_max = 0;
  1833.             bool firstValue = true;
  1834.  
  1835.             for (int j = 0; j < k_syn.size(); j++)
  1836.             {
  1837.                 if (k_syn_time[j] >= Ca_intervals[i].T_start && k_syn_time[j] <= Ca_intervals[i].T_end)
  1838.                 {
  1839.                     if (firstValue)
  1840.                     {
  1841.                         k_syn_local_max = k_syn[j];
  1842.                         firstValue = false;
  1843.                     }
  1844.                     else if (k_syn[j] > k_syn_local_max)
  1845.                     {
  1846.                         k_syn_local_max = k_syn[j];
  1847.                     }
  1848.                 }
  1849.             }
  1850.  
  1851.             k_syn_local_max_mean += k_syn_local_max;
  1852.             k_syn_local_max_mean_count++;
  1853.  
  1854.             fprintf(fp0, "%f\n", k_syn_local_max);
  1855.         }
  1856.         else if (g_astro < 0)
  1857.         {
  1858.             double k_syn_local_min = 1;
  1859.             bool firstValue = true;
  1860.  
  1861.             for (int j = 0; j < k_syn.size(); j++)
  1862.             {
  1863.                 if (k_syn_time[j] >= Ca_intervals[i].T_start && k_syn_time[j] <= Ca_intervals[i].T_end)
  1864.                 {
  1865.                     if (firstValue)
  1866.                     {
  1867.                         k_syn_local_min = k_syn[j];
  1868.                         firstValue = false;
  1869.                     }
  1870.                     else if (k_syn[j] < k_syn_local_min)
  1871.                     {
  1872.                         k_syn_local_min = k_syn[j];
  1873.                     }
  1874.                 }
  1875.             }
  1876.  
  1877.             k_syn_local_min_mean += k_syn_local_min;
  1878.             k_syn_local_min_mean_count++;
  1879.  
  1880.             fprintf(fp0, "%f\n", k_syn_local_min);
  1881.         }
  1882.     }
  1883.  
  1884.     fclose(fp0);
  1885.  
  1886.     k_syn_local_max_mean /= k_syn_local_max_mean_count;
  1887.     k_syn_local_min_mean /= k_syn_local_min_mean_count;
  1888.  
  1889.     //fp0 = fopen("corr_aver_mean.txt", "w+");
  1890.     //fprintf(fp0, "%f\n", corr_aver_mean);
  1891.     //fclose(fp0);
  1892.  
  1893.     //fp0 = fopen("corr_aver_mean_not_zero.txt", "w+");
  1894.     //fprintf(fp0, "%f\n", corr_aver_mean_not_zero);
  1895.     //fclose(fp0);
  1896.     ////// CHUNKS END
  1897.  
  1898.     //2
  1899.     double* V_freq_sync_time = new double[Node_count];
  1900.     double* V_freq_sync_time_relative = new double[Node_count];
  1901.  
  1902. #pragma omp parallel for
  1903.     for (int i = 0; i < Node_count; i++)
  1904.     {
  1905.         list<double>::iterator it_V_spikes = V_spikes[i].begin();
  1906.  
  1907.         while (it_V_spikes != V_spikes[i].end() && *it_V_spikes < 0.25 * t_max)
  1908.         {
  1909.             V_spikes[i].pop_front();
  1910.             it_V_spikes = V_spikes[i].begin();
  1911.         }
  1912.  
  1913.         list<double> V_freqs_normalized;
  1914.         list<double> V_freqs_time;
  1915.  
  1916.         it_V_spikes = V_spikes[i].begin();
  1917.  
  1918.         if (V_spikes[i].size() <= 1)
  1919.         {
  1920.             continue;
  1921.         }
  1922.         else
  1923.         {
  1924.             for (int j = 1; j < V_spikes[i].size(); j++)
  1925.             {
  1926.                 double first = *it_V_spikes;
  1927.                 advance(it_V_spikes, 1);
  1928.                 double next = *it_V_spikes;
  1929.  
  1930.                 double T = next - first;
  1931.                 V_freqs_normalized.push_back(1 / T - V_mean_mean_Freq);
  1932.                 V_freqs_time.push_back(next);
  1933.             }
  1934.         }
  1935.  
  1936.         list<double>::iterator it_Freq_normalized = V_freqs_normalized.begin();
  1937.         list<double>::iterator it_Freq_time = V_freqs_time.begin();
  1938.  
  1939.         V_freq_sync_time[i] = 0;
  1940.  
  1941.         for (int j = 1; j < V_freqs_normalized.size(); j++)
  1942.         {
  1943.             double Freq_normalized_last = *it_Freq_normalized;
  1944.             advance(it_Freq_normalized, 1);
  1945.             double Freq_normalized_next = *it_Freq_normalized;
  1946.  
  1947.             double Freq_time_last = *it_Freq_time;
  1948.             advance(it_Freq_time, 1);
  1949.             double Freq_time_next = *it_Freq_time;
  1950.  
  1951.             if (abs(Freq_normalized_last) <= 0.5 && abs(Freq_normalized_next) <= 0.5 && (Freq_time_next - Freq_time_last) <= 0.035)
  1952.                 V_freq_sync_time[i] += Freq_time_next - Freq_time_last;
  1953.         }
  1954.  
  1955.         V_freq_sync_time_relative[i] = V_freq_sync_time[i] / (t_max - (0.25 * t_max));
  1956.     }
  1957.  
  1958.     double V_mean_freq_sync_time_relative = 0;
  1959.  
  1960.     for (int j = 0; j < Node_count; j++)
  1961.     {
  1962.         V_mean_freq_sync_time_relative += V_freq_sync_time_relative[j];
  1963.     }
  1964.  
  1965.     V_mean_freq_sync_time_relative /= Node_count;
  1966.  
  1967.     fp0 = fopen("V_freq_sync_time_relative.txt", "w+");
  1968.     for (int i = 0; i < Node_count; i++)
  1969.     {
  1970.         fprintf(fp0, "%f\t", V_freq_sync_time_relative[i]);
  1971.     }
  1972.     fclose(fp0);
  1973.  
  1974.     fp0 = fopen("V_mean_freq_sync_time_relative.txt", "w+");
  1975.     fprintf(fp0, "%f\n", V_mean_freq_sync_time_relative);
  1976.     //fprintf(fp_res, "%f\n", V_mean_freq_sync_time_relative);
  1977.     fclose(fp0);
  1978.  
  1979.     //fclose(fp_res);
  1980.     fclose(fp_Max_magnitude);
  1981.     fclose(fp_I_stim);
  1982.     ///fclose(fp_I_syn);
  1983.     fclose(fp_Ca);
  1984.     //fclose(fp_IP3);
  1985.     //fclose(fp_z);
  1986.     fclose(fp_V);
  1987.     //fclose(fp_m);
  1988.     //fclose(fp_n);
  1989.     //fclose(fp_h);
  1990.     //fclose(fp_G_P);
  1991.     fclose(fp_V_P);
  1992.     //fclose(fp_m_P);
  1993.     //fclose(fp_n_P);
  1994.     //fclose(fp_h_P);
  1995.  
  1996.     fclose(fp_V_spikes);
  1997.  
  1998.     //end_rk4 = omp_get_wtime();
  1999.     end_rk4 = clock();
  2000.     double extime_rk4 = (double)(end_rk4 - start_rk4);// / CLOCKS_PER_SEC;
  2001.     int minutes = (int)extime_rk4 / 60;
  2002.     int seconds = (int)extime_rk4 % 60;
  2003.     printf("\nExecution time is: %d minutes %d seconds\n ", minutes, seconds);
  2004.  
  2005.     /*int nth;
  2006.     #pragma omp parallel
  2007.     {
  2008.         #pragma omp master
  2009.         nth = omp_get_num_threads();
  2010.     }*/
  2011.  
  2012.     fp0 = fopen("time_exec.txt", "a");
  2013.     //fprintf(fp0, "%d %lf\n", nth, extime_rk4);
  2014.     fprintf(fp0, "%lf\n", extime_rk4);
  2015.     fclose(fp0);
  2016.  
  2017.     if (g_astro > 0)
  2018.     {
  2019.         fp0 = fopen("results_Poisson_Freq_g_syn_k_syn.txt", "a");
  2020.         fprintf(fp0, "%f\t%f\t%f\n", Poisson_Freq, g_syn, k_syn_local_max_mean);
  2021.         fclose(fp0);
  2022.     }
  2023.     else if (g_astro < 0)
  2024.     {
  2025.         fp0 = fopen("results_Poisson_Freq_g_syn_k_syn.txt", "a");
  2026.         fprintf(fp0, "%f\t%f\t%f\n", Poisson_Freq, g_syn, k_syn_local_min_mean);
  2027.         fclose(fp0);
  2028.     }
  2029.     else
  2030.     {
  2031.         fp0 = fopen("results_Poisson_Freq_g_syn_k_syn.txt", "a");
  2032.         fprintf(fp0, "%f\t%f\t%f\n", Poisson_Freq, g_syn, corr_aver_mean_not_zero);
  2033.         fclose(fp0);
  2034.     }
  2035.  
  2036.     for (int i = 0; i < Node_count; i++)
  2037.         delete[] A_A[i];
  2038.  
  2039.     delete[] A_A;
  2040.  
  2041.     for (int i = 0; i < Node_count; i++)
  2042.         delete[] B_A[i];
  2043.  
  2044.     delete[] B_A;
  2045.  
  2046.     delete[] C_A;
  2047.  
  2048.     for (int i = 0; i < Node_count; i++)
  2049.         delete[] A_N[i];
  2050.  
  2051.     delete[] A_N;
  2052.  
  2053.     for (int i = 0; i < Node_count; i++)
  2054.         delete[] B_N[i];
  2055.  
  2056.     delete[] B_N;
  2057.  
  2058.     delete[] C_N;
  2059.  
  2060.     delete[] A_N_P;
  2061.  
  2062.     for (int i = 0; i < Node_count; i++)
  2063.         delete[] B_N_P[i];
  2064.  
  2065.     delete[] B_N_P;
  2066.  
  2067.     delete[] C_N_P;
  2068.  
  2069.     delete[] f;
  2070.     delete[] f_diff;
  2071.     delete[] v_4;
  2072.     delete[] E_syn;
  2073.     delete[] I_app;
  2074.     delete[] V_spikes;
  2075.     delete[] V_spikes_Freq;
  2076.  
  2077.     delete[] Meander_start_from_zero;
  2078.     delete[] Meander_width;
  2079.     delete[] Meander_height;
  2080.     delete[] Meander_interval;
  2081.     delete[] last_meander_end;
  2082.     delete[] V_freq_sync_time;
  2083.  
  2084.     for (int i = 0; i < Node_count; i++)
  2085.         delete[] tau[i];
  2086.  
  2087.     delete[] tau;
  2088.  
  2089.     for (int i = 0; i < Node_count; i++)
  2090.         delete[] V_old_array[i];
  2091.  
  2092.     delete[] V_old_array;
  2093.  
  2094.     for (int i = 0; i < Equations_count; i++)
  2095.         delete[] k[i];
  2096.  
  2097.     delete[] k;
  2098.  
  2099.     delete[] phi_k1;
  2100.     delete[] phi_k2;
  2101.     delete[] phi_k3;
  2102. }
  2103.  
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!