evolve.h

Go to the documentation of this file.

//
// Created by Giuliano Iorio on 2020-12-26.
//
 
#ifndef SEVN_EVOLVE_H
#define SEVN_EVOLVE_H
 
#ifdef _OPENMP
#include <omp.h>
#endif
 
#include <star.h>
#include <binstar.h>
#include <IO.h>
#include <vector>
#include <random>
#include <errhand.h>
#include <sevnlog.h>
using sevnstd::SevnLogging;
 
//TODO With the new functor implementation, the evolve function can be passed as parameters as used for other options. TO be implemented
//TODO the evolve function should not be put in the SEVN libraries, so we should move it in the main folder
 
namespace evolve_utility{
 
    struct Options{
        std::vector<double> num_parameters;
        std::vector<std::string> str_parameters;
    };
 
 
    class EvolveFunctor {
    public:
        explicit EvolveFunctor(SevnLogging& svlog, bool record_state=true, bool include_failed=false, Options* options=nullptr)
                : _svlog(svlog), _record_state(record_state), _include_failed(include_failed), _options{options} {};
        virtual ~EvolveFunctor() = default;
        virtual inline int operator() (_UNUSED Binstar& binary){ return EXIT_SUCCESS;}
        virtual inline int operator() (_UNUSED Star& star){ return EXIT_SUCCESS;}
        virtual inline std::string name() { return "Pure Virtual Base EvolveFunctor";}
 
        template <typename System>
        static inline bool check_stalling(System& system) {
            static long dt_threshold = long(system.get_svpar_num("check_stalling_time")); //Threshold to check for stalling system in seconds
            if (system.get_svpar_bool("check_stalling")){
 
                static std::chrono::steady_clock::time_point clock_begin;
                static size_t _id{0};
                static std::string _name;
 
                //Set value and start new clock for a new system
                if (system.get_name()!=_name or  system.get_ID()!=_id){
                    clock_begin = std::chrono::steady_clock::now();
                    _id         = system.get_ID();
                    _name       = system.get_name();
                }
 
                if (elapsed_time(clock_begin)>dt_threshold){
                    SevnLogging svlog;
                    svlog.critical("System with name " + system.get_name() +
                                   " and ID " + utilities::n2s(system.get_ID(),__FILE__,__LINE__) +
                                   " evolved for more than " + utilities::n2s(dt_threshold,__FILE__,__LINE__) +
                                   "s this is not expected."
                                   "The evolution is skipped and the system flagged ad failed.",
                                   __FILE__,
                                   __LINE__,
                                   sevnstd::sevnerr(""));
                }
 
            }
 
            return EXIT_SUCCESS;
        }
 
    protected:
        SevnLogging& _svlog;
        const bool _record_state;
        const bool _include_failed;
        const Options* _options;
 
        template <class System>
        inline void final_print(System& system, bool evolution_result){
 
            if (evolution_result==EXIT_SUCCESS)
                system.print();
            else
                system.print_failed(_include_failed);
        }
 
        template <class System>
        static void catched_error_message(System& system, const sevnstd::sevnerr &e){
            std::string errmedd= "Evolution of system "+system.get_name() + "(ID "+ std::to_string(system.get_ID())+
                                 ") broken by an error: "+e.what();
            std::cerr<<errmedd<<std::endl;
        }
 
        static inline long elapsed_time(std::chrono::steady_clock::time_point start_clock){
            return std::chrono::duration_cast<std::chrono::seconds>(std::chrono::steady_clock::now()-start_clock).count();
        }
    };
 
    class EvolveDebug : public EvolveFunctor {
    public:
        explicit EvolveDebug(SevnLogging& svlog, bool record_state=true, bool include_failed=false)
                : EvolveFunctor(svlog,record_state,include_failed) {};
        inline int operator() (Binstar& binary) override  {
            return evolve_debug(binary);
        }
        inline int operator() (Star& star) override {
            return evolve_debug(star);
        }
        inline std::string name() override { return "EvolveDefault";}
 
