using System; using System.Collections.Generic; using System.Linq; using System.Text; using System.Threading.Tasks; using System.Threading; namespace TwoU { class Cholesky { //main method for Cholesky decomposition. //input n size of matrix //input/output a Symmetric positive def. matrix //output p vector of resulting diag of a //author: //Converted to c# by Andrew Polar private bool choldc1(int n, double[,] a, double[] p) { int i, j, k; double sum; for (i = 0; i < n; i++) { for (j = i; j < n; j++) { sum = a[i, j]; for (k = i - 1; k >= 0; k--) { sum -= a[i, k] * a[j, k]; } if (i == j) { if (sum <= 0) { Console.WriteLine("Matrix is not positive definite!"); return false; } p[i] = Math.Sqrt(sum); } else { a[j, i] = sum / p[i]; } } } return true; } private bool choldc(int n, double[,] A, double[,] a) { int i, j; double[] p = new double[n]; for (i = 0; i < n; i++) for (j = 0; j < n; j++) a[i, j] = A[i, j]; if (!choldc1(n, a, p)) { return false; } for (i = 0; i < n; i++) { a[i, i] = p[i]; for (j = i + 1; j < n; j++) { a[i, j] = 0; } } p = null; return true; } /////////////////////////////////////////////////////////////////// //solution for M * x = y //M must be positive definite public bool CholeskySolution(double[,] M, double[] x, double[] y, int size) { //regularization //double trace = 0.0; //for (int i = 0; i < size; ++i) //{ // trace += M[i, i]; //} //trace /= size; //trace *= 0.001; //for (int i = 0; i < size; ++i) //{ // M[i, i] += trace; //} //end regularization double[,] T = new double[size, size]; //find lower triangular if (!choldc(size, M, T)) { return false; } //make is symmetric for (int i = 0; i < size; ++i) for (int j = i + 1; j < size; ++j) T[i, j] = T[j, i]; double[] v = new double[size]; //solve lower triangular for (int i = 0; i < size; ++i) { double sum = 0.0; for (int j = 0; j < i; ++j) { sum += T[i, j] * v[j]; } v[i] = (y[i] - sum) / T[i, i]; } //solve upper triangular for (int i = size - 1; i >= 0; --i) { double sum = 0.0; for (int j = size - 1; j > i; --j) { sum += T[i, j] * x[j]; } x[i] = (v[i] - sum) / T[i, i]; } v = null; T = null; return true; } public void SelfTest() { const int size = 10; double[,] H = new double[size, size]; for (int i = 0; i < size; ++i) { for (int j = 0; j < size; ++j) { H[i, j] = 1.0 / (double)(i + j + 1); } } double[] x = new double[size]; double[] y = new double[size]; for (int i = 0; i < size; ++i) { x[i] = (double)(i + 1); y[i] = 0.0; } for (int i = 0; i < size; ++i) { for (int j = 0; j < size; ++j) { y[i] += H[i, j] * x[j]; } } for (int i = 0; i < size; ++i) { x[i] = 0.0; } Cholesky ch = new Cholesky(); if (ch.CholeskySolution(H, x, y, size)) { for (int i = 0; i < size; ++i) { Console.WriteLine(x[i]); } } else { Console.WriteLine("fatal: failed to obtain cholesky solution"); Environment.Exit(1); } } } class Program { static int _nt = 30; static int _nx = 40; static int _N = 100000; static int _NTraining = _N / 2; static int _NLinear = _nt / 2; static double noise_suppress = 0.01; static double[,] _U = new double[_nx, _nt]; static double[,] _Ui = new double[_nx, _nt]; static int[] _x = new int[_N]; static double[] _y = new double[_N]; static void GenerateX() { Thread.Sleep(200); Random rnd = new Random(); for (int i = 0; i <_N; ++i) { _x[i] = rnd.Next() % 100; } for (int k = 0; k < 4; ++k) { for (int i = 1; i < _N; ++i) { _x[i] = (_x[i] + _x[i - 1]) / 2; } } int min = _x.Min(); int max = _x.Max(); for (int i=0; i < _N; ++i) { _x[i] = (int)((_x[i] - min) / (double)(max - min) * (double)(_nx)); if (_x[i] < 0) _x[i] = 0; if (_x[i] > _nx - 1) _x[i] = _nx - 1; } } static void ShowX() { for (int i = 0; i < _N; ++i) { Console.WriteLine("{0} ", _x[i]); } Console.WriteLine(); } static void ShowY() { for (int i = 0; i < _N; ++i) { Console.Write("{0} ", _y[i]); } Console.WriteLine(); } static void ShowOperators() { for (int i = 0; i <_nx; ++i) { for (int j = 1; j < _nt - 1; ++j) { Console.Write("{0:0.00}, ", _U[i, j]); } Console.Write("{0:0.00}", _U[i, _nt - 1]); Console.WriteLine(); } Console.WriteLine("----------------------------"); for (int i = 0; i < _nx; ++i) { for (int j = 1; j < _nt - 1; ++j) { Console.Write("{0:0.00}, ", _Ui[i, j]); } Console.Write("{0:0.00}", _Ui[i, _nt - 1]); Console.WriteLine(); } } static void GenerateOperator() { for (int i = 0; i < _nt; ++i) { for (int j = 0; j < _nx; ++j) { _U[j, i] = 0.0; } } Thread.Sleep(200); Random rnd = new Random(); for (int k = 0; k < 10; ++k) { for (int i = k; i < _nt - k; ++i) { for (int j = k; j < _nx - k; ++j) { _U[j, i] += rnd.Next() % 10; } } } for (int counter = 0; counter < 2; ++counter) { for (int i = 0; i < _nt; ++i) { for (int j = 1; j < _nx; ++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 = 1; i < _nt; ++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; } } } for (int i = 0; i < _nx; ++i) { double k = 1.0; for (int j = 0; j < _nt; ++j) { k /= 1.2; _U[i, j] *= k; } } } static void ComputeY() { for (int i = 0; i < _nt; ++i) _y[i] = 0; for (int i = _nt - 1; i <_N; ++i) { double y = 0.0; for (int j = 0; j <_nt; ++j) { y += _U[_x[i - j], j]; } _y[i] = y; } } static double ComputePearson(double[] a, double[] b) { double ma = 0.0; double mb = 0.0; double mamb = 0.0; double ma2 = 0.0; double mb2 = 0.0; int len = a.Length; for (int i = 0; i < len; ++i) { ma += a[i]; mb += b[i]; mamb += a[i] * b[i]; ma2 += a[i] * a[i]; mb2 += b[i] * b[i]; } ma /= len; mb /= len; mamb /= len; ma2 /= len; mb2 /= len; double pearson = mamb - ma * mb; pearson /= Math.Sqrt(ma2 - ma * ma); pearson /= Math.Sqrt(mb2 - mb * mb); return pearson; } static void IdentifyTraining() { int nSteps = 100; for (int step = 0; step < nSteps; ++step) { double accuracy = 0.0; int cnt = 0; double min = _y[_nt]; double max = _y[_nt]; for (int i = _nt; i < _NTraining; ++i) { double delta = 0.0; for (int j = 0; j < _nt; ++j) { delta += _Ui[_x[i - j], j]; } delta = _y[i] - delta; accuracy += Math.Abs(delta); ++cnt; delta /= _nt; for (int j = 0; j < _nt; ++j) { _Ui[_x[i - j], j] += delta * noise_suppress; } if (_y[i] < min) min = _y[i]; if (_y[i] > max) max = _y[i]; } if (step >= nSteps - 1) { Console.WriteLine("Accuracy for training data = {0:0.0000}", accuracy / (max - min) / (double)(cnt)); } } } static void Validate() { double accuracy = 0.0; double min = _y[_NTraining]; double max = _y[_NTraining]; int cnt = 0; for (int i = _NTraining; i < _N; ++i) { double delta = 0.0; for (int j = 0; j < _nt; ++j) { delta += _Ui[_x[i - j], j]; } delta = _y[i] - delta; accuracy += Math.Abs(delta); ++cnt; if (_y[i] < min) min = _y[i]; if (_y[i] > max) max = _y[i]; } Console.WriteLine("Accuracy for unseen data without an error = {0}", accuracy / (max - min) / (double)(cnt)); } static void AddCorrelatedError(int p) { Thread.Sleep(200); Random rnd = new Random(); double minx = _x.Min(); double maxx = _x.Max(); double[] err = new double[_NTraining]; err[0] = 0.0; err[1] = 0.0; err[2] = 0.0; for (int i = 3; i < _NTraining; ++i) { err[i] = 0.2 * _x[i] + 0.2 * _x[i - 1] + 0.8 * _x[i - 2] + rnd.Next() % maxx * p / 100.0; } double min = err.Min(); double max = err.Max(); double av = err.Average(); for (int i = 0; i < _NTraining; ++i) { err[i] = (err[i] - av) / (max - min) * (maxx - minx) * p / 100.0; } double[] xd = new double[_NTraining]; for (int i = 0; i < xd.Length; ++i) { xd[i] = _x[i]; } double pearson = ComputePearson(err, xd); Console.WriteLine("Pearson corr. coef. for input and error: {0:0.0000}", pearson); for (int i = 0; i < _NTraining; ++i) { _x[i] += (int)err[i]; if (_x[i] < 0) _x[i] = 0; if (_x[i] >= _nx) _x[i] = _nx - 1; } min = err.Min(); max = err.Max(); minx = _x.Min(); maxx = _x.Max(); Console.WriteLine("Correlated error is added:"); Console.WriteLine("Ranges for x[{0}, {1}] and for error [{2}, {3}]", minx, maxx, min, max); } static void TestDifferentSize() { int R1 = 17; int R2 = 10; int nx = _nx * R1 / R2; int[] x = new int[_N]; for (int i = 0; i < _N; ++i) { x[i] = _x[i] * R1 / R2; if (x[i] > nx - 1) x[i] = nx - 1; } double[,] _U2 = new double[nx, _nt]; int nSteps = 100; for (int step = 0; step < nSteps; ++step) { double accuracy = 0.0; int cnt = 0; double min = _y[_nt]; double max = _y[_nt]; for (int i = _nt; i <_NTraining; ++i) { double delta = 0.0; for (int j = 0; j < _nt; ++j) { delta += _U2[x[i - j], j]; } delta = _y[i] - delta; accuracy += Math.Abs(delta); ++cnt; delta /= _nt; for (int j = 0; j < _nt; ++j) { _U2[x[i - j], j] += delta * noise_suppress; } if (_y[i] < min) min = _y[i]; if (_y[i] > max) max = _y[i]; } if (step >= nSteps - 1) { Console.WriteLine("Accuracy for resized operator (training data) = {0}", accuracy / (max - min) / (double)(cnt)); } } } static void TestLinearModelInput() { int NLinear = 2 * _NLinear; int size = _NTraining - NLinear; double[] vector = new double[size]; for (int i = _nt - 1; i < size; ++i) { vector[i] = _y[i + NLinear]; } double[,] matrix = new double[size, NLinear]; for (int k = 0; k < NLinear; ++k) { for (int i = _nt - 1; i < size; ++i) { matrix[i, k] = _x[i + k]; } } int nCols = NLinear; double[] MTV = new double[nCols]; for (int i = 0; i < nCols; ++i) { MTV[i] = 0.0; for (int j = _nt - 1; j < size; ++j) { MTV[i] += vector[j] * matrix[j, i]; } } double[,] MTM = new double[nCols, nCols]; for (int i = 0; i < nCols; ++i) { for (int j = 0; j < nCols; ++j) { MTM[i, j] = 0.0; for (int k = _nt - 1; k < size; ++k) { MTM[i, j] += matrix[k, i] * matrix[k, j]; } } } double[] z = new double[nCols]; Cholesky ch = new Cholesky(); ch.CholeskySolution(MTM, z, MTV, nCols); //for (int i = 0; i < z.Length; ++i) //{ // Console.WriteLine(z[i]); //} //Accuracy for training data double accuracy = 0.0; int count = 0; double min = _y[_nt - 1]; double max = _y[_nt - 1]; for (int i = _nt - 1; i < size; ++i) { double y = 0; for (int k = 0; k < NLinear; ++k) { y += _x[i + k] * z[k]; } accuracy += (_y[i + NLinear] - y) * (_y[i + NLinear] - y); ++count; if (_y[i] < min) min = _y[i]; if (_y[i] > max) max = _y[i]; } accuracy /= (double)(count); accuracy = Math.Sqrt(accuracy); Console.WriteLine("Accuracy for training data (linear model, input only) {0:0.0000}", accuracy / (max - min)); //Accuracy for unseen data size = _N - NLinear; accuracy = 0.0; count = 0; min = _y[_nt]; max = _y[_nt]; for (int i = _NTraining; i < size; ++i) { double y = 0; for (int k = 0; k < NLinear; ++k) { y += _x[i + k] * z[k]; } accuracy += (_y[i + NLinear] - y) * (_y[i + NLinear] - y); ++count; if (_y[i] < min) min = _y[i]; if (_y[i] > max) max = _y[i]; } accuracy /= (double)(count); accuracy = Math.Sqrt(accuracy); Console.WriteLine("Accuracy for unseen data (linear model, input only) {0:0.0000}", accuracy / (max - min)); } static void Main(string[] args) { GenerateOperator(); GenerateX(); ComputeY(); AddCorrelatedError(20); IdentifyTraining(); Validate(); TestDifferentSize(); TestLinearModelInput(); } } }