File: gradient.cpp

package info (click to toggle)
fityk 1.3.1-6
  • links: PTS, VCS
  • area: main
  • in suites: bookworm, bullseye, sid
  • size: 3,784 kB
  • sloc: cpp: 34,396; ansic: 4,673; python: 971; makefile: 366; sh: 117; java: 31; ruby: 27; perl: 25; xml: 16
file content (131 lines) | stat: -rw-r--r-- 4,548 bytes parent folder | download | duplicates (3)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131

#include <stdio.h>
#include <boost/scoped_ptr.hpp>
#include "fityk/logic.h"
#include "fityk/data.h"
#include "fityk/fit.h"

#include "catch.hpp"

using namespace std;
using namespace fityk;


// Box and Betts exponential quadratic sum, equivalent to least-squares of
// f(x) = exp(-0.1*a0*x) - exp(-0.1*a1*x) - a2 * (exp(-0.1*x) - exp(-x))
// with points (1,0), (2,0), ... (10,0)
// 
// domains of a0, a1, a2 are, respectively, (0.9, 1.2), (9, 11.2), (0.9, 1.2)
// minimum: (1,10,1) -> 0

// modified boxbetts_f() from nlopt-2.3/test/testfuncs.c
static double boxbetts_f(const double *a, double *grad)
{
    double f = 0;
    if (grad)
        grad[0] = grad[1] = grad[2] = 0;
    for (int x = 1; x <= 10; ++x) {
        double e0 = exp(-0.1*x*a[0]);
        double e1 = exp(-0.1*x*a[1]);
        double e2 = exp(-0.1*x) - exp(-x);
        double g = e0 - e1 - e2 * a[2];
        f += g*g;
        if (grad) {
            grad[0] += (2 * g) * (-0.1*x*e0);
            grad[1] += (2 * g) * (0.1*x*e1);
            grad[2] += -(2 * g) * e2;
        }
    }
    return f;
}

static double boxbetts_in_fityk(const double *a, double *grad)
{
    boost::scoped_ptr<Fityk> ftk(new Fityk);
    Full* priv = ftk->priv();
    ftk->set_option_as_number("verbosity", -1);
    for (int i = 1; i <= 10; ++i)
        priv->dk.data(0)->add_one_point(i, 0, 1);
    ftk->execute("define BoxBetts(a0,a1,a2) = "
            "exp(-0.1*a0*x) - exp(-0.1*a1*x) - a2 * (exp(-0.1*x) - exp(-x))");
    ftk->execute("F = BoxBetts(~0.9, ~11.8, ~1.08)");

    priv->get_fit()->get_dof(priv->dk.datas()); // to invoke update_par_usage()
    vector<realt> avec(a, a+3);
    return priv->get_fit()->compute_wssr_gradient(avec, priv->dk.datas(), grad);
}

#if 0
static int test_gradient()
{
    const double a[3] = { 0.9, 11.8, 1.08 };
    double grad[3], grad_again[3];
    double ssr = boxbetts_f(a, grad);
    printf("ssr=%.10g  grad: %.10g %.10g %.10g\n",
            ssr, grad[0], grad[1], grad[2]);

    double ssr_again = boxbetts_in_fityk(a, grad_again);
    printf("ssr=%.10g  grad: %.10g %.10g %.10g\n",
            ssr_again, grad_again[0], grad_again[1], grad_again[2]);
    return 0;
}
int main() { return test_gradient(); }
#endif


TEST_CASE("gradient", "test Fit::compute_wssr_gradient()") {
    const double a[3] = { 0.9, 11.8, 1.08 };
    double grad[3], grad_again[3];
    double ssr = boxbetts_f(a, grad);
    double ssr_again = boxbetts_in_fityk(a, grad_again);
    REQUIRE(ssr == Approx(ssr_again));
    REQUIRE(grad[0] == Approx(grad_again[0]));
    REQUIRE(grad[1] == Approx(grad_again[1]));
    REQUIRE(grad[2] == Approx(grad_again[2]));
}

//----------- + some unrelated random tests

TEST_CASE("set-throws", "test Fityk::set_throws()") {
    boost::scoped_ptr<Fityk> fik(new Fityk);
    REQUIRE(fik->get_throws() == true);
    REQUIRE_THROWS_AS(fik->execute("hi"), SyntaxError);
    REQUIRE_THROWS_AS(fik->execute("fit"), ExecuteError);
    fik->set_throws(false);
    REQUIRE_NOTHROW(fik->execute("hi"));
    fik->execute("hi");
    REQUIRE(fik->last_error() != "");
    REQUIRE(fik->last_error().substr(0,12) == "SyntaxError:");
    fik->clear_last_error();
    REQUIRE(fik->last_error() == "");
}


TEST_CASE("set-get-option", "test Fityk::[gs]et_option_as_*()") {
    boost::scoped_ptr<Fityk> fik(new Fityk);

    REQUIRE(fik->get_option_as_string("numeric_format") == "%g"); // string
    REQUIRE(fik->get_option_as_string("on_error") == "stop"); // enum
    REQUIRE(fik->get_option_as_number("verbosity") == 0.); // integer
    REQUIRE(fik->get_option_as_number("fit_replot") == 0.); // bool
    REQUIRE(fik->get_option_as_number("epsilon") == 1e-12); // double

    // no such option
    REQUIRE_THROWS_AS(fik->get_option_as_string("hi"), ExecuteError);
    // non-numeric option
    REQUIRE_THROWS_AS(fik->get_option_as_number("on_error"), ExecuteError);
    REQUIRE_THROWS_AS(fik->get_option_as_number("cwd"), ExecuteError);

    fik->set_option_as_string("numeric_format", "%.8E");
    fik->set_option_as_string("on_error", "nothing");
    fik->set_option_as_number("verbosity", -1);
    fik->set_option_as_number("fit_replot", 1);
    fik->set_option_as_number("epsilon", 1e-10);

    REQUIRE(fik->get_option_as_string("numeric_format") == "%.8E"); // string
    REQUIRE(fik->get_option_as_string("on_error") == "nothing"); // enum
    REQUIRE(fik->get_option_as_number("verbosity") == -1); // integer
    REQUIRE(fik->get_option_as_number("fit_replot") == 1.); // bool
    REQUIRE(fik->get_option_as_number("epsilon") == 1e-10); // double
}