tlsfit.cpp - numeric - C++ library with numerical algorithms
 (HTM) git clone git://src.adamsgaard.dk/numeric
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) LICENSE
       ---
       tlsfit.cpp (1246B)
       ---
            1 #include <iostream>
            2 #include <cmath> // NaN is here
            3 #include <armadillo>
            4 #include "header.h"
            5 #include "lsfit.h"
            6 #include "qrfunc.h"
            7 
            8 // A number of i fitting functions, evaluated at x
            9 Floattype LSfit::fitfuncts(char i, Floattype x)
           10 {
           11   // Choose function i
           12   switch (i) {
           13     case 0:
           14       return 1.0f;
           15       break;
           16     case 1:
           17       return x;
           18       break;
           19     case 2:
           20       return x*x;
           21       break;
           22     default:
           23       std::cout << "Wrong function (i = "
           24                 << i << ") specified in fitfuncts, x = " 
           25                 << x << '\n';
           26       return NAN;
           27   }
           28 }
           29 
           30 // Constructor
           31 LSfit::LSfit(arma::Col<Floattype> &x,
           32                  arma::Col<Floattype> &y,
           33              arma::Col<Floattype> &delta_y) : n(x.n_rows), m(3), A(n,m), b(n), c(m)
           34 {
           35 
           36   // Initialize b values
           37   b = y/delta_y;
           38 
           39   // Initialize A values
           40   for (Lengthtype i=0; i<n; ++i) {
           41     for (Lengthtype k=0; k<m; ++k) {
           42       A(i,k) = fitfuncts(k,x(i)) / delta_y(i);
           43     }
           44   }
           45 
           46   // Perform QR decomposition
           47   QR qr(A);
           48 
           49   // Calculate the least-squares solution
           50   c = qr.backsub(b);
           51 }
           52 
           53 // Evaluate function at x
           54 Floattype LSfit::eval(Floattype x)
           55 {
           56   // Find the linear combination of m functions to calculate F(x)
           57   Floattype sum = 0.0f;
           58   for (char k=0; k<m; ++k)
           59     sum += c(k) * fitfuncts(k, x);
           60   return sum;
           61 }