#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 5*5*2
#define Node_count_half Node_count / 2

#define Equations_per_node 8
#define Equations_count Node_count * Equations_per_node

#define WireWidth (int)sqrt(Node_count_half)

double f[Equations_count];

double c_0 = 2; // uM
double c_1 = 0.185;
double v_1 = 6; // s^-1
double v_2 = 0.11; // s^-1
double v_3 = 2.2; // uM/s
double v_4[Equations_count]; // uM/s - Controling parameter
double v_5 = 0.025; // uM/s
double v_6 = 0.2;  // uM/s
double k_1 = 0.5; // s^-1
double k_2 = 1; // uM
double k_3 = 0.1;
double k_4 = 1.1; // uM/s
double a_2 = 0.14; // uM/s
double d_1 = 0.13; // uM
double d_2 = 1.049; // uM
double d_3 = 0.9434; // uM
double d_5 = 0.082; // uM
double alpha = 0.8;
double tau_IP3 = 7.143; // s
double IP3_star = 0.16; // uM
double d_Ca = 0.001;
double d_IP3 = 0.12;

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]; // pA - Controling neuron parameter
double g_syn = 0.1;
double g_astro = 3;
double k_syn = 0.2;

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

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

#define V_old_length 350
double** V_old_array;

int V_old_offset = 0;

double Ca(int i)
{
	return f[i * 8];
}

double IP3(int i)
{
	return f[i * 8 + 1];
}

double z(int i)
{
	return f[i * 8 + 2];
}

double G(int i)
{
	return f[i * 8 + 3];
}

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

double m(int i)
{
	return f[i * 8 + 5];
}

double n(int i)
{
	return f[i * 8 + 6];
}

double h(int i)
{
	return f[i * 8 + 7];
}

double V_old(int i)
{
	if (V_old_offset == 0)
		return V_old_array[0][i];

	if (V_old_offset < V_old_length)
		return V_old_array[V_old_length - V_old_offset][i];

	return V_old_array[0][i];
}

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 J_channel(double* f, int i)
{
	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);
}

double J_PLC(double* f, int i)
{
	return v_4[i] * (Ca(i) + (1 - alpha) * k_4) / (Ca(i) + k_4);
}

double J_leak(double* f, int i)
{
	return c_1 * v_2 * (c_0 / c_1 - (1 + 1 / c_1) * Ca(i));
}

double J_pump(double* f, int i)
{
	return v_3 * pow(Ca(i), 2) / (pow(k_3, 2) + pow(Ca(i), 2));
}

double J_in(double* f, int i)
{
	return v_5 + v_6 * pow(IP3(i), 2) / (pow(k_2, 2) + pow(IP3(i), 2));
}

double J_out(double* f, int i)
{
	return k_1 * Ca(i);
}

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 UllahJung_HodgkinHuxley(int i, double* f)
{
  int in = i / 8;
	int il = i % 8;

	switch (il)
	{
	  case 0: // Ca
    {
      double sum_1 = 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_1 += sigma[in][(int)B[in][j]] * d_Ca * (Ca((int)B[in][j]) - Ca(in));
      }

  		return J_channel(f, in) - J_pump(f, in) + J_leak(f, in) + J_in(f, in) - J_out(f, in) + sum_1; // Ca
    }

  	case 1: // IP3
    {
      double sum_2 = 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_2 += sigma[in][(int)B[in][j]] * d_IP3 * (IP3((int)B[in][j]) - IP3(in));
      }

  		return (IP3_star - IP3(in)) / tau_IP3 + J_PLC(f, in) + sum_2; // IP3
    }

  	case 2: // z
    {
  		return a_2 * (d_2 * (IP3(in) + d_1) / (IP3(in) + d_3) * (1 - z(in)) - Ca(in) * z(in)); // z
    }

    case 3: // G
    {
      return 1; // G
    }

    case 4: // V
  	{
  		double sum_3 = 0;

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

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

  		for (int j = 0; j < C[in]; j++)
  		{
  			// только тут подставил V_old
  			sum_3 += sigma[in][(int)B[in][j]] * g_syn * (1 + g_astro * Ca(in)) * V((int)B[in][j]) / (1 + exp(-V_old(in) / k_syn));
  		}

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

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

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

  	case 7: // h
  	{
  		return 1000 * (alpha_h(f, in) * (1 - h(in)) - beta_h(f, in) * h(in)); // h
  	}
	}
  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] = UllahJung_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] = UllahJung_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] = UllahJung_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] = UllahJung_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;
}

bool CheckSameLine(int i, int j)
{
	return i / WireWidth == j / WireWidth;
}

bool IsWireNeighbors(int i, int j)
{
	if (CheckSameLine(i, j) && (i == j - 1 || i == j + 1))
		return true;

	if (i == j - WireWidth || i == j + WireWidth)
		return true;

	return false;
}

