integrale.cpp 9.32 KB
#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;
}