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

#define N 10*10*2 // должно быть таким, что бы из этого числа, разделенного на два, извлекался корень
#define NHalf N / 2
#define WireWidth (int)sqrt(NHalf) // число осциляторов в строке/столбце одного слоя
//#define K 1

const double omega_rand_glial_min = 0.5;
const double omega_rand_glial_max = 1.5;

const double omega_rand_neuron_min = 9.5;
const double omega_rand_neuron_max = 10.5;

const double theta_rand_min = 0;
const double theta_rand_max = 2 * M_PI;

const double sigma_G  = 0.5;
double sigma_GN;
const double sigma_N  = 0.5;

double** A;
double** B;
double*  C;
double*  omega;
double*  theta;
double** sigma;

/*double sinux_x[N];
double cosinus_x[N];
double sinus_y[N];
double cosinus_y[N];*/

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 kuramoto(int i, double theta[N]) // Без оптимизации
{
	double sum = 0;

	for (int j = 0; j < N; j++)
	{
		sum += A[i][j] * sin(theta[j] - theta[i]);
	}

	return omega[i] + (double)sigma[i] * sum;
}*/

double kuramoto(int i, double theta[N]) // С оптимизацией 1
{
	double sum = 0;

	for (int j = 0; j < C[i]; j++)
	{
		sum += sigma[i][(int)B[i][j]] * sin(theta[(int)B[i][j]] - theta[i]);
	}

	return omega[i] + sum;
}

/*double kuramoto(int i, double theta[N]) // С оптимизацией 2
{
	double sum = 0;
	for (int j = 0; j < C[i]; j++)
	{
		sinux_x[j] = sin(theta[(int)B[i][j]]);
		cosinus_x[j] = cos(theta[(int)B[i][j]]);
		sinus_y[j] = sin(theta[i]);
		cosinus_y[j] = cos(theta[i]);
	}

	for (int j = 0; j < C[i]; j++)
	{
		sum += sigma[i][(int)B[i][j]] * ( sinux_x[j] * cosinus_y[j] - cosinus_x[j] * sinus_y[j] );
	}

	return omega[i] + sum;
}*/

void RungeKutta(double dt, double theta[N], double theta_next[N])
{
	double k[N][4];

	// k1
	for (int i = 0; i < N; i++)
		k[i][0] = kuramoto(i, theta) * dt;

	double theta_k1[N];
	for (int i = 0; i < N; i++)
		theta_k1[i] = theta[i] + k[i][0] / 2;

	// k2
	for (int i = 0; i < N; i++)
		k[i][1] = kuramoto(i, theta_k1) * dt;

	double theta_k2[N];
	for (int i = 0; i < N; i++)
		theta_k2[i] = theta[i] + k[i][1] / 2;

	// k3
	for (int i = 0; i < N; i++)
		k[i][2] = kuramoto(i, theta_k2) * dt;

	double theta_k3[N];
	for (int i = 0; i < N; i++)
		theta_k3[i] = theta[i] + k[i][2] / 2;

	// k4
	for (int i = 0; i < N; i++)
		k[i][3] = kuramoto(i, theta_k3) * dt;

	for (int i = 0; i < N; i++)
		theta_next[i] = theta[i] + (k[i][0] + 2 * k[i][1] + 2 * k[i][2] + k[i][3]) / 6;
}

void CopyArray(double source[N], double target[N])
{
	for (int i = 0; i < N; i++)
		target[i] = source[i];
}

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

