You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
1029 lines
34 KiB
1029 lines
34 KiB
// This file is part of Bertini 2.0.
|
|
//
|
|
// trackers/explicit_predictors.hpp 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.
|
|
//
|
|
// trackers/explicit_predictors.hpp 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 trackers/explicit_predictors.hpp. If not, see
|
|
// <http://www.gnu.org/licenses/>.
|
|
//
|
|
// copyright 2015
|
|
// James B. Collins
|
|
// West Texas A&M University
|
|
// Department of Mathematics
|
|
// Spring 2016
|
|
//
|
|
// silviana amethyst, university of wisconsin-eau claire
|
|
|
|
/**
|
|
\file explicit_predictors.hpp
|
|
|
|
\brief Contains a base class for all ODE predictors.
|
|
*/
|
|
|
|
#ifndef BERTINI_EXPLICIT_PREDICTORS_HPP
|
|
#define BERTINI_EXPLICIT_PREDICTORS_HPP
|
|
|
|
#include <Eigen/LU>
|
|
#include <boost/type_index.hpp>
|
|
|
|
#include "eigen_extensions.hpp"
|
|
#include "mpfr_extensions.hpp"
|
|
#include "system/system.hpp"
|
|
#include "trackers/amp_criteria.hpp"
|
|
|
|
namespace bertini {
|
|
namespace tracking {
|
|
namespace predict {
|
|
|
|
// Constant,
|
|
// Euler,
|
|
// Heun,
|
|
// RK4,
|
|
// HeunEuler,
|
|
// RKNorsett34,
|
|
// RKF45,
|
|
// RKCashKarp45,
|
|
// RKDormandPrince56,
|
|
// RKVerner67
|
|
|
|
/**
|
|
\brief Get the Bertini2 default predictor.
|
|
|
|
Currently set to RKF45.
|
|
|
|
\return The default predictor method to use.
|
|
*/
|
|
inline Predictor DefaultPredictor() { return Predictor::RKF45; }
|
|
|
|
/**
|
|
\brief The order of the predictor.
|
|
|
|
The order of the error estimate is this plus one.
|
|
|
|
\return The order of the predictor method.
|
|
\param predictor_choice The predictor method to query.
|
|
*/
|
|
inline unsigned Order(Predictor predictor_choice) {
|
|
switch (predictor_choice) {
|
|
case (Predictor::Constant):
|
|
return 0;
|
|
case (Predictor::Euler):
|
|
return 1;
|
|
case (Predictor::Heun):
|
|
return 2;
|
|
case (Predictor::HeunEuler):
|
|
return 2;
|
|
case (Predictor::RKNorsett34):
|
|
return 3;
|
|
case (Predictor::RK4):
|
|
return 4;
|
|
case (Predictor::RKF45):
|
|
return 4;
|
|
case (Predictor::RKCashKarp45):
|
|
return 4;
|
|
case (Predictor::RKDormandPrince56):
|
|
return 5;
|
|
case (Predictor::RKVerner67):
|
|
return 6;
|
|
default: {
|
|
throw std::runtime_error("incompatible predictor choice in Order");
|
|
}
|
|
}
|
|
}
|
|
|
|
/**
|
|
\brief Ask whether a predictor method provides an error estimate.
|
|
|
|
\return Yes or no, does it or does it not.
|
|
\param predictor_choice The predictor method to query.
|
|
*/
|
|
inline bool HasErrorEstimate(Predictor predictor_choice) {
|
|
switch (predictor_choice) {
|
|
case (Predictor::Constant):
|
|
return false;
|
|
case (Predictor::Euler):
|
|
return false;
|
|
case (Predictor::Heun):
|
|
return false;
|
|
case (Predictor::HeunEuler):
|
|
return true;
|
|
case (Predictor::RKNorsett34):
|
|
return true;
|
|
case (Predictor::RK4):
|
|
return false;
|
|
case (Predictor::RKF45):
|
|
return true;
|
|
case (Predictor::RKCashKarp45):
|
|
return true;
|
|
case (Predictor::RKDormandPrince56):
|
|
return true;
|
|
case (Predictor::RKVerner67):
|
|
return true;
|
|
default: {
|
|
throw std::runtime_error(
|
|
"incompatible predictor choice in HasErrorEstimate");
|
|
}
|
|
}
|
|
}
|
|
|
|
namespace {
|
|
template <typename T>
|
|
struct LUSelector {};
|
|
|
|
template <>
|
|
struct LUSelector<dbl> {
|
|
template <typename N>
|
|
static Eigen::PartialPivLU<Mat<dbl>>& Run(N& n) {
|
|
return n.GetLU_d();
|
|
}
|
|
};
|
|
|
|
template <>
|
|
struct LUSelector<mpfr_complex> {
|
|
template <typename N>
|
|
static Eigen::PartialPivLU<Mat<mpfr_complex>>& Run(N& n) {
|
|
return n.GetLU_mp();
|
|
}
|
|
};
|
|
} // namespace
|
|
|
|
/**
|
|
\class ExplicitRKPredictor
|
|
|
|
\brief A class which stores all the explicit single-step multi-stage ODE
|
|
predictor methods.
|
|
|
|
## Purpose
|
|
Stores all the information needed to implement the predictor method
|
|
- Butcher Table
|
|
- Number of Stages
|
|
- Order of the method.
|
|
|
|
Also stores information computed during implementation of the method.
|
|
|
|
|
|
## Use
|
|
Each predictor method is stored as a static Butcher table. To perform a
|
|
predict step, you must instantiate an object with a particular predictor method
|
|
and call Predict:
|
|
|
|
\code
|
|
ExplicitRKPredictors<Complex,Real> euler(Predictor::Euler, sys)
|
|
success_code = euler.Predict( ... )
|
|
\endcode
|
|
*/
|
|
class ExplicitRKPredictor {
|
|
friend LUSelector<dbl>;
|
|
friend LUSelector<mpfr_complex>;
|
|
|
|
public:
|
|
/**
|
|
\brief Construct a predictor to work on a system.
|
|
|
|
\param S the system the predictor will be predicting on.
|
|
*/
|
|
ExplicitRKPredictor(const System& S)
|
|
: current_precision_(DefaultPrecision()), s_(0) {
|
|
ChangeSystem(S);
|
|
PredictorMethod(DefaultPredictor());
|
|
}
|
|
|
|
/**
|
|
\brief Constructor for a particular predictor method
|
|
|
|
\param method The predictor method to be implemented.
|
|
\param S the system to be predicting on.
|
|
*/
|
|
ExplicitRKPredictor(Predictor method, const System& S)
|
|
: current_precision_(DefaultPrecision()), s_(0) {
|
|
ChangeSystem(S);
|
|
PredictorMethod(method);
|
|
}
|
|
|
|
/**
|
|
/brief Sets the local variables to correspond to a particular predictor
|
|
method
|
|
|
|
\param method Enum class that determines the predictor method
|
|
*/
|
|
void PredictorMethod(Predictor method) {
|
|
predictor_ = method;
|
|
p_ = predict::Order(method);
|
|
switch (method) {
|
|
case Predictor::Constant: {
|
|
s_ = 1;
|
|
Mat<double>& arefd = std::get<Mat<double>>(a_);
|
|
Vec<double>& brefd = std::get<Vec<double>>(b_);
|
|
Vec<double>& crefd = std::get<Vec<double>>(c_);
|
|
crefd.resize(s_);
|
|
crefd(0) = 0;
|
|
arefd.resize(s_, s_);
|
|
arefd(0, 0) = 0;
|
|
brefd.resize(s_);
|
|
brefd(0) = 0;
|
|
Mat<mpfr_float>& arefmp = std::get<Mat<mpfr_float>>(a_);
|
|
Vec<mpfr_float>& brefmp = std::get<Vec<mpfr_float>>(b_);
|
|
Vec<mpfr_float>& crefmp = std::get<Vec<mpfr_float>>(c_);
|
|
crefmp.resize(s_);
|
|
crefmp(0) = 0;
|
|
arefmp.resize(s_, s_);
|
|
arefmp(0, 0) = 0;
|
|
brefmp.resize(s_);
|
|
brefmp(0) = 0;
|
|
uses_embedded_ = false;
|
|
|
|
break;
|
|
}
|
|
case Predictor::Euler: {
|
|
s_ = 1;
|
|
Mat<double>& arefd = std::get<Mat<double>>(a_);
|
|
Vec<double>& brefd = std::get<Vec<double>>(b_);
|
|
Vec<double>& crefd = std::get<Vec<double>>(c_);
|
|
crefd.resize(s_);
|
|
crefd(0) = static_cast<double>(cEuler_(0));
|
|
arefd.resize(s_, s_);
|
|
arefd(0, 0) = static_cast<double>(aEuler_(0, 0));
|
|
brefd.resize(s_);
|
|
brefd(0) = static_cast<double>(bEuler_(0));
|
|
Mat<mpfr_float>& arefmp = std::get<Mat<mpfr_float>>(a_);
|
|
Vec<mpfr_float>& brefmp = std::get<Vec<mpfr_float>>(b_);
|
|
Vec<mpfr_float>& crefmp = std::get<Vec<mpfr_float>>(c_);
|
|
crefmp.resize(s_);
|
|
crefmp(0) = static_cast<mpfr_float>(cEuler_(0));
|
|
arefmp.resize(s_, s_);
|
|
arefmp(0, 0) = static_cast<mpfr_float>(aEuler_(0, 0));
|
|
brefmp.resize(s_);
|
|
brefmp(0) = static_cast<mpfr_float>(bEuler_(0));
|
|
uses_embedded_ = false;
|
|
break;
|
|
}
|
|
case Predictor::HeunEuler: {
|
|
s_ = 2;
|
|
|
|
FillButcherTable<double>(s_, aHeunEuler_, bHeunEuler_,
|
|
b_minus_bstarHeunEuler_, cHeunEuler_);
|
|
FillButcherTable<mpfr_float>(s_, aHeunEuler_, bHeunEuler_,
|
|
b_minus_bstarHeunEuler_, cHeunEuler_);
|
|
|
|
break;
|
|
}
|
|
case Predictor::RK4: {
|
|
s_ = 4;
|
|
|
|
FillButcherTable<double>(s_, aRK4_, bRK4_, cRK4_);
|
|
FillButcherTable<mpfr_float>(s_, aRK4_, bRK4_, cRK4_);
|
|
|
|
break;
|
|
}
|
|
|
|
case Predictor::RKF45: {
|
|
s_ = 6;
|
|
|
|
FillButcherTable<double>(s_, aRKF45_, bRKF45_, b_minus_bstarRKF45_,
|
|
cRKF45_);
|
|
FillButcherTable<mpfr_float>(s_, aRKF45_, bRKF45_, b_minus_bstarRKF45_,
|
|
cRKF45_);
|
|
|
|
break;
|
|
}
|
|
|
|
case Predictor::RKCashKarp45: {
|
|
s_ = 6;
|
|
|
|
FillButcherTable<double>(s_, aRKCK45_, bRKCK45_, b_minus_bstarRKCK45_,
|
|
cRKCK45_);
|
|
FillButcherTable<mpfr_float>(s_, aRKCK45_, bRKCK45_,
|
|
b_minus_bstarRKCK45_, cRKCK45_);
|
|
|
|
break;
|
|
}
|
|
|
|
case Predictor::RKDormandPrince56: {
|
|
s_ = 8;
|
|
|
|
FillButcherTable<double>(s_, aRKDP56_, bRKDP56_, b_minus_bstarRKDP56_,
|
|
cRKDP56_);
|
|
FillButcherTable<mpfr_float>(s_, aRKDP56_, bRKDP56_,
|
|
b_minus_bstarRKDP56_, cRKDP56_);
|
|
|
|
break;
|
|
}
|
|
|
|
case Predictor::RKVerner67: {
|
|
s_ = 10;
|
|
|
|
FillButcherTable<double>(s_, aRKV67_, bRKV67_, b_minus_bstarRKV67_,
|
|
cRKV67_);
|
|
FillButcherTable<mpfr_float>(s_, aRKV67_, bRKV67_, b_minus_bstarRKV67_,
|
|
cRKV67_);
|
|
|
|
break;
|
|
}
|
|
|
|
default: {
|
|
throw std::runtime_error(
|
|
"incompatible predictor choice in ExplicitPredict");
|
|
}
|
|
}
|
|
ResizeK();
|
|
}; // re: PredictorMethod
|
|
|
|
/**
|
|
\brief Change the system (number of total functions) that the predictor uses.
|
|
|
|
\param S New system to switch to.
|
|
*/
|
|
void ChangeSystem(const System& S) {
|
|
numTotalFunctions_ = S.NumTotalFunctions();
|
|
numVariables_ = S.NumVariables();
|
|
// you cannot set K_ here, because s_ may not have been set
|
|
std::get<Mat<dbl>>(dh_dx_0_).resize(numTotalFunctions_, numVariables_);
|
|
std::get<Mat<mpfr_complex>>(dh_dx_0_).resize(numTotalFunctions_,
|
|
numVariables_);
|
|
std::get<Mat<dbl>>(dh_dx_temp_).resize(numTotalFunctions_, numVariables_);
|
|
std::get<Mat<mpfr_complex>>(dh_dx_temp_)
|
|
.resize(numTotalFunctions_, numVariables_);
|
|
std::get<Vec<dbl>>(dh_dt_temp_).resize(numTotalFunctions_);
|
|
std::get<Vec<mpfr_complex>>(dh_dt_temp_).resize(numTotalFunctions_);
|
|
|
|
ResizeK();
|
|
}
|
|
|
|
void ResizeK() {
|
|
std::get<Mat<dbl>>(K_).resize(numTotalFunctions_, s_);
|
|
std::get<Mat<mpfr_complex>>(K_).resize(numTotalFunctions_, s_);
|
|
}
|
|
|
|
/**
|
|
\brief get the current precision of the predictor
|
|
*/
|
|
unsigned precision() const { return current_precision_; }
|
|
|
|
/**
|
|
/brief Change the precision of the predictor variables and reassign the
|
|
Butcher table variables.
|
|
|
|
\param new_precision The new precision.
|
|
|
|
*/
|
|
void ChangePrecision(unsigned new_precision) {
|
|
Precision(std::get<Mat<mpfr_complex>>(K_), new_precision);
|
|
|
|
Precision(std::get<Vec<mpfr_complex>>(dh_dt_temp_), new_precision);
|
|
Precision(std::get<Mat<mpfr_complex>>(dh_dx_0_), new_precision);
|
|
Precision(std::get<Mat<mpfr_complex>>(dh_dx_temp_), new_precision);
|
|
|
|
Precision(std::get<Mat<mpfr_float>>(a_), new_precision);
|
|
Precision(std::get<Vec<mpfr_float>>(b_), new_precision);
|
|
Precision(std::get<Vec<mpfr_float>>(b_minus_bstar_), new_precision);
|
|
Precision(std::get<Vec<mpfr_float>>(c_), new_precision);
|
|
|
|
PredictorMethod(predictor_);
|
|
|
|
current_precision_ = new_precision;
|
|
|
|
PrecisionSanityCheck();
|
|
}
|
|
|
|
void PrecisionSanityCheck() const {
|
|
assert(current_precision_ == DefaultPrecision());
|
|
|
|
Vec<mpfr_complex>& dhdttemp = std::get<Vec<mpfr_complex>>(dh_dt_temp_);
|
|
Mat<mpfr_complex>& dhdx0 = std::get<Mat<mpfr_complex>>(dh_dx_0_);
|
|
Mat<mpfr_complex>& dhdxtemp = std::get<Mat<mpfr_complex>>(dh_dx_temp_);
|
|
|
|
Mat<mpfr_float>& a = std::get<Mat<mpfr_float>>(a_);
|
|
Vec<mpfr_float>& b = std::get<Vec<mpfr_float>>(b_);
|
|
Vec<mpfr_float>& bstar = std::get<Vec<mpfr_float>>(b_minus_bstar_);
|
|
Vec<mpfr_float>& c = std::get<Vec<mpfr_float>>(c_);
|
|
|
|
assert(Precision(dhdttemp) == current_precision_);
|
|
assert(Precision(dhdx0) == current_precision_);
|
|
assert(Precision(dhdxtemp) == current_precision_);
|
|
|
|
assert(Precision(a) == current_precision_);
|
|
assert(Precision(b) == current_precision_);
|
|
if (uses_embedded_) assert(Precision(bstar) == current_precision_);
|
|
assert(Precision(c) == current_precision_);
|
|
}
|
|
|
|
/**
|
|
\brief Perform a generic predictor step.
|
|
|
|
\param next_space The computed prediction.
|
|
\param method An enum class selecting the predictor method to use.
|
|
\param S The system being solved.
|
|
\param current_space The current space variable vector.
|
|
\param current_time The current time.
|
|
\param delta_t The size of the time step.
|
|
\param condition_number_estimate The computed estimate of the condition
|
|
number of the Jacobian. \param
|
|
num_steps_since_last_condition_number_computation. Updated in this function.
|
|
\param frequency_of_CN_estimation How many steps to take between condition
|
|
number estimates. \param prec_type The operating precision type. \param
|
|
tracking_tolerance How tightly to track the path.
|
|
|
|
\return SuccessCode indicating how the prediction went.
|
|
*/
|
|
|
|
template <typename ComplexType>
|
|
SuccessCode Predict(
|
|
Vec<ComplexType>& next_space, System const& S,
|
|
const Vec<ComplexType>& current_space, ComplexType current_time,
|
|
ComplexType const& delta_t, NumErrorT& condition_number_estimate,
|
|
unsigned& num_steps_since_last_condition_number_computation,
|
|
unsigned frequency_of_CN_estimation,
|
|
NumErrorT const& tracking_tolerance) {
|
|
auto step_success =
|
|
FullStep(next_space, S, current_space, current_time, delta_t);
|
|
|
|
NumErrorT norm_J, norm_J_inverse;
|
|
SetNormsCond<ComplexType>(norm_J, norm_J_inverse, condition_number_estimate,
|
|
num_steps_since_last_condition_number_computation,
|
|
frequency_of_CN_estimation);
|
|
|
|
return step_success;
|
|
}
|
|
|
|
/**
|
|
\brief Perform a generic predictor step and return size_proportion and
|
|
condition number information
|
|
|
|
\param next_space The computed prediction.
|
|
\param method An enum class selecting the predictor method to use.
|
|
\param size_proportion $a$ in AMP2 paper.
|
|
\param norm_J The computed estimate of the norm of the Jacobian matrix.
|
|
\param norm_J_inverse The computed estimate of the norm of the inverse of the
|
|
Jacobian matrix. \param S The system being solved. \param current_space The
|
|
current space variable vector. \param current_time The current time. \param
|
|
delta_t The size of the time step. \param condition_number_estimate The
|
|
computed estimate of the condition number of the Jacobian. \param
|
|
num_steps_since_last_condition_number_computation. Updated in this function.
|
|
\param frequency_of_CN_estimation How many steps to take between condition
|
|
number estimates. \param prec_type The operating precision type. \param
|
|
tracking_tolerance How tightly to track the path. \param AMP_config The
|
|
settings for adaptive multiple precision.
|
|
|
|
\return SuccessCode indicating how the prediction went.
|
|
*/
|
|
|
|
template <typename ComplexType>
|
|
SuccessCode Predict(
|
|
Vec<ComplexType>& next_space, NumErrorT& size_proportion,
|
|
NumErrorT& norm_J, NumErrorT& norm_J_inverse, System const& S,
|
|
const Vec<ComplexType>& current_space, ComplexType current_time,
|
|
ComplexType const& delta_t, NumErrorT& condition_number_estimate,
|
|
unsigned& num_steps_since_last_condition_number_computation,
|
|
unsigned frequency_of_CN_estimation, NumErrorT const& tracking_tolerance,
|
|
AdaptiveMultiplePrecisionConfig const& AMP_config) {
|
|
auto success_code =
|
|
Predict<ComplexType>(next_space, S, current_space, current_time,
|
|
delta_t, condition_number_estimate,
|
|
num_steps_since_last_condition_number_computation,
|
|
frequency_of_CN_estimation, tracking_tolerance);
|
|
|
|
SetNormsCond<ComplexType>(norm_J, norm_J_inverse, condition_number_estimate,
|
|
num_steps_since_last_condition_number_computation,
|
|
frequency_of_CN_estimation);
|
|
|
|
// Set size_proportion
|
|
SetSizeProportion(size_proportion, delta_t);
|
|
|
|
if (success_code != SuccessCode::Success) return success_code;
|
|
|
|
// AMP Criteria
|
|
if (!amp::CriterionA<ComplexType>(norm_J, norm_J_inverse,
|
|
AMP_config)) // AMP_criterion_A != ok
|
|
{
|
|
return SuccessCode::HigherPrecisionNecessary;
|
|
} else if (!amp::CriterionC<ComplexType>(
|
|
norm_J_inverse, current_space, tracking_tolerance,
|
|
AMP_config)) // AMP_criterion_C != ok
|
|
{
|
|
return SuccessCode::HigherPrecisionNecessary;
|
|
}
|
|
|
|
return success_code;
|
|
}
|
|
|
|
/**
|
|
\brief Perform a generic predictor step and return error estimate,
|
|
size_proportion and condition number information
|
|
|
|
\param next_space The computed prediction.
|
|
\param method An enum class selecting the predictor method to use.
|
|
\param error_estimate Estimate of the error from an embedded method.
|
|
\param size_proportion $a$ in AMP2 paper.
|
|
\param norm_J The computed estimate of the norm of the Jacobian matrix.
|
|
\param norm_J_inverse The computed estimate of the norm of the inverse of the
|
|
Jacobian matrix. \param S The system being solved. \param current_space The
|
|
current space variable vector. \param current_time The current time. \param
|
|
delta_t The size of the time step. \param condition_number_estimate The
|
|
computed estimate of the condition number of the Jacobian. \param
|
|
num_steps_since_last_condition_number_computation. Updated in this function.
|
|
\param frequency_of_CN_estimation How many steps to take between condition
|
|
number estimates. \param prec_type The operating precision type. \param
|
|
tracking_tolerance How tightly to track the path. \param AMP_config The
|
|
settings for adaptive multiple precision.
|
|
|
|
\return SuccessCode indicating how the prediction went.
|
|
*/
|
|
|
|
template <typename ComplexType>
|
|
SuccessCode Predict(
|
|
Vec<ComplexType>& next_space, NumErrorT& error_estimate,
|
|
NumErrorT& size_proportion, NumErrorT& norm_J, NumErrorT& norm_J_inverse,
|
|
System const& S, const Vec<ComplexType>& current_space,
|
|
ComplexType current_time, ComplexType const& delta_t,
|
|
NumErrorT& condition_number_estimate,
|
|
unsigned& num_steps_since_last_condition_number_computation,
|
|
unsigned frequency_of_CN_estimation, NumErrorT const& tracking_tolerance,
|
|
AdaptiveMultiplePrecisionConfig const& AMP_config) {
|
|
// If this is a method without an error estimator, then can't calculate size
|
|
// proportion and should throw an error
|
|
|
|
if (!predict::HasErrorEstimate(predictor_))
|
|
throw std::runtime_error(
|
|
"incompatible predictor choice in ExplicitPredict, no error "
|
|
"estimator");
|
|
|
|
auto success_code =
|
|
Predict(next_space, size_proportion, norm_J, norm_J_inverse, S,
|
|
current_space, current_time, delta_t, condition_number_estimate,
|
|
num_steps_since_last_condition_number_computation,
|
|
frequency_of_CN_estimation, tracking_tolerance, AMP_config);
|
|
|
|
SetErrorEstimate(error_estimate, delta_t);
|
|
|
|
return success_code;
|
|
}
|
|
|
|
/**
|
|
\brief Get the currently used prediction method.
|
|
|
|
\return The current method.
|
|
*/
|
|
Predictor PredictorMethod() { return predictor_; }
|
|
|
|
/**
|
|
\brief Get the order of the currently used prediction method.
|
|
|
|
This is the lowest order of the predictor. The order of the error estimate is
|
|
this plus one.
|
|
|
|
\return The aforementioned order.
|
|
*/
|
|
inline unsigned Order() { return p_; }
|
|
|
|
/**
|
|
\brief Get whether the current prediction method provides an error estimate.
|
|
|
|
\return Yes or no.
|
|
*/
|
|
inline bool HasErrorEstimate() {
|
|
return predict::HasErrorEstimate(predictor_);
|
|
}
|
|
|
|
protected:
|
|
///////////////////////////
|
|
//
|
|
// Protected Methods
|
|
//
|
|
////////////////////
|
|
|
|
template <typename T>
|
|
Eigen::PartialPivLU<Mat<T>>& GetLU() {
|
|
return LUSelector<T>::Run(*this);
|
|
}
|
|
|
|
Eigen::PartialPivLU<Mat<dbl>>& GetLU_d() { return LU_d_; }
|
|
|
|
Eigen::PartialPivLU<Mat<mpfr_complex>>& GetLU_mp() {
|
|
assert(current_precision_ == DefaultPrecision());
|
|
return LU_mp_[current_precision_];
|
|
}
|
|
|
|
/**
|
|
\brief Performs a full prediction step from current_time to current_time +
|
|
delta_t
|
|
|
|
\param next_space The computed prediction space
|
|
\param S The homotopy system
|
|
\param current_space The current space values
|
|
\param current_time The current time values
|
|
\param delta_t The time step
|
|
|
|
\return SuccessCode determining result of the computation
|
|
*/
|
|
|
|
template <typename ComplexType>
|
|
SuccessCode FullStep(Vec<ComplexType>& next_space, System const& S,
|
|
Vec<ComplexType> const& current_space,
|
|
ComplexType const& current_time,
|
|
ComplexType const& delta_t) {
|
|
// If using constant predictor
|
|
if (s_ == 0) {
|
|
next_space = current_space;
|
|
return SuccessCode::Success;
|
|
}
|
|
|
|
using RealType = typename Eigen::NumTraits<ComplexType>::Real;
|
|
|
|
Mat<ComplexType>& Kref = std::get<Mat<ComplexType>>(K_);
|
|
Mat<RealType>& aref = std::get<Mat<RealType>>(a_);
|
|
Vec<RealType>& bref = std::get<Vec<RealType>>(b_);
|
|
Vec<RealType>& cref = std::get<Vec<RealType>>(c_);
|
|
Kref.fill(ComplexType(0));
|
|
Vec<ComplexType> temp(S.NumTotalFunctions());
|
|
|
|
if (EvalRHS(S, current_space, current_time, Kref, 0) !=
|
|
SuccessCode::Success) {
|
|
return SuccessCode::MatrixSolveFailureFirstPartOfPrediction;
|
|
}
|
|
|
|
for (int ii = 1; ii < s_; ++ii) {
|
|
temp.setZero(); // see https://github.com/bertiniteam/b2/issues/198
|
|
for (int jj = 0; jj < ii; ++jj) temp += aref(ii, jj) * Kref.col(jj);
|
|
|
|
// Vec<ComplexType> wfp =
|
|
if (EvalRHS<ComplexType>(S, current_space + delta_t * temp,
|
|
current_time + cref(ii) * delta_t, Kref,
|
|
ii) != SuccessCode::Success)
|
|
return SuccessCode::MatrixSolveFailure;
|
|
}
|
|
|
|
temp.setZero();
|
|
for (int ii = 0; ii < s_; ++ii) temp += bref(ii) * Kref.col(ii);
|
|
|
|
next_space = current_space + delta_t * temp;
|
|
|
|
return SuccessCode::Success;
|
|
};
|
|
|
|
template <typename ComplexType>
|
|
void SetNormsCond(NumErrorT& norm_J, NumErrorT& norm_J_inverse,
|
|
NumErrorT& condition_number_estimate,
|
|
unsigned num_steps_since_last_condition_number_computation,
|
|
unsigned frequency_of_CN_estimation) {
|
|
// Calculate condition number and update if needed
|
|
Eigen::PartialPivLU<Mat<ComplexType>>& LUref = GetLU<ComplexType>();
|
|
Mat<ComplexType>& dhdxref = std::get<Mat<ComplexType>>(dh_dx_0_);
|
|
|
|
// TODO this random vector should not be made fresh every time. especiallyif
|
|
// the numeric type is mpfr_complex!
|
|
Vec<ComplexType> randy = RandomOfUnits<ComplexType>(numVariables_);
|
|
Vec<ComplexType> temp_soln = LUref.solve(randy);
|
|
|
|
norm_J = NumErrorT(dhdxref.norm());
|
|
norm_J_inverse = NumErrorT(temp_soln.norm());
|
|
|
|
if (num_steps_since_last_condition_number_computation >=
|
|
frequency_of_CN_estimation) {
|
|
condition_number_estimate = NumErrorT(norm_J * norm_J_inverse);
|
|
num_steps_since_last_condition_number_computation =
|
|
1; // reset the counter to 1
|
|
} else // no need to compute the condition number
|
|
num_steps_since_last_condition_number_computation++;
|
|
}
|
|
|
|
/**
|
|
\brief Computes the error estimate of this prediction step.
|
|
|
|
\param error_estimate Computed error estimate
|
|
\param delta_t The time step
|
|
|
|
\return Success code or the computation
|
|
|
|
*/
|
|
|
|
template <typename ComplexType>
|
|
SuccessCode SetErrorEstimate(NumErrorT& error_estimate,
|
|
ComplexType const& delta_t) {
|
|
using RealType = typename Eigen::NumTraits<ComplexType>::Real;
|
|
|
|
Mat<ComplexType>& Kref = std::get<Mat<ComplexType>>(K_);
|
|
Vec<RealType>& b_minus_bstar_ref = std::get<Vec<RealType>>(b_minus_bstar_);
|
|
|
|
auto numFuncs = Kref.rows();
|
|
Vec<ComplexType> err(numFuncs);
|
|
|
|
err.setZero();
|
|
for (int ii = 0; ii < s_; ++ii) {
|
|
err += (b_minus_bstar_ref(ii)) * Kref.col(ii);
|
|
}
|
|
|
|
err *= delta_t;
|
|
|
|
error_estimate = NumErrorT(err.norm());
|
|
|
|
return SuccessCode::Success;
|
|
};
|
|
|
|
/**
|
|
\brief Compute the size proportion variable for AMP computation
|
|
|
|
\param size_proportion Computed size proportion
|
|
\param delta_t The time step
|
|
|
|
\return Success code of the computation
|
|
|
|
*/
|
|
|
|
template <typename ComplexType>
|
|
SuccessCode SetSizeProportion(NumErrorT& size_proportion,
|
|
ComplexType const& delta_t) {
|
|
if (predict::HasErrorEstimate(predictor_)) {
|
|
NumErrorT err_est;
|
|
SetErrorEstimate(err_est, delta_t);
|
|
|
|
using std::pow;
|
|
size_proportion = err_est / NumErrorT(pow(abs(delta_t), p_ + 1));
|
|
|
|
return SuccessCode::Success;
|
|
} else {
|
|
Mat<ComplexType>& Kref = std::get<Mat<ComplexType>>(K_);
|
|
using std::pow;
|
|
size_proportion =
|
|
NumErrorT(Kref.array().abs().maxCoeff() / (pow(abs(delta_t), p_)));
|
|
return SuccessCode::Success;
|
|
}
|
|
};
|
|
|
|
/**
|
|
\brief Evaluates the RHS of the Davidenko differential equation at a
|
|
particular time and space
|
|
|
|
\param S The homotopy system
|
|
\param space The space variable used to evaluate RHS
|
|
\param time The time variable used to evaluate RHS
|
|
\param K Matrix of stage variables
|
|
\param stage Which stage variable(column of K) should be filled by this
|
|
computation
|
|
|
|
\return Success code of this computation
|
|
*/
|
|
|
|
template <typename ComplexType>
|
|
SuccessCode EvalRHS(System const& S, const Vec<ComplexType>& space,
|
|
const ComplexType& time, Mat<ComplexType>& K,
|
|
unsigned stage) {
|
|
if (std::is_same<ComplexType, mpfr_complex>::value) PrecisionSanityCheck();
|
|
|
|
if (stage == 0) {
|
|
Eigen::PartialPivLU<Mat<ComplexType>>& LUref = GetLU<ComplexType>();
|
|
Mat<ComplexType>& dhdxref = std::get<Mat<ComplexType>>(dh_dx_0_);
|
|
|
|
if (!std::is_same<ComplexType, dbl>::value) {
|
|
assert(DefaultPrecision() == current_precision_);
|
|
|
|
assert(Precision(space) == current_precision_);
|
|
assert(Precision(time) == current_precision_);
|
|
assert(Precision(dhdxref) == current_precision_);
|
|
assert(Precision(K) == current_precision_);
|
|
}
|
|
S.SetAndReset<ComplexType>(space, time);
|
|
S.JacobianInPlace(dhdxref);
|
|
LUref = dhdxref.lu();
|
|
if (!std::is_same<ComplexType, dbl>::value) {
|
|
assert(Precision(dhdxref) == current_precision_);
|
|
assert(Precision(LUref.matrixLU()) == current_precision_);
|
|
}
|
|
|
|
if (LUPartialPivotDecompositionSuccessful(LUref.matrixLU()) !=
|
|
MatrixSuccessCode::Success)
|
|
return SuccessCode::MatrixSolveFailureFirstPartOfPrediction;
|
|
|
|
Vec<ComplexType>& dhdtref = std::get<Vec<ComplexType>>(dh_dt_temp_);
|
|
S.TimeDerivativeInPlace(dhdtref);
|
|
K.col(stage) = LUref.solve(-dhdtref);
|
|
|
|
return SuccessCode::Success;
|
|
|
|
} else {
|
|
S.SetAndReset<ComplexType>(space, time);
|
|
|
|
Mat<ComplexType>& dhdxtempref = std::get<Mat<ComplexType>>(dh_dx_temp_);
|
|
S.JacobianInPlace(dhdxtempref);
|
|
auto LU = dhdxtempref.lu();
|
|
|
|
if (LUPartialPivotDecompositionSuccessful(LU.matrixLU()) !=
|
|
MatrixSuccessCode::Success)
|
|
return SuccessCode::MatrixSolveFailure;
|
|
|
|
Vec<ComplexType>& dhdtref = std::get<Vec<ComplexType>>(dh_dt_temp_);
|
|
S.TimeDerivativeInPlace(dhdtref);
|
|
K.col(stage) = LU.solve(-dhdtref);
|
|
|
|
return SuccessCode::Success;
|
|
}
|
|
}
|
|
|
|
private:
|
|
///////////////////////////
|
|
//
|
|
// Private Methods
|
|
//
|
|
////////////////////
|
|
|
|
/**
|
|
/brief Fills the local embedded butcher table variables a,b,bstar and c with
|
|
the constant static values stored in the class.
|
|
|
|
\param stages Number of stages. Used to create correct size on variables
|
|
\param a Matrix of the Butcher table
|
|
\param b Weights in the Butcher table
|
|
\param bstar Weights of embedded method in Butcher table
|
|
\param c Time offsets in Butcher table
|
|
|
|
*/
|
|
|
|
template <typename RealType>
|
|
void FillButcherTable(int stages, const Mat<mpq_rational>& a,
|
|
const Mat<mpq_rational>& b,
|
|
const Mat<mpq_rational>& b_minus_bstar,
|
|
const Mat<mpq_rational>& c) {
|
|
Mat<RealType>& aref = std::get<Mat<RealType>>(a_);
|
|
aref.resize(stages, stages);
|
|
for (int ii = 0; ii < stages; ++ii) {
|
|
for (int jj = 0; jj < s_; ++jj) {
|
|
aref(ii, jj) = static_cast<RealType>(a(ii, jj));
|
|
}
|
|
}
|
|
|
|
Vec<RealType>& bref = std::get<Vec<RealType>>(b_);
|
|
bref.resize(stages);
|
|
for (int ii = 0; ii < stages; ++ii) {
|
|
bref(ii) = static_cast<RealType>(b(ii));
|
|
}
|
|
|
|
Vec<RealType>& b_minus_bstar_ref = std::get<Vec<RealType>>(b_minus_bstar_);
|
|
b_minus_bstar_ref.resize(stages);
|
|
for (int ii = 0; ii < stages; ++ii) {
|
|
b_minus_bstar_ref(ii) = static_cast<RealType>(b_minus_bstar(ii));
|
|
}
|
|
|
|
Vec<RealType>& cref = std::get<Vec<RealType>>(c_);
|
|
cref.resize(stages);
|
|
for (int ii = 0; ii < stages; ++ii) {
|
|
cref(ii) = static_cast<RealType>(c(ii));
|
|
}
|
|
uses_embedded_ = true;
|
|
}
|
|
|
|
/**
|
|
/brief Fills the local butcher table variables a,b and c with the constant
|
|
static values stored in the class.
|
|
|
|
\param stages Number of stages. Used to create correct size on variables
|
|
\param a Matrix of the Butcher table
|
|
\param b Weights in the Butcher table
|
|
\param c Time offsets in Butcher table
|
|
|
|
*/
|
|
|
|
template <typename RealType>
|
|
void FillButcherTable(int stages, const Mat<mpq_rational>& a,
|
|
const Mat<mpq_rational>& b,
|
|
const Mat<mpq_rational>& c) {
|
|
Mat<RealType>& aref = std::get<Mat<RealType>>(a_);
|
|
aref.resize(stages, stages);
|
|
for (int ii = 0; ii < stages; ++ii) {
|
|
for (int jj = 0; jj < s_; ++jj) {
|
|
aref(ii, jj) = static_cast<RealType>(a(ii, jj));
|
|
}
|
|
}
|
|
|
|
Vec<RealType>& bref = std::get<Vec<RealType>>(b_);
|
|
bref.resize(stages);
|
|
for (int ii = 0; ii < stages; ++ii) {
|
|
bref(ii) = static_cast<RealType>(b(ii));
|
|
}
|
|
|
|
Vec<RealType>& cref = std::get<Vec<RealType>>(c_);
|
|
cref.resize(stages);
|
|
for (int ii = 0; ii < stages; ++ii) {
|
|
cref(ii) = static_cast<RealType>(c(ii));
|
|
}
|
|
uses_embedded_ = false;
|
|
}
|
|
|
|
///////////////////////////
|
|
//
|
|
// Private Data Members
|
|
//
|
|
////////////////////
|
|
|
|
unsigned
|
|
numTotalFunctions_; // Number of total functions for the current system
|
|
unsigned numVariables_; // Number of variables for the current system
|
|
mutable std::tuple<Mat<dbl>, Mat<mpfr_complex>>
|
|
K_; // All the stage variables. Each column represents a different
|
|
// stage.
|
|
Predictor predictor_; // Method for prediction
|
|
unsigned p_; // Order of the prediction method
|
|
mutable std::tuple<Mat<dbl>, Mat<mpfr_complex>>
|
|
dh_dx_0_; // Jacobian for the initial stage. Use for AMP testing
|
|
mutable std::tuple<Mat<dbl>, Mat<mpfr_complex>>
|
|
dh_dx_temp_; // Temporary jacobian for all other stages
|
|
mutable std::tuple<Vec<dbl>, Vec<mpfr_complex>>
|
|
dh_dt_temp_; // Temporary time derivative used for all stages
|
|
// std::tuple< Eigen::PartialPivLU<Mat<dbl>>,
|
|
// Eigen::PartialPivLU<Mat<mpfr_complex>> > LU_0_; // LU from the intial
|
|
// stage used for AMP testing
|
|
|
|
mutable Eigen::PartialPivLU<Mat<dbl>> LU_d_;
|
|
mutable std::map<unsigned, Eigen::PartialPivLU<Mat<mpfr_complex>>> LU_mp_;
|
|
|
|
// Butcher Table (notation from
|
|
// https://en.wikipedia.org/wiki/List_of_Runge%E2%80%93Kutta_methods )
|
|
mutable unsigned s_; // Number of stages
|
|
mutable std::tuple<Mat<double>, Mat<mpfr_float>> a_;
|
|
mutable std::tuple<Vec<double>, Vec<mpfr_float>> b_;
|
|
mutable std::tuple<Vec<double>, Vec<mpfr_float>> b_minus_bstar_;
|
|
mutable std::tuple<Vec<double>, Vec<mpfr_float>> c_;
|
|
|
|
mutable bool uses_embedded_;
|
|
mutable unsigned current_precision_;
|
|
|
|
// static const variables that store the Butcher table in mpq_rational form
|
|
// Euler
|
|
static const mpq_rational aEulerPtr_[];
|
|
static const Eigen::Matrix<mpq_rational, 1, 1> aEuler_;
|
|
static const mpq_rational bEulerPtr_[];
|
|
static const Eigen::Matrix<mpq_rational, 1, 1> bEuler_;
|
|
static const mpq_rational cEulerPtr_[];
|
|
static const Eigen::Matrix<mpq_rational, 1, 1> cEuler_;
|
|
|
|
// Heun-Euler
|
|
static const mpq_rational aHeunEulerPtr_[];
|
|
static const Eigen::Matrix<mpq_rational, 2, 2> aHeunEuler_;
|
|
static const mpq_rational bHeunEulerPtr_[];
|
|
static const Eigen::Matrix<mpq_rational, 2, 1> bHeunEuler_;
|
|
static const mpq_rational b_minus_bstarHeunEulerPtr_[];
|
|
static const Eigen::Matrix<mpq_rational, 2, 1> b_minus_bstarHeunEuler_;
|
|
static const mpq_rational cHeunEulerPtr_[];
|
|
static const Eigen::Matrix<mpq_rational, 2, 1> cHeunEuler_;
|
|
|
|
// RK4
|
|
static const mpq_rational aRK4Ptr_[];
|
|
static const Eigen::Matrix<mpq_rational, 4, 4> aRK4_;
|
|
static const mpq_rational bRK4Ptr_[];
|
|
static const Eigen::Matrix<mpq_rational, 4, 1> bRK4_;
|
|
static const mpq_rational cRK4Ptr_[];
|
|
static const Eigen::Matrix<mpq_rational, 4, 1> cRK4_;
|
|
|
|
// RKF45
|
|
static const mpq_rational aRKF45Ptr_[];
|
|
static const Eigen::Matrix<mpq_rational, 6, 6> aRKF45_;
|
|
static const mpq_rational bRKF45Ptr_[];
|
|
static const Eigen::Matrix<mpq_rational, 6, 1> bRKF45_;
|
|
static const mpq_rational b_minus_bstarRKF45Ptr_[];
|
|
static const Eigen::Matrix<mpq_rational, 6, 1> b_minus_bstarRKF45_;
|
|
static const mpq_rational cRKF45Ptr_[];
|
|
static const Eigen::Matrix<mpq_rational, 6, 1> cRKF45_;
|
|
|
|
// RK Cash-Karp45
|
|
static const mpq_rational aRKCK45Ptr_[];
|
|
static const Eigen::Matrix<mpq_rational, 6, 6> aRKCK45_;
|
|
static const mpq_rational bRKCK45Ptr_[];
|
|
static const Eigen::Matrix<mpq_rational, 6, 1> bRKCK45_;
|
|
static const mpq_rational b_minus_bstarRKCK45Ptr_[];
|
|
static const Eigen::Matrix<mpq_rational, 6, 1> b_minus_bstarRKCK45_;
|
|
static const mpq_rational cRKCK45Ptr_[];
|
|
static const Eigen::Matrix<mpq_rational, 6, 1> cRKCK45_;
|
|
|
|
// RK Dormand-Prince 56
|
|
static const mpq_rational aRKDP56Ptr_[];
|
|
static const Eigen::Matrix<mpq_rational, 8, 8> aRKDP56_;
|
|
static const mpq_rational bRKDP56Ptr_[];
|
|
static const Eigen::Matrix<mpq_rational, 8, 1> bRKDP56_;
|
|
static const mpq_rational b_minus_bstarRKDP56Ptr_[];
|
|
static const Eigen::Matrix<mpq_rational, 8, 1> b_minus_bstarRKDP56_;
|
|
static const mpq_rational cRKDP56Ptr_[];
|
|
static const Eigen::Matrix<mpq_rational, 8, 1> cRKDP56_;
|
|
|
|
// RK Verner 67
|
|
static const mpq_rational aRKV67Ptr_[];
|
|
static const Eigen::Matrix<mpq_rational, 10, 10> aRKV67_;
|
|
static const mpq_rational bRKV67Ptr_[];
|
|
static const Eigen::Matrix<mpq_rational, 10, 1> bRKV67_;
|
|
static const mpq_rational b_minus_bstarRKV67Ptr_[];
|
|
static const Eigen::Matrix<mpq_rational, 10, 1> b_minus_bstarRKV67_;
|
|
static const mpq_rational cRKV67Ptr_[];
|
|
static const Eigen::Matrix<mpq_rational, 10, 1> cRKV67_;
|
|
|
|
}; // re: ExplicitRKPredictor class
|
|
|
|
} // namespace predict
|
|
} // namespace tracking
|
|
} // namespace bertini
|
|
|
|
#endif
|
|
|