#define _USE_MATH_DEFINES
#include "math.h"
#include <stdlib.h>
#include <stdio.h>
#include <locale.h>
#include <time.h>
#include <stdbool.h>

#define Node_count 10*10*2 // должно быть таким, что бы из этого числа, разделенного на два, извлекался корень
#define Node_count_half Node_count / 2

#define Equations_per_node 4
#define Equations_count Node_count * Equations_per_node

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 g_syn = 0.05; // 0.1
double k_syn = 0.2; // 0.2
double E_syn[Node_count];

double I_app[Node_count];

const double sigma_G  = 0.0;
const double sigma_GN = 0.0;
const double sigma_N  = 1.0;

double** A;
double** B;
double*  C;
double** sigma;

int V_old_min = 0;
int V_old_max = 2000; // 700 steps = 7 ms delay;

int V_old_length[Node_count];
double** V_old_array;

int V_old_offset = 0;

const double Min_freq      = 50;
const double Max_freq      = 500;
const double Min_magintude = -1.5;
const double Max_magintude =  1.5;
const double Duration      = 0.001;

double Meander_start_from_zero[Node_count];
double Meander_width[Node_count];
double Meander_height[Node_count];
double Meander_interval[Node_count];

double I_stim(double t, int i)
{
    if (t < Meander_start_from_zero[i])
        return 0;

    t -= Meander_start_from_zero[i];
    t = fmod(t, Meander_width[i] + Meander_interval[i]);

    return t < Meander_width[i] ? Meander_height[i] : 0;
}

double V(int i)
{
	return f[i * 4];
}

void SetV(int i, double value)
{
	f[i * 4] = value;
}

double m(int i)
{
	return f[i * 4 + 1];
}

void Setm(int i, double value)
{
	f[i * 4 + 1] = value;
}

double n(int i)
{
	return f[i * 4 + 2];
}

void Setn(int i, double value)
{
	f[i * 4 + 2] = value;
}

double h(int i)
{
	return f[i * 4 + 3];
}

void Seth(int i, double value)
{
	f[i * 4 + 3] = value;
}

double V_old(int i)
{
	if (V_old_offset < V_old_length[i])
		return V(i);

	return V_old_array[i][0];
}

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_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_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_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 < Node_count; j++)
		  //{
		  //sum += /*sigma[in][j] * */A[in][j] * g_syn * (V(in) - V(j));
		  //sum += /*sigma[in][j] * */A[in][j] * g_syn * (V(in) - E_syn[in]) / (1 + exp(-V_old(j) / k_syn));
		  //}

		  /*for (int j = 0; j < C[in]; j++)
		  {
		    sum += sigma[in][(int)B[in][j]] * (V((int)B[in][j]) - V(in));
		  }*/

		  for (int j = 0; j < C[in]; j++)
		  {
			  // только тут подставил V_old
			  sum += /*sigma[in][(int)B[in][j]] * */ g_syn * (V(in) - E_syn[in]) / (1 + exp(-V_old((int)B[in][j]) / k_syn));
			  //sum += g_syn * (V((int)B[in][j]) - V(in));
		  }
		  //printf("i = %d\t sum = %f\n", in, sum);

		  return 1000 * ((g_Na * pow(m(in), 3) * h(in) * (E_Na - V(in)) + g_K * pow(n(in), 4) * (E_K - V(in)) + g_L * (E_L - V(in)) + I_app[in] + sum) / C_m); // V
	  }

	  case 1:
	  {
		  return 1000 * (alpha_m(f, in) * (1 - m(in)) - beta_m(f, in) * m(in));
	  }

	  case 2:
	  {
		  return 1000 * (alpha_n(f, in) * (1 - n(in)) - beta_n(f, in) * n(in));
	  }

	  case 3:
	  {
		  return 1000 * (alpha_h(f, in) * (1 - h(in)) - beta_h(f, in) * h(in));
	  }
	}

	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 < Node_count; i++)
	{
		for (int j = 0; j < Node_count; j++)
		{
			A[i][j] = 0;
		}
	}
}

void FillBCMatrix()
{
	for (int i = 0; i < Node_count; i++)
	{
		int bIndex = 0;
		C[i] = 0;
		for (int j = 0; j < Node_count; j++)
		{
			if (A[i][j] == 1)
			{
				B[i][bIndex] = j;
				bIndex++;
				C[i]++;
			}
		}
	}
}

void FillSigmaMatrix()
{
	for (int i = 0; i < Node_count; i++)
	{
		for (int j = 0; j < Node_count; j++)
		{
			if ((i >= 0 && i < Node_count_half) && (j >= 0 || j < Node_count_half))
			{
				sigma[i][j] = sigma_G;
			}
			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))))
			{
				sigma[i][j] = sigma_GN;
			}
			if ((i >= Node_count_half && i < Node_count) && (j >= Node_count_half && j < Node_count))
			{
				sigma[i][j] = sigma_N;
			}
		}
	}
}

void FillRandomMatrix()
{
	for (int k = 0; k < 2 * Node_count_half; k++)
	{
		int i, j;

		do
		{
			i = RandomI(Node_count_half, Node_count);
			j = RandomI(Node_count_half, Node_count);
		} while ((i == j) || (A[i][j] == 1));

		A[i][j] = 1;
		A[j][i] = 1;
	}
}

