// TrfVector_LA.h

#ifndef TrfVector_LA_H
#define TrfVector_LA_H

// Version of TrfVector using FNAL ZOOM linearAlgebra package.

// Contains headers for underlying TrfVector, TrfMatrix and TrfSMatrix
// classes and corresponding typedefs.

#include <cassert>
#include <cmath>
#include <vector>
#include "LinearAlgebra/ColumnVector.h"
#include "LinearAlgebra/Matrix.h"
#include "trfutil/TRFMath.h"

//**********************************************************************

class TrfVector : public ColumnVector {

public:  // methods

  // Constructor.
  explicit TrfVector(int len) : ColumnVector(len) { }

  // Constructor from subclass.
  TrfVector(const ColumnVector& vec)
  : ColumnVector(vec) { }

  // Constructor from subclass.
  TrfVector(const MatrixD& mtx) : ColumnVector(mtx.rows()) {
    assert( mtx.columns() == 1 );
    MatrixD::operator=(mtx);
  }

  // Return the length.
  int length() const { return rows(); }

  // Return the length.
  int nrow() const { return rows(); }

  // Fill the vector.
  void fill(double val) {
    for ( int i=0; i<rows(); ++i ) operator()(i) = val;
  }

  // Return the minimum value in the vector.
  double min() const {
    double val = operator()(0);
    for ( int i=1; i<rows(); ++i ) {
      double newval = operator()(i);
      if ( newval < val ) val = newval;
    }
    return val;
  }
    
  // Return the maximum value in the vector.
  double max() const {
    double val = operator()(0);
    for ( int i=1; i<rows(); ++i ) {
      double newval = operator()(i);
      if ( newval > val ) val = newval;
    }
    return val;
  }
    
  // Return the minimum absolute value in the vector.
  double amin() const {
    double val = std::abs(operator()(0));
    for ( int i=1; i<rows(); ++i ) {
      double newval = std::abs(operator()(i));
      if ( newval < val ) val = newval;
    }
    return val;
  }
    
  // Return the minimum absolute value in the vector.
  double amax() const {
    double val = std::abs(operator()(0));
    for ( int i=1; i<rows(); ++i ) {
      double newval = std::abs(operator()(i));
      if ( newval > val ) val = newval;
    }
    return val;
  }
    
  // Addition.
  TrfVector& operator+=(const TrfVector& rhs) {
    TrfVector& lhs = *this;
    assert( lhs.nrow() == rhs.nrow() );
    for ( int i=0; i<lhs.nrow(); ++i ) {
      lhs(i) += rhs(i);
    }
    return lhs;
  }

  // Subtraction.
  TrfVector& operator-=(const TrfVector& rhs) {
    TrfVector& lhs = *this;
    assert( lhs.nrow() == rhs.nrow() );
    for ( int i=0; i<lhs.nrow(); ++i ) {
      lhs(i) -= rhs(i);
    }
    return lhs;
  }

};

// Addition.
inline
TrfVector operator+(const TrfVector& rhs1, const TrfVector& rhs2) {
  TrfVector lhs = rhs1;
  lhs += rhs2;
  return lhs;
}

// Subtraction.
inline
TrfVector operator-(const TrfVector& rhs1, const TrfVector& rhs2) {
  TrfVector lhs = rhs1;
  lhs -= rhs2;
  return lhs;
}

// Equality.
inline
bool operator==(const TrfVector lhs, const TrfVector& rhs) {
  if ( lhs.nrow() != rhs.nrow() ) {
    return false;
  }
  for ( int i=0; i<lhs.nrow(); ++i ) {
    if ( lhs(i) != rhs(i) ) {
      return false;
    }
  }
  return true;
}

// Near equality.
inline
bool is_equal(const TrfVector lhs, const TrfVector& rhs) {
  if ( lhs.nrow() != rhs.nrow() ) {
    return false;
  }
  for ( int i=0; i<lhs.nrow(); ++i ) {
    if ( ! is_equal(lhs(i), rhs(i)) ) {
      return false;
    }
  }
  return true;
}

//**********************************************************************

class TrfMatrix : public MatrixD {

public:  // methods

  // constructor
  explicit TrfMatrix(int ncol, int nrow)
  : MatrixD(ncol,nrow) { }

  // Constructor from subclass.
  TrfMatrix(const MatrixD& mtx)
  : MatrixD(mtx) { }

  // Return the length.
  int nrow() const { return rows(); }
  int ncol() const { return columns(); }

  // Fill the matrix.
  void fill(double val) {
    for ( int i=0; i<rows(); ++i ) {
      for ( int j=0; j<columns(); ++j ) {
        operator()(i,j) = val;
      }
    }
  }

  // Return the minimum value in the matrix.
  double min() const {
    double val = operator()(0,0);
    for ( int i=0; i<rows(); ++i ) {
      for ( int j=0; j<columns(); ++j ) {
        double newval = operator()(i,j);
        if ( newval < val ) val = newval;
      }
    }
    return val;
  }

  // Return the maximum value in the matrix.
  double max() const {
    double val = operator()(0,0);
    for ( int i=0; i<rows(); ++i ) {
      for ( int j=0; j<columns(); ++j ) {
        double newval = operator()(i,j);
        if ( newval > val ) val = newval;
      }
    }
    return val;
  }

  // Return the minimum absolute value in the matrix.
  double amin() const {
    double val = std::abs(operator()(0,0));
    for ( int i=0; i<rows(); ++i ) {
      for ( int j=0; j<columns(); ++j ) {
        double newval = std::abs(operator()(i,j));
        if ( newval < val ) val = newval;
      }
    }
    return val;
  }

