File: NLfit.cpp

package info (click to toggle)
fityk 1.3.1-6
  • links: PTS, VCS
  • area: main
  • in suites: bullseye
  • 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 (119 lines) | stat: -rw-r--r-- 3,393 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
// This file is part of fityk program. Copyright 2001-2013 Marcin Wojdyr.
// Licence: GNU General Public License ver. 2+

#define BUILDING_LIBFITYK
#include "NLfit.h"
#include "logic.h"
#include "settings.h"
#include "var.h"

#if HAVE_LIBNLOPT

using namespace std;

namespace fityk {

NLfit::NLfit(Full* F, const char* name, nlopt_algorithm algorithm)
    : Fit(F, name), algorithm_(algorithm), opt_(NULL)
{
}

NLfit::~NLfit()
{
    if (opt_ != NULL)
        nlopt_destroy(opt_);
}

static
double calculate_for_nlopt(unsigned n, const double* x,
                           double* grad, void* f_data)
{
    return static_cast<NLfit*>(f_data)->calculate(n, x, grad);
}


double NLfit::calculate(int n, const double* par, double* grad)
{
    assert(n == na_);
    vector<realt> A(par, par+n);
    if (F_->get_verbosity() >= 1)
        output_tried_parameters(A);
    bool stop = common_termination_criteria();
    if (stop)
        nlopt_force_stop(opt_);

    double wssr;
    if (!grad || stop)
        wssr = compute_wssr(A, fitted_datas_);
    else
        wssr = compute_wssr_gradient(A, fitted_datas_, grad);
    if (F_->get_verbosity() >= 1)
        F_->ui()->mesg(iteration_info(wssr));
    return wssr;
}

static
const char* nlresult_to_string(nlopt_result r)
{
    switch (r) {
        case NLOPT_FAILURE: return "failure";
        case NLOPT_INVALID_ARGS: return "invalid arguments";
        case NLOPT_OUT_OF_MEMORY: return "out of memory";
        case NLOPT_ROUNDOFF_LIMITED: return "roundoff errors limit progress";
        case NLOPT_FORCED_STOP: return "interrupted";
        case NLOPT_SUCCESS: return "success";
        case NLOPT_STOPVAL_REACHED: return "stop-value reached";
        case NLOPT_FTOL_REACHED: return "ftol-value reached";
        case NLOPT_XTOL_REACHED: return "xtol-value reached";
        case NLOPT_MAXEVAL_REACHED: return "max. evaluation number reached";
        case NLOPT_MAXTIME_REACHED: return "max. time reached";
    }
    return NULL;
}

double NLfit::run_method(vector<realt>* best_a)
{
    if (opt_ != NULL && na_ != (int) nlopt_get_dimension(opt_)) {
        nlopt_destroy(opt_);
        opt_ = NULL;
    }

    if (opt_ == NULL) {
        opt_ = nlopt_create(algorithm_, na_);
        nlopt_set_min_objective(opt_, calculate_for_nlopt, this);
    }

    // this is also handled in Fit::common_termination_criteria()
    nlopt_set_maxtime(opt_, F_->get_settings()->max_fitting_time);
    nlopt_set_maxeval(opt_, max_eval() - 1); // save 1 eval for final calc.
    nlopt_set_ftol_rel(opt_, F_->get_settings()->ftol_rel);
    nlopt_set_xtol_rel(opt_, F_->get_settings()->xtol_rel);

    if (F_->get_settings()->box_constraints) {
        double *lb = new double[na_];
        double *ub = new double[na_];
        for (int i = 0; i < na_; ++i) {
            const RealRange& d = F_->mgr.get_variable(i)->domain;
            lb[i] = d.lo;
            ub[i] = d.hi;
        }
        nlopt_set_lower_bounds(opt_, lb);
        nlopt_set_upper_bounds(opt_, ub);
        delete [] lb;
        delete [] ub;
    }

    double opt_f;
    double *a = new double[na_];
    for (int i = 0; i < na_; ++i)
        a[i] = a_orig_[i];
    nlopt_result r = nlopt_optimize(opt_, a, &opt_f);
    F_->msg("NLopt says: " + S(nlresult_to_string(r)));
    best_a->assign(a, a+na_);
    delete [] a;
    return opt_f;
}

} // namespace fityk

#endif //HAVE_LIBNLOPT