void FillVOldFromCurrent()
{
	for (int i = 0; i < Node_count; i++)
		for (int j = 0; j < V_old_length[i]; j++)
			V_old_array[i][j] = V(i);
}

void UpdateVOld()
{
	for (int i = 0; i < Node_count; i++)
	{
		for (int j = 1; j < V_old_length[i]; j++)
			V_old_array[i][j - 1] = V_old_array[i][j];

		V_old_array[i][V_old_length[i] - 1] = V(i);
	}

	V_old_offset++;
}

int main(int argc, char *argv[])
{
	FILE *fp0;
	srand(time(NULL));

	for (int i = 0; i < Node_count; i++)
		V_old_length[i] = RandomI(V_old_min, V_old_max);

	V_old_array = malloc(Node_count * sizeof(double*));
	for (int i = 0; i < Node_count; i++)
		V_old_array[i] = malloc(V_old_length[i] * sizeof(double));

	A = malloc(Node_count * sizeof(double));
	for (int i = 0; i < Node_count; i++)
		A[i] = malloc(Node_count * sizeof(double));

	B = malloc(Node_count * sizeof(double));
	for (int i = 0; i < Node_count; i++)
		B[i] = malloc(Node_count * sizeof(double));

	C = malloc(Node_count * sizeof(double));

	sigma = malloc(Node_count * sizeof(double));
	for (int i = 0; i < Node_count; i++)
		sigma[i] = malloc(Node_count * sizeof(double));

	FillAMatrixZero();
	FillRandomMatrix();
	FillBCMatrix();
	FillSigmaMatrix();

	fp0 = fopen("A.txt", "w+");
	for (int i = 0; i < Node_count; i++)
	{
		for (int j = 0; j < Node_count; j++)
		{
			fprintf(fp0, "%d\t", (int)A[i][j]);
		}
		fprintf(fp0, "\n");
	}
	fclose(fp0);

	//Пишем в файл число связей у каждого осциллятора в четвертом квадранте
	fp0 = fopen("links.txt", "w+");
	for (int i = Node_count_half; i < Node_count; i++)
	{
		int links_count = 0;
		for (int j = Node_count_half; j < Node_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 < Node_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 < Node_count; i++)
	{
		fprintf(fp0, "%d\n", (int)C[i]);
	}
	fclose(fp0);

	fp0 = fopen("sigma.txt", "w+");
	for (int i = 0; i < Node_count; i++)
	{
		for (int j = 0; j < Node_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(9, 40);
	}*/

	double percent_stable_state = 0.45;

	double V0 = -61.5364;
	double m0 = 0.0789;
	double n0 = 0.3718;
	double h0 = 0.4723;

	double V1 = 34.3334;
	double m1 = 0.9165;
	double n1 = 0.5626;
	double h1 = 0.2440;

	for (int i = 0; i < Node_count; i++) // init array for all neurons
	{
		SetV(i, 0); // V
		Setm(i, 0); // m
		Setn(i, 0); // n
		Seth(i, 0); // h
	}
	for (int i = Node_count_half; i < Node_count; i++) // init only neuron nodes
	{
		double random = RandomD(0, 1);

		SetV(i, random < percent_stable_state ? V0 : V1); // V
		Setm(i, random < percent_stable_state ? m0 : m1); // m
		Setn(i, random < percent_stable_state ? n0 : n1); // n
		Seth(i, random < percent_stable_state ? h0 : h1); // h
	}

	for (int i = 0; i < Node_count; i++)
	{
		I_app[i] = 0; // init for neurons all array
	}
	for (int i = Node_count_half; i < Node_count; i++)
	{
		I_app[i] = 5.27; //RandomD(5.20, 5.34); // init for neurons; Bifurcation point: I_app = 5.27
	}

	double percent_excitable = 1.0; // 0.8
	double E_syn0 = 0;
	double E_syn1 = -90;

	for (int i = 0; i < Node_count; i++)
	{
		E_syn[i] = 0; // init for neurons all array
	}
	for (int i = Node_count_half; i < Node_count; i++)
	{
		double random = RandomD(0, 1);
		E_syn[i] = random < percent_excitable ? E_syn0 : E_syn1;
		//printf("i = %d\t E_syn = %f\n", i, E_syn[i]);
	}

	const double t_start = 0;
	const double t_max   = 10; // 100 msec = 0.1 sec
	const double dt      = 0.00001; // 0.01 msec = 0.00001 sec; 1 msec = 0.001 sec

	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;

	FillVOldFromCurrent();

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

		//printf("V(24) = %f\t V_old(24) = %f\n", f[24*4], V_old(24));
		UpdateVOld();
	}

	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 < Node_count; i++)
		free(V_old_array[i]);
	free(V_old_array);

	for (int i = 0; i < Node_count; i++)
		free(A[i]);
	free(A);

	for (int i = 0; i < Node_count; i++)
		free(B[i]);
	free(B);

	for (int i = 0; i < Node_count; i++)
		free(sigma[i]);
	free(sigma);

	free(C);
}