// $Id$

#ifndef _utl_CovarianceMatrix_h_
#define _utl_CovarianceMatrix_h_

#include <vector>
#include <cmath>


namespace utl {

  class CovarianceMatrix {
  public:
    CovarianceMatrix(const unsigned int n = 1)
    {
      SetExtent(n);
    }

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

    unsigned int
    GetExtent()
      const
    {
      return fExtent;
    }

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

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

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

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

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

    /// filling matrix with x in diagonal and zero everywhere else
    CovarianceMatrix&
    operator=(const 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:
    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