    private:
        template <typename System>
        inline int evolve_debug(System& system){
 
            bool evolution_result = EXIT_SUCCESS; //Reset to true
            final_print(system,evolution_result); //print out all the recorded states for the star
            /********************************/
 
            return evolution_result;
        }
    };
 
 
    class EvolveDefault : public EvolveFunctor {
    public:
        explicit EvolveDefault(SevnLogging& svlog, bool record_state=true, bool include_failed=false)
                : EvolveFunctor(svlog,record_state,include_failed) {};
        inline int operator() (Binstar& binary) override  {
            return evolve_default(binary);
        }
        inline int operator() (Star& star) override {
            return evolve_default(star);
        }
        inline std::string name() override { return "EvolveDefault";}
 
    private:
        template <typename System>
        inline int evolve_default(System& system){
 
            bool evolution_result = EXIT_SUCCESS; //Reset to true
 
            if (_record_state) system.recordstate(); //always record the initial state
 
            for(;;) {
                //TODO Myabe is better to use try..catch.. outside the for?
                try {
                    check_stalling(system); //Check for stalling systems
                    system.evolve();
                    if (system.breaktrigger()) //stopping condition for evolving the star
                        break; //go to evolve the next star
 
                    if ( (system.isoutputtime() or system.printall()) and _record_state) {
                        system.recordstate_w_timeupdate();
                    }
                }
                catch(sevnstd::sevnerr& e){
 
                    catched_error_message(system,e);
 
                    evolution_result=EXIT_FAILURE;
                    break;
                }
            }
            if (_record_state) {
                system.recordstate(); //always record the final state
            }
 
 
            //print out all the recorded states for the star
            final_print(system,evolution_result); //print out all the recorded states for the star
            /********************************/
 
            return evolution_result;
        }
    };
 
    class EvolveStopCondition : public EvolveFunctor{
    public:
        explicit EvolveStopCondition(SevnLogging& svlog, bool record_state=true, bool include_failed=false)
        : EvolveFunctor(svlog,record_state,include_failed) {};
        virtual inline int operator() (Binstar& binary) override {
 
 
            bool evolution_result = EXIT_SUCCESS; //Reset to true
 
            for(;;) {
 
                try {
                    check_stalling(binary); //Check for stalling systems
                    binary.evolve();
 
                    //If special break condition is flagged and _record_state is true:
                    //  record_state, print it to the file and exit.
                    if (stop_condition(binary) and _record_state){
                        binary.recordstate(); //always record the final state
                        break; //Exit
                    }
                    //If special break condition or natural breaktrigger  (controlled by input options) just break
                    else if (stop_condition(binary) or binary.breaktrigger())
                        break;
                }
                catch(sevnstd::sevnerr& e){
 
                    catched_error_message(binary,e);
 
                    evolution_result = EXIT_FAILURE;
                    break;
 
                }
            }
            /********************************/
            final_print(binary,evolution_result); //print out all the recorded states for the star
            return evolution_result;
        }
        virtual inline int operator() (Star& star) override {
 
 
            bool evolution_result = EXIT_SUCCESS; //Reset to true
 
            for(;;) {
 
                try {
                    check_stalling(star); //Check for stalling systems
                    star.evolve();
 
                    //If special break condition is flagged and _record_state is true:
                    //  record_state, print it to the file and exit.
                    if (stop_condition(star) and _record_state){
                        star.recordstate(); //always record the final state
                        break; //Exit
                    }
                        //If special break condition or natural breaktrigger  (controlled by input options) just break
                    else if (stop_condition(star) or star.breaktrigger())
                        break;
                }
                catch(sevnstd::sevnerr& e){
 
                    catched_error_message(star,e);
 
                    evolution_result = EXIT_FAILURE;
                    break;
 
                }
            }
            /********************************/
            final_print(star,evolution_result); //print out all the recorded states for the star
            return evolution_result;
        }
        inline std::string name() override { return "EvolveStopCondition";}
    protected:
        virtual inline bool stop_condition(_UNUSED Binstar& binstar){
            return false;
        }
 
        virtual inline bool stop_condition(_UNUSED Star& star){
            return false;
        }
    };
 
