src/include/PolynomialSystem.h

/*
 * PolynomialSystem.h
 *
 *  Created on: Oct 14, 2010
 *      Author: romain
 */
#include "time.h"
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "ostream"
#include "vector"
#include "iostream"
#include "assert.h"
#include "utils.h"
//#include "Bbox.h"
#include "Box.h"
#include "DynamicalSystem.h"

#ifndef POLYNOMIALSYSTEM_H_
#define POLYNOMIALSYSTEM_H_

namespace nltool {

using namespace std;

typedef std::vector<int> multiIndice;

struct Monom {
        Number scalar;
        multiIndice coef;

        Monom();
        Monom(int dim, Number s, multiIndice c);
        Monom(int dim, string strEq, bool sign);
        virtual ~Monom();
        int dimension;
        Number eval(const Vector &x);
        Monom derivate(const int dx);

};

struct Polynom {
        std::vector<Monom> monomials;

        Polynom();
        Polynom(int dim, std::string strEq);
        Polynom(int dim, const std::vector<Monom> monoms);

        virtual ~Polynom();

        void clean();

        virtual Number eval(const Vector &x);

        void linearInterpollation(stdVector point, stdVector &a, Number b);

        Polynom derivate(const int dx);

        Polynom replaceVariable(std::vector<Polynom> expression);

        multiIndice computeMaxDegrees();

        bool containsMultiIndice(multiIndice i, int& position);

        Number computeBernsteinCoefficient(multiIndice i, multiIndice l);

        stdVector getBernsteinCoefficients();

        void computeAffinebound(stdVector al, Number bl, stdVector au, Number bu);

        int dimension;

};

class PolynomialSystem: public DynamicalSystem {
public:
        PolynomialSystem(int dim, int deg, int nbMon, const bool stable, int s);

        PolynomialSystem(int dim, std::vector<Polynom> polys);

        PolynomialSystem(int dim, string polysys);

        PolynomialSystem();

        virtual ~PolynomialSystem();

        Vector eval(const Vector &x);

        DynamicalSystem& derivate(const unsigned dx);

        Number maximize(const unsigned indice, const Template& constraintsA,
                        const Vector& constraintsB, Vector& pt);

        Number minimize(const unsigned indice, const Template& constraintsA,
                        const Vector& constraintsB, Vector& pt);

        void linearInterpollation(Vector point, Matrix &A, Vector &b);

        Number eval(const Vector &x, const int indice);

        Vector evalDiscretised(const Vector &x, Number dt);

        Number curvature(const Vector &x, const int indice, const Vector &d1,
                        const Vector &d2);

        PolynomialSystem getDiscretizedSystem(Number dt);

        PolynomialSystem replaceVariable(std::vector<Polynom> expression);

        PolynomialSystem derivate(const int dx);

        //#################################
        // Bernstein part
        //#################################

        //Matrix computeBernsteinCoefficient();
//      std::vector<Polynom> computeUnitBoxMap(Bbox box);
        multiIndice getMaxDegree();

        std::vector<Polynom> polynomials;
        multiIndice maxDegrees;

private:

        const char* toString();
        bool inline equal(PolynomialSystem s2) {
                return (std::string(toString()) == std::string(s2.toString()));
        }

        static Number evalFunc(const int N, const Vector &x,
                        const int indiceParameter, const Number sign, void *data) {
                PolynomialSystem *t = static_cast<PolynomialSystem*>(data);
                return t->eval(x, indiceParameter);
        }

        static Number evalNumberDerivatedFunc(stdVector &pt, int d, int x, int y,
                        void *data) {
                return 0.0;
                // TODO
        }

        static Number curvatureFunc(const int systemId, const Vector &x,
                        const Vector &v1, const Vector &v2, void* data = NULL) {
                PolynomialSystem *t = static_cast<PolynomialSystem*>(data);
                return t->curvature(x, systemId, v1, v2);
        }

        SystemType findSystemType();

        int degree;
        int nbMonom;
        std::vector<multiIndice> multiIndices;
        //Matrix coefficients;
        stdVector domainU;
        stdVector domainL;

};

//#################################
// Operator
//#################################

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

inline std::ostream& operator<<(std::ostream& out, const Monom &monom) {
        bool prev = false;
        //if(fabs(monom.scalar)>1)
        out << monom.scalar << "*";
        //else return out<<0;
        for (unsigned int i = 0; i < monom.coef.size(); i++) {
                if (prev && monom.coef[i] >= 1)
                        out << "*";
                if (monom.coef[i] > 1)
                        out << "x" << i << "^" << monom.coef[i];
                if (monom.coef[i] == 1)
                        out << "x" << i;
                if (monom.coef[i] >= 1)
                        prev = true;
        }
        return out;
}

inline std::ostream& operator<<(std::ostream& out, const Polynom &poly) {
        for (unsigned int i = 0; i < poly.monomials.size(); i++) {
                out << poly.monomials[i];
                if (i != poly.monomials.size() - 1)
                        out << " + ";
        }
        return out;
}

inline std::ostream& operator<<(std::ostream& out, const PolynomialSystem sys) {
        for (unsigned int i = 0; i < (sys.polynomials).size(); i++) {
                out << "x" << i << "': " << (sys.polynomials[i]) << std::endl;
        }
        return out;
}

// Monom and polynomilas operations

inline Monom operator*(const Monom m1, const Monom m2) {
        std::vector<int> coef(m1.coef.size(), 0);
        Number scalar = 0.0;
        for (unsigned int i = 0; i < coef.size(); i++) {
                coef[i] = m1.coef[i] + m2.coef[i];
        }
        scalar = m1.scalar * m2.scalar;
        Monom res(m1.dimension, scalar, coef);
        return res;
}

inline Polynom operator*(const Polynom poly1, const Polynom poly2) {
        std::vector<Monom> monoms;
        for (unsigned int i = 0; i < poly1.monomials.size(); i++) {
                for (unsigned int j = 0; j < poly2.monomials.size(); j++) {
                        monoms.push_back(poly1.monomials[i] * poly2.monomials[j]);
                }
        }
        Polynom res(poly1.dimension, monoms);
        res.clean();
        return res;
}

inline Polynom operator*(Number d, const Polynom poly) {
        std::vector<Monom> monoms = poly.monomials;

        for (unsigned int i = 0; i < poly.monomials.size(); i++) {
                monoms[i].scalar *= d;
        }
        Polynom res(poly.dimension, monoms);
        res.clean();
        return res;
}

inline Polynom operator+(const Polynom p1, const Polynom p2) {
        Polynom res = p1;
        for (unsigned int i = 0; i < p2.monomials.size(); i++) {
                bool notPresent = true;
                for (unsigned int j = 0; j < res.monomials.size(); j++) {
                        if (res.monomials[j].coef == p2.monomials[i].coef) {
                                notPresent = false;
                                res.monomials[j].scalar += p2.monomials[i].scalar;
                        }
                }
                if (notPresent)
                        res.monomials.push_back(p2.monomials[i]);
        }
        res.clean();
        return res;
}

inline Polynom operator*(stdVector v, const PolynomialSystem polySys) {
        assert(v.size() == polySys.getDimension());
        //        cout << "vector:" << v << endl;
        Polynom res = v[0] * polySys.polynomials[0];
        //cout << "i=" << 0 << " res="<<res<<endl;
        for (unsigned int i = 1; i < polySys.polynomials.size(); i++) {
                res = res + (v[i] * polySys.polynomials[i]);
                //      cout << "i=" << i << " res="<<res<<endl;
        }
        res.clean();
        return res;
}

inline Polynom operator*(Vector v, const PolynomialSystem polySys) {
        assert(v.n_elem == polySys.getDimension());
        //        cout << "vector:" << v << endl;
        Polynom res = v[0] * polySys.polynomials[0];
        //cout << "i=" << 0 << " res="<<res<<endl;
        for (unsigned int i = 1; i < polySys.polynomials.size(); i++) {
                res = res + (v[i] * polySys.polynomials[i]);
                //      cout << "i=" << i << " res="<<res<<endl;
        }
        res.clean();
        return res;
}

inline bool operator==(const multiIndice ind1, const multiIndice ind2) {
        assert(ind1.size() == ind2.size());
        for (unsigned int i = 0; i < ind1.size(); i++)
                if (!(ind1[i] == ind2[i]))
                        return false;
        return true;
}

inline stdVector operator/(const multiIndice ind1, const multiIndice ind2) {
        assert(ind1.size() == ind2.size());
        stdVector result(ind1.size());
        for (unsigned int i = 0; i < ind1.size(); i++) {
                if (ind2[i] != 0)
                        result[i] = ((Number) ind1[i] / (Number) ind2[i]);
                else
                        result[i] = 0;
        }

        return result;
}

inline bool operator<=(const multiIndice ind1, const multiIndice ind2) {
        assert(ind1.size() == ind2.size());
        for (unsigned int i = 0; i < ind1.size(); i++)
                if (ind1[i] > ind2[i])
                        return false;
        return true;
}

inline int sum(const multiIndice ind) {
        int res = 0;
        for (unsigned int i = 0; i < ind.size(); i++) {
                res += ind[i];
        }
        return res;
}

} // namespace nltool

#endif /* POLYNOMIALSYSTEM_H_ */