src/include/ApproximationBuilder.h

#ifndef HYBRIDIZATIONSIMPLEX_H_
#define HYBRIDIZATIONSIMPLEX_H_

#include "utils.h"
#include "DynamicalSystem.h"
#include "LinearSystem.h"
#include "PointerSystem.h"
#include "PolynomialSystem.h"
#include "Hpolyhedron.h"
#include "Box.h"
//#include "CurvatureTensor.h"

namespace nltool {

using namespace std;

enum system_type {
        NOT_DEFINED = 1, HESSIAN_CONSTANT = 2, HESSIAN_NOT_CONSTANT = 3
};

class ApproximationBuilder {
public:

        ApproximationBuilder();

        ApproximationBuilder(int dimension);

        ApproximationBuilder(int dimension, DynamicalSystem* sys,
                        bool useT = false);

        virtual ~ApproximationBuilder();

        LinearSystem createApproximation(Hpolyhedron& initialSet,
                        Hpolyhedron& approximationDomain, Number errorBound,
                        bool& requiresSplitting);

        Simplex createSimplex(Hpolyhedron initialSet, Number error);

        Number computeMaximalCurvature();

        Matrix computeTransformationMatrix();

        /*
         * ##########################################
         * GETTTERS and SETTERS
         * ##########################################
         */

        Matrix getTransformation() const;
        void setTransformation(Matrix transformation);

        bool isUseTransformation() const;
        void setUseTransformation(bool useTransformation);

        double getMaxCurvature() {
                return maxCurvature;
        }
        void setMaxCurvature(double maxCurvature);

        void setUseDisplacement(bool useDisplacement) {
                this->useDisplacement = useDisplacement;
        }

        Box getCurvatureEstimationDomain() const {
                return curvatureEstimationDomain;
        }

        void setCurvatureEstimationDomain(Box curvatureEstimationDomain) {
                this->curvatureEstimationDomain = curvatureEstimationDomain;
        }

        unsigned getDimension() {
                return dimension;
        }

        void setNonlinearSystem(DynamicalSystem* nlsys) {
                if (nonlinearSystem != NULL)
                        delete &nonlinearSystem;
                nonlinearSystem = nlsys;
                cout << "dim: " << nonlinearSystem->getDimension() << endl;
                nonlinearSystem->printType();
                nlsys->printType();
        }

private:

        /*
         * ##########################################
         * TRANSFORMATION MATRIX COMPUTATION
         * ##########################################
         */

        Matrix computeTransformationMatrix(Matrix &curvature);

        Matrix computeTransformationMatrix(vector<Matrix>& hessian);

        vector<Matrix> computeHessianMatrix();

        Matrix computeCurvatureMatrix(vector<Matrix>& hessian);

        /*
         * ##########################################
         * SIMPLEX CREATION
         * ##########################################
         */

        Simplex createRegularSimplex(Number radius, int dimension);

        void alignToDirection(Vertices &points, Vector direction);

        bool placeSimplex(Vertices& points, Hpolyhedron& init,
                        Number setMaxSize = 1);

        Vector computeDirectionalVector(Hpolyhedron &initial);

        Vector computeCentroid(Hpolyhedron &poly);

        Vector computeCentroid(Vertices &points);

        void computeCurvatureDomain(Hpolyhedron set);

        void testCurvature(Matrix& curvature, int nbTest = 100);

        // use or not the isotropic transformation to increase the approximation domain volume
        bool useTransformation;

        unsigned dimension; // dimension of the system
        int type; // information on the system
        const Number minCurvature;

        vector<Matrix> cstHessian; //constant hessian matrix list
        vector<vector<DynamicalSystem*> > hessian; //not constant hessian matrix
        Matrix curvature; //curvature matrix
        Matrix transformation; //transformation matrix

        Simplex initialSimplex; //initial equilateral simplex
        bool useDisplacement; // use additional displacement when positionning the
                                                  //approximation domain around the initial set

        Number maxCurvature; // maximal curvature of the nonlinear system
        //TODO: for now we recompute the transformation matrix at each
        //              approximation computation for non quadratic system,
        //              this should be optimized

        Box curvatureEstimationDomain; // domain where the curvature tensor is computed
        Number curvatureEstimationCoef; // a coefficient to determine the size of the curvature estimation domain

        DynamicalSystem* nonlinearSystem;

};

// ###############################
// Error Approximation computation
// ###############################

struct ErrorEstimationParameters {
        DynamicalSystem* originalSystem;
        DynamicalSystem* approximateSystem;
        Hpolyhedron domain;
        int indice;
};

Box estimateApproximationerror(DynamicalSystem* original,
                DynamicalSystem* approximation, Hpolyhedron& domain, Number initialerror);

}
 // namespace nltool

#endif