ProteoWizard
Functions | Variables
IsotopeCalculatorTest.cpp File Reference
#include "IsotopeCalculator.hpp"
#include "pwiz/utility/misc/unit.hpp"
#include "pwiz/utility/misc/Std.hpp"
#include <cstring>

Go to the source code of this file.

Functions

bool nonzero (int a)
 
double probability (const vector< double > &p, const vector< int > &i)
 
void testUsage (const IsotopeCalculator &calc)
 
void testProbabilites (const IsotopeCalculator &calc)
 
void testNormalization (const IsotopeCalculator &calc)
 
int main (int argc, char *argv[])
 

Variables

ostream * os_ = 0
 

Function Documentation

§ nonzero()

bool nonzero ( int  a)
inline

Definition at line 37 of file IsotopeCalculatorTest.cpp.

Referenced by probability().

38 {
39  return a != 0;
40 }

§ probability()

double probability ( const vector< double > &  p,
const vector< int > &  i 
)

Definition at line 44 of file IsotopeCalculatorTest.cpp.

References nonzero().

Referenced by testProbabilites().

45 {
46  if (p.size() < i.size())
47  throw runtime_error("p not big enough");
48 
49  // p = probability distribution
50  // i = partition, sum(i) == n
51 
52  const int n = accumulate(i.begin(), i.end(), 0);
53 
54  // calculate p0^i0 * p1^i1 * ... * pr^ir * C(n; i0, ... , ir)
55 
56  vector<int> p_count = i, C_count = i;
57  int n_count = n;
58 
59  double result = 1;
60 
61  for (int count=0; count<3*n; )
62  {
63  if (n_count && result<=1)
64  {
65  //cout << count << ") multiply: " << n_count << endl;
66  result *= n_count--;
67  count++;
68  }
69 
70  while ((result>=1 || n_count==0) && accumulate(C_count.begin(), C_count.end(), 0))
71  {
72  vector<int>::iterator it = find_if(C_count.begin(), C_count.end(), nonzero);
73  if (it == C_count.end()) throw runtime_error("blech");
74  //cout << count << ") divide: " << *it << endl;
75  result /= (*it)--;
76  count++;
77  }
78 
79  while ((result>=1 || n_count==0) && accumulate(p_count.begin(), p_count.end(), 0))
80  {
81  vector<int>::iterator it = find_if(p_count.begin(), p_count.end(), nonzero);
82  if (it == p_count.end()) throw runtime_error("blech2");
83  size_t index = it - p_count.begin();
84  //cout << count << ") multiply: " << p[index] << endl;
85  result *= p[index];
86  (*it)--;
87  count++;
88  }
89  }
90 
91  return result;
92 }
bool nonzero(int a)

§ testUsage()

void testUsage ( const IsotopeCalculator calc)

Definition at line 95 of file IsotopeCalculatorTest.cpp.

References pwiz::chemistry::IsotopeCalculator::distribution(), and os_.

Referenced by main().

96 {
97  Formula angiotensin("C50 H71 N13 O12 ");
98  Formula bombesin("C71 H110 N24 O18 S1");
99  Formula substanceP("C63 H98 N18 O13 S1");
100  Formula neurotensin("C78 H121 N21 O20");
101  Formula alpha1_6("C45 H59 N11 O8");
102 
103  MassDistribution md = calc.distribution(neurotensin, 2);
104  if (os_) *os_ << "MassDistribution:\n" << md << endl;
105 }
MassDistribution distribution(const Formula &formula, int chargeState=0, int normalization=0) const
class to represent a chemical formula
Definition: Chemistry.hpp:133
ostream * os_
std::vector< MassAbundance > MassDistribution
struct for holding isotope distribution
Definition: Chemistry.hpp:66

§ testProbabilites()

void testProbabilites ( const IsotopeCalculator calc)

Definition at line 108 of file IsotopeCalculatorTest.cpp.

References C, pwiz::chemistry::IsotopeCalculator::distribution(), os_, probability(), pwiz::proteome::AminoAcid::Info::record(), and unit_assert_equal.

Referenced by main().

