Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #define _USE_MATH_DEFINES
- #include "math.h"
- #include <stdlib.h>
- #include <stdio.h>
- #include <locale.h>
- #include <time.h>
- #include <stdbool.h>
- #define Neuron_count 20
- #define Neuron_count_half Neuron_count / 2
- #define Equations_per_neuron 4
- #define Equations_count Neuron_count * Equations_per_neuron
- double f[Equations_count];
- double C_m = 1;
- double g_K = 36;
- double g_Na = 120;
- double g_L = 0.3;
- double E_K = -77;
- double E_Na = 55;
- double E_L = -54.4;
- double I_app[Equations_count];
- const double sigma_G = 0.3;
- const double sigma_GN = 0.5;
- const double sigma_N = 0.0;
- double** A;
- double** B;
- double* C;
- double** sigma;
- double V(int i)
- {
- return f[i * 4];
- }
- double m(int i)
- {
- return f[i * 4 + 1];
- }
- double n(int i)
- {
- return f[i * 4 + 2];
- }
- double h(int i)
- {
- return f[i * 4 + 3];
- }
- int RandomI(int min, int max)
- {
- return ((double)rand() / (RAND_MAX - 1)) * (max - min) + min;
- }
- double RandomD(double min, double max)
- {
- return ((double)rand() / RAND_MAX) * (max - min) + min;
- }
- double alpha_n(double* f, int i)
- {
- return 0.01 * (V(i) + 55) / (1 - exp(-(V(i) + 55) / 10));
- }
- double beta_n(double* f, int i)
- {
- return 0.125 * exp(-(V(i) + 65) / 80);
- }
- double alpha_m(double* f, int i)
- {
- return 0.1 * (V(i) + 40) / (1 - exp(-(V(i) + 40) / 10));
- }
- double beta_m(double* f, int i)
- {
- return 4 * exp(-(V(i) + 65) / 18);
- }
- double alpha_h(double* f, int i)
- {
- return 0.07 * exp(-(V(i) + 65) / 20);
- }
- double beta_h(double* f, int i)
- {
- return 1 / (exp(-(V(i) + 35) / 10) + 1);
- }
- double HodgkinHuxley(int i, double* f)
- {
- int in = i / 4;
- int il = i % 4;
- switch (il)
- {
- case 0:
- {
- double sum = 0;
- /*for (int j = 0; j < Neuron_count; j++)
- {
- sum += sigma[in][j] * A[in][j] * (V(j) - V(in));
- }*/
- for (int j = 0; j < C[in]; j++)
- {
- sum += sigma[in][(int)B[in][j]] * (V((int)B[in][j]) - V(in));
- }
- return (g_Na * m(il) * m(il) * m(il) * h(il) * (E_Na - V(il)) + g_K * n(il) * n(il) * n(il) * n(il) * (E_K - V(il)) + g_L * (E_L - V(il)) + I_app[il] + sum) / C_m;
- }
- case 1:
- return alpha_m(f, il) * (1 - m(il)) - beta_m(f, il) * m(il);
- case 2:
- return alpha_n(f, il) * (1 - n(il)) - beta_n(f, il) * n(il);
- case 3:
- return alpha_h(f, il) * (1 - h(il)) - beta_h(f, il) * h(il);
- }
- return 0;
- }
- void RungeKutta(double dt, double* f, double* f_next)
- {
- double k[Equations_count][4];
- // k1
- for (int i = 0; i < Equations_count; i++)
- k[i][0] = HodgkinHuxley(i, f) * dt;
- double phi_k1[Equations_count];
- for (int i = 0; i < Equations_count; i++)
- phi_k1[i] = f[i] + k[i][0] / 2;
- // k2
- for (int i = 0; i < Equations_count; i++)
- k[i][1] = HodgkinHuxley(i, phi_k1) * dt;
- double phi_k2[Equations_count];
- for (int i = 0; i < Equations_count; i++)
- phi_k2[i] = f[i] + k[i][1] / 2;
- // k3
- for (int i = 0; i < Equations_count; i++)
- k[i][2] = HodgkinHuxley(i, phi_k2) * dt;
- double phi_k3[Equations_count];
- for (int i = 0; i < Equations_count; i++)
- phi_k3[i] = f[i] + k[i][2] / 2;
- // k4
- for (int i = 0; i < Equations_count; i++)
- k[i][3] = HodgkinHuxley(i, phi_k3) * dt;
- for (int i = 0; i < Equations_count; i++)
- f_next[i] = f[i] + (k[i][0] + 2 * k[i][1] + 2 * k[i][2] + k[i][3]) / 6;
- }
- void CopyArray(double* source, double* target, int N)
- {
- for (int i = 0; i < N; i++)
- target[i] = source[i];
- }
- bool Approximately(double a, double b)
- {
- if (a < 0)
- a = -a;
- if (b < 0)
- b = -b;
- return a - b <= 0.000001;
- }
- void FillAMatrixZero()
- {
- for (int i = 0; i < Neuron_count; i++)
- {
- for (int j = 0; j < Neuron_count; j++)
- {
- A[i][j] = 0;
- }
- }
- }
- void FillBCMatrix()
- {
- for (int i = 0; i < Neuron_count; i++)
- {
- int bIndex = 0;
- C[i] = 0;
- for (int j = 0; j < Neuron_count; j++)
- {
- if (A[i][j] == 1)
- {
- B[i][bIndex] = j;
- bIndex++;
- C[i]++;
- }
- }
- }
- }
- void FillSigmaMatrix()
- {
- for (int i = 0; i < Neuron_count; i++)
- {
- for (int j = 0; j < Neuron_count; j++)
- {
- if ((i >= 0 && i < Neuron_count_half) && (j >= 0 || j < Neuron_count_half))
- {
- sigma[i][j] = sigma_G;
- }
- if ((((i >= Neuron_count_half && i < Neuron_count) && (j >= 0 && j < Neuron_count_half)) || ((i >= 0 && i < Neuron_count_half) && (j >= Neuron_count_half && j < Neuron_count))))
- {
- sigma[i][j] = sigma_GN;
- }
- if ((i >= Neuron_count_half && i < Neuron_count) && (j >= Neuron_count_half && j < Neuron_count))
- {
- sigma[i][j] = sigma_N;
- }
- }
- }
- }
- void FillRandomMatrix()
- {
- for (int k = 0; k < 2 * Neuron_count_half; k++)
- {
- int i, j;
- do
- {
- i = RandomI(Neuron_count_half, Neuron_count);
- j = RandomI(Neuron_count_half, Neuron_count);
- } while ((i == j) || (A[i][j] == 1));
- A[i][j] = 1;
- A[j][i] = 1;
- }
- }
- int main(int argc, char *argv[])
- {
- FILE *fp0;
- srand(time(NULL));
- A = malloc(Neuron_count * sizeof(double));
- for (int i = 0; i < Neuron_count; i++)
- A[i] = malloc(Neuron_count * sizeof(double));
- B = malloc(Neuron_count * sizeof(double));
- for (int i = 0; i < Neuron_count; i++)
- B[i] = malloc(Neuron_count * sizeof(double));
- C = malloc(Neuron_count * sizeof(double));
- sigma = malloc(Neuron_count * sizeof(double));
- for (int i = 0; i < Neuron_count; i++)
- sigma[i] = malloc(Neuron_count * sizeof(double));
- FillAMatrixZero();
- FillRandomMatrix();
- FillBCMatrix();
- FillSigmaMatrix();
- fp0 = fopen("A.txt", "w+");
- for (int i = 0; i < Neuron_count; i++)
- {
- for (int j = 0; j < Neuron_count; j++)
- {
- fprintf(fp0, "%d\t", (int)A[i][j]);
- }
- fprintf(fp0, "\n");
- }
- fclose(fp0);
- //Пишем в файл число связей у каждого осциллятора в четвертом квадранте
- fp0 = fopen("links.txt", "w+");
- for (int i = Neuron_count_half; i < Neuron_count; i++)
- {
- int links_count = 0;
- for (int j = Neuron_count_half; j < Neuron_count; j++)
- {
- if (A[i][j] == 1)
- {
- links_count++;
- }
- }
- fprintf(fp0, "%d\n", (int)links_count);
- }
- fclose(fp0);
- fp0 = fopen("B.txt", "w+");
- for (int i = 0; i < Neuron_count; i++)
- {
- for (int j = 0; j < C[i]; j++)
- {
- fprintf(fp0, "%d\t", (int)B[i][j]);
- }
- fprintf(fp0, "\n");
- }
- fclose(fp0);
- fp0 = fopen("C.txt", "w+");
- for (int i = 0; i < Neuron_count; i++)
- {
- fprintf(fp0, "%d\n", (int)C[i]);
- }
- fclose(fp0);
- fp0 = fopen("sigma.txt", "w+");
- for (int i = 0; i < Neuron_count; i++)
- {
- for (int j = 0; j < Neuron_count; j++)
- {
- fprintf(fp0, "%f\t", sigma[i][j]);
- }
- fprintf(fp0, "\n");
- }
- fclose(fp0);
- // Initial values
- for (int i = 0; i < Equations_count; i++)
- f[i] = 0;
- for (int i = 0; i < Equations_count; i++)
- I_app[i] = RandomD(10, 50);
- const double t_start = 0;
- const double t_max = 500;
- const double dt = 0.01;
- double t = t_start;
- fp0 = fopen("results.txt", "w+");
- //setlocale(LC_NUMERIC, "French_Canada.1252");
- clock_t start_rk4, end_rk4;
- start_rk4 = clock();
- int lastPercent = -1;
- while (t < t_max || Approximately(t, t_max))
- {
- fprintf(fp0, "%f\t", t);
- for (int i = 0; i < Equations_count; i += 4)
- fprintf(fp0, i == Equations_count - 1 ? "%f" : "%f\t", f[i]);
- fprintf(fp0, "\n");
- double f_next[Equations_count];
- RungeKutta(dt, f, f_next);
- CopyArray(f_next, f, Equations_count);
- t += dt;
- int percent = (int)(100 * (t - t_start) / (t_max - t_start));
- if (percent != lastPercent)
- {
- printf("Progress: %d%%\n", percent);
- lastPercent = percent;
- }
- }
- fclose(fp0);
- end_rk4 = clock();
- double extime_rk4 = (double)(end_rk4 - start_rk4) / CLOCKS_PER_SEC;
- int minutes = (int)extime_rk4 / 60;
- int seconds = (int)extime_rk4 % 60;
- printf("\nExecution time is: %d minutes %d seconds\n ", minutes, seconds);
- fp0 = fopen("time_exec.txt", "w+");
- fprintf(fp0, "%f\n", extime_rk4);
- fclose(fp0);
- for (int i = 0; i < Neuron_count; i++)
- free(A[i]);
- free(A);
- for (int i = 0; i < Neuron_count; i++)
- free(B[i]);
- free(B);
- for (int i = 0; i < Neuron_count; i++)
- free(sigma[i]);
- free(sigma);
- free(C);
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement