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 
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));
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 
PeakDetectorMatchedFilterTestData.hpp
PeakDetectorMatchedFilter.hpp
pwiz::frequency::PeakDetectorMatchedFilter::Score::monoisotopicIntensity
std::complex< double > monoisotopicIntensity
Definition: PeakDetectorMatchedFilter.hpp:137
pwiz::data::peakdata::PeakData
Definition: PeakData.hpp:190
pwiz::data::CalibrationParameters
Definition: CalibrationParameters.hpp:45
pwiz::data::FrequencyData::data
const container & data() const
const access to underlying data
pwiz::frequency::PeakDetectorMatchedFilter::Config::log
std::ostream * log
log stream (0 == no logging)
Definition: PeakDetectorMatchedFilter.hpp:85
pwiz::data::peakdata::Peak
Definition: PeakData.hpp:52
pwiz::frequency
Definition: DerivativeTest.hpp:37
unit_assert_equal
#define unit_assert_equal(x, y, epsilon)
Definition: unit.hpp:99
testDataCalibrationB_
const double testDataCalibrationB_
Definition: PeakDetectorMatchedFilterTestData.hpp:44
pwiz::data::FrequencyData::observationDuration
double observationDuration() const
pwiz::frequency::PeakDetectorMatchedFilter::findPeaks
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.
pwiz::frequency::PeakDetectorMatchedFilter::Config::collapseRadius
double collapseRadius
multiple peaks within this radius (Hz) are reported as single peak
Definition: PeakDetectorMatchedFilter.hpp:76
testCreation
void testCreation(const IsotopeEnvelopeEstimator &isotopeEnvelopeEstimator)
Definition: PeakDetectorMatchedFilterTest.cpp:53
createIsotopeEnvelopeEstimator
auto_ptr< IsotopeEnvelopeEstimator > createIsotopeEnvelopeEstimator()
Definition: PeakDetectorMatchedFilterTest.cpp:142
pwiz::data::peakdata::PeakFamily
Definition: PeakData.hpp:111
pwiz::data::peakdata
Definition: PeakData.hpp:40
pwiz::data::peakdata::Scan
Definition: PeakData.hpp:134
pwiz::data
Definition: BinaryIndexStream.hpp:31
pwiz::frequency::PeakDetectorMatchedFilter::Config::isotopeMaxChargeState
int isotopeMaxChargeState
isotope filter maximum charge state to score
Definition: PeakDetectorMatchedFilter.hpp:70
pwiz::data::FrequencyData
Class for binary storage of complex frequency data.
Definition: FrequencyData.hpp:48
pwiz::frequency::PeakDetectorMatchedFilter::Config::logDetailLevel
int logDetailLevel
log detail level (0 == normal, 1 == extra)
Definition: PeakDetectorMatchedFilter.hpp:82
testDataObservationDuration_
const double testDataObservationDuration_
Definition: PeakDetectorMatchedFilterTestData.hpp:42
pwiz::chemistry::IsotopeEnvelopeEstimator
Class used for calculating a theoretical isotope envelope for a given mass, based on an estimate of t...
Definition: IsotopeEnvelopeEstimator.hpp:42
testDataSize_
const unsigned int testDataSize_
Definition: PeakExtractorTest.cpp:167
pwiz::frequency::PeakDetectorMatchedFilter::Config::isotopeMaxNeutronCount
int isotopeMaxNeutronCount
isotope filter maximum number of neutrons to score
Definition: PeakDetectorMatchedFilter.hpp:73
pwiz::frequency::PeakDetectorMatchedFilter::Config::peakMaxCorrelationAngle
double peakMaxCorrelationAngle
maximum correlation angle (degrees) for initial peak reporting
Definition: PeakDetectorMatchedFilter.hpp:61
pwiz::frequency::PeakDetectorMatchedFilter::Config::peakThresholdFactor
double peakThresholdFactor
noise floor multiple for initial peak reporting threshold
Definition: PeakDetectorMatchedFilter.hpp:58
pwiz::frequency::PeakDetectorMatchedFilter::Config::filterMatchRate
int filterMatchRate
number of filter correlations computed per frequency step
Definition: PeakDetectorMatchedFilter.hpp:52
pwiz::util
Definition: almost_equal.hpp:33
pwiz::data::peakdata::PeakFamily::peaks
std::vector< Peak > peaks
Definition: PeakData.hpp:116
TEST_EPILOG
#define TEST_EPILOG
Definition: unit.hpp:183
pwiz::data::peakdata::Peak::intensity
double intensity
Definition: PeakData.hpp:59
pwiz::chemistry::IsotopeCalculator
Definition: IsotopeCalculator.hpp:37
pwiz::frequency::PeakDetectorMatchedFilter::Config::filterSampleRadius
int filterSampleRadius
number of filter samples taken on either side of 0
Definition: PeakDetectorMatchedFilter.hpp:55
Std.hpp
pwiz::chemistry::IsotopeEnvelopeEstimator::Config::isotopeCalculator
const IsotopeCalculator * isotopeCalculator
Definition: IsotopeEnvelopeEstimator.hpp:54
pwiz::frequency::PeakDetectorMatchedFilter::Score::charge
int charge
Definition: PeakDetectorMatchedFilter.hpp:132
testDataCalibrationA_
const double testDataCalibrationA_
Definition: PeakDetectorMatchedFilterTestData.hpp:43
pwiz::chemistry
Definition: Chemistry.hpp:37
os_
ostream * os_
Definition: PeakDetectorMatchedFilterTest.cpp:38
TEST_FAILED
#define TEST_FAILED(x)
Definition: unit.hpp:177
pwiz::chemistry::IsotopeEnvelopeEstimator::Config
Definition: IsotopeEnvelopeEstimator.hpp:46
testFind
void testFind(FrequencyData &fd, const IsotopeEnvelopeEstimator &isotopeEnvelopeEstimator)
Definition: PeakDetectorMatchedFilterTest.cpp:83
pwiz::frequency::PeakDetectorMatchedFilter::config
virtual const Config & config() const =0
access to the configuration
TEST_PROLOG
#define TEST_PROLOG(argc, argv)
Definition: unit.hpp:175
pwiz::data::peakdata::PeakFamily::charge
int charge
Definition: PeakData.hpp:114
testData_
TestDatum testData_[]
Definition: PeakExtractorTest.cpp:39
pwiz::data::peakdata::Peak::getAttribute
double getAttribute(Attribute attribute) const
initializeWithTestData
void initializeWithTestData(FrequencyData &fd)
Definition: PeakDetectorMatchedFilterTest.cpp:41
pwiz::frequency::PeakDetectorMatchedFilter::Score::monoisotopicFrequency
double monoisotopicFrequency
Definition: PeakDetectorMatchedFilter.hpp:136
test
void test()
Definition: PeakDetectorMatchedFilterTest.cpp:155
pwiz::frequency::PeakDetectorMatchedFilter::Score
structure for holding the matched filter calculation results
Definition: PeakDetectorMatchedFilter.hpp:129
pwiz::data::FrequencyDatum
SampleDatum< double, std::complex< double > > FrequencyDatum
Definition: FrequencyData.hpp:40
TestDatum
Definition: PeakDetectorMatchedFilterTestData.hpp:28
pwiz::data::peakdata::PeakData::scans
std::vector< Scan > scans
Definition: PeakData.hpp:194
pwiz::frequency::PeakDetectorMatchedFilter::Config::isotopeEnvelopeEstimator
const chemistry::IsotopeEnvelopeEstimator * isotopeEnvelopeEstimator
IsotopeEnvelopeEstimator pointer, must be valid for PeakDetector lifetime.
Definition: PeakDetectorMatchedFilter.hpp:49
main
int main(int argc, char *argv[])
Definition: PeakDetectorMatchedFilterTest.cpp:170
pwiz::data::FrequencyData::calibrationParameters
const CalibrationParameters & calibrationParameters() const
pwiz::data::FrequencyData::analyze
void analyze()
recache statistics calculations after any direct data changes via non-const data()
unit.hpp
unit_assert
#define unit_assert(x)
Definition: unit.hpp:85
pwiz::frequency::PeakDetectorMatchedFilter::Config::useMagnitudeFilter
bool useMagnitudeFilter
use the magnitude of the peak shape filter kernel for finding peaks
Definition: PeakDetectorMatchedFilter.hpp:79
pwiz::frequency::PeakDetectorMatchedFilter::Config::monoisotopicPeakThresholdFactor
double monoisotopicPeakThresholdFactor
noise floor multiple for monoisotopic peak threshold
Definition: PeakDetectorMatchedFilter.hpp:67
pwiz::frequency::PeakDetectorMatchedFilter::Config::isotopeThresholdFactor
double isotopeThresholdFactor
noise floor multiple for isotope filter threshold
Definition: PeakDetectorMatchedFilter.hpp:64
pwiz::math::MatchedFilter::details::norm
double norm(double d)
Definition: MatchedFilter.hpp:221
pwiz::frequency::PeakDetectorMatchedFilter::Config
structure for holding configuration
Definition: PeakDetectorMatchedFilter.hpp:46