109 {
110  const MassDistribution& md = Element::Info::record(Element::C).isotopes;
111 
112  if (os_) *os_ << "C distribution: " << md << endl;
113 
114  vector<double> p;
115  for (MassDistribution::const_iterator it=md.begin(); it!=md.end(); ++it)
116  p.push_back(it->abundance);
117 
118  vector<int> neutron0; neutron0.push_back(100);
119  vector<int> neutron1; neutron1.push_back(99); neutron1.push_back(1);
120  vector<int> neutron2; neutron2.push_back(98); neutron2.push_back(2);
121  vector<int> neutron3; neutron3.push_back(97); neutron3.push_back(3);
122  vector<int> neutron4; neutron4.push_back(96); neutron4.push_back(4);
123 
124  vector<double> abundances;
125  abundances.push_back(probability(p, neutron0));
126  abundances.push_back(probability(p, neutron1));
127  abundances.push_back(probability(p, neutron2));
128  abundances.push_back(probability(p, neutron3));
129  abundances.push_back(probability(p, neutron4));
130 
131  if (os_)
132  {
133  *os_ << "C100 abundances (manually calculated):\n";
134  copy(abundances.begin(), abundances.end(), ostream_iterator<double>(*os_, "\n"));
135  *os_ << endl;
136  }
137 
138  MassDistribution md_C100 = calc.distribution(Formula("C100"));
139  if (os_) *os_ << "C100 distribution (from IsotopeCalculator):\n"
140  << md_C100 << endl;
141 
142  // compare manually calculated probabilities to those returned by IsotopeCalculator
143  for (unsigned int i=0; i<abundances.size(); i++)
144  unit_assert_equal(abundances[i], md_C100[i].abundance, 1e-10);
145 }
MassDistribution distribution(const Formula &formula, int chargeState=0, int normalization=0) const
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
#define unit_assert_equal(x, y, epsilon)
Definition: unit.hpp:99
double probability(const vector< double > &p, const vector< int > &i)
ostream * os_
C
Definition: Chemistry.hpp:80
std::vector< MassAbundance > MassDistribution
struct for holding isotope distribution
Definition: Chemistry.hpp:66

§ testNormalization()

void testNormalization ( const IsotopeCalculator calc)

Definition at line 148 of file IsotopeCalculatorTest.cpp.

References pwiz::chemistry::IsotopeCalculator::distribution(), NormalizeMass, os_, unit_assert, and unit_assert_equal.

Referenced by main().

149 {
150  if (os_) *os_ << "mass normalized:\n"
151  << calc.distribution(Formula("C100"), 0,
153 
154  if (os_) *os_ << "abundance normalized:\n"
155  << calc.distribution(Formula("C100"), 0,
156  IsotopeCalculator::NormalizeAbundance) << endl;
157 
158  MassDistribution md = calc.distribution(Formula("C100"), 0,
160  IsotopeCalculator::NormalizeAbundance);
161  if (os_) *os_ << "both normalized:\n" << md << endl;
162 
163  double sumSquares = 0;
164  for (MassDistribution::iterator it=md.begin(); it!=md.end(); ++it)
165  sumSquares += it->abundance * it->abundance;
166  if (os_) *os_ << "sumSquares: " << sumSquares << endl;
167 
168  unit_assert_equal(sumSquares, 1, 1e-12);
169  unit_assert(md[0].mass == 0);
170 }
MassDistribution distribution(const Formula &formula, int chargeState=0, int normalization=0) const
class to represent a chemical formula
Definition: Chemistry.hpp:133
#define unit_assert_equal(x, y, epsilon)
Definition: unit.hpp:99
NormalizeMass
ostream * os_
std::vector< MassAbundance > MassDistribution
struct for holding isotope distribution
Definition: Chemistry.hpp:66
#define unit_assert(x)
Definition: unit.hpp:85

§ main()

int main ( int  argc,
char *  argv[] 
)

Definition at line 173 of file IsotopeCalculatorTest.cpp.

References os_, TEST_EPILOG, TEST_FAILED, TEST_PROLOG, testNormalization(), testProbabilites(), and testUsage().

174 {
175  TEST_PROLOG(argc, argv)
176 
177  try
178  {
179  if (argc>1 && !strcmp(argv[1],"-v")) os_ = &cout;
180  if (os_) *os_ << "IsotopeCalculatorTest\n" << setprecision(12);
181 
182  //IsotopeCalculator calc(1e-10, .2);
183  IsotopeCalculator calc(1e-3, .2);
184  testUsage(calc);
185  testProbabilites(calc);
186  testNormalization(calc);
187  }
188  catch (exception& e)
189  {
190  TEST_FAILED(e.what())
191  }
192  catch (...)
193  {
194  TEST_FAILED("Caught unknown exception.")
195  }
196 
198 }
#define TEST_EPILOG
Definition: unit.hpp:182
void testProbabilites(const IsotopeCalculator &calc)
#define TEST_FAILED(x)
Definition: unit.hpp:176
void testUsage(const IsotopeCalculator &calc)
ostream * os_
#define TEST_PROLOG(argc, argv)
Definition: unit.hpp:174
void testNormalization(const IsotopeCalculator &calc)

Variable Documentation

§ os_

ostream* os_ = 0

Definition at line 34 of file IsotopeCalculatorTest.cpp.

Referenced by main(), testNormalization(), testProbabilites(), and testUsage().