///////////////////////////////////////////////////////////////////////////
// ASVR: Active Support Vector Regression                                //
// Copyright 2002 David R. Musicant and Alexander Feinberg.              //
// This software is free for academic use only.                          //
// For commercial use, contact dmusican@carleton.edu.                    //
///////////////////////////////////////////////////////////////////////////
// DataSource.H: Header information for the objects dealing with data    //
// storage.                                                              //
///////////////////////////////////////////////////////////////////////////

#ifndef DATASOURCE_H
#define DATASOURCE_H

#include <stdio.h>
#include <iostream.h>
#include <fstream.h>
#include <assert.h>
#include <stdlib.h>
#include <set>
#include "math.h"
#include "f2c.h"

// Visual C++ is finicky about wanting this for some things
using namespace std ;

// Indicate CLAPACK routines as C external so as to avoid name
// mangling problems
extern "C" {
  int dgetrf_(int *m, int *n, double *A, int *lda, int *ipiv, int*info);
  int dgetri_(int *n, double *A, int *lda, int *ipiv, double *work,
        int *lwork, int*info);
}

#define MAXFILENAMELENGTH 500 // Maximum filename length
#define MAXROWLENGTH 5000   // Maximum text length of a row in data file
#define TRUE 1
#define FALSE 0

// An inheritance type structure is used here to allow future development in:
// - optimizations in sparse matrices
// - data which is resident on disk (and not in memory)

// A template structure is used to allow the class to function for
// storage of doubles, integers, shorts, etc.

class Vector;
class DenseVector;

// Abstract Matrix class
class MatrixType {
 public:
  FILE *filePtr;
  int rows, cols;
  MatrixType(int rows, int cols);
  virtual ~MatrixType(){};
  virtual double Show(int row, int col)=0;
  virtual void Display()=0;
  virtual void DiagAddInPlace(double addOn)=0;
  virtual void MultiplyInPlace(double factor)=0;
  virtual void GaussianInverse()=0;
  virtual void AddInPlace(MatrixType *m)=0;
  virtual void Predict(Vector *w, double gamma, DenseVector *results)=0;
  virtual double Accuracy(Vector *w, double gamma, Vector *y)=0;
  virtual void GetMeanStd(double *mean, double *std)=0;
  virtual void Normalize(double *mean, double *std)=0;
  virtual double Sum()=0;
  virtual double Square()=0;
  virtual double Norm1()=0;
  virtual double NormT1()=0;
};

// Abstract matrix class for a specific data type
template<class T> class MatrixTemplate : public MatrixType {
 public:
  MatrixTemplate(int rows, int cols) : MatrixType(rows,cols){};
  virtual ~MatrixTemplate(){};
};

// Concrete matrix class for dense data
template<class T> class DenseMatrixType : public MatrixTemplate<T> {
 public:
  T **dataMatrix;
  DenseMatrixType(int rows, int cols);
  ~DenseMatrixType();
  double Show(int row, int col);
  void Display();
  void DiagAddInPlace(double addOn);
  void MultiplyInPlace(double factor);
  void GaussianInverse();
  void AddInPlace(MatrixType *m);
  void Predict(Vector *w, double gamma, DenseVector *results);
  double Accuracy(Vector *w, double gamma, Vector *y);
  void GetMeanStd(double *mean, double *std);
  void Normalize(double *mean, double *std);
  double Sum();
  double Square();
  double Norm1();
  double NormT1();
};

// Abstract vector class
class Vector{
 public:
  int cols;
  Vector(int cols);
  virtual ~Vector(){};
  virtual double operator[](int col)=0;
  virtual bool operator>=(double value)=0;
  virtual bool operator<=(double value)=0;
  virtual int Size()=0;
  virtual void Assign(int col,double value)=0;
  virtual void Display()=0;
  virtual void Display(int m)=0;
  virtual set<int> findTolNonBasic(double tol)=0;
  virtual set<int> findTolBasic(double tol)=0;
  virtual void SubtractInPlace(Vector *m2)=0;
  virtual void AddInPlace(Vector *m2)=0;
  virtual void AddInPlace(double value)=0;
  virtual void MultiplyInPlace(double factor)=0;
  virtual void PlusFunctionInPlace()=0;
  virtual void Copy(Vector *m)=0;
  virtual void SubtractInPlaceFrom(double value)=0;
  virtual void Update(Vector *m, set<int> indices)=0;
  virtual void Zero(set<int> indices)=0;
  virtual void Strip(Vector *m, set<int> indices)=0;
  virtual double Square()=0;
  virtual double Sum()=0;
  virtual double ScalarProduct(Vector *m2) = 0;

