File: LMfit.cpp

package info (click to toggle)
fityk 0.4.4-1
  • links: PTS
  • area: main
  • in suites: sarge
  • size: 2,472 kB
  • ctags: 2,617
  • sloc: cpp: 19,705; sh: 5,965; xml: 2,325; yacc: 356; makefile: 183; lex: 178
file content (150 lines) | stat: -rw-r--r-- 4,738 bytes parent folder | download
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
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
// This file is part of fityk program. Copyright (C) Marcin Wojdyr
#include "common.h"
RCSID ("$Id: LMfit.cpp,v 1.6 2004/06/17 19:01:12 wojdyr Exp $")

#include "LMfit.h"
#include "ui.h"
#include <math.h>
#include <vector>
#include <algorithm>

using namespace std;

LMfit::LMfit () 
    : v_fit ('m', "Levenberg-Marquardt"),
      lambda_starting_value (0.001),
      lambda_up_factor (10), lambda_down_factor (10),
      stop_rel (1e-4), shake_before (0), shake_type ('u'),
      alpha(0), alpha_(0), beta(0), beta_(0)
{
    fpar ["lambda-starting-value"] = &lambda_starting_value;
    fpar ["lambda-up-factor"] = &lambda_up_factor;
    fpar ["lambda-down-factor"] = &lambda_down_factor;
    fpar ["stop-rel-change"] = &stop_rel;
    fpar ["shake-before"] = &shake_before;
    epar.insert(pair<string, Enum_string>("shake-type", 
                               Enum_string (Distrib_enum, &shake_type)));
}    
    
LMfit::~LMfit () {}

// WSSR is also called chi2
fp LMfit::init ()   
{
    alpha.resize (na*na);
    alpha_.resize (na*na);
    beta.resize (na);
    beta_.resize (na);
    if (na < 1 ) {
        warn ("What should I fit ?");
        return -1;
    }
    lambda = lambda_starting_value;
    if (shake_before > 0.) {
        for (int i = 0; i < na; i++) 
            a[i] = draw_a_from_distribution (i, shake_type, shake_before);
    }
    else
        a = a_orig; 

    info (print_matrix (a, 1, na, "Initial A"));
    //no need to optimise it (and compute chi2 and derivatives together)
    chi2 = compute_wssr (a);
    compute_derivatives(a, alpha, beta);
    return chi2;
}

int LMfit::autoiter () 
{
    wssr_before = (shake_before > 0. ? compute_wssr() : chi2);
    fp prev_chi2 = chi2;
    verbose("\t === Levenberg-Marquardt method ===");
    info ("Initial values:  lambda=" + S(lambda) + "  WSSR=" + S(chi2));
    verbose ("Max. number of iterations: " + max_iterations);
    if (stop_rel > 0) {
        verbose ("Stopping when relative change of WSSR is "
                  "second time below " + S (stop_rel * 100.) + "%");
    }
    bool converged = false;
    int small_change_counter = 0;
    for (int iter = 0; !common_termination_criteria(iter); iter++) {
        int result = do_iteration();
        if (result < 0) {
            warn ("Error when processing iteration " + S(iter+1) + ".");
            return result;
        }
        if (result == 1) { //better fit
            fp d = prev_chi2 - chi2;
            if (iter % output_one_of == 0)
                info ("#" + S(iter_nr) + ":  WSSR=" + S(chi2) 
                        + "  lambda=" + S(lambda) + "  d(WSSR)=" +  S(-d) 
                        + "  (" + S (d / prev_chi2 * 100) + "%)");  
            if (d / prev_chi2 < stop_rel || chi2 == 0) { //another termination
                small_change_counter++;                  // criterium:
                if (small_change_counter >= 2 || chi2 == 0) { //second time
                    info("Fit converged.");              // neglectable change 
                    converged = true;                    // of chi2; or chi2==0
                    break;
                }
            }
            else
                small_change_counter = 0;
            prev_chi2 = chi2;
            iteration_plot(a);
        }
        else { // result == 0, worse fit
            info ("#" + S(iter_nr) + ": (WSSR=" + S(chi2_) 
                    + ")  lambda=" + S(lambda));
        }
    }
    post_fit (a, chi2);
    return 1;
}

int LMfit::do_iteration()
    //pre: init() callled
{
    if (na < 1) {
        warn ("What am I to fit ?");
        return -1;
    }
    iter_nr++;
    alpha_ = alpha;
    for (int j = 0; j < na; j++) 
        alpha_[na * j + j] *= (1.0 + lambda);
    beta_ = beta;
#ifdef debug
    info (print_matrix (beta_, 1, na, "beta"));
    info (print_matrix (alpha_, na, na, "alpha'"));
#endif /*debug*/

    // Matrix solution (Ax=b)  alpha_ * da == beta_
    if (Jordan (alpha_, beta_, na) < 0)
        return -1;

    // da is in beta_  
    if (getUI()->getVerbosity() >= 4) {
        vector<fp> rel (na);
        for (int q = 0; q < na; q++)
            rel[q] = beta_[q] / a[q] * 100;
        verbose (print_matrix (rel, 1, na, "delta(A)/A[%]"));
    }
    for (int i = 0; i < na; i++) 
        beta_[i] = a[i] + beta_[i];   // and now there is new a[] in beta_[] 
    verbose_lazy (print_matrix (beta_, 1, na, "Trying A"));
    //  compute chi2_
    chi2_ = compute_wssr (beta_);

    if (chi2_ < chi2) { // better fitting
        chi2 = chi2_; 
        a = beta_;
        compute_derivatives(a, alpha, beta);
        lambda /= lambda_down_factor;
        return 1;
    }
    else {// worse fitting
        lambda *= lambda_up_factor;
        return 0;
    }
}