using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace U_invert
{
    partial class Program
    {
        static Random rnd = new Random();

        static int[] GenerateRandomInput(int levels, int len)
        {
            if (levels < 2 || len < 1000)
            {
                Console.WriteLine("Wrong parameters for input process");
                Environment.Exit(0);
            }

            double[] x = new double[len + 10];
            for (int i = 0; i < x.Length; ++i)
            {
                x[i] = (rnd.Next() % 100) / 100.0;
            }

            for (int step = 0; step < 4; ++step)
            {
                for (int i=0; i < x.Length - 1; ++i)
                {
                    x[i] = x[i] + x[i + 1];
                }
            }

            double min = x.Min();
            double max = x.Max();
            double delta = (max - min) / (levels);

            int[] xx = new int[len];
            for (int i = 5; i < len + 5; ++i)
            {
                xx[i - 5] = (int)Math.Round((x[i] - min) / delta) - 1;
                if (xx[i - 5] < 0) xx[i - 5] = 0;
                if (xx[i - 5] > levels - 1) xx[i - 5] = levels - 1;
            }

            return xx;
        }

        static double[,] GenerateOperator(int nx, int nt)
        {
            double[,] U = new double[nx, nt];
            for (int i = 0; i < nt; ++i)
            {
                for (int j = 0; j < nx; ++j)
                {
                    U[j, i] = rnd.Next() % 100;
                }
            }

            double average = 0.0;
            for (int i = 0; i < nt; ++i)
            {
                for (int j = 0; j < nx; ++j)
                {
                    average += U[j, i];
                }
            }
            average /= nt;
            average /= nx;

            for (int i = 0; i < nt; ++i)
            {
                for (int j = 0; j < nx; ++j)
                {
                    U[j, i] -= average;
                }
            }

            for (int counter = 0; counter < 3; ++counter)
            {
                for (int i = 0; i < nt; ++i)
                {
                    for (int j = 0; j < nx - 1; ++j)
                    {
                        U[j, i] = (U[j, i] + U[j + 1, i]) / 2.0;
                    }
                    for (int j = nx - 1; j > 0; --j)
                    {
                        U[j, i] = (U[j, i] + U[j - 1, i]) / 2.0;
                    }
                }
                for (int j = 0; j < nx; ++j)
                {
                    for (int i = 0; i < nt - 1; ++i)
                    {
                        U[j, i] = (U[j, i] + U[j, i + 1]) / 2.0;
                    }
                    for (int i = nt - 1; i > 0; --i)
                    {
                        U[j, i] = (U[j, i] + U[j, i - 1]) / 2.0;
                    }
                }
            }

            double k = Math.Exp(4.0 / nt);
            for (int i = 0; i < nx; ++i)
            {
                double p = 1.0;
                for (int j = 0; j < nt; ++j)
                {
                    p /= k;
                    U[i, j] *= p;
                }
            }

            return U;
        }

        static double[] ComputeOutput(int[] x, double[,] U)
        {
            double[] z = new double[x.Length];
            int T = U.GetLength(1);
 
            for (int i = T - 1; i < z.Length; ++i)
            {
                z[i] = 0.0;
                for (int j = 0; j < T; ++j)
                {
                    z[i] += U[x[i - j], j];
                }
            }

            double av = z.Average();
            for (int i = 0; i < T - 1; ++i) z[i] = av;

            return z;
        }

        static double[,] IdentifyOperator(int[] x, double[] z, int nx, int nt)
        {
            double learning_rate = 0.05;
            double[,] U = new double[nx, nt];
            double zmin = z.Min();
            double zmax = z.Max();

            double lastError = 0.0;
            for (int step = 0; step < 8; ++step)
            {
                double accuracy = 0.0;
                int cnt = 0;
                for (int i = nt - 1; i < z.Length; ++i)
                {
                    double predicted = 0.0;
                    for (int j = 0; j < nt; ++j)
                    {
                        predicted += U[x[i - j], j];
                    }
                    double diff = z[i] - predicted;
                    accuracy += diff * diff;
                    ++cnt;
                    double error = diff / nt * learning_rate;
                    for (int j = 0; j < nt; ++j)
                    {
                        U[x[i - j], j] += error;
                    }
                }
                accuracy /= cnt;
                accuracy = Math.Sqrt(accuracy) / (zmax - zmin);
                lastError = accuracy;                
            }
            Console.WriteLine("Relative error for identification {0:0.0000}", lastError);

            return U;
        }

        static int[] AdjustInput(double[] z, int[] x, double[,] U)
        {
            int[] tmpX = new int[x.Length];
            for (int i = 0; i < tmpX.Length; ++i)
            {
                tmpX[i] = x[i];
            }

            int nx = U.GetLength(0);
            int nt = U.GetLength(1);
            for (int i = nt - 1; i < z.Length - nt; ++i)
            {
                int selectedPos = 0;
                double min = 0.0;
                for (int j = 0; j < nx; ++j)
                {
                    tmpX[i] = j;
                    double abs = GetFragmentAccuracy(tmpX, z, U, i);
                    if (0 == j)
                    {
                        selectedPos = 0;
                        min = abs;
                    }
                    else
                    {
                        if (abs < min)
                        {
                            min = abs;
                            selectedPos = j;
                        }
                    }
                }
                tmpX[i] = selectedPos;
            }

            return tmpX;
        }

        static double GetFragmentAccuracy(int[] x, double[] z, double[,] U, int currentPos)
        {
            int nt = U.GetLength(1);
            double abs = 0.0;
            for (int i = 0; i < nt; ++i)
            {
                double output = 0.0;
                for (int j = 0; j < nt; ++j)
                {
                    output += U[x[currentPos - j + i], j];
                }
                abs += Math.Abs(z[currentPos + i] - output);
            }

            return abs;
        }

        static void ShowOperator(double[,] U)
        {
            int nt = U.GetLength(1);
            int nx = U.GetLength(0);
            Console.WriteLine("-------------- Operator, each row is time layer --------------");
            for (int i = 0; i < nt; ++i)
            {
                for (int j = 0; j < nx - 1; ++j)
                {
                    Console.Write("{0:0.00}, ", U[j, i]);
                }
                Console.Write("{0:0.00}", U[nx - 1, i]);
                Console.WriteLine();
            }
            Console.WriteLine("----------------------");
        }

        static void Main(string[] args)
        {
            int N = 32000;
            int levels = 20;
            int timeLayers = 25;

            //Generate data
            double[,] U = GenerateOperator(levels, timeLayers);
            int[] originalInput = GenerateRandomInput(levels, N);
            double[] z = ComputeOutput(originalInput, U);
            //double[,] UI = IdentifyOperator(originalInput, z, levels, timeLayers);
            //
            //Adjust input
            int[] adjustedInput = GenerateRandomInput(levels, N);
            int steps = 16;
            for (int step = 0; step < steps; ++step)
            {
                adjustedInput = AdjustInput(z, adjustedInput, U);
                Console.WriteLine("adistment loop {0} of {1} lools", step, steps);
            }

            //Verify accuracy
            double[] zvalidate = ComputeOutput(adjustedInput, U);
            double min = z.Min();
            double max = z.Max();
            double diff = 0.0;
            int cnt = 0;
            for (int i = timeLayers - 1; i < N - timeLayers; ++i)
            {
                diff += (z[i] - zvalidate[i]) * (z[i] - zvalidate[i]);
                ++cnt;
            }
            diff /= cnt;
            diff = Math.Sqrt(diff);
            diff /= (max - min);
            Console.WriteLine("Relative difference for outputs {0:0.0000}", diff);
        }
    }
}