ProteoWizard
IsotopeCalculatorTest.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
24#include "IsotopeCalculator.hpp"
27#include <cstring>
28
29
30using namespace pwiz::util;
31using namespace pwiz::chemistry;
32
33
34ostream* os_ = 0;
35
36
37inline bool nonzero(int a)
38{
39 return a != 0;
40}
41
42
43// calculate multinomial probabilities manually
44double probability(const vector<double>& p, const vector<int>& i)
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}
93
94
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}
106
107
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}
146
147
149{
150 if (os_) *os_ << "mass normalized:\n"
151 << calc.distribution(Formula("C100"), 0,
152 IsotopeCalculator::NormalizeMass) << endl;
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,
159 IsotopeCalculator::NormalizeMass |
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}
171
172
173int main(int argc, char* argv[])
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}
bool nonzero(int a)
int main(int argc, char *argv[])
void testProbabilites(const IsotopeCalculator &calc)
double probability(const vector< double > &p, const vector< int > &i)
ostream * os_
void testNormalization()
class to represent a chemical formula
MassDistribution distribution(const Formula &formula, int chargeState=0, int normalization=0) const
PWIZ_API_DECL const Record & record(Type type)
retrieve the record for an element
std::vector< MassAbundance > MassDistribution
struct for holding isotope distribution
Definition Chemistry.hpp:66
MassDistribution isotopes
the most abundant isotope
#define unit_assert(x)
Definition unit.hpp:85
#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