← Back to branch summary

~ubuntu-branches/ubuntu/wily/openms/wily

« back to all changes in this revision

Viewing changes to source/COMPARISON/SPECTRA/SpectrumAlignmentScore.C

  • Committer: Package Import Robot
  • Author(s): Filippo Rusconi
  • Date: 2012-11-12 15:58:12 UTC
  • Revision ID: package-import@ubuntu.com-20121112155812-vr15wtg9b50cuesg
Tags: upstream-1.9.0
Import upstream version 1.9.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
source/COMPARISON/SPECTRA/SpectrumAlignmentScore.C
 
1
// -*- mode: C++; tab-width: 2; -*-
 
2
// vi: set ts=2:
 
3
//
 
4
// --------------------------------------------------------------------------
 
5
//                   OpenMS Mass Spectrometry Framework
 
6
// --------------------------------------------------------------------------
 
7
//  Copyright (C) 2003-2011 -- Oliver Kohlbacher, Knut Reinert
 
8
//
 
9
//  This library is free software; you can redistribute it and/or
 
10
//  modify it under the terms of the GNU Lesser General Public
 
11
//  License as published by the Free Software Foundation; either
 
12
//  version 2.1 of the License, or (at your option) any later version.
 
13
//
 
14
//  This library is distributed in the hope that it will be useful,
 
15
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
 
16
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 
17
//  Lesser General Public License for more details.
 
18
//
 
19
//  You should have received a copy of the GNU Lesser General Public
 
20
//  License along with this library; if not, write to the Free Software
 
21
//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
22
//
 
23
// --------------------------------------------------------------------------
 
24
// $Maintainer: Andreas Bertsch $
 
25
// $Authors: Andreas Bertsch $
 
26
// --------------------------------------------------------------------------
 
27
 
 
28
#include <OpenMS/COMPARISON/SPECTRA/SpectrumAlignmentScore.h>
 
29
#include <OpenMS/COMPARISON/SPECTRA/SpectrumAlignment.h>
 
30
#include <cmath>
 
31
 
 
32
#include <boost/math/special_functions/erf.hpp>
 
33
 
 
34
using namespace std;
 
35
 
 
36
namespace OpenMS
 
