using System;

namespace Science.Mathematics.LinearAlgebra
{
    public class Projection
    {
        public Projection(Subspace subspace, Vector b)
        {
            Matrix A = new Matrix(subspace.Basis[0].Size, subspace.Basis.Length);
            for (int i = 0; i < A.SizeOfColumn; i++)
                for (int j = 0; j < A.SizeOfRow; j++)
                    A[j, i] = subspace.Basis[i][j];

            Matrix AtA = A.Transpose*A;
            Matrix AA = new Matrix(AtA.SizeOfRow, AtA.SizeOfColumn);
            for (int i = 0; i < AA.SizeOfRow; i++)
                for (int j = 0; j < AA.SizeOfColumn; j++)
                    AA[i, j] = AtA[i, j];
            
            P = A*AA.Inverse*A.Transpose;
            p = P * b;
            e = b - p;
        }
        public Projection(Vector a, Vector b)
        {
            double ata = a * a;
            Matrix aat = new Matrix(a.Size,a.Size);
            for (int i = 0; i < aat.SizeOfRow; i++)
                for (int j = 0; j < aat.SizeOfColumn; j++)
                    aat[i, j] = a[i]*a[j];

            P = aat *(1.0/ata);
            p = P * b;
            e = b - p;
        }
        public Projection(Vector a)
        {
            double ata = a * a;
            Matrix aat = new Matrix(a.Size, a.Size);
            for (int i = 0; i < aat.SizeOfRow; i++)
                for (int j = 0; j < aat.SizeOfColumn; j++)
                    aat[i, j] = a[i] * a[j];

            P = aat * (1.0 / ata);
        }
        private Vector p;
        public Vector OntoSubspace
        {
            get { return p; }
        }
        private Vector e;
        public Vector Error
        {
            get { return e; }
        }
        private Matrix P;
        public Matrix ProjectionMatrix
        {
            get { return P; }
        }
    }
}