ProteoWizard
MagnitudeLorentzianTest.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#include <boost/filesystem/operations.hpp>
30#include <cstring>
31
32
33using namespace pwiz::util;
34using namespace pwiz::frequency;
35using namespace pwiz::data;
36
37
38ostream* os_ = 0;
39double epsilon_ = numeric_limits<double>::epsilon();
40
41
43{
44 MagnitudeLorentzian m(1,0,1); // m(x) = 1/sqrt(x^2+1)
46 unit_assert_equal(m(1), 1/sqrt(2.), epsilon_);
47 unit_assert_equal(m(2), 1/sqrt(5.), epsilon_);
48 unit_assert_equal(m(3), 1/sqrt(10.), epsilon_);
49 unit_assert_equal(m(4), 1/sqrt(17.), epsilon_);
50
51 // center == 0, alpha == 2*pi, tau == 1/(2*pi)
53 unit_assert_equal(m.alpha(), 2*M_PI, epsilon_);
54 unit_assert_equal(m.tau(), 1/(2*M_PI), epsilon_);
55
56 if (os_) *os_ << "testBasic(): success!\n";
57}
58
59
60void testFit()
61{
62 MagnitudeLorentzian ref(1,0,1);
63
64 // choose sample values near 1!
65 // weighting pow(y,6) gives big roundoff errors
66
67 vector< pair<double,double> > samples;
68 for (int i=-2; i<3; i++)
69 samples.push_back(make_pair(i/10.,ref(i/10.)));
70
71 MagnitudeLorentzian m(samples);
72
73 if (os_)
74 {
75 *os_ << "coefficients: " << setprecision(14);
76 copy(m.coefficients().begin(), m.coefficients().end(), ostream_iterator<double>(*os_, " "));
77 *os_ << endl;
78
79 *os_ << "error: " << m(0)-1 << endl;
80
81 for (int i=0; i<5; i++)
82 *os_ << i << ", " << m(i) << endl;
83 }
84
85 unit_assert_equal(m(0), 1, epsilon_*100);
86 unit_assert_equal(m(1), 1/sqrt(2.), epsilon_*100);
87 unit_assert_equal(m(2), 1/sqrt(5.), epsilon_*100);
88 unit_assert_equal(m(3), 1/sqrt(10.), epsilon_*100);
89 unit_assert_equal(m(4), 1/sqrt(17.), epsilon_*100);
90 if (os_) *os_ << "testFit(): success!\n";
91}
92
93
95{
96 string filename = "MagnitudeLorentizianTest.cfd.temp.txt";
97 ofstream temp(filename.c_str());
98 temp << sampleData_;
99 temp.close();
100
101 FrequencyData fd(filename);
102 boost::filesystem::remove(filename);
103
105 if (os_) *os_ << "max: (" << max->x << ", " << abs(max->y) << ")\n";
106
107 // fit MagnitudeLorentzian to 3 points on unnormalized data
108
109 vector< pair<double,double> > samples1;
110 transform(fd.max()-1, fd.max()+2, back_inserter(samples1), FrequencyData::magnitudeSample);
111
112 if (os_)
113 {
114 *os_ << "raw data:\n";
115 for (unsigned int i=0; i<samples1.size(); i++)
116 *os_ << "sample " << i << ": (" << samples1[i].first << ", " << samples1[i].second << ")\n";
117 }
118
119 const MagnitudeLorentzian m1(samples1);
120
121 if (os_)
122 {
123 *os_ << "m1: ";
124 copy(m1.coefficients().begin(), m1.coefficients().end(), ostream_iterator<double>(*os_, " "));
125 *os_ << endl;
126 *os_ << "error: " << scientific << m1.leastSquaresError() << endl;
127
128 for (unsigned int i=0; i<samples1.size(); i++)
129 *os_ << "m1(" << i << ") == " << m1(samples1[i].first) << endl;
130 }
131
132 // now on normalized data
133
134
135 fd.normalize();
136
137 vector< pair<double,double> > samples2;
138 transform(fd.max()-1, fd.max()+2, back_inserter(samples2), FrequencyData::magnitudeSample);
139
140 if (os_)
141 {
142 *os_ << "normalized: \n";
143 for (unsigned int i=0; i<samples2.size(); i++)
144 *os_ << "sample " << i << ": (" << samples2[i].first << ", " << samples2[i].second << ")\n";
145 }
146
147 const MagnitudeLorentzian m2(samples2);
148
149 if (os_)
150 {
151 *os_ << "m2: ";
152 copy(m2.coefficients().begin(), m2.coefficients().end(), ostream_iterator<double>(*os_, " "));
153 *os_ << endl;
154 *os_ << "error: " << scientific << m2.leastSquaresError() << endl;
155
156 for (unsigned int i=0; i<samples2.size(); i++)
157 *os_ << "m2(" << i << ") == " << m2(samples2[i].first) << " [" << fd.scale()*m2(samples2[i].first) << "]\n";
158 }
159
160 unit_assert_equal(m2.leastSquaresError(), 0, 1e-15);
161}
162
163
164int main(int argc, char* argv[])
165{
166 TEST_PROLOG(argc, argv)
167
168 try
169 {
170 if (argc>1 && !strcmp(argv[1],"-v")) os_ = &cout;
171 if (os_) *os_ << "MagnitudeLorentzianTest\n";
172 testBasic();
173 testFit();
174 testData();
175 }
176 catch (exception& e)
177 {
178 TEST_FAILED(e.what())
179 }
180 catch (...)
181 {
182 TEST_FAILED("Caught unknown exception.")
183 }
184
186}
187
int main(int argc, char *argv[])
double epsilon_
void testData()
void testBasic()
void testFit()
ostream * os_
const char * sampleData_
Class for binary storage of complex frequency data.
const_iterator max() const
returns an iterator to FrequencyDatum with highest magnitude
container::const_iterator const_iterator
std::complex< double > scale() const
return current scale of data (compared to original)
static std::pair< double, double > magnitudeSample(const FrequencyDatum &datum)
Returns a <frequency,magnitude> pair.
void normalize()
normalize by transform( -max.x, 1/abs(max.y) )
std::vector< double > & coefficients()
#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