//This program is written by Andrew Polar and Mike Poluektov.
//The goal is to demonstrate identification of Wiener type dynamic system with with two inputs.
//The identified model is two input Urysohn operator that is capable to model Wiener object when nonlinearity is weak.
//The validation of identified operator is tested on unseen data,
//not on data used for building the models.

//Patent pending "Method for identifying discrete urysohn models", 15/998381, filed 08/11/2018.

//C# compiler is free download now, but
//it is not necessary to have C# compiler to build this self-test if you run WIN10.
//Find MSBuild.exe in C:\Windows\Microsoft.NET\Framework64\v4.0.30319 or other version of your framework
//Make ASCII file with extension *.csproj and copy/paste following script:

//<Project xmlns="http://schemas.microsoft.com/developer/msbuild/2003" DefaultTargets="Build" ToolsVersion="12">
//    <ItemGroup>
//        <FilesToCompile Include="Program.cs" />
//    </ItemGroup>
//    <Target Name="Build">
//        <Csc Sources="@(FilesToCompile)" OutputAssembly="$(OutputAssembly)" Optimize="$(Optimize)" TargetType="exe" Resources="@(CompiledResources)"/>
//    </Target>
// </Project>


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

namespace Wiener
{
    static class Generic
    {
        static Random _rnd = new Random();
        static public double[] GetImpulseFunction(int size)
        {
            double[] h = new double[size];
            for (int i = 0; i < size; ++i)
            {
                h[i] = _rnd.Next() % 100.0;
            }

            for (int k = 0; k < 2; ++k)
            {
                for (int i = 0; i < size - 1; ++i)
                {
                    h[i] = (h[i] + h[i + 1]) / 2.0;
                }
                for (int i = size - 1; i > 0.0; --i)
                {
                    h[i] = (h[i] + h[i - 1]) / 2.0;
                }
            }

            double av = h.Average();
            for (int i = 0; i < size; ++i)
            {
                h[i] -= av;
            }

            double coeff = 1.0;
            double decay = 0.85;
            for (int i = 0; i < size; ++i)
            {
                h[i] *= coeff;
                coeff *= decay;
            }

            double norm = 0.0;
            for (int i = 0; i < size; ++i)
            {
                norm += h[i] * h[i];
            }
            norm = Math.Sqrt(norm);
            for (int i = 0; i < size; ++i)
            {
                h[i] /= norm;
            }
            return h;
        }

        static public double[] GenerateSignal(int N)
        {
            double[] x = new double[N];
            for (int i = 0; i < N; ++i)
            {
                x[i] = (_rnd.Next() % 1000 - 500) / 11.0;
            }
            for (int k = 0; k < 4; ++k)
            {
                for (int i = 0; i < N - 1; ++i)
                {
                    x[i] = (x[i] + x[i + 1]) / 2.0;
                }
            }

            double min = x.Min();
            double max = x.Max();
            for (int i = 0; i < N; ++i)
            {
                x[i] = (x[i] - min) / (max - min) * 2.0 - 1.0;
            }

            return x;
        }

        public static double[] ComputerLinear(double[] signal, double[] h)
        {
            int N = signal.Length;
            int size = h.Length;
            double[] reaction = new double[N];
            for (int i = size - 1; i < N; ++i)
            {
                reaction[i] = 0.0;
                for (int j = 0; j < size; ++j)
                {
                    reaction[i] += signal[i - j] * h[j];
                }
            }

            double av = reaction.Average();
            for (int i = 0; i < size - 1; ++i)
            {
                reaction[i] = av;
            }

            return reaction;
        }
    }

    class Program
    {
        private static double GetPoint(double x, double y, double NW, double NE, double SW, double SE)
        {
            return NW + (NE - NW) * x + (SW - NW) * y + ((SE + NW) - (NE + SW)) * x * y;
        }

        private static double GetStaticNonlinearValue(double y)
        {
            return Math.Sign(y) * y * y + 0.2 * y * y * y;
        }

