Commit 82ea64f151f82d12381a7468077338745b34dc32

Authored by Olivier Mullier
1 parent d186d731

ajout des fichiers source

CMakeLists.txt 0 → 100644
... ... @@ -0,0 +1,41 @@
  1 +cmake_minimum_required(VERSION 2.6)
  2 +
  3 +project(integrale)
  4 +
  5 +#Directory of executables
  6 +set(EXECUTABLE_OUTPUT_PATH bin)
  7 +
  8 +set(CMAKE_CXX_FLAGS "-Wall -O2 -frounding-math ")#-pg -fprofile-arcs -ftest-coverage
  9 +
  10 +
  11 +# include_directories(/usr/local/include_test/)
  12 +#
  13 +# link_directories(~/ibex/lib)
  14 +
  15 +include_directories(/usr/local/include/ibex)
  16 +# include_directories(/usr/include/graphviz)
  17 +
  18 +# link_directories(/usr/local/lib)
  19 +link_directories(/usr/local/lib)
  20 +
  21 +
  22 +
  23 +#Source file list generation
  24 +file(
  25 + GLOB_RECURSE
  26 + source_files
  27 + src/*
  28 +)
  29 +
  30 +
  31 +#Executable declaration
  32 +add_executable(
  33 + integrale
  34 + ${source_files}
  35 +)
  36 +
  37 +#Libraries to add for specified target
  38 +target_link_libraries(
  39 + integrale ibex prim z m stdc++
  40 +)
  41 +
... ...
src/integrale.cpp 0 → 100644
... ... @@ -0,0 +1,221 @@
  1 +#include "integrale.h"
  2 +
  3 +integrale::integrale(ibex::Interval inf, ibex::Interval sup, ibex::Function* f, ibex::Variable* var) : _inf_bound(inf), _sup_bound(sup), _f(f)
  4 +{
  5 + _zeros_inf = nullptr;
  6 + _zeros_sup = nullptr;
  7 +}
  8 +
  9 +
  10 + std::list<zero> * integrale::construct_zeros_intern(const ibex::Interval & inf_bound)
  11 +{
  12 +#ifdef __DEBUG__
  13 + __SHOW_FCT_NAME__
  14 +#endif
  15 + //ATTENTION pas de preuve d'unicité
  16 + std::list<zero> * res = new std::list<zero>();
  17 +
  18 + ibex::CtcNewton newton(*_f);
  19 +
  20 + ibex::CtcFwdBwd fb (*_f);
  21 + ibex::RoundRobin bisector(0);
  22 +
  23 + ibex::CellStack buff;
  24 +
  25 + ibex::CtcCompo comp( fb, newton);
  26 +
  27 + ibex::CtcFixPoint fix(comp, 1e-7);
  28 +
  29 + ibex::Solver s(fix, bisector, buff);
  30 +
  31 +
  32 + ibex::IntervalVector inf(1);
  33 + inf[0] = inf_bound;
  34 + std::vector<ibex::IntervalVector> sols = s.solve(inf);
  35 +
  36 + unsigned int i;
  37 + for (i = 0; i < sols.size(); i++){
  38 + //ATTENTION detection foireuse de la pente dans la cas ou [f(temp)] contient 0
  39 + ibex::IntervalVector temp(sols[i].diam());
  40 + if (i < sols.size() - 1){
  41 + for (int j = 0; j < sols[i].size(); j++){
  42 + temp[j] = (sols[i][j] + sols[i+1][j]) / 2.;
  43 + }
  44 + } else {
  45 + for (int j = 0; j < sols[i].size(); j++){
  46 + temp[j] = (sols[i][j] + inf_bound.ub()) / 2.;
  47 + }
  48 + }
  49 + if( (_f->eval_vector(temp)[0].ub() < 0.) ){
  50 + res->push_back(zero(sols[i][0], -1));
  51 + } else {
  52 + res->push_back(zero(sols[i][0], 1));
  53 + }
  54 + }
  55 + return res;
  56 +}
  57 +
  58 +unsigned int integrale::construct_zeros_inf()
  59 +{
  60 + _zeros_inf = construct_zeros_intern(_inf_bound);
  61 + _zeros_inf->sort();
  62 + return _zeros_inf->size();
  63 +}
  64 +
  65 +unsigned int integrale::construct_zeros_sup()
  66 +{
  67 + _zeros_sup = construct_zeros_intern(_sup_bound);
  68 + return _zeros_sup->size();
  69 +}
  70 +
  71 +ibex::Interval integrale::compute_integral(const double& inf_bound, const double& sup_bound)
  72 +{
  73 + ibex::Interval res(0.0);
  74 + double current = inf_bound;
  75 + ibex::IntervalVector current_vec(1);
  76 + while (current < sup_bound){
  77 + current_vec[0] = ibex::Interval(current, current + __EPS__);
  78 + res += _f->eval_vector(current_vec)[0] *__EPS__;
  79 + current += __EPS__;
  80 + }
  81 + return res;
  82 +}
  83 +
  84 +unsigned int integrale::construct_zeros()
  85 +{
  86 + return construct_zeros_inf() + construct_zeros_sup();
  87 +}
  88 +
  89 +integrale_sol integrale::compute_integral(const zero& inf, const zero& sup)
  90 +{
  91 + ibex::Interval res_int = compute_integral(inf.get_x().lb(), sup.get_x().ub());
  92 + return integrale_sol(inf, sup, res_int);
  93 +}
  94 +
  95 +
  96 +std::pair<integrale_sol, integrale_sol> integrale::compute_integral_inf(const ibex::Interval& range)
  97 +{
  98 + std::list<integrale_sol> res_tmp;
  99 + if (!_zeros_inf){
  100 + construct_zeros_inf();
  101 + }
  102 +
  103 + std::list<zero>::iterator it = _zeros_inf->begin();
  104 + zero bound_inf(ibex::Interval(range.lb()), 0);
  105 + integrale_sol integ0 = compute_integral(bound_inf, *(it));
  106 + res_tmp.push_back(integ0);
  107 +
  108 + while (it != _zeros_inf->end()){
  109 + if (*it != _zeros_inf->back()){
  110 + std::list<zero>::iterator it_first = it;
  111 + it++;
  112 + std::list<zero>::iterator it_second = it;
  113 + integrale_sol integ = compute_integral(*it_first, *it_second);
  114 + for(auto elem : res_tmp){
  115 + elem+=integ;
  116 + }
  117 + res_tmp.push_back(integ);
  118 +
  119 + } else {
  120 + zero bound(ibex::Interval(range.ub()), 0);
  121 + integrale_sol integ = compute_integral(*it, bound);
  122 + for(auto elem : res_tmp){
  123 + elem+=integ;
  124 + }
  125 + it++;
  126 + res_tmp.push_back(integ);
  127 + }
  128 + }
  129 + if (res_tmp.size() == 0){
  130 + cout << "Erreur, il devrait y avoir au moins une solution" << endl;
  131 + exit(EXIT_FAILURE);
  132 + }
  133 + integrale_sol min = *(res_tmp.begin());
  134 + integrale_sol max = *(res_tmp.begin());
  135 + for (auto elem : res_tmp){
  136 + if (elem.get_res_integrale().lb() < min.get_res_integrale().lb()){
  137 + min = elem;
  138 + }
  139 + if (elem.get_res_integrale().ub() > max.get_res_integrale().ub()){
  140 + max = elem;
  141 + }
  142 + }
  143 + return pair<integrale_sol, integrale_sol>(min, max);
  144 +}
  145 +
  146 +
  147 +std::pair<integrale_sol, integrale_sol> integrale::compute_integral_sup(const ibex::Interval& range)
  148 +{
  149 + std::list<integrale_sol> res_temp;
  150 + if (!_zeros_sup){
  151 + construct_zeros_sup();
  152 + }
  153 +
  154 + std::list<zero>::iterator it = _zeros_sup->begin();
  155 + zero bound_sup(ibex::Interval(range.ub()), 0);
  156 + std::cout << "inf bound: " << *it << endl;
  157 + std::cout << "sup bound: " << bound_sup << endl;
  158 + integrale_sol integ0 = compute_integral((*it), bound_sup);
  159 + res_temp.push_back(integ0);
  160 +
  161 + while (it != _zeros_sup->end()){
  162 + if (*it != _zeros_sup->back()){
  163 + std::list<zero>::iterator it_second = it;
  164 + it++;
  165 + std::list<zero>::iterator it_first = it;
  166 +// std::cout << "inf bound: " << *it_first << endl;
  167 +// std::cout << "sup bound: " << *it_second << endl;
  168 + integrale_sol integ = compute_integral(*it_first, *it_second);
  169 + for(auto elem : res_temp){
  170 +// cout << "elem : " << elem << endl;
  171 +// cout << "integ : " << integ << endl;
  172 + elem+=integ;
  173 +// cout << "new elem : " << elem << endl;
  174 + }
  175 + res_temp.push_back(integ);
  176 +
  177 + } else {
  178 + zero bound(ibex::Interval(range.lb()), 0);
  179 +// std::cout << "inf bound: " << bound << endl;
  180 +// std::cout << "sup bound: " << *it << endl;
  181 + integrale_sol integ = compute_integral(bound, *it);
  182 + for(auto elem : res_temp){
  183 +// cout << "elem : " << elem << endl;
  184 +// cout << "integ : " << integ << endl;
  185 + elem+=integ;
  186 +// cout << "new elem : " << elem << endl;
  187 + }
  188 + it++;
  189 + res_temp.push_back(integ);
  190 + }
  191 + }
  192 + if (res_temp.size() == 0){
  193 + cout << "Erreur, il devrait y avoir au moins une solution" << endl;
  194 + exit(EXIT_FAILURE);
  195 + }
  196 + integrale_sol min = *(res_temp.begin());
  197 + integrale_sol max = *(res_temp.begin());
  198 + for (auto elem : res_temp){
  199 + cout << "\tcandidate : " << elem << endl;
  200 + if (elem.get_res_integrale().lb() < min.get_res_integrale().lb()){
  201 + min = elem;
  202 + }
  203 + if (elem.get_res_integrale().ub() > max.get_res_integrale().ub()){
  204 + max = elem;
  205 + }
  206 + }
  207 + return pair<integrale_sol, integrale_sol>(min, max);
  208 +}
  209 +
  210 +
  211 +
  212 +
  213 +
  214 +//========================================================================
  215 +std::ostream& operator<<(std::ostream& os, const integrale& inte)
  216 +{
  217 + os << "Integral from " << inte.get_inf_bound() << " to " << inte.get_sup_bound() << " of "
  218 + << *inte.get_f();
  219 + return os;
  220 +}
  221 +
... ...
src/integrale.h 0 → 100644
... ... @@ -0,0 +1,66 @@
  1 +#ifndef _INTEGRALE_
  2 +#define _INTEGRALE_
  3 +#include <ibex.h>
  4 +#include <vector>
  5 +#include "zero.h"
  6 +#include "integrale_sol.h"
  7 +
  8 +// #define __DEBUG__
  9 +
  10 +
  11 +
  12 +#define __PREC__ (1e-3)
  13 +#define __EPS__ (1e-1)
  14 +#define __SHOW_FCT_NAME__ {std::cout << __func__ << std::endl;}
  15 +
  16 +
  17 +
  18 +class integrale{
  19 +private:
  20 +
  21 + ibex::Interval _inf_bound;
  22 + ibex::Interval _sup_bound;
  23 +
  24 + ibex::Function * _f;
  25 + ibex::Variable * _x_var;
  26 +
  27 + std::list<zero> * _zeros_inf;
  28 + std::list<zero> * _zeros_sup;
  29 +
  30 +
  31 +
  32 + std::list<zero> * construct_zeros_intern(const ibex::Interval & inf_bound);
  33 +
  34 +public:
  35 + //CONSTRUCTORS
  36 + integrale(ibex::Interval inf, ibex::Interval sup, ibex::Function * f, ibex::Variable * var);
  37 +
  38 + //GETTERS
  39 + inline ibex::Interval get_inf_bound() const {return _inf_bound;};
  40 + inline ibex::Interval get_sup_bound() const {return _sup_bound;};
  41 + inline ibex::Function * get_f() const {return _f;};
  42 + inline ibex::Variable * get_x_var() const {return _x_var;};
  43 + inline std::list<zero> * get_zeros_inf(){return _zeros_inf;};
  44 + inline std::list<zero> * get_zeros_sup(){return _zeros_sup;};
  45 +
  46 + unsigned int construct_zeros_inf();
  47 +
  48 + unsigned int construct_zeros_sup();
  49 +
  50 + unsigned int construct_zeros();
  51 +
  52 + integrale_sol compute_integral(const zero & inf, const zero & sup);
  53 +
  54 + std::list<ibex::Interval> compute_integral(const double & inf_bound, const ibex::Interval & sup_bound);
  55 +
  56 + ibex::Interval compute_integral(const double & inf_bound, const double & sup_bound);
  57 +
  58 + std::pair<integrale_sol, integrale_sol> compute_integral_inf(const ibex::Interval & range);
  59 +
  60 + std::pair<integrale_sol, integrale_sol> compute_integral_sup(const ibex::Interval & range);
  61 +};
  62 +
  63 +
  64 +std::ostream& operator<<(std::ostream& os, const integrale& inte);
  65 +
  66 +#endif //_INTEGRALE_
... ...
src/integrale_sol.cpp 0 → 100644
... ... @@ -0,0 +1,56 @@
  1 +/*
  2 + * <one line to give the library's name and an idea of what it does.>
  3 + * Copyright (C) 2019 <copyright holder> <email>
  4 + *
  5 + * This program is free software: you can redistribute it and/or modify
  6 + * it under the terms of the GNU General Public License as published by
  7 + * the Free Software Foundation, either version 3 of the License, or
  8 + * (at your option) any later version.
  9 + *
  10 + * This program is distributed in the hope that it will be useful,
  11 + * but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  13 + * GNU General Public License for more details.
  14 + *
  15 + * You should have received a copy of the GNU General Public License
  16 + * along with this program. If not, see <http://www.gnu.org/licenses/>.
  17 + */
  18 +
  19 +#include "integrale_sol.h"
  20 +
  21 +integrale_sol::integrale_sol(const zero& inf, const zero& sup, const ibex::Interval& res) : _res_integrale(res)
  22 +{
  23 + _inf = new zero(inf);
  24 + _sup = new zero(sup);
  25 + _min_candidate = false;
  26 + _max_candidate = false;
  27 + if (res.ub() < 0.){
  28 + _min_candidate = true;
  29 + }
  30 + if (res.lb() > 0.){
  31 + _max_candidate = true;
  32 + }
  33 +}
  34 +
  35 +integrale_sol & integrale_sol::operator+=(const integrale_sol &rhs){
  36 + *this = *this + rhs;
  37 + return *this;
  38 +}
  39 +
  40 +integrale_sol operator+(const integrale_sol& lhs, const integrale_sol& rhs)
  41 +{
  42 + if (lhs.get_sup_bound()->get_x().lb() == rhs.get_sup_bound()->get_x().lb()){
  43 + cout << "ici" << endl;
  44 + return integrale_sol(*lhs.get_inf_bound(), *rhs.get_sup_bound(), (lhs.get_res_integrale() + rhs.get_res_integrale()));
  45 + } else {
  46 + cout << "la" << endl;
  47 + return integrale_sol(*rhs.get_inf_bound(), *lhs.get_sup_bound(), (lhs.get_res_integrale() + rhs.get_res_integrale()));
  48 + }
  49 +
  50 +}
  51 +
  52 +std::ostream & operator<<(std::ostream& os, const integrale_sol& inte)
  53 +{
  54 + os << "Integral from (" << *(inte.get_inf_bound()) << ") to (" << *inte.get_sup_bound() << ") : " << inte.get_res_integrale();
  55 + return os;
  56 +}
... ...
src/integrale_sol.h 0 → 100644
... ... @@ -0,0 +1,54 @@
  1 +/*
  2 + * <one line to give the library's name and an idea of what it does.>
  3 + * Copyright (C) 2019 <copyright holder> <email>
  4 + *
  5 + * This program is free software: you can redistribute it and/or modify
  6 + * it under the terms of the GNU General Public License as published by
  7 + * the Free Software Foundation, either version 3 of the License, or
  8 + * (at your option) any later version.
  9 + *
  10 + * This program is distributed in the hope that it will be useful,
  11 + * but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  13 + * GNU General Public License for more details.
  14 + *
  15 + * You should have received a copy of the GNU General Public License
  16 + * along with this program. If not, see <http://www.gnu.org/licenses/>.
  17 + */
  18 +
  19 +#ifndef INTEGRALE_SOL_H
  20 +#define INTEGRALE_SOL_H
  21 +
  22 +#include "zero.h"
  23 +
  24 +/**
  25 + * @todo write docs
  26 + */
  27 +class integrale_sol
  28 +{
  29 +private:
  30 + zero * _inf;
  31 + zero * _sup;
  32 + ibex::Interval _res_integrale;
  33 + bool _min_candidate;
  34 + bool _max_candidate;
  35 +
  36 +public:
  37 + integrale_sol(const zero & inf, const zero & sup, const ibex::Interval & res);
  38 +
  39 + inline zero * get_inf_bound() const {return _inf;};
  40 + inline zero * get_sup_bound() const {return _sup;};
  41 + inline ibex::Interval get_res_integrale() const {return _res_integrale;};
  42 + inline bool is_min_candidate(){return _min_candidate;};
  43 + inline bool is_max_candidate(){return _max_candidate;};
  44 + void set_min_candidate(bool min){_min_candidate = min;};
  45 + void set_max_candidate(bool max){_max_candidate = max;};
  46 +
  47 +
  48 + integrale_sol & operator+=(const integrale_sol &rhs);
  49 +};
  50 +
  51 +integrale_sol operator+(const integrale_sol &, const integrale_sol &);
  52 +std::ostream& operator<<(std::ostream& os, const integrale_sol& inte_sol);
  53 +
  54 +#endif // INTEGRALE_SOL_H
... ...
src/main.cpp 0 → 100644
... ... @@ -0,0 +1,39 @@
  1 +#include <ibex.h>
  2 +#include "integrale.h"
  3 +
  4 +using namespace std;
  5 +
  6 +
  7 +
  8 +void example(){
  9 + cout << setprecision(15);
  10 + ibex::Variable x(1);
  11 + ibex::Function f(x, ibex::sin(x));
  12 + ibex::Interval inf(0., 5.5);
  13 + ibex::Interval sup(10, 15.);
  14 +
  15 + integrale example(inf, sup, &f, &x);
  16 +
  17 + cout << example << endl;
  18 +
  19 + example.construct_zeros();
  20 +// pair<integrale_sol, integrale_sol> res = example.compute_integral_inf(inf);
  21 +// cout << "minimum: " << res.first << endl;
  22 +// cout << "maximum: " << res.second << endl;
  23 +// for (auto elem : *(example.get_zeros_inf())){
  24 +// cout << elem << endl;
  25 +// }
  26 + cout << "======================================================" << endl;
  27 + pair<integrale_sol, integrale_sol> res2 = example.compute_integral_sup(sup);
  28 + cout << "minimum: " << res2.first << endl;
  29 + cout << "maximum: " << res2.second << endl;
  30 + for (auto elem : *(example.get_zeros_sup())){
  31 + cout << elem << endl;
  32 + }
  33 +}
  34 +
  35 +
  36 +int main(/*int argc, char * argv[]*/){
  37 + example();
  38 + return EXIT_SUCCESS;
  39 +}
... ...
src/zero.cpp 0 → 100644
... ... @@ -0,0 +1,33 @@
  1 +/*
  2 + * <one line to give the library's name and an idea of what it does.>
  3 + * Copyright (C) 2019 <copyright holder> <email>
  4 + *
  5 + * This program is free software: you can redistribute it and/or modify
  6 + * it under the terms of the GNU General Public License as published by
  7 + * the Free Software Foundation, either version 3 of the License, or
  8 + * (at your option) any later version.
  9 + *
  10 + * This program is distributed in the hope that it will be useful,
  11 + * but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  13 + * GNU General Public License for more details.
  14 + *
  15 + * You should have received a copy of the GNU General Public License
  16 + * along with this program. If not, see <http://www.gnu.org/licenses/>.
  17 + */
  18 +
  19 +#include "zero.h"
  20 +
  21 +ostream& operator<<(ostream& os, const zero& z)
  22 +{
  23 + if (z.is_non_zero())
  24 + os << "boundary: " << z.get_x();
  25 + else {
  26 + os << "Zero of the function: " << z.get_x() << ", derivative is ";
  27 + if (z.is_neg())
  28 + os << " negative";
  29 + else
  30 + os << " positive";
  31 + }
  32 + return os;
  33 +}
... ...
src/zero.h 0 → 100644
... ... @@ -0,0 +1,55 @@
  1 +/*
  2 + * <one line to give the library's name and an idea of what it does.>
  3 + * Copyright (C) 2019 <copyright holder> <email>
  4 + *
  5 + * This program is free software: you can redistribute it and/or modify
  6 + * it under the terms of the GNU General Public License as published by
  7 + * the Free Software Foundation, either version 3 of the License, or
  8 + * (at your option) any later version.
  9 + *
  10 + * This program is distributed in the hope that it will be useful,
  11 + * but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  13 + * GNU General Public License for more details.
  14 + *
  15 + * You should have received a copy of the GNU General Public License
  16 + * along with this program. If not, see <http://www.gnu.org/licenses/>.
  17 + */
  18 +
  19 +#ifndef ZERO_H
  20 +#define ZERO_H
  21 +
  22 +#include <ibex.h>
  23 +
  24 +/**
  25 + * @todo write docs
  26 + */
  27 +class zero
  28 +{
  29 +private:
  30 + ibex::Interval _x;
  31 + int _slope;
  32 +
  33 +public:
  34 + zero(const zero &z) : _x(z.get_x()), _slope(z.get_slope()){};
  35 + zero(ibex::Interval x, int slope) : _x(x), _slope(slope) {};
  36 +
  37 + inline ibex::Interval get_x() const {return _x;};
  38 + inline int get_slope() const {return _slope;};
  39 +
  40 + inline bool is_neg() const {return _slope < 0;};
  41 + inline bool is_pos() const {return _slope > 0;};
  42 + inline bool is_non_zero() const {return _slope == 0;};
  43 +};
  44 +
  45 +inline bool operator< (const zero& lhs, const zero& rhs){ return lhs.get_x().mid() < rhs.get_x().mid(); }
  46 +inline bool operator> (const zero& lhs, const zero& rhs){ return rhs < lhs; }
  47 +inline bool operator<=(const zero& lhs, const zero& rhs){ return !(lhs > rhs); }
  48 +inline bool operator>=(const zero& lhs, const zero& rhs){ return !(lhs < rhs); }
  49 +
  50 +inline bool operator==(const zero& lhs, const zero& rhs){ return lhs.get_x() == rhs.get_x(); }
  51 +inline bool operator!=(const zero& lhs, const zero& rhs){ return !(lhs == rhs); }
  52 +
  53 +ostream& operator<<(ostream& os, const zero& z);
  54 +
  55 +#endif // ZERO_H
... ...