newton_corrector.hpp

//
//  newton_corrector.hpp
//  Xcode_b2
//
//  Created by Collins, James B. on 4/27/16.
//  Copyright (c) 2016 West Texas A&M University. All rights reserved.
//
 
#ifndef BERTINI_NEWTON_CORRECTOR_HPP
#define BERTINI_NEWTON_CORRECTOR_HPP
 
#include "bertini2/trackers/amp_criteria.hpp"
#include "bertini2/trackers/config.hpp"
#include "bertini2/system/system.hpp"
 
 
namespace bertini{
    namespace tracking{
        namespace correct{
 

 
            class NewtonCorrector
            {
            public:
                
                
                NewtonCorrector(const System& S) : current_precision_(DefaultPrecision())
                {
                    ChangeSystem(S);
                }
 
                
                
                
                
                
                
                
                void Settings(const NewtonConfig& newton_settings)
                {
                    newton_config_ = newton_settings;
                }
                
                
    
                
                
                void ChangePrecision(unsigned new_precision)
                {
                    Precision(std::get< Vec<mpfr_complex> >(f_temp_), new_precision);
                    Precision(std::get< Vec<mpfr_complex> >(step_temp_), new_precision);
                    Precision(std::get< Mat<mpfr_complex> >(J_temp_), new_precision);
 
                    std::get< Eigen::PartialPivLU<Mat<mpfr_complex>> >(LU_) = Eigen::PartialPivLU<Mat<mpfr_complex>>(numTotalFunctions_);
 
                    current_precision_ = new_precision;             
                }
 
 
                unsigned precision() const
                {
                    return current_precision_;
                }
                
                void ChangeSystem(const System& S)
                {
                    numTotalFunctions_ = S.NumTotalFunctions();
                    numVariables_ = S.NumVariables();
                    std::get< Mat<dbl> >(J_temp_).resize(numTotalFunctions_, numVariables_);
                    std::get< Mat<mpfr_complex> >(J_temp_).resize(numTotalFunctions_, numVariables_);
                    std::get< Vec<dbl> >(f_temp_).resize(numTotalFunctions_);
                    std::get< Vec<mpfr_complex> >(f_temp_).resize(numTotalFunctions_);
                    std::get< Vec<dbl> >(step_temp_).resize(numTotalFunctions_);
                    std::get< Vec<mpfr_complex> >(step_temp_).resize(numTotalFunctions_);
                }
 
                
                
                
                
