fit_linear.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2013-2017, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 /* multifit/multilinear.c
24  *
25  * Copyright (C) 2000, 2007, 2010 Brian Gough
26  * Copyright (C) 2013, 2015 Patrick Alken
27  *
28  * This program is free software; you can redistribute it and/or modify
29  * it under the terms of the GNU General Public License as published by
30  * the Free Software Foundation; either version 3 of the License, or (at
31  * your option) any later version.
32  *
33  * This program is distributed in the hope that it will be useful, but
34  * WITHOUT ANY WARRANTY; without even the implied warranty of
35  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
36  * General Public License for more details.
37  *
38  * You should have received a copy of the GNU General Public License
39  * along with this program; if not, write to the Free Software
40  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
41  * 02110-1301, USA.
42  */
43 #ifndef O2SCL_FIT_LINEAR_H
44 #define O2SCL_FIT_LINEAR_H
45 
46 /** \file fit_linear.h
47  \brief File defining \ref o2scl::fit_linear
48 */
49 
50 #include <boost/numeric/ublas/vector.hpp>
51 #include <boost/numeric/ublas/matrix.hpp>
52 
53 #include <o2scl/vector.h>
54 #include <o2scl/fit_base.h>
55 #include <o2scl/permutation.h>
56 #include <o2scl/cblas.h>
57 #include <o2scl/svd.h>
58 
59 #ifndef DOXYGEN_NO_O2NS
60 namespace o2scl {
61 #endif
62 
63  /** \brief Linear least-squares fitting class (GSL)
64  */
65  template<class vec_t=boost::numeric::ublas::vector<double>,
66  class mat_t=boost::numeric::ublas::matrix<double> >
67  class fit_linear {
68 
69  protected:
70 
71  /// \name Storage
72  //@{
73  /** \brief Local copy of <tt>xpred</tt>
74  (and used as workspace by the SV decomposition)
75  */
76  mat_t A;
77  /// The first unitary matrix from the SV decomposition
78  mat_t Q;
79  /// Workspace for the SV decomposition and storage for \f$ Q S^{-1} \f$
80  mat_t QSI;
81 
82  /// The singular values from the SV decomposition
83  vec_t S;
84  /// SV decomposition workspace and also used to store new predicted values
85  vec_t xt;
86  /// Balancing factors for A
87  vec_t D;
88  /// Only used for the weighted fit (not yet implemented)
89  vec_t t;
90  //@}
91 
92  /// Number of parameters
93  size_t size_par;
94  /// Number of data points
95  size_t size_dat;
96 
97  public:
98 
99  fit_linear() {
100  column_scaling=true;
101  tol=2.2204460492503131e-16;
102  size_par=0;
103  size_dat=0;
104  }
105 
106  virtual ~fit_linear() {
107  }
108 
109  /** \brief If true, discard fit components if the associated singular
110  value becomes too small (default true)
111  */
113 
114  /// Tolerance (default \f$ \sim 2.22\times 10^{-16} \f$)
115  double tol;
116 
117  /** \brief The rank of the linear system from the last call to
118  <tt>fit_linear()</tt>
119  */
120  size_t rank;
121 
122  /** \brief Perform the SV decomposition
123  */
124  virtual void fit_svd(size_t ndat, size_t npar) {
125  o2scl_linalg::SV_decomp_mod(ndat,npar,A,QSI,Q,S,xt);
126  return;
127  }
128 
129  /** \brief Perform a least-squares fit of a linear system
130 
131  This function performs a least-squares fit of the
132  system
133  \f[
134  \mathrm{xpred} \cdot \mathrm{parms} = \mathrm{ydat}
135  \f]
136 
137  The variance-covariance matrix for the parameters is returned in
138  \c covar and the value of \f$ \chi^2 \f$ is returned in \c chi2.
139  */
140  virtual void fit(size_t npar, size_t ndat, const vec_t &ydat,
141  const mat_t &xpred, vec_t &parms,
142  mat_t &covar, double &chi2) {
143 
144  if (npar>ndat) {
145  O2SCL_ERR2("Insufficient data points in ",
146  "fit_linear::fit_linear().",exc_einval);
147  }
148 
149  if (npar!=size_par || ndat!=size_dat) {
150  A.resize(ndat,npar);
151  Q.resize(npar,npar);
152  QSI.resize(npar,npar);
153  S.resize(npar);
154  t.resize(ndat);
155  xt.resize(npar);
156  D.resize(npar);
157  size_par=npar;
158  size_dat=ndat;
159  }
160 
161  // necessary?
162  for(size_t i=0;i<npar;i++) {
163  xt[i]=0.0;
164  D[i]=0.0;
165  }
166 
167  for(size_t i=0;i<ndat;i++) {
168  for(size_t j=0;j<npar;j++) {
169  A(i,j)=xpred(i,j);
170  }
171  }
172 
173  if (column_scaling) {
174  o2scl_linalg::balance_columns(ndat,npar,A,D);
175  } else {
176  for(size_t i=0;i<npar;i++) {
177  D[i]=1.0;
178  }
179  }
180 
181  // [GSL] Decompose A into U S Q^T
182  fit_svd(ndat,npar);
183 
184  // [GSL] Solve y = A c for c
186  (o2scl_cblas::o2cblas_RowMajor,
187  o2scl_cblas::o2cblas_Trans,ndat,npar,1.0,A,ydat,0.0,xt);
188 
189  // [GSL] Scale the matrix Q, Q' = Q S^-1
190  for(size_t i=0;i<npar;i++) {
191  for(size_t j=0;j<npar;j++) {
192  QSI(i,j)=Q(i,j);
193  }
194  }
195 
196  double alpha0=S[0];
197  size_t p_eff=0;
198 
199  for(size_t j=0;j<npar;j++) {
200  double alpha=S[j];
201  if (alpha<=tol*alpha0) {
202  alpha=0.0;
203  } else {
204  alpha=1.0/alpha;
205  p_eff++;
206  }
207  o2scl_cblas::dscal_subcol(QSI,0,j,npar,alpha);
208  }
209 
210  rank=p_eff;
211  for(size_t i=0;i<npar;i++) {
212  parms[i]=0.0;
213  }
214 
216  (o2scl_cblas::o2cblas_RowMajor,
217  o2scl_cblas::o2cblas_NoTrans,npar,npar,
218  1.0,QSI,xt,0.0,parms);
219 
220  // [GSL] Unscale the balancing factors
221 
222  for(size_t i=0;i<npar;i++) {
223  parms[i]/=D[i];
224  }
225 
226  // [GSL] Compute chi-squared from residual, r = y - X c
227 
228  double s2=0.0, r2=0.0;
229  for(size_t i=0;i<ndat;i++) {
230  double yi=ydat[i];
231  double y_est=o2scl_cblas::ddot_subrow(npar,xpred,i,0,parms);
232  double ri=yi-y_est;
233  r2+=ri*ri;
234  }
235  s2=r2/(ndat-p_eff);
236  chi2=r2;
237 
238  // [GSL] Form variance-covariance matrix, cov = s2 * (Q S^-1) (Q S^-1)^T
239 
240  for(size_t i=0;i<npar;i++) {
241  double d_i=D[i];
242  for(size_t j=i;j<npar;j++) {
243  double d_j=D[j];
244  double s=0.0;
245  for(size_t k=0;k<npar;k++) {
246  s+=QSI(i,k)*QSI(j,k);
247  }
248  covar(i,j)=s*s2/(d_i*d_j);
249  covar(j,i)=s*s2/(d_i*d_j);
250  }
251  }
252 
253  return;
254  }
255 
256  /// Return string denoting type ("fit_linear")
257  virtual const char *type() { return "fit_linear"; }
258 
259  };
260 
261 #ifndef DOXYGEN_NO_O2NS
262 }
263 #endif
264 
265 #endif
size_t rank
The rank of the linear system from the last call to fit_linear()
Definition: fit_linear.h:120
The main O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl names...
Definition: anneal.h:42
double tol
Tolerance (default )
Definition: fit_linear.h:115
Linear least-squares fitting class (GSL)
Definition: fit_linear.h:67
size_t size_par
Number of parameters.
Definition: fit_linear.h:93
invalid argument supplied by user
Definition: err_hnd.h:59
void balance_columns(size_t M, size_t N, mat_t &A, vec_t &D)
Balance a general matrix A by scaling the columns by the diagonal matrix D.
Definition: svd_base.h:615
double ddot_subrow(const size_t N, const mat_t &X, const size_t ir, const size_t ic, const vec_t &Y)
Compute for a subrow of a matrix.
Definition: cblas_base.h:1223
size_t size_dat
Number of data points.
Definition: fit_linear.h:95
vec_t D
Balancing factors for A.
Definition: fit_linear.h:87
virtual void fit_svd(size_t ndat, size_t npar)
Perform the SV decomposition.
Definition: fit_linear.h:124
virtual void fit(size_t npar, size_t ndat, const vec_t &ydat, const mat_t &xpred, vec_t &parms, mat_t &covar, double &chi2)
Perform a least-squares fit of a linear system.
Definition: fit_linear.h:140
void dgemv(const enum o2cblas_order order, const enum o2cblas_transpose TransA, const size_t M, const size_t N, const double alpha, const mat_t &A, const vec_t &X, const double beta, vec2_t &Y)
Compute .
Definition: cblas_base.h:218
mat_t QSI
Workspace for the SV decomposition and storage for .
Definition: fit_linear.h:80
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
vec_t S
The singular values from the SV decomposition.
Definition: fit_linear.h:83
mat_t A
Local copy of xpred (and used as workspace by the SV decomposition)
Definition: fit_linear.h:76
mat_t Q
The first unitary matrix from the SV decomposition.
Definition: fit_linear.h:78
virtual const char * type()
Return string denoting type ("fit_linear")
Definition: fit_linear.h:257
void SV_decomp_mod(size_t M, size_t N, mat_t &A, mat2_t &X, mat3_t &V, vec_t &S, vec2_t &work)
SV decomposition by the modified Golub-Reinsch algorithm which is better for .
Definition: svd_base.h:290
vec_t t
Only used for the weighted fit (not yet implemented)
Definition: fit_linear.h:89
bool column_scaling
If true, discard fit components if the associated singular value becomes too small (default true) ...
Definition: fit_linear.h:112
void dscal_subcol(mat_t &A, const size_t ir, const size_t ic, const size_t M, const double alpha)
Compute for a subcolumn of a matrix.
Definition: cblas_base.h:1093
vec_t xt
SV decomposition workspace and also used to store new predicted values.
Definition: fit_linear.h:85

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).