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

#define N 50*50*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;

double sigma_G;
double sigma_GN;
double sigma_N;

double sigma_G_1;
double sigma_N_1;

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

double p = 0.3; // 0.1 rewiring neurons parameter
double p_in = 0.2; // 0.2 inhibitory neurons percent

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

// возвращает целое случайное число в интервале от min до max включительно
int RandomI(int min, int max)
{
    return round(((double)rand() / (RAND_MAX)) * (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++)
    {
        if (A[i][j] != 0)
            sum += A[i][j] * sigma[i][j] * sin(theta[j] - theta[i]);
    }

    return omega[i] + sum;
}*/

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

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

int GetLineIndex(int i)
{
    return i / WireWidth;
}

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

bool CheckAdjacentLine(int i, int j)
{
    return abs(GetLineIndex(i) - GetLineIndex(j)) == 1;
}

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;

    if (CheckAdjacentLine(i, j))
    {
        if (i == j - WireWidth - 1)
            return true;

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

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

        if (i == j + WireWidth + 1)
            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 FillNetworkWireMatrix()
{
    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 RandomizeWireMatrix()
{
    srand(time(NULL));

    for (int i = NHalf; i < N; i++)
    {
        if (i == N - NHalf / 2)
            srand(time(NULL));

        for (int j = NHalf; j < N; j++)
        {
            if (A[i][j] == 0)
                continue;

            if (i > j)
                continue;

            double x = RandomD(0, 1);

            if (x > p)
                continue;

            int rndJ;
            int rndI;

            do
            {
                rndJ = RandomI(NHalf, N - 1);
                rndI = RandomI(NHalf, N - 1);
            } while (rndJ == rndI || IsWireNeighbors(rndI, rndJ) || j == rndJ || i == rndI || A[rndI][rndJ] == 1);

            A[i][j] = 0;
            A[j][i] = 0;

            A[rndI][rndJ] = 1;
            A[rndJ][rndI] = 1;
        }
    }
}

void RandomizeInhibitoryWireMatrix()
{
    for (int i = NHalf; i < N; i++)
    {
        //if (A[i][0] == 0)
        //  continue;

        double x = RandomD(0, 1);

        if (x > p_in)
            continue;

        for (int j = NHalf; j < N; j++)
        {
            if (A[i][j] == 1)
                A[i][j] = -1;
        }
    }
}

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

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

        int leftUpIndex = i - WireWidth - 1;
        if (leftUpIndex >= min && CheckAdjacentLine(leftUpIndex, i))
            Connect(leftUpIndex, nextLayerIndex);

        int rightUpIndex = i - WireWidth + 1;
        if (rightUpIndex >= min && CheckAdjacentLine(rightUpIndex, i))
            Connect(rightUpIndex, nextLayerIndex);

        int leftDownIndex = i + WireWidth - 1;
        if (leftDownIndex < max && CheckAdjacentLine(leftDownIndex, i))
            Connect(leftDownIndex, nextLayerIndex);

        int rightDownIndex = i + WireWidth + 1;
        if (rightDownIndex < max && CheckAdjacentLine(rightDownIndex, i))
            Connect(rightDownIndex, 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 || A[i][j] == -1)
            {
                B[i][bIndex] = j;
                bIndex++;
                C[i]++;
            }
        }
    }
}

int main(int argc, char *argv[])
{
    sscanf(argv[1], "%lf", &sigma_G_1);
    sscanf(argv[2], "%lf", &sigma_N_1);

    sigma_G = sigma_G_1 / 17;
    sigma_GN = sigma_G;
    sigma_N = sigma_N_1 / 17;

    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();
    FillNetworkWireMatrix();
    FillLayerConnectionMatrix();
    //FillRandomMatrix();
    RandomizeWireMatrix();
    RandomizeInhibitoryWireMatrix();

    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 || A[i][j] == -1)
            {
                links_count++;
            }
        }
        fprintf(fp0, "%d\n", (int)links_count);

    }
    fclose(fp0);

    fp0 = fopen("linksInhib.txt", "w+");
    for (int i = NHalf; i < N; i++)
    {
        bool connect = false;

        for (int j = NHalf; j < N; j++)
        {
            if (A[i][j] != 0)
            {
                fprintf(fp0, "%d\n", (int)A[i][j]);
                connect = true;
                break;
            }
        }

        if (!connect)
            fprintf(fp0, "%d\n", (int)0);
    }
    fclose(fp0);


    const double t_start = 0;
    const double t_max = 2000; //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_1, sigma_N_1);
    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);

}