Loading src/integrale.cppdeleted 100644 → 0 +0 −289 Original line number Diff line number Diff line #include "integrale.h" integrale::integrale(const ibex::Interval & inf, const ibex::Interval & sup, ibex::Function* f) : _inf_bound(inf), _sup_bound(sup), _f(f) { _zeros_inf = nullptr; _zeros_sup = nullptr; _result = nullptr; } integrale::~integrale() { if (_zeros_inf) delete _zeros_inf; if (_zeros_sup) delete _zeros_sup; if (_result) delete _result; // if (_f) // delete _f; } std::list<zero*> * integrale::construct_zeros_intern(const ibex::Interval & inf_bound) { #ifdef __DEBUG__ __SHOW_FCT_NAME__ #endif //ATTENTION pas de cout << elem << endl;preuve d'unicité std::list<zero*> * res = new std::list<zero*>(); ibex::CtcNewton newton(*_f); ibex::CtcFwdBwd fb (*_f); ibex::RoundRobin bisector(0); ibex::CellStack buff; ibex::CtcCompo comp( fb, newton); ibex::CtcFixPoint fix(comp, 1e-7); ibex::Solver s(fix, bisector, buff); ibex::IntervalVector inf(1); inf[0] = inf_bound; std::vector<ibex::IntervalVector> sols = s.solve(inf); unsigned int i; for (i = 0; i < sols.size(); i++){ //ATTENTION detection foireuse de la pente dans la cas ou [f(temp)] contient 0 ibex::IntervalVector temp(sols[i].diam()); if (i < sols.size() - 1){ for (int j = 0; j < sols[i].size(); j++){ temp[j] = (sols[i][j] + sols[i+1][j]) / 2.; } } else { for (int j = 0; j < sols[i].size(); j++){ temp[j] = (sols[i][j] + inf_bound.ub()) / 2.; } } if( (_f->eval_vector(temp)[0].ub() < 0.) ){ zero *current = new zero(sols[i][0], -1); res->push_back(current); } else { zero *current = new zero(sols[i][0], 1); res->push_back(current); } } return res; } unsigned int integrale::construct_zeros_inf() { _zeros_inf = construct_zeros_intern(_inf_bound); _zeros_inf->sort(); return _zeros_inf->size(); } unsigned int integrale::construct_zeros_sup() { _zeros_sup = construct_zeros_intern(_sup_bound); return _zeros_sup->size(); } ibex::Interval integrale::compute_integral(const double inf, const double sup) { ibex::Interval res(0.0); double current = inf; ibex::IntervalVector current_vec(1); while (current < sup){ current_vec[0] = ibex::Interval(current, current + __EPS__); res += _f->eval_vector(current_vec)[0] *__EPS__; current += __EPS__; } return res; } unsigned int integrale::construct_zeros() { return construct_zeros_inf() + construct_zeros_sup(); } integrale_sol integrale::compute_integral(zero *inf, zero *sup) { ibex::Interval res_inter = compute_integral(inf->get_x().lb(), sup->get_x().ub()); integrale_sol res(inf, sup, res_inter); return res; } std::pair<integrale_sol, integrale_sol> integrale::compute_integral_inf(const ibex::Interval& range) { std::list<integrale_sol> res_tmp; if (!_zeros_inf){ construct_zeros_inf(); } if (_zeros_inf->size()){ std::list<zero*>::iterator it = _zeros_inf->begin(); zero *bound_inf = new zero(ibex::Interval(range.lb()), 0); integrale_sol integ0 = compute_integral(bound_inf, *(it)); res_tmp.push_back(integ0); while (it != _zeros_inf->end()){ if (*it != _zeros_inf->back()){ std::list<zero*>::iterator it_first = it; it++; std::list<zero*>::iterator it_second = it; integrale_sol integ = compute_integral(*it_first, *it_second); std::list<integrale_sol>::iterator it_tmp = res_tmp.begin(); while (it_tmp != res_tmp.end()){ *it_tmp = *it_tmp + integ; it_tmp++; } res_tmp.push_back(integ); } else { zero *bound = new zero(ibex::Interval(range.ub()), 0); integrale_sol integ = compute_integral(*it, bound); std::list<integrale_sol>::iterator it_tmp = res_tmp.begin(); while (it_tmp != res_tmp.end()){ (*it_tmp) = (*it_tmp) + integ; it_tmp++; } it++; res_tmp.push_back(integ); } } } else { zero * bound_sup = new zero(ibex::Interval(range.ub()), 0); zero * bound_inf = new zero(ibex::Interval(range.lb()), 0); integrale_sol integ = compute_integral(bound_inf, bound_sup); res_tmp.push_back(integ); ibex::Interval nul_inter(0.); integrale_sol nul(bound_sup, bound_sup, nul_inter); res_tmp.push_back(nul); } if (res_tmp.size() == 0){ cout << "Erreur, il devrait y avoir au moins une solution" << endl; exit(EXIT_FAILURE); } integrale_sol min = *(res_tmp.begin()); integrale_sol max = *(res_tmp.begin()); for (auto elem : res_tmp){ if (elem.get_res_integrale().lb() < min.get_res_integrale().lb()){ min = elem; } if (elem.get_res_integrale().ub() > max.get_res_integrale().ub()){ max = elem; } } for (auto elem : res_tmp){ if (elem != min || elem != max){ if (elem.get_inf_bound() && !elem.get_inf_bound()->is_non_zero()) delete elem.get_inf_bound(); if (elem.get_sup_bound() && !elem.get_sup_bound()->is_non_zero()) delete elem.get_sup_bound(); } } return pair<integrale_sol, integrale_sol>(min, max); } std::pair<integrale_sol, integrale_sol> integrale::compute_integral_sup(const ibex::Interval& range) { std::list<integrale_sol> res_temp; if (!_zeros_sup){ construct_zeros_sup(); } if (_zeros_sup->size() > 0){ std::list<zero*>::iterator it = _zeros_sup->begin(); zero * bound_sup = new zero(ibex::Interval(range.ub()), 0); integrale_sol integ0 = compute_integral((*it), bound_sup); res_temp.push_back(integ0); while (it != _zeros_sup->end()){ if (*it != _zeros_sup->back()){ std::list<zero*>::iterator it_second = it; it++; std::list<zero*>::iterator it_first = it; integrale_sol integ = compute_integral(*it_first, *it_second); std::list<integrale_sol>::iterator it_tmp = res_temp.begin(); while (it_tmp != res_temp.end()){ *it_tmp = *it_tmp + integ; it_tmp++; } res_temp.push_back(integ); } else { zero * bound = new zero(ibex::Interval(range.lb()), 0); integrale_sol integ = compute_integral(bound, *it); std::list<integrale_sol>::iterator it_tmp = res_temp.begin(); while (it_tmp != res_temp.end()){ *it_tmp = *it_tmp + integ; it_tmp++; } it++; res_temp.push_back(integ); } } } else { zero * bound_sup = new zero(ibex::Interval(range.ub()), 0); zero * bound_inf = new zero(ibex::Interval(range.lb()), 0); integrale_sol integ = compute_integral(bound_inf, bound_sup); res_temp.push_back(integ); ibex::Interval nul_inter(0.); integrale_sol nul(bound_inf, bound_inf, nul_inter); res_temp.push_back(nul); } if (res_temp.size() == 0){ cout << "Erreur, il devrait y avoir au moins une solution" << endl; exit(EXIT_FAILURE); } integrale_sol min = *(res_temp.begin()); integrale_sol max = *(res_temp.begin()); for (auto elem : res_temp){ if (elem.get_res_integrale().lb() < min.get_res_integrale().lb()){ min = elem; } if (elem.get_res_integrale().ub() > max.get_res_integrale().ub()){ max = elem; } } for (auto elem : res_temp){ if (elem != min || elem != max){ if (elem.get_inf_bound() && !elem.get_inf_bound()->is_non_zero()) delete elem.get_inf_bound(); if (elem.get_sup_bound() && !elem.get_sup_bound()->is_non_zero()) delete elem.get_sup_bound(); } } return pair<integrale_sol, integrale_sol>(min, max); } std::pair<integrale_sol, integrale_sol> integrale::compute_integral() { std::pair<integrale_sol, integrale_sol> res = compute_integral_inf(_inf_bound); zero * inf = res.first.get_sup_bound(); std::pair<integrale_sol, integrale_sol> res_sup = compute_integral_sup(_sup_bound); zero * sup = res_sup.first.get_inf_bound(); ibex::Interval mid = compute_integral(_inf_bound.ub(), _sup_bound.lb()); integrale_sol integ_mid(inf, sup, mid); res.first+= integ_mid; res.second+= integ_mid; res.first+= res_sup.first; res.second+= res_sup.second; if (_result == nullptr){ _result = new std::pair<integrale_sol, integrale_sol>(res); } else { *_result = res; } return res; } //======================================================================== std::ostream& operator<<(std::ostream& os, const integrale& inte) { os << "Integral from " << inte.get_inf_bound() << " to " << inte.get_sup_bound() << " of " << *inte.get_f(); return os; } src/integrale.hdeleted 100644 → 0 +0 −70 Original line number Diff line number Diff line #ifndef _INTEGRALE_ #define _INTEGRALE_ #include <ibex.h> #include <vector> #include "zero.h" #include "integrale_sol.h" // #define __DEBUG__ #define __PREC__ (1e-5) #define __EPS__ (1e-1) #define __SHOW_FCT_NAME__ {std::cout << __func__ << std::endl;} class integrale{ private: ibex::Interval _inf_bound; ibex::Interval _sup_bound; ibex::Function * _f; // ibex::Variable * _x_var; std::list<zero*> * _zeros_inf; std::list<zero*> * _zeros_sup; std::pair<integrale_sol, integrale_sol> * _result; std::list<zero*> * construct_zeros_intern(const ibex::Interval & inf_bound); integrale_sol compute_integral(zero* inf, zero* sup); std::list<ibex::Interval> compute_integral(const double & inf_bound, const ibex::Interval & sup_bound); ibex::Interval compute_integral(const double inf_bound, const double sup_bound); std::pair<integrale_sol, integrale_sol> compute_integral_inf(const ibex::Interval & range); std::pair<integrale_sol, integrale_sol> compute_integral_sup(const ibex::Interval & range); public: //CONSTRUCTORS integrale(const ibex::Interval& inf, const ibex::Interval& sup, ibex::Function * f); ~integrale(); //GETTERS inline ibex::Interval get_inf_bound() const {return _inf_bound;}; inline ibex::Interval get_sup_bound() const {return _sup_bound;}; inline ibex::Function * get_f() const {return _f;}; // inline ibex::Variable * get_x_var() const {return _x_var;}; inline std::list<zero*> * get_zeros_inf(){return _zeros_inf;}; inline std::list<zero*> * get_zeros_sup(){return _zeros_sup;}; inline std::pair<integrale_sol, integrale_sol> * get_result(){return _result;}; unsigned int construct_zeros_inf(); unsigned int construct_zeros_sup(); unsigned int construct_zeros(); std::pair<integrale_sol, integrale_sol> compute_integral(); }; std::ostream& operator<<(std::ostream& os, const integrale& inte); #endif //_INTEGRALE_ src/integralii.cpp 0 → 100644 +176 −0 Original line number Diff line number Diff line /* * <one line to give the library's name and an idea of what it does.> * Copyright (C) 2019 mullier <email> * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program. If not, see <http://www.gnu.org/licenses/>. */ #include "integralii.h" std::list<ibex::Interval *>& IntegralII::compute_list_of_zeros(ibex::Interval searchdomain) { //ATTENTION pas de preuve d'unicité std::list<ibex::Interval *>* res = new std::list<ibex::Interval *>(); ibex::CtcNewton newton(*_f); ibex::CtcFwdBwd fb (*_f); ibex::RoundRobin bisector(0); ibex::CellStack buff; ibex::CtcCompo comp( fb, newton); ibex::CtcFixPoint fix(comp, 1e-7); ibex::Solver s(fix, bisector, buff); ibex::IntervalVector inf(1); inf[0] = searchdomain; std::vector<ibex::IntervalVector> sols = s.solve(inf); unsigned int i; for (i = 0; i < sols.size(); i++){ //ATTENTION detection foireuse de la pente dans la cas ou [f(temp)] contient 0 ibex::Interval * current = new ibex::Interval(sols[i][0]); res->push_back(current); } return *res; } std::list<ibex::Interval *>& IntegralII::compute_list_of_zeros_inf() { std::list<ibex::Interval *>& res = compute_list_of_zeros(_inf_endpoint); res.sort(); return res; } std::list<ibex::Interval *>& IntegralII::compute_list_of_zeros_sup() { return compute_list_of_zeros(_sup_endpoint); } void IntegralII::compute_integral() { std::list<IntegralRR> res_temp; std::list<IntegralRR>::iterator it_sols; ibex::Interval *current = new ibex::Interval(_inf_endpoint.lb()); IntegralRR min, max; //======================================================================= //Computation of the integrals in the _inf_endpoint std::list<ibex::Interval *> list_of_zeros_inf = compute_list_of_zeros_inf(); if (list_of_zeros_inf.size()>0){ std::list<ibex::Interval *>::iterator it_zeros_inf = list_of_zeros_inf.begin(); while (it_zeros_inf != list_of_zeros_inf.end()){//calcul des intégrales entre chaque zéro de f IntegralRR current_integral(current->mid(), (*it_zeros_inf)->mid(), _f); it_sols = res_temp.begin(); while (it_sols != res_temp.end()){//On ajoute l'intégrale courante aux autres intégrales *it_sols += current_integral; it_sols++; } res_temp.push_back(current_integral); current = *it_zeros_inf; it_zeros_inf++; } //Dernière intégrale du dernier zéro vers la borne inf de l'intervalle IntegralRR current_integral(current->mid(), _inf_endpoint.ub(), _f); it_sols = res_temp.begin(); while (it_sols != res_temp.end()){//On ajoute l'intégrale courante aux autres intégrales *it_sols += current_integral; it_sols++; } res_temp.push_back(current_integral); } else { IntegralRR whole_integral(_inf_endpoint.lb(), _inf_endpoint.ub(), _f); res_temp.push_back(whole_integral); res_temp.push_back(IntegralRR(_inf_endpoint.ub(), _inf_endpoint.ub(), _f, 0.)); } it_sols = res_temp.begin(); min = *it_sols; max = *it_sols; it_sols++; while (it_sols != res_temp.end()){//On recherche le min et le max if (max < *it_sols){ max = *it_sols; } if (*it_sols < min){ min = *it_sols; } it_sols ++; } delete current; //Calcul de l'intégrale entre les intervalles IntegralRR integrale_milieu(_inf_endpoint.ub(), _sup_endpoint.lb(), _f); _solution.first = min + integrale_milieu; _solution.second = max + integrale_milieu; //======================================================================= current = new ibex::Interval(_sup_endpoint.ub()); res_temp.clear(); std::list<ibex::Interval *> list_of_zeros_sup = compute_list_of_zeros_sup(); if (list_of_zeros_sup.size()>0){ std::list<ibex::Interval *>::iterator it_zeros_sup = list_of_zeros_sup.begin(); while (it_zeros_sup != list_of_zeros_sup.end()){//calcul des intégrales entre chaque zéro de f IntegralRR current_integral((*it_zeros_sup)->mid(), current->mid(), _f); it_sols = res_temp.begin(); while (it_sols != res_temp.end()){//On ajoute l'intégrale courante aux autres intégrales *it_sols += current_integral; it_sols++; } res_temp.push_back(current_integral); current = *it_zeros_sup; it_zeros_sup++; } //Dernière intégrale du dernier zéro vers la borne inf de l'intervalle IntegralRR current_integral(_sup_endpoint.lb(), current->mid(), _f); it_sols = res_temp.begin(); while (it_sols != res_temp.end()){//On ajoute l'intégrale courante aux autres intégrales *it_sols += current_integral; it_sols++; } res_temp.push_back(current_integral); } else { IntegralRR whole_integral(_sup_endpoint.lb(), _sup_endpoint.ub(), _f); res_temp.push_back(whole_integral); res_temp.push_back(IntegralRR(_sup_endpoint.lb(), _sup_endpoint.lb(), _f, 0.)); } it_sols = res_temp.begin(); min = *it_sols; max = *it_sols; it_sols++; while (it_sols != res_temp.end()){//On recherche le min et le max if (max < *it_sols){ max = *it_sols; } if (*it_sols < min){ min = *it_sols; } it_sols ++; } _solution.first += min; _solution.second += max; } std::ostream& operator<<(std::ostream& os, const IntegralII& integral){ os << "Integrale from " << integral.get_inf_endpoint() << " to " << integral.get_sup_endpoint(); os << " of the function " << *integral.get_f() << " : " << endl << "minimum : " << integral.get_solution().first << endl << "maximum : " << integral.get_solution().second; return os; } src/zero.h→src/integralii.h +27 −22 Original line number Diff line number Diff line /* * <one line to give the library's name and an idea of what it does.> * Copyright (C) 2019 <copyright holder> <email> * Copyright (C) 2019 mullier <email> * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by Loading @@ -16,40 +16,45 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. */ #ifndef ZERO_H #define ZERO_H #ifndef INTEGRALII_H #define INTEGRALII_H #include <ibex.h> #include "integralrr.h" /** * @todo write docs */ class zero class IntegralII { private: ibex::Interval _x; int _slope; ibex::Interval _inf_endpoint; ibex::Interval _sup_endpoint; ibex::Function * _f; std::pair<IntegralRR, IntegralRR> _solution; std::list<ibex::Interval*>& compute_list_of_zeros(ibex::Interval searchdomain); public: zero(const zero &z) : _x(z.get_x()), _slope(z.get_slope()){}; zero(ibex::Interval x, int slope) : _x(x), _slope(slope) {}; IntegralII(const ibex::Interval & inf, const ibex::Interval & sup, ibex::Function * f) : _inf_endpoint(inf), _sup_endpoint(sup), _f(f){}; inline ibex::Interval get_x() const {return _x;}; inline int get_slope() const {return _slope;}; ibex::Interval get_inf_endpoint() const {return _inf_endpoint;}; ibex::Interval get_sup_endpoint() const {return _sup_endpoint;}; inline bool is_neg() const {return _slope < 0;}; inline bool is_pos() const {return _slope > 0;}; inline bool is_non_zero() const {return _slope == 0;}; }; ibex::Function * get_f() const {return _f;}; inline bool operator< (const zero& lhs, const zero& rhs){ return lhs.get_x().mid() < rhs.get_x().mid(); } inline bool operator> (const zero& lhs, const zero& rhs){ return rhs < lhs; } inline bool operator<=(const zero& lhs, const zero& rhs){ return !(lhs > rhs); } inline bool operator>=(const zero& lhs, const zero& rhs){ return !(lhs < rhs); } std::pair<IntegralRR, IntegralRR> get_solution() const {return _solution;}; inline bool operator==(const zero& lhs, const zero& rhs){ return lhs.get_x() == rhs.get_x(); } inline bool operator!=(const zero& lhs, const zero& rhs){ return !(lhs == rhs); } std::list<ibex::Interval*>& compute_list_of_zeros_inf(); std::list<ibex::Interval*>& compute_list_of_zeros_sup(); void compute_integral(); }; ostream& operator<<(ostream& os, const zero& z); std::ostream& operator<<(std::ostream& os, const IntegralII& integral); #endif // ZERO_H #endif // INTEGRALEII_H src/integrale_sol.cpp→src/integralrr.cpp +29 −36 Original line number Diff line number Diff line /* * <one line to give the library's name and an idea of what it does.> * Copyright (C) 2019 <copyright holder> <email> * Copyright (C) 2019 mullier <email> * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by Loading @@ -16,61 +16,54 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. */ #include "integrale_sol.h" #include "integralrr.h" integrale_sol::integrale_sol(zero* inf, zero* sup, const ibex::Interval& res) : _res_integrale(res) IntegralRR::IntegralRR(double inf_endpoint, double sup_endpoint, ibex::Function* f) : _inf_endpoint(inf_endpoint), _sup_endpoint(sup_endpoint), _f(f), _solution(0.) { _inf = (inf); _sup = (sup); _min_candidate = false; _max_candidate = false; if (res.ub() < 0.){ _min_candidate = true; } if (res.lb() > 0.){ _max_candidate = true; } compute_integral(); } integrale_sol::~integrale_sol() void IntegralRR::compute_integral() { // if (_inf) // delete _inf; // if (_sup) // delete _sup; ibex::Interval res(0.0); double current = _inf_endpoint; ibex::IntervalVector current_vec(1); while (current < _sup_endpoint){ current_vec[0] = ibex::Interval(current, current + __EPS__); res += _f->eval_vector(current_vec)[0] *__EPS__; current += __EPS__; } current_vec[0] = ibex::Interval(current - __EPS__, _sup_endpoint); res += _f->eval_vector(current_vec)[0] * (_sup_endpoint - current + __EPS__); _solution = res; } integrale_sol & integrale_sol::operator+=(const integrale_sol &rhs){ *this = *this + rhs; IntegralRR& IntegralRR::operator+=(const IntegralRR& other) { *this = *this + other; return *this; } integrale_sol operator+(const integrale_sol& lhs, const integrale_sol& rhs) { if (lhs.get_sup_bound() == rhs.get_inf_bound()){ return integrale_sol(lhs.get_inf_bound(), rhs.get_sup_bound(), (lhs.get_res_integrale() + rhs.get_res_integrale())); IntegralRR operator+(const IntegralRR& lhs, const IntegralRR& rhs) { if (lhs.get_sup_endpoint() == rhs.get_inf_endpoint()){ return IntegralRR(lhs.get_inf_endpoint(), rhs.get_sup_endpoint(), lhs.get_f(), lhs.get_solution() + rhs.get_solution()); } else { return integrale_sol(rhs.get_inf_bound(), lhs.get_sup_bound(), (lhs.get_res_integrale() + rhs.get_res_integrale())); return IntegralRR(rhs.get_inf_endpoint(), lhs.get_sup_endpoint(), lhs.get_f(), lhs.get_solution() + rhs.get_solution()); } } std::ostream & operator<<(std::ostream& os, const integrale_sol& inte) std::ostream & operator<<(std::ostream& os, const IntegralRR& integral) { os << "Integral from (" << *(inte.get_inf_bound()) << ") to (" << *inte.get_sup_bound() << ") : " << inte.get_res_integrale(); os << "Integrale from " << integral.get_inf_endpoint() << " to " << integral.get_sup_endpoint(); os << " of the function " << *integral.get_f() << " : " << integral.get_solution(); return os; } bool operator==(const integrale_sol& lhs, const integrale_sol& rhs) bool operator<(const IntegralRR& lhs, const IntegralRR& rhs) { return (lhs.get_inf_bound() == rhs.get_inf_bound() && lhs.get_sup_bound() == rhs.get_sup_bound()); return lhs.get_solution().ub() < rhs.get_solution().lb(); } bool operator!=(const integrale_sol& lhs, const integrale_sol& rhs) { return !(lhs == rhs); } Loading
src/integrale.cppdeleted 100644 → 0 +0 −289 Original line number Diff line number Diff line #include "integrale.h" integrale::integrale(const ibex::Interval & inf, const ibex::Interval & sup, ibex::Function* f) : _inf_bound(inf), _sup_bound(sup), _f(f) { _zeros_inf = nullptr; _zeros_sup = nullptr; _result = nullptr; } integrale::~integrale() { if (_zeros_inf) delete _zeros_inf; if (_zeros_sup) delete _zeros_sup; if (_result) delete _result; // if (_f) // delete _f; } std::list<zero*> * integrale::construct_zeros_intern(const ibex::Interval & inf_bound) { #ifdef __DEBUG__ __SHOW_FCT_NAME__ #endif //ATTENTION pas de cout << elem << endl;preuve d'unicité std::list<zero*> * res = new std::list<zero*>(); ibex::CtcNewton newton(*_f); ibex::CtcFwdBwd fb (*_f); ibex::RoundRobin bisector(0); ibex::CellStack buff; ibex::CtcCompo comp( fb, newton); ibex::CtcFixPoint fix(comp, 1e-7); ibex::Solver s(fix, bisector, buff); ibex::IntervalVector inf(1); inf[0] = inf_bound; std::vector<ibex::IntervalVector> sols = s.solve(inf); unsigned int i; for (i = 0; i < sols.size(); i++){ //ATTENTION detection foireuse de la pente dans la cas ou [f(temp)] contient 0 ibex::IntervalVector temp(sols[i].diam()); if (i < sols.size() - 1){ for (int j = 0; j < sols[i].size(); j++){ temp[j] = (sols[i][j] + sols[i+1][j]) / 2.; } } else { for (int j = 0; j < sols[i].size(); j++){ temp[j] = (sols[i][j] + inf_bound.ub()) / 2.; } } if( (_f->eval_vector(temp)[0].ub() < 0.) ){ zero *current = new zero(sols[i][0], -1); res->push_back(current); } else { zero *current = new zero(sols[i][0], 1); res->push_back(current); } } return res; } unsigned int integrale::construct_zeros_inf() { _zeros_inf = construct_zeros_intern(_inf_bound); _zeros_inf->sort(); return _zeros_inf->size(); } unsigned int integrale::construct_zeros_sup() { _zeros_sup = construct_zeros_intern(_sup_bound); return _zeros_sup->size(); } ibex::Interval integrale::compute_integral(const double inf, const double sup) { ibex::Interval res(0.0); double current = inf; ibex::IntervalVector current_vec(1); while (current < sup){ current_vec[0] = ibex::Interval(current, current + __EPS__); res += _f->eval_vector(current_vec)[0] *__EPS__; current += __EPS__; } return res; } unsigned int integrale::construct_zeros() { return construct_zeros_inf() + construct_zeros_sup(); } integrale_sol integrale::compute_integral(zero *inf, zero *sup) { ibex::Interval res_inter = compute_integral(inf->get_x().lb(), sup->get_x().ub()); integrale_sol res(inf, sup, res_inter); return res; } std::pair<integrale_sol, integrale_sol> integrale::compute_integral_inf(const ibex::Interval& range) { std::list<integrale_sol> res_tmp; if (!_zeros_inf){ construct_zeros_inf(); } if (_zeros_inf->size()){ std::list<zero*>::iterator it = _zeros_inf->begin(); zero *bound_inf = new zero(ibex::Interval(range.lb()), 0); integrale_sol integ0 = compute_integral(bound_inf, *(it)); res_tmp.push_back(integ0); while (it != _zeros_inf->end()){ if (*it != _zeros_inf->back()){ std::list<zero*>::iterator it_first = it; it++; std::list<zero*>::iterator it_second = it; integrale_sol integ = compute_integral(*it_first, *it_second); std::list<integrale_sol>::iterator it_tmp = res_tmp.begin(); while (it_tmp != res_tmp.end()){ *it_tmp = *it_tmp + integ; it_tmp++; } res_tmp.push_back(integ); } else { zero *bound = new zero(ibex::Interval(range.ub()), 0); integrale_sol integ = compute_integral(*it, bound); std::list<integrale_sol>::iterator it_tmp = res_tmp.begin(); while (it_tmp != res_tmp.end()){ (*it_tmp) = (*it_tmp) + integ; it_tmp++; } it++; res_tmp.push_back(integ); } } } else { zero * bound_sup = new zero(ibex::Interval(range.ub()), 0); zero * bound_inf = new zero(ibex::Interval(range.lb()), 0); integrale_sol integ = compute_integral(bound_inf, bound_sup); res_tmp.push_back(integ); ibex::Interval nul_inter(0.); integrale_sol nul(bound_sup, bound_sup, nul_inter); res_tmp.push_back(nul); } if (res_tmp.size() == 0){ cout << "Erreur, il devrait y avoir au moins une solution" << endl; exit(EXIT_FAILURE); } integrale_sol min = *(res_tmp.begin()); integrale_sol max = *(res_tmp.begin()); for (auto elem : res_tmp){ if (elem.get_res_integrale().lb() < min.get_res_integrale().lb()){ min = elem; } if (elem.get_res_integrale().ub() > max.get_res_integrale().ub()){ max = elem; } } for (auto elem : res_tmp){ if (elem != min || elem != max){ if (elem.get_inf_bound() && !elem.get_inf_bound()->is_non_zero()) delete elem.get_inf_bound(); if (elem.get_sup_bound() && !elem.get_sup_bound()->is_non_zero()) delete elem.get_sup_bound(); } } return pair<integrale_sol, integrale_sol>(min, max); } std::pair<integrale_sol, integrale_sol> integrale::compute_integral_sup(const ibex::Interval& range) { std::list<integrale_sol> res_temp; if (!_zeros_sup){ construct_zeros_sup(); } if (_zeros_sup->size() > 0){ std::list<zero*>::iterator it = _zeros_sup->begin(); zero * bound_sup = new zero(ibex::Interval(range.ub()), 0); integrale_sol integ0 = compute_integral((*it), bound_sup); res_temp.push_back(integ0); while (it != _zeros_sup->end()){ if (*it != _zeros_sup->back()){ std::list<zero*>::iterator it_second = it; it++; std::list<zero*>::iterator it_first = it; integrale_sol integ = compute_integral(*it_first, *it_second); std::list<integrale_sol>::iterator it_tmp = res_temp.begin(); while (it_tmp != res_temp.end()){ *it_tmp = *it_tmp + integ; it_tmp++; } res_temp.push_back(integ); } else { zero * bound = new zero(ibex::Interval(range.lb()), 0); integrale_sol integ = compute_integral(bound, *it); std::list<integrale_sol>::iterator it_tmp = res_temp.begin(); while (it_tmp != res_temp.end()){ *it_tmp = *it_tmp + integ; it_tmp++; } it++; res_temp.push_back(integ); } } } else { zero * bound_sup = new zero(ibex::Interval(range.ub()), 0); zero * bound_inf = new zero(ibex::Interval(range.lb()), 0); integrale_sol integ = compute_integral(bound_inf, bound_sup); res_temp.push_back(integ); ibex::Interval nul_inter(0.); integrale_sol nul(bound_inf, bound_inf, nul_inter); res_temp.push_back(nul); } if (res_temp.size() == 0){ cout << "Erreur, il devrait y avoir au moins une solution" << endl; exit(EXIT_FAILURE); } integrale_sol min = *(res_temp.begin()); integrale_sol max = *(res_temp.begin()); for (auto elem : res_temp){ if (elem.get_res_integrale().lb() < min.get_res_integrale().lb()){ min = elem; } if (elem.get_res_integrale().ub() > max.get_res_integrale().ub()){ max = elem; } } for (auto elem : res_temp){ if (elem != min || elem != max){ if (elem.get_inf_bound() && !elem.get_inf_bound()->is_non_zero()) delete elem.get_inf_bound(); if (elem.get_sup_bound() && !elem.get_sup_bound()->is_non_zero()) delete elem.get_sup_bound(); } } return pair<integrale_sol, integrale_sol>(min, max); } std::pair<integrale_sol, integrale_sol> integrale::compute_integral() { std::pair<integrale_sol, integrale_sol> res = compute_integral_inf(_inf_bound); zero * inf = res.first.get_sup_bound(); std::pair<integrale_sol, integrale_sol> res_sup = compute_integral_sup(_sup_bound); zero * sup = res_sup.first.get_inf_bound(); ibex::Interval mid = compute_integral(_inf_bound.ub(), _sup_bound.lb()); integrale_sol integ_mid(inf, sup, mid); res.first+= integ_mid; res.second+= integ_mid; res.first+= res_sup.first; res.second+= res_sup.second; if (_result == nullptr){ _result = new std::pair<integrale_sol, integrale_sol>(res); } else { *_result = res; } return res; } //======================================================================== std::ostream& operator<<(std::ostream& os, const integrale& inte) { os << "Integral from " << inte.get_inf_bound() << " to " << inte.get_sup_bound() << " of " << *inte.get_f(); return os; }
src/integrale.hdeleted 100644 → 0 +0 −70 Original line number Diff line number Diff line #ifndef _INTEGRALE_ #define _INTEGRALE_ #include <ibex.h> #include <vector> #include "zero.h" #include "integrale_sol.h" // #define __DEBUG__ #define __PREC__ (1e-5) #define __EPS__ (1e-1) #define __SHOW_FCT_NAME__ {std::cout << __func__ << std::endl;} class integrale{ private: ibex::Interval _inf_bound; ibex::Interval _sup_bound; ibex::Function * _f; // ibex::Variable * _x_var; std::list<zero*> * _zeros_inf; std::list<zero*> * _zeros_sup; std::pair<integrale_sol, integrale_sol> * _result; std::list<zero*> * construct_zeros_intern(const ibex::Interval & inf_bound); integrale_sol compute_integral(zero* inf, zero* sup); std::list<ibex::Interval> compute_integral(const double & inf_bound, const ibex::Interval & sup_bound); ibex::Interval compute_integral(const double inf_bound, const double sup_bound); std::pair<integrale_sol, integrale_sol> compute_integral_inf(const ibex::Interval & range); std::pair<integrale_sol, integrale_sol> compute_integral_sup(const ibex::Interval & range); public: //CONSTRUCTORS integrale(const ibex::Interval& inf, const ibex::Interval& sup, ibex::Function * f); ~integrale(); //GETTERS inline ibex::Interval get_inf_bound() const {return _inf_bound;}; inline ibex::Interval get_sup_bound() const {return _sup_bound;}; inline ibex::Function * get_f() const {return _f;}; // inline ibex::Variable * get_x_var() const {return _x_var;}; inline std::list<zero*> * get_zeros_inf(){return _zeros_inf;}; inline std::list<zero*> * get_zeros_sup(){return _zeros_sup;}; inline std::pair<integrale_sol, integrale_sol> * get_result(){return _result;}; unsigned int construct_zeros_inf(); unsigned int construct_zeros_sup(); unsigned int construct_zeros(); std::pair<integrale_sol, integrale_sol> compute_integral(); }; std::ostream& operator<<(std::ostream& os, const integrale& inte); #endif //_INTEGRALE_
src/integralii.cpp 0 → 100644 +176 −0 Original line number Diff line number Diff line /* * <one line to give the library's name and an idea of what it does.> * Copyright (C) 2019 mullier <email> * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program. If not, see <http://www.gnu.org/licenses/>. */ #include "integralii.h" std::list<ibex::Interval *>& IntegralII::compute_list_of_zeros(ibex::Interval searchdomain) { //ATTENTION pas de preuve d'unicité std::list<ibex::Interval *>* res = new std::list<ibex::Interval *>(); ibex::CtcNewton newton(*_f); ibex::CtcFwdBwd fb (*_f); ibex::RoundRobin bisector(0); ibex::CellStack buff; ibex::CtcCompo comp( fb, newton); ibex::CtcFixPoint fix(comp, 1e-7); ibex::Solver s(fix, bisector, buff); ibex::IntervalVector inf(1); inf[0] = searchdomain; std::vector<ibex::IntervalVector> sols = s.solve(inf); unsigned int i; for (i = 0; i < sols.size(); i++){ //ATTENTION detection foireuse de la pente dans la cas ou [f(temp)] contient 0 ibex::Interval * current = new ibex::Interval(sols[i][0]); res->push_back(current); } return *res; } std::list<ibex::Interval *>& IntegralII::compute_list_of_zeros_inf() { std::list<ibex::Interval *>& res = compute_list_of_zeros(_inf_endpoint); res.sort(); return res; } std::list<ibex::Interval *>& IntegralII::compute_list_of_zeros_sup() { return compute_list_of_zeros(_sup_endpoint); } void IntegralII::compute_integral() { std::list<IntegralRR> res_temp; std::list<IntegralRR>::iterator it_sols; ibex::Interval *current = new ibex::Interval(_inf_endpoint.lb()); IntegralRR min, max; //======================================================================= //Computation of the integrals in the _inf_endpoint std::list<ibex::Interval *> list_of_zeros_inf = compute_list_of_zeros_inf(); if (list_of_zeros_inf.size()>0){ std::list<ibex::Interval *>::iterator it_zeros_inf = list_of_zeros_inf.begin(); while (it_zeros_inf != list_of_zeros_inf.end()){//calcul des intégrales entre chaque zéro de f IntegralRR current_integral(current->mid(), (*it_zeros_inf)->mid(), _f); it_sols = res_temp.begin(); while (it_sols != res_temp.end()){//On ajoute l'intégrale courante aux autres intégrales *it_sols += current_integral; it_sols++; } res_temp.push_back(current_integral); current = *it_zeros_inf; it_zeros_inf++; } //Dernière intégrale du dernier zéro vers la borne inf de l'intervalle IntegralRR current_integral(current->mid(), _inf_endpoint.ub(), _f); it_sols = res_temp.begin(); while (it_sols != res_temp.end()){//On ajoute l'intégrale courante aux autres intégrales *it_sols += current_integral; it_sols++; } res_temp.push_back(current_integral); } else { IntegralRR whole_integral(_inf_endpoint.lb(), _inf_endpoint.ub(), _f); res_temp.push_back(whole_integral); res_temp.push_back(IntegralRR(_inf_endpoint.ub(), _inf_endpoint.ub(), _f, 0.)); } it_sols = res_temp.begin(); min = *it_sols; max = *it_sols; it_sols++; while (it_sols != res_temp.end()){//On recherche le min et le max if (max < *it_sols){ max = *it_sols; } if (*it_sols < min){ min = *it_sols; } it_sols ++; } delete current; //Calcul de l'intégrale entre les intervalles IntegralRR integrale_milieu(_inf_endpoint.ub(), _sup_endpoint.lb(), _f); _solution.first = min + integrale_milieu; _solution.second = max + integrale_milieu; //======================================================================= current = new ibex::Interval(_sup_endpoint.ub()); res_temp.clear(); std::list<ibex::Interval *> list_of_zeros_sup = compute_list_of_zeros_sup(); if (list_of_zeros_sup.size()>0){ std::list<ibex::Interval *>::iterator it_zeros_sup = list_of_zeros_sup.begin(); while (it_zeros_sup != list_of_zeros_sup.end()){//calcul des intégrales entre chaque zéro de f IntegralRR current_integral((*it_zeros_sup)->mid(), current->mid(), _f); it_sols = res_temp.begin(); while (it_sols != res_temp.end()){//On ajoute l'intégrale courante aux autres intégrales *it_sols += current_integral; it_sols++; } res_temp.push_back(current_integral); current = *it_zeros_sup; it_zeros_sup++; } //Dernière intégrale du dernier zéro vers la borne inf de l'intervalle IntegralRR current_integral(_sup_endpoint.lb(), current->mid(), _f); it_sols = res_temp.begin(); while (it_sols != res_temp.end()){//On ajoute l'intégrale courante aux autres intégrales *it_sols += current_integral; it_sols++; } res_temp.push_back(current_integral); } else { IntegralRR whole_integral(_sup_endpoint.lb(), _sup_endpoint.ub(), _f); res_temp.push_back(whole_integral); res_temp.push_back(IntegralRR(_sup_endpoint.lb(), _sup_endpoint.lb(), _f, 0.)); } it_sols = res_temp.begin(); min = *it_sols; max = *it_sols; it_sols++; while (it_sols != res_temp.end()){//On recherche le min et le max if (max < *it_sols){ max = *it_sols; } if (*it_sols < min){ min = *it_sols; } it_sols ++; } _solution.first += min; _solution.second += max; } std::ostream& operator<<(std::ostream& os, const IntegralII& integral){ os << "Integrale from " << integral.get_inf_endpoint() << " to " << integral.get_sup_endpoint(); os << " of the function " << *integral.get_f() << " : " << endl << "minimum : " << integral.get_solution().first << endl << "maximum : " << integral.get_solution().second; return os; }
src/zero.h→src/integralii.h +27 −22 Original line number Diff line number Diff line /* * <one line to give the library's name and an idea of what it does.> * Copyright (C) 2019 <copyright holder> <email> * Copyright (C) 2019 mullier <email> * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by Loading @@ -16,40 +16,45 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. */ #ifndef ZERO_H #define ZERO_H #ifndef INTEGRALII_H #define INTEGRALII_H #include <ibex.h> #include "integralrr.h" /** * @todo write docs */ class zero class IntegralII { private: ibex::Interval _x; int _slope; ibex::Interval _inf_endpoint; ibex::Interval _sup_endpoint; ibex::Function * _f; std::pair<IntegralRR, IntegralRR> _solution; std::list<ibex::Interval*>& compute_list_of_zeros(ibex::Interval searchdomain); public: zero(const zero &z) : _x(z.get_x()), _slope(z.get_slope()){}; zero(ibex::Interval x, int slope) : _x(x), _slope(slope) {}; IntegralII(const ibex::Interval & inf, const ibex::Interval & sup, ibex::Function * f) : _inf_endpoint(inf), _sup_endpoint(sup), _f(f){}; inline ibex::Interval get_x() const {return _x;}; inline int get_slope() const {return _slope;}; ibex::Interval get_inf_endpoint() const {return _inf_endpoint;}; ibex::Interval get_sup_endpoint() const {return _sup_endpoint;}; inline bool is_neg() const {return _slope < 0;}; inline bool is_pos() const {return _slope > 0;}; inline bool is_non_zero() const {return _slope == 0;}; }; ibex::Function * get_f() const {return _f;}; inline bool operator< (const zero& lhs, const zero& rhs){ return lhs.get_x().mid() < rhs.get_x().mid(); } inline bool operator> (const zero& lhs, const zero& rhs){ return rhs < lhs; } inline bool operator<=(const zero& lhs, const zero& rhs){ return !(lhs > rhs); } inline bool operator>=(const zero& lhs, const zero& rhs){ return !(lhs < rhs); } std::pair<IntegralRR, IntegralRR> get_solution() const {return _solution;}; inline bool operator==(const zero& lhs, const zero& rhs){ return lhs.get_x() == rhs.get_x(); } inline bool operator!=(const zero& lhs, const zero& rhs){ return !(lhs == rhs); } std::list<ibex::Interval*>& compute_list_of_zeros_inf(); std::list<ibex::Interval*>& compute_list_of_zeros_sup(); void compute_integral(); }; ostream& operator<<(ostream& os, const zero& z); std::ostream& operator<<(std::ostream& os, const IntegralII& integral); #endif // ZERO_H #endif // INTEGRALEII_H
src/integrale_sol.cpp→src/integralrr.cpp +29 −36 Original line number Diff line number Diff line /* * <one line to give the library's name and an idea of what it does.> * Copyright (C) 2019 <copyright holder> <email> * Copyright (C) 2019 mullier <email> * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by Loading @@ -16,61 +16,54 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. */ #include "integrale_sol.h" #include "integralrr.h" integrale_sol::integrale_sol(zero* inf, zero* sup, const ibex::Interval& res) : _res_integrale(res) IntegralRR::IntegralRR(double inf_endpoint, double sup_endpoint, ibex::Function* f) : _inf_endpoint(inf_endpoint), _sup_endpoint(sup_endpoint), _f(f), _solution(0.) { _inf = (inf); _sup = (sup); _min_candidate = false; _max_candidate = false; if (res.ub() < 0.){ _min_candidate = true; } if (res.lb() > 0.){ _max_candidate = true; } compute_integral(); } integrale_sol::~integrale_sol() void IntegralRR::compute_integral() { // if (_inf) // delete _inf; // if (_sup) // delete _sup; ibex::Interval res(0.0); double current = _inf_endpoint; ibex::IntervalVector current_vec(1); while (current < _sup_endpoint){ current_vec[0] = ibex::Interval(current, current + __EPS__); res += _f->eval_vector(current_vec)[0] *__EPS__; current += __EPS__; } current_vec[0] = ibex::Interval(current - __EPS__, _sup_endpoint); res += _f->eval_vector(current_vec)[0] * (_sup_endpoint - current + __EPS__); _solution = res; } integrale_sol & integrale_sol::operator+=(const integrale_sol &rhs){ *this = *this + rhs; IntegralRR& IntegralRR::operator+=(const IntegralRR& other) { *this = *this + other; return *this; } integrale_sol operator+(const integrale_sol& lhs, const integrale_sol& rhs) { if (lhs.get_sup_bound() == rhs.get_inf_bound()){ return integrale_sol(lhs.get_inf_bound(), rhs.get_sup_bound(), (lhs.get_res_integrale() + rhs.get_res_integrale())); IntegralRR operator+(const IntegralRR& lhs, const IntegralRR& rhs) { if (lhs.get_sup_endpoint() == rhs.get_inf_endpoint()){ return IntegralRR(lhs.get_inf_endpoint(), rhs.get_sup_endpoint(), lhs.get_f(), lhs.get_solution() + rhs.get_solution()); } else { return integrale_sol(rhs.get_inf_bound(), lhs.get_sup_bound(), (lhs.get_res_integrale() + rhs.get_res_integrale())); return IntegralRR(rhs.get_inf_endpoint(), lhs.get_sup_endpoint(), lhs.get_f(), lhs.get_solution() + rhs.get_solution()); } } std::ostream & operator<<(std::ostream& os, const integrale_sol& inte) std::ostream & operator<<(std::ostream& os, const IntegralRR& integral) { os << "Integral from (" << *(inte.get_inf_bound()) << ") to (" << *inte.get_sup_bound() << ") : " << inte.get_res_integrale(); os << "Integrale from " << integral.get_inf_endpoint() << " to " << integral.get_sup_endpoint(); os << " of the function " << *integral.get_f() << " : " << integral.get_solution(); return os; } bool operator==(const integrale_sol& lhs, const integrale_sol& rhs) bool operator<(const IntegralRR& lhs, const IntegralRR& rhs) { return (lhs.get_inf_bound() == rhs.get_inf_bound() && lhs.get_sup_bound() == rhs.get_sup_bound()); return lhs.get_solution().ub() < rhs.get_solution().lb(); } bool operator!=(const integrale_sol& lhs, const integrale_sol& rhs) { return !(lhs == rhs); }