//It is not necessary to have CS 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 commands:

//<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>

//File name "Program.cs" is this code

//Using DOS prompt navigate to folder with MSBuild.exe and execute command:
//MSBuild C:\YourFolder\YourProjectFile.csproj

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

namespace UCon
{
    class Data
    {
        const int nx = 15;
        const int ny = 10;
        const int nt = 17;
        const int xysize = 5000;

        public double[, ,] U = new double[nx, ny, nt];
        public double[, ,] U1 = new double[nx, ny, nt];
        public double[, ,] U2 = new double[nx, ny, nt];

        public int[] x1 = new int[xysize];
        public int[] y1 = new int[xysize];
        public double[] z1 = new double[xysize];

        public int[] x2 = new int[xysize];
        public int[] y2 = new int[xysize];
        public double[] z2 = new double[xysize];

        public double[] ztmp = new double[xysize];

        public int getnx()
        {
            return nx;
        }

        public int getny()
        {
            return ny;
        }

        public int getnt()
        {
            return nt;
        }

        public int getsize()
        {
            return xysize;
        }
    }

    class Model
    {
        private Data data = null;

        private void InitializeInput(int[] X, int limit)
        {
            Thread.Sleep(20);
            Random rnd = new Random();
            for (int i = 0; i < X.Length; ++i)
            {
                X[i] = rnd.Next() % limit;
            }
        }

        private void ComputeOutput(int[] X, int[] Y, double[] Z, double[, ,] G, int nT)
        {
            for (int i = 0; i < nT + 1; ++i)
            {
                Z[i] = 0.0;
            }

            for (int i = nT - 1; i < Z.Length; ++i)
            {
                Z[i] = 0.0;
                for (int j = 0; j < nT; ++j)
                {
                    Z[i] += G[X[i - j], Y[i - j], j];
                }
            }
        }

        private double ComputeError(double[] X, double[] Y)
        {
            int size = X.Length;
            if (size < Y.Length) size = Y.Length;

            double error = 0.0;
            for (int i = 0; i < size; ++i)
            {
                double diff = X[i] - Y[i];
                if (diff >= 0) error += diff;
                else error -= diff;
            }

            if (size > 0) error /= (double)(size);
            return error;
        }

        private void InitializeOperator(double[, ,] G, int nX, int nY, int nT)
        {
            Thread.Sleep(20);
            Random rnd = new Random();
            for (int i = 0; i < nX; ++i)
            {
                for (int j = 0; j < nY; ++j)
                {
                    for (int k = 0; k < nT; ++k)
                    {
                        G[i, j, k] = (double)(rnd.Next() % 100) / 10.0;
                    }
                }
            }
        }

        public void InitializeData()
        {
            if (null != data) return;

            data = new Data();

            //inputs for training sample
            InitializeInput(data.x1, data.getnx());
            InitializeInput(data.y1, data.getny());

            //inputs for control sample
            InitializeInput(data.x2, data.getnx());
            InitializeInput(data.y2, data.getny());

            //exact operator
            InitializeOperator(data.U, data.getnx(), data.getny(), data.getnt());

            //outputs for training and control samples and exact operator
            ComputeOutput(data.x1, data.y1, data.z1, data.U, data.getnt());
            ComputeOutput(data.x2, data.y2, data.z2, data.U, data.getnt());
        }

        public void Identification()
        {
            //assign all zeros as initial values for estimated operators
            for (int i = 0; i < data.getnx(); ++i)
            {
                for (int j = 0; j < data.getny(); ++j)
                {
                    for (int k = 0; k < data.getnt(); ++k)
                    {
                        data.U1[i, j, k] = 0.0;
                    }
                }
            }

            for (int i = 0; i < data.getnx(); ++i)
            {
                for (int j = 0; j < data.getny(); ++j)
                {
                    for (int k = 0; k < data.getnt(); ++k)
                    {
                        data.U2[i, j, k] = 0.0;
                    }
                }
            }

            Console.WriteLine("Identification for first implementation");
            for (int cnt = 0; cnt < 100; ++cnt)
            {
                //identification step
                for (int j = data.getnt(); j < data.getsize(); ++j)
                {
                    double value = 0.0;
                    //compute single output value using training data
                    for (int i = 0; i < data.getnt(); ++i)
                    {
                        value += data.U1[data.x1[j - i], data.y1[j - i], i];
                    }
                    double delta = data.z1[j] - value;
                    delta /= data.getnt();
                    //update estimated operator using trainging data
                    for (int i = 0; i < data.getnt(); ++i)
                    {
                        data.U1[data.x1[j - i], data.y1[j - i], i] += delta;
                    }
                }
                //estimate error on control data
                //we estimate operator on sample x1,y1,z1 and compute discrepancy for x2,y2,z2
                if (cnt % 10 == 0)
                {
                    ComputeOutput(data.x2, data.y2, data.ztmp, data.U1, data.getnt());
                    Console.WriteLine("Residual error = " + ComputeError(data.ztmp, data.z2));
                }
            }
            Console.WriteLine();

            Console.WriteLine("Identification for second implementation");
            for (int cnt = 0; cnt < 100; ++cnt)
            {
                //identification step
                for (int j = data.getnt(); j < data.getsize(); ++j)
                {
                    double value = 0.0;
                    //compute single output value using training data
                    for (int i = 0; i < data.getnt(); ++i)
                    {
                        value += data.U2[data.x2[j - i], data.y2[j - i], i];
                    }
                    double delta = data.z2[j] - value;
                    delta /= data.getnt();
                    //update estimated operator using trainging data
                    for (int i = 0; i < data.getnt(); ++i)
                    {
                        data.U2[data.x2[j - i], data.y2[j - i], i] += delta;
                    }
                }
                //estimate error on control data
                //we estimate operator on sample x2,y2,z2 and compute discrepancy for x1,y1,z1
                if (cnt % 10 == 0)
                {
                    ComputeOutput(data.x1, data.y1, data.ztmp, data.U2, data.getnt());
                    Console.WriteLine("Residual error = " + ComputeError(data.ztmp, data.z1));
                }
            }
            Console.WriteLine();

            //compare two operators
            double operatorDiff = 0.0;
            int counter = 0;
            for (int i = 0; i < data.getnx(); ++i)
            {
                for (int j = 0; j < data.getny(); ++j)
                {
                    for (int k = 0; k < data.getnt(); ++k)
                    {
                        double diff = data.U1[i, j, k] - data.U2[i, j, k];
                        if (diff >= 0) operatorDiff += diff;
                        else operatorDiff -= diff;
                        ++counter;
                    }
                }
            }
            Console.WriteLine("Operators' discrepancy " + operatorDiff / (double)counter);
        }
    }

    class Program
    {
        static void Main(string[] args)
        {
            Model model = new Model();
            model.InitializeData();
            model.Identification();
        }
    }
}