  // Return the maximum absolute value in the matrix.
  double amax() const {
    double val = std::abs(operator()(0,0));
    for ( int i=0; i<rows(); ++i ) {
      for ( int j=0; j<columns(); ++j ) {
        double newval = std::abs(operator()(i,j));
        if ( newval > val ) val = newval;
      }
    }
    return val;
  }

  // Transpose.
  TrfMatrix transpose() const {
    return MatrixD::transpose();
  }

};

//**********************************************************************

class TrfSMatrix : public MatrixD {

public:  // methods

  // constructor
  explicit TrfSMatrix(int dim) : MatrixD(dim,dim) {
    declareSymmetric();
  }

  // Constructor from subclass.
  TrfSMatrix(const MatrixD& mtx)
  : MatrixD(mtx) { }

  // Return the length.
  int nrow() const { return rows(); }
  int ncol() const { return columns(); }

  // Fill the matrix.
  void fill(double val) {
    for ( int i=0; i<rows(); ++i ) {
      for ( int j=0; j<=i; ++j ) {
        operator()(i,j) = val;
      }
    }
  }

  // Return the minimum value in the matrix.
  double min() const {
    double val = operator()(0,0);
    for ( int i=1; i<rows(); ++i ) {
      for ( int j=0; j<=i; ++j ) {
        double newval = operator()(i,j);
        if ( newval < val ) val = newval;
      }
    }
    return val;
  }

  // Return the maximum value in the matrix.
  double max() const {
    double val = operator()(0,0);
    for ( int i=1; i<rows(); ++i ) {
      for ( int j=0; j<i; ++j ) {
        double newval = operator()(i,j);
        if ( newval > val ) val = newval;
      }
    }
    return val;
  }
    
  // Return the minimum absolute value in the matrix.
  double amin() const {
    double val = std::abs(operator()(0,0));
    for ( int i=1; i<rows(); ++i ) {
      for ( int j=0; j<=i; ++j ) {
        double newval = std::abs(operator()(i,j));
        if ( newval > val ) val = newval;
      }
    }
    return val;
  }

  // Return the maximum absolute value in the matrix.
  double amax() const {
    double val = std::abs(operator()(0,0));
    for ( int i=1; i<rows(); ++i ) {
      for ( int j=0; j<=i; ++j ) {
        double newval = std::abs(operator()(i,j));
        if ( newval < val ) val = newval;
      }
    }
    return val;
  }

  // Xform  S --> M * S * MT
  TrfSMatrix xform(const TrfMatrix& mtx) const {
    return xsxt(mtx);
  }

  // Xform  S --> MT * S * M
  TrfSMatrix xformt(const TrfMatrix& mtx) const {
    return xtsx(mtx);
  }

};

// Normalize a matrix.
// Transformation which makes all diagonal elements unity.
inline
void normalize(TrfSMatrix& sma) {
  std::vector<double> vec(sma.nrow());
  for ( int irow=0; irow<sma.nrow(); irow++ ) {
    vec[irow] = 1.0/std::sqrt( sma(irow,irow) );
    double rowfac = sqrt( vec[irow] );
    for ( int icol=0; icol<=irow; icol++ ) {
      sma(irow,icol) *= vec[irow]*vec[icol];
    }
  }
}

//**********************************************************************

// Chi-square difference.
inline
double chisq_diff(const TrfVector& vecdiff, const TrfSMatrix& inverr) {
  MatrixD mtx = inverr.xtsx(vecdiff);
  assert( mtx.rows() == 1 );
  assert( mtx.columns() == 1 );
  return mtx(0,0);
}

// Matrix Equality.
inline
bool operator==(const TrfMatrix& lhs, const TrfMatrix& rhs) {
  if ( lhs.nrow() != rhs.nrow() ) {
    return false;
  }
  if ( lhs.ncol() != rhs.ncol() ) {
    return false;
  }
  for ( int i=0; i<lhs.nrow(); ++i ) {
    for ( int j=0; j<=lhs.ncol(); ++j ) {
      if ( lhs(i,j) != rhs(i,j) ) {
        return false;
      }
    }
  }
  return true;
}

// Symmetric matrix Equality.
inline
bool operator==(const TrfSMatrix& lhs, const TrfSMatrix& rhs) {
  if ( lhs.nrow() != rhs.nrow() ) {
    return false;
  }
  for ( int i=0; i<lhs.nrow(); ++i ) {
    for ( int j=0; j<=i; ++j ) {
      if ( lhs(i,j) != rhs(i,j) ) {
        return false;
      }
    }
  }
  return true;
}

// Matrix near equality.
inline
bool is_equal(const TrfMatrix& lhs, const TrfMatrix& rhs) {
  if ( lhs.nrow() != rhs.nrow() ) {
    return false;
  }
  if ( lhs.ncol() != rhs.ncol() ) {
    return false;
  }
  for ( int i=0; i<lhs.nrow(); ++i ) {
    for ( int j=0; j<=lhs.ncol(); ++j ) {
      if ( ! is_equal(lhs(i,j), rhs(i,j)) ) {
        return false;
      }
    }
  }
  return true;
}

// Near equality.
inline
bool is_equal(const TrfSMatrix& lhs, const TrfSMatrix& rhs) {
  if ( lhs.nrow() != rhs.nrow() ) {
    return false;
  }
  for ( int i=0; i<lhs.nrow(); ++i ) {
    for ( int j=0; j<=i; ++j ) {
      if ( ! is_equal(lhs(i,j), rhs(i,j)) ) {
        return false;
      }
    }
  }
  return true;
}

// Invert.
// Return 0 for success.
inline
int invert(TrfSMatrix& mtx) {
  mtx = mtx.inverse();
  return 0;
}

//**********************************************************************

#endif