//
// Libregresion -- A library to help you with
// your linear regresions, interpolations, systems solving, etc...
//
// (c) 2003, Jaime Anguiano Olarra <jaime@gnome.org>
//

using System;
using System.IO;
using System.Collections;

namespace Libregresions {

public class Libregresion {
	static double[] xarr; // x data 
	static double[] yarr; // y data
	static double[] x_sums; // the x^j sums
	static double y_sum; // the y sum
	static double[] xy_sums; // the y*x^j sums
	static double[,] coeffs; // coeffs matrix for the linear systems methods
	static int order; // order of the polinomy for the regresions
	static bool saved = true; // data saved?
	static readonly double PI = 3.141592653589;
	
	// Do a simple linear regresion
	public static double[] Regres (double[] xarr, double[] yarr)
	{
		order = 1;
		return Regres (xarr, yarr, 1, 15000);
	}

	// Libregresion for a user choosen number of iterations and order of the polinomy
	public static double[] Regres (double[] xarr, double[] yarr, int n, long iter)
	{
		order = n;
		if (xarr.Length != yarr.Length)
			throw new Exception ();

		// We do the sums for x^j, y and y*x^j
		x_sums = new double[2*n];
		y_sum = Sum(yarr);
		xy_sums = new double[n];

		for (int i=0; i < 2*n; i++) 
			x_sums[i] = Sum(ToNth(xarr, i));
			
		for (int j=0; j < n; j++) 
			xy_sums[j] = Sum (Multiply (ToNth(xarr, j), yarr));

		// We set the coefficients matrix for the linear system obtained by
		// the regresion method
		coeffs = new double[order, order];
		
		for (int i=0; i < order; i++) {
			for (int j=0; j < order; j++) {
				coeffs[i,j] = x_sums[i+j];
			}
		}

		return GaussSeidel (coeffs, xy_sums, iter);
	}

	// Returns the sum of all the elements of the array passed to it
	public static double Sum (double[] arr)
	{
		double sum = 0.0;
		
	        for (int i=0; i < arr.Length; i++) {
			sum += arr[i];
		}
		return sum;
	}

	// This method takes the array passed and it powers each of its
	// elements to the power passed, n, (counting from 1). Then it
	// returns a new array with the new values.
	// F.ex. if you pass [2,3] you get [4,9].
	public static double[] ToNth (double[] arr, int n)
	{
		double[] elems = new double[arr.Length];
		double pwr;
		for (int i=0; i < arr.Length; i++) {
			pwr = 1;
			for (int j=0; j < n; j++) {
				pwr *= arr[i];
			}
			elems[i] = pwr;
		}

		return elems;
	}

	// Returns an array which each element is calculated by
	// multiplying each the corresponding elements from the
	// two arrays passed to it.
	public static double[] Multiply (double[] xarr, double[] yarr)
	{
		double[] zarr = new double[xarr.Length];
		
		for (int i=0; i < xarr.Length; i++) {
			zarr[i] = xarr[i] * yarr[i];
		}

		return zarr;
	}

	// Returns the matrix generated by multiplying a * b
	public static double[,] Multiply (double[,] a, double[,] b, uint n)
	{
		double[,] c;
		double y;

		c = new double[n,n];
			
		for (uint i=0; i < n; ++i) {
			for (uint j=0; j < n; ++j) {
				y = 0;
				for (uint k=0; k < n; ++k)
					y += a[i,k]*b[k,j];
				c[i,j] = y;
			}
		}

		return c;
	}

	// Returns the vector obtained by multiplying m (matrix) for v (vector)
	public static double[] Multiply (double[,] m, double[] v)
	{
		double y;
		double[] s;

		s = new double[v.Length];
		
		for (uint i=0; i < v.Length; ++i) {
			y = 0;
			for (uint j=0; j < v.Length; ++j)
				y += m[i,j]*v[j];
			s[i] = y;
		}

		return s;
	}