    class EvolveRecordCondition : public  EvolveFunctor{
    public:
        explicit EvolveRecordCondition(SevnLogging& svlog, bool record_state=true, bool include_failed=false)
                : EvolveFunctor(svlog,record_state,include_failed) {};
        virtual inline int operator() (Binstar& binary) override {
 
 
            bool evolution_result = EXIT_SUCCESS; //Reset to true
            unsigned int Record_counter=0;
 
            for(;;) {
 
                try {
                    check_stalling(binary); //Check for stalling systems
                    binary.evolve();
 
                    //If record condition is true, update the BLC counter and record state (if allowed by other flags)
                    if (record_condition(binary)){
                        Record_counter++;
                        //Record only if it is time (or print only end and this is the start of the Record condition Record_counter=1) and if _record state is true
                        if ( ( (binary.isoutputtime() or binary.printall()) or Record_counter==1 ) and _record_state ) {
                            binary.recordstate_w_timeupdate();
                        }
                    }
                        //If record_condition is false, but Record_counter>0 we just exit from the record condition, so record state (if allowed by other flags) ans reset the Record_counter
                    else if(Record_counter>0){
                        //Always print the end
                        if (_record_state){
                            binary.recordstate_w_timeupdate();
                        }
                        //reset the Record_counter counter
                        Record_counter=0;
                    }
 
                    if (binary.breaktrigger()) //stopping condition for evolving the star
                        break; //go to evolve the next star
 
                }
                catch(sevnstd::sevnerr& e){
 
 
                    catched_error_message(binary,e);
 
                    evolution_result = EXIT_FAILURE;
                    break;
 
                }
            }
            /********************************/
            final_print(binary,evolution_result); //print out all the recorded states for the star
            return evolution_result;
        }
        virtual inline int operator() (Star& star) override {
 
 
            bool evolution_result = EXIT_SUCCESS; //Reset to true
            unsigned int Record_counter=0;
 
            for(;;) {
 
                try {
                    check_stalling(star); //Check for stalling systems
                    star.evolve();
 
                    //If record condition is true, update the BLC counter and record state (if allowed by other flags)
                    if (record_condition(star)){
                        Record_counter++;
                        //Record only if it is time (or print only end and this is the start of the Record condition Record_counter=1) and if _record state is true
                        if ( ( (star.isoutputtime() or star.printall()) or Record_counter==1 ) and _record_state ) {
                            star.recordstate_w_timeupdate();
                        }
                    }
                        //If record_condition is false, but Record_counter>0 we just exit from the record condition, so record state (if allowed by other flags) ans reset the Record_counter
                    else if(Record_counter>0){
                        //Always print the end
                        if (_record_state){
                            star.recordstate_w_timeupdate();
                        }
                        //reset the Record_counter counter
                        Record_counter=0;
                    }
 
                    if (star.breaktrigger()) //stopping condition for evolving the star
                        break; //go to evolve the next star
 
                }
                catch(sevnstd::sevnerr& e){
 
                    catched_error_message(star,e);
 
                    evolution_result = EXIT_FAILURE;
                    break;
 
                }
            }
            /********************************/
            final_print(star,evolution_result); //print out all the recorded states for the star
            return evolution_result;
        }
        inline std::string name() override { return "EvolveRecordCondition";}
    protected:
 
 
        virtual inline bool record_condition(_UNUSED Binstar& binstar){
            return false;
        }
 
        virtual inline bool record_condition(_UNUSED Star& star){
            return false;
        }
    };
 
 
    class EvolveBinaryCompact : public EvolveStopCondition{
    public:
        explicit EvolveBinaryCompact(SevnLogging& svlog, bool record_state=true, bool include_failed=false)
        : EvolveStopCondition(svlog,record_state,include_failed) {};
        inline int operator() (_UNUSED Star& star) override {
            _svlog.critical("Evolve for star is not implemented for functor EvolveBinaryCompact",__FILE__,__LINE__,
                            sevnstd::notimplemented_error());
            return -1;
        }
 
        inline std::string name() override { return "EvolveBinaryCompact";}
 
