ProteoWizard
Main Page
Namespaces
Classes
Files
File List
File Members
pwiz
utility
math
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
"
23
#include "
pwiz/utility/misc/unit.hpp
"
24
#include "
pwiz/utility/misc/Std.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
= 1
e
-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., 1
e
-13);
60
unit_assert_equal
(
x
(1), 2., 1
e
-13);
61
}
62
63
void
testDoubleQR
()
64
{
65
if
(
os_
) *
os_
<<
"testDoubleQR()\n"
;
66
67
LinearLeastSquares<LinearLeastSquaresType_QR>
lls;
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
87
void
testExactFitQR
()
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
103
LinearLeastSquares<LinearLeastSquaresType_QR>
lls;
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
113
void
testSimpleRectangleQR
()
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
136
LinearLeastSquares<LinearLeastSquaresType_QR>
lls;
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
146
void
testLeastSquaresQR
()
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
172
LinearLeastSquares<LinearLeastSquaresType_QR>
lls;
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
();
193
testSimpleRectangleQR
();
194
testLeastSquaresQR
();
195
}
196
catch
(exception&
e
)
197
{
198
TEST_FAILED
(e.what())
199
}
200
catch
(...)
201
{
202
TEST_FAILED
(
"Caught unknown exception."
)
203
}
204
205
TEST_EPILOG
206
}
Generated on Mon Nov 26 2012 18:05:49 for ProteoWizard by
1.8.1.1