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

namespace Operator
{
    class Program
    {
        private const int N = 10000;
        private const double decline = 0.8;
        private const int smoothLevel = 2;
        private const int _xn = 25;
        private const int _tn = 21;

        private static int[,] map = new int[_xn, _tn];
        private static double[,] U = new double[_xn, _tn];
        private static double[,] U1 = new double[_xn, _tn];
        private static double[,] U2 = new double[_xn, _tn];
        private static double[,] U3 = new double[_xn, _tn];
        private static double[,] U4 = new double[_xn, _tn];
        private static int[] x1 = new int[N];
        private static double[] y1 = new double[N];
        private static int[] x2 = new int[N];
        private static double[] y2 = new double[N];

        private static void MakeOperator()
        {
            Random rnd = new Random();
            for (int i = 0; i < _xn; ++i)
            {
                for (int j = 0; j < _tn; ++j)
                {
                    U[i, j] = (rnd.Next() % 100) / 100.0;
                }
            }
            for (int counter = 0; counter < smoothLevel; ++counter)
            {
                for (int i = 0; i < _xn; ++i)
                {
                    for (int j = 1; j < _tn; ++j)
                    {
                        U[i, j] = (U[i, j] + U[i, j - 1]) / 2.0;
                    }
                }
                for (int j = 0; j < _tn; ++j)
                {
                    for (int i = 1; i < _xn; ++i)
                    {
                        U[i, j] = (U[i, j] + U[i - 1, j]) / 2.0;
                    }
                }
                for (int i = _xn - 1; i >= 0; --i)
                {
                    double rate = 1.0;
                    for (int j = 0; j < _tn; ++j)
                    {
                        if (j < _tn / 3) rate /= decline;
                        else rate *= decline;
                        U[i, j] *= rate;
                    }
                }
                for (int j = _tn - 1; j >= 0; --j)
                {
                    double rate = 1.0;
                    for (int i = _xn - 1; i >= 0; --i)
                    {
                        if (i < _xn / 3) rate /= decline;
                        else rate *= decline;
                        U[i, j] *= rate;
                    }
                }
                for (int i = 0; i < _xn; ++i)
                {
                    for (int j = 1; j < _tn; ++j)
                    {
                        U[i, j] = (U[i, j] + U[i, j - 1]) / 2.0;
                    }
                }
                for (int j = 0; j < _tn; ++j)
                {
                    for (int i = 1; i < _xn; ++i)
                    {
                        U[i, j] = (U[i, j] + U[i - 1, j]) / 2.0;
                    }
                }
            }
        }

        private static void ShowOperator<T>(string label, T[,] U)
        {
            Console.WriteLine(label);
            string format = "{0:0} ";
            if (typeof(T).Equals(typeof(double)))
            {
                format = "{0:0.000} ";
            }
            for (int i = 0; i < _xn; ++i)
            {
                for (int j = 0; j < _tn; ++j)
                {
                    Console.Write(format, U[i, j]);
                }
                Console.WriteLine();
            }
            Console.WriteLine();
        }

        private static void MakeX(int[] x)
        {
            Random rnd = new Random();
            for (int i = 0; i < N; ++i)
            {
                x[i] = rnd.Next() % _xn;
            }
            for (int i = 1; i < N; ++i)
            {
                x[i] = (x[i] + x[i - 1]) / 2;
            }
        }

        private static void ComputeY(int[] x, double[,] U, double[] y)
        {
            for (int i = 0; i < _tn; ++i)
            {
                y[i] = 0.0;
            }
            for (int i = _tn - 1; i < N; ++i)
            {
                y[i] = 0.0;
                for (int j = 0; j < _tn; ++j)
                {
                    y[i] += U[x[i - j], j];
                }
            }
        }

        static double AddError(int[] x)
        {
            Random rnd = new Random();
            int total = 0;
            int totalX = 0;
            for (int i = 0; i < x.Length; ++i)
            {
                int xold = x[i];
                int level = rnd.Next() % 100;
                if (level < 20) x[i] -= 2;
                else if (level < 40) x[i] -= 1;
                else if (level < 60) x[i] = x[i];
                else if (level < 80) x[i] += 1;
                else x[i] += 2;
                if (x[i] >= _xn) x[i] = _xn - 1;
                if (x[i] <= 0) x[i] = 0;
                int xnew = x[i];
                total += Math.Abs(xold - xnew);
                totalX += Math.Abs(xold);
            }
            return (double)(total) / (double)(totalX);
        }

        static void IdentificationClassic(int[] x, double[] y, double[,] U)
        {
            for (int i = 0; i < _xn; ++i)
            {
                for (int j = 0; j < _tn; ++j)
                {
                    U[i, j] = 0.0;
                }
            }
            for (int counter = 0; counter < 10; ++counter)
            {
                for (int i = _tn - 1; i < N; ++i)
                {
                    double delta = 0.0;
                    for (int j = 0; j < _tn; ++j)
                    {
                        delta += U[x[i - j], j];
                    }
                    delta = (y[i] - delta) / _tn;
                    for (int j = 0; j < _tn; ++j)
                    {
                        U[x[i - j], j] += delta;
                    }
                }
            }
        }

        // Next block is identification with linear interpolation of the kernel
        private static void MakeMap(int steps)
        {
            for (int i = 0; i < _xn; ++i)
            {
                for (int j = 0; j < _tn; ++j)
                {
                    map[i, j] = 0;
                }
            }
            int step = _xn / steps;
            if (step < 1) step = 1;
            for (int j = 0; j < _tn; ++j)
            {
                for (int i = 0; i < _xn; ++i)
                {
                    if (i % step == 0) map[i, j] = 1;
                }
                if (map[_xn - 1, j] == 0) map[_xn - 1, j] = 1;
            }
        }