37
{
 
38
  SpectrumAlignmentScore::SpectrumAlignmentScore()
 
39
    : PeakSpectrumCompareFunctor()
 
40
  {
 
41
                setName(SpectrumAlignmentScore::getProductName());
 
42
                defaults_.setValue("tolerance", 0.3, "Defines the absolut (in Da) or relative (in ppm) tolerance");
 
43
                defaults_.setValue("is_relative_tolerance", "false", "if true, the tolerance value is interpreted as ppm");
 
44
                defaults_.setValidStrings("is_relative_tolerance", StringList::create("true,false"));
 
45
                defaults_.setValue("use_linear_factor", "false", "if true, the intensities are weighted with the relative m/z difference");
 
46
                defaults_.setValidStrings("use_linear_factor", StringList::create("true,false"));
 
47
                defaults_.setValue("use_gaussian_factor", "false", "if true, the intensities are weighted with the relative m/z difference using a gaussian");
 
48
                defaults_.setValidStrings("use_gaussian_factor", StringList::create("true,false"));
 
49
                defaultsToParam_();
 
50
  }
 
51
 
 
52
  SpectrumAlignmentScore::SpectrumAlignmentScore(const SpectrumAlignmentScore& source)
 
53
    : PeakSpectrumCompareFunctor(source)
 
54
  {
 
55
  }
 
56
 
 
57
  SpectrumAlignmentScore::~SpectrumAlignmentScore()
 
58
  {
 
59
  }
 
60
 
 
61
  SpectrumAlignmentScore& SpectrumAlignmentScore::operator = (const SpectrumAlignmentScore& source)
 
62
  {
 
63
                if (this != &source)
 
64
                {
 
65
        PeakSpectrumCompareFunctor::operator = (source);
 
66
                }
 
67
    return *this;
 
68
  }
 
69
 
 
70
        DoubleReal SpectrumAlignmentScore::operator () (const PeakSpectrum& spec) const
 
71
        {
 
72
                return operator () (spec, spec);
 
73
        }
 
74
 
 
75
  DoubleReal SpectrumAlignmentScore::operator () (const PeakSpectrum& s1, const PeakSpectrum& s2) const
 
76
  {                     
 
77
                const DoubleReal tolerance = (DoubleReal)param_.getValue("tolerance");
 
78
                bool is_relative_tolerance = param_.getValue("is_relative_tolerance").toBool();
 
79
                bool use_linear_factor = param_.getValue("use_linear_factor").toBool();
 
80
                bool use_gaussian_factor = param_.getValue("use_gaussian_factor").toBool();
 
81
 
 
82
                if (use_linear_factor && use_gaussian_factor)
 
83
                {
 
84
                        cerr << "Warning: SpectrumAlignmentScore, use either 'use_linear_factor' or 'use_gaussian_factor'!" << endl;
 
85
                }
 
86
 
 
87
                SpectrumAlignment aligner;
 
88
                Param p;
 
89
                p.setValue("tolerance", tolerance);
 
90
                p.setValue("is_relative_tolerance", (String)param_.getValue("is_relative_tolerance"));
 
91
                aligner.setParameters(p);
 
92
 
 
93
                vector<pair<Size, Size> > alignment;
 
94
                aligner.getSpectrumAlignment(alignment, s1, s2);
 
95
 
 
96
                DoubleReal score(0), sum(0), sum1(0), sum2(0);
 
97
                for (PeakSpectrum::ConstIterator it1 = s1.begin(); it1 != s1.end(); ++it1)
 
98
                {
 
99
                        sum1 += it1->getIntensity() * it1->getIntensity();
 
100
                }
 
101
                
 
102
                for (PeakSpectrum::ConstIterator it1 = s2.begin(); it1 != s2.end(); ++it1)
 
103
                {
 
104
                        sum2 += it1->getIntensity() * it1->getIntensity();
 
105
                }
 
106
                
 
107
                for (vector<pair<Size, Size> >::const_iterator it = alignment.begin(); it != alignment.end(); ++it)
 
108
                {
 
109
                        //DoubleReal factor(0.0);
 
110
                        //factor = (epsilon - fabs(s1[it->first].getPosition()[0] - s2[it->second].getPosition()[0])) / epsilon;
 
111
                        DoubleReal mz_tolerance(tolerance);
 
112
 
 
113
                        if (is_relative_tolerance)
 
114
                        {
 
115
                                mz_tolerance = mz_tolerance * s1[it->first].getPosition()[0] / 1e6;
 
116
                        }
 
117
        
 
118
                        DoubleReal mz_difference(fabs(s1[it->first].getPosition()[0] - s2[it->second].getPosition()[0]));
 
119
                        DoubleReal factor = 1.0;
 
120
                        
 
121
                        if (use_linear_factor || use_gaussian_factor)
 
122
                        {
 
123
                                factor = getFactor_(mz_tolerance, mz_difference, use_gaussian_factor);
 
124
                        }
 
125
                        sum += sqrt(s1[it->first].getIntensity() * s2[it->second].getIntensity() * factor);
 
126
                }
 
127
 
 
128
    score = sum / (sqrt(sum1 * sum2));
 
129
 
 
130
    return score;
 
131
        }
 
132
 
 
133
        DoubleReal SpectrumAlignmentScore::getFactor_(DoubleReal mz_tolerance, DoubleReal mz_difference, bool is_gaussian) const
 
134
        {
 
135
                DoubleReal factor(0.0);
 
136
 
 
137
                if (is_gaussian)
 
138
                {
 
139
                        static const DoubleReal denominator = mz_tolerance * 3.0 * sqrt(2.0);
 
140
                        factor = boost::math::erfc(mz_difference / denominator);
 
141
                        //cerr << "Factor: " << factor << " " << mz_tolerance << " " << mz_difference << endl;
 
142
                }
 
143
                else
 
144
                {
 
145
                        factor = (mz_tolerance - mz_difference) / mz_tolerance;
 
146
                }
 
147
                return factor;
 
148
        }
 
149
 
 
150
}