    protected:
        inline bool stop_condition(Binstar& binstar) override {
            return binstar.getstar(0)->amiCompact() and binstar.getstar(1)->amiCompact() and !binstar.broken;
        }
 
    };
 
     //TODO Deprecated to be removed, new implementation is EvolveBinaryCompact
    class EvolveBinaryCompactOld : public EvolveFunctor {
    public:
        explicit EvolveBinaryCompactOld(SevnLogging& svlog, bool record_state=true, bool include_failed=false)
                : EvolveFunctor(svlog,record_state,include_failed) {};
        inline int operator() (Binstar& binary) override {
 
            bool special_break_condition = false;
            bool evolution_result = EXIT_SUCCESS; //Reset to true
 
            for(;;) {
 
                try {
                    check_stalling(binary); //Check for stalling systems
 
                    binary.evolve();
 
                    special_break_condition  = binary.getstar(0)->amiCompact() and binary.getstar(1)->amiCompact() and !binary.broken;
 
                    //If special break condition is flagged and _record_state is true:
                    //  record_state, print it to the file and exit.
                    if (special_break_condition and _record_state){
                        binary.recordstate(); //always record the final state
                        break; //Exit
                    }
                    //If special break condition or natural breaktrigger  (controlled by input options) just break
                    else if (special_break_condition or binary.breaktrigger())
                        break;
                }
                catch(sevnstd::sevnerr& e){
 
                    catched_error_message(binary,e);
 
                    evolution_result = EXIT_FAILURE;
                    break;
 
                }
            }
            /********************************/
            final_print(binary,evolution_result); //print out all the recorded states for the star
            return evolution_result;
        }
        inline int operator() (_UNUSED Star& star) override {
            _svlog.critical("Evolve for star is not implemented for functor EvolveBinaryCompact",__FILE__,__LINE__,
                            sevnstd::notimplemented_error());
            return -1;
        }
 
        inline std::string name() override { return "EvolveBinaryCompact";}
 
    };
 
    class EvolveBBH :  public EvolveStopCondition{
    public:
        explicit EvolveBBH(SevnLogging& svlog, bool record_state=true, bool include_failed=false)
        : EvolveStopCondition(svlog,record_state,include_failed) {};
        inline int operator() (_UNUSED Star& star) override {
            _svlog.critical("Evolve for star is not implemented for functor EvolveBinaryCompact",__FILE__,__LINE__,
                            sevnstd::notimplemented_error());
            return -1;
        }
 
        inline std::string name() override { return "EvolveBBH";}
 
