mtupca.cpp

Go to the documentation of this file.

#include "PCA.h"
#include "statutils.h"
#include "VecMat.h"
#include "TimeSeriesData.h"
#include <iostream>
#include <fstream>
#include <vector>

int main()
{
        TimeSeriesData TsData;
        std::string mtuname;
        
        std::cout << "Data file name: ";
        std::cin >> mtuname;
        
        
        TsData.GetData(mtuname);
        
        
        const size_t nobs = TsData.GetData().Size();
        const size_t nchan = 4;
        gplib::rmat input(nchan,nobs);
        gplib::cmat evec(nchan,nchan);
        gplib::cvec eval(nchan);
        
        for (unsigned int i = 0; i < nobs; ++i)
        {
                input(0,i) = TsData.GetData().GetEx().GetData().at(i);
                input(1,i) =    TsData.GetData().GetEy().GetData().at(i);
                input(2,i) = TsData.GetData().GetHx().GetData().at(i);
                input(3,i) = TsData.GetData().GetHy().GetData().at(i);
        }
    
    PCA(input,evec,eval);
        
        gplib::cmat wmat(WhiteMat(evec, eval));
        std::cout << "pca evec: " << evec << std::endl;
        std::cout << "pca eval: " << eval << std::endl;
        std::cout << "pca WhM: " << wmat << std::endl;
        std::cout << "pca DeWhM: " << DeWhiteMat(evec, eval) << std::endl;
        
        gplib::cmat output(nchan,nobs);
        noalias(output) = prec_prod(wmat,input);
        for (unsigned int i = 0; i < nobs; ++i)
        {
                TsData.GetData().GetEx().GetData().at(i) = output(0,i).real();
                TsData.GetData().GetEy().GetData().at(i) = output(1,i).real();
                TsData.GetData().GetHx().GetData().at(i) = output(2,i).real();
                TsData.GetData().GetHy().GetData().at(i) = output(3,i).real();
        }
        TsData.WriteAsBirrp(mtuname+".pca");
}

---