ProteoWizard
MatchedFilter.hpp
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#ifndef _MATCHEDFILTER_HPP_
25#define _MATCHEDFILTER_HPP_
26
27
29#include <iostream>
30#include <vector>
31#include <iterator>
32#include <algorithm>
33#include <complex>
34#include <limits>
35
36
37namespace pwiz {
38namespace math {
39namespace MatchedFilter {
40
41
42template <typename X, typename Y>
44{
45 typedef X abscissa_type;
47};
48
49
52
53
54template <typename space_type>
56{
57 typedef typename space_type::abscissa_type abscissa_type;
58 typedef typename space_type::ordinate_type ordinate_type;
59 typedef std::pair<abscissa_type,abscissa_type> domain_type;
60 typedef std::vector<ordinate_type> samples_type;
61
64
66 {
67 return domain.second - domain.first;
68 }
69
71 {
72 return samples.empty() ? 0 : domainWidth()/(samples.size()-1);
73 }
74
75 abscissa_type x(typename samples_type::size_type index) const
76 {
77 return domain.first + domainWidth() * index / (samples.size()-1);
78 }
79
80 typename samples_type::size_type sampleIndex(const abscissa_type& x) const
81 {
82 typename samples_type::size_type sampleCount = samples.size();
83 int result = (int)round((sampleCount-1)*(x-domain.first)/domainWidth());
84 if (result < 0) result = 0;
85 if (result > (int)(sampleCount-1)) result = sampleCount-1;
86 return (typename samples_type::size_type)(result);
87 }
88
90 {
91 return samples[sampleIndex(x)];
92 }
93};
94
95
96template <typename space_type>
97std::ostream& operator<<(std::ostream& os, const SampledData<space_type>& data)
98{
99 os << "[" << data.domain.first << "," << data.domain.second << "] "
100 << "(" << data.samples.size() << " samples)\n";
101
103 for (unsigned int index=0; index!=data.samples.size(); ++index, ++it)
104 os << data.x(index) << "\t" << *it << std::endl;
105
106 return os;
107}
108
109
110template <typename Kernel>
112{
113 // When using a kernel function of type Kernel,
114 // KernelTraitsBase<Kernel> must define space_type,
115 // which in turn must define abscissa_type and ordinate_type, e.g:
116 // typedef ProductSpace<X,Y> space_type;
117
118 // As a shortcut, the following default typedef allows a client to
119 // define space_type in the definition of Kernel:
120 typedef typename Kernel::space_type space_type;
121};
122
123
124// partial specialization of KernelTraitsBase for function pointers
125template <typename X, typename Y>
126struct KernelTraitsBase<Y(*)(X)>
127{
129};
130
131
132namespace {
133
134template <typename Kernel>
135struct KernelConcept
136{
137 void check()
138 {
139 y = k(x);
140 }
141
142 typename KernelTraitsBase<Kernel>::space_type::abscissa_type x;
143 typename KernelTraitsBase<Kernel>::space_type::ordinate_type y;
144 Kernel k;
145};
146
147}
148
149template <typename Kernel>
151{
152 // force compile of KernelConcept::check()
153 void (KernelConcept<Kernel>::*dummy)() = &KernelConcept<Kernel>::check;
154 (void)dummy;
155}
156
157
158template <typename Y>
160{
162 double e2;
163 double tan2angle;
164
165 Correlation(Y _dot = 0, double _e2 = 0, double _tan2angle = 0)
166 : dot(_dot), e2(_e2), tan2angle(_tan2angle)
167 {}
168
169 double angle() const {return atan(sqrt(tan2angle))*180/M_PI;}
170};
171
172
173template<typename Y>
174std::ostream& operator<<(std::ostream& os, const Correlation<Y>& c)
175{
176 os << "<" << c.dot << ", " << c.e2 << ", " << c.angle() << ">";
177 return os;
178}
179
180
181template <typename Kernel>
183{
185 typedef typename space_type::abscissa_type abscissa_type;
186 typedef typename space_type::ordinate_type ordinate_type;
187
189 typedef typename sampled_data_type::samples_type samples_type;
194
195 // verify Kernel concept at compile time
196 template <void(*T)()> struct Dummy;
198};
199
200
201namespace details {
202
203
204template <typename Kernel>
206createFilter(const Kernel& kernel,
207 int sampleRadius,
210{
211 checkKernelConcept<Kernel>();
212
213 typename KernelTraits<Kernel>::filter_type filter;
214 for (int i=-sampleRadius; i<=sampleRadius; i++)
215 filter.push_back(kernel(i*dx - shift));
216 return filter;
217}
218
219
220// mimic complex<> functions
221inline double norm(double d) {return d*d;}
222inline double conj(double d) {return d;}
223
224
225template <typename Filter>
226void normalizeFilter(Filter& filter)
227{
228 double normalization = 0;
229 for (typename Filter::const_iterator it=filter.begin(); it!=filter.end(); ++it)
230 normalization += norm(*it);
231 normalization = sqrt(normalization);
232
233 for (typename Filter::iterator it=filter.begin(); it!=filter.end(); ++it)
234 *it /= normalization;
235}
236
237
238template <typename Kernel>
239std::vector<typename KernelTraits<Kernel>::filter_type>
240createFilters(const Kernel& kernel,
241 int sampleRadius,
242 int subsampleFactor,
244{
245 checkKernelConcept<Kernel>();
246
247 typedef typename KernelTraits<Kernel>::filter_type filter_type;
248 std::vector<filter_type> filters;
249
250 for (int i=0; i<subsampleFactor; i++)
251 filters.push_back(createFilter(kernel, sampleRadius, dx, dx*i/subsampleFactor));
252
253 for_each(filters.begin(), filters.end(), normalizeFilter<filter_type>);
254
255 return filters;
256}
257
258
259template<typename Kernel>
260void
265{
266 checkKernelConcept<Kernel>();
267 result.dot = 0;
268 double normData = 0;
269
270 for (; samples!=samplesEnd; ++samples, ++filter)
271 {
272 result.dot += (*samples) * conj(*filter);
273 normData += norm(*samples);
274 }
275
276 double normDot = norm(result.dot);
277 result.e2 = (std::max)(normData - normDot, 0.);
278 result.tan2angle = normDot>0 ? result.e2/normDot : std::numeric_limits<double>::infinity();
279}
280
281
282} // namespace details
283
284
285template<typename Kernel>
288 const Kernel& kernel,
289 int sampleRadius,
290 int subsampleFactor)
291{
292 checkKernelConcept<Kernel>();
293
294 typedef typename KernelTraits<Kernel>::correlation_data_type result_type;
295 result_type result;
296
297 result.domain = data.domain;
298 if (data.samples.empty()) return result;
299 result.samples.resize((data.samples.size()-1) * subsampleFactor + 1);
300
301 typedef typename KernelTraits<Kernel>::filter_type filter_type;
302 std::vector<filter_type> filters = details::createFilters(kernel,
303 sampleRadius,
304 subsampleFactor,
305 data.dx());
306
307 typedef typename KernelTraits<Kernel>::samples_type samples_type;
308
309 unsigned int sampleIndex = sampleRadius;
310 for (typename samples_type::const_iterator itData = data.samples.begin() + sampleRadius;
311 itData + sampleRadius != data.samples.end(); ++itData, ++sampleIndex)
312 for (unsigned int filterIndex=0; filterIndex<filters.size(); ++filterIndex)
313 {
314 unsigned int index = sampleIndex * filters.size() + filterIndex;
315
316 if (index >= result.samples.size()) // only when sampleRadius==0, filterIndex>0
317 break;
318
319 details::computeCorrelation<Kernel>(itData-sampleRadius, itData+sampleRadius+1,
320 filters[filterIndex].begin(),
321 result.samples[index]);
322 }
323
324 return result;
325}
326
327
328} // namespace MatchedFilter
329} // namespace math
330} // namespace pwiz
331
332
333#endif // _MATCHEDFILTER_HPP_
334
Y
Definition Chemistry.hpp:84
KernelTraitsBase< Kernel >::space_type::abscissa_type x
Kernel k
KernelTraitsBase< Kernel >::space_type::ordinate_type y
KernelTraits< Kernel >::filter_type createFilter(const Kernel &kernel, int sampleRadius, typename KernelTraits< Kernel >::abscissa_type dx, typename KernelTraits< Kernel >::abscissa_type shift)
void computeCorrelation(typename KernelTraits< Kernel >::samples_type::const_iterator samples, typename KernelTraits< Kernel >::samples_type::const_iterator samplesEnd, typename KernelTraits< Kernel >::samples_type::const_iterator filter, typename KernelTraits< Kernel >::correlation_type &result)
std::vector< typename KernelTraits< Kernel >::filter_type > createFilters(const Kernel &kernel, int sampleRadius, int subsampleFactor, typename KernelTraits< Kernel >::abscissa_type dx)
ProductSpace< double, double > DxD
std::ostream & operator<<(std::ostream &os, const SampledData< space_type > &data)
ProductSpace< double, std::complex< double > > DxCD
KernelTraits< Kernel >::correlation_data_type computeCorrelationData(const typename KernelTraits< Kernel >::sampled_data_type &data, const Kernel &kernel, int sampleRadius, int subsampleFactor)
Correlation(Y _dot=0, double _e2=0, double _tan2angle=0)
Correlation< ordinate_type > correlation_type
sampled_data_type::samples_type samples_type
Dummy< &checkKernelConcept< Kernel > > dummy
ProductSpace< abscissa_type, correlation_type > correlation_space_type
SampledData< space_type > sampled_data_type
space_type::ordinate_type ordinate_type
SampledData< correlation_space_type > correlation_data_type
space_type::abscissa_type abscissa_type
KernelTraitsBase< Kernel >::space_type space_type
space_type::abscissa_type abscissa_type
samples_type::size_type sampleIndex(const abscissa_type &x) const
std::pair< abscissa_type, abscissa_type > domain_type
const ordinate_type & sample(abscissa_type x) const
abscissa_type x(typename samples_type::size_type index) const
space_type::ordinate_type ordinate_type
std::vector< ordinate_type > samples_type