    protected:
        inline bool stop_condition(Binstar& binstar) override {
            return binstar.getstar(0)->amiBH() and binstar.getstar(1)->amiBH() and !binstar.broken;
        }
 
    };
 
 
    class EvolveBLC : public EvolveFunctor {
    public:
        explicit EvolveBLC(SevnLogging& svlog, bool record_state=true, bool include_failed=false)
                : EvolveFunctor(svlog,record_state,include_failed) {};
        inline int operator() (Binstar& binary) override {
 
            bool BLC_flag = false, BLC_case1=false, BLC_case2=false;
            //The BLC counter counts the evolution steps the systems has been in the BLC configuration
            //It is increased by one if the BLC_flag is zero, while it is reset to zero when BLC_flag is false
            unsigned int BLC_counter=0;
            bool evolution_result = EXIT_SUCCESS; //Reset to true
 
            for(;;) {
                check_stalling(binary); //Check for stalling systems
 
                try {
 
                    binary.evolve();
 
                    BLC_case1 = binary.getstar(0)->amiCompact() and !binary.getstar(1)->amiCompact();
                    BLC_case2 = binary.getstar(1)->amiCompact() and !binary.getstar(0)->amiCompact();
 
                    BLC_flag  = (BLC_case1 or BLC_case2) and !binary.broken;
 
                    //If BLC_flag is true, update the BLC counter and record state (if allowed by other flags)
                    if (BLC_flag){
                        BLC_counter++;
                        //Record only if it is time (or print only end and this is the start of the BLC conditions BLC_counter=) and if _record state is true
                        if ( ( (binary.isoutputtime() or binary.printall()) or BLC_counter==1 ) and _record_state ) {
                            binary.recordstate_w_timeupdate();
                        }
                    }
                    //If BLC_flag is false, but BLC_counter>0 we are at the end of a BLC phase, so record state (if allowed by other flags) ans reset the BLC counter
                    //FInally break, in fact if this system was in a BLC phase and not it is not, there is no way it will enter again in a BLC phase
                    else if(BLC_counter>0){
 
                        //Always print the end
                        if (_record_state){
                            binary.recordstate_w_timeupdate();
                        }
                        //reset the BLC counter
                        BLC_counter=0;
                        break;
                    }
 
                    if (binary.breaktrigger()) //stopping condition for evolving the star
                        break; //go to evolve the next star
 
                }
                catch(sevnstd::sevnerr& e){
 
                    catched_error_message(binary,e);
 
                    evolution_result=EXIT_FAILURE;
                    break;
                }
            }
            /********************************/
            final_print(binary,evolution_result); //print out all the recorded states for the star
            return evolution_result;
        }
        inline int operator() (_UNUSED Star& star) override {
            _svlog.critical("Evolve for star is not implemented for functor EvolveBinaryCompact",__FILE__,__LINE__,
                            sevnstd::notimplemented_error());
            return -1;
        }
        inline std::string name() override { return "EvolveBLC";}
 
 
    };
 
 
 
 
    class EvolveBCX1 : public EvolveStopCondition {
    public:
        explicit EvolveBCX1(SevnLogging& svlog, bool record_state=true, bool include_failed=false)
                : EvolveStopCondition(svlog,record_state,include_failed) {};
        inline int operator() (_UNUSED Star& star) override {
            _svlog.critical("Evolve for star is not implemented for functor EvolveBinaryCompact",__FILE__,__LINE__,
                            sevnstd::notimplemented_error());
            return -1;
        }
 
        inline std::string name() override { return "EvolveBCX1";}
 
    protected:
        inline bool stop_condition(Binstar& binstar) override {
            bool BH_cond = binstar.getstar(0)->amiBH() and binstar.getstar(0)->getp(Mass::ID)>=14 and binstar.getstar(0)->getp(Mass::ID)<=24;
            bool  MS_cond = binstar.getstar(1)->getp(Phase::ID)==1 and binstar.getstar(1)->getp(Mass::ID)>=30 and binstar.getstar(1)->getp(Mass::ID)<=50;
            bool Semimajor_cond = binstar.getp(Semimajor::ID)<=1000000000000.0;
 
            return BH_cond and MS_cond and Semimajor_cond;
        }
 
    };
 
    class EvolveBHMSTarantula : public EvolveStopCondition {
    public:
        explicit EvolveBHMSTarantula(SevnLogging& svlog, bool record_state=true, bool include_failed=false)
                : EvolveStopCondition(svlog,record_state,include_failed) {};
        inline int operator() (_UNUSED Star& star) override {
            _svlog.critical("Evolve for star is not implemented for functor EvolveBinaryCompact",__FILE__,__LINE__,
                            sevnstd::notimplemented_error());
            return -1;
        }
 
        inline std::string name() override { return "EvolveBHMSTarantula";}
 
    protected:
        inline bool stop_condition(Binstar& binstar) override {
 
            constexpr double  yr_to_day = 265.24;
 
            bool BH_cond = binstar.getstar(0)->amiBH() and binstar.getstar(0)->getp(Mass::ID)>=6 and binstar.getstar(0)->getp(Mass::ID)<=16;
            bool MS_cond = binstar.getstar(1)->getp(Phase::ID)==1 and binstar.getstar(1)->getp(Mass::ID)>=15 and binstar.getstar(1)->getp(Mass::ID)<=32;
            bool P_cond  = binstar.getp(Period::ID)*yr_to_day < 100;
            bool e_cond  = binstar.getp(Eccentricity::ID) < 0.3;
 
            return BH_cond and MS_cond and P_cond and e_cond;
        }
 
    };
 
