250 lines
8 KiB
C++
250 lines
8 KiB
C++
#include <catch2/catch.hpp>
|
|
#include <iostream>
|
|
|
|
#include "../utils.hpp"
|
|
#include "../TOMS748.hpp"
|
|
#include "utils_test.hpp"
|
|
|
|
#define print(x) PRINT_VAR(x);
|
|
#define printvec(x) PRINT_VEC(x);
|
|
#define printstr(x) PRINT_STR(x);
|
|
|
|
template<typename T> T f(T E) { return E - 0.5*sin(E) - 0.3; }
|
|
template<typename T> T f2(T x) { return exp(x)-T(2); }
|
|
template<typename T> T f3(T x) { return (x*x)-T(1); }
|
|
|
|
using namespace TOMS748;
|
|
using namespace TOMS748::internal;
|
|
|
|
TEST_CASE( "sureAbs", "[TOMS748]" )
|
|
{
|
|
REQUIRE(TOMS748::internal::sureAbs(-1) == 1);
|
|
REQUIRE(TOMS748::internal::sureAbs(-1.) == 1.);
|
|
}
|
|
|
|
TEST_CASE( "checkTwoValuesIdentical", "[TOMS748]" )
|
|
{
|
|
std::cout.precision(16);
|
|
|
|
SECTION( "checkTwoValuesIdentical" ) {
|
|
REQUIRE(checkTwoValuesIdentical(1,2,3,4) == false);
|
|
REQUIRE(checkTwoValuesIdentical(1,2,3,1) == true);
|
|
REQUIRE(checkTwoValuesIdentical(1,1,3,4) == true);
|
|
REQUIRE(checkTwoValuesIdentical(1,2,2,4) == true);
|
|
REQUIRE(checkTwoValuesIdentical(1,2,3,3) == true);
|
|
REQUIRE(checkTwoValuesIdentical(1,2,1,4) == true);
|
|
REQUIRE(checkTwoValuesIdentical(1,2,3,2) == true);
|
|
}
|
|
}
|
|
|
|
TEST_CASE( "swap_ab_fafb_if_not_ascending", "[TOMS748]" )
|
|
{
|
|
std::cout.precision(16);
|
|
|
|
SECTION( "ascending (do nothing)" ) {
|
|
double a = 1., b = 2.,
|
|
fa = f(a), fb = f(b);
|
|
|
|
swap_ab_fafb_if_not_ascending(a,b,fa,fb);
|
|
|
|
REQUIRE(a == 1.);
|
|
REQUIRE(b == 2.);
|
|
REQUIRE(fa == f(1.));
|
|
REQUIRE(fb == f(2.));
|
|
}
|
|
|
|
SECTION( "descending (swap)" ) {
|
|
double a = 2., b = 1.,
|
|
fa = f(a), fb = f(b);
|
|
|
|
swap_ab_fafb_if_not_ascending(a,b,fa,fb);
|
|
|
|
REQUIRE(a == 1.);
|
|
REQUIRE(b == 2.);
|
|
REQUIRE(fa == f(1.));
|
|
REQUIRE(fb == f(2.));
|
|
}
|
|
}
|
|
|
|
TEST_CASE( "Bracket", "[TOMS748]" )
|
|
{
|
|
std::cout.precision(16);
|
|
|
|
double a = 0, b = 1, c = 0.3, d = 1.2;
|
|
double fa = f(a), fb = f(b), fc = f(c), fd = f(d);
|
|
double tol = 1e-12;
|
|
bool ok;
|
|
|
|
SECTION( "checkTwoValuesIdentical" ) {
|
|
REQUIRE(check_almost_equal(fbracket1(a, b, fa, fb), 0.5792645075960517, tol));
|
|
REQUIRE(check_almost_equal(fbracket2(a, b, d, fa, fb, fd), 0.1619293662546865, tol));
|
|
}
|
|
|
|
SECTION( "bracket normal operation [c b]" ) {
|
|
std::tie(a, b, d, fa, fb, fd) = bracket(a, b, c, fa, fb, fc);
|
|
REQUIRE(check_almost_equal(a, 0.3, tol));
|
|
REQUIRE(check_almost_equal(b, 1., tol));
|
|
REQUIRE(check_almost_equal(d, 0., tol));
|
|
REQUIRE(check_almost_equal(fa, f(0.3), tol));
|
|
REQUIRE(check_almost_equal(fb, f(1.), tol));
|
|
REQUIRE(check_almost_equal(fd, f(0.), tol));
|
|
}
|
|
|
|
SECTION( "bracket normal operation [a c]" ) {
|
|
std::tie(a, b, d, fa, fb, fd) = bracket(a, b, 0.7, fa, fb, f(0.7));
|
|
REQUIRE(check_almost_equal(a, 0., tol));
|
|
REQUIRE(check_almost_equal(b, 0.7, tol));
|
|
REQUIRE(check_almost_equal(d, 1., tol));
|
|
REQUIRE(check_almost_equal(fa, f(0.), tol));
|
|
REQUIRE(check_almost_equal(fb, f(0.7), tol));
|
|
REQUIRE(check_almost_equal(fd, f(1.), tol));
|
|
}
|
|
|
|
SECTION( "bracket and check normal operation [c b]" ) {
|
|
std::tie(a, b, d, fa, fb, fd, ok) = bracketAndCheckConvergence(a, b, c, fa, fb, fc, tol);
|
|
REQUIRE(check_almost_equal(a, 0.3, tol));
|
|
REQUIRE(check_almost_equal(b, 1., tol));
|
|
REQUIRE(check_almost_equal(d, 0., tol));
|
|
REQUIRE(check_almost_equal(fa, f(0.3), tol));
|
|
REQUIRE(check_almost_equal(fb, f(1.), tol));
|
|
REQUIRE(check_almost_equal(fd, f(0.), tol));
|
|
REQUIRE_FALSE(ok);
|
|
}
|
|
|
|
SECTION( "bracket and check normal operation [a c]" ) {
|
|
std::tie(a, b, d, fa, fb, fd, ok) = bracketAndCheckConvergence(a, b, 0.7, fa, fb, f(0.7), tol);
|
|
REQUIRE(check_almost_equal(a, 0., tol));
|
|
REQUIRE(check_almost_equal(b, 0.7, tol));
|
|
REQUIRE(check_almost_equal(d, 1., tol));
|
|
REQUIRE(check_almost_equal(fa, f(0.), tol));
|
|
REQUIRE(check_almost_equal(fb, f(0.7), tol));
|
|
REQUIRE(check_almost_equal(fd, f(1.), tol));
|
|
REQUIRE_FALSE(ok);
|
|
}
|
|
|
|
SECTION( "bracket and check normal operation c is root" ) {
|
|
std::tie(a, b, d, fa, fb, fd, ok) = bracketAndCheckConvergence(a, b, 0.569682256443945, fa, fb, f(0.569682256443945), tol);
|
|
REQUIRE(check_almost_equal(a, 0., tol));
|
|
REQUIRE(check_almost_equal(b, 0.569682256443945, tol));
|
|
REQUIRE(check_almost_equal(d, 1., tol));
|
|
REQUIRE(check_almost_equal(fa, f(0.), tol));
|
|
REQUIRE(check_almost_equal(fb, f(0.569682256443945), tol));
|
|
REQUIRE(check_almost_equal(fd, f(1.), tol));
|
|
REQUIRE(ok);
|
|
}
|
|
}
|
|
|
|
TEST_CASE( "TOMS748_solve1", "[TOMS748]")
|
|
{
|
|
double a = 0, b = 1, tol = 1e-12;
|
|
|
|
SECTION( "TOMS748_solve1 good" ) {
|
|
auto [x, fx, ok] = TOMS748_solve1(f<double>, a, b, tol);
|
|
CHECK(check_almost_equal(fx, 0., tol));
|
|
CHECK(ok == true);
|
|
}
|
|
|
|
SECTION( "TOMS748_solve1 wrong initial bracket" ) {
|
|
auto [x, fx, ok] = TOMS748_solve1(f<double>, 1., 2., tol);// still converges !
|
|
CHECK(check_almost_equal(fx, 0., tol));
|
|
CHECK(ok == true);
|
|
}
|
|
}
|
|
|
|
TEST_CASE( "findEnclosingBracket", "[TOMS748]")
|
|
{
|
|
double tol = 1e-12, factor = 2;
|
|
|
|
SECTION( "findEnclosingBracket start at 0" ) {
|
|
double x = 0;
|
|
auto [a, b, fa, fb, ok] = findEnclosingBracket(f2<double>, x, factor);
|
|
CHECK(a == x);
|
|
CHECK(a == b);
|
|
CHECK_FALSE(fa*fb <= 0.);
|
|
CHECK_FALSE(ok);
|
|
}
|
|
|
|
SECTION( "findEnclosingBracket sol below" ) {
|
|
double x = 2.;
|
|
auto [a, b, fa, fb, ok] = findEnclosingBracket(f2<double>, x, factor);
|
|
CHECK(check_almost_equal(b, x, tol));
|
|
CHECK(a < b);
|
|
CHECK(fa*fb <= 0.);
|
|
CHECK(ok);
|
|
}
|
|
|
|
SECTION( "findEnclosingBracket sol above" ) {
|
|
double x = .1;
|
|
auto [a, b, fa, fb, ok] = findEnclosingBracket(f2<double>, x, factor);
|
|
CHECK(check_almost_equal(a, x, tol));
|
|
CHECK(a < b);
|
|
CHECK(fa*fb <= 0.);
|
|
CHECK(ok);
|
|
}
|
|
|
|
SECTION( "findEnclosingBracket initial bracket good" ) {
|
|
double x = .5;
|
|
auto [a, b, fa, fb, ok] = findEnclosingBracket(f2<double>, x, factor);
|
|
CHECK(check_almost_equal(a, x, tol));
|
|
CHECK(a < b);
|
|
CHECK(fa*fb <= 0.);
|
|
CHECK(ok);
|
|
}
|
|
|
|
SECTION( "findEnclosingBracket sol above far far away" ) {
|
|
double x = 1e-10;
|
|
auto [a, b, fa, fb, ok] = findEnclosingBracket(f2<double>, x, factor);
|
|
CHECK(check_almost_equal(a, x, tol));
|
|
CHECK(a < b);
|
|
CHECK(fa*fb <= 0.);
|
|
CHECK(ok);
|
|
}
|
|
|
|
SECTION( "findEnclosingBracket sol on the other side of 0 (does not converge)" ) {
|
|
double x = -1.;
|
|
auto [a, b, fa, fb, ok] = findEnclosingBracket(f2<double>, x, factor);
|
|
CHECK(a < b);
|
|
CHECK_FALSE(fa*fb <= 0.);
|
|
CHECK_FALSE(ok);
|
|
}
|
|
|
|
SECTION( "findEnclosingBracket sol below (negative)" ) {
|
|
double x = -.1;
|
|
auto [a, b, fa, fb, ok] = findEnclosingBracket(f3<double>, x, factor);
|
|
CHECK(check_almost_equal(b, x, tol));
|
|
CHECK(a < b);
|
|
CHECK(fa*fb <= 0.);
|
|
CHECK(ok);
|
|
}
|
|
|
|
SECTION( "findEnclosingBracket sol above (negative)" ) {
|
|
double x = -2.;
|
|
auto [a, b, fa, fb, ok] = findEnclosingBracket(f3<double>, x, factor);
|
|
CHECK(check_almost_equal(a, x, tol));
|
|
CHECK(a < b);
|
|
CHECK(fa*fb <= 0.);
|
|
CHECK(ok);
|
|
}
|
|
}
|
|
|
|
TEST_CASE( "bracket_and_solve", "[TOMS748]" )
|
|
{
|
|
double tol = 1e-12;
|
|
SECTION( "initial guess close to solution" ) {
|
|
auto [x, fx, ok] = bracket_and_solve(f<double>, 0.5, tol);
|
|
CHECK(check_almost_equal(fx, 0., tol));
|
|
CHECK(ok == true);
|
|
}
|
|
|
|
SECTION( "initial guess far from solution" ) {
|
|
auto [x, fx, ok] = bracket_and_solve(f<double>, 100., tol);
|
|
CHECK(check_almost_equal(fx, 0., tol));
|
|
CHECK(ok == true);
|
|
}
|
|
|
|
SECTION( "initial guess on other side of 0 compared to solution" ) {
|
|
auto [x, fx, ok] = bracket_and_solve(f<double>, -1., tol);
|
|
CHECK_FALSE(check_almost_equal(fx, 0., tol));
|
|
CHECK_FALSE(ok);
|
|
}
|
|
}
|