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