  virtual set<int> findTolNonBasic(double tol, set<int> indices)=0;
  virtual double findMin(set<int> &indices)=0;
  virtual double ScalarProduct(Vector *m2, set<int> &b) = 0;
  virtual void SubtractInPlaceFrom(Vector *vec)=0;
  virtual set<int> findTolNegative(double tol, set<int> &indices)=0;
  virtual set<int> findTolNegative(double tol)=0;
  virtual void ZeroInactivesInPlace(const set<int> &b)=0;
  virtual void MultiplyThenAddInPlace(double factor, Vector *addMe, const set<int> &b)=0;
  virtual void MultiplyThenAddInPlace(double factor, Vector *addMe)=0;
};

// Concrete vector class for dense data
class DenseVector : public Vector {
 public:
  double *dataMatrix;
  DenseVector(int cols);
  ~DenseVector();
  double operator[](int col);
  bool operator>=(double value);
  bool operator<=(double value);
  void Assign(int col,double value);
  int Size();
  void Display();
  void Display(int m);
  set<int> findTolNonBasic(double tol);
  set<int> findTolBasic(double tol);
  void SubtractInPlace(Vector *m2);
  void AddInPlace(Vector *m2);
  void AddInPlace(double value);
  void MultiplyInPlace(double factor);
  void PlusFunctionInPlace();
  void Copy(Vector *m);
  void SubtractInPlaceFrom(double value);
  void Update(Vector *m, set<int> indices);
  void Zero(set<int> indices);
  void Strip(Vector *m, set<int> indices);
  double Square();
  double Sum();
  double ScalarProduct(Vector *m2);

  set<int> findTolNonBasic(double tol, set<int> indices);
  double findMin(set<int> &indices);
  double ScalarProduct(Vector *m2, set<int> &b);
  void SubtractInPlaceFrom(Vector *vec);
  set<int> findTolNegative(double tol, set<int> &indices);
  set<int> findTolNegative(double tol);
  void ZeroInactivesInPlace(const set<int> &b);
  void MultiplyThenAddInPlace(double factor, Vector *addMe, const set<int> &b);
  void MultiplyThenAddInPlace(double factor, Vector *addMe);
};

// Matrix and vector binary operations
template<class T1, class T2, class T3>
void MultAB(MatrixTemplate<T1> *m1, MatrixTemplate<T2> *m2,
	    MatrixTemplate<T3> *output);
template<class T1, class T2, class T3>
void MultAB(MatrixTemplate<T1> *m1, set<int> i1, 
	    MatrixTemplate<T2> *m2, MatrixTemplate<T3> *output);
template<class T1>
void MultAB(MatrixTemplate<T1> *m1, set<int> i1, Vector *m2, 
	    Vector *output);
template<class T1>
void MultAB(MatrixTemplate<T1> *m1, Vector *m2, Vector *output);
template<class T1, class T2, class T3>
void MultAtB(MatrixTemplate<T1> *m1, MatrixTemplate<T2> *m2,
	     MatrixTemplate<T3> *output);
template<class T1>
void MultAtB(MatrixTemplate<T1> *m1, Vector *m2, Vector *output);
template<class T1, class T2, class T3>
void MultAtB(MatrixTemplate<T1> *m1, set<int> i1,
	     MatrixTemplate<T2> *m2,
	     set<int> i2, MatrixTemplate<T3> *output);
template<class T1>
void MultAtB(MatrixTemplate<T1>*m1, set<int> i1, Vector *m2, 
	     set<int> i2, Vector *output);
template<class T1>
void MultAtB(MatrixTemplate<T1>*m1, set<int> set1, Vector *m2, 
	     Vector *output);
template<class T1, class T2, class T3>
void MultAtB(MatrixTemplate<T1> *m1, set<int> i1,
	     MatrixTemplate<T2> *m2, MatrixTemplate<T3> *output);
template<class T1, class T2, class T3>
void MultABt(MatrixTemplate<T1> *m1, MatrixTemplate<T2> *m2,
	     set<int> i2, MatrixTemplate<T3> *output);

// Miscellaneous functions
void ThrowIt(const char *text);



// Functions added for regression.
template <class T1, class T2>
void AppendConstCol(MatrixTemplate<T1> *in, double val, MatrixTemplate<T2> *out);

void ElDivAB(Vector *vIn1, Vector *vIn2, set<int> indices, Vector *vOut);

template <class T1>
void MultAtB(Vector *in1, MatrixTemplate<T1> *in2, Vector *out);

template<class T1, class T2, class T3>
void MultABt(MatrixTemplate<T1> *A, MatrixTemplate<T2> *B, MatrixTemplate<T3> *output);

template <class T1>
void MultAB(Vector *in1, set<int> &b, MatrixTemplate<T1> *in2, Vector *out);

template <class T1>
void MultAB(Vector *in1, MatrixTemplate<T1> *in2, Vector *out);

void SubAB(Vector *A, Vector *B, const set<int> &b, Vector *C);

#endif
