Bertini 2 - C++ Core
Software for Numerical Algebraic Geometry
Loading...
Searching...
No Matches
explicit_predictors.hpp
Go to the documentation of this file.
1//This file is part of Bertini 2.0.
2//
3//bertini2/trackers/explicit_predictors.hpp is free software: you can redistribute it and/or modify
4//it under the terms of the GNU General Public License as published by
5//the Free Software Foundation, either version 3 of the License, or
6//(at your option) any later version.
7//
8//bertini2/trackers/explicit_predictors.hpp is distributed in the hope that it will be useful,
9//but WITHOUT ANY WARRANTY; without even the implied warranty of
10//MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11//GNU General Public License for more details.
12//
13//You should have received a copy of the GNU General Public License
14//along with bertini2/trackers/explicit_predictors.hpp. If not, see <http://www.gnu.org/licenses/>.
15//
16// copyright 2015
17// James B. Collins
18// West Texas A&M University
19// Department of Mathematics
20// Spring 2016
21//
22// silviana amethyst, university of wisconsin-eau claire
23
24
25
31
32#ifndef BERTINI_EXPLICIT_PREDICTORS_HPP
33#define BERTINI_EXPLICIT_PREDICTORS_HPP
34
35#include "bertini2/trackers/amp_criteria.hpp"
36
37#include "bertini2/system/system.hpp"
38#include "bertini2/mpfr_extensions.hpp"
39#include <Eigen/LU>
40
41#include <boost/type_index.hpp>
42
43#include "bertini2/eigen_extensions.hpp"
44
45namespace bertini{
46 namespace tracking{
47 namespace predict{
48
49
50 // Constant,
51 // Euler,
52 // Heun,
53 // RK4,
54 // HeunEuler,
55 // RKNorsett34,
56 // RKF45,
57 // RKCashKarp45,
58 // RKDormandPrince56,
59 // RKVerner67
60
68 inline
69 Predictor DefaultPredictor()
70 {
71 return Predictor::RKF45;
72 }
73
74
83 inline
84 unsigned Order(Predictor predictor_choice)
85 {
86 switch (predictor_choice)
87 {
88 case (Predictor::Constant):
89 return 0;
90 case (Predictor::Euler):
91 return 1;
92 case (Predictor::Heun):
93 return 2;
94 case (Predictor::HeunEuler):
95 return 2;
96 case (Predictor::RKNorsett34):
97 return 3;
98 case (Predictor::RK4):
99 return 4;
100 case (Predictor::RKF45):
101 return 4;
102 case (Predictor::RKCashKarp45):
103 return 4;
104 case (Predictor::RKDormandPrince56):
105 return 5;
106 case (Predictor::RKVerner67):
107 return 6;
108 default:
109 {
110 throw std::runtime_error("incompatible predictor choice in Order");
111 }
112 }
113 }
114
121 inline bool HasErrorEstimate(Predictor predictor_choice)
122 {
123 switch (predictor_choice)
124 {
125 case (Predictor::Constant):
126 return false;
127 case (Predictor::Euler):
128 return false;
129 case (Predictor::Heun):
130 return false;
131 case (Predictor::HeunEuler):
132 return true;
133 case (Predictor::RKNorsett34):
134 return true;
135 case (Predictor::RK4):
136 return false;
137 case (Predictor::RKF45):
138 return true;
139 case (Predictor::RKCashKarp45):
140 return true;
141 case (Predictor::RKDormandPrince56):
142 return true;
143 case (Predictor::RKVerner67):
144 return true;
145 default:
146 {
147 throw std::runtime_error("incompatible predictor choice in HasErrorEstimate");
148 }
149 }
150 }
151
152
153 namespace {
154 template<typename T>
155 struct LUSelector
156 {};
157
158 template<>
159 struct LUSelector<dbl>
160 {
161 template<typename N>
162 static Eigen::PartialPivLU<Mat<dbl>>& Run(N & n)
163 {
164 return n.GetLU_d();
165 }
166 };
167
168 template<>
169 struct LUSelector<mpfr_complex>
170 {
171 template<typename N>
172 static Eigen::PartialPivLU<Mat<mpfr_complex>>& Run(N & n)
173 {
174 return n.GetLU_mp();
175 }
176 };
177 }
178
201 class ExplicitRKPredictor
202 {
203 friend LUSelector<dbl>;
204 friend LUSelector<mpfr_complex>;
205 public:
206
212 ExplicitRKPredictor(const System& S) : current_precision_(DefaultPrecision()), s_(0)
213 {
214 ChangeSystem(S);
215 PredictorMethod(DefaultPredictor());
216 }
217
224 ExplicitRKPredictor(Predictor method, const System& S) : current_precision_(DefaultPrecision()), s_(0)
225 {
226 ChangeSystem(S);
227 PredictorMethod(method);
228 }
229
230
231
232
238 void PredictorMethod(Predictor method)
239 {
240 predictor_ = method;
241 p_ = predict::Order(method);
242 switch(method)
243 {
244 case Predictor::Constant:
245 {
246 s_ = 1;
247 Mat<double>& arefd = std::get< Mat<double> >(a_);
248 Vec<double>& brefd = std::get< Vec<double> >(b_);
249 Vec<double>& crefd = std::get< Vec<double> >(c_);
250 crefd.resize(s_); crefd(0) = 0;
251 arefd.resize(s_,s_); arefd(0,0) = 0;
252 brefd.resize(s_); brefd(0) = 0;
253 Mat<mpfr_float>& arefmp = std::get< Mat<mpfr_float> >(a_);
254 Vec<mpfr_float>& brefmp = std::get< Vec<mpfr_float> >(b_);
255 Vec<mpfr_float>& crefmp = std::get< Vec<mpfr_float> >(c_);
256 crefmp.resize(s_); crefmp(0) = 0;
257 arefmp.resize(s_,s_); arefmp(0,0) = 0;
258 brefmp.resize(s_); brefmp(0) = 0;
259 uses_embedded_ = false;
260
261 break;
262 }
263 case Predictor::Euler:
264 {
265 s_ = 1;
266 Mat<double>& arefd = std::get< Mat<double> >(a_);
267 Vec<double>& brefd = std::get< Vec<double> >(b_);
268 Vec<double>& crefd = std::get< Vec<double> >(c_);
269 crefd.resize(s_); crefd(0) = static_cast<double>(cEuler_(0));
270 arefd.resize(s_,s_); arefd(0,0) = static_cast<double>(aEuler_(0,0));
271 brefd.resize(s_); brefd(0) = static_cast<double>(bEuler_(0));
272 Mat<mpfr_float>& arefmp = std::get< Mat<mpfr_float> >(a_);
273 Vec<mpfr_float>& brefmp = std::get< Vec<mpfr_float> >(b_);
274 Vec<mpfr_float>& crefmp = std::get< Vec<mpfr_float> >(c_);
275 crefmp.resize(s_); crefmp(0) = static_cast<mpfr_float>(cEuler_(0));
276 arefmp.resize(s_,s_); arefmp(0,0) = static_cast<mpfr_float>(aEuler_(0,0));
277 brefmp.resize(s_); brefmp(0) = static_cast<mpfr_float>(bEuler_(0));
278 uses_embedded_ = false;
279 break;
280 }
281 case Predictor::HeunEuler:
282 {
283 s_ = 2;
284
285 FillButcherTable<double>(s_, aHeunEuler_, bHeunEuler_, b_minus_bstarHeunEuler_, cHeunEuler_);
286 FillButcherTable<mpfr_float>(s_, aHeunEuler_, bHeunEuler_, b_minus_bstarHeunEuler_, cHeunEuler_);
287
288 break;
289 }
290 case Predictor::RK4:
291 {
292 s_ = 4;
293
294 FillButcherTable<double>(s_, aRK4_, bRK4_, cRK4_);
295 FillButcherTable<mpfr_float>(s_, aRK4_, bRK4_, cRK4_);
296
297 break;
298 }
299
300 case Predictor::RKF45:
301 {
302 s_ = 6;
303
304 FillButcherTable<double>(s_, aRKF45_, bRKF45_, b_minus_bstarRKF45_, cRKF45_);
305 FillButcherTable<mpfr_float>(s_, aRKF45_, bRKF45_, b_minus_bstarRKF45_, cRKF45_);
306
307 break;
308 }
309
310 case Predictor::RKCashKarp45:
311 {
312 s_ = 6;
313
314 FillButcherTable<double>(s_, aRKCK45_, bRKCK45_, b_minus_bstarRKCK45_, cRKCK45_);
315 FillButcherTable<mpfr_float>(s_, aRKCK45_, bRKCK45_, b_minus_bstarRKCK45_, cRKCK45_);
316
317 break;
318 }
319
320 case Predictor::RKDormandPrince56:
321 {
322 s_ = 8;
323
324 FillButcherTable<double>(s_, aRKDP56_, bRKDP56_, b_minus_bstarRKDP56_, cRKDP56_);
325 FillButcherTable<mpfr_float>(s_, aRKDP56_, bRKDP56_, b_minus_bstarRKDP56_, cRKDP56_);
326
327 break;
328 }
329
330 case Predictor::RKVerner67:
331 {
332 s_ = 10;
333
334 FillButcherTable<double>(s_, aRKV67_, bRKV67_, b_minus_bstarRKV67_, cRKV67_);
335 FillButcherTable<mpfr_float>(s_, aRKV67_, bRKV67_, b_minus_bstarRKV67_, cRKV67_);
336
337 break;
338 }
339
340 default:
341 {
342 throw std::runtime_error("incompatible predictor choice in ExplicitPredict");
343 }
344 }
345 ResizeK();
346 }; // re: PredictorMethod
347
348
349
350
351
357 void ChangeSystem(const System& S)
358 {
359 numTotalFunctions_ = S.NumTotalFunctions();
360 numVariables_ = S.NumVariables();
361 // you cannot set K_ here, because s_ may not have been set
362 std::get< Mat<dbl> >(dh_dx_0_).resize(numTotalFunctions_, numVariables_);
363 std::get< Mat<mpfr_complex> >(dh_dx_0_).resize(numTotalFunctions_, numVariables_);
364 std::get< Mat<dbl> >(dh_dx_temp_).resize(numTotalFunctions_, numVariables_);
365 std::get< Mat<mpfr_complex> >(dh_dx_temp_).resize(numTotalFunctions_, numVariables_);
366 std::get< Vec<dbl> >(dh_dt_temp_).resize(numTotalFunctions_);
367 std::get< Vec<mpfr_complex> >(dh_dt_temp_).resize(numTotalFunctions_);
368
369 ResizeK();
370 }
371
372
373 void ResizeK()
374 {
375 std::get< Mat<dbl> >(K_).resize(numTotalFunctions_, s_);
376 std::get< Mat<mpfr_complex> >(K_).resize(numTotalFunctions_, s_);
377 }
378
379
383 unsigned precision() const
384 {
385 return current_precision_;
386 }
387
394 void ChangePrecision(unsigned new_precision)
395 {
396 Precision(std::get< Mat<mpfr_complex> >(K_),new_precision);
397
398 Precision(std::get< Vec<mpfr_complex> >(dh_dt_temp_),new_precision);
399 Precision(std::get< Mat<mpfr_complex> >(dh_dx_0_),new_precision);
400 Precision(std::get< Mat<mpfr_complex> >(dh_dx_temp_),new_precision);
401
402 Precision(std::get< Mat<mpfr_float> >(a_),new_precision);
403 Precision(std::get< Vec<mpfr_float> >(b_),new_precision);
404 Precision(std::get< Vec<mpfr_float> >(b_minus_bstar_),new_precision);
405 Precision(std::get< Vec<mpfr_float> >(c_),new_precision);
406
407 PredictorMethod(predictor_);
408
409 current_precision_ = new_precision;
410
411 PrecisionSanityCheck();
412 }
413
414 void PrecisionSanityCheck() const
415 {
416 assert(current_precision_==DefaultPrecision());
417
418 Vec<mpfr_complex>& dhdttemp = std::get< Vec<mpfr_complex> >(dh_dt_temp_);
419 Mat<mpfr_complex>& dhdx0 = std::get< Mat<mpfr_complex> >(dh_dx_0_);
420 Mat<mpfr_complex>& dhdxtemp = std::get< Mat<mpfr_complex> >(dh_dx_temp_);
421
422 Mat<mpfr_float>& a = std::get< Mat<mpfr_float> >(a_);
423 Vec<mpfr_float>& b = std::get< Vec<mpfr_float> >(b_);
424 Vec<mpfr_float>& bstar = std::get< Vec<mpfr_float> >(b_minus_bstar_);
425 Vec<mpfr_float>& c = std::get< Vec<mpfr_float> >(c_);
426
427
428
429 assert(Precision(dhdttemp)==current_precision_);
430 assert(Precision(dhdx0)==current_precision_);
431 assert(Precision(dhdxtemp)==current_precision_);
432
433 assert(Precision(a)==current_precision_);
434 assert(Precision(b)==current_precision_);
435 if (uses_embedded_)
436 assert(Precision(bstar)==current_precision_);
437 assert(Precision(c)==current_precision_);
438 }
439
440
458
459 template<typename ComplexType>
460 SuccessCode Predict(Vec<ComplexType> & next_space,
461 System const& S,
462 const Vec<ComplexType>& current_space, ComplexType current_time,
463 ComplexType const& delta_t,
464 NumErrorT & condition_number_estimate,
465 unsigned & num_steps_since_last_condition_number_computation,
466 unsigned frequency_of_CN_estimation,
467 NumErrorT const& tracking_tolerance)
468 {
469
470 auto step_success = FullStep(next_space, S, current_space, current_time, delta_t);
471
472 NumErrorT norm_J, norm_J_inverse;
473 SetNormsCond<ComplexType>(norm_J, norm_J_inverse, condition_number_estimate, num_steps_since_last_condition_number_computation, frequency_of_CN_estimation);
474
475 return step_success;
476 }
477
478
479
501
502 template<typename ComplexType>
503 SuccessCode Predict(Vec<ComplexType> & next_space,
504 NumErrorT & size_proportion,
505 NumErrorT & norm_J,
506 NumErrorT & norm_J_inverse,
507 System const& S,
508 const Vec<ComplexType>& current_space, ComplexType current_time,
509 ComplexType const& delta_t,
510 NumErrorT & condition_number_estimate,
511 unsigned & num_steps_since_last_condition_number_computation,
512 unsigned frequency_of_CN_estimation,
513 NumErrorT const& tracking_tolerance,
514 AdaptiveMultiplePrecisionConfig const& AMP_config)
515 {
516
517
518 auto success_code = Predict<ComplexType>(next_space, S, current_space, current_time, delta_t,
519 condition_number_estimate, num_steps_since_last_condition_number_computation,
520 frequency_of_CN_estimation, tracking_tolerance);
521
522 SetNormsCond<ComplexType>(norm_J, norm_J_inverse, condition_number_estimate, num_steps_since_last_condition_number_computation, frequency_of_CN_estimation);
523
524 // Set size_proportion
525 SetSizeProportion(size_proportion, delta_t);
526
527 if(success_code != SuccessCode::Success)
528 return success_code;
529
530
531
532 //AMP Criteria
533 if (!amp::CriterionA<ComplexType>(norm_J, norm_J_inverse, AMP_config)) // AMP_criterion_A != ok
534 {
535 return SuccessCode::HigherPrecisionNecessary;
536 }
537 else if (!amp::CriterionC<ComplexType>(norm_J_inverse, current_space, tracking_tolerance, AMP_config)) // AMP_criterion_C != ok
538 {
539 return SuccessCode::HigherPrecisionNecessary;
540 }
541
542
543 return success_code;
544 }
545
546
569
570 template<typename ComplexType>
571 SuccessCode Predict(Vec<ComplexType> & next_space,
572 NumErrorT & error_estimate,
573 NumErrorT & size_proportion,
574 NumErrorT & norm_J,
575 NumErrorT & norm_J_inverse,
576 System const& S,
577 const Vec<ComplexType>& current_space, ComplexType current_time,
578 ComplexType const& delta_t,
579 NumErrorT & condition_number_estimate,
580 unsigned & num_steps_since_last_condition_number_computation,
581 unsigned frequency_of_CN_estimation,
582 NumErrorT const& tracking_tolerance,
583 AdaptiveMultiplePrecisionConfig const& AMP_config)
584 {
585
586 // If this is a method without an error estimator, then can't calculate size proportion and should throw an error
587
588 if(!predict::HasErrorEstimate(predictor_))
589 throw std::runtime_error("incompatible predictor choice in ExplicitPredict, no error estimator");
590
591
592 auto success_code = Predict(next_space, size_proportion, norm_J, norm_J_inverse,
593 S, current_space, current_time, delta_t,
594 condition_number_estimate, num_steps_since_last_condition_number_computation,
595 frequency_of_CN_estimation, tracking_tolerance, AMP_config);
596
597 SetErrorEstimate(error_estimate, delta_t);
598
599 return success_code;
600 }
601
602
603
604
605
611 Predictor PredictorMethod()
612 {
613 return predictor_;
614 }
615
616
617
618
626 inline
627 unsigned Order()
628 {
629 return p_;
630 }
631
637 inline bool HasErrorEstimate()
638 {
639 return predict::HasErrorEstimate(predictor_);
640 }
641
642
643
644
645
646
647 protected:
649 //
650 // Protected Methods
651 //
653
654 template <typename T>
655 Eigen::PartialPivLU<Mat<T>>& GetLU()
656 {
657 return LUSelector<T>::Run(*this);
658 }
659
660
661 Eigen::PartialPivLU<Mat<dbl>>& GetLU_d()
662 {
663 return LU_d_;
664 }
665
666 Eigen::PartialPivLU<Mat<mpfr_complex>>& GetLU_mp()
667 {
668 assert(current_precision_==DefaultPrecision());
669 return LU_mp_[current_precision_];
670 }
671
683
684 template<typename ComplexType>
685 SuccessCode FullStep(Vec<ComplexType> & next_space,
686 System const& S,
687 Vec<ComplexType> const& current_space, ComplexType const& current_time,
688 ComplexType const& delta_t)
689 {
690
691 // If using constant predictor
692 if(s_ == 0)
693 {
694 next_space = current_space;
695 return SuccessCode::Success;
696 }
697
698 using RealType = typename Eigen::NumTraits<ComplexType>::Real;
699
700 Mat<ComplexType>& Kref = std::get< Mat<ComplexType> >(K_);
701 Mat<RealType>& aref = std::get< Mat<RealType> >(a_);
702 Vec<RealType>& bref = std::get< Vec<RealType> >(b_);
703 Vec<RealType>& cref = std::get< Vec<RealType> >(c_);
704 Kref.fill(ComplexType(0));
705 Vec<ComplexType> temp(S.NumTotalFunctions());
706
707 if(EvalRHS(S, current_space, current_time, Kref, 0) != SuccessCode::Success)
708 {
709 return SuccessCode::MatrixSolveFailureFirstPartOfPrediction;
710 }
711
712 for(int ii = 1; ii < s_; ++ii)
713 {
714 temp.setZero(); // see https://github.com/bertiniteam/b2/issues/198
715 for(int jj = 0; jj < ii; ++jj)
716 temp += aref(ii,jj)*Kref.col(jj);
717
718 // Vec<ComplexType> wfp =
719 if(EvalRHS<ComplexType>(S, current_space + delta_t*temp, current_time + cref(ii)*delta_t, Kref, ii) != SuccessCode::Success)
720 return SuccessCode::MatrixSolveFailure;
721 }
722
723
724 temp.setZero();
725 for(int ii = 0; ii < s_; ++ii)
726 temp += bref(ii)*Kref.col(ii);
727
728 next_space = current_space + delta_t*temp;
729
730 return SuccessCode::Success;
731 };
732
733
734 template<typename ComplexType>
735 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)
736 {
737 // Calculate condition number and update if needed
738 Eigen::PartialPivLU<Mat<ComplexType>>& LUref = GetLU<ComplexType>();
739 Mat<ComplexType>& dhdxref = std::get< Mat<ComplexType> >(dh_dx_0_);
740
741 // TODO this random vector should not be made fresh every time. especiallyif the numeric type is mpfr_complex!
742 Vec<ComplexType> randy = RandomOfUnits<ComplexType>(numVariables_);
743 Vec<ComplexType> temp_soln = LUref.solve(randy);
744
745 norm_J = NumErrorT(dhdxref.norm());
746 norm_J_inverse = NumErrorT(temp_soln.norm());
747
748 if (num_steps_since_last_condition_number_computation >= frequency_of_CN_estimation)
749 {
750 condition_number_estimate = NumErrorT(norm_J * norm_J_inverse);
751 num_steps_since_last_condition_number_computation = 1; // reset the counter to 1
752 }
753 else // no need to compute the condition number
754 num_steps_since_last_condition_number_computation++;
755 }
756
757
767
768 template<typename ComplexType>
769 SuccessCode SetErrorEstimate(NumErrorT & error_estimate, ComplexType const& delta_t)
770 {
771 using RealType = typename Eigen::NumTraits<ComplexType>::Real;
772
773 Mat<ComplexType>& Kref = std::get< Mat<ComplexType> >(K_);
774 Vec<RealType>& b_minus_bstar_ref = std::get< Vec<RealType> >(b_minus_bstar_);
775
776 auto numFuncs = Kref.rows();
777 Vec<ComplexType> err(numFuncs);
778
779 err.setZero();
780 for(int ii = 0; ii < s_; ++ii)
781 {
782 err += (b_minus_bstar_ref(ii))*Kref.col(ii);
783 }
784
785 err *= delta_t;
786
787 error_estimate = NumErrorT(err.norm());
788
789 return SuccessCode::Success;
790 };
791
792
793
794
795
796
806
807 template<typename ComplexType>
808 SuccessCode SetSizeProportion(NumErrorT & size_proportion, ComplexType const& delta_t)
809 {
810 if(predict::HasErrorEstimate(predictor_))
811 {
812 NumErrorT err_est;
813 SetErrorEstimate(err_est, delta_t);
814
815 using std::pow;
816 size_proportion = err_est/NumErrorT(pow(abs(delta_t), p_+1));
817
818 return SuccessCode::Success;
819 }
820 else
821 {
822 Mat<ComplexType>& Kref = std::get< Mat<ComplexType> >(K_);
823 using std::pow;
824 size_proportion = NumErrorT(Kref.array().abs().maxCoeff()/(pow(abs(delta_t), p_)));
825 return SuccessCode::Success;
826 }
827 };
828
829
830
842
843 template<typename ComplexType>
844 SuccessCode EvalRHS(System const& S,
845 const Vec<ComplexType>& space, const ComplexType& time, Mat<ComplexType> & K, unsigned stage)
846 {
847
848 if (std::is_same<ComplexType, mpfr_complex>::value)
849 PrecisionSanityCheck();
850
851 if(stage == 0)
852 {
853 Eigen::PartialPivLU<Mat<ComplexType>>& LUref = GetLU<ComplexType>();
854 Mat<ComplexType>& dhdxref = std::get< Mat<ComplexType> >(dh_dx_0_);
855
856 if (!std::is_same<ComplexType,dbl>::value)
857 {
858 assert(DefaultPrecision()==current_precision_);
859
860 assert(Precision(space)==current_precision_);
861 assert(Precision(time)==current_precision_);
862 assert(Precision(dhdxref)==current_precision_);
863 assert(Precision(K)==current_precision_);
864 }
865 S.SetAndReset<ComplexType>(space, time);
866 S.JacobianInPlace(dhdxref);
867 LUref = dhdxref.lu();
868 if (!std::is_same<ComplexType,dbl>::value)
869 {
870 assert(Precision(dhdxref)==current_precision_);
871 assert(Precision(LUref.matrixLU())==current_precision_);
872 }
873
874 if (LUPartialPivotDecompositionSuccessful(LUref.matrixLU())!=MatrixSuccessCode::Success)
875 return SuccessCode::MatrixSolveFailureFirstPartOfPrediction;
876
877 Vec<ComplexType>& dhdtref = std::get< Vec<ComplexType> >(dh_dt_temp_);
878 S.TimeDerivativeInPlace(dhdtref);
879 K.col(stage) = LUref.solve(-dhdtref);
880
881 return SuccessCode::Success;
882
883 }
884 else
885 {
886 S.SetAndReset<ComplexType>(space, time);
887
888 Mat<ComplexType>& dhdxtempref = std::get< Mat<ComplexType> >(dh_dx_temp_);
889 S.JacobianInPlace(dhdxtempref);
890 auto LU = dhdxtempref.lu();
891
892 if (LUPartialPivotDecompositionSuccessful(LU.matrixLU())!=MatrixSuccessCode::Success)
893 return SuccessCode::MatrixSolveFailure;
894
895 Vec<ComplexType>& dhdtref = std::get< Vec<ComplexType> >(dh_dt_temp_);
896 S.TimeDerivativeInPlace(dhdtref);
897 K.col(stage) = LU.solve(-dhdtref);
898
899 return SuccessCode::Success;
900 }
901 }
902
903
904
905
906
907
908
909
910 private:
912 //
913 // Private Methods
914 //
916
917
918
929
930 template<typename RealType>
931 void FillButcherTable(int stages, const Mat<mpq_rational>& a,
932 const Mat<mpq_rational> & b,
933 const Mat<mpq_rational> & b_minus_bstar,
934 const Mat<mpq_rational> & c)
935 {
936 Mat<RealType>& aref = std::get< Mat<RealType> >(a_);
937 aref.resize(stages, stages);
938 for(int ii = 0; ii < stages; ++ii)
939 {
940 for(int jj = 0; jj < s_; ++jj)
941 {
942 aref(ii,jj) = static_cast<RealType>(a(ii,jj));
943 }
944 }
945
946 Vec<RealType>& bref = std::get< Vec<RealType> >(b_);
947 bref.resize(stages);
948 for(int ii = 0; ii < stages; ++ii)
949 {
950 bref(ii) = static_cast<RealType>(b(ii));
951 }
952
953 Vec<RealType>& b_minus_bstar_ref = std::get< Vec<RealType> >(b_minus_bstar_);
954 b_minus_bstar_ref.resize(stages);
955 for(int ii = 0; ii < stages; ++ii)
956 {
957 b_minus_bstar_ref(ii) = static_cast<RealType>(b_minus_bstar(ii));
958 }
959
960 Vec<RealType>& cref = std::get< Vec<RealType> >(c_);
961 cref.resize(stages);
962 for(int ii = 0; ii < stages; ++ii)
963 {
964 cref(ii) = static_cast<RealType>(c(ii));
965
966 }
967 uses_embedded_ = true;
968 }
969
970
971
972
973
983
984 template<typename RealType>
985 void FillButcherTable(int stages, const Mat<mpq_rational>& a,
986 const Mat<mpq_rational> & b,
987 const Mat<mpq_rational> & c)
988 {
989 Mat<RealType>& aref = std::get< Mat<RealType> >(a_);
990 aref.resize(stages, stages);
991 for(int ii = 0; ii < stages; ++ii)
992 {
993 for(int jj = 0; jj < s_; ++jj)
994 {
995 aref(ii,jj) = static_cast<RealType>(a(ii,jj));
996 }
997 }
998
999 Vec<RealType>& bref = std::get< Vec<RealType> >(b_);
1000 bref.resize(stages);
1001 for(int ii = 0; ii < stages; ++ii)
1002 {
1003 bref(ii) = static_cast<RealType>(b(ii));
1004 }
1005
1006 Vec<RealType>& cref = std::get< Vec<RealType> >(c_);
1007 cref.resize(stages);
1008 for(int ii = 0; ii < stages; ++ii)
1009 {
1010 cref(ii) = static_cast<RealType>(c(ii));
1011
1012 }
1013 uses_embedded_ = false;
1014 }
1015
1016
1017
1018
1019
1020
1022 //
1023 // Private Data Members
1024 //
1026
1027 unsigned numTotalFunctions_; // Number of total functions for the current system
1028 unsigned numVariables_; // Number of variables for the current system
1029 mutable std::tuple< Mat<dbl>, Mat<mpfr_complex> > K_; // All the stage variables. Each column represents a different stage.
1030 Predictor predictor_; // Method for prediction
1031 unsigned p_; //Order of the prediction method
1032 mutable std::tuple< Mat<dbl>, Mat<mpfr_complex> > dh_dx_0_; // Jacobian for the initial stage. Use for AMP testing
1033 mutable std::tuple< Mat<dbl>, Mat<mpfr_complex> > dh_dx_temp_; // Temporary jacobian for all other stages
1034 mutable std::tuple< Vec<dbl>, Vec<mpfr_complex> > dh_dt_temp_; // Temporary time derivative used for all stages
1035 // std::tuple< Eigen::PartialPivLU<Mat<dbl>>, Eigen::PartialPivLU<Mat<mpfr_complex>> > LU_0_; // LU from the intial stage used for AMP testing
1036
1037 mutable Eigen::PartialPivLU<Mat<dbl>> LU_d_;
1038 mutable std::map<unsigned,Eigen::PartialPivLU<Mat<mpfr_complex>>> LU_mp_;
1039
1040
1041 // Butcher Table (notation from https://en.wikipedia.org/wiki/List_of_Runge%E2%80%93Kutta_methods )
1042 mutable unsigned s_; // Number of stages
1043 mutable std::tuple< Mat<double>, Mat<mpfr_float> > a_;
1044 mutable std::tuple< Vec<double>, Vec<mpfr_float> > b_;
1045 mutable std::tuple< Vec<double>, Vec<mpfr_float> > b_minus_bstar_;
1046 mutable std::tuple< Vec<double>, Vec<mpfr_float> > c_;
1047
1048 mutable bool uses_embedded_;
1049 mutable unsigned current_precision_;
1050
1051
1052
1053
1054
1055 // static const variables that store the Butcher table in mpq_rational form
1056 // Euler
1057 static const mpq_rational aEulerPtr_[];
1058 static const Eigen::Matrix<mpq_rational,1,1> aEuler_;
1059 static const mpq_rational bEulerPtr_[];
1060 static const Eigen::Matrix<mpq_rational,1,1> bEuler_;
1061 static const mpq_rational cEulerPtr_[];
1062 static const Eigen::Matrix<mpq_rational,1,1> cEuler_;
1063
1064 // Heun-Euler
1065 static const mpq_rational aHeunEulerPtr_[];
1066 static const Eigen::Matrix<mpq_rational,2,2> aHeunEuler_;
1067 static const mpq_rational bHeunEulerPtr_[];
1068 static const Eigen::Matrix<mpq_rational,2,1> bHeunEuler_;
1069 static const mpq_rational b_minus_bstarHeunEulerPtr_[];
1070 static const Eigen::Matrix<mpq_rational,2,1> b_minus_bstarHeunEuler_;
1071 static const mpq_rational cHeunEulerPtr_[];
1072 static const Eigen::Matrix<mpq_rational,2,1> cHeunEuler_;
1073
1074 // RK4
1075 static const mpq_rational aRK4Ptr_[];
1076 static const Eigen::Matrix<mpq_rational,4,4> aRK4_;
1077 static const mpq_rational bRK4Ptr_[];
1078 static const Eigen::Matrix<mpq_rational,4,1> bRK4_;
1079 static const mpq_rational cRK4Ptr_[];
1080 static const Eigen::Matrix<mpq_rational,4,1> cRK4_;
1081
1082 // RKF45
1083 static const mpq_rational aRKF45Ptr_[];
1084 static const Eigen::Matrix<mpq_rational,6,6> aRKF45_;
1085 static const mpq_rational bRKF45Ptr_[];
1086 static const Eigen::Matrix<mpq_rational,6,1> bRKF45_;
1087 static const mpq_rational b_minus_bstarRKF45Ptr_[];
1088 static const Eigen::Matrix<mpq_rational,6,1> b_minus_bstarRKF45_;
1089 static const mpq_rational cRKF45Ptr_[];
1090 static const Eigen::Matrix<mpq_rational,6,1> cRKF45_;
1091
1092 // RK Cash-Karp45
1093 static const mpq_rational aRKCK45Ptr_[];
1094 static const Eigen::Matrix<mpq_rational,6,6> aRKCK45_;
1095 static const mpq_rational bRKCK45Ptr_[];
1096 static const Eigen::Matrix<mpq_rational,6,1> bRKCK45_;
1097 static const mpq_rational b_minus_bstarRKCK45Ptr_[];
1098 static const Eigen::Matrix<mpq_rational,6,1> b_minus_bstarRKCK45_;
1099 static const mpq_rational cRKCK45Ptr_[];
1100 static const Eigen::Matrix<mpq_rational,6,1> cRKCK45_;
1101
1102 // RK Dormand-Prince 56
1103 static const mpq_rational aRKDP56Ptr_[];
1104 static const Eigen::Matrix<mpq_rational,8,8> aRKDP56_;
1105 static const mpq_rational bRKDP56Ptr_[];
1106 static const Eigen::Matrix<mpq_rational,8,1> bRKDP56_;
1107 static const mpq_rational b_minus_bstarRKDP56Ptr_[];
1108 static const Eigen::Matrix<mpq_rational,8,1> b_minus_bstarRKDP56_;
1109 static const mpq_rational cRKDP56Ptr_[];
1110 static const Eigen::Matrix<mpq_rational,8,1> cRKDP56_;
1111
1112 // RK Verner 67
1113 static const mpq_rational aRKV67Ptr_[];
1114 static const Eigen::Matrix<mpq_rational,10,10> aRKV67_;
1115 static const mpq_rational bRKV67Ptr_[];
1116 static const Eigen::Matrix<mpq_rational,10,1> bRKV67_;
1117 static const mpq_rational b_minus_bstarRKV67Ptr_[];
1118 static const Eigen::Matrix<mpq_rational,10,1> b_minus_bstarRKV67_;
1119 static const mpq_rational cRKV67Ptr_[];
1120 static const Eigen::Matrix<mpq_rational,10,1> cRKV67_;
1121
1122
1123
1124 }; //re: ExplicitRKPredictor class
1125
1126
1127
1128
1129
1130
1131
1132 } // re: predict
1133 }// re: tracking
1134}// re: bertini
1135
1136#endif
amp_criteria.hpp
Provides the Adaptive Multiple Precision criteria functions.
bertini::tracking::amp::CriterionC
bool CriterionC(double const &norm_J_inverse, const Eigen::MatrixBase< Derived > &z, double tracking_tolerance, AdaptiveMultiplePrecisionConfig const &AMP_config)
Check AMP Criterion C.
Definition amp_criteria.hpp:206
bertini::tracking::amp::CriterionA
bool CriterionA(double const &norm_J, double const &norm_J_inverse, AdaptiveMultiplePrecisionConfig const &AMP_config)
Check AMP Criterion A.
Definition amp_criteria.hpp:71
bertini::tracking::predict::HasErrorEstimate
bool HasErrorEstimate(Predictor predictor_choice)
Ask whether a predictor method provides an error estimate.
Definition base_predictor.hpp:94
bertini::tracking::predict::DefaultPredictor
Predictor DefaultPredictor()
Get the Bertini2 default predictor.
Definition base_predictor.hpp:58
bertini::tracking::predict::Order
unsigned Order(Predictor predictor_choice)
The order of the predictor.
Definition base_predictor.hpp:68
bertini::System
The fundamental polynomial system class for Bertini2.
Definition system/system.hpp:99
bertini::System::TimeDerivativeInPlace
void TimeDerivativeInPlace(Vec< T > &ds_dt, const Eigen::MatrixBase< Derived > &variable_values, const T &path_variable_value) const
Compute the time-derivative of a system.
Definition system/system.hpp:680
bertini::System::JacobianInPlace
void JacobianInPlace(Mat< T > &J) const
Evaluate the Jacobian matrix of the system, using the previous space and time values,...
Definition system/system.hpp:470
bertini::System::NumVariables
size_t NumVariables() const
bertini::System::NumTotalFunctions
size_t NumTotalFunctions() const
bertini::tracking::predict::ExplicitRKPredictor::ChangePrecision
void ChangePrecision(unsigned new_precision)
Definition explicit_predictors.hpp:394
bertini::tracking::predict::ExplicitRKPredictor::precision
unsigned precision() const
get the current precision of the predictor
Definition explicit_predictors.hpp:383
bertini::tracking::predict::ExplicitRKPredictor::EvalRHS
SuccessCode EvalRHS(System const &S, const Vec< ComplexType > &space, const ComplexType &time, Mat< ComplexType > &K, unsigned stage)
Evaluates the RHS of the Davidenko differential equation at a particular time and space.
Definition explicit_predictors.hpp:844
bertini::tracking::predict::ExplicitRKPredictor::PredictorMethod
Predictor PredictorMethod()
Get the currently used prediction method.
Definition explicit_predictors.hpp:611
bertini::tracking::predict::ExplicitRKPredictor::HasErrorEstimate
bool HasErrorEstimate()
Get whether the current prediction method provides an error estimate.
Definition explicit_predictors.hpp:637
bertini::tracking::predict::ExplicitRKPredictor::ExplicitRKPredictor
ExplicitRKPredictor(Predictor method, const System &S)
Constructor for a particular predictor method.
Definition explicit_predictors.hpp:224
bertini::tracking::predict::ExplicitRKPredictor::Predict
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)
Perform a generic predictor step.
Definition explicit_predictors.hpp:460
bertini::tracking::predict::ExplicitRKPredictor::SetSizeProportion
SuccessCode SetSizeProportion(NumErrorT &size_proportion, ComplexType const &delta_t)
Compute the size proportion variable for AMP computation.
Definition explicit_predictors.hpp:808
bertini::tracking::predict::ExplicitRKPredictor::Order
unsigned Order()
Get the order of the currently used prediction method.
Definition explicit_predictors.hpp:627
bertini::tracking::predict::ExplicitRKPredictor::Predict
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)
Perform a generic predictor step and return size_proportion and condition number information.
Definition explicit_predictors.hpp:503
bertini::tracking::predict::ExplicitRKPredictor::ExplicitRKPredictor
ExplicitRKPredictor(const System &S)
Construct a predictor to work on a system.
Definition explicit_predictors.hpp:212
bertini::tracking::predict::ExplicitRKPredictor::PredictorMethod
void PredictorMethod(Predictor method)
Definition explicit_predictors.hpp:238
bertini::tracking::predict::ExplicitRKPredictor::Predict
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)
Perform a generic predictor step and return error estimate, size_proportion and condition number info...
Definition explicit_predictors.hpp:571
bertini::tracking::predict::ExplicitRKPredictor::FullStep
SuccessCode FullStep(Vec< ComplexType > &next_space, System const &S, Vec< ComplexType > const &current_space, ComplexType const &current_time, ComplexType const &delta_t)
Performs a full prediction step from current_time to current_time + delta_t.
Definition explicit_predictors.hpp:685
bertini::tracking::predict::ExplicitRKPredictor::ChangeSystem
void ChangeSystem(const System &S)
Change the system (number of total functions) that the predictor uses.
Definition explicit_predictors.hpp:357
bertini::tracking::predict::ExplicitRKPredictor::SetErrorEstimate
SuccessCode SetErrorEstimate(NumErrorT &error_estimate, ComplexType const &delta_t)
Computes the error estimate of this prediction step.
Definition explicit_predictors.hpp:769
bertini::pow
dbl pow(const dbl &z, int power)
Definition double_extensions.hpp:88
eigen_extensions.hpp
Bertini extensions of the Eigen linear algebra library.
bertini::RandomOfUnits
Mat< NumberType > RandomOfUnits(uint rows, uint cols)
Make a random matrix of units (numbers with norm 1).
Definition eigen_extensions.hpp:502
bertini::LUPartialPivotDecompositionSuccessful
MatrixSuccessCode LUPartialPivotDecompositionSuccessful(Eigen::MatrixBase< Derived > const &LU)
Check the diagonal elements of an LU decomposition for small values and large ratios.
Definition eigen_extensions.hpp:417
bertini::Precision
unsigned Precision(Eigen::MatrixBase< Derived > const &v)
Get the precision of an Eigen object. If the object is empty, it's the precision of a default-constru...
Definition eigen_extensions.hpp:292
mpfr_extensions.hpp
Extensions to the Boost.Multiprecision library.
bertini::tracking::AdaptiveMultiplePrecisionConfig
Definition trackers/config.hpp:154
system.hpp
Provides the bertini::System class.