                template <typename ComplexType>
                SuccessCode Correct(Vec<ComplexType> & next_space,
                                       System const& S,
                                       Vec<ComplexType> const& current_space, // pass by value to get a copy of it
                                       ComplexType const& current_time,
                                       NumErrorT const& tracking_tolerance,
                                       unsigned min_num_newton_iterations,
                                       unsigned max_num_newton_iterations)
                {
                    #ifndef BERTINI_DISABLE_ASSERTS
                    assert(max_num_newton_iterations >= min_num_newton_iterations && "max number newton iterations must be at least the min.");
                    #endif
                    
                    Vec<ComplexType>& step_ref = std::get< Vec<ComplexType> >(step_temp_);
                    
                    next_space = current_space;
                    for (unsigned ii = 0; ii < max_num_newton_iterations; ++ii)
                    {
                        //Update the newton iterate by one iteration
                        auto success_code = EvalIterationStep(step_ref, S, next_space, current_time);
                        if(success_code != SuccessCode::Success)
                            return success_code;
                        
                        next_space += step_ref;
                        
                        if ( (step_ref.template lpNorm<Eigen::Infinity>() < tracking_tolerance) && (ii >= (min_num_newton_iterations-1)) )
                            return SuccessCode::Success;
                    }
                    
                    return SuccessCode::FailedToConverge;
 
                    
 
                }
 
                
                
                
                
                
                template <typename ComplexType>
                SuccessCode Correct(Vec<ComplexType> & next_space,
                                       System const& S,
                                       Vec<ComplexType> const& current_space, // pass by value to get a copy of it
                                       ComplexType const& current_time,
                                       NumErrorT const& tracking_tolerance,
                                       unsigned min_num_newton_iterations,
                                       unsigned max_num_newton_iterations,
                                       AdaptiveMultiplePrecisionConfig const& AMP_config)
                {
                    #ifndef BERTINI_DISABLE_ASSERTS
                    assert(max_num_newton_iterations >= min_num_newton_iterations && "max number newton iterations must be at least the min.");
                    #endif
 
                    Vec<ComplexType>& step_ref = std::get< Vec<ComplexType> >(step_temp_);
                    
                    next_space = current_space;
                    for (unsigned ii = 0; ii < max_num_newton_iterations; ++ii)
                    {
                        //Update the newton iterate by one iteration
                        auto success_code = EvalIterationStep(step_ref, S, next_space, current_time);
                        if(success_code != SuccessCode::Success)
                            return success_code;
                        
                        next_space += step_ref;
                        
                        Mat<ComplexType>& J_temp_ref = std::get< Mat<ComplexType> >(J_temp_);
                        Eigen::PartialPivLU< Mat<ComplexType> >& LU_ref = std::get< Eigen::PartialPivLU< Mat<ComplexType> > >(LU_);
                        
                        if ( (step_ref.template lpNorm<Eigen::Infinity>() < tracking_tolerance) && (ii >= (min_num_newton_iterations-1)) )
                            return SuccessCode::Success;
                        
                        NumErrorT norm_J_inverse(LU_ref.solve(RandomOfUnits<ComplexType>(S.NumVariables())).norm());
 
                        if (!amp::CriterionB<ComplexType>(NumErrorT(J_temp_ref.norm()), norm_J_inverse, max_num_newton_iterations - ii, tracking_tolerance, NumErrorT(step_ref.template lpNorm<Eigen::Infinity>()), AMP_config))
                            return SuccessCode::HigherPrecisionNecessary;
                        
                        if (!amp::CriterionC<ComplexType>(norm_J_inverse, next_space, tracking_tolerance, AMP_config))
                            return SuccessCode::HigherPrecisionNecessary;
                    }
                    
                    return SuccessCode::FailedToConverge;
                }
 
                
                
                
                
                
                template <typename ComplexType>
                SuccessCode Correct(Vec<ComplexType> & next_space,
                                       NumErrorT & norm_delta_z,
                                       NumErrorT & norm_J,
                                       NumErrorT & norm_J_inverse,
                                       NumErrorT & condition_number_estimate,
                                       System const& S,
                                       Vec<ComplexType> const& current_space, // pass by value to get a copy of it
                                       ComplexType const& current_time,
                                       NumErrorT const& tracking_tolerance,
                                       unsigned min_num_newton_iterations,
                                       unsigned max_num_newton_iterations,
                                       AdaptiveMultiplePrecisionConfig const& AMP_config)
                {
                    #ifndef BERTINI_DISABLE_ASSERTS
                    assert(max_num_newton_iterations >= min_num_newton_iterations && "max number newton iterations must be at least the min.");
                    #endif
                    
                    Vec<ComplexType>& step_ref = std::get< Vec<ComplexType> >(step_temp_);
                    
                    next_space = current_space;
                    for (unsigned ii = 0; ii < max_num_newton_iterations; ++ii)
                    {
                        //Update the newton iterate by one iteration
                        auto success_code = EvalIterationStep(step_ref, S, next_space, current_time);
                        if(success_code != SuccessCode::Success)
                            return success_code;
                        
                        next_space += step_ref;
                        
                        Mat<ComplexType>& J_temp_ref = std::get< Mat<ComplexType> >(J_temp_);
                        Eigen::PartialPivLU< Mat<ComplexType> >& LU_ref = std::get< Eigen::PartialPivLU< Mat<ComplexType> > >(LU_);
                        
                        
                        norm_delta_z = NumErrorT(step_ref.template lpNorm<Eigen::Infinity>());
                        norm_J = NumErrorT(J_temp_ref.norm());
                        norm_J_inverse = NumErrorT(LU_ref.solve(RandomOfUnits<ComplexType>(S.NumVariables())).norm());
                        condition_number_estimate = NumErrorT(norm_J*norm_J_inverse);
                                                
                        if ( (norm_delta_z < tracking_tolerance) && (ii >= (min_num_newton_iterations-1)) )
                            return SuccessCode::Success;
                        
                        if (!amp::CriterionB<ComplexType>(norm_J, norm_J_inverse, max_num_newton_iterations - ii, tracking_tolerance, norm_delta_z, AMP_config))
                            return SuccessCode::HigherPrecisionNecessary;
                        
                        if (!amp::CriterionC<ComplexType>(norm_J_inverse, next_space, tracking_tolerance, AMP_config))
                            return SuccessCode::HigherPrecisionNecessary;
                    }
                    
                    return SuccessCode::FailedToConverge;
                }
 
                
            private:

                //
                // Private Data Methods
                //
                
                
                
                template<typename ComplexType, typename Derived>
                SuccessCode EvalIterationStep(Vec<ComplexType> & newton_step,
                                              const System& S,
                                              const Eigen::MatrixBase<Derived>& current_space, const ComplexType& current_time)
                {
                    Vec<ComplexType>& f_temp_ref = std::get< Vec<ComplexType> >(f_temp_);
                    Mat<ComplexType>& J_temp_ref = std::get< Mat<ComplexType> >(J_temp_);
                    
                    Eigen::PartialPivLU< Mat<ComplexType> >& LU_ref = std::get< Eigen::PartialPivLU< Mat<ComplexType> > >(LU_);
 
                    S.SetAndReset<ComplexType>(current_space, current_time);
                    S.EvalInPlace(f_temp_ref);
                    S.JacobianInPlace(J_temp_ref);
                    LU_ref = J_temp_ref.lu();
                    
                    if (LUPartialPivotDecompositionSuccessful(LU_ref.matrixLU())!=MatrixSuccessCode::Success)
                        return SuccessCode::MatrixSolveFailure;
                    
                    newton_step = LU_ref.solve(-f_temp_ref);
                    
                    return SuccessCode::Success;
                    
                }
                
 
                
                //
                // Private Data Members
                //
                
                unsigned numTotalFunctions_; // Number of total functions for the current system
                unsigned numVariables_;  // Number of variables for the current system
                
                std::tuple< Vec<dbl>, Vec<mpfr_complex> > f_temp_; // Variable to hold temporary evaluation of the system
                std::tuple< Vec<dbl>, Vec<mpfr_complex> > step_temp_; // Variable to hold temporary evaluation of the newton step
                std::tuple< Mat<dbl>, Mat<mpfr_complex> > J_temp_; // Variable to hold temporary evaluation of the Jacobian
                
                std::tuple< Eigen::PartialPivLU<Mat<dbl>>, Eigen::PartialPivLU<Mat<mpfr_complex>> > LU_; // The LU factorization from the Newton iterates
                
                unsigned current_precision_;
 
                NewtonConfig newton_config_; // Hold the settings of the Newton iteration
 
                
            }; //re: class NewtonCorrector
            
        } //re: namespace correct
    }// re: namespace tracking
}// re: namespace bertini
 
#endif