/*
 * <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;
}