    class EvolveBHMS : public EvolveRecordCondition {
    public:
        explicit EvolveBHMS(SevnLogging& svlog, bool record_state=true, bool include_failed=false)
                : EvolveRecordCondition(svlog,record_state,include_failed) {};
        inline int operator() (_UNUSED Star& star) override {
            _svlog.critical("Evolve for star is not implemented for functor EvolveBinaryCompact",__FILE__,__LINE__,
                            sevnstd::notimplemented_error());
            return -1;
        }
 
        inline std::string name() override { return "EvolveBHMSTarantula";}
 
    protected:
        inline bool record_condition(Binstar& binstar) override {
            bool cond_1      = binstar.getstar(0)->amiBH() and (binstar.getstar(1)->getp(Phase::ID)<7);
            bool cond_2      = binstar.getstar(1)->amiBH() and (binstar.getstar(0)->getp(Phase::ID)<7);
            bool binary_cond = !binstar.broken;
 
            return (cond_1 or cond_2) and binary_cond;
        }
 
    };
 
 
    class EvolveW1 : public EvolveRecordCondition {
 
    public:
        explicit EvolveW1(SevnLogging& svlog, bool record_state=true, bool include_failed=false)
                : EvolveRecordCondition(svlog,record_state,include_failed) {};
        inline int operator() (_UNUSED Star& star) override {
            _svlog.critical("Evolve for star is not implemented for functor EvolveBinaryCompact",__FILE__,__LINE__,
                            sevnstd::notimplemented_error());
            return -1;
        }
 
        inline std::string name() override { return "EvolveW1";}
 
    protected:
        inline bool record_condition(Binstar& binstar) override {
            bool BHNS_cond = (binstar.getstar(0)->amiCompact() and !binstar.getstar(1)->amiCompact()) or (binstar.getstar(1)->amiCompact() and !binstar.getstar(0)->amiCompact());
            bool Period_cond = !binstar.broken and binstar.getp(Period::ID)>=Pmin and binstar.getp(Period::ID)<=Pmax;
 
            return BHNS_cond and Period_cond;
        }
    private:
        static constexpr double hours_to_yr = 1.0/8766.0; //1 hours to year, considering 1 year=365.25 days
        static constexpr double Pmin = 1*hours_to_yr; //Minimum Period 1h
        static constexpr double Pmax = 5*hours_to_yr; //Maximum Period 5h
 
    };
 
 
    class EvolveRRL : public EvolveRecordCondition {
 
    public:
        explicit EvolveRRL(SevnLogging& svlog, bool record_state=true, bool include_failed=false)
                : EvolveRecordCondition(svlog,record_state,include_failed) {};
        inline std::string name() override { return "EvolveRRL";}
 
    protected:
        inline bool record_condition(Binstar& binstar) override {
            return binstar.getstar(0)->amiRRL() or binstar.getstar(1)->amiRRL();
        }
        inline bool record_condition(Star& star) override {
            return star.amiRRL() or star.amiRRL();
        }
    };
 
 
    inline int evolve_single(Binstar& binary, SevnLogging& svlog, bool record_state=true){
 
        if (record_state) binary.recordstate(); //always record the initial state
        for(;;) {
 
            try {
                EvolveFunctor::check_stalling(binary); //Check stalling
                binary.evolve();
 
 
                if (binary.breaktrigger()) //stopping condition for evolving the star
                    break; //go to evolve the next star
 
                if ( (binary.isoutputtime() or binary.printall()) and record_state ) {
                    binary.recordstate_w_timeupdate();
                }
            }
            catch(sevnstd::sevnerr& e){
 
 
                svlog.error("Evolution of binary "+binary.get_name() + "(ID "+ utilities::n2s(binary.get_ID(),__FILE__,__LINE__)+
                            ") broken by an error: "+e.what(),__FILE__,__LINE__,false);
 
                //std::cerr<< binaries[i].get_name()  <<" " <<e.what()<<std::endl << std::flush;
                return EXIT_FAILURE;
            }
 
            //utilities::wait();
 
 
        }
        if (record_state) {
            binary.recordstate(); //always record the final state
            binary.print(); //print out all the recorded states for the star
        }
        /********************************/
 
        return EXIT_SUCCESS;
    }
 
