ProteoWizard
LinearLeastSquaresTest.cpp
Go to the documentation of this file.
1 //
2 // $Id: LinearLeastSquaresTest.cpp 4140 2012-11-22 01:07:16Z pcbrefugee $
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 #include "LinearLeastSquares.hpp"
25 #include <cstring>
26 
27 
28 using namespace pwiz::util;
29 using namespace pwiz::math;
30 using namespace pwiz::math::types;
31 
32 
33 namespace ublas = boost::numeric::ublas;
34 
35 
36 const double epsilon = 1e-16;
37 ostream* os_ = 0;
38 
39 
40 void testDouble()
41 {
42  if (os_) *os_ << "testDouble()\n";
43 
44  LinearLeastSquares<> lls;
45  ublas::matrix<double> A(2, 2);
46  A(0,0) = 1; A(0,1) = 2;
47  A(1,0) = 3; A(1,1) = 4;
48 
49  ublas::vector<double> y(2);
50  y(0) = 5;
51  y(1) = 11;
52 
53  ublas::vector<double> x = lls.solve(A, y);
54 
55  if (os_) *os_ << "A: " << A << endl;
56  if (os_) *os_ << "y: " << y << endl;
57  if (os_) *os_ << "x: " << x << endl;
58 
59  unit_assert_equal(x(0), 1., 1e-13);
60  unit_assert_equal(x(1), 2., 1e-13);
61 }
62 
64 {
65  if (os_) *os_ << "testDoubleQR()\n";
66 
68  ublas::matrix<double> A(2, 2);
69  A(0,0) = 1; A(0,1) = 2;
70  A(1,0) = 3; A(1,1) = 4;
71 
72  ublas::vector<double> y(2);
73  y(0) = 5;
74  y(1) = 11;
75 
76  ublas::vector<double> x = lls.solve(A, y);
77 
78  if (os_) *os_ << "A: " << A << endl;
79  if (os_) *os_ << "y: " << y << endl;
80  if (os_) *os_ << "x: " << x << endl;
81 
82  if (os_) *os_ << "x(0) = " << x(0) - 1. << ", x(1) = " << x(1) - 2. << endl;
83  unit_assert_equal(x(0), 1., 100*epsilon);
84  unit_assert_equal(x(1), 2., 100*epsilon);
85 }
86 
88 {
89  if (os_) *os_ << "***************************\n";
90  if (os_) *os_ << "testExactFit()\n";
91 
92  dmatrix m(2,2);
93  dvector obs(2);
94 
95  m(0, 0) = 1;
96  m(1, 0) = 0;
97  m(0, 1) = 0;
98  m(1, 1) = 1;
99 
100  obs(0) = 1;
101  obs(1) = 1;
102 
104  const dvector result = lls.solve(m,obs);
105 
106  unit_assert_equal(obs(0), 1, epsilon);
107  unit_assert_equal(obs(1), 1, epsilon);
108 
109  if (os_) *os_ << "testExactFit(): success\n";
110 }
111 
112 
114 {
115  if (os_) *os_ << "***************************\n";
116  if (os_) *os_ << "testSimpleRectangleQR()\n";
117 
118  dmatrix samples(4, 2);
119 
120  samples.clear();
121  samples(0,0) = 1;
122  samples(0,1) = 0;
123  samples(1,0) = 0;
124  samples(1,1) = 1;
125  samples(2,0) = 1;
126  samples(2,1) = 0;
127  samples(3,0) = 0;
128  samples(3,1) = 1;
129 
130  dvector obs(4);
131  obs(0) = 2;
132  obs(1) = 2;
133  obs(2) = 2;
134  obs(3) = 2;
135 
137 
138  dvector a = lls.solve(samples, obs);
139 
140  unit_assert_equal(a(0), 2, epsilon);
141  unit_assert_equal(a(1), 2, epsilon);
142 
143  if (os_) *os_ << "testSimpleRectangleQR(): success\n";
144 }
145 
147 {
148  if (os_) *os_ << "***************************\n";
149  if (os_) *os_ << "testLeastSquaresQR()\n";
150 
151  dmatrix samples(5, 2);
152 
153  samples.clear();
154  samples(0,0) = 1;
155  samples(0,1) = 1;
156  samples(1,0) = 2;
157  samples(1,1) = 2;
158  samples(2,0) = 3;
159  samples(2,1) = 3;
160  samples(3,0) = 0;
161  samples(3,1) = 4;
162  samples(4,0) = -1;
163  samples(4,1) = 5;
164 
165  dvector obs(5);
166  obs(0) = 1;
167  obs(1) = 3;
168  obs(2) = 9;
169  obs(3) = 3;
170  obs(4) = -9;
171 
173 
174  dvector a = lls.solve(samples, obs);
175 
176  unit_assert_equal(a(0), 3.16666666666667, epsilon*100);
177  unit_assert_equal(a(1), -0.5, epsilon*100);
178 
179  if (os_) *os_ << "testLeastSquaresQR(): success\n";
180 }
181 
182 int main(int argc, char* argv[])
183 {
184  TEST_PROLOG(argc, argv)
185 
186  try
187  {
188  if (argc>1 && !strcmp(argv[1],"-v")) os_ = &cout;
189  if (os_) *os_ << "LinearLeastSquaresTest\n";
190  testDouble();
191  testDoubleQR();
192  testExactFitQR();
195  }
196  catch (exception& e)
197  {
198  TEST_FAILED(e.what())
199  }
200  catch (...)
201  {
202  TEST_FAILED("Caught unknown exception.")
203  }
204 
206 }