ProteoWizard
PeakDetectorMatchedFilterTest.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
28
29using namespace pwiz::util;
30using namespace pwiz::frequency;
31using namespace pwiz::chemistry;
32using namespace pwiz::data;
33using namespace pwiz::data::peakdata;
34using std::norm;
35using std::polar;
36
37
38ostream* os_ = 0;
39
40
42{
43 for (TestDatum* datum=testData_; datum<testData_+testDataSize_; datum++)
44 fd.data().push_back(FrequencyDatum(datum->frequency,
45 complex<double>(datum->real, datum->imaginary)));
46
49 fd.analyze();
50}
51
52
53void testCreation(const IsotopeEnvelopeEstimator& isotopeEnvelopeEstimator)
54{
55 if (os_) *os_ << "testCreation()\n";
56 const int filterMatchRate = 4;
57 const int filterSampleRadius = 2;
58 const double peakThresholdFactor = 0;
59 const double peakMaxCorrelationAngle = 5;
60 const double isotopeThresholdFactor = 666;
61 const double monoisotopicPeakThresholdFactor = 777;
62
64 config.isotopeEnvelopeEstimator = &isotopeEnvelopeEstimator;
65 config.filterMatchRate = filterMatchRate;
66 config.filterSampleRadius = filterSampleRadius;
67 config.peakThresholdFactor = peakThresholdFactor;
68 config.peakMaxCorrelationAngle = peakMaxCorrelationAngle;
69 config.isotopeThresholdFactor = isotopeThresholdFactor;
70 config.monoisotopicPeakThresholdFactor = monoisotopicPeakThresholdFactor;
71
72 auto_ptr<PeakDetectorMatchedFilter> pd = PeakDetectorMatchedFilter::create(config);
73
74 unit_assert(pd->config().filterMatchRate == filterMatchRate);
75 unit_assert(pd->config().filterSampleRadius == filterSampleRadius);
76 unit_assert(pd->config().peakThresholdFactor == peakThresholdFactor);
77 unit_assert(pd->config().peakMaxCorrelationAngle == peakMaxCorrelationAngle);
78 unit_assert(pd->config().isotopeThresholdFactor == isotopeThresholdFactor);
79 unit_assert(pd->config().monoisotopicPeakThresholdFactor == monoisotopicPeakThresholdFactor);
80}
81
82
83void testFind(FrequencyData& fd, const IsotopeEnvelopeEstimator& isotopeEnvelopeEstimator)
84{
85 if (os_) *os_ << "testFind()\n";
86
87 // fill in config structure
89 config.isotopeEnvelopeEstimator = &isotopeEnvelopeEstimator;
90 config.filterMatchRate = 4;
91 config.filterSampleRadius = 2;
92 config.peakThresholdFactor = 2;
93 config.peakMaxCorrelationAngle = 30;
94 config.isotopeThresholdFactor = 2;
96 config.isotopeMaxChargeState = 6;
97 config.isotopeMaxNeutronCount = 4;
98 config.collapseRadius = 15;
99 config.useMagnitudeFilter = false;
100 config.logDetailLevel = 1;
101 config.log = os_;
102
103 // instantiate
104 auto_ptr<PeakDetectorMatchedFilter> pd = PeakDetectorMatchedFilter::create(config);
105
106 // find peaks
107 PeakData data;
108 data.scans.push_back(Scan());
109 vector<PeakDetectorMatchedFilter::Score> scores;
110 pd->findPeaks(fd, data.scans[0], scores);
111
112 // report results
113 if (os_)
114 {
115 *os_ << "peaks found: " << data.scans[0].peakFamilies.size() << endl;
116 *os_ << data.scans[0];
117 *os_ << "scores: " << scores.size() << endl;
118 copy(scores.begin(), scores.end(),
119 ostream_iterator<PeakDetectorMatchedFilter::Score>(*os_, "\n"));
120 }
121
122 // assertions
123 unit_assert(data.scans[0].peakFamilies.size() == 1);
124 const PeakFamily& peakFamily = data.scans[0].peakFamilies.back();
125
126 if (os_) *os_ << "peakFamily: " << peakFamily << endl;
127 unit_assert(peakFamily.peaks.size() > 1);
128 const Peak& peak = peakFamily.peaks[0];
130 unit_assert(peakFamily.charge == 2);
131
132 unit_assert(scores.size() == 1);
133 const PeakDetectorMatchedFilter::Score& score = scores.back();
134 unit_assert(score.charge == peakFamily.charge);
137 polar(peak.intensity, peak.getAttribute(Peak::Attribute_Phase))),
138 0, 1e-14);
139}
140
141
142auto_ptr<IsotopeEnvelopeEstimator> createIsotopeEnvelopeEstimator()
143{
144 const double abundanceCutoff = .01;
145 const double massPrecision = .1;
146 IsotopeCalculator isotopeCalculator(abundanceCutoff, massPrecision);
147
149 config.isotopeCalculator = &isotopeCalculator;
150
151 return auto_ptr<IsotopeEnvelopeEstimator>(new IsotopeEnvelopeEstimator(config));
152}
153
154
155void test()
156{
157 if (os_) *os_ << setprecision(12);
158
159 auto_ptr<IsotopeEnvelopeEstimator> isotopeEnvelopeEstimator = createIsotopeEnvelopeEstimator();
160
161 testCreation(*isotopeEnvelopeEstimator);
162
163 FrequencyData fd;
165
166 testFind(fd, *isotopeEnvelopeEstimator);
167}
168
169
170int main(int argc, char* argv[])
171{
172 TEST_PROLOG(argc, argv)
173
174 try
175 {
176 if (argc>1 && !strcmp(argv[1],"-v")) os_ = &cout;
177 if (os_) *os_ << "PeakDetectorMatchedFilterTest\n";
178 test();
179 }
180 catch (exception& e)
181 {
182 TEST_FAILED(e.what())
183 }
184 catch (...)
185 {
186 TEST_FAILED("Caught unknown exception.")
187 }
188
190}
191
int main(int argc, char *argv[])
auto_ptr< IsotopeEnvelopeEstimator > createIsotopeEnvelopeEstimator()
void initializeWithTestData(FrequencyData &fd)
TestDatum testData_[]
const double testDataObservationDuration_
const unsigned int testDataSize_
const double testDataCalibrationB_
const double testDataCalibrationA_
void testCreation()
void testFind()
Class used for calculating a theoretical isotope envelope for a given mass, based on an estimate of t...
Class for binary storage of complex frequency data.
double observationDuration() const
const CalibrationParameters & calibrationParameters() const
void analyze()
recache statistics calculations after any direct data changes via non-const data()
const container & data() const
const access to underlying data
static std::auto_ptr< PeakDetectorMatchedFilter > create(const Config &config)
create an instance.
SampleDatum< double, std::complex< double > > FrequencyDatum
std::vector< Scan > scans
Definition PeakData.hpp:194
double getAttribute(Attribute attribute) const
int logDetailLevel
log detail level (0 == normal, 1 == extra)
double peakMaxCorrelationAngle
maximum correlation angle (degrees) for initial peak reporting
int isotopeMaxChargeState
isotope filter maximum charge state to score
const chemistry::IsotopeEnvelopeEstimator * isotopeEnvelopeEstimator
IsotopeEnvelopeEstimator pointer, must be valid for PeakDetector lifetime.
int filterMatchRate
number of filter correlations computed per frequency step
double monoisotopicPeakThresholdFactor
noise floor multiple for monoisotopic peak threshold
double peakThresholdFactor
noise floor multiple for initial peak reporting threshold
double collapseRadius
multiple peaks within this radius (Hz) are reported as single peak
int filterSampleRadius
number of filter samples taken on either side of 0
double isotopeThresholdFactor
noise floor multiple for isotope filter threshold
bool useMagnitudeFilter
use the magnitude of the peak shape filter kernel for finding peaks
int isotopeMaxNeutronCount
isotope filter maximum number of neutrons to score
structure for holding the matched filter calculation results
#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