void FillGlialWireMatrix()
{
	for (int i = 0; i < NHalf; i++)
	{
		for (int j = 0; j < NHalf; 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 FilNetworklWireMatrix()
{
	for (int i = NHalf; i < N; i++)
	{
		for (int j = NHalf; j < N; 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 * NHalf; k++)
	{
		int i, j;

		do
		{
			i = RandomI(NHalf, N);
			j = RandomI(NHalf, N);
		} 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 = NHalf;
	int max = N;
	int layerOffset = -NHalf;

	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 FillFullNetworkNeuronMatrix()
{
	for (int i = NHalf; i < N; i++)
	{
		for (int j = NHalf; j < N; j++)
		{
			A[i][j] = 1;
			if (i == j)
			{
				A[i][j] = 0;
			}
		}
	}
}

void FillsigmaMatrix()
{
	for (int i = 0; i < N; i++)
	{
		for (int j = 0; j < N; j++)
		{
			if ((i >= 0 && i < NHalf) && (j >= 0 || j < NHalf))
			{
				sigma[i][j] = sigma_G;
			}
			if ((((i >= NHalf && i < N) && (j >= 0 && j < NHalf)) || ((i >= 0 && i < NHalf) && (j >= NHalf && j < N))))
			{
				sigma[i][j] = sigma_GN;
			}
      if ((i >= NHalf && i < N) && (j >= NHalf && j < N))
			{
				sigma[i][j] = sigma_N;
			}
		}
	}
}

bool Approximately(double a, double b)
{
	if (a < 0)
		a = -a;

	if (b < 0)
		b = -b;

	return a - b <= 0.000001;
}

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

int main()
{
	sigma_GN = sigma_G;

	FILE *fp0;
	srand(time(NULL));

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

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

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

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

	omega = malloc(N * sizeof(double));

	theta = malloc(N * sizeof(double));

	// omega_init open for write. begin
	for (int i = 0; i < NHalf; i++)
		omega[i] = RandomD(omega_rand_glial_min, omega_rand_glial_max);

	for (int i = NHalf; i < N; i++)
		omega[i] = RandomD(omega_rand_neuron_min, omega_rand_neuron_max);

	fp0 = fopen("omega_init.txt", "w+");
	for (int i = 0; i < N; i++)
	{
		fprintf(fp0, "%f\n", omega[i]);
	}
	fclose(fp0);
	// omega_init open for write. end

	// omega_init open for read. begin
	/*fp0 = fopen("omega_init.txt", "r");
	for (int i = 0; i < N; i++)
	{
		fscanf(fp0, "%lf", &omega[i]);
	}*/
	// omega_init open for read. end

	// theta_init open for write. begin
	for (int i = 0; i < N; i++)
		theta[i] = RandomD(theta_rand_min, theta_rand_max);

	fp0 = fopen("theta_init.txt", "w+");
	for (int i = 0; i < N; i++)
	{
		fprintf(fp0, "%f\n", theta[i]);
	}
	fclose(fp0);
	// theta_init open for write.end

	// theta_init open for read. begin
	/*fp0 = fopen("theta_init.txt", "r");
	for (int i = 0; i < N; i++)
	{
		fscanf(fp0, "%lf", &theta[i]);
	}
	fclose(fp0);*/
	// theta_init open for read.end

	fp0 = fopen("A.txt", "w+");

	FillAMatrixZero();
	FillGlialWireMatrix();
	FilNetworklWireMatrix();
	FillLayerConnectionMatrix();

	FillBCMatrix();

	for (int i = 0; i < N; i++)
	{
		for (int j = 0; j < N; j++)
		{
			/*printf("i: %d\n", i);
			printf("j: %d\n", j);
			printf("A: %d\n", A[i][j]);*/
			fprintf(fp0, "%d\t", (int)A[i][j]);
		}
		fprintf(fp0, "\n");
	}
	fclose(fp0);

	fp0 = fopen("B.txt", "w+");

	for (int i = 0; i < N; 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 < N; i++)
	{
		fprintf(fp0, "%d\n", (int)C[i]);
	}

	fclose(fp0);

  FillsigmaMatrix();

  fp0 = fopen("sigma.txt", "w+");
  for (int i = 0; i < N; i++)
	{
		for (int j = 0; j < N; j++)
		{
			fprintf(fp0, "%f\t", sigma[i][j]);
		}
		fprintf(fp0, "\n");
	}

	fclose(fp0);

	//Пишем в файл число связей у каждого осциллятора в четвертом квадранте
	/*fp0 = fopen("links.txt", "w+");
	for (int i = NHalf; i < N; i++)
	{
		int links_count = 0;
		for (int j = NHalf; j < N; j++)
		{
			if (A[i][j] == 1)
			{
				links_count++;
			}
		}
		fprintf(fp0, "%d\n", (int)links_count);

	}
	fclose(fp0);*/

	const double t_start = 0;
	const double t_max = 200; //2000
	const double dt = 0.05; //0.05

	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 < N; i++)
		{
			fprintf(fp0, i == N - 1 ? "%f" : "%f\t", theta[i]);
		}
		fprintf(fp0, "\n");

		double theta_next[N];

		RungeKutta(dt, theta, theta_next);
		CopyArray(theta_next, theta);

		t += dt;

		int percent = (int)(100 * (t - t_start) / (t_max - t_start));
		if (percent != lastPercent)
		{
			printf("Progress: %d%%\n", percent);
			lastPercent = percent;
		}
	}

	fp0 = fopen("s_G_s_N_Delta_rho.txt", "a");
  fprintf(fp0, "%f\t%f\t", sigma_G, sigma_N);
  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);

	fclose(fp0);

	fp0 = fopen("time_exec.txt", "w+");
	fprintf(fp0, "%f\n", extime_rk4);
	fclose(fp0);

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

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

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

	free(C);

	free(omega);

	free(theta);

}