// $Id: CovarianceMatrix.h 18212 2010-12-21 14:16:51Z darko $

#ifndef _utl_CovarianceMatrix_h_
#define _utl_CovarianceMatrix_h_

#include <vector>
#include <cmath>

namespace utl {

  class CovarianceMatrix
  {
  public:

    CovarianceMatrix(unsigned int n=1)
    {
      SetExtent(n);
    }

    inline
    void
    SetExtent(unsigned int n)
    {
      fExtent = n;
      fData.resize((n*(n+1))/2);
    }

    inline
    unsigned int
    GetExtent()
      const
    {
      return fExtent;
    }

    inline
    double&
    operator()(unsigned int i, unsigned int j)
    {
      return fData[Internal(i,j)];
    }

    inline
    const double&
    operator()(unsigned int i, unsigned int j)
      const
    {
      return fData[Internal(i,j)];
    }

    inline
    double&
    operator[](unsigned int i)
    {
      return fData[Internal(i,i)];
    }

    inline
    const double&
    operator[](unsigned int i)
      const
    {
      return fData[Internal(i,i)];
    }

    /// convenience: return std. deviation of parameter i
    inline
    double
    Std(unsigned int i)
      const
    {
      return std::sqrt(fData[Internal(i,i)]);
    }

    /// filling matrix with x in diagonal and zero everywhere else
    inline
    CovarianceMatrix&
    operator=(double x)
    {
      for (unsigned int i=0;i<fExtent;++i)
        for (unsigned int j=i;j<fExtent;++j)
          fData[Internal(i,j)] = i==j ? x : 0;
      return *this;
    }

  protected:

    inline
    unsigned int
    Internal(unsigned int i, unsigned int j)
      const
    {
      if (i>j) std::swap(i,j);
      return j + i*fExtent - (i*(i+1))/2;
    }
  
    unsigned int fExtent;
    std::vector<double> fData;
  };
  
} // NS utl


#endif