flexiblesusy is hosted by Hepforge, IPPP Durham
FlexibleSUSY
fixed_point_iterator.hpp
Go to the documentation of this file.
1// ====================================================================
2// This file is part of FlexibleSUSY.
3//
4// FlexibleSUSY is free software: you can redistribute it and/or modify
5// it under the terms of the GNU General Public License as published
6// by the Free Software Foundation, either version 3 of the License,
7// or (at your option) any later version.
8//
9// FlexibleSUSY is distributed in the hope that it will be useful, but
10// WITHOUT ANY WARRANTY; without even the implied warranty of
11// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12// General Public License for more details.
13//
14// You should have received a copy of the GNU General Public License
15// along with FlexibleSUSY. If not, see
16// <http://www.gnu.org/licenses/>.
17// ====================================================================
18
19#ifndef FIXED_POINT_ITERATOR_H
20#define FIXED_POINT_ITERATOR_H
21
22#include <algorithm>
23#include <cmath>
24#include <limits>
25#include <string>
26#include <utility>
27#include <Eigen/Core>
28
29#include "logger.hpp"
30#include "eigen_utils.hpp"
31#include "error.hpp"
32#include "ewsb_solver.hpp"
33
34namespace flexiblesusy {
35
36namespace fixed_point_iterator {
37
38template <std::size_t dimension>
40public:
41 using Vector_t = Eigen::Matrix<double,dimension,1>;
42 enum Status : int { SUCCESS, CONTINUE };
43
44 explicit Convergence_tester_absolute(double precision_ = 1.0e-2)
45 : precision(precision_)
46 {}
47
48 const char* name() const { return "Convergence_tester_absolute"; }
49
59 int operator()(const Vector_t& a, const Vector_t& b) const {
60 const double res = std::sqrt((a - b).array().square().sum());
61 return res < precision ? SUCCESS : CONTINUE;
62 }
63
64private:
65 double precision;
66};
67
68template <std::size_t dimension>
70public:
71 using Vector_t = Eigen::Matrix<double,dimension,1>;
72 enum Status : int { SUCCESS, CONTINUE };
73
74 explicit Convergence_tester_relative(double precision_ = 1.0e-2)
75 : precision(precision_)
76 {}
77
78 const char* name() const { return "Convergence_tester_relative"; }
79
89 int operator()(const Vector_t& a, const Vector_t& b) const {
90 return is_equal_rel(a, b, precision) ? SUCCESS : CONTINUE;
91 }
92
93private:
94 double precision;
95};
96
97template <std::size_t dimension>
99public:
100 using Vector_t = Eigen::Matrix<double,dimension,1>;
101 using Function_t = std::function<Vector_t(const Vector_t&)>;
102 enum Status : int { SUCCESS, CONTINUE };
103
104 Convergence_tester_tadpole(double precision_,
105 const Function_t& tadpole_function_)
106 : precision(precision_)
107 , tadpole_function(tadpole_function_)
108 {}
109
110 const char* name() const { return "Convergence_tester_tadpole"; }
111
124 int operator()(const Vector_t& a, const Vector_t& b) const
125 {
126 if (!is_equal_rel(a, b, precision)) {
127 return CONTINUE;
128 }
129
130 static const double eps = 10*std::pow(10., -std::numeric_limits<double>::digits10);
131
132 if (is_equal_rel(a, b, eps)) {
133 return SUCCESS;
134 }
135
136 return check_tadpoles(a);
137 }
138
139private:
140 double precision;
142
143 int check_tadpoles(const Vector_t& x) const {
144 const double res = x.cwiseAbs().sum();
145 return res < precision ? SUCCESS : CONTINUE;
146 }
147};
148
149} // namespace fixed_point_iterator
150
172template <std::size_t dimension, class Convergence_tester = fixed_point_iterator::Convergence_tester_relative<dimension>>
174public:
175 using Vector_t = Eigen::Matrix<double,dimension,1>;
176 using Function_t = std::function<Vector_t(const Vector_t&)>;
177
179 Fixed_point_iterator(const Function_t&, std::size_t, const Convergence_tester&);
180 virtual ~Fixed_point_iterator() = default;
181
182 void set_function(const Function_t& f) { function = f; }
183 void set_max_iterations(std::size_t n) { max_iterations = n; }
184 int find_fixed_point(const Eigen::VectorXd&);
185
186 // EWSB_solver interface methods
187 virtual std::string name() const override { return std::string("Fixed_point_iterator<") + convergence_tester.name() + ">"; }
188 virtual int solve(const Eigen::VectorXd& start) override { return find_fixed_point(start); }
189 virtual Eigen::VectorXd get_solution() const override { return fixed_point; }
190
191private:
192 std::size_t max_iterations{100};
193 Vector_t xn{Vector_t::Zero()};
194 Vector_t fixed_point{Vector_t::Zero()};
197
199 void print_state(std::size_t) const;
200
201 static bool is_finite(const Vector_t& v) {
202 return std::any_of(v.data(), v.data() + v.size(),
203 [](double x) { return std::isfinite(x); });
204 }
205};
206
214template <std::size_t dimension, class Convergence_tester>
216 const Function_t& function_,
217 std::size_t max_iterations_,
218 const Convergence_tester& convergence_tester_
219)
220 : max_iterations(max_iterations_)
221 , function(function_)
222 , convergence_tester(convergence_tester_)
223{
224}
225
233template <std::size_t dimension, class Convergence_tester>
235 const Eigen::VectorXd& start
236)
237{
238 if (!function) {
239 throw SetupError("Fixed_point_iterator: function not callable");
240 }
241
242 int status;
243 std::size_t iter = 0;
244
245 fixed_point = xn = start;
246
247 print_state(iter);
248
249 do {
250 iter++;
251 status = fixed_point_iterator_iterate();
252
253 print_state(iter);
254
255 if (status == EWSB_solver::FAIL) {
256 break;
257 }
258
259 status = convergence_tester(fixed_point, xn);
260
261 } while (status == Convergence_tester::CONTINUE && iter < max_iterations);
262
263 VERBOSE_MSG("\t\t\tFixed_point_iterator status = " << status);
264
265 return status;
266}
267
273template <std::size_t dimension, class Convergence_tester>
275{
276 int status = EWSB_solver::SUCCESS;
277 xn = fixed_point;
278
279 try {
280 fixed_point = function(xn);
281 } catch (const flexiblesusy::Error&) {
282 status = EWSB_solver::FAIL;
283 }
284
285 // For safety, include a check for nans or infs
286 if (!is_finite(fixed_point)) {
287 status = EWSB_solver::FAIL;
288 }
289
290 return status;
291}
292
298template <std::size_t dimension, class Convergence_tester>
300{
301 VERBOSE_MSG("\t\t\tIteration n = " << iteration
302 << ": x_{n} = [" << xn.transpose()
303 << "], x_{n+1} = [" << fixed_point.transpose() << "]");
304}
305
306} // namespace flexiblesusy
307
308#endif
interface for numeric EWSB solvers
Definition: ewsb_solver.hpp:31
Eigen::Matrix< double, dimension, 1 > Vector_t
Vector_t fixed_point
vector of fixed point estimate
int find_fixed_point(const Eigen::VectorXd &)
Convergence_tester convergence_tester
convergence tester
Vector_t xn
current iteration point
std::size_t max_iterations
maximum number of iterations
Function_t function
function defining fixed point
virtual std::string name() const override
virtual ~Fixed_point_iterator()=default
virtual int solve(const Eigen::VectorXd &start) override
std::function< Vector_t(const Vector_t &)> Function_t
static bool is_finite(const Vector_t &v)
virtual Eigen::VectorXd get_solution() const override
Spectrum generator was not setup correctly.
Definition: error.hpp:46
int operator()(const Vector_t &a, const Vector_t &b) const
int operator()(const Vector_t &a, const Vector_t &b) const
const Function_t tadpole_function
function to calculate tadpole
Convergence_tester_tadpole(double precision_, const Function_t &tadpole_function_)
int operator()(const Vector_t &a, const Vector_t &b) const
#define VERBOSE_MSG(msg)
Definition: logger.hpp:57
std::complex< double > f(double tau) noexcept
bool is_finite(const gsl_vector *x)
Returns true if GSL vector contains only finite elements, false otherwise.
Definition: gsl_utils.cpp:32
bool is_equal_rel(const Eigen::PlainObjectBase< Derived > &a, const Eigen::PlainObjectBase< Derived > &b, double eps)
Definition: eigen_utils.hpp:89
auto sum(Idx ini, Idx fin, Function f) -> decltype(EvalEigenXpr< Idx >(ini, f))
Definition: sum.hpp:91