#pragma once

template<class T>
void PCA<T>::InitFromPointMatrix(DenseMatrix<T> &B)
{
    const UINT PointCount = B.ColCount();
    const UINT Dimension = B.RowCount();

    Console::WriteLine(String("Initializing PCA, ") + String(B.ColCount()) + String(" points, ") + String(B.RowCount()) + String(" dimensions"));
    
    _Means.Allocate(Dimension);
    _Means.Clear(0.0);
    for(UINT PointIndex = 0; PointIndex < PointCount; PointIndex++)
    {
        for(UINT DimensionIndex = 0; DimensionIndex < Dimension; DimensionIndex++)
        {
            _Means[DimensionIndex] += B.Cell(DimensionIndex, PointIndex);
        }
    }
    for(UINT DimensionIndex = 0; DimensionIndex < Dimension; DimensionIndex++)
    {
        _Means[DimensionIndex] /= PointCount;
    }
    for(UINT PointIndex = 0; PointIndex < PointCount; PointIndex++)
    {
        for(UINT DimensionIndex = 0; DimensionIndex < Dimension; DimensionIndex++)
        {
            B.Cell(DimensionIndex, PointIndex) -= _Means[DimensionIndex];
        }
    }

    DenseMatrix<T> C;
    Console::WriteLine("Building correlation matrix...");
    DenseMatrix<T>::MultiplyMMTranspose(C, B);
    DenseMatrix<T>::MultiplyInPlace(C, T(1.0) / T(PointCount));
    InitFromCorrelationMatrix(C);
}

template<class T>
void PCA<T>::InitFromDensePoints(const Vector<const T*> &Points, UINT Dimension)
{
    Console::WriteLine(String("Initializing PCA, ") + String(Points.Length()) + String(" points, ") + String(Dimension) + String(" dimensions"));
    
    const UINT PointCount = Points.Length();
    DenseMatrix<T> B(Dimension, PointCount);
    
    _Means.Allocate(Dimension);
    _Means.Clear(0.0);
    for(UINT PointIndex = 0; PointIndex < PointCount; PointIndex++)
    {
        const T *CurPoint = Points[PointIndex];
        for(UINT DimensionIndex = 0; DimensionIndex < Dimension; DimensionIndex++)
        {
            _Means[DimensionIndex] += CurPoint[DimensionIndex];
        }
    }
    for(UINT DimensionIndex = 0; DimensionIndex < Dimension; DimensionIndex++)
    {
        _Means[DimensionIndex] /= PointCount;
    }
    for(UINT PointIndex = 0; PointIndex < PointCount; PointIndex++)
    {
        const T *CurPoint = Points[PointIndex];
        for(UINT DimensionIndex = 0; DimensionIndex < Dimension; DimensionIndex++)
        {
            B[DimensionIndex][PointIndex] = CurPoint[DimensionIndex] - _Means[DimensionIndex];
        }
    }

    DenseMatrix<T> C;
    Console::WriteLine("Building correlation matrix...");
    DenseMatrix<T>::MultiplyMMTranspose(C, B);
    DenseMatrix<T>::MultiplyInPlace(C, T(1.0) / T(PointCount));
    InitFromCorrelationMatrix(C);
}

template<class T>
void PCA<T>::InitFromCorrelationMatrix(const DenseMatrix<T> &C)
{
    const UINT Dimension = C.RowCount();
    const bool OutputMatlabAndExit = false;
    if(OutputMatlabAndExit)
    {
        Console::WriteLine("Outputting MATLAB matrix and exiting...");
        C.SaveToMATLAB("C:\\MATLAB7\\work\\PCA.txt");
        _Means.SaveToASCIIFile("C:\\JReaderData\\RawMeanVector.txt");
        exit(0);
    }
    else
    {
        Console::WriteLine("Computing eigensystem...");
        C.EigenSystem(_Eigenvalues, _Eigenvectors);
        Console::WriteLine(C.EigenTest(_Eigenvalues, _Eigenvectors));

        //_Eigenvectors.SaveToMATLAB("C:\\MATLAB7\\work\\MyEigenvectors.txt");
        //_Eigenvalues.SaveToASCIIFile("C:\\MATLAB7\\work\\MyEigenvalues.txt");
    }

    FinalizeFromEigenSystem();
}

