//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: // // // // // // // // 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); } } }