1 /*
2 * ===========================================================================
3 *
4 * PUBLIC DOMAIN NOTICE
5 * National Center for Biotechnology Information
6 *
7 * This software/database is a "United States Government Work" under the
8 * terms of the United States Copyright Act. It was written as part of
9 * the author's official duties as a United States Government employee and
10 * thus cannot be copyrighted. This software/database is freely available
11 * to the public for use. The National Library of Medicine and the U.S.
12 * Government have not placed any restriction on its use or reproduction.
13 *
14 * Although all reasonable efforts have been taken to ensure the accuracy
15 * and reliability of the software and data, the NLM and the U.S.
16 * Government do not and cannot warrant the performance or results that
17 * may be obtained by using this software or data. The NLM and the U.S.
18 * Government disclaim all warranties, express or implied, including
19 * warranties of performance, merchantability or fitness for any particular
20 * purpose.
21 *
22 * Please cite the authors in any work or product based on this material.
23 *
24 * ===========================================================================
25 *
26 * Author: Lewis Y. Geer
27 *
28 * File Description:
29 * code to do mass spec scoring
30 *
31 * ===========================================================================
32 */
33
34 //
35 //
36 //
37 // Please keep this file clear of NCBI toolkit dependent code as it
38 // is meant to be used externally from the toolkit.
39 //
40 // NCBI dependent code can go into omssascore.cpp
41 //
42 //
43 //
44
45 // precompiled headers #include found in msscore.cpp
46 #if 0
47 #include <ncbi_pch.hpp>
48 #endif
49
50 #ifdef HAVE_LIMITS
51 #include <limits>
52 #endif
53 #include <math.h>
54 #include "msscore.hpp"
55
56 #ifdef USING_OMSSA_SCOPE
57 USING_OMSSA_SCOPE
58 #endif
59
60 // define log of gamma fcn if not defined already
61 #ifndef MSLNGAMMA
62 #define MSLNGAMMA lgamma
63 #endif
64
65 // define error function
66 #ifndef MSERF
67 #define MSERF erf
68 #endif
69
70
IsForwardSeries(EMSIonSeries Series)71 const int CMSBasicMatchedPeak::IsForwardSeries(EMSIonSeries Series)
72 {
73 if ( Series == eMSIonTypeA ||
74 Series == eMSIonTypeB ||
75 Series == eMSIonTypeC)
76 return 1;
77 if ( Series == eMSIonTypeX ||
78 Series == eMSIonTypeY ||
79 Series == eMSIonTypeZ)
80 return 0;
81 return -1;
82 }
83
CMSMatchedPeak(void)84 CMSMatchedPeak::CMSMatchedPeak(void):
85 MassTolerance(-1),
86 ExpIons(-1),
87 MatchType(eMSMatchTypeUnknown)
88 {
89 }
90
Assign(CMSMatchedPeak * in)91 void CMSMatchedPeak::Assign(CMSMatchedPeak *in)
92 {
93 if (!in || in == this) return;
94 SetCharge() = in->GetCharge();
95 SetExpIons() = in->GetExpIons();
96 SetIntensity() = in->GetIntensity();
97 SetIonSeries() = in->GetIonSeries();
98 SetMassTolerance() = in->GetMassTolerance();
99 SetMatchType() = in->GetMatchType();
100 SetMZ() = in->GetMZ();
101 SetNumber() = in->GetNumber();
102 }
103
Assign(CMSBasicMatchedPeak * in)104 void CMSMatchedPeak::Assign(CMSBasicMatchedPeak *in)
105 {
106 if (!in || in == this) return;
107 SetCharge() = in->GetCharge();
108 SetIntensity() = in->GetIntensity();
109 SetIonSeries() = in->GetIonSeries();
110 SetMZ() = in->GetMZ();
111 SetNumber() = in->GetNumber();
112
113 SetMassTolerance() = -1;
114 SetExpIons() = -1;
115 SetMatchType() = eMSMatchTypeUnknown;
116 }
117
118
CMSMatchedPeakSet(void)119 CMSMatchedPeakSet::CMSMatchedPeakSet(void)
120 {
121 }
122
~CMSMatchedPeakSet()123 CMSMatchedPeakSet::~CMSMatchedPeakSet()
124 {
125 DeleteMatchedPeakSet();
126 }
127
CreateMatchedPeakSet(int SizeIn)128 void CMSMatchedPeakSet::CreateMatchedPeakSet(int SizeIn)
129 {
130 DeleteMatchedPeakSet();
131 int i;
132 for (i = 0; i < SizeIn; ++i)
133 SetMatchedPeakSet().push_back(new CMSMatchedPeak);
134 }
135
DeleteMatchedPeakSet(void)136 void CMSMatchedPeakSet::DeleteMatchedPeakSet(void)
137 {
138 while (!GetMatchedPeakSet().empty()) {
139 delete SetMatchedPeakSet().back();
140 SetMatchedPeakSet().pop_back();
141 }
142 }
143
144 void
Compare(CMSMatchedPeakSet * Other,bool SameDirection)145 CMSMatchedPeakSet::Compare(CMSMatchedPeakSet *Other, bool SameDirection)
146 {
147 if(!Other ||
148 (GetMatchedPeakSet().size() != Other->GetMatchedPeakSet().size()) ) return;
149
150 unsigned i, j;
151
152 for(i = 0; i < GetMatchedPeakSet().size(); ++i) {
153 if(SameDirection) j = i;
154 else j = Other->GetSize() - i - 1;
155 // don't do the comparison if there are no peaks to compare!
156 if(j >= Other->GetMatchedPeakSet().size()) break;
157
158 if ((GetMatchedPeakSet()[i]->GetMatchType() == eMSMatchTypeIndependent ||
159 GetMatchedPeakSet()[i]->GetMatchType() == eMSMatchTypeSemiIndependent) &&
160 (Other->GetMatchedPeakSet()[j]->GetMatchType() == eMSMatchTypeIndependent ||
161 Other->GetMatchedPeakSet()[j]->GetMatchType() == eMSMatchTypeSemiIndependent)) {
162 Other->GetMatchedPeakSet()[j]->SetMatchType() = eMSMatchTypeDependent;
163 }
164 }
165 }
166
CMSMatchedPeakSetMap(void)167 CMSMatchedPeakSetMap::CMSMatchedPeakSetMap(void)
168 {
169 }
170
~CMSMatchedPeakSetMap()171 CMSMatchedPeakSetMap::~CMSMatchedPeakSetMap()
172 {
173 TIonSeriesMatchMap::iterator i;
174 for ( i = SetMatchMap().begin(); i != SetMatchMap().end(); ++i)
175 delete i->second;
176 SetMatchMap().clear();
177 }
178
179 CMSMatchedPeakSet *
CreateSeries(TMSCharge Charge,TMSIonSeries Series,int Size,int Maxproductions)180 CMSMatchedPeakSetMap::CreateSeries(
181 TMSCharge Charge,
182 TMSIonSeries Series,
183 int Size,
184 int Maxproductions)
185 {
186 CMSMatchedPeakSet *retval;
187 int Key = ChargeSeries2Key(Charge, Series);
188 unsigned Realsize = min(Size, Maxproductions);
189
190 // find old matches
191 if (GetMatchMap().find(Key) != GetMatchMap().end()) {
192 // if size is same, don't bother reallocating
193 if (SetMatchMap()[Key]->SetMatchedPeakSet().size() == Realsize) {
194 return SetMatchMap()[Key];
195 }
196 // otherwise delete and reallocate
197 else {
198 delete SetMatchMap()[Key];
199 SetMatchMap().erase(Key);
200 }
201 }
202
203 retval = new CMSMatchedPeakSet;
204 if (retval) {
205 retval->CreateMatchedPeakSet(Realsize);
206 retval->SetSize() = Size;
207 SetMatchMap().insert(std::pair <int, CMSMatchedPeakSet * > (Key, retval));
208 }
209
210 return retval;
211 }
212
SetSeries(TMSCharge Charge,TMSIonSeries Series)213 CMSMatchedPeakSet * CMSMatchedPeakSetMap::SetSeries(TMSCharge Charge, TMSIonSeries Series)
214 {
215 CMSMatchedPeakSet *retval;
216 const int Key = ChargeSeries2Key(Charge, Series);
217 if (GetMatchMap().find(Key) != GetMatchMap().end())
218 retval = SetMatchMap()[Key];
219 else
220 retval = 0;
221
222 return retval;
223 }
224
225 //typedef std::map <int, CMSMatchedPeakSet *> TIonSeriesMatchMap;
226
227 const int
ChargeSeries2Key(TMSCharge Charge,TMSIonSeries Series)228 CMSMatchedPeakSetMap::ChargeSeries2Key(TMSCharge Charge, TMSIonSeries Series)
229 {
230 return eMSIonTypeMax * Charge + Series;
231 }
232
233 const TMSCharge
Key2Charge(int Key)234 CMSMatchedPeakSetMap::Key2Charge(int Key)
235 {
236 return static_cast <TMSCharge> (Key/eMSIonTypeMax);
237 }
238
239 const TMSIonSeries
Key2Series(int Key)240 CMSMatchedPeakSetMap::Key2Series(int Key)
241 {
242 return static_cast <TMSIonSeries> (Key - Key2Charge(Key)*eMSIonTypeMax);
243 }
244
245
246
CMSSpectrumMatch(void)247 CMSSpectrumMatch::CMSSpectrumMatch(void):
248 HitInfo(0),
249 Hits(0),
250 ExpMass(-1),
251 TheoreticalMass(-1),
252 Sum(0),
253 M(0),
254 N(0)
255 {
256 }
257
~CMSSpectrumMatch()258 CMSSpectrumMatch::~CMSSpectrumMatch()
259 {
260 delete [] HitInfo;
261 }
262
CreateHitInfo(void)263 void CMSSpectrumMatch::CreateHitInfo(void)
264 {
265 delete [] HitInfo;
266 HitInfo = 0;
267 if (GetHits() > 0)
268 HitInfo = new CMSBasicMatchedPeak[GetHits()];
269 }
270
Find(TMSNumber Number,TMSCharge ChargeIn,TMSIonSeries Series)271 CMSBasicMatchedPeak * CMSSpectrumMatch::Find(
272 TMSNumber Number,
273 TMSCharge ChargeIn,
274 TMSIonSeries Series)
275 {
276 int i;
277 for (i = 0; i < GetHits(); ++i) {
278 if (GetHitInfo(i).GetNumber() == Number &&
279 GetHitInfo(i).GetCharge() == ChargeIn &&
280 GetHitInfo(i).GetIonSeries() == Series)
281 return &SetHitInfo(i);
282 }
283 return 0;
284 }
285
286
FillMatchedPeaks(TMSCharge ChargeIn,TMSIonSeries Series,unsigned Size,TMSIntensity MinIntensity,bool Skipb1,EMSTerminalBias TerminalIon,int Maxproductions,string & Sequence,bool NoProline)287 void CMSSpectrumMatch::FillMatchedPeaks(
288 TMSCharge ChargeIn,
289 TMSIonSeries Series,
290 unsigned Size,
291 TMSIntensity MinIntensity,
292 bool Skipb1,
293 EMSTerminalBias TerminalIon,
294 int Maxproductions
295 //#if 0
296 ,
297 string &Sequence,
298 bool NoProline
299 //#endif
300 )
301 {
302 // create a new match set
303 CMSMatchedPeakSet *MatchPeakSet =
304 SetIonSeriesMatchMap().CreateSeries(ChargeIn, Series, Size, Maxproductions);
305
306 unsigned i;
307
308 // iterate over series and look for matching hits
309 for (i = 0; i < MatchPeakSet->GetMatchedPeakSet().size(); ++i) {
310 CMSBasicMatchedPeak * FoundPeak = Find(i, ChargeIn, Series);
311 if (FoundPeak && FoundPeak->GetIntensity() > MinIntensity) {
312 // copy hit to match
313 MatchPeakSet->SetMatchedPeakSet()[i]->Assign(FoundPeak);
314 // set as untyped hit
315 MatchPeakSet->SetMatchedPeakSet()[i]->SetMatchType() =
316 eMSMatchTypeNotTyped;
317
318 // check for terminal match
319 // n-term
320 if (TerminalIon == eMSNTerminalBias || TerminalIon == eMSBothTerminalBias) {
321 if ( (i==0 && CMSBasicMatchedPeak::IsForwardSeries(static_cast <EMSIonSeries> (Series))) ||
322 (i==Size-1 && !CMSBasicMatchedPeak::IsForwardSeries(static_cast <EMSIonSeries> (Series)))) {
323 MatchPeakSet->SetMatchedPeakSet()[i]->SetMatchType() = eMSMatchTypeTerminus;
324 }
325 }
326 // c-term
327 else if (TerminalIon == eMSCTerminalBias || TerminalIon == eMSBothTerminalBias) {
328 if ( (i==0 && !CMSBasicMatchedPeak::IsForwardSeries(static_cast <EMSIonSeries> (Series))) ||
329 (i==Size-1 && CMSBasicMatchedPeak::IsForwardSeries(static_cast <EMSIonSeries> (Series)))) {
330 MatchPeakSet->SetMatchedPeakSet()[i]->SetMatchType() = eMSMatchTypeTerminus;
331 }
332 }
333 }
334 else {
335 MatchPeakSet->SetMatchedPeakSet()[i]->SetCharge() = ChargeIn;
336 MatchPeakSet->SetMatchedPeakSet()[i]->SetIonSeries() = Series;
337 MatchPeakSet->SetMatchedPeakSet()[i]->SetMatchType() = eMSMatchTypeUnknown;
338 MatchPeakSet->SetMatchedPeakSet()[i]->SetMZ() = 0;
339 MatchPeakSet->SetMatchedPeakSet()[i]->SetIntensity() = 0;
340 MatchPeakSet->SetMatchedPeakSet()[i]->SetNumber() = i;
341 }
342 // if b1 should not be searched, mark it as so.
343 if (i == 0 && Skipb1 && CMSBasicMatchedPeak::IsForwardSeries(static_cast <EMSIonSeries> (Series)))
344 MatchPeakSet->SetMatchedPeakSet()[i]->SetMatchType() = eMSMatchTypeNoSearch;
345 //#if 0
346 // apply rule that no
347 if (NoProline &&
348 CMSBasicMatchedPeak::IsForwardSeries(static_cast <EMSIonSeries> (Series)) &&
349 i+1 < Sequence.size() &&
350 Sequence[i+1] == 'P') {
351 MatchPeakSet->SetMatchedPeakSet()[i]->SetMatchType() = eMSMatchTypeNoSearch;
352 }
353
354 if (NoProline &&
355 !CMSBasicMatchedPeak::IsForwardSeries(static_cast <EMSIonSeries> (Series)) &&
356 i < Sequence.size() &&
357 Sequence[Sequence.size() - i - 1] == 'P') {
358 MatchPeakSet->SetMatchedPeakSet()[i]->SetMatchType() = eMSMatchTypeNoSearch;
359 }
360 //#endif
361
362 }
363
364
365 // iterate over array, look for unset m/z. keep track of index of last set m/z
366 // if unset, iterate until first set, then fill in mz
367
368 // find first -1
369 unsigned j;
370 for (i = 0; i < MatchPeakSet->GetMatchedPeakSet().size(); ++i) {
371 if (MatchPeakSet->GetMatchedPeakSet()[i]->GetMZ() != 0) break;
372 }
373 // if first filled hit is not the first one, fill in array
374 for (j = 0; j < i && i < MatchPeakSet->GetMatchedPeakSet().size(); ++j) {
375 MatchPeakSet->SetMatchedPeakSet()[j]->SetMZ() =
376 (MatchPeakSet->GetMatchedPeakSet()[i]->GetMZ()/ (i+1)) * (j+1);
377 }
378 // now fill out sections in the center
379 unsigned LastIndex(i);
380 for (; i < MatchPeakSet->GetMatchedPeakSet().size(); ++i) {
381 if (MatchPeakSet->GetMatchedPeakSet()[i]->GetMZ() != 0) {
382 if (i != LastIndex) {
383 // fill in the blank spaces
384 for (j = LastIndex+1; j < i; j++) {
385 MatchPeakSet->SetMatchedPeakSet()[j]->SetMZ() =
386 MatchPeakSet->GetMatchedPeakSet()[LastIndex]->GetMZ() +
387 ((MatchPeakSet->GetMatchedPeakSet()[i]->GetMZ() -
388 MatchPeakSet->GetMatchedPeakSet()[LastIndex]->GetMZ()) /
389 (i - LastIndex )) * (j - LastIndex) ;
390 }
391 }
392 LastIndex = i;
393 }
394 }
395
396 // now fill out last section
397 if (MatchPeakSet->GetMatchedPeakSet()[MatchPeakSet->GetMatchedPeakSet().size()-1]->GetMZ() == 0) {
398
399 int StartIndex;
400 TMSMZ StartMass;
401 TMSMZ MassIncrement;
402 // if completely empty, force the array to be completely filled out
403 if (MatchPeakSet->GetMatchedPeakSet()[0]->GetMZ() == 0) {
404 MassIncrement = (GetExpMass()/ChargeIn)/(Size+1);
405 StartIndex = 0;
406 StartMass = (GetExpMass()/ChargeIn)/(Size+1);
407 }
408 else {
409 MassIncrement = ((GetExpMass()/ChargeIn) - MatchPeakSet->GetMatchedPeakSet()[LastIndex]->GetMZ()) /
410 (Size - LastIndex);
411 StartIndex = LastIndex + 1;
412 StartMass = MatchPeakSet->GetMatchedPeakSet()[LastIndex]->GetMZ() +
413 MassIncrement;
414 }
415
416 for (j = StartIndex; j < MatchPeakSet->GetMatchedPeakSet().size(); ++j) {
417 MatchPeakSet->SetMatchedPeakSet()[j]->SetMZ() =
418 StartMass + MassIncrement * (j - StartIndex);
419 }
420 }
421
422
423 // first determine if the first peak is independent. This simplifies the loop
424 if (MatchPeakSet->GetMatchedPeakSet()[0]->GetMatchType() == eMSMatchTypeNotTyped)
425 MatchPeakSet->SetMatchedPeakSet()[0]->SetMatchType() = eMSMatchTypeIndependent;
426
427 // now characterize all of the matches depending on the preceeding match
428 for (i = 1; i < MatchPeakSet->GetMatchedPeakSet().size(); ++i) {
429 // only type matched peaks
430 if (MatchPeakSet->GetMatchedPeakSet()[i]->GetMatchType() == eMSMatchTypeNotTyped) {
431 if (MatchPeakSet->GetMatchedPeakSet()[i-1]->GetMatchType() == eMSMatchTypeIndependent ||
432 MatchPeakSet->GetMatchedPeakSet()[i-1]->GetMatchType() == eMSMatchTypeSemiIndependent) {
433 MatchPeakSet->SetMatchedPeakSet()[i]->SetMatchType() = eMSMatchTypeDependent;
434 }
435 else if (MatchPeakSet->GetMatchedPeakSet()[i-1]->GetMatchType() == eMSMatchTypeNoSearch ||
436 MatchPeakSet->GetMatchedPeakSet()[i-1]->GetMatchType() == eMSMatchTypeUnknown ||
437 MatchPeakSet->GetMatchedPeakSet()[i-1]->GetMatchType() == eMSMatchTypeTerminus) {
438 MatchPeakSet->SetMatchedPeakSet()[i]->SetMatchType() = eMSMatchTypeIndependent;
439 }
440 else if (MatchPeakSet->GetMatchedPeakSet()[i-1]->GetMatchType() == eMSMatchTypeDependent) {
441 MatchPeakSet->SetMatchedPeakSet()[i]->SetMatchType() = eMSMatchTypeSemiIndependent;
442 }
443 }
444 }
445
446 return;
447 }
448
449
450
451 const double
CalcPoissonMean(double ProbTerminal,int NumTerminalMasses,double ProbDependent,int NumUniqueMasses,double ToleranceAdjust) const452 CMSSpectrumMatch::CalcPoissonMean(double ProbTerminal,
453 int NumTerminalMasses,
454 double ProbDependent,
455 int NumUniqueMasses,
456 double ToleranceAdjust
457 ) const
458 {
459 double Mean(0.0L);
460 // scan thru entire map
461 TIonSeriesMatchMap::const_iterator i;
462 for (i = GetIonSeriesMatchMap().GetMatchMap().begin();
463 i != GetIonSeriesMatchMap().GetMatchMap().end();
464 ++i) {
465 CMSMatchedPeakSet *PeakSet = i->second;
466 // scan thru each match set
467 TMatchedPeakSet::const_iterator j;
468 for (j = PeakSet->GetMatchedPeakSet().begin();
469 j != PeakSet->GetMatchedPeakSet().end();
470 ++j) {
471 CMSMatchedPeak *Peak = *j;
472
473 // skip unsearched peaks
474 if (Peak->GetMatchType() == eMSMatchTypeNoSearch) continue;
475
476 // calculate single poisson
477 // add it to total poisson
478 Mean += 2.0 * Peak->GetMassTolerance() * Peak->GetExpIons() * ToleranceAdjust;
479 _TRACE( "mean=" << Mean << " masstol=" << Peak->GetMassTolerance() << " expions=" << Peak->GetExpIons());
480
481 // if a statistically dependent peak, add in dependent match probability
482 if (Peak->GetMatchType() == eMSMatchTypeDependent && NumUniqueMasses != 0) {
483 Mean += ProbDependent / NumUniqueMasses;
484 _TRACE( "mean=" << Mean << " probdep=" << ProbDependent << " numunique=" << NumUniqueMasses);
485 }
486
487 // if a statistically biased terminal peak, add in biased match probability
488 if (Peak->GetMatchType() == eMSMatchTypeTerminus && NumTerminalMasses != 0) {
489 Mean += ProbTerminal / NumTerminalMasses;
490 }
491 }
492 }
493 return Mean;
494 }
495
CalcPoisson(double Mean,int i) const496 const double CMSSpectrumMatch::CalcPoisson(double Mean, int i) const
497 {
498 return exp(-Mean) * pow(Mean, i) / exp(MSLNGAMMA(i+1));
499 }
500
CalcPoissonTopHit(double Mean,int i,double TopHitProb) const501 const double CMSSpectrumMatch::CalcPoissonTopHit(double Mean, int i, double TopHitProb) const
502 {
503 double retval;
504 retval = exp(-Mean) * pow(Mean, i) / exp(MSLNGAMMA(i+1));
505 if (TopHitProb != 1.0L) retval *= 1.0L - pow((1.0L-TopHitProb), i);
506 return retval;
507 }
508
CalcNormalTopHit(double Mean,double TopHitProb) const509 const double CMSSpectrumMatch::CalcNormalTopHit(double Mean, double TopHitProb) const
510 {
511 int i;
512 double retval(0.0L), before(-1.0L), increment;
513
514 for (i = 1; i < 1000; i++) {
515 increment = CalcPoissonTopHit(Mean, i, TopHitProb);
516 // convergence hack -- on linux (at least) convergence test doesn't work
517 // for gcc release build
518 if (increment <= MSDOUBLELIMIT && i > 10) break;
519 retval += increment;
520 if (retval == before) break; // convergence
521 before = retval;
522 }
523 return retval;
524 }
525
CalcPvalue(double Mean,int HitsIn) const526 const double CMSSpectrumMatch::CalcPvalue(double Mean, int HitsIn) const
527 {
528 if (HitsIn <= 0) return 1.0L;
529
530 int i;
531 double retval(0.0L), increment, beforeincrement(0.0L);
532
533 for (i = 0; i < HitsIn; i++) {
534 increment = CalcPoisson(Mean, i);
535 #ifdef HAVE_LIMITS
536 // detect limit of increment if numeric_limits available
537 if(increment < retval * std::numeric_limits <double> ::epsilon()) break;
538 #endif
539 retval += increment;
540 if(retval == 1.0L) break;
541 beforeincrement = increment;
542 }
543
544 // if retval is calculated to the maximum precision of double
545 // back off one increment to give a lower bound of the e-value
546 if(retval == 1.0L)
547 retval -= beforeincrement;
548
549 retval = 1.0L - retval;
550 if (retval <= 0.0L) retval = MSDOUBLELIMIT;
551 return retval;
552 }
553
CalcPvalueTopHit(double Mean,int HitsIn,double Normal,double TopHitProb) const554 const double CMSSpectrumMatch::CalcPvalueTopHit(double Mean, int HitsIn, double Normal, double TopHitProb) const
555 {
556 if (HitsIn <= 0) return 1.0L;
557
558 int i;
559 double retval(0.0L), increment, beforeretval(-1.0), beforeincrement(0.0L);
560
561 for (i = 1; i < HitsIn; i++) {
562 increment = CalcPoissonTopHit(Mean, i, TopHitProb);
563 #ifdef HAVE_LIMITS
564 // detect limit of increment if numeric_limits available
565 if(increment < retval * std::numeric_limits <double> ::epsilon()) break;
566 #endif
567 retval += increment;
568 // note that in gcc 3.4.2 on linux with -O9 -ffast-math, the following statement doesn't appear to work
569 // trying to fix this with retval - beforeretval < increment also does not work -- it breaks too soon!
570 if(retval == beforeretval) break;
571 beforeretval = retval;
572 beforeincrement = increment;
573 }
574
575 // if retval and Normal are calculated to the maximum precision of double
576 // back off one increment to give a lower bound of the e-value
577 if(retval == Normal) retval -= beforeincrement;
578
579 retval /= Normal;
580 retval = 1.0L - retval;
581 if (retval <= 0.0L) retval = MSDOUBLELIMIT;
582 return retval;
583 }
584
CalcRankProb(void) const585 const double CMSSpectrumMatch::CalcRankProb(void) const
586 {
587 double Average, StdDev, Perf(1.0L);
588 if (GetM() != 0.0) {
589 Average = (GetM() *(GetN()+1))/2.0;
590 StdDev = sqrt(Average * (GetN() - GetM()) / 6.0);
591 if ( StdDev != 0.0)
592 Perf = 0.5*(1.0 + MSERF((GetSum() - Average)/(StdDev*sqrt(2.0))));
593 }
594 return Perf;
595 }
596
597
598 const TMSMZ
GetMeanDelta(void) const599 CMSSpectrumMatch::GetMeanDelta(void) const
600 {
601 TMSMZ Mean(0);
602
603 int i;
604 for (i = 0; i < GetHits(); ++i) {
605 Mean += GetHitInfo(i).GetDelta();
606 }
607 return static_cast <int> (Mean/(static_cast <double> (GetHits())));
608 }
609
610
611 const TMSMZ
GetStdDevDelta(void) const612 CMSSpectrumMatch::GetStdDevDelta(void) const
613 {
614 TMSMZ Mean;
615 double StdDev(0.0L);
616
617 Mean = GetMeanDelta();
618
619 int i;
620 for (i = 0; i < GetHits(); ++i) {
621 StdDev += pow((double)GetHitInfo(i).GetDelta() - Mean, 2.0);
622 }
623 StdDev = pow((double)StdDev/GetHits(), 0.5);
624 return (TMSMZ)StdDev;
625 }
626
627
628 const TMSMZ
GetMaxDelta(void) const629 CMSSpectrumMatch::GetMaxDelta(void) const
630 {
631 TMSMZ Max(0);
632 int i;
633 for (i = 0; i < GetHits(); ++i) {
634 if(abs(GetHitInfo(i).GetDelta()) > Max)
635 Max = abs(GetHitInfo(i).GetDelta());
636 }
637 return Max;
638 }
639
640