	// Returns the matrix passed multiplied by the constant passed
	public static double[,] Multiply (double[,] x, double c, uint order)
	{
		for (uint i=0; i < order; ++i)
			for (uint j=0; j < order; ++j)
				x[i,j] *= c;
		return x;
	}

	// Returns the matrix passed multiplied by the constant passed
	public static double[,] Multiply (double c, double[,] x, uint order)
	{
		return Multiply (x, c, order);
	}

	// Returns the matrix obtained by multiplying v (vector) for m (matrix)
	public static double[,] Multiply (double[] v, double[,] m)
	{
		// Not implemented yet
		return m;
	}
	
	// Returns the traspose of the matrix passed as argument
	public static double[,] Traspose (double[,] m, uint n)
	{
		double tmp;
		double[,] x;

		x = new double[n,n];
		
		for (uint j=0; j < n; ++j) 
			for (uint i=0; i < n; ++i) {
				if (i != j)
					x[i,j] = m[j,i];
				else
					x[i,i] = m[i,i];
			}
		return x;
	}

	// Returns the adjunt of the element given by (n,m) in the matrix x
	public static double[,] AdjointMatrix (double[,] x, int n, int m, uint order)
	{
		double[,] adj;
		uint r = 0, c = 0;
		
		adj = new double[order - 1, order - 1];
		for (uint i=0; i < order; ++i) {
			if (i == n)
				continue;
			for (uint j=0; j < order; ++j) {
				if (c == (order - 1)) { 
					++r;
					c = 0;
				}
				if (j == m) 
					continue;
				adj[r,c] = x[i,j];
				++c;
			}
		}
		return adj;
	}

	// Returns the determinant of an adjoint matrix element
	public static double Adjoint (double[,] x, int n, int m, uint order)
	{
		return (Sign(n+m) * Determinant (AdjointMatrix(x, n, m, order), order - 1));
	}

	// Returns the adjoint matrix of a complete matrix
	public static double[,] Adjoint (double[,] x, uint order)
	{
		double[,] matrix;

		matrix = new double[order, order];

		for (int i=0; i < order; ++i)
			for (int j=0; j < order; ++j) 
				matrix[i,j] = Adjoint(x, i, j, order);
		return matrix;
	}
	
	// Returns the inverse of the matrix passed
	public static double[,] Inverse (double[,] x, uint order)
	{
		double det;
		if ((det = Determinant(x, order)) == 0)
			ThrowException ("Error: matrix has no inverse.", x);
		return (Multiply (Traspose(Adjoint(x, order), order), det, order));
	}
	
	// Returns the determinant of the matrix passed
	public static double Determinant (double[,] x, uint order)
	{
		double d;
		double[,] matrix;
		double[,] tmp;
		int morder;

		morder = (int) order;
		matrix = Copy (x, (uint) morder);

		if (morder == 2)
			return ((matrix[0,0]*matrix[1,1]) - (matrix[0,1]*matrix[1,0]));
		d = 0;
		try {
			for (int i=0; i < morder; ++i)
					d += Sign(i) * matrix[0,i] * Determinant(tmp = AdjointMatrix(matrix, 0, i, (uint) morder), (uint) morder - 1);
		} catch (Exception e) {
			Console.WriteLine (e);
		}
		return d;
	}

	// Returns the trace of the matrix
	public static double Trace (double[,] x, uint order)
	{
		double t;
		t = 0;
		for (uint i=0; i < order; ++i)
			t += x[i,i];
		return t;
	}

	// Returns the sign of the integer passed
	public static int Sign (int s)
	{
		int result;
		result = ((s%2) == 0) ? 1 : -1;
		return result;
	}

	// Returns a copy of the matrix passed
	public static double[,] Copy (double[,] x, uint order)
	{
		double[,] matrix;

		matrix = new double[order, order];
		for (uint i=0; i < order; ++i)
			for (uint j=0; j < order; ++j)
				matrix[i,j] = x[i,j];
		return matrix;
	}
	
