ProteoWizard
ChemistryTest.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 
25 #include "Chemistry.hpp"
26 #include "Ion.hpp"
28 #include <cstring>
30 #include "boost/thread/thread.hpp"
31 #include "boost/thread/barrier.hpp"
32 
33 
34 using namespace pwiz::util;
35 using namespace pwiz::chemistry;
36 
37 
38 ostream* os_ = 0;
39 
40 
42 {
43  MassAbundance ma(1, 2);
44  MassAbundance ma2(1, 4);
45  unit_assert(ma != ma2);
46  ma2.abundance = 2;
47  unit_assert(ma == ma2);
48 }
49 
50 
52 {
53  const char* formula;
54  int numC, numH, numN, numO, numS, num13C, num2H, num15N, num18O;
55  double monoMass;
56  double avgMass;
57 };
58 
60 {
61  { "C1H2N3O4S5", 1, 2, 3, 4, 5, 0, 0, 0, 0, 279.864884, 280.36928 },
62  { "CH2N3O4S5", 1, 2, 3, 4, 5, 0, 0, 0, 0, 279.864884, 280.36928 },
63  { "H2N3O4S5C", 1, 2, 3, 4, 5, 0, 0, 0, 0, 279.864884, 280.36928 },
64  { "C H2N3O4S5", 1, 2, 3, 4, 5, 0, 0, 0, 0, 279.864884, 280.36928 },
65  { "H2N3O4S5C ", 1, 2, 3, 4, 5, 0, 0, 0, 0, 279.864884, 280.36928 },
66  { "H2N3O4S5 C ", 1, 2, 3, 4, 5, 0, 0, 0, 0, 279.864884, 280.36928 },
67  { "C1H 2N3O4 S5", 1, 2, 3, 4, 5, 0, 0, 0, 0, 279.864884, 280.36928 },
68  { "H-42", 0, -42, 0, 0, 0, 0, 0, 0, 0, -42.328651, -42.333512 },
69  { "N2C-1", -1, 0, 2, 0, 0, 0, 0, 0, 0, 28.006148 - 12, 28.013486 - 12.0107 },
70  { "C39H67N11O10", 39, 67, 11, 10, 0, 0, 0, 0, 0, 849.507238, 850.01698 },
71  { "C3H7N1O2Se1", 3, 7, 1, 2, 0, 0, 0, 0, 0, 168.9642, 168.0532 },
72  { "_13C1_2H2_15N3_18O4", 0, 0, 0, 0, 0, 1, 2, 3, 4, 134.0285, 134.0285 },
73  { "_13C1D2_15N3_18O4", 0, 0, 0, 0, 0, 1, 2, 3, 4, 134.0285, 134.0285 },
74  { "_13C1_3H2_15N3_18O4", 0, 0, 0, 0, 0, 1, 2, 3, 4, 136.0324, 136.0324 },
75  { "_13C_3H2_15N3_18O4", 0, 0, 0, 0, 0, 1, 2, 3, 4, 136.0324, 136.0324 },
76  { "_3H2_15N3_18O4_13C", 0, 0, 0, 0, 0, 1, 2, 3, 4, 136.0324, 136.0324 },
77  { "_13C1T2_15N3_18O4", 0, 0, 0, 0, 0, 1, 2, 3, 4, 136.0324, 136.0324 },
78  { "C-1 _13C1 _2H2 H-2 N-3 _15N3 _18O4 O-4", -1, -2, -3, -4, 0, 1, 2, 3, 4, 14.024, 13.984 },
79  { "C-1_13C1_2H2H-2N-3_15N3_18O4O-4", -1, -2, -3, -4, 0, 1, 2, 3, 4, 14.024, 13.984 },
80 };
81 
82 const int testFormulaDataSize = sizeof(testFormulaData)/sizeof(TestFormula);
83 
85 {
86  for (int i=0; i < testFormulaDataSize; ++i)
87  {
88  const TestFormula& testFormula = testFormulaData[i];
89  Formula formula(testFormula.formula);
90 
91  const double EPSILON = 0.001;
92 
93  if (os_) *os_ << formula << " " << formula.monoisotopicMass() << " " << formula.molecularWeight() << endl;
94  unit_assert_equal(formula.monoisotopicMass(), testFormula.monoMass, EPSILON);
95  unit_assert_equal(formula.molecularWeight(), testFormula.avgMass, EPSILON);
96 
97  formula[Element::C] += 2;
98  if (os_) *os_ << formula << " " << formula.monoisotopicMass() << " " << formula.molecularWeight() << endl;
99  unit_assert_equal(formula.monoisotopicMass(), testFormula.monoMass+24, EPSILON);
100  unit_assert_equal(formula.molecularWeight(), testFormula.avgMass+12.0107*2, EPSILON);
101 
102  //const Formula& constFormula = formula;
103  //constFormula[Element::C] = 1; // this won't compile, by design
104 
105  // test copy constructor
106  Formula formula2 = formula;
107  formula2[Element::C] -= 2;
108  if (os_) *os_ << "formula: " << formula << endl;
109  if (os_) *os_ << "formula2: " << formula2 << endl;
110  unit_assert_equal(formula.monoisotopicMass(), testFormula.monoMass+24, EPSILON);
111  unit_assert_equal(formula2.monoisotopicMass(), testFormula.monoMass, EPSILON);
112 
113  // test operator=
114  formula = formula2;
115  if (os_) *os_ << "formula: " << formula << endl;
116  if (os_) *os_ << "formula2: " << formula2 << endl;
117  unit_assert_equal(formula.monoisotopicMass(), formula2.monoisotopicMass(), EPSILON);
118 
119  // test operator==
120  unit_assert(formula == testFormula.formula); // implicit construction from string
121  unit_assert(formula == formula2);
122  formula2[Element::C] += 4; // test difference in CHONSP
123  unit_assert(formula != formula2);
124  formula2[Element::C] -= 4;
125  unit_assert(formula == formula2);
126  formula2[Element::U] += 2; // test difference outside CHONSP
127  unit_assert(formula != formula2);
128  formula2[Element::U] -= 2;
129  unit_assert(formula == formula2);
130 
131  // test data()
132  Formula::Map data = formula.data();
133  if (os_)
134  {
135  *os_ << "map: ";
136  for (Formula::Map::iterator it=data.begin(), end=data.end(); it!=end; ++it)
137  *os_ << it->first << it->second << " ";
138  *os_ << "\n";
139  }
140  }
141 }
142 
143 
145 {
146  using namespace Element;
147 
148  Formula water("H2O1");
149  Formula a("C1 H2 N3 O4 S5");
150 
151  a *= 2;
152  unit_assert(a[C]==2 && a[H]==4 && a[N]==6 && a[O]==8 && a[S]==10);
153  a += water;
154  unit_assert(a[H]==6 && a[O]==9);
155  a -= water;
156  unit_assert(a[C]==2 && a[H]==4 && a[N]==6 && a[O]==8 && a[S]==10);
157  a += 2*water;
158  unit_assert(a[H]==8 && a[O]==10);
159  a = (a - water*2);
160  unit_assert(a[C]==2 && a[H]==4 && a[N]==6 && a[O]==8 && a[S]==10);
161  a = water + water;
162  unit_assert(a[H]==4 && a[O]==2);
163  if (os_) *os_ << "water: " << a-water << endl;
164 }
165 
166 
167 void testInfo()
168 {
169  if (os_)
170  {
171  for (Element::Type e=Element::H; e<=Element::Ca; e=(Element::Type)(e+1))
172  *os_ << Element::Info::record(e).symbol << " " << Element::Info::record(e).atomicNumber << endl;
173  *os_ << endl;
174  }
175 
176  unit_assert(Element::Info::record(Element::C).atomicNumber == 6);
177  unit_assert(Element::Info::record("C").atomicNumber == 6);
178  unit_assert(Element::Info::record(Element::U).atomicNumber == 92);
179  unit_assert(Element::Info::record("U").atomicNumber == 92);
180  unit_assert(Element::Info::record(Element::Uuh).atomicNumber == 116);
181  unit_assert(Element::Info::record("Uuh").atomicNumber == 116);
182 
184  runtime_error,
185  "[chemistry::text2enum()] Error translating symbol foo");
186 }
187 
188 
190 {
191  if (os_)
192  {
193  *os_ << "Sulfur isotopes: " << Element::Info::record(Element::S).isotopes.size() << endl
194  << Element::Info::record(Element::S).isotopes;
195  }
196 }
197 
198 
200 {
201  Formula formula("Si6C12H36O6");
202 
203  if (os_)
204  {
205  *os_ << "polysiloxane:\n"
206  << formula << " "
207  << formula.monoisotopicMass() << " "
208  << formula.molecularWeight() << endl
209  << "ion: " << Ion::mz(formula.monoisotopicMass(), 1) << endl;
210  }
211 }
212 
213 
214 void testThreadSafetyWorker(boost::barrier* testBarrier, int& result)
215 {
216  testBarrier->wait(); // wait until all threads have started
217 
218  try
219  {
221  testFormula();
223  testInfo();
224  infoExample();
226  result = 0;
227  return;
228  }
229  catch (exception& e)
230  {
231  cerr << "Exception in worker thread: " << e.what() << endl;
232  }
233  catch (...)
234  {
235  cerr << "Unhandled exception in worker thread." << endl;
236  }
237  result = 1;
238 }
239 
240 void testThreadSafety(const int& testThreadCount)
241 {
242  boost::barrier testBarrier(testThreadCount);
243  boost::thread_group testThreadGroup;
244  vector<int> results(testThreadCount, 0);
245  for (int i=0; i < testThreadCount; ++i)
246  testThreadGroup.add_thread(new boost::thread(&testThreadSafetyWorker, &testBarrier, boost::ref(results[i])));
247  testThreadGroup.join_all();
248 
249  int failedThreads = std::accumulate(results.begin(), results.end(), 0);
250  if (failedThreads > 0)
251  throw runtime_error(lexical_cast<string>(failedThreads) + " thread failed");
252 }
253 
254 
255 int main(int argc, char* argv[])
256 {
257  TEST_PROLOG(argc, argv)
258 
259  if (argc>1 && !strcmp(argv[1],"-v")) os_ = &cout;
260  if (os_) *os_ << "ChemistryTest\n" << setprecision(12);
261 
262  try
263  {
264  //testThreadSafety(1); // does not test thread-safety of singleton initialization
265  testThreadSafety(2);
266  testThreadSafety(4);
267  testThreadSafety(8);
268  testThreadSafety(16);
269  }
270  catch (exception& e)
271  {
272  TEST_FAILED(e.what())
273  }
274  catch (...)
275  {
276  TEST_FAILED("Caught unknown exception.")
277  }
278 
280 }
281 
282 
void testMassAbundance()
#define unit_assert_throws_what(x, exception, whatStr)
Definition: unit.hpp:119
const TestFormula testFormulaData[]
class to represent a chemical formula
Definition: Chemistry.hpp:133
PWIZ_API_DECL const Record & record(Type type)
returns the amino acid&#39;s Record by type
void testFormula()
void testFormulaOperations()
#define TEST_EPILOG
Definition: unit.hpp:183
#define unit_assert_equal(x, y, epsilon)
Definition: unit.hpp:99
double molecularWeight() const
const int testFormulaDataSize
double monoisotopicMass() const
U
Definition: Chemistry.hpp:80
void testThreadSafetyWorker(boost::barrier *testBarrier, int &result)
int main(int argc, char *argv[])
ostream * os_
Uuh
Definition: Chemistry.hpp:80
void infoExample()
N
Definition: Chemistry.hpp:80
Ca
Definition: Chemistry.hpp:80
double mz(double neutralMass, int protonDelta, int electronDelta=0, int neutronDelta=0)
Definition: Ion.hpp:78
void testPolysiloxane()
void testThreadSafety(const int &testThreadCount)
const char * formula
H
Definition: Chemistry.hpp:80
#define TEST_FAILED(x)
Definition: unit.hpp:177
S
Definition: Chemistry.hpp:80
#define TEST_PROLOG(argc, argv)
Definition: unit.hpp:175
pwiz::chemistry::Element::Type Type
O
Definition: Chemistry.hpp:80
C
Definition: Chemistry.hpp:80
struct for holding isotope information
Definition: Chemistry.hpp:51
#define unit_assert(x)
Definition: unit.hpp:85
std::map< Element::Type, int > Map
Definition: Chemistry.hpp:153
void testInfo()