ProteoWizard
ParametrizedFunctionTest.cpp
Go to the documentation of this file.
1 //
2 // $Id: ParametrizedFunctionTest.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 
24 #include "DerivativeTest.hpp"
25 #include "TruncatedLorentzian.hpp"
28 #include <cstring>
29 
30 
31 using namespace pwiz::util;
32 using namespace pwiz::frequency;
33 
34 
35 ostream* os_ = 0;
36 
37 
39 {
40  if (os_) *os_ << "****************************************************\n";
41  if (os_) *os_ << "testDifferenceQuotient()\n";
42 
43  using namespace DerivativeTest;
44 
45  class TestFunction : public VectorFunction<double>
46  {
47  public:
48 
49  // f(x,y) = (x^2, xy, y^2)
50 
51  virtual unsigned int argumentCount() const {return 2;}
52  virtual unsigned int valueCount() const {return 3;}
53 
54  virtual ublas::vector<double> operator()(ublas::vector<double> x) const
55  {
56  if (x.size() != argumentCount())
57  throw logic_error("[TestFunction::()] Wrong argument count.");
58 
59  ublas::vector<double> result(3);
60  result(0) = x(0)*x(0);
61  result(1) = x(0)*x(1);
62  result(2) = x(1)*x(1);
63  return result;
64  }
65  };
66 
67 
68  TestFunction f;
69  ublas::vector<double> args(2);
70  args(0) = 5; args(1) = 7;
71  if (os_) *os_ << "f(5,7): " << f(args) << endl;
72 
73  if (os_) f.printDifferenceQuotientSequence(args, *os_);
74 
75  // f'(x,y) = ((2x, y, 0), (0, x, 2y))
76  // f'(5,7) = ((10, 7, 0), (0, 5, 14))
77 
78  ublas::matrix<double> d(2,3);
79  d(0,0) = 10;
80  d(0,1) = 7;
81  d(0,2) = 0;
82  d(1,0) = 0;
83  d(1,1) = 5;
84  d(1,2) = 14;
85 
86  const double delta = 1e-9;
87  const double epsilon = 1e-5;
88  unit_assert_matrices_equal(d, f.differenceQuotient(args,delta), epsilon);
89 }
90 
91 
93 {
94  // F(x) = Acos(Bx), p = <A, B>
95 
96  public:
97 
98  virtual unsigned int parameterCount() const {return 2;}
99 
100  virtual double operator()(double x, const ublas::vector<double>& p) const
101  {
102  preprocess(x,p);
103  return A_*cosBx_;
104  }
105 
106 
107  virtual ublas::vector<double> dp(double x, const ublas::vector<double>& p) const
108  {
109  preprocess(x,p);
110  ublas::vector<double> v(2);
111  v(0) = cosBx_; // dF/dA
112  v(1) = -A_*x_*sinBx_; // dF/dB
113  return v;
114  }
115 
116 
117  virtual ublas::matrix<double> dp2(double x, const ublas::vector<double>& p) const
118  {
119  preprocess(x,p);
120  ublas::matrix<double> m(2,2);
121  m(0,0) = 0; // d2F/dA2
122  m(1,0) = m(0,1) = -x_*sinBx_; // d2F/dAdB
123  m(1,1) = -A_*x_*x_*cosBx_; // d2F/dB2
124  return m;
125  }
126 
127 
128  private:
129 
130  void preprocess(double x, const ublas::vector<double>& p) const
131  {
132  // check parameter size
133  if (p.size() != parameterCount())
134  throw logic_error("[Parabola] Wrong parameter size.");
135 
136  // cache arguments and do expensive calculations
137  if (x!=x_ || p(0)!=A_ || p(1)!=B_)
138  {
139  x_ = x;
140  A_ = p(0);
141  B_ = p(1);
142  sinBx_ = sin(B_*x);
143  cosBx_ = cos(B_*x);
144  }
145  else
146  {
147  //if (os_) *os_ << "cache hit!\n";
148  }
149  }
150 
151  // cached values
152  mutable double x_;
153  mutable double A_;
154  mutable double B_;
155  mutable double sinBx_;
156  mutable double cosBx_;
157 };
158 
159 
161 {
162  if (os_) *os_ << "****************************************************\n";
163  if (os_) *os_ << "testDerivatives()\n";
164 
166 
167  ublas::vector<double> p(2);
168  p(0) = 5;
169  p(1) = M_PI/4;
170 
171  for (int i=0; i<8; i++)
173 }
174 
175 
177 {
178  if (os_) *os_ << "****************************************************\n";
179  if (os_) *os_ << "testErrorFunction()\n";
180 
182 
183  ublas::vector<double> p(2);
184  p(0) = 4;
185  p(1) = 30;
186 
187  ParametrizedCosine::ErrorFunction::Data data;
188  typedef ParametrizedCosine::ErrorFunction::Datum Datum;
189  data.push_back(Datum(0,3));
190  data.push_back(Datum(M_PI/2,0));
191 
192  ParametrizedCosine::ErrorFunction e(f, data);
193  if (os_) *os_ << "error: " << e(p) << endl;
194 
195  DerivativeTest::testDerivatives<double>(e, p, os_);
196 
197  if (os_) *os_ << "8*pi^2: " << 8*M_PI*M_PI << endl;
198 }
199 
200 
202 {
203  if (os_) *os_ << "****************************************************\n";
204  if (os_) *os_ << "testErrorLorentzian()\n";
205 
206  TruncatedLorentzian f(1);
207 
208  ublas::vector<double> p(4);
212  p(TruncatedLorentzian::F0) = 0;
213 
214  TruncatedLorentzian::ErrorFunction::Data data;
215  typedef TruncatedLorentzian::ErrorFunction::Datum Datum;
216  data.push_back(Datum(0,3));
217  data.push_back(Datum(M_PI/2,0));
218 
219  TruncatedLorentzian::ErrorFunction e(f, data);
220  if (os_) *os_ << "error: " << e(p) << endl;
221 
222  DerivativeTest::testDerivatives< complex<double> >(e, p, os_);
223 }
224 
225 
226 int main(int argc, char* argv[])
227 {
228  TEST_PROLOG(argc, argv)
229 
230  try
231  {
232  if (argc>1 && !strcmp(argv[1],"-v")) os_ = &cout;
233  if (os_) *os_ << "ParametrizedFunctionTest\n";
235  testDerivatives();
238  }
239  catch (exception& e)
240  {
241  TEST_FAILED(e.what())
242  }
243  catch (...)
244  {
245  TEST_FAILED("Caught unknown exception.")
246  }
247 
249 }
250