	// Prints the content of an array to the console.
	public static void Show (double[] arr)
	{
		for (int i=0; i < arr.Length; i++) 
			Console.Write (arr[i] + " ");
		Console.WriteLine ();
	}

	// Converts an array to a single string form
	public static string ArrToString (double[] arr)
	{
		string s = null;
		for (int i=0; i < arr.Length; ++i)
			s += (arr[i]).ToString() + " ";
		return s;
	}

	// Prints the content of a matrix of n rows, m columns to the
	// console.
	public static void Show (double[,] matrix, int n, int m)
	{
		for (int i=0; i < n; i++) {
			for (int j=0; j < m; j++) {
				Console.Write (matrix[i,j] + " ");
			}
			Console.WriteLine ();
		}
		Console.WriteLine ();
	}

	// Prints the content of a solution array with an improved
	// look.
	public static void ShowSols (double[] arr)
	{
		for (int i=0; i < arr.Length; i++)
			Console.WriteLine ("x" + i + ": {0}", arr[i]);
	}

	public static string SolsToString (double[] arr)
	{
		string s = null;
		for (int i=0; i < arr.Length; ++i)
			s += "x" + (i).ToString() + ": " + (arr[i]).ToString() + "  ";
		return s;
	}
	
	// Resolvs a linear system by Gauss-Seidel's method given the
	// coefficient matrix for the unknown variables, a solution vector
	// and the number of iterations it should use.
	public static double[] GaussSeidel (double[,] matrix, double[] solVctr, long iter)
	{
		double q = 0;
		double[] sols = new double[solVctr.Length];

		for (int t=0; t < sols.Length; t++)
			sols[t] = 0;

		for (int k=0; k < iter; k++) {
			for (int i=0; i < solVctr.Length; i++) {
				for (int j=0; j < solVctr.Length; j++) {
					if (i != j) 
						q += matrix[i,j]*sols[j];
				}
				sols[i] = (solVctr[i] - q)/matrix[i,i];
				q = 0;
			}
		}
		return sols;
	}

	// This needs a second look. Something is going wrong.
	public static double[] GaussJordan (double[,] matrix, double[] solVctr)
	{
		// Old code erased. You better start rewritting it with the new
		// set of functions (Unitary, etc)...
		return solVctr;
	}

	// Pivoting of matrix and returns the order vector
	public static double[] Pivoting (ref double[,] m, uint order)
	{
		Console.WriteLine ("Starting pivoting");
		double[] ord;
		uint row;
		double temp;
		
		ord = new double[order];
		
		for (uint i=0; i < order; ++i)
			ord[i] = i;

		for (uint i=0; i < order; ++i) {
			row = i;
			for (uint p=i; p < order; ++p) {
				if (Math.Abs(m[row, i]) < Math.Abs(m[p, i]))
					row = p;
			}

			if (row != i) {
				for (uint j=0; j < order; ++j) {
					temp = m[i,j];
					m[i,j] = m[row,j];
					m[row,j] = temp;
				}
				temp = ord[i];
				ord[i] = ord[row];
				ord[row] = temp;
			}
			Show (m, (int) order, (int) order);
		}
		Console.WriteLine ("Finished pivoting");
		return ord;
	}

	// Check for singularity of the inverse of matrix
	public static bool HasInverse (double[,] matrix, uint n)
	{
		return HasInverse (matrix, n, 0);
	}
		
	// Check for singularity of the inverse of matrix for minimum value epsil
	public static bool HasInverse (double[,] matrix, uint n, double epsil)
	{
		bool response = false;
		for (uint i=0; i < n; ++i) 
			response = (Math.Abs(matrix[i,i]) < epsil) ? false : true;
		return response;
	}
	
	// Return a unitary matrix of order n
	public static double[,] Unitary (uint n)
	{
		double[,] m;
		m = new double[n,n];

		for (int i=0; i < n; ++i)
			for (int j=0; j < n; ++j) {
				if (i != j)
					m[i,j] = 0;
				else
					m[i,i] = 1; 
			}
		return m;
	}

