ProteoWizard
FrequencyEstimatorSimpleTest.cpp
Go to the documentation of this file.
1//
2// $Id$
3//
4//
5// Original author: Darren Kessner <darren@proteowizard.org>
6//
7// Copyright 2006 Louis Warschaw Prostate Cancer Center
8// Cedars Sinai Medical Center, Los Angeles, California 90048
9//
10// Licensed under the Apache License, Version 2.0 (the "License");
11// you may not use this file except in compliance with the License.
12// You may obtain a copy of the License at
13//
14// http://www.apache.org/licenses/LICENSE-2.0
15//
16// Unless required by applicable law or agreed to in writing, software
17// distributed under the License is distributed on an "AS IS" BASIS,
18// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
19// See the License for the specific language governing permissions and
20// limitations under the License.
21//
22
23
27
28using namespace pwiz::util;
29using namespace pwiz::frequency;
30using namespace pwiz::data;
31using namespace pwiz::data::peakdata;
32
33
34ostream* os_ = 0;
35
36
37void test()
38{
39 if (os_) *os_ << "**************************************************\ntest()\n";
40
41 // create data with Lorentzian peak at 0, parabolic peak at 7.5
43 for (int i=-5; i<=5; i++)
44 fd.data().push_back(FrequencyDatum(i, 3/sqrt(i*i+1.)));
45 fd.data().push_back(FrequencyDatum(6, 1));
46 fd.data().push_back(FrequencyDatum(7, 9));
47 fd.data().push_back(FrequencyDatum(8, 9));
48 fd.data().push_back(FrequencyDatum(9, 1));
49
50 // create initial peak list
51 vector<Peak> peaks(2);
52 peaks[0].attributes[Peak::Attribute_Frequency] = .1;
53 peaks[1].attributes[Peak::Attribute_Frequency] = 7.4;
54
55 // storage for results
56
57 vector<Peak> estimatedPeaks;
58
59 // run Parabola estimator
60
61 auto_ptr<FrequencyEstimatorSimple>
62 fe(FrequencyEstimatorSimple::create(FrequencyEstimatorSimple::Parabola));
63
64 for (vector<Peak>::const_iterator it=peaks.begin(); it!=peaks.end(); ++it)
65 estimatedPeaks.push_back(fe->estimate(fd, *it));
66
67 // check results
68
69 if (os_)
70 {
71 *os_ << setprecision(10);
72
73 *os_ << "Initial peaks:\n";
74 copy(peaks.begin(), peaks.end(), ostream_iterator<Peak>(*os_, "\n"));
75 *os_ << endl;
76
77 *os_ << "Parabola estimated peaks:\n";
78 copy(estimatedPeaks.begin(), estimatedPeaks.end(), ostream_iterator<Peak>(*os_, "\n"));
79 *os_ << endl;
80 }
81
82 unit_assert(estimatedPeaks.size() == 2);
83
84 const Peak& pi0 = estimatedPeaks[0];
86 unit_assert(pi0.intensity == 3.);
87
88 const Peak& pi1 = estimatedPeaks[1];
90 unit_assert_equal(abs(pi1.intensity-10.), 0., 1e-10);
91
92
93 // run Lorentzian estimator
94
95 estimatedPeaks.clear();
96 fe = FrequencyEstimatorSimple::create(FrequencyEstimatorSimple::Lorentzian);
97
98 for (vector<Peak>::const_iterator it=peaks.begin(); it!=peaks.end(); ++it)
99 estimatedPeaks.push_back(fe->estimate(fd, *it));
100
101 // check results
102
103 if (os_)
104 {
105 *os_ << "Lorentzian estimated peaks:\n";
106 copy(estimatedPeaks.begin(), estimatedPeaks.end(), ostream_iterator<Peak>(*os_, "\n"));
107 *os_ << endl;
108 }
109
110 unit_assert(estimatedPeaks.size() == 2);
111
112 const Peak& pil0 = estimatedPeaks[0];
114 unit_assert_equal(abs(pil0.intensity-3.), 0, 1e-10);
115
116 const Peak& pil1 = estimatedPeaks[1];
118 unit_assert(pil1.intensity == 0.); // intensity is nan
119}
120
121
123{
124 if (os_) *os_ << "**************************************************\ntestData()\n";
125
126 FrequencyData fd;
127 fd.data().push_back(FrequencyDatum(28558.59375, complex<double>(25243.032972361, -2820.6360692452)));
128 fd.data().push_back(FrequencyDatum(28559.895833333, complex<double>(39978.141686921, 291.1363106641)));
129 fd.data().push_back(FrequencyDatum(28561.197916667, complex<double>(189200.35822792, -2636.9254689346)));
130 fd.data().push_back(FrequencyDatum(28562.5, complex<double>(-62230.480432624, -2546.1033855971)));
131 fd.data().push_back(FrequencyDatum(28563.802083333, complex<double>(-32263.08735743, -2769.7946573836)));
132
133 vector<Peak> peaks(1);
134 peaks[0].attributes[Peak::Attribute_Frequency] = 28561.2;
135
136 if (os_)
137 {
138 *os_ << "Initial peaks:\n";
139 copy(peaks.begin(), peaks.end(), ostream_iterator<Peak>(*os_, "\n"));
140 *os_ << endl;
141 }
142
143 vector<Peak> estimatedPeaks;
144
145 auto_ptr<FrequencyEstimatorSimple>
146 fe(FrequencyEstimatorSimple::create(FrequencyEstimatorSimple::Parabola));
147
148 for (vector<Peak>::const_iterator it=peaks.begin(); it!=peaks.end(); ++it)
149 estimatedPeaks.push_back(fe->estimate(fd, *it));
150
151 if (os_)
152 {
153 *os_ << setprecision(10);
154
155 *os_ << "Parabola estimated peaks:\n";
156 copy(estimatedPeaks.begin(), estimatedPeaks.end(), ostream_iterator<Peak>(*os_, "\n"));
157 *os_ << endl;
158 }
159
160 unit_assert(estimatedPeaks.size() == 1);
161 unit_assert_equal(estimatedPeaks[0].getAttribute(Peak::Attribute_Frequency), 28561.25049, 1e-5);
162}
163
164
165struct Datum
166{
167 double frequency;
168 complex<double> intensity;
169};
170
171
173{
174 {10981.7708333, complex<double>(1494.77333435,-827.565405375)},
175 {10982.421875, complex<double>(-522.951336943,-2933.77290646)},
176 {10983.0729167, complex<double>(1026.58070488,-2790.70883417)},
177 {10983.7239583, complex<double>(1002.03072708,-1020.67745139)},
178 {10984.375, complex<double>(-567.573503924,-2220.62261993)},
179 {10985.0260417, complex<double>(-6322.94426498,-3013.78424791)},
180 {10985.6770833, complex<double>(430.465272274,502.150355144)},
181 {10986.328125, complex<double>(2578.81032322,795.379729653)},
182 {10986.9791667, complex<double>(2864.69277204,-470.140696311)},
183 {10987.6302083, complex<double>(2788.00641762,4788.24971282)},
184 {10988.28125, complex<double>(-366.077646703,-6084.91428783)},
185 {10988.9322917, complex<double>(1220.81029308,-3297.88016503)},
186 {10989.5833333, complex<double>(2268.72858986,-646.091997391)},
187 {10990.234375, complex<double>(4681.74708664,2313.31976782)},
188 {10990.8854167, complex<double>(-955.7765424,5903.76925847)},
189 {10991.5364583, complex<double>(3957.71316667,-225.389114599)},
190 {10992.1875, complex<double>(-2159.60123121,3597.12682291)},
191 {10992.8385417, complex<double>(653.493128029,2229.46593497)},
192 {10993.4895833, complex<double>(6037.67518189,-4347.22639235)},
193 {10994.140625, complex<double>(-167.455004321,1848.75455373)}
194};
195
196
197const int dataSize_ = sizeof(data_)/sizeof(Datum);
198
199
201{
202 if (os_) *os_ << "**************************************************\ntestData2_LocalMax()\n";
203
204 FrequencyData fd;
205
206 for (const Datum* p=data_; p!=data_+dataSize_; ++p)
207 fd.data().push_back(FrequencyDatum(p->frequency, p->intensity));
208
209 vector<Peak> peaks(1);
210 peaks[0].attributes[Peak::Attribute_Frequency] = 10983.74;
211
212 if (os_)
213 {
214 *os_ << "Initial peaks:\n";
215 copy(peaks.begin(), peaks.end(), ostream_iterator<Peak>(*os_, "\n"));
216 *os_ << endl;
217 }
218
219 vector<Peak> estimatedPeaks;
220
221 auto_ptr<FrequencyEstimatorSimple>
222 fe(FrequencyEstimatorSimple::create(FrequencyEstimatorSimple::LocalMax));
223
224 for (vector<Peak>::const_iterator it=peaks.begin(); it!=peaks.end(); ++it)
225 estimatedPeaks.push_back(fe->estimate(fd, *it));
226
227 if (os_)
228 {
229 *os_ << setprecision(10);
230
231 *os_ << "Local max peaks:\n";
232 copy(estimatedPeaks.begin(), estimatedPeaks.end(), ostream_iterator<Peak>(*os_, "\n"));
233 *os_ << endl;
234 }
235
236 unit_assert(estimatedPeaks.size() == 1);
237 unit_assert_equal(estimatedPeaks[0].getAttribute(Peak::Attribute_Frequency), 10985.02604, 1e-5);
238}
239
240
242{
243 if (os_) *os_ << "**************************************************\ntestData2()\n";
244
245 FrequencyData fd;
246
247 for (const Datum* p=data_; p!=data_+dataSize_; ++p)
248 fd.data().push_back(FrequencyDatum(p->frequency, p->intensity));
249
250 vector<Peak> peaks(1);
251 peaks[0].attributes[Peak::Attribute_Frequency] = 10987.6;
252
253 if (os_)
254 {
255 *os_ << "Initial peaks:\n";
256 copy(peaks.begin(), peaks.end(), ostream_iterator<Peak>(*os_, "\n"));
257 *os_ << endl;
258 }
259
260 vector<Peak> estimatedPeaks;
261
262 auto_ptr<FrequencyEstimatorSimple>
263 fe(FrequencyEstimatorSimple::create(FrequencyEstimatorSimple::Parabola));
264
265 for (vector<Peak>::const_iterator it=peaks.begin(); it!=peaks.end(); ++it)
266 estimatedPeaks.push_back(fe->estimate(fd, *it));
267
268 if (os_)
269 {
270 *os_ << setprecision(10);
271
272 *os_ << "Parabola peaks:\n";
273 copy(estimatedPeaks.begin(), estimatedPeaks.end(), ostream_iterator<Peak>(*os_, "\n"));
274 *os_ << endl;
275 }
276
277 unit_assert(estimatedPeaks.size() == 1);
278 unit_assert_equal(estimatedPeaks[0].getAttribute(Peak::Attribute_Frequency), 10988.07103, 1e-5);
279}
280
281
282int main(int argc, char* argv[])
283{
284 TEST_PROLOG(argc, argv)
285
286 try
287 {
288 if (argc>1 && !strcmp(argv[1],"-v")) os_ = &cout;
289 if (os_) *os_ << "FrequencyEstimatorSimpleTest\n";
290 test();
291 testData();
294 }
295 catch (exception& e)
296 {
297 TEST_FAILED(e.what())
298 }
299 catch (...)
300 {
301 TEST_FAILED("Caught unknown exception.")
302 }
303
305}
306
int main(int argc, char *argv[])
void testData2_LocalMax()
void testData2_Parabola()
const int dataSize_
Class for binary storage of complex frequency data.
const container & data() const
const access to underlying data
static std::auto_ptr< FrequencyEstimatorSimple > create(Type type=Parabola, unsigned int windowRadius=1)
create an instance
SampleDatum< double, std::complex< double > > FrequencyDatum
double getAttribute(Attribute attribute) const
#define unit_assert(x)
Definition unit.hpp:85
#define TEST_EPILOG
Definition unit.hpp:183
#define TEST_FAILED(x)
Definition unit.hpp:177
#define unit_assert_equal(x, y, epsilon)
Definition unit.hpp:99
#define TEST_PROLOG(argc, argv)
Definition unit.hpp:175