src/include/utils.h

/*
 * utils.h
 *
 *  Created on:
 *      Author: romain
 */

#ifndef UTILS_H_
#define UTILS_H_

#include "iostream"
#include "string"
#include "fstream"
#include "iomanip"
#include "stdio.h"
#include "vector"
#include "math.h"
#include "stdlib.h"
#include "assert.h"
#include <gmpxx.h>
#include <armadillo>


namespace nltool {
#define PI 3.14159265

#ifdef GMP_RAT
typedef mpq_class Number;
#else
typedef double Number;
#endif

const Number EPSILON=1e-10;

typedef arma::Mat<Number> Matrix;
typedef arma::Col<Number> Vector;
typedef std::vector<Vector> Vertices;
typedef std::vector<Vector> Template;
typedef std::vector<double> doubleVector;
typedef std::vector<Number> stdVector;
typedef std::vector<std::vector<Number> > stdMatrix;

extern bool useBinomTable;
extern int** binomialCoefficientTable;



//class Simplex:public Set{
//public:
//      Simplex(int dim, stdMatrix& vert):dimension(dim),verticesData(vert){}
//      Simplex():dimension(0),verticesData(stdMatrix(0)){};
//      int dim(){return dimension;}
//      stdMatrix vertices(){
//               return getVertices();
//      }
//      stdMatrix getVertices(){return verticesData;}
//      unsigned int nb_vertices(){return verticesData.size();}
//
//private:
//      stdMatrix verticesData;
//      int dimension;
//};

// ###############################
// Matrix and vector operations
// ###############################

inline arma::Mat<Number> template2Matrix(const Template& temp) {
        arma::Mat<Number> result(temp.size(), (temp[0].n_elem));
        for (unsigned int i = 0; i < temp.size(); ++i) {
                for (unsigned int j = 0; j < temp[0].n_elem; ++j) {
                        result(i, j) = temp[i][j];
                }
        }
        return result;
}

inline Template matrixToTemplate(const Matrix& mat) {
        Template result(mat.n_rows, Vector(mat.n_cols));
        for (unsigned int i = 0; i < mat.n_rows; ++i) {
                for (unsigned int j = 0; j < mat.n_cols; ++j) {
                        result[i][j] = mat(i,j);
                }
        }
        return result;
}

inline Number infinityNorm(stdVector v) {
        Number res = fabs(v[0]);
        for(unsigned i=1; i<v.size();i++)
                if (res<fabs(v[i]))
                        res = fabs(v[i]);
        return res;
}

stdMatrix matrixTranspose(const stdMatrix m1);

inline stdVector operator*(stdMatrix A, const stdVector v) {
        assert((A.size()>0) && (A[0].size()==v.size()));
        stdVector result(A.size(), 0);
        for (unsigned i = 0; i < A.size(); i++)
                for (unsigned j = 0; j < v.size(); j++)
                        result[i] += A[i][j] * v[j];
        return result;
}

inline stdVector operator*(Number n, const stdVector v) {
        stdVector result(v.size());
        for (unsigned i = 0; i < v.size(); i++)
                result[i] = n * v[i];
        return result;
}

inline stdMatrix operator*(Number n, const stdMatrix m) {
        stdMatrix result(m.size());
        for (unsigned i = 0; i < m.size(); i++)
                result[i] = n * m[i];
        return result;
}

inline stdVector operator+(stdVector v1, const stdVector v2) {
        assert(v1.size()==v2.size());
        stdVector result(v1.size());
        for (unsigned i = 0; i < v1.size(); i++)
                result[i] = v1[i] + v2[i];
        return result;
}

inline stdVector operator-(stdVector v1, const stdVector v2) {
        assert(v1.size()==v2.size());
        stdVector result(v1.size());
        for (unsigned i = 0; i < v1.size(); i++)
                result[i] = v1[i] - v2[i];
        return result;
}

inline stdMatrix operator+(stdMatrix m1, const stdMatrix m2) {
        assert(m1.size()==m2.size());
        stdMatrix result(m1.size());
        for (unsigned i = 0; i < m1.size(); i++) {
                assert(m1[i].size()==m2[i].size());
                result[i] = m1[i] + m2[i];
        }
        return result;
}

inline bool operator<=(stdVector v1, const stdVector v2) {
        assert(v1.size()==v2.size());
        for (unsigned i = 0; i < v1.size(); i++)
                if (v1[i] > v2[i])
                        return false;
        return true;
}

inline Number dot(stdVector v1, stdVector v2) {
        Number result = 0.0;
        for (unsigned int i = 0; i < v1.size(); ++i) {
                result += v1[i] * v2[i];
        }
        return result;
}

inline bool notNul(stdVector v) {
        for (unsigned int i = 0; i < v.size(); i++)
                if (v[i] != 0)
                        return true;
        return false;
}

inline int fact(const int n) {
        int result = 1;
        for (int i = 2; i <= n; i++)
                result *= i;
        return result;
}

inline Number max(const stdVector& v){
        assert(v.size()>0);
        Number maximum=v[0];
        for(unsigned i=1; i<v.size(); i++)
                if(maximum<v[i])
                        maximum = v[i];
        return maximum;
}

inline Number min(const stdVector& v){
        assert(v.size()>0);
        Number maximum=v[0];
        for(unsigned i=1; i<v.size(); i++)
                if(maximum>v[i])
                        maximum = v[i];
        return maximum;
}

inline bool closeToZero(Number nb){
        return (fabs(nb)<EPSILON);
}


double binomialCoef(const int n, const int k);

void createBinomialCoefficientTable(const int max);

inline std::ostream& operator<<(std::ostream& out, stdVector v) {
        for (unsigned int i = 0; i < v.size(); i++)
                out << v[i] << " ";
        out << std::endl;
        return out;
}

inline std::ostream& operator<<(std::ostream& out, stdMatrix t) {
        for (unsigned int i = 0; i < t.size(); i++)
                out << t[i];
        return out;
}

inline std::ostream& operator<<(std::ostream& out, Template t) {
        for (unsigned int i = 0; i < t.size(); i++)
                out << t[i].t()<<std::endl;
        return out;
}


stdVector setUnitVector(const Vector v);

Matrix matrixExponential(Matrix A, Number dt);



// ###############################
// LP solver
// ###############################

Number simplex(Vector obj, Template M, Vector N, Vector& solution, bool max = true);

Number simplex(Vector obj, Template M, Vector N, bool max = true);

Number simplex(stdVector obj, stdMatrix Aineq, stdVector Bineq, stdMatrix Aeq,
                stdVector Beq, stdVector bl, stdVector bu, int &err);

Number simplify(const Number& num, const Number& precision);

Vector simplify(const Vector& vec, const Number& precision);

//Number simplex(stdVector obj, Template M, stdVector N, stdVector bl,
//              stdVector bu, int &err);
//
//Number simplex(stdVector obj, stdMatrix Aineq, stdVector Bineq, stdMatrix Aeq,
//              stdVector Beq, stdVector bl, stdVector bu, int &err);

//bool adjacencyTest(doubleMatrix M, doubleVector N, int face1, int face2);

inline void wait(){
        getchar();
}

// ###############################
// Multi index operation
// ###############################

bool increaseMultiIndice(std::vector<int>& m, const std::vector<int> d);


class Chrono{
public:
        Chrono();
        ~Chrono();
         void start();
         Number stop();
private :
         clock_t startClock;
};

}  // namespace nltool

//using namespace std;

#endif /* UTILS_H_ */