ProteoWizard
LinearLeastSquares.hpp
Go to the documentation of this file.
1//
2// $Id$
3//
4// Original author: Robert Burke <robert.burke@cshs.org>
5//
6// Copyright 2006 Louis Warschaw Prostate Cancer Center
7// Cedars Sinai Medical Center, Los Angeles, California 90048
8//
9// Licensed under the Apache License, Version 2.0 (the "License");
10// you may not use this file except in compliance with the License.
11// You may obtain a copy of the License at
12//
13// http://www.apache.org/licenses/LICENSE-2.0
14//
15// Unless required by applicable law or agreed to in writing, software
16// distributed under the License is distributed on an "AS IS" BASIS,
17// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
18// See the License for the specific language governing permissions and
19// limitations under the License.
20//
21
22#ifndef _LINEARLEASTSQUARES_HPP_
23#define _LINEARLEASTSQUARES_HPP_
24
25
27#include <iostream>
28#include "LinearSolver.hpp"
29#include "Types.hpp"
30
31namespace pwiz {
32namespace math {
33
34enum PWIZ_API_DECL LinearLeastSquaresType {LinearLeastSquaresType_LU, LinearLeastSquaresType_QR};
35
36template <LinearLeastSquaresType lls_type = LinearLeastSquaresType_LU>
38
39template<>
41{
42public:
43 template<typename T>
44 boost::numeric::ublas::vector<T> solve(const boost::numeric::ublas::matrix<T>& A,
45 const boost::numeric::ublas::vector<T>& y)
46 {
47 boost::numeric::ublas::permutation_matrix<std::size_t> m(A.size1());
48 boost::numeric::ublas::matrix<T> AtA = prod(trans(A), A);
49 boost::numeric::ublas::vector<T> b = y;
50 boost::numeric::ublas::vector<T> r;
51
52 // This serves as a sanity check. Note that an exception here
53 // probably indicates a data file error.
54 if (boost::numeric::ublas::lu_factorize(AtA, m) == 0.)
55 {
56 r = prod(trans(A), b);
57
58 boost::numeric::ublas::lu_substitute(AtA, m, r);
59 }
60
61 return r;
62 }
63};
64
65template<>
66class LinearLeastSquares<LinearLeastSquaresType_QR>
67{
68public:
69
70 template<typename T>
71 boost::numeric::ublas::vector<T> solve(
72 const boost::numeric::ublas::matrix<T>& A,
73 const boost::numeric::ublas::vector<T>& x)
74 {
76
77 boost::numeric::ublas::vector<T> y = solver.solve(A, x);
78
79 return y;
80 }
81};
82
83}
84}
85
86#endif // _LINEARLEASTSQUARES_HPP_
#define PWIZ_API_DECL
Definition Export.hpp:32
#define A
LinearLeastSquaresType_LU
KernelTraitsBase< Kernel >::space_type::abscissa_type x
KernelTraitsBase< Kernel >::space_type::ordinate_type y
boost::numeric::ublas::vector< T > solve(const boost::numeric::ublas::matrix< T > &A, const boost::numeric::ublas::vector< T > &y)
boost::numeric::ublas::vector< T > solve(const boost::numeric::ublas::matrix< T > &A, const boost::numeric::ublas::vector< T > &x)