	// Returns the 'y' interpolated value for the table given by the arrays
	// passed, as well as the iteration and the error
	public static double[] Interpolation (double[] xarr, double[] yarr, double xval)
	{
		int m; // m is the number of intervals (n of points - 1)
		double y, er = 0;
		double n = xarr.Length - 1; // Grade of the interpolation polinomy
		double[] l = new double[xarr.Length];

		m = (int) n;
		for (int i=0; i < n; ++i) {
			y = 1;
			for (int j=0; j < n; ++j) {
				if ((i != j) && ((xarr[j] - xarr[i]) != 0)) {
					y *= ((xval - xarr[j]) / (xarr[i] - xarr[j]));
				}
			}
			l[i] = y;
		}
		y = 0;
		for (int i=0; i < n; ++i)
			y += yarr[i]*l[i];
		
		return new double[3] {y, m, er};
	}

	// Get the factorial of the argument
	public static long Factorial (uint n)
	{
		if (n == 1 || n == 0)
			return 1;

		return (n * Factorial(n - 1));
	}

        // Calculates the binomial coefficient (combinatory n,m)
	public static double Combination (uint n, uint m)
	{
		return (Factorial(n)/(Factorial(m)*Factorial(n - m)));
	}

	// Calculates the variation of n elements taken by m
	public static long Variation (uint n, uint m)
	{
		return (Factorial(n)/Factorial(n - m));
	}

	//*** WARNING ***
	//**
	//** Look for a better, general definition of Gamma if any
	// Calculates the gamma of (2*n - 1) , n => 1
	public static double Gamma (uint n)
	{
		double g;
		g = 1;
		for (uint i=1; i <= n; ++i)
			g *= (2*i - 1);
		return (Math.Sqrt(PI)*(1/Math.Pow(2,n))*g);
	}

	// Calcultates the beta of n, m
	public static double Beta (uint n, uint m)
	{
		return (Gamma(n)*Gamma(m)/Gamma(n + m));
	}
	
	// Launches a LibregresionException instance with the message passed
	public static void ThrowException (string msg)
	{
		throw new LibregresionException (msg);
	}

	// Launches a LibregresionException instance with the message passed and the matrix that caused it
	public static void ThrowException (string msg, double[,] x)
	{
		throw new LibregresionException (msg, x);
	}
	
	// Asks for X/Y style input, that is, it allows you to retrieve the x_data
	// followed by the y_data
	public static void AskForInput ()
	{
		ArrayList xa = new ArrayList();
		saved = false;
		string s = null;
		bool yep = true;
		int i = 0;
		Console.WriteLine ("Enter the independent variable (X) points: ");
		Console.WriteLine ("(enter 'q' to stop)");
		while (yep) {
			Console.Write ("x" + i + ": ");
			s = Console.ReadLine ();
			if (String.Compare (s, "q") == 0) {
				yep = false;
				continue;
			}
			xa.Add (Convert.ToDouble(s));
			i++;
		}

		xarr = new double[xa.Count];
		yarr = new double[xa.Count];
		
		Console.WriteLine ("Enter the dependent variable (Y) points: ");
		for (int j=0; j < xa.Count; j++) {
			Console.Write ("y" + j + ": ");
			s = Console.ReadLine ();
			yarr[j] = Convert.ToDouble(s);
		}

		for (int j=0; j < xa.Count; j++) {
			xarr[j] = (double) xa[j];
		}
	}