        private static int GetTop(int x, int j)
        {
            if (map[x, j] == 1) return x;
            while (true)
            {
                if (map[++x, j] == 1) return x;
            }
        }

        private static int GetBottom(int x, int j)
        {
            if (map[x, j] == 1) return x;
            while (true)
            {
                if (map[--x, j] == 1) return x;
            }
        }

        private static double GetLinear(int pos1, int pos2, double U1, double U2, int pos)
        {
            if (pos == pos1) return U1;
            if (pos == pos2) return U2;
            return (U2 - U1) / (double)(pos2 - pos1) * (double)(pos - pos1) + U1;
        }

        private static void IdentificationLinear(int[] x, double[] y, double[,] U)
        {
            MakeMap(5);
            //set all zeros
            for (int i = 0; i < _xn; ++i)
            {
                for (int j = 0; j < _tn; ++j)
                {
                    U[i, j] = 0.0;
                }
            }
            for (int counter = 0; counter < 50; ++counter)
            {
                double avdelta = 0.0;
                double avoutput = 0.0;
                for (int i = _tn - 1; i < N; ++i)
                {
                    //get difference
                    double delta = 0.0;
                    for (int j = 0; j < _tn; ++j)
                    {
                        int pos = (int)(x[i - j]);
                        int top = GetTop(x[i - j], j);
                        int bottom = GetBottom(x[i - j], j);
                        double u = GetLinear(top, bottom, U[top, j], U[bottom, j], pos);
                        delta += u;
                    }
                    delta = (y[i] - delta) / _tn;
                    avdelta += delta * delta;
                    avoutput += y[i] * y[i];
                    //update
                    for (int j = 0; j < _tn; ++j)
                    {
                        int pos = (int)(x[i - j]);
                        int top = GetTop(x[i - j], j);
                        int bottom = GetBottom(x[i - j], j);
                        double deltaTop = 0.0;
                        double deltaBottom = 0.0;
                        if (pos == top) deltaTop = delta;
                        else if (pos == bottom) deltaBottom = delta;
                        else
                        {
                            deltaTop = delta / (double)(top - bottom) * (double)(pos - bottom);
                            deltaBottom = delta / (double)(top - bottom) * (double)(top - pos);
                        }
                        U[top, j] += deltaTop;
                        U[bottom, j] += deltaBottom;
                    }
                }
                avdelta /= (double)(N - _tn);
                avdelta = Math.Sqrt(avdelta);
                avoutput /= (double)(N - _tn);
                avoutput = Math.Sqrt(avoutput);
                //Console.WriteLine("Precision {0}", avdelta / avoutput);
            }
            //Filling operator gap
            for (int j = 0; j < _tn; ++j)
            {
                int pos1 = 0;
                for (int i = 1; i < _xn; ++i)
                {
                    if (map[i, j] > 0)
                    {
                        int pos2 = i;
                        if (pos2 - pos1 > 1)
                        {
                            double u1 = U[pos1, j];
                            double u2 = U[pos2, j];
                            for (int k = pos1 + 1; k < pos2; ++k)
                            {
                                U[k, j] = (u2 - u1) / (double)(pos2 - pos1) * (double)(k - pos1) + u1;
                            }
                        }
                        pos1 = pos2;
                    }
                }
            }
        }
        //end identification with linearization of kernel

        static double ComputeAccuracy(int[] x, double[] y, double[,] U)
        {
            double totalErr = 0.0;
            double totalY = 0.0;
            for (int i = _tn - 1; i < N; ++i)
            {
                double result = 0.0;
                for (int j = 0; j < _tn; ++j)
                {
                    result += U[x[i - j], j];
                }
                double err = y[i] - result;
                err *= err;
                totalY += y[i] * y[i];
                totalErr += err;
            }
            double accuracy = totalErr / totalY / (double)(N - _tn);
            return accuracy;
        }

        static void Main(string[] args)
        {
            Console.WriteLine("Test started");
            MakeOperator();
            MakeX(x1);
            Thread.Sleep(1000);
            MakeX(x2);
            ComputeY(x1, U, y1);
            ComputeY(x2, U, y2);
            Console.WriteLine("Relative error is added to x1 {0}", AddError(x1));
            Thread.Sleep(100);
            Console.WriteLine("Relative error is added to x2 {0}", AddError(x2));
            IdentificationClassic(x1, y1, U1);
            IdentificationClassic(x2, y2, U2);
            IdentificationLinear(x1, y1, U3);
            IdentificationLinear(x2, y2, U4);

            ShowOperator("Original operator", U);
            ShowOperator("Classic identified operator U1", U1);
            ShowOperator("Classic identified operator U2", U2);
            ShowOperator("Linear", U3);
            ShowOperator("Linear", U4);
            ShowOperator("Map", map);

            Console.WriteLine("Test completed, accuracies are below");
            Console.WriteLine("Accuracy of classic identification, conducted on sample 1, tested on sample2 {0}", ComputeAccuracy(x2, y2, U1));
            Console.WriteLine("Accuracy of classic identification, conducted on sample 2, tested on sample1 {0}", ComputeAccuracy(x1, y1, U2));
            Console.WriteLine("Accuracy of linear identification, conducted on sample 1, tested on sample2 {0}", ComputeAccuracy(x2, y2, U3));
            Console.WriteLine("Accuracy of linear identification, conducted on sample 2, tested on sample1 {0}", ComputeAccuracy(x1, y1, U4));
        }
    }
}