Line data Source code
1 : #include <catch2/catch_get_random_seed.hpp>
2 : #include <catch2/catch_test_macros.hpp>
3 : #include <catch2/matchers/catch_matchers_floating_point.hpp>
4 : #include <cmath>
5 :
6 : #include "Opt/GeneticIntervalSolver.hpp"
7 : #include "Opt/GoldenSectionSolver.hpp"
8 : #include "Opt/QuadraticSolver.hpp"
9 : #include "Opt/UnivariateSolver.hpp"
10 :
11 : using namespace std;
12 :
13 : using Catch::Matchers::WithinAbs;
14 :
15 6 : void testInterval(
16 : const Opt::IntervalSolver::Interval &p,
17 : const Opt::IntervalSolver::Interval &psol
18 : ) {
19 : // Check consistency of solution
20 12 : REQUIRE(p.first <= p.second);
21 :
22 : // Check range is completely included in range allowed by error around true solution
23 12 : REQUIRE(psol.first <= p.first);
24 12 : REQUIRE(p.second <= psol.second);
25 6 : }
26 :
27 5 : void testIntervalSolverQuadratic(
28 : Opt::IntervalSolver &solver,
29 : const double &a,
30 : const double &b,
31 : const double &c,
32 : const double &lmargin,
33 : const double &rmargin,
34 : const double &e
35 : ) {
36 5 : double sol = -b / (a * 2);
37 5 : Opt::IntervalSolver::Interval psol(sol - e, sol + e);
38 5 : double l = sol - lmargin;
39 5 : double r = sol + rmargin;
40 185 : Opt::UnivariateSolver::Problem prob = [a, b, c](double x) {
41 180 : return (((a)*x + b) * x + c);
42 5 : };
43 5 : solver.setInterval(l, r);
44 5 : solver.setStopCriteria(e);
45 5 : Opt::IntervalSolver::Interval p = solver.solveInterval(prob);
46 :
47 5 : testInterval(p, psol);
48 5 : }
49 :
50 0 : void testSolverQuadratic(
51 : Opt::UnivariateSolver &solver,
52 : const double &a,
53 : const double &b,
54 : const double &c,
55 : const double &e
56 : ) {
57 0 : double sol = -b / (a * 2);
58 0 : Opt::IntervalSolver::Interval psol(sol - e, sol + e);
59 0 : Opt::UnivariateSolver::Problem prob = [a, b, c](double x) {
60 0 : return (((a)*x + b) * x + c);
61 0 : };
62 0 : solver.setStopCriteria(e);
63 0 : double x = solver.solve(prob);
64 :
65 0 : REQUIRE_THAT(x, WithinAbs(sol, e));
66 0 : }
67 :
68 6 : TEST_CASE("Golden Section solver", "[golden-section]") {
69 6 : Catch::StringMaker<float>::precision = 40;
70 :
71 6 : Opt::GoldenSectionSolver solver;
72 :
73 7 : SECTION("Quadratic 1,0,0, error 1e-3, margins 1,1") {
74 1 : double a = 1, b = 0, c = 0, e = 1e-3, lmargin = 1, rmargin = 1;
75 1 : testIntervalSolverQuadratic(solver, a, b, c, lmargin, rmargin, e);
76 : }
77 :
78 7 : SECTION("Quadratic 1,1,1, error 1e-3, margins 10,100") {
79 1 : double a = 1, b = 1, c = 1, e = 1e-3, lmargin = 10, rmargin = 100;
80 1 : testIntervalSolverQuadratic(solver, a, b, c, lmargin, rmargin, e);
81 : }
82 :
83 7 : SECTION("Quadratic 1,0,0, error 1e-7, margins 1,1") {
84 1 : double a = 1, b = 0, c = 0, e = 1e-7, lmargin = 1, rmargin = 1;
85 1 : testIntervalSolverQuadratic(solver, a, b, c, lmargin, rmargin, e);
86 : }
87 :
88 7 : SECTION("Quadratic 1,1,1, error 1e-7, margins 10,100") {
89 1 : double a = 1, b = 1, c = 1, e = 1e-7, lmargin = 10, rmargin = 100;
90 1 : testIntervalSolverQuadratic(solver, a, b, c, lmargin, rmargin, e);
91 : }
92 :
93 7 : SECTION("Quadratic 123,456,789, error 1e-7, margins 10000,10000") {
94 1 : double a = 123, b = 456, c = 789, e = 1e-7, lmargin = 10000, rmargin = 10000;
95 1 : testIntervalSolverQuadratic(solver, a, b, c, lmargin, rmargin, e);
96 : }
97 :
98 7 : SECTION("Absolute value, sol 2, error 1e-7, interval -3,14") {
99 1 : double sol = 2, e = 1e-15, l = -3, r = 14;
100 80 : Opt::UnivariateSolver::Problem prob = [sol](double x) {
101 79 : return fabs(x - sol);
102 2 : };
103 1 : solver.setInterval(l, r);
104 1 : solver.setStopCriteria(e);
105 1 : Opt::IntervalSolver::Interval p = solver.solveInterval(prob);
106 :
107 1 : testInterval(p, make_pair(sol - e, sol + e));
108 : }
109 6 : }
110 :
111 5 : void testIntervalSolverGenetic(
112 : Opt::GeneticIntervalSolver &solver,
113 : const double &a,
114 : const double &b,
115 : const double &c,
116 : const double &lmargin,
117 : const double &rmargin,
118 : const double &e
119 : ) {
120 5 : double sol = -b / (a * 2);
121 5 : Opt::IntervalSolver::Interval psol(sol - e, sol + e);
122 5 : double l = sol - lmargin;
123 5 : double r = sol + rmargin;
124 36270 : Opt::UnivariateSolver::Problem prob = [a, b, c](double x) {
125 36265 : return (((a)*x + b) * x + c);
126 5 : };
127 5 : solver.setInterval(l, r);
128 5 : solver.setStopCriteria(e);
129 5 : Opt::IntervalSolver::Var x = solver.solve(prob);
130 :
131 20 : REQUIRE_THAT(x, WithinAbs(sol, e));
132 5 : }
133 :
134 7 : TEST_CASE("Genetic solver", "[genetic]") {
135 7 : Catch::StringMaker<float>::precision = 40;
136 :
137 7 : Opt::GeneticIntervalSolver solver(
138 : 100,
139 : 1000,
140 : 0.5,
141 : 1000,
142 : 8,
143 7 : std::make_shared<std::mt19937>(0)
144 7 : );
145 :
146 8 : SECTION("Quadratic 1,0,0, error 1e-3, margins 1,1") {
147 1 : double a = 1, b = 0, c = 0, e = 1e-3, lmargin = 1, rmargin = 1;
148 1 : testIntervalSolverGenetic(solver, a, b, c, lmargin, rmargin, e);
149 : }
150 :
151 8 : SECTION("Quadratic 1,1,1, error 1e-3, margins 10,100") {
152 1 : double a = 1, b = 1, c = 1, e = 1e-3, lmargin = 10, rmargin = 100;
153 1 : testIntervalSolverGenetic(solver, a, b, c, lmargin, rmargin, e);
154 : }
155 :
156 8 : SECTION("Quadratic 1,0,0, error 1e-7, margins 1,1") {
157 1 : double a = 1, b = 0, c = 0, e = 1e-7, lmargin = 1, rmargin = 1;
158 1 : testIntervalSolverGenetic(solver, a, b, c, lmargin, rmargin, e);
159 : }
160 :
161 8 : SECTION("Quadratic 1,1,1, error 1e-7, margins 10,100") {
162 1 : double a = 1, b = 1, c = 1, e = 1e-7, lmargin = 10, rmargin = 100;
163 1 : testIntervalSolverGenetic(solver, a, b, c, lmargin, rmargin, e);
164 : }
165 :
166 8 : SECTION("Quadratic 123,456,789, error 1e-7, margins 10000,10000") {
167 1 : double a = 123, b = 456, c = 789, e = 1e-7, lmargin = 10000, rmargin = 10000;
168 1 : testIntervalSolverGenetic(solver, a, b, c, lmargin, rmargin, e);
169 : }
170 :
171 8 : SECTION("Absolute value, sol 2, error 1e-7, interval -3,14") {
172 1 : double sol = 2, e = 1e-15, l = -3, r = 14;
173 16073 : Opt::UnivariateSolver::Problem prob = [sol](double x) {
174 16072 : return fabs(x - sol);
175 2 : };
176 1 : solver.setInterval(l, r);
177 1 : solver.setStopCriteria(e);
178 1 : Opt::IntervalSolver::Var x = solver.solve(prob);
179 :
180 4 : REQUIRE_THAT(x, WithinAbs(sol, e));
181 : }
182 :
183 8 : SECTION("Quadratic formula with local minima") {
184 8091 : Opt::UnivariateSolver::Problem prob = [](double x) {
185 8091 : return pow(x, 4) - pow(x, 3) - 8 * pow(x, 2) + 30;
186 2 : };
187 1 : double l = -5, r = 3;
188 1 : double e = 1e-7;
189 :
190 1 : double sol = +2.40985257;
191 1 : double localMinimum = -1.65985257;
192 :
193 : // To prove the golden section method cannot solve this case because
194 : // the function is not convex
195 2 : Opt::GoldenSectionSolver goldenSolver;
196 1 : goldenSolver.setInterval(l, r);
197 1 : goldenSolver.setStopCriteria(e);
198 1 : Opt::IntervalSolver::Var xGS = goldenSolver.solve(prob);
199 3 : REQUIRE_THAT(xGS, WithinAbs(localMinimum, e));
200 :
201 1 : solver.setInterval(l, r);
202 1 : solver.setStopCriteria(e);
203 1 : Opt::IntervalSolver::Var x = solver.solve(prob);
204 4 : REQUIRE_THAT(x, WithinAbs(sol, e));
205 : }
206 7 : }
|
