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