PCATool.cxx

Go to the documentation of this file.

//***************************************************************************************
//***************************************************************************************
 
#include "NAGASH/PCATool.h"
 
using namespace std;
using namespace NAGASH;
using namespace Eigen;
 
namespace NAGASH
{
    PCATool::PCATool(std::shared_ptr<MSGTool> MSG, int size) : Tool(MSG)
    {
        Size = size;
        Covariance_Matrix = MatrixXd(size, size);
        Observable_Vector = VectorXd(size);
    }
 
    void PCATool::SetObservable(int index, double val)
    {
        auto titleguard = MSGUser()->StartTitleWithGuard("PCATool::SetObservable");
        if (index < 0 || index >= Size)
        {
            MSGUser()->MSG_ERROR("Observable index out of range: ", index);
            return;
        }
        Observable_Vector(index) = val;
    }
 
    void PCATool::SetCovariance(int index_i, int index_j, double cov)
    {
        auto titleguard = MSGUser()->StartTitleWithGuard("PCATool::SetCovariance");
        if (index_i < 0 || index_i >= Size)
        {
            MSGUser()->MSG_ERROR("Observable index out of range: ", index_i);
            return;
        }
        if (index_j < 0 || index_j >= Size)
        {
            MSGUser()->MSG_ERROR("Observable index out of range: ", index_j);
            return;
        }
        if (index_i > index_j)
            SetCovariance(index_j, index_i, cov);
 
        Covariance_Matrix(index_i, index_j) = cov;
        Covariance_Matrix(index_j, index_i) = cov;
        isProcessed = false;
    }
 
    void PCATool::SetUncertainty(int index, double unc)
    {
        auto titleguard = MSGUser()->StartTitleWithGuard("PCATool::SetUncertainty");
        SetCovariance(index, index, unc * unc);
    }
 
    void PCATool::Process()
    {
        auto titleguard = MSGUser()->StartTitleWithGuard("PCATool::Process");
 
        MSGUser()->MSG_INFO("Calculate Eigen Values");
        Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eigensolver(Covariance_Matrix);
        if (eigensolver.info() != Eigen::Success)
            abort();
        EigenValue = eigensolver.eigenvalues();
        MSGUser()->MSG_INFO("Calculate Eigen Vectors");
        Transform_Matrix = eigensolver.eigenvectors();
        MSGUser()->MSG_INFO("Calculate New Correlation Matrix");
        New_Covariance_Matrix = Transform_Matrix.transpose() * Covariance_Matrix * Transform_Matrix;
        MSGUser()->MSG_INFO("Calculate Rotated Observables");
        New_Observable_Vector = Transform_Matrix.transpose() * Observable_Vector;
 
        MSGUser()->MSG_INFO("Calculate Inverse Transform Matrix");
        MatrixXd TransMatrix_Inv = Transform_Matrix.transpose().inverse();
        MSGUser()->MSG_INFO("Calculate New Variations");
        MSGUser()->MSG_INFO("Decide NVariations");
        MatrixXd m_Covariance_Matrix(Size, Size);
        for (int i = 0; i < Size; i++)
            for (int j = 0; j < Size; j++)
                m_Covariance_Matrix(i, j) = 0;
 
        for (int i = 0; i < Size; i++)
        {
            for (int j = 0; j < 2; j++)
            {
                double factor = 0;
                if (j == 0)
                    factor = 1;
                else
                    factor = -1;
 
                VectorXd RotatedNewVariation = New_Observable_Vector;
                RotatedNewVariation(Size - 1 - i) = New_Observable_Vector(Size - 1 - i) +
                                                    factor * sqrt(New_Covariance_Matrix(Size - 1 - i, Size - 1 - i));
                VectorXd newvar = TransMatrix_Inv * RotatedNewVariation;
                if (j == 0)
                    New_Observable_UpVariations.push_back(newvar);
                else
                    New_Observable_DownVariations.push_back(newvar);
            }
 
            for (int j = 0; j < Size; j++)
            {
                for (int k = 0; k < Size; k++)
                {
                    double diff_j = (New_Observable_UpVariations[i])(j)-Observable_Vector(j);
                    double diff_k = (New_Observable_UpVariations[i])(k)-Observable_Vector(k);
 
                    m_Covariance_Matrix(j, k) = m_Covariance_Matrix(j, k) + diff_j * diff_k;
                }
            }
 
            // Check Unc
            bool isOK = true;
            for (int j = 0; j < Size; j++)
            {
                double unc = sqrt(Covariance_Matrix(j, j));
                if (unc < 1e-8)
                    continue;
                double m_unc = sqrt(m_Covariance_Matrix(j, j));
 
                double factor = 1.0;
                int iunc = factor * unc;
                while (iunc < 10 || iunc >= 100)
                {
                    if (iunc < 10)
                    {
                        factor = factor * 10;
                    }
                    if (iunc >= 100)
                    {
                        factor = factor / 10;
                    }
                    iunc = factor * unc;
                }
 
                int imunc = factor * m_unc;
                if (imunc != iunc)
                {
                    isOK = false;
                    break;
                }
            }
 
            if (isOK)
            {
                NVariation = i + 1;
                break;
            }
        }
 
        if (NVariation < 0)
            NVariation = Size;
 
        isProcessed = true;
    }
 
    int PCATool::GetNPCAVariation()
    {
        if (!isProcessed)
            Process();
        return NVariation;
    }
 
    double PCATool::GetPCAResult(int ivar, int index, double nsigma)
    {
        auto titleguard = MSGUser()->StartTitleWithGuard("PCATool::GetPCAResult");
        if (ivar < 0 || ivar >= Size)
        {
            MSGUser()->MSG_ERROR("Variation Index not in range: ", ivar);
            return 0;
        }
        if (index < 0 || index >= Size)
        {
            MSGUser()->MSG_ERROR("Observable Index not in range: ", index);
            return 0;
        }
 
        if (!isProcessed)
            Process();
 
        if (ivar >= NVariation)
        {
            MSGUser()->MSG_ERROR("Variation Index larger than NVariation: ", ivar);
            return 0;
        }
 
        double diff;
        if (nsigma > 0)
        {
            diff = (New_Observable_UpVariations[ivar])(index)-Observable_Vector(index);
        }
        else
        {
            diff = (New_Observable_DownVariations[ivar])(index)-Observable_Vector(index);
        }
        return fabs(nsigma) * diff + Observable_Vector(index);
    }
 
    double PCATool::GetPCAResult(int index)
    {
        auto titleguard = MSGUser()->StartTitleWithGuard("PCATool::GetPCAResult");
        if (index < 0 || index >= Size)
        {
            MSGUser()->MSG_ERROR("Observable Index not in range: ", index);
            return 0;
        }
        if (!isProcessed)
            Process();
        return New_Observable_Vector(index);
    }
 
    double PCATool::GetPCAUpVariation(int ivar, int index)
    {
        return GetPCAResult(ivar, index, 1.0);
    }
 
    double PCATool::GetPCADownVariation(int ivar, int index)
    {
        return GetPCAResult(ivar, index, -1.0);
    }
 
    double PCATool::GetEigenValue(int ivar)
    {
        if (ivar < 0 || ivar >= Size)
        {
            MSGUser()->MSG_ERROR("Variation Index not in range: ", ivar);
            return 0;
        }
        if (!isProcessed)
            Process();
        return EigenValue(Size - 1 - ivar);
    }
} // namespace NAGASH