void FillAMatrixZero()
{
	for (int i = 0; i < Node_count; i++)
	{
		for (int j = 0; j < Node_count; j++)
		{
			A[i][j] = 0;
		}
	}
}

void FillWireMatrix()
{
	for (int i = 0; i < Node_count_half; i++)
	{
		for (int j = 0; j < Node_count_half; j++)
		{
			if (i == j)
			{
				A[i][j] = 0;
				continue;
			}

			if (i > j)
			{
				A[i][j] = A[j][i];
				continue;
			}

			A[i][j] = IsWireNeighbors(i, j) ? 1 : 0;
		}
	}
}

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 Connect(int i, int j)
{
	A[i][j] = 1;
	A[j][i] = 1;
}

void FillLayerConnectionMatrix()
{
	//int min = 0;
	//int max = NHalf;
	//int layerOffset = NHalf;

	int min = Node_count_half;
	int max = Node_count;
	int layerOffset = -Node_count_half;

	for (int i = min; i < max; i++)
	{
		int nextLayerIndex = i + layerOffset;
		Connect(i, nextLayerIndex);

		int leftIndex = i - 1;
		if (leftIndex >= min && CheckSameLine(leftIndex, i))
			Connect(leftIndex, nextLayerIndex);

		int rightIndex = i + 1;
		if (rightIndex < max && CheckSameLine(rightIndex, i))
			Connect(rightIndex, nextLayerIndex);

		int upIndex = i - WireWidth;
		if (upIndex >= min)
			Connect(upIndex, nextLayerIndex);

		int downIndex = i + WireWidth;
		if (downIndex < max)
			Connect(downIndex, nextLayerIndex);
	}
}

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

void UpdateVOld()
{
	double* tmp = V_old_array[0];

	for (int i = 1; i < V_old_length; i++)
		V_old_array[i - 1] = V_old_array[i];

	for (int i = 0; i < Node_count; i++)
		tmp[i] = V(i);

	V_old_array[V_old_length - 1] = tmp;

	V_old_offset++;
}

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

  V_old_array = malloc(V_old_length * sizeof(double*));
  for (int i = 0; i < V_old_length; i++)
    V_old_array[i] = malloc(Node_count * 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();
	FillWireMatrix();
	FillRandomMatrix();
  FillLayerConnectionMatrix();
	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("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 at t = 0
  for (int i = 0; i < Equations_count; i += 8)
  {
    f[i] = 0.07; // Ca
  }

  for (int i = 1; i < Equations_count; i += 8)
  {
    f[i] = 0.16; // IP3
  }

  for (int i = 2; i < Equations_count; i += 8)
  {
    f[i] = 0.67; // z
  }

  for (int i = 3; i < Equations_count; i += 8)
	{
		f[i] = 0.0; // G - ????
	}

	for (int i = 2; i < Equations_count; i += 8)
	{
		v_4[i] = RandomD(0.5, 1.0); // init for astrocites
	}

	for (int i = 4; i < Equations_count; i += 8)
	{
		f[i] = 0; // V
	}

  for (int i = 5; i < Equations_count; i += 8)
	{
		f[i] = 0; // m
	}

  for (int i = 6; i < Equations_count; i += 8)
	{
		f[i] = 0; // n
	}

  for (int i = 7; i < Equations_count; i += 8)
	{
		f[i] = 0; // h
	}

	for (int i = 0; i < Equations_count; i++)
	{
		I_app[i] = RandomD(9, 20); // init for neurons
	}

	const double t_start = 0;
	const double t_max = 50; // 100 msec = 0.1 sec
	const double dt = 0.00001; // 0.01 msec = 0.00001 sec

	double t = t_start;

	/*fp0 = fopen("I_stim_height.txt", "a");
	fprintf(fp0, "%f\t", I_app);
	fclose(fp0);*/

	//fp0 = fopen("results.txt", "w+");
	//setlocale(LC_NUMERIC, "French_Canada.1252");

	clock_t start_rk4, end_rk4;
	start_rk4 = clock();
	int lastPercent = -1;

  FillVOldFromCurrent();

  fp_Ca = fopen("results_Ca.txt", "w+");
  fp_G = fopen("results_G.txt", "w+");
  fp_V = fopen("results_V.txt", "w+");

  while (t < t_max || Approximately(t, t_max))
	{
    fprintf(fp_Ca, "%f\t", t);
    fprintf(fp_G, "%f\t", t);
    fprintf(fp_V, "%f\t", t);

    for (int i = 0; i < Equations_count; i += 8)
			fprintf(fp_Ca, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // Ca

    for (int i = 3; i < Equations_count; i += 8)
    	fprintf(fp_G, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // G

    for (int i = 4; i < Equations_count; i += 8)
  		fprintf(fp_V, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // V

		fprintf(fp_Ca, "\n");
    fprintf(fp_G, "\n");
    fprintf(fp_V, "\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;
		}

    UpdateVOld();
	}

  fclose(fp_Ca);
  fclose(fp_G);
  fclose(fp_V);

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