ProteoWizard
MatchedFilterTest.cpp
Go to the documentation of this file.
1//
2// $Id$
3//
4//
5// Original author: Darren Kessner <darren@proteowizard.org>
6//
7// Copyright 2007 Spielberg Family Center for Applied Proteomics
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
24#include "MatchedFilter.hpp"
27#include <complex>
28#include <cstring>
29#include <typeinfo>
30
31using namespace pwiz::math;
32using namespace pwiz::util;
33
34
35ostream* os_ = 0;
36
37
39{
40 if (os_) *os_ << "test_SampledData()\n";
41 using namespace MatchedFilter;
42
43 SampledData<DxD> sd;
44 sd.domain = make_pair(0.0, 2.0);
45 sd.samples.push_back(10);
46 sd.samples.push_back(11);
47 sd.samples.push_back(12);
48 sd.samples.push_back(13);
49 sd.samples.push_back(14);
50
51 if (os_) *os_ << sd << endl;
52
53 if (os_) *os_ << "domainWidth: " << sd.domainWidth() << endl;
54 unit_assert(sd.domainWidth() == 2.0);
55
56 if (os_) *os_ << "dx: " << sd.dx() << endl;
57 unit_assert(sd.dx() == .5);
58
59 for (unsigned int i=0; i<sd.samples.size(); i++)
60 {
61 if (os_) *os_ << "x(" << i << "): " << sd.x(i) << endl;
62 unit_assert(sd.x(i) == i * .5);
63 }
64
65 unsigned int count = 0;
66 for (double x=-.2; x<2.3; x+=.1, count++)
67 {
68 //if (os_) *os_ << "sampleIndex(" << x << "): " << sd.sampleIndex(x) << endl;
69 unit_assert(sd.sampleIndex(x) == count/5);
70 //if (os_) *os_ << "sample(" << x << "): " << sd.sample(x) << endl;
71 unit_assert(sd.sample(x) == 10 + count/5);
72 }
73}
74
75
76// Kernel types defined as objects must define space_type in the class definition,
77// or in a specialization of KernelTraitsBase<>.
79{
81 double operator()(double d) const {return (d>=-1 && d<=1) ? 1 - abs(d) : 0;}
82};
83
84
85// KernelTraitsBase<> has default specialization for function pointer types
86complex<double> OneMinusAbsComplex(double d)
87{
88 return (d>=-1 && d<=1) ? 1 - abs(d) : 0;
89}
90
91
92template <typename Kernel>
93void test_createFilter(const Kernel& f)
94{
95 using namespace MatchedFilter;
96
97 if (os_) *os_ << "test_createFilter() " << typeid(f).name() << endl;
98
99 int sampleRadius = 4;
100 double dx = .25;
101 double shift = 0;
102
103 typedef typename KernelTraits<Kernel>::filter_type filter_type;
104 typedef typename KernelTraits<Kernel>::abscissa_type abscissa_type;
105 typedef typename KernelTraits<Kernel>::ordinate_type ordinate_type;
106
107 filter_type filter = details::createFilter(f, sampleRadius, dx, shift);
108
109 if (os_)
110 {
111 copy(filter.begin(), filter.end(), ostream_iterator<ordinate_type>(*os_, " "));
112 *os_ << endl;
113 }
114
115 unit_assert((int)filter.size() == sampleRadius*2 + 1);
116 for (int i=-sampleRadius; i<=sampleRadius; ++i)
117 unit_assert(filter[sampleRadius+i] == f(i*dx - shift));
118
119 if (os_) *os_ << endl;
120}
121
122
123template <typename Kernel>
124void test_createFilters(const Kernel& f)
125{
126 using namespace MatchedFilter;
127
128 if (os_) *os_ << "test_createFilters() " << typeid(f).name() << endl;
129
130 int sampleRadius = 2;
131 int subsampleFactor = 4;
132 double dx = 1;
133
134 typedef typename KernelTraits<Kernel>::filter_type filter_type;
135 typedef typename KernelTraits<Kernel>::ordinate_type ordinate_type;
136 vector<filter_type> filters = details::createFilters(f,
137 sampleRadius,
138 subsampleFactor,
139 dx);
140
141 // verify filter count
142 unit_assert((int)filters.size() == subsampleFactor);
143
144 for (typename vector<filter_type>::const_iterator it=filters.begin(); it!=filters.end(); ++it)
145 {
146 if (os_)
147 {
148 copy(it->begin(), it->end(), ostream_iterator<ordinate_type>(*os_, " "));
149 *os_ << endl;
150 }
151
152 // verify filter size
153 unit_assert((int)it->size() == sampleRadius*2 + 1);
154
155 // verify filter normalization
156 double sum = 0;
157 for (typename filter_type::const_iterator jt=it->begin(); jt!=it->end(); ++jt)
158 sum += norm(complex<double>(*jt));
159 unit_assert_equal(sum, 1, 1e-14);
160 }
161
162 if (os_) *os_ << endl;
163}
164
165
166template <typename Kernel>
167void test_compute(const Kernel& f)
168{
169 using namespace MatchedFilter;
170
171 if (os_) *os_ << "test_compute() " << typeid(f).name() << endl;
172
173 typename KernelTraits<Kernel>::sampled_data_type data;
174 data.domain = make_pair(0, 10);
175 data.samples.resize(11);
176 data.samples[5] = 1.;
177
178 if (os_) *os_ << "data: " << data << endl;
179
180 int sampleRadius = 2;
181 int sampleFactor = 4;
182
183 typedef typename KernelTraits<Kernel>::correlation_data_type CorrelationData;
184
185 CorrelationData correlationData =
186 computeCorrelationData(data, f, sampleRadius, sampleFactor);
187
188 if (os_) *os_ << "correlationData: " << correlationData << endl;
189
190 unit_assert(correlationData.samples.size() == 41);
191 unit_assert(abs(correlationData.samples[20].dot - 1.) < 1e-12);
192}
193
194
195template <typename Kernel>
196void test_kernel(const Kernel& kernel)
197{
198 if (os_) *os_ << "***************************************************************\n";
199 if (os_) *os_ << "test_kernel() " << typeid(kernel).name() << endl;
200 if (os_) *os_ << "***************************************************************\n";
201 if (os_) *os_ << endl;
202
203 test_createFilter(kernel);
204 test_createFilters(kernel);
205 test_compute(kernel);
206}
207
208
209int main(int argc, char* argv[])
210{
211 TEST_PROLOG(argc, argv)
212
213 try
214 {
215 if (argc>1 && !strcmp(argv[1],"-v")) os_ = &cout;
216 if (os_) *os_ << "MatchedFilterTest\n";
217
221
222 }
223 catch (exception& e)
224 {
225 TEST_FAILED(e.what())
226 }
227 catch (...)
228 {
229 TEST_FAILED("Caught unknown exception.")
230 }
231
233}
234
KernelTraitsBase< Kernel >::space_type::abscissa_type x
int main(int argc, char *argv[])
void test_SampledData()
complex< double > OneMinusAbsComplex(double d)
void test_kernel(const Kernel &kernel)
ostream * os_
void test_createFilter(const Kernel &f)
void test_createFilters(const Kernel &f)
void test_compute(const Kernel &f)
ProductSpace< double, double > DxD
MatchedFilter::DxD space_type
double operator()(double d) 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