        private static void IdentificationTISO(double[] x, double[] y, double[] z, int X_size, int Y_size, int T_size)
        {
            int N = x.Length;
            int M = y.Length;
            if (M < N) N = M;

            double mu = 0.1;
            double xmin = x.Min();
            double xmax = x.Max();
            double ymin = y.Min();
            double ymax = y.Max();
            double zmin = z.Min();
            double zmax = z.Max();
            double deltaX = (xmax - xmin) * 0.999 / (double)(X_size - 2);
            double deltaY = (ymax - ymin) * 0.999 / (double)(Y_size - 2);

            double[,,] Model = new double[X_size, Y_size, T_size];

            for (int counter = 0; counter < 25; ++counter)
            {
                double accuracy = 0.0;
                int cnt = 0;
                for (int i = T_size - 1; i < N / 2; ++i)
                {
                    double signal = 0.0;
                    for (int j = 0; j < T_size; ++j)
                    {
                        int X = (int)((x[i - j] - xmin) / deltaX);
                        int Y = (int)((y[i - j] - ymin) / deltaY);

                        double xr = ((X + 1) * deltaX - (x[i - j] - xmin)) / deltaX;
                        double yr = ((Y + 1) * deltaY - (y[i - j] - ymin)) / deltaY;

                        signal += GetPoint(xr, yr, Model[X + 1, Y, j], Model[X + 1, Y + 1, j], Model[X, Y, j], Model[X, Y + 1, j]);
                    }
                    double error = z[i] - signal;

                    accuracy += error * error;
                    ++cnt;

                    error /= T_size;
                    error *= mu;

                    for (int j = 0; j < T_size; ++j)
                    {
                        int X = (int)((x[i - j] - xmin) / deltaX);
                        int Y = (int)((y[i - j] - ymin) / deltaY);

                        double xr = ((X + 1) * deltaX - (x[i - j] - xmin)) / deltaX;
                        double yr = ((Y + 1) * deltaY - (y[i - j] - ymin)) / deltaY;

                        double wNW = (1.0 - xr) * (1.0 - yr);
                        double wNE = xr * (1.0 - yr);
                        double wSW = (1.0 - xr) * yr;
                        double wSE = xr * yr;

                        Model[X + 1, Y, j] += error * wNW;
                        Model[X + 1, Y + 1, j] += error * wNE;
                        Model[X, Y, j] += error * wSW;
                        Model[X, Y + 1, j] += error * wSE;
                    }
                }
                accuracy = Math.Sqrt(accuracy / cnt) / (zmax - zmin);
                Console.WriteLine("Relative error on self data {0:0.0000}", accuracy);
            }

            //Validation
            double accuracy2 = 0.0;
            int cnt2 = 0;
            for (int i = N / 2 + T_size - 1; i < N; ++i)
            {
                double signal = 0.0;
                for (int j = 0; j < T_size; ++j)
                {
                    int X = (int)((x[i - j] - xmin) / deltaX);
                    int Y = (int)((y[i - j] - ymin) / deltaY);

                    double xr = ((X + 1) * deltaX - (x[i - j] - xmin)) / deltaX;
                    double yr = ((Y + 1) * deltaY - (y[i - j] - ymin)) / deltaY;

                    signal += GetPoint(xr, yr, Model[X + 1, Y, j], Model[X + 1, Y + 1, j], Model[X, Y, j], Model[X, Y + 1, j]);
                }
                double error = z[i] - signal;
                //Console.WriteLine(z[i] + " " + signal);

                accuracy2 += error * error;
                ++cnt2;
            }
            accuracy2 = Math.Sqrt(accuracy2 / cnt2) / (zmax - zmin);
            Console.WriteLine("------------- Validation --------------");
            Console.WriteLine("Relative error on unseen data {0:0.0000}", accuracy2);
        }

        private static void ShowStaticNonlinearity(double min, double max)
        {
            int steps = 10;
            double delta = (max - min) / (steps - 1);
            double x = min;
            for (int i = 0; i < steps; ++i)
            {
                Console.WriteLine("{0:0.0000}, {1:0.0000}", x, GetStaticNonlinearValue(x));
                x += delta;
            }
        }

        static void Main(string[] args)
        {
            //Generate data
            int N = 100000;
            int t = 20;
            double[] h1 = Generic.GetImpulseFunction(t);
            double[] h2 = Generic.GetImpulseFunction(t);
            double[] x1 = Generic.GenerateSignal(N);
            double[] x2 = Generic.GenerateSignal(N);
            double[] y1 = Generic.ComputerLinear(x1, h1);
            double[] y2 = Generic.ComputerLinear(x2, h2);
            double[] y = new double[N];
            for (int i = 0; i < N; ++i)
            {
                y[i] = y1[i] + y2[i];
                y[i] = GetStaticNonlinearValue(y[i]);
            }
            double ymin = y.Min();
            double ymax = y.Max();
            Console.WriteLine("Data is ready output min = {0:0.0000}, max = {1:0.0000}", ymin, ymax);

            //Identification 2-dim Urysohn
            IdentificationTISO(x1, x2, y, 25, 30, t);
        }
    }
}