template<class T>
void PCA<T>::InitFromMATLAB(UINT Dimension, const String &EigenvectorFilename, const String &EigenvalueFilename, const String &MeanFilename)
{
    Console::WriteLine(String("Initializing PCA from MATLAB eigensystem files..."));
    //M = load('PCA.txt');
    //[Eigenvectors, Eigenvalues] = eig(M);
    //dlmwrite('Eigenvectors.txt', Eigenvectors, ' ');
    //dlmwrite('Eigenvalues.txt', diag(Eigenvalues), ' ');

    _Means.LoadFromASCIIFile("C:\\JReaderData\\RawMeanVector.txt");
    _Eigenvectors.Allocate(Dimension, Dimension);
    _Eigenvalues.Allocate(Dimension);

    ifstream EigenvectorFile(EigenvectorFilename.CString());
    for(UINT Row = 0; Row < Dimension; Row++)
    {
        for(UINT Col = 0; Col < Dimension; Col++)
        {
            EigenvectorFile >> _Eigenvectors.Cell(Row, Dimension - 1 - Col);
        }
    }

    ifstream EigenvalueFile(EigenvalueFilename.CString());
    for(UINT Index = 0; Index < Dimension; Index++)
    {
        EigenvalueFile >> _Eigenvalues[Dimension - 1 - Index];
    }

    //_Eigenvectors.SaveToMATLAB("C:\\MATLAB7\\work\\RepeatedEigenvectors.txt");
    //_Eigenvalues.SaveToASCIIFile("C:\\MATLAB7\\work\\RepeatedEigenvalues.txt");

    FinalizeFromEigenSystem();
}

template<class T>
void PCA<T>::FinalizeFromEigenSystem()
{
    const UINT Dimension = _Eigenvalues.Length();
    double EigenvalueSum = 0.0;
    for(UINT DimensionIndex = 0; DimensionIndex < Dimension; DimensionIndex++)
    {
        EigenvalueSum += Math::Max(T(0.0), _Eigenvalues[DimensionIndex]);
    }
    double CumulativeEnergy = 0.0;
    for(UINT DimensionIndex = 0; DimensionIndex < Dimension; DimensionIndex++)
    {
        CumulativeEnergy += _Eigenvalues[DimensionIndex];
        Console::WriteLine(String("Energy at ") + String(DimensionIndex + 1) + String(" terms: ") + String(CumulativeEnergy / EigenvalueSum * 100.0f) + String("%"));
    }
}

template<class T>
UINT PCA<T>::ReducedDimension(double EnergyPercent)
{
    double EigenvalueSum = 0.0;
    for(UINT DimensionIndex = 0; DimensionIndex < _Eigenvalues.Length(); DimensionIndex++)
    {
        EigenvalueSum += _Eigenvalues[DimensionIndex];
    }
    double CumulativeEnergy = 0.0;
    for(UINT DimensionIndex = 0; DimensionIndex < _Eigenvalues.Length(); DimensionIndex++)
    {
        CumulativeEnergy += _Eigenvalues[DimensionIndex];
        if(CumulativeEnergy / EigenvalueSum >= EnergyPercent)
        {
            return DimensionIndex + 1;
        }
    }
    return _Eigenvalues.Length();
}

template<class T>
void PCA<T>::Transform(Vector<T> &Result, const Vector<T> &Input, UINT ReducedDimension)
{
    if(Result.Length() != ReducedDimension)
    {
        Result.Allocate(ReducedDimension);
    }
    Transform(Result.CArray(), Input.CArray(), ReducedDimension);
}

template<class T>
void PCA<T>::InverseTransform(Vector<T> &Result, const Vector<T> &Input)
{
    const UINT Dimension = _Means.Length();
    Result.Allocate(Dimension);
    InverseTransform(Result.CArray(), Input.CArray(), Input.Length());
}

template<class T>
void PCA<T>::Transform(T *Result, const T *Input, UINT ReducedDimension)
{
    const UINT Dimension = _Means.Length();
    for(UINT Row = 0; Row < ReducedDimension; Row++)
    {
        T Total = 0.0;
        for(UINT Index = 0; Index < Dimension; Index++)
        {
            Total += _Eigenvectors[Row][Index] * (Input[Index] - _Means[Index]);
        }
        Result[Row] = Total;
    }
}

template<class T>
void PCA<T>::InverseTransform(T *Result, const T *Input, UINT ReducedDimension)
{
    UINT Dimension = _Means.Length();
    for(UINT Col = 0; Col < Dimension; Col++)
    {
        T Total = 0.0;
        for(UINT Index = 0; Index < ReducedDimension; Index++)
        {
            Total += _Eigenvectors[Index][Col] * Input[Index];
        }
        Result[Col] = Total + _Means[Col];
    }
}