ProteoWizard
DerivativeTest.hpp
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
24#ifndef _DERIVATIVETEST_HPP_
25#define _DERIVATIVETEST_HPP_
26
27
31#include <iomanip>
32#include <exception>
33#include <stdexcept>
34
35
36namespace pwiz {
37namespace frequency {
38namespace DerivativeTest {
39
40
41using namespace pwiz::util;
42
43
44// VectorFunction is a generic interface used for numeric derivative calculations.
45// Clients implement the pure virtual functions, and the calculations are handled automatically
46
47
48template<typename value_type>
50{
51 public:
52
53 virtual unsigned int argumentCount() const = 0;
54 virtual unsigned int valueCount() const = 0;
55 virtual ublas::vector<value_type> operator()(ublas::vector<double> x) const = 0;
56
57 ublas::matrix<value_type> differenceQuotient(ublas::vector<double> x, double delta) const;
58 void printDifferenceQuotientSequence(ublas::vector<double> x,
59 std::ostream& os) const;
60
61 protected:
62 virtual ~VectorFunction(){}
63};
64
65
66template<typename value_type>
67ublas::matrix<value_type> VectorFunction<value_type>::differenceQuotient(ublas::vector<double> x, double delta) const
68{
69 ublas::matrix<value_type> result(argumentCount(), valueCount());
70 result.clear();
71 for (unsigned int i=0; i<argumentCount(); i++)
72 {
73 ublas::vector<double> x2(x);
74 x2(i) += delta;
75 row(result, i) = ((*this)(x2)-(*this)(x))/delta;
76 }
77 return result;
78}
79
80
81template<typename value_type>
83 std::ostream& os) const
84{
85 using namespace std;
86
87 for (double delta=.1; delta>1e-9; delta/=10)
88 {
89 os << scientific << setprecision(1) << "[delta: " << delta << "] ";
90 os.unsetf(std::ios::scientific);
91 os << setprecision(8) << differenceQuotient(x, delta) << endl;
92 }
93}
94
95
96// Adaptors from ParametrizedFunction and its derivative, error, and error-derivative functions to
97// implement the VectorFunction interface.
98
99
100template<typename value_type>
102{
103 public:
107
108 virtual unsigned int argumentCount() const {return f_.parameterCount();}
109 virtual unsigned int valueCount() const {return 1;}
110
111 virtual ublas::vector<value_type> operator()(ublas::vector<double> p) const
112 {
113 ublas::vector<value_type> result(1);
114 result(0) = f_(x_, p);
115 return result;
116 }
117
118 private:
120 double x_;
121};
122
123
124template<typename value_type>
126{
127 public:
131
132 virtual unsigned int argumentCount() const {return f_.parameterCount();}
133 virtual unsigned int valueCount() const {return f_.parameterCount();}
134
135 virtual ublas::vector<value_type> operator()(ublas::vector<double> p) const
136 {
137 return f_.dp(x_,p);
138 }
139
140 private:
142 double x_;
143};
144
145
146template<typename value_type>
147class AdaptedErrorFunction : public VectorFunction<double> // double: error functions are real-valued
148{
149 public:
153
154 virtual unsigned int argumentCount() const {return e_.parameterCount();}
155 virtual unsigned int valueCount() const {return 1;}
156
157 virtual ublas::vector<double> operator()(ublas::vector<double> p) const
158 {
159 ublas::vector<double> result(1);
160 result(0) = e_(p);
161 return result;
162 }
163
164 private:
166};
167
168
169template<typename value_type>
171{
172 public:
176
177 virtual unsigned int argumentCount() const {return e_.parameterCount();}
178 virtual unsigned int valueCount() const {return e_.parameterCount();}
179
180 virtual ublas::vector<double> operator()(ublas::vector<double> p) const
181 {
182 return e_.dp(p);
183 }
184
185 private:
187};
188
189
190// Numeric derivative calculations for ParametrizedFunction and associated ErrorFunction.
191
192
193template<typename value_type>
195 double x,
196 const ublas::vector<double>& p,
197 std::ostream* os = 0,
198 double delta = 1e-7,
199 double epsilon = 1e-4)
200{
201 using namespace std;
202
203 if (os)
204 {
205 *os << "x: " << x << endl;
206 *os << "p: " << p << endl;
207 }
208
209 if (os) *os << "f.dp: " << f.dp(x,p) << endl;
211 if (os) slice.printDifferenceQuotientSequence(p, *os);
212
213 ublas::matrix<value_type> dp(f.dp(x,p).size(),1);
214 column(dp,0) = f.dp(x,p);
216
217 if (os) *os << "f.dp2: " << f.dp2(x,p) << endl;
219 if (os) derivativeSlice.printDifferenceQuotientSequence(p, *os);
220
221 unit_assert_matrices_equal(f.dp2(x,p), derivativeSlice.differenceQuotient(p,delta), epsilon);
222}
223
224
225template<typename value_type>
227 const ublas::vector<double>& p,
228 std::ostream* os = 0,
229 double delta = 1e-7,
230 double epsilon = 1e-4)
231{
232 using namespace std;
233
234 if (os) *os << "p: " << p << endl;
235
236 if (os) *os << "e.dp: " << e.dp(p) << endl;
238 if (os) adapted.printDifferenceQuotientSequence(p, *os);
239
240 ublas::matrix<value_type> dp(e.dp(p).size(), 1);
241 column(dp,0) = e.dp(p);
243
244 if (os) *os << "e.dp2: " << e.dp2(p) << endl;
245 AdaptedErrorDerivative<value_type> adaptedDerivative(e);
246 if (os) adaptedDerivative.printDifferenceQuotientSequence(p, *os);
247
248 unit_assert_matrices_equal(e.dp2(p), adaptedDerivative.differenceQuotient(p,delta), epsilon);
249}
250
251
252} // namespace DerivativeTest
253} // namespace frequency
254} // namespace pwiz
255
256
257#endif // _DERIVATIVETEST_HPP_
258
KernelTraitsBase< Kernel >::space_type::abscissa_type x
void testDerivatives()
const ParametrizedFunction< value_type >::ErrorFunction & e_
AdaptedErrorDerivative(const typename ParametrizedFunction< value_type >::ErrorFunction &e)
virtual ublas::vector< double > operator()(ublas::vector< double > p) const
virtual ublas::vector< double > operator()(ublas::vector< double > p) const
AdaptedErrorFunction(const typename ParametrizedFunction< value_type >::ErrorFunction &e)
const ParametrizedFunction< value_type >::ErrorFunction & e_
ParametrizedDerivativeSlice(const ParametrizedFunction< value_type > &f, double x)
virtual ublas::vector< value_type > operator()(ublas::vector< double > p) const
ParametrizedFunctionSlice(const ParametrizedFunction< value_type > &f, double x)
virtual ublas::vector< value_type > operator()(ublas::vector< double > p) const
const ParametrizedFunction< value_type > & f_
virtual ublas::vector< value_type > operator()(ublas::vector< double > x) const =0
virtual unsigned int valueCount() const =0
virtual unsigned int argumentCount() const =0
void printDifferenceQuotientSequence(ublas::vector< double > x, std::ostream &os) const
ublas::matrix< value_type > differenceQuotient(ublas::vector< double > x, double delta) const
ublas::vector< double > dp(const ublas::vector< double > &p) const
ublas::matrix< double > dp2(const ublas::vector< double > &p) const
virtual ublas::matrix< value_type > dp2(double x, const ublas::vector< double > &p) const =0
virtual ublas::vector< value_type > dp(double x, const ublas::vector< double > &p) const =0
const double epsilon
Definition DiffTest.cpp:41
STL namespace.
#define unit_assert_matrices_equal(A, B, epsilon)
Definition unit.hpp:135