	// Main method for Console versions of the program.
#if (CONSOLE)
	public static void Main ()
	{
		Console.WriteLine ("\n\n\t\t\tReAgresion 1.0\n\n");
		Console.WriteLine ("\tThe definitive calculator for physicists and\n" +
				   "\tthe people that need to do lots of regressions,\n" +
				   "\tinterpolations, solving systems...\n");
		Console.WriteLine ("\tDistributed under the GNU General Public License v 2.0\n");
		Console.WriteLine ("\t(c) 2003, Jaime Anguiano Olarra <jaime@gnome.org>\n\n");

		bool yep = true;
		
		while (yep) {
			PrintMenu ();
			switch (Console.ReadLine()) {
			case "1":
				DoLibregresion ();
				break;
			case "2":
				DoInterpolation ();
				break;
			case "3":
				DoGaussSeidel ();
				break;
			case "4":
				Save ();
				break;
			case "5":
				Modify ();
				break;
			case "6":
				Show ();
				break;
			case "7":
				if (!saved) {
					Console.WriteLine ("Your data hasn't been saved!\n" +
							   "Do you want to save it now? (y/n)");
					if (String.Compare ((Console.ReadLine()).ToUpper(),"Y") == 0) {
						Save ();
						saved = true;
					}
				}
					
				yep = false;
				break;
			}
		}
	}
#endif 
	
	public static void PrintMenu ()
	{
		Console.WriteLine ("\n\nEnter your choice: ");
		Console.WriteLine ("\t1. Libregresion");
		Console.WriteLine ("\t2. Interpolation");
		Console.WriteLine ("\t3. Gauss-Seidel");
		Console.WriteLine ("\t4. Save data");
		Console.WriteLine ("\t5. Modify data");
		Console.WriteLine ("\t6. Show data");
		Console.WriteLine ("\t7. Quit\n\n");
	}
	
	public static void DoLibregresion ()
	{
		int ord;
		long itr;
		double[] sols;
		
		AskForInput ();
		Console.Write ("Enter the order of the polinomy to use: ");
		ord = Convert.ToInt32 (Console.ReadLine ()) + 1;
		Console.WriteLine ();
		Console.WriteLine ("Note: This program uses the Gauss-Seidel method by default to calculate\n" +
				   "the polinomy coefficients. Normally you should get a nice convergency and\n" +
				   "you shouldn't enter numbers above 15000.\n\n");
		Console.Write ("Enter the number of iterations: ");
		itr = (long) Convert.ToDouble (Console.ReadLine ());

		sols = Regres (xarr, yarr, ord, itr);
		
		ShowSols (sols);
		Console.Write ("Save results? (y/n): ");
		if (String.Compare ((Console.ReadLine ()).ToUpper(), "Y") == 0) {
			Console.Write ("Enter the file to use: ");
			StreamWriter sw = new StreamWriter (Console.ReadLine());
			sw.WriteLine ("Results obtained by regresion:\n\nThe vectors entered:\n(x data)");
			WriteVector (sw, xarr);
			sw.WriteLine ("(y data):");
			WriteVector (sw, yarr);
			sw.WriteLine ("The Gauss-Seidel solutions for them after " + itr + " iterations:");
			WriteVector (sw, sols);
			sw.Flush ();
			sw.Close ();
			sw = null;
			saved = true;
		}
	}

	public static void DoInterpolation ()
	{
		double xval;
		double[] results;
		AskForInput ();
		Console.Write ("Enter the value for the one you want to interpolate (x): ");
		xval = Convert.ToDouble (Console.ReadLine());
		results = Interpolation (xarr, yarr, xval);
		Show (results);
	}

	public static void DoGaussSeidel ()
	{
		string s = null;
		int order;
		double[,] matrix;
		double[] sols;
		double[] gs_sols;
		int iter;
		
		Console.Write ("Enter the number of rows: ");
		order = Convert.ToInt32 (Console.ReadLine ());
		matrix = new double[order, order];
		Console.WriteLine ("\nEnter the matrix elements:\n");
		for (int i=0; i < order; i++) {
			for (int j=0; j < order; j++) {
				Console.Write ("a" + i + j + ": ");
				matrix[i,j] = Convert.ToDouble (Console.ReadLine());
			}
		}

		Console.WriteLine ("\nEnter the solution vector:");
		sols = new double[order];
		for (int k=0; k < order; k++) {
			Console.Write ("s" + k + ": ");
			sols[k] = Convert.ToDouble (Console.ReadLine ());
		}

		Console.Write ("Enter the number of iterations to use: ");
		iter = Convert.ToInt32 (Console.ReadLine ());
		gs_sols = GaussSeidel (matrix, sols, iter);
		ShowSols (gs_sols);
		saved = false;
		Console.Write ("Save results? (y/n): ");
		if (String.Compare ((Console.ReadLine ()).ToUpper(), "Y") == 0) {
			Console.Write ("Enter the file to use: ");
			StreamWriter sw = new StreamWriter (Console.ReadLine());
			sw.WriteLine ("Results obtained by Gauss-Seidel:\n\nThe matrix entered:");
			WriteMatrix (sw, matrix, order);
			sw.WriteLine ("The solution vector entered:");
			WriteVector (sw, sols);
			sw.WriteLine ("The Gauss-Seidel solutions for them after " + iter + " iterations:");
			WriteVector (sw, gs_sols);
			sw.Flush ();
			sw.Close ();
			sw = null;
			saved = true;
		}
	}

