ProteoWizard
PeakDetectorMatchedFilterTest.cpp
Go to the documentation of this file.
1 //
2 // $Id: PeakDetectorMatchedFilterTest.cpp 4129 2012-11-20 00:05:37Z chambm $
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 
29 using namespace pwiz::util;
30 using namespace pwiz::frequency;
31 using namespace pwiz::chemistry;
32 using namespace pwiz::data;
33 using namespace pwiz::data::peakdata;
34 using std::norm;
35 using std::polar;
36 
37 
38 ostream* 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 
53 void 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 
83 void 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];
129  unit_assert_equal(peak.getAttribute(Peak::Attribute_Frequency), 159455, 1);
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);
135  unit_assert(score.monoisotopicFrequency == peak.getAttribute(Peak::Attribute_Frequency));
136  unit_assert_equal(norm(score.monoisotopicIntensity -
137  polar(peak.intensity, peak.getAttribute(Peak::Attribute_Phase))),
138  0, 1e-14);
139 }
140 
141 
142 auto_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 
155 void 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 
170 int 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 
double collapseRadius
multiple peaks within this radius (Hz) are reported as single peak
TestDatum testData_[]
int filterSampleRadius
number of filter samples taken on either side of 0
Class for binary storage of complex frequency data.
int isotopeMaxNeutronCount
isotope filter maximum number of neutrons to score
int logDetailLevel
log detail level (0 == normal, 1 == extra)
Class used for calculating a theoretical isotope envelope for a given mass, based on an estimate of t...
void testFind(FrequencyData &fd, const IsotopeEnvelopeEstimator &isotopeEnvelopeEstimator)
#define TEST_EPILOG
Definition: unit.hpp:182
#define unit_assert_equal(x, y, epsilon)
Definition: unit.hpp:99
void initializeWithTestData(FrequencyData &fd)
const chemistry::IsotopeEnvelopeEstimator * isotopeEnvelopeEstimator
IsotopeEnvelopeEstimator pointer, must be valid for PeakDetector lifetime.
int main(int argc, char *argv[])
void analyze()
recache statistics calculations after any direct data changes via non-const data() ...
const double testDataCalibrationB_
virtual const Config & config() const =0
access to the configuration
const container & data() const
const access to underlying data
double monoisotopicPeakThresholdFactor
noise floor multiple for monoisotopic peak threshold
std::vector< Scan > scans
Definition: PeakData.hpp:194
const double testDataCalibrationA_
double observationDuration() const
std::ostream * log
log stream (0 == no logging)
void testCreation(const IsotopeEnvelopeEstimator &isotopeEnvelopeEstimator)
auto_ptr< IsotopeEnvelopeEstimator > createIsotopeEnvelopeEstimator()
double peakThresholdFactor
noise floor multiple for initial peak reporting threshold
int filterMatchRate
number of filter correlations computed per frequency step
double peakMaxCorrelationAngle
maximum correlation angle (degrees) for initial peak reporting
int isotopeMaxChargeState
isotope filter maximum charge state to score
SampleDatum< double, std::complex< double > > FrequencyDatum
const CalibrationParameters & calibrationParameters() const
#define TEST_FAILED(x)
Definition: unit.hpp:176
const double testDataObservationDuration_
virtual void findPeaks(const pwiz::data::FrequencyData &fd, pwiz::data::peakdata::Scan &result) const =0
Find the peaks in the frequency data, filling in Scan structure.
double isotopeThresholdFactor
noise floor multiple for isotope filter threshold
bool useMagnitudeFilter
use the magnitude of the peak shape filter kernel for finding peaks
#define TEST_PROLOG(argc, argv)
Definition: unit.hpp:174
structure for holding the matched filter calculation results
#define unit_assert(x)
Definition: unit.hpp:85
const unsigned int testDataSize_