using System;

namespace Science.Mathematics.LinearAlgebra
{
    public class Nullspace : VectorSpace
    { 
        private Matrix AA;
        public Nullspace(Matrix A)
        {
            AA = A;
            FindSolution();
            this.Basis = s;
            this.Dimension = AA.SizeOfColumn - nofpivot;
        }
        public Matrix NullspaceMatrix
        {
            get
            {
                Matrix A = new Matrix(s);
                return A;
            }
        }
        private Vector[] s;
        private int nofpivot;
        private bool[] pivot;
        public bool[] PivotColumnQ
        {
            get { return pivot; }
        }
        private int[] pivotrow;
        public int[] PivotRow
        {
            get { return pivotrow; }
        }
        private void FindSolution()
        {
            FactorizationEAeqR obj = new FactorizationEAeqR(AA);
            obj.Compute();
            nofpivot = 0;
            pivot = new bool[AA.SizeOfColumn];
            pivotrow = new int[AA.SizeOfColumn];
            for(int i = 0; i < AA.SizeOfRow; i++)
                for(int j = 0; j < AA.SizeOfColumn; j++)
                    if (obj.ReducedRowEchelonForm[i, j] > 1.0E-12)
                    {
                        nofpivot++;
                        pivot[j] = true;
                        pivotrow[j] = i;
                        break;
                    }
            if (AA.SizeOfColumn == nofpivot)
            {
                s = new Vector[1];
                s[0] = new Vector(AA.SizeOfColumn);
            }
            else
            {
                s = new Vector[AA.SizeOfColumn - nofpivot];
                int nn = 0;
                for (int j = 0; j < AA.SizeOfColumn; j++)
                {
                    if (!pivot[j])
                    {
                        double[] x = new double[AA.SizeOfColumn];
                        x[j] = 1.0;
                        for (int i = 0; i < AA.SizeOfColumn; i++)
                        {
                            if (pivot[i])
                            {
                                x[i] = -obj.ReducedRowEchelonForm[pivotrow[i], j];
                            }
                        }
                        s[nn] = new Vector(x);
                        nn++;
                    }
                }
            }
        }
        public Vector LinearCombination(double[] x)
        {
            Vector xx = new Vector(this.Basis[0].Size);
            for(int i = 0; i < x.Length; i++)
                xx += x[i]*this.Basis[i];
            return xx;
        }
    }
}