#ifndef __MATRIX_H__
#define __MATRIX_H__

#ifndef PI
#define PI 3.1415
#endif

#include <math.h>
#include <iostream>
#include <stdlib.h>
#include <time.h>
#include <vector>
#include <string.h>
#include <pthread.h>
#include "timer.h"

///class for a R2 matrix
/**
 * has some useful overloaded function for matrix
 * multiply and addition.  The number of overloads 
 * need to be expanded to include more cases.
 */
class matrix
{
 public:
  matrix();
  matrix( const matrix &o );
  ///constructore to allocate 
  matrix(int r, int c);
  ~matrix();
  int del();
  ///initialize a matrix from double array
  int init(double **md, int r, int c);
  ///initialize a matrix from array reshape to size
  int init(double *md, int r, int c);
  ///initialize a column matrix
  int init(double *md, int r);
  ///initialize a blank matrix
  int init(int r, int c);
  /// initialize a matrix fill with some data
  int initrand(int r, int c);
  ///initialize a identity matrix with v down the column
  /**
   * @param v value to place on diagnal
   * @param s size of matrix
   */
  int initidentity(double v, int s);

  double maxdiag();

  ///zeros out the matrix
  int zero();

  int notgood();

  int serialize( char **buf, int &size );
  int deserialize( char *buf );

  matrix &operator=(const matrix &);
  const matrix operator*(const matrix &);
  const matrix operator*(const double);
  const matrix operator+(const matrix &);
  const matrix operator-(const matrix &);
  const matrix operator/(const double );

  /// multiply a vector of matrices
  /**
   * results of the multiplication is stored in the matrix
   */
  int multiply( std::vector< matrix *> &ml );

  ///transpose a matrix
  /**
   * @return a copy of the transposed matrix
   */
  matrix transpose();
  int transpose( matrix & ) const;

  ///invert a matrix
  /**
   * @return inverted matrix
   * @warning matrix invertion does not check for invertable matrix
   */
  matrix invert();

  matrix invertU();
  matrix invertL();

  /// returns a matrix of same size with only diagnal values
  int diag(matrix &c);

  ///set a value
  /**
   * @param v value
   * @param r row #
   * @param c column #
   */
  int set(double v, int r, int c);
  ///gets a value
  /**
   * @param r row #
   * @param c column #
   */
  double get(int r, int c) const;
  
  ///return the number of rows
  int rows() const;
  ///return the number of columns
  int cols() const;
  ///craete a minor matrix used to find determinant
  /**
   * @param m original matrix
   * @param i row to skip
   * @param j col to skip
   */
  int initminor(matrix &m, int i, int j);
  ///return the determinate of the matrix
  double determinant();
  ///checks if matrix is a identity matrix
  int isidentity();
  ///calculate the Cholesky factorization of the matrix
  /**
   * implemenataion as seen in Numberical Recipes
   * @param c place to store the matrix
   * @return 0 on success
   */
  int cholesky(matrix &c);
  ///returns the sum of the diagonal
  double diagsum();
  ///returns the sum of the log diagonal
  double logdiagsum();

  int add(const matrix &m);
  int divide(double m);
  int multiply(double m);
  int addtranspose(matrix &m);

  int addmult(const matrix &m,
              const double c);
  int addtransposemult(const matrix &m,
                       const double c);

  /// this = this*b'
  int transposemult( const matrix &m );

  /// transpose this = a*b'
  int transposemult( const matrix &a, const matrix &b );

  /// sum of all values in matrix
  double sum();

  double *dataptr()
    {
      return d;
    }

 private:
  /// matrix multiply
  /**
   * a*b -> c
   * use a buffer if provided
   * create buffer if it does not exist or is too small
   */
  int multiply( matrix *a, matrix *b,
                matrix *c,
                double **buffer,
                int &bufsize );

  ///matrix pivot function used to invert a matrix
  int pivot(double **d, int pi, int pj,
            int sr, int sc);
  int col;
  int row;
  int detflag;
  double cdeterminant;
  double *d;

  //threads locking stuff
#ifdef __THREADING__
  pthread_mutex_t lock;   ///
#endif
  int transpose_nolock( matrix &c ) const;
  ///private functions that do not lock object
  int init_nolock(double **md, int r, int c);
  ///initialize a matrix from array reshape to size
  int init_nolock(double *md, int r, int c);
  ///initialize a column matrix
  int init_nolock(double *md, int r);
  ///initialize a blank matrix
  int init_nolock(int r, int c);
  double determinant_nolock();
  int initminor_nolock(matrix &m, int i, int j);
  int zero_nolock(); 
};

std::ostream& operator <<(std::ostream &o, const matrix &m);

double point_dis(matrix &p1, matrix &p2 );

#endif
