pca.h

Go to the documentation of this file.

#pragma once

/*+
________________________________________________________________________

 (C) dGB Beheer B.V.; (LICENSE) http://opendtect.org/OpendTect_license.txt
 Author:        Kristofer Tingdahl
 Date:          4-11-2002
________________________________________________________________________


-*/

#include "algomod.h"
#include "arrayndimpl.h"
#include "sets.h"
#include "trigonometry.h"
#include "threadwork.h"

template <class T> class Array2D;
class SequentialTask;
class PCACovarianceCalculator;

mExpClass(Algo) PCA
{
public:
                                        PCA( int nrvars );
    virtual                             ~PCA();

    void                                clearAllSamples();
    template <class IDXABL> void        addSample( const IDXABL& sample );
    bool                                calculate();
    float                               getEigenValue(int idx) const;
    template <class IDXABL> void        getEigenVector(int idx,
                                                       IDXABL& vec) const;
protected:
    bool                        tqli( float[], float[], int, ObjectSet<float>&);
    void                        tred2( ObjectSet<float>&, int, float[],float[]);

    const int                   nrvars_;
    Array2DImpl<float>          covariancematrix_;
    ObjectSet<TypeSet<float> >  samples_;
    TypeSet<float>              samplesums_;
    TypeSet<Threads::Work>      workload_;
    ObjectSet<SequentialTask>   tasks_;
    float*                      eigenvalues_;
    int*                        eigenvecindexes_;
};


template <class IDXABL> inline
void PCA::addSample( const IDXABL& sample )
{
    TypeSet<float>& ownsample = *new TypeSet<float>;
    for ( int idx=0; idx<nrvars_; idx++ )
    {
        const float val = (float) sample[idx];
        ownsample += val;
        samplesums_[idx] += val;
    }

    samples_ += &ownsample;
}

template <class IDXABL> inline
void PCA::getEigenVector(int idy, IDXABL& vec ) const
{
    for ( int idx=0; idx<nrvars_; idx++ )
        vec[idx] = covariancematrix_.get(idx, eigenvecindexes_[idy] );
}