    inline int evolve_single(Star& star, SevnLogging& svlog, bool record_state=true){
 
        if (record_state) star.recordstate(); //always record the initial state
        for(;;) {
 
            try {
 
                EvolveFunctor::check_stalling(star); //Check stalling
                star.evolve();
 
 
                if (star.breaktrigger()) //stopping condition for evolving the star
                    break; //go to evolve the next star
 
                if ( (star.isoutputtime() or star.printall()) and record_state ) {
 
                    star.recordstate_w_timeupdate();
                }
            }
            catch(sevnstd::sevnerr& e){
 
                svlog.error("Evolution of binary "+star.get_name() + "(ID "+ utilities::n2s(star.get_ID(),__FILE__,__LINE__)+
                            ") broken by an error: "+e.what(),__FILE__,__LINE__,false);
 
                return EXIT_FAILURE;
            }
 
            //utilities::wait();
 
 
        }
        if (record_state) {
            star.recordstate(); //always record the final state
            star.print(); //print out all the recorded states for the star
        }
        /********************************/
 
        return EXIT_SUCCESS;
    }
 
 
    template <typename System>
    inline int evolve_list(EvolveFunctor* evolve_function, std::vector<System>& systems, _UNUSED IO& sevnio, int Nevolve=-1){
        SevnLogging svlog;
 
        if (Nevolve==-1) Nevolve=systems.size();
        unsigned Nfailed =0;
 
#ifdef _OPENMP
#pragma omp parallel num_threads(sevnio.nthreads)
#endif
        {
 
#ifdef _OPENMP
#pragma omp for schedule(static)  reduction(+: Nfailed)
#endif
            for (size_t i = 0; i < (size_t)Nevolve; i++) {
 
                //NOTICE: THE FOLLOWING IS EXTREMELY IMPORTANT
                //mtrand is a static threadlocal variaible, i.e. it is a global variable for each thread
                //it is "seeded" each time we initiliase a system (binary or star), but we initiliase them
                //before to evolve them, therefore is you not "re-seed" mtrand before the evolution
                //mtrand will be seeded with the last loaded system and this will make impossibile to reproduce our results
                //unless we use exactly the same list of systems (at the same relative position in the input file)
                //Notice also that it is important that mtrand IS NEVER RE-SEEDED during the stellar/binary evolution
                //(for this reason when we create an auxiliary star we save the current mtrand status, then we initiliase the
                //auxiliary star where mtrand is reseed and finally we restore the old mtrand)
                //TODO We should surely improve this behaviour (see issue 71 on gitlab: https://gitlab.com/sevn/SEVN/-/issues/71 )
                utilities::mtrand.seed(systems[i].get_rseed());//Set random state for riproducibility, DO NOT MOVE OR REMOVE THIS LINE
 
 
                /****** EVOLVE  ******/
                if((*evolve_function)(systems[i])==EXIT_FAILURE)
                    Nfailed++;
 
 
            }
 
 
        }
 
        return Nfailed;
    }
 
 
    //[[deprecated("Replaced by functor evolve implementation")]] This attribute can be used only from C++14
    template<typename T>
    inline int chunk_dispatcher(unsigned int Nchunk, IO& sevnio, std::vector<T>& systems, bool record_state=true, bool progress=true){
 
        //If systems is not an empty vector clear it
        if (!systems.empty())
            systems.clear();
 
        //Max size
        unsigned int Ntot = sevnio.STARS_MATRIX.size();
        //Reset Nchunk if needed
        Nchunk = std::min(Nchunk,Ntot);
        //Reserve space for systems
        systems.reserve(Nchunk);
 
 
        if (progress)
            std::cout<<"Evolving systems:"<<std::endl;
 
        unsigned int Ndone=0;
        size_t current_idx=0;
        int Ntodo=0;
        int Nfailed=0;
        while (Ndone<Ntot){
 
 
            //Assign Ntodo
            Ntodo = std::min(Ntot-Ndone, Nchunk);
 
            //Fill vector
            for (size_t i = 0; i < Ntodo; i++){
                systems.emplace_back(&sevnio, sevnio.STARS_MATRIX[current_idx], current_idx);
                current_idx++;
            }
 
            //MEMORY DEBUG
            //utilities::hardwait("After filling systems",__FILE__,__LINE__);
 
            //Evolve and update Nfailed
            Nfailed+=evolve_list(systems,sevnio,Ntodo,record_state);
 
            //MEMORY DEBUG
            //utilities::hardwait("After evolving systems",__FILE__,__LINE__);
 
            //Clear
            systems.clear();
 
            //MEMORY DEBUG
            //utilities::hardwait("After clearing systems",__FILE__,__LINE__);
 
            //Update Ndone
            Ndone+=Ntodo;
 
            //TODO 1: Here we can estimate time needed to perform a chunk run and estimate the time to the end
            //TODO 2: Maybe Nfailed is not needed because svlog has already static counters to warning, error and critical messages
            if (progress){
                std::cout<<"\r"<<Ndone<<"/"<<Ntot<<" (Nfailed:"<<Nfailed<<")";
                std::cout<<std::flush;
            }
 
 
 
        }
 
        if (progress)
            std::cout<<std::endl;
 
        return EXIT_SUCCESS;
    }
 