	public static void Save ()
	{
		StreamWriter sw = new StreamWriter ("data.txt");
		sw.WriteLine ("X Data:\n");
		for (int i=0; i < xarr.Length; i++)
			sw.WriteLine ("x" + i + ": " + xarr[i]);
		sw.WriteLine ("\nY Data:\n");
		for (int j=0; j < yarr.Length; j++)
			sw.WriteLine ("y" + j + ": " + yarr[j]);
		sw.Flush ();
		sw.Close ();
		sw = null;
		saved = true;
	}

	public static void WriteVector (StreamWriter sw, double[] vctr)
	{
		for (int i=0; i < vctr.Length; i++)
			sw.WriteLine ("s" + i + ": " + vctr[i]);
		sw.WriteLine ();
	}

	public static void WriteMatrix (StreamWriter sw, double[,] mtrx, int order)
	{
		for (int i=0; i < order; i++) {
			for (int j=0; j < order; j++) {
				sw.Write (mtrx[i,j] + " ");
				if (j == order - 1)
					sw.WriteLine ();
			}
		}
		sw.WriteLine ();
	}
	
	public static void Modify ()
	{
		Console.Write ("\n\nDo you want to modify an 'X' or an 'Y' (x/y)?: ");
		if (String.Compare ((Console.ReadLine ()).ToUpper(), "X") == 0) {
			modifyx:
			Console.Write ("\nEnter the index (starting at 0): ");
			int i = Convert.ToInt32 (Console.ReadLine());
			Console.WriteLine ("The current value of x" + i + " is: " + xarr[i]);
			Console.Write ("\nEnter the new value for x" + i + ": ");
			xarr[i] = Convert.ToDouble (Console.ReadLine());
			Console.WriteLine ("Modify another value?: (y/n)");
			if (String.Compare ((Console.ReadLine ()).ToUpper(), "Y") == 0)
				goto modifyx;
		} else {
			modifyy:
			Console.Write ("\nEnter the index (starting at 0): ");
			int j = Convert.ToInt32 (Console.ReadLine());
			Console.WriteLine ("The current value of y" + j + " is: " + yarr[j]);
			Console.Write ("\nEnter the new value for y" + j + ": ");
			yarr[j] = Convert.ToDouble (Console.ReadLine());
			Console.WriteLine ("Modify another value?: (y/n)");
			if (String.Compare ((Console.ReadLine ()).ToUpper(), "Y") == 0)
				goto modifyy;
		}

	}

	public static void Show ()
	{
		Console.WriteLine ("\nThe current data values are:");
		Console.WriteLine ("\nX\t\tY\n");
		for (int i=0; i < xarr.Length; i++) {
			Console.WriteLine (xarr[i] + "\t\t" + yarr[i]);
		}
	}
}

[Serializable]
public class LibregresionException : System.Exception {
	private double[,] matrix;
	
	public LibregresionException () {}

	public LibregresionException (string msg) : base (msg)
	{
	}

	public LibregresionException (string msg, double[,] matrix) : base (msg)
	{
		this.matrix = matrix;
	}

	public double[,] Matrix { get { return matrix; } }
}
}