// 
//  MinVariance.h
//  Qdos2 - Minimum Variance Algorithm 
//
//  Janet Barnes
//  May 2000
//

#ifndef QDOS_MINVARIANCE_H
#define QDOS_MINVARIANCE_H

#include "Algorithm.h"
#include "Matrix.cc"

namespace QDOS{

class MinVariance
{
    Matrix<double> _eigvals;       // Matrix Dim(3) of eigenvalues
    Matrix<double> _eigvecs;       // Matrix Dim(3,3) of eigenvectors (in rows)
    Matrix<double> _eigqual;       // Matrix Dim(3) eigenvalue quality
    double _orgqual;               // orthogonality quality
    Matrix<double> _covariance;    // covariance matrix
    
public:
    MinVariance(const MatrixSeq<double>& ms):_eigvals(Matrix<double>(Dim(3))),\
        _eigvecs(Matrix<double>(Dim(3,3))),_covariance(Matrix<double>(Dim(3,3))),\
        _eigqual(Matrix<double>(Dim(3)))
    {
        if(ms.empty())
            return;
        if (!ms.is_homogeneous_seq())
            throw MatrixDomainError();
        if(ms.elt(0).rank()!=1)
            throw MatrixRankError();
        if(ms.elt(0).dimdata()!=Dim(3))
             throw MatrixLengthError();
             
        double sz=ms.size();
        // compute the variance matrix 
for (int k=0; k<sz; k++) cout << ms[k] << endl;       
        ScalarSeq<double> x=ms.get_element_seq(vector<int>((size_t)1,0));
        ScalarSeq<double> y=ms.get_element_seq(vector<int>((size_t)1,1));
        ScalarSeq<double> z=ms.get_element_seq(vector<int>((size_t)1,2));

        double sum_x=x.sum();
        double sum_y=y.sum();
        double sum_z=z.sum();

        _covariance.resize(Dim(3,3));
   
        _covariance[0][0]=((x*x).sum() - sum_x*sum_x/sz)/sz;
        _covariance[0][1]=((x*y).sum() - sum_x*sum_y/sz)/sz;
        _covariance[0][2]=((x*z).sum() - sum_x*sum_z/sz)/sz;
        _covariance[1][1]=((y*y).sum() - sum_y*sum_y/sz)/sz;
        _covariance[1][2]=((y*z).sum() - sum_y*sum_z/sz)/sz;
        _covariance[2][2]=((z*z).sum() - sum_z*sum_z/sz)/sz;
   
        _covariance[1][0] = _covariance[0][1];
        _covariance[2][0] = _covariance[0][2];
        _covariance[2][1] = _covariance[1][2];

        // get eigenvectors and values  
        MatrixEig eig(_covariance);
        _eigvals=eig.eigenvalues().seq_to_matrix();
        _eigvecs=eig.eigenvectors().seq_to_matrix();
		
		// check the positive sense of the eigenvector associated with
        // the minimum eigenvalue and if necessary reverse the signs
        // of all three components i.e. vint x vmax = vector perpendicular
        // to both vint and vmax. If vmin is perpendicular to vint 
        // and vmax (i.e. system is ortogonal), then (vinx x vmax).vmin=1
        // If < 0 then not a right-handed system, so reverse signs of vmin
        double rh_factor=(_eigvecs[1].cross(_eigvecs[2])).dot(_eigvecs[0]);
        if(rh_factor<0.)
            _eigvecs[0]*=-1;
        _orgqual=std::abs(std::abs(rh_factor)-1.0);
                
        // ensure that the sign of the average minimum variance component
        // is positive by computing average of sum of B-field data.min
        // eigenvector. If less than zero negate both vmin and vint (so
        // preserving right-handed system)
        double avg_minvar_comp=(ms*_eigvecs[0]).sum();
        if(avg_minvar_comp<0) 
        {
            _eigvecs[0]*=-1;
            _eigvecs[1]*=-1;
        }
        _eigqual=_eigvals/_eigvals[2];
            
    }
    
    // methods for returning members and calculating variance
    Matrix<double>& covariance() 
    { return _covariance; }
    
    Matrix<double>& eigenvalues() 
    { return _eigvals;}
     
    Matrix<double>& eigenvectors()
    { return _eigvecs;}   
    
    Matrix<double>& eigenquality() 
    { return _eigqual;}
    
    double& orthogonality_quality()
    { return _orgqual;}    // orthogonality quality     
        
};
    
}  // end namespace QDOS
#endif