    template<typename T>
    inline int chunk_dispatcher(EvolveFunctor* evolve_function,unsigned int Nchunk, IO& sevnio, std::vector<T>& systems, bool progress=true){
 
        //If systems is not an empty vector clear it
        if (!systems.empty())
            systems.clear();
 
        SevnLogging sevnlog;
 
        //Max size
        unsigned int Ntot = sevnio.STARS_MATRIX.size();
        //Reset Nchunk if needed
        Nchunk = std::min(Nchunk,Ntot);
        //Reserve space for systems
        systems.reserve(Nchunk);
 
        if (progress)
            std::cout<<"Evolving systems:"<<std::endl;
 
        unsigned int Ndone=0;
        size_t current_idx=0;
        int Ntodo=0;
        int Nfailed=0;
        while (Ndone<Ntot){
 
 
            //Assign Ntodo
            Ntodo = std::min(Ntot-Ndone, Nchunk);
 
            //Fill vector
            for (size_t i = 0; i < (size_t)Ntodo; i++){
 
                try{
                    systems.emplace_back(&sevnio, sevnio.STARS_MATRIX[current_idx], current_idx);
                }
                catch(sevnstd::sevnio_error& e){ //sevnio error contains initialisation errors
                    sevnio.print_failed_initilisation_summary(current_idx);
                    sevnlog.error("Failed initilisation for System with ID="+utilities::n2s(current_idx,__FILE__,__LINE__)+
                    " with message:\n"+e.what(),__FILE__,__LINE__,sevnio.svpar.get_bool("initerror_stop"));
                }
                current_idx++;
 
 
            }
 
            //MEMORY DEBUG
            //utilities::hardwait("After filling systems",__FILE__,__LINE__);
            //Evolve and update Nfailed
            //Nfailed+=evolve_list(evolve_function,systems,sevnio,Ntodo);
            //We use system.size instead of Ntodo because wome of the systema could be not initliased due to initilisation errors
            Nfailed+=evolve_list(evolve_function,systems,sevnio,systems.size());
            //MEMORY DEBUG
            //utilities::hardwait("After evolving systems",__FILE__,__LINE__);
            //Clear
            systems.clear();
 
            //MEMORY DEBUG
            //utilities::hardwait("After clearing systems",__FILE__,__LINE__);
 
            //Update Ndone
            Ndone+=Ntodo;
 
            //TODO 1: Here we can estimate time needed to perform a chunk run and estimate the time to the end
            //TODO 2: Maybe Nfailed is not needed because svlog has already static counters to warning, error and critical messages
            if (progress){
                std::cout<<"\r"<<Ndone<<"/"<<Ntot<<" (Nfailed:"<<Nfailed<<")";
                std::cout<<std::flush;
            }
 
 
 
        }
 
        if (progress)
            std::cout<<std::endl;
 
        return EXIT_SUCCESS;
    }
 
 
 
}
 
 
#endif //SEVN_EVOLVE_H