← Back to branch summary

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

« back to all changes in this revision

Viewing changes to source/ANALYSIS/ID/HiddenMarkovModel.C

  • Committer: Package Import Robot
  • Author(s): Filippo Rusconi
  • Date: 2013-12-20 11:30:16 UTC
  • mfrom: (5.1.2 sid)
  • Revision ID: package-import@ubuntu.com-20131220113016-wre5g9bteeheq6he
Tags: 1.11.1-3
* remove version number from libbost development package names;
* ensure that AUTHORS is correctly shipped in all packages.

Show diffs side-by-side

added added

removed removed

Lines of Context:
source/ANALYSIS/ID/HiddenMarkovModel.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
 
1
// --------------------------------------------------------------------------
 
2
//                   OpenMS -- Open-Source Mass Spectrometry
 
3
// --------------------------------------------------------------------------
 
4
// Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
 
5
// ETH Zurich, and Freie Universitaet Berlin 2002-2013.
 
6
//
 
7
// This software is released under a three-clause BSD license:
 
8
//  * Redistributions of source code must retain the above copyright
 
9
//    notice, this list of conditions and the following disclaimer.
 
10
//  * Redistributions in binary form must reproduce the above copyright
 
11
//    notice, this list of conditions and the following disclaimer in the
 
12
//    documentation and/or other materials provided with the distribution.
 
13
//  * Neither the name of any author or any participating institution
 
14
//    may be used to endorse or promote products derived from this software
 
15
//    without specific prior written permission.
 
16
// For a full list of authors, refer to the file AUTHORS.
 
17
// --------------------------------------------------------------------------
 
18
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 
19
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 
20
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 
21
// ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING
 
22
// INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
 
23
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
 
24
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
 
25
// OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
 
26
// WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
 
27
// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
 
28
// ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
22
29
//
23
30
// --------------------------------------------------------------------------
24
31
// $Maintainer: Andreas Bertsch $
54
61
 
55
62
using namespace std;
56
63
 
57
 
namespace OpenMS 
 
64
namespace OpenMS
58
65
{
59
 
        // HMMState implementation
60
 
        HMMState::HMMState()
61
 
                : hidden_(true)
62
 
        {
63
 
        }
64
 
 
65
 
        HMMState::HMMState(const String& name, bool hidden)
66
 
                : hidden_(hidden),
67
 
                        name_(name)
68
 
        {
69
 
                #ifdef STATE_DEBUG
70
 
                cerr << "State: " << name << ", hidden=" << hidden << endl;
71
 
                #endif
72
 
        }
73
 
 
74
 
        HMMState::HMMState(const HMMState& state)
75
 
                : hidden_(state.hidden_),
76
 
                        name_(state.name_)
77
 
        {
78
 
        }
79
 
        /*              pre_states_(state.pre_states_),
80
 
                        succ_states_(state.succ_states_)
81
 
        {
82
 
                
83
 
        }*/
84
 
        
85
 
        HMMState::~HMMState()
86
 
        {
87
 
        }
88
 
        
89
 
        void HMMState::setName(const String& name)
90
 
        {
91
 
                name_ = name;
92
 
        }
93
 
 
94
 
        const String& HMMState::getName() const
95
 
        {
96
 
                return name_;
97
 
        }
98
 
 
99
 
        HMMState& HMMState::operator = (const HMMState& state)
100
 
        {
101
 
                hidden_ = state.hidden_;
102
 
                name_ = state.name_;
103
 
                pre_states_.clear();
104
 
                succ_states_.clear();
105
 
                //pre_states_ = state.pre_states_;
106
 
                //succ_states_ = state.succ_states_;
107
 
                return *this;
108
 
        }
109
 
 
110
 
        void HMMState::addSuccessorState(HMMState* state)
111
 
        {
112
 
#ifdef SIMPLE_DEBUG
113
 
                cerr << "'" << state->getName() << "'" << endl;
114
 
#endif
115
 
                succ_states_.insert(state);
116
 
        }
117
 
 
118
 
        void HMMState::deleteSuccessorState(HMMState* state)
119
 
        {
120
 
#ifdef SIMPLE_DEBUG
121
 
    cerr << "'" << state->getName() << "'" << endl;
122
 
#endif
123
 
                succ_states_.erase(state);
124
 
        }
125
 
 
126
 
        void HMMState::addPredecessorState(HMMState* state)
127
 
        {
128
 
#ifdef SIMPLE_DEBUG
129
 
    cerr << "'" << state->getName() << "'" << endl;
130
 
#endif
131
 
                pre_states_.insert(state);
132
 
        }
133
 
 
134
 
        void HMMState::deletePredecessorState(HMMState* state)
135
 
        {
136
 
#ifdef SIMPLE_DEBUG
137
 
    cerr << "'" << state->getName() << "'" << endl;
138
 
#endif
139
 
                pre_states_.erase(state);
140
 
        }
141
 
 
142
 
        const set<HMMState*>& HMMState::getPredecessorStates() const
143
 
        {
144
 
                return pre_states_;
145
 
        }
146
 
 
147
 
        const set<HMMState*>& HMMState::getSuccessorStates() const
148
 
        {
149
 
                return succ_states_;
150
 
        }
151
 
 
152
 
        bool HMMState::isHidden() const
153
 
        {
154
 
                return hidden_;
155
 
        }
156
 
 
157
 
        void HMMState::setHidden(bool hidden)
158
 
        {
159
 
                hidden_ = hidden;
160
 
        }
161
 
 
162
 
        // The hidden markov model
163
 
        HiddenMarkovModel::HiddenMarkovModel() 
164
 
                : pseudo_counts_(0.0)
165
 
        {
166
 
        }
167
 
 
168
 
        HiddenMarkovModel::HiddenMarkovModel(const HiddenMarkovModel& hmm)
169
 
        {
170
 
                copy_(hmm);
171
 
        }
172
 
        
173
 
        HiddenMarkovModel::~HiddenMarkovModel()
174
 
        {
175
 
                clear();
176
 
        }
177
 
 
178
 
        HiddenMarkovModel& HiddenMarkovModel::operator = (const HiddenMarkovModel& hmm)
179
 
        {
180
 
                if (&hmm != this)
181
 
                {
182
 
                        clear();
183
 
                        copy_(hmm);
184
 
                }
185
 
                return *this;
186
 
        }
187
 
 
188
 
        HMMState* HiddenMarkovModel::getState(const String& name)
189
 
        {
190
 
                if (name_to_state_.find(name) == name_to_state_.end())
191
 
                {
192
 
                        throw Exception::ElementNotFound(__FILE__, __LINE__, __PRETTY_FUNCTION__, name);
193
 
                }
194
 
                return name_to_state_[name];
195
 
        }
196
 
 
197
 
        const HMMState* HiddenMarkovModel::getState(const String& name) const
198
 
        {
199
 
                if (name_to_state_.find(name) == name_to_state_.end())
200
 
                {
201
 
                        throw Exception::ElementNotFound(__FILE__, __LINE__, __PRETTY_FUNCTION__, name);
202
 
                }
203
 
                return name_to_state_.find(name)->second;
204
 
        }
205
 
        
206
 
        void HiddenMarkovModel::writeGraphMLFile(const String& filename)
207
 
        {
208
 
                set<HMMState*> states = states_;
209
 
                Map<HMMState*, vector<HMMState*> > transitions;
210
 
                
211
 
                ofstream out(filename.c_str());
 
66
  // HMMState implementation
 
67
  HMMState::HMMState() :
 
68
    hidden_(true)
 
69
  {
 
70
  }
 
71
 
 
72
  HMMState::HMMState(const String & name, bool hidden) :
 
73
    hidden_(hidden),
 
74
    name_(name)
 
75
  {
 
76
#ifdef STATE_DEBUG
 
77
    cerr << "State: " << name << ", hidden=" << hidden << endl;
 
78
#endif
 
79
  }
 
80
 
 
81
  HMMState::HMMState(const HMMState & state) :
 
82
    hidden_(state.hidden_),
 
83
    name_(state.name_)
 
84
  {
 
85
  }
 
86
 
 
87
  /*            pre_states_(state.pre_states_),
 
88
          succ_states_(state.succ_states_)
 
89
  {
 
90
 
 
91
  }*/
 
92
 
 
93
  HMMState::~HMMState()
 
94
  {
 
95
  }
 
96
 
 
97
  void HMMState::setName(const String & name)
 
98
  {
 
99
    name_ = name;
 
100
  }
 
101
 
 
102
  const String & HMMState::getName() const
 
103
  {
 
104
    return name_;
 
105
  }
 
106
 
 
107
  HMMState & HMMState::operator=(const HMMState & state)
 
108
  {
 
109
    hidden_ = state.hidden_;
 
110
    name_ = state.name_;
 
111
    pre_states_.clear();
 
112
    succ_states_.clear();
 
113
    //pre_states_ = state.pre_states_;
 
114
    //succ_states_ = state.succ_states_;
 
115
    return *this;
 
116
  }
 
117
 
 
118
  void HMMState::addSuccessorState(HMMState * state)
 
119
  {
 
120
#ifdef SIMPLE_DEBUG
 
121
    cerr << "'" << state->getName() << "'" << endl;
 
122
#endif
 
123
    succ_states_.insert(state);
 
124
  }
 
125
 
 
126
  void HMMState::deleteSuccessorState(HMMState * state)
 
127
  {
 
128
#ifdef SIMPLE_DEBUG
 
129
    cerr << "'" << state->getName() << "'" << endl;
 
130
#endif
 
131
    succ_states_.erase(state);
 
132
  }
 
133
 
 
134
  void HMMState::addPredecessorState(HMMState * state)
 
135
  {
 
136
#ifdef SIMPLE_DEBUG
 
137
    cerr << "'" << state->getName() << "'" << endl;
 
138
#endif
 
139
    pre_states_.insert(state);
 
140
  }
 
141
 
 
142
  void HMMState::deletePredecessorState(HMMState * state)
 
143
  {
 
144
#ifdef SIMPLE_DEBUG
 
145
    cerr << "'" << state->getName() << "'" << endl;
 
146
#endif
 
147
    pre_states_.erase(state);
 
148
  }
 
149
 
 
150
  const set<HMMState *> & HMMState::getPredecessorStates() const
 
151
  {
 
152
    return pre_states_;
 
153
  }
 
154
 
 
155
  const set<HMMState *> & HMMState::getSuccessorStates() const
 
156
  {
 
157
    return succ_states_;
 
158
  }
 
159
 
 
160
  bool HMMState::isHidden() const
 
161
  {
 
162
    return hidden_;
 
163
  }
 
164
 
 
165
  void HMMState::setHidden(bool hidden)
 
166
  {
 
167
    hidden_ = hidden;
 
168
  }
 
169
 
 
170
  // The hidden markov model
 
171
  HiddenMarkovModel::HiddenMarkovModel() :
 
172
    pseudo_counts_(0.0)
 
173
  {
 
174
  }
 
175
 
 
176
  HiddenMarkovModel::HiddenMarkovModel(const HiddenMarkovModel & hmm)
 
177
  {
 
178
    copy_(hmm);
 
179
  }
 
180
 
 
181
  HiddenMarkovModel::~HiddenMarkovModel()
 
182
  {
 
183
    clear();
 
184
  }
 
185
 
 
186
  HiddenMarkovModel & HiddenMarkovModel::operator=(const HiddenMarkovModel & hmm)
 
187
  {
 
188
    if (&hmm != this)
 
189
    {
 
190
      clear();
 
191
      copy_(hmm);
 
192
    }
 
193
    return *this;
 
194
  }
 
195
 
 
196
  HMMState * HiddenMarkovModel::getState(const String & name)
 
197
  {
 
198
    if (name_to_state_.find(name) == name_to_state_.end())
 
199
    {
 
200
      throw Exception::ElementNotFound(__FILE__, __LINE__, __PRETTY_FUNCTION__, name);
 
201
    }
 
202
    return name_to_state_[name];
 
203
  }
 
204
 
 
205
  const HMMState * HiddenMarkovModel::getState(const String & name) const
 
206
  {
 
207
    if (name_to_state_.find(name) == name_to_state_.end())
 
208
    {
 
209
      throw Exception::ElementNotFound(__FILE__, __LINE__, __PRETTY_FUNCTION__, name);
 
210
    }
 
211
    return name_to_state_.find(name)->second;
 
212
  }
 
213
 
 
214
  void HiddenMarkovModel::writeGraphMLFile(const String & filename)
 
215
  {
 
216
    set<HMMState *> states = states_;
 
217
    Map<HMMState *, vector<HMMState *> > transitions;
 
218
 
 
219
    ofstream out(filename.c_str());
212
220
    out << "<?xml version=\"1.0\" encoding=\"UTF-8\"?>" << endl;
213
221
 
214
 
                out << "<graphml xmlns=\"http://graphml.graphdrawing.org/xmlns/graphml\" xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\" " <<
215
 
         "xsi:schemaLocation=\"http://graphml.graphdrawing.org/xmlns/graphml http://www.yworks.com/xml/schema/graphml/1.0/ygraphml.xsd\" " <<
216
 
                                 "xmlns:y=\"http://www.yworks.com/xml/graphml\">" << endl;
 
222
    out << "<graphml xmlns=\"http://graphml.graphdrawing.org/xmlns/graphml\" xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\" " <<
 
223
    "xsi:schemaLocation=\"http://graphml.graphdrawing.org/xmlns/graphml http://www.yworks.com/xml/schema/graphml/1.0/ygraphml.xsd\" " <<
 
224
    "xmlns:y=\"http://www.yworks.com/xml/graphml\">" << endl;
217
225
 
218
 
                out << "<key id=\"d0\" for=\"node\" yfiles.type=\"nodegraphics\"/>" << endl;
219
 
        out << "<key id=\"d1\" for=\"edge\" yfiles.type=\"edgegraphics\"/>" << endl;
 
226
    out << "<key id=\"d0\" for=\"node\" yfiles.type=\"nodegraphics\"/>" << endl;
 
227
    out << "<key id=\"d1\" for=\"edge\" yfiles.type=\"edgegraphics\"/>" << endl;
220
228
    out << "  <graph id=\"G\" edgedefault=\"directed\">" << endl;
221
 
                for (set<HMMState*>::const_iterator it = states.begin(); it != states.end(); ++it)
222
 
                {
223
 
                        out << "    <node id=\"" << (*it)->getName() << "\">" << endl;
224
 
                        out << "      <data key=\"d0\">" << endl;
225
 
                        out << "        <y:ShapeNode>" << endl;
226
 
                        out << "          <y:NodeLabel>" << (*it)->getName() << "</y:NodeLabel>" << endl;
227
 
                        out << "        </y:ShapeNode>" << endl;
228
 
                        out << "      </data>" << endl;
229
 
                        out << "    </node>" << endl;
230
 
 
231
 
                        set<HMMState*> succ = (*it)->getSuccessorStates();
232
 
                        for (set<HMMState*>::const_iterator sit = succ.begin(); sit != succ.end(); ++sit)
233
 
                        {
234
 
                                transitions[*it].push_back(*sit);
235
 
                        }
236
 
                }
237
 
 
238
 
 
239
 
                
240
 
                for (Map<HMMState*, vector<HMMState*> >::const_iterator it = transitions.begin(); it != transitions.end(); ++it)
241
 
                {
242
 
                        for (vector<HMMState*>::const_iterator it1 = it->second.begin(); it1 != it->second.end(); ++it1)
243
 
                        {
244
 
                                out << "    <edge source=\"" << it->first->getName() << "\" target=\"" << (*it1)->getName() << "\" directed=\"true\">" << endl;
245
 
                                out << "      <data key=\"d1\">" << endl;
246
 
                                out << "        <y:PolyLineEdge>" << endl;
247
 
                                out << "          <y:EdgeLabel>" << getTransitionProbability_(it->first, *it1) << "</y:EdgeLabel>" << endl;
248
 
                                out << "        </y:PolyLineEdge>" << endl;
249
 
                                out << "      </data>" << endl;
250
 
                                out << "    </edge>" << endl;
251
 
                        }
252
 
                }
253
 
 
254
 
                out << "  </graph>" << endl;
255
 
                out << "</graphml>" << endl;
256
 
                
257
 
        }
258
 
 
259
 
        void HiddenMarkovModel::write(ostream& out) const
260
 
        {
261
 
                // write states
262
 
                for (set<HMMState*>::const_iterator it = states_.begin(); it != states_.end(); ++it)
263
 
                {
264
 
                        out << "State " << (*it)->getName();
265
 
 
266
 
                        if (!(*it)->isHidden())
267
 
                        {
268
 
                                out << " false";
269
 
                        }
270
 
                        out << endl;
271
 
                }
272
 
 
273
 
                // write transitions
274
 
                for (Map<HMMState*, Map<HMMState*, DoubleReal> >::const_iterator it1 = trans_.begin(); it1 != trans_.end(); ++it1)
275
 
                {
276
 
                        for (Map<HMMState*, DoubleReal>::const_iterator it2 = it1->second.begin(); it2 != it1->second.end(); ++it2)
277
 
                        {
278
 
                                out << "Transition " << it1->first->getName() << " " << it2->first->getName() << " " << it2->second << endl;
279
 
                        }
280
 
                }
281
 
        
282
 
                // write synonyms
283
 
                /*
284
 
                for (Map<String, Map<String, pair<String, String> > >::const_iterator it = synonym_trans_names_.begin(); it != synonym_trans_names_.end(); ++it)
285
 
                {
286
 
                        for (Map<String, pair<String, String> >::const_iterator it2 = it->second.begin(); it2 != it->second.end(); ++it2)
287
 
                        {
288
 
                                out << "Synonym " << it->first << " " << it2->first << " " << it2->second.first << " " << it2->second.second << endl;
289
 
                        }
290
 
                }*/
291
 
 
292
 
    for (Map<HMMState*, Map<HMMState*, std::pair<HMMState*, HMMState*> > >::const_iterator it1 = synonym_trans_.begin(); it1 != synonym_trans_.end(); ++it1)
293
 
    {
294
 
      for (Map<HMMState*, std::pair<HMMState*, HMMState*> >::const_iterator it2 = it1->second.begin(); it2 != it1->second.end(); ++it2)
295
 
      {
296
 
                                //cerr << "Synonym " << it1->first->getName() << " " << it2->first->getName() << " " << it2->second.first->getName() << " " << it2->second.second->getName() << endl;
 
229
    for (set<HMMState *>::const_iterator it = states.begin(); it != states.end(); ++it)
 
230
    {
 
231
      out << "    <node id=\"" << (*it)->getName() << "\">" << endl;
 
232
      out << "      <data key=\"d0\">" << endl;
 
233
      out << "        <y:ShapeNode>" << endl;
 
234
      out << "          <y:NodeLabel>" << (*it)->getName() << "</y:NodeLabel>" << endl;
 
235
      out << "        </y:ShapeNode>" << endl;
 
236
      out << "      </data>" << endl;
 
237
      out << "    </node>" << endl;
 
238
 
 
239
      set<HMMState *> succ = (*it)->getSuccessorStates();
 
240
      for (set<HMMState *>::const_iterator sit = succ.begin(); sit != succ.end(); ++sit)
 
241
      {
 
242
        transitions[*it].push_back(*sit);
 
243
      }
 
244
    }
 
245
 
 
246
 
 
247
 
 
248
    for (Map<HMMState *, vector<HMMState *> >::const_iterator it = transitions.begin(); it != transitions.end(); ++it)
 
249
    {
 
250
      for (vector<HMMState *>::const_iterator it1 = it->second.begin(); it1 != it->second.end(); ++it1)
 
251
      {
 
252
        out << "    <edge source=\"" << it->first->getName() << "\" target=\"" << (*it1)->getName() << "\" directed=\"true\">" << endl;
 
253
        out << "      <data key=\"d1\">" << endl;
 
254
        out << "        <y:PolyLineEdge>" << endl;
 
255
        out << "          <y:EdgeLabel>" << getTransitionProbability_(it->first, *it1) << "</y:EdgeLabel>" << endl;
 
256
        out << "        </y:PolyLineEdge>" << endl;
 
257
        out << "      </data>" << endl;
 
258
        out << "    </edge>" << endl;
 
259
      }
 
260
    }
 
261
 
 
262
    out << "  </graph>" << endl;
 
263
    out << "</graphml>" << endl;
 
264
 
 
265
  }
 
266
 
 
267
  void HiddenMarkovModel::write(ostream & out) const
 
268
  {
 
269
    // write states
 
270
    for (set<HMMState *>::const_iterator it = states_.begin(); it != states_.end(); ++it)
 
271
    {
 
272
      out << "State " << (*it)->getName();
 
273
 
 
274
      if (!(*it)->isHidden())
 
275
      {
 
276
        out << " false";
 
277
      }
 
278
      out << endl;
 
279
    }
 
280
 
 
281
    // write transitions
 
282
    for (Map<HMMState *, Map<HMMState *, DoubleReal> >::const_iterator it1 = trans_.begin(); it1 != trans_.end(); ++it1)
 
283
    {
 
284
      for (Map<HMMState *, DoubleReal>::const_iterator it2 = it1->second.begin(); it2 != it1->second.end(); ++it2)
 
285
      {
 
286
        out << "Transition " << it1->first->getName() << " " << it2->first->getName() << " " << it2->second << endl;
 
287
      }
 
288
    }
 
289
 
 
290
    // write synonyms
 
291
    /*
 
292
    for (Map<String, Map<String, pair<String, String> > >::const_iterator it = synonym_trans_names_.begin(); it != synonym_trans_names_.end(); ++it)
 
293
    {
 
294
        for (Map<String, pair<String, String> >::const_iterator it2 = it->second.begin(); it2 != it->second.end(); ++it2)
 
295
        {
 
296
            out << "Synonym " << it->first << " " << it2->first << " " << it2->second.first << " " << it2->second.second << endl;
 
297
        }
 
298
    }*/
 
299
 
 
300
    for (Map<HMMState *, Map<HMMState *, std::pair<HMMState *, HMMState *> > >::const_iterator it1 = synonym_trans_.begin(); it1 != synonym_trans_.end(); ++it1)
 
301
    {
 
302
      for (Map<HMMState *, std::pair<HMMState *, HMMState *> >::const_iterator it2 = it1->second.begin(); it2 != it1->second.end(); ++it2)
 
303
      {
 
304
        //cerr << "Synonym " << it1->first->getName() << " " << it2->first->getName() << " " << it2->second.first->getName() << " " << it2->second.second->getName() << endl;
297
305
        out << "Synonym " << it1->first->getName() << " " << it2->first->getName() << " " << it2->second.first->getName() << " " << it2->second.second->getName() << endl;
298
306
      }
299
307
    }
300
308
 
301
 
                return;
302
 
        }
 
309
    return;
 
310
  }
303
311
 
304
 
        void HiddenMarkovModel::clear()
305
 
        {
306
 
                for (set<HMMState*>::const_iterator it = states_.begin(); it != states_.end(); ++it)
 
312
  void HiddenMarkovModel::clear()
 
313
  {
 
314
    for (set<HMMState *>::const_iterator it = states_.begin(); it != states_.end(); ++it)
307
315
    {
308
316
      delete *it;
309
317
    }
310
318
    trans_.clear();
311
319
    count_trans_.clear();
312
 
                train_count_trans_all_.clear();
 
320
    train_count_trans_all_.clear();
313
321
    forward_.clear();
314
322
    backward_.clear();
315
323
    name_to_state_.clear();
318
326
    states_.clear();
319
327
    trained_trans_.clear();
320
328
    synonym_trans_.clear();
321
 
    synonym_trans_names_.clear();       
322
 
                var_modifications_ = StringList::create("");
323
 
                return;
324
 
        }
325
 
        
326
 
        Size HiddenMarkovModel::getNumberOfStates() const
327
 
        {
328
 
                return states_.size();
329
 
        }
330
 
 
331
 
        void HiddenMarkovModel::addNewState(HMMState* s)
332
 
        {
333
 
                states_.insert(s);
334
 
                if (name_to_state_.find(s->getName()) == name_to_state_.end())
335
 
                {
336
 
                        name_to_state_[s->getName()] = s;
337
 
                }
338
 
                else
339
 
                {
340
 
                        cerr << "HiddenMarkovModel: state name '" << s->getName() << "' (" << s << ") already used!" << endl;
341
 
                }
342
 
        }
343
 
 
344
 
        void HiddenMarkovModel::addNewState(const String& name)
345
 
        {
346
 
                HMMState* s = new HMMState(name);
347
 
                states_.insert(s);
348
 
                if (name_to_state_.find(name) == name_to_state_.end())
349
 
                {
350
 
                        name_to_state_[name] = s;
351
 
                }
352
 
                else
353
 
                {
354
 
                        cerr << "HiddenMarkovModel: state name '" << name << "' (" << name_to_state_[name] << ") already used!" << endl;
355
 
                }
356
 
                return;
357
 
        }
358
 
 
359
 
        void HiddenMarkovModel::addSynonymTransition(const String& name1, const String& name2, const String& synonym1, const String& synonym2)
360
 
        {
 
329
    synonym_trans_names_.clear();
 
330
    var_modifications_ = StringList::create("");
 
331
    return;
 
332
  }
 
333
 
 
334
  Size HiddenMarkovModel::getNumberOfStates() const
 
335
  {
 
336
    return states_.size();
 
337
  }
 
338
 
 
339
  void HiddenMarkovModel::addNewState(HMMState * s)
 
340
  {
 
341
    states_.insert(s);
 
342
    if (name_to_state_.find(s->getName()) == name_to_state_.end())
 
343
    {
 
344
      name_to_state_[s->getName()] = s;
 
345
    }
 
346
    else
 
347
    {
 
348
      cerr << "HiddenMarkovModel: state name '" << s->getName() << "' (" << s << ") already used!" << endl;
 
349
    }
 
350
  }
 
351
 
 
352
  void HiddenMarkovModel::addNewState(const String & name)
 
353
  {
 
354
    HMMState * s = new HMMState(name);
 
355
    states_.insert(s);
 
356
    if (name_to_state_.find(name) == name_to_state_.end())
 
357
    {
 
358
      name_to_state_[name] = s;
 
359
    }
 
360
    else
 
361
    {
 
362
      cerr << "HiddenMarkovModel: state name '" << name << "' (" << name_to_state_[name] << ") already used!" << endl;
 
363
    }
 
364
    return;
 
365
  }
 
366
 
 
367
  void HiddenMarkovModel::addSynonymTransition(const String & name1, const String & name2, const String & synonym1, const String & synonym2)
 
368
  {
361
369
#ifdef SIMPLE_DEBUG2
362
 
                cerr << "addSynonymTransition: '" << name1 << "' -> '" << name2 << "' :: '" << synonym1 << "' -> '" << synonym2 << "'" << endl;
 
370
    cerr << "addSynonymTransition: '" << name1 << "' -> '" << name2 << "' :: '" << synonym1 << "' -> '" << synonym2 << "'" << endl;
363
371
#endif
364
 
                if (name_to_state_.find(name1) == name_to_state_.end())
365
 
                {
366
 
                        cerr << "state '" << name1 << "' unknown" << endl;
367
 
                }
 
372
    if (name_to_state_.find(name1) == name_to_state_.end())
 
373
    {
 
374
      cerr << "state '" << name1 << "' unknown" << endl;
 
375
    }
368
376
    if (name_to_state_.find(name2) == name_to_state_.end())
369
377
    {
370
378
      cerr << "state '" << name2 << "' unknown" << endl;
377
385
    {
378
386
      cerr << "state '" << synonym2 << "' unknown" << endl;
379
387
    }
380
 
                synonym_trans_names_[synonym1][synonym2] = make_pair(name1, name2);
 
388
    synonym_trans_names_[synonym1][synonym2] = make_pair(name1, name2);
381
389
 
382
 
                synonym_trans_[name_to_state_[synonym1]][name_to_state_[synonym2]] = make_pair(name_to_state_[name1], name_to_state_[name2]);
 
390
    synonym_trans_[name_to_state_[synonym1]][name_to_state_[synonym2]] = make_pair(name_to_state_[name1], name_to_state_[name2]);
383
391
 
384
392
/*
385
 
                for (Map<String, Map<String, std::pair<String, String> > >::const_iterator it = synonym_trans_names_.begin(); it != synonym_trans_names_.end(); ++it)
 
393
        for (Map<String, Map<String, std::pair<String, String> > >::const_iterator it = synonym_trans_names_.begin(); it != synonym_trans_names_.end(); ++it)
386
394
    {
387
395
      for (Map<String, pair<String, String> >::const_iterator it2 = it->second.begin(); it2 != it->second.end(); ++it2)
388
396
      {
390
398
            make_pair(name_to_state_[it2->second.first], name_to_state_[it2->second.second]);
391
399
      }
392
400
    }
393
 
        */      
394
 
        }
395
 
 
396
 
        void HiddenMarkovModel::setTransitionProbability_(HMMState * s1, HMMState * s2, DoubleReal trans_prob)
397
 
        {
398
 
                trans_[s1][s2] = trans_prob;
399
 
                s1->addSuccessorState(s2);
400
 
                s2->addPredecessorState(s1);
401
 
                enabled_trans_[s1].insert(s2);
402
 
                training_steps_count_[s1][s2] = 0;
403
 
        }
404
 
 
405
 
        void HiddenMarkovModel::setTransitionProbability(const String& s1, const String& s2, DoubleReal trans_prob)
406
 
        {
 
401
    */
 
402
  }
 
403
 
 
404
  void HiddenMarkovModel::setTransitionProbability_(HMMState * s1, HMMState * s2, DoubleReal trans_prob)
 
405
  {
 
406
    trans_[s1][s2] = trans_prob;
 
407
    s1->addSuccessorState(s2);
 
408
    s2->addPredecessorState(s1);
 
409
    enabled_trans_[s1].insert(s2);
 
410
    training_steps_count_[s1][s2] = 0;
 
411
  }
 
412
 
 
413
  void HiddenMarkovModel::setTransitionProbability(const String & s1, const String & s2, DoubleReal trans_prob)
 
414
  {
407
415
    OPENMS_PRECONDITION(
408
 
                name_to_state_.find(s1) != name_to_state_.end() && name_to_state_.find(s2) != name_to_state_.end(),
409
 
                String("HiddenMarkovModel::setTransitionProbability(" + String(s1) + ", " + String(s2) + ", " +  String(trans_prob) + "), no suchstate!").c_str());
 
416
      name_to_state_.find(s1) != name_to_state_.end() && name_to_state_.find(s2) != name_to_state_.end(),
 
417
      String("HiddenMarkovModel::setTransitionProbability(" + String(s1) + ", " + String(s2) + ", " +  String(trans_prob) + "), no suchstate!").c_str());
410
418
#ifdef SIMPLE_DEBUG2
411
419
    cerr << "setTransitionProbability: '" << s1 << "' -> '" << s2 << "'" << " " << trans_prob << endl;
412
420
#endif
413
421
 
414
 
                trans_[name_to_state_[s1]][name_to_state_[s2]] = trans_prob;
415
 
                name_to_state_[s1]->addSuccessorState(name_to_state_[s2]);
416
 
                name_to_state_[s2]->addPredecessorState(name_to_state_[s1]);
417
 
                enabled_trans_[name_to_state_[s1]].insert(name_to_state_[s2]);
418
 
                training_steps_count_[name_to_state_[s1]][name_to_state_[s2]] = 0;
419
 
        }
420
 
 
421
 
        DoubleReal HiddenMarkovModel::getTransitionProbability(const String& s1, const String& s2) const
422
 
        {
423
 
                if (name_to_state_.find(s1) == name_to_state_.end())
424
 
                {
425
 
                        throw Exception::ElementNotFound(__FILE__, __LINE__, __PRETTY_FUNCTION__, s1);
426
 
                }
427
 
                HMMState* state1 = name_to_state_[s1];
428
 
                if (name_to_state_.find(s2) == name_to_state_.end())
429
 
                {
430
 
                        throw Exception::ElementNotFound(__FILE__, __LINE__, __PRETTY_FUNCTION__, s2);
431
 
                }
432
 
                HMMState* state2 = name_to_state_[s2];
433
 
                return getTransitionProbability_(state1, state2);
434
 
        }
435
 
        
436
 
        DoubleReal HiddenMarkovModel::getTransitionProbability_(HMMState* s1, HMMState* s2) const
437
 
        {
438
 
                HMMState* state1 = s1;
439
 
                HMMState* state2 = s2;
440
 
 
441
 
#ifdef SIMPLE_DEBUG2
442
 
                cerr << "getTransitionProbability(" << s1->getName() << ", " << s2->getName() << ")" << endl;
443
 
#endif
444
 
                
445
 
                // check if transition is a synonym transition
446
 
                if (synonym_trans_.find(state1) != synonym_trans_.end() && 
447
 
                                synonym_trans_.find(state1)->second.find(state2) != synonym_trans_.find(state1)->second.end())
448
 
                {
449
 
                        HMMState* tmp = synonym_trans_.find(state1)->second.find(state2)->second.first;
450
 
                        state2 = synonym_trans_.find(state1)->second.find(state2)->second.second; //name_to_state_[p.second];
451
 
                        state1 = tmp;
452
 
                }
453
 
                
454
 
#ifdef SIMPLE_DEBUG2
455
 
                cerr << "getTransitionProbability: " << state1->getName() << " " << state2->getName() << endl;
456
 
#endif
457
 
                
458
 
                // get the transition prob 
459
 
                if (trans_.find(state1) != trans_.end() && trans_.find(state1)->second.find(state2) != trans_.find(state1)->second.end())
460
 
                {
461
 
                /*
462
 
                        // TODO !!!! ???? path length correction
463
 
                        if (state2->getName().hasSubstring("next"))
464
 
                        {
465
 
                                return sqrt(trans_[state1][state2]);
466
 
                        }
467
 
                        else
468
 
                        {*/
469
 
                                return trans_.find(state1)->second.find(state2)->second;
470
 
                                /*
471
 
                        }
472
 
                        */
473
 
                }
474
 
                else
475
 
                {
476
 
                        return 0;
477
 
                }
478
 
        }
479
 
 
480
 
        void HiddenMarkovModel::train()
481
 
        {
482
 
#ifdef HIDDEN_MARKOV_MODEL_DEBUG
483
 
                cerr << "HiddenMarkovModel::train()" << endl;
484
 
#endif
485
 
                trained_trans_.clear();
486
 
 
487
 
                #ifdef HIDDEN_MARKOV_MODEL_DEBUG
488
 
                cerr << "calc forward part" << endl;
489
 
                #endif
490
 
                
491
 
                calculateForwardPart_();
492
 
 
493
 
                #ifdef HIDDEN_MARKOV_MODEL_DEBUG
494
 
                cerr << "calc backward part" << endl;
495
 
                #endif
496
 
                
497
 
                calculateBackwardPart_();
498
 
        
499
 
 
500
 
                #ifdef HIDDEN_MARKOV_MODEL_DEBUG
501
 
                cerr << "calc px" << endl;
502
 
                #endif
503
 
        
504
 
                // calc p_x from forward part
505
 
                DoubleReal px(0);
506
 
 
507
 
                for (Map<HMMState*, DoubleReal>::const_iterator it1 = train_emission_prob_.begin(); it1 != train_emission_prob_.end(); ++it1)
508
 
                {
509
 
                        for (set<HMMState*>::const_iterator it2 = it1->first->getPredecessorStates().begin(); it2 != it1->first->getPredecessorStates().end(); ++it2)
510
 
                        {
511
 
                                px += getForwardVariable_(*it2);
512
 
                        }
513
 
                }
514
 
 
515
 
                DoubleReal num_px(0);
516
 
                if (px != 0)
517
 
                {
518
 
                        num_px = 1.0/px;
519
 
                }
520
 
 
521
 
#ifdef HIDDEN_MARKOV_MODEL_DEBUG
522
 
                cerr << "px=" << num_px << endl;
523
 
                cerr << "add contributions to count_trans" << endl;
524
 
#endif
525
 
 
526
 
                // add contributions to count_trans_
527
 
                for (set<pair<HMMState*, HMMState*> >::const_iterator it = trained_trans_.begin(); it != trained_trans_.end(); ++it)
528
 
                {
529
 
                        DoubleReal tmp(0);
530
 
                        tmp = num_px * getForwardVariable_(it->first) * getBackwardVariable_(it->second) * getTransitionProbability_(it->first, it->second);
531
 
                        tmp += pseudo_counts_;
532
 
                        HMMState* s1 = it->first;
533
 
                        HMMState* s2 = it->second;
534
 
 
535
 
                        if (synonym_trans_.find(s1) != synonym_trans_.end() && synonym_trans_[s1].find(s2) != synonym_trans_[s1].end())
536
 
                        {
537
 
                                HMMState* tmp = synonym_trans_[s1][s2].first;
538
 
                                s2 = synonym_trans_[s1][s2].second;
539
 
                                s1 = tmp;
540
 
                        }
541
 
                        
542
 
                        train_count_trans_all_[s1][s2].push_back(tmp);
543
 
                        if (count_trans_.find(s1) != count_trans_.end() && count_trans_[s1].find(s2) != count_trans_[s1].end())
544
 
                        {
545
 
                                count_trans_[s1][s2] += tmp;
546
 
                        }
547
 
                        else
548
 
                        {
549
 
                                count_trans_[s1][s2] = tmp;
550
 
                        }
551
 
                        training_steps_count_[s1][s2]++;
552
 
                }
553
 
        }
554
 
 
555
 
        void HiddenMarkovModel::calculateForwardPart_()
556
 
        {
557
 
                forward_.clear();
558
 
                set<HMMState*> succ;
559
 
                for (Map<HMMState*, DoubleReal>::iterator it = init_prob_.begin(); it != init_prob_.end(); ++it)
560
 
                {
561
 
                        //cerr << it->first << " " << it->second << endl;
562
 
                        //cerr << it->first->getName() << endl;
563
 
                        forward_[it->first] = it->second;
564
 
                }
565
 
                
566
 
                for (Map<HMMState*, DoubleReal>::iterator it = init_prob_.begin(); it != init_prob_.end(); ++it)
567
 
                {
568
 
                        succ.insert(it->first->getSuccessorStates().begin(), it->first->getSuccessorStates().end());
569
 
                
570
 
      while ( !succ.empty() )
571
 
                        {
572
 
                                set<HMMState*> succ_new;
573
 
                                for (set<HMMState*>::const_iterator it = succ.begin(); it != succ.end(); ++it)
574
 
                                {
575
 
                                        set<HMMState*> pre = (*it)->getPredecessorStates();
576
 
                                        DoubleReal sum(0);
577
 
                                        for (set<HMMState*>::const_iterator it2 = pre.begin(); it2 != pre.end(); ++it2)
578
 
                                        {
579
 
                                                #ifdef HIDDEN_MARKOV_MODEL_DEBUG
580
 
                                                cerr << "\tadding from " << (*it2)->getName() << " " << getForwardVariable_(*it2) << " * " << getTransitionProbability(*it2, *it) << endl;
581
 
                                                #endif
582
 
                                                sum += getForwardVariable_(*it2) * getTransitionProbability_(*it2, *it);
583
 
                                                trained_trans_.insert(make_pair(*it2, *it));
584
 
                                        }
585
 
                                        forward_[*it] = sum;
586
 
                                        succ_new.insert((*it)->getSuccessorStates().begin(), (*it)->getSuccessorStates().end());
587
 
                                        #ifdef HIDDEN_MARKOV_MODEL_DEBUG
588
 
                                        cerr << "f: " << (*it)->getName() << "\t" << sum << endl;
589
 
                                        #endif
590
 
                                }
591
 
                                succ = succ_new;
592
 
                        }
593
 
                }
594
 
        }
595
 
 
596
 
        void HiddenMarkovModel::calculateBackwardPart_()
597
 
        {
598
 
                backward_.clear();
599
 
                set<HMMState*> pre;
600
 
                for (Map<HMMState*, DoubleReal>::iterator it = train_emission_prob_.begin(); it != train_emission_prob_.end(); ++it)
601
 
                {
602
 
                        backward_[it->first] = it->second;
603
 
                }
604
 
 
605
 
                for (Map<HMMState*, DoubleReal>::iterator it = train_emission_prob_.begin(); it != train_emission_prob_.end(); ++it)
606
 
                {
607
 
                        #ifdef HIDDEN_MARKOV_MODEL_DEBUG
608
 
                        cerr << "b:" << it->first << " " <<  it->first->getName() << " " << it->second << endl;
609
 
                        #endif
610
 
                        pre.insert(it->first->getPredecessorStates().begin(), it->first->getPredecessorStates().end());
611
 
                
612
 
      while ( !pre.empty())
613
 
                        {
614
 
                                set <HMMState*> pre_new;
615
 
                                for (set<HMMState*>::const_iterator it = pre.begin(); it != pre.end(); ++it)
616
 
                                {
617
 
                                        set<HMMState*> succ = (*it)->getSuccessorStates();
618
 
                                        DoubleReal sum(0);
619
 
                                        for (set<HMMState*>::const_iterator it2 = succ.begin(); it2 != succ.end(); ++it2)
620
 
                                        {
621
 
                                                #ifdef HIDDEN_MARKOV_MODEL_DEBUG
622
 
                                                cerr << "\tadding from " << (*it2)->getName() << " " << getBackwardVariable_(*it2) << " * " << getTransitionProbability(*it, *it2) << endl;
623
 
                                                #endif
624
 
                                                sum += getBackwardVariable_(*it2) * getTransitionProbability_(*it, *it2);
625
 
                                                trained_trans_.insert(make_pair(*it, *it2));
626
 
                                        }
627
 
                                        backward_[*it] = sum;
628
 
                                        pre_new.insert((*it)->getPredecessorStates().begin(), (*it)->getPredecessorStates().end());
629
 
                                        #ifdef HIDDEN_MARKOV_MODEL_DEBUG
630
 
                                        cerr << "b: " << (*it)->getName() << "\t" << sum << endl;
631
 
                                        #endif
632
 
                                }
633
 
                                pre = pre_new;
634
 
                        }
635
 
                }
636
 
        }
637
 
 
638
 
        DoubleReal HiddenMarkovModel::getForwardVariable_(HMMState* state)
639
 
        {
640
 
                return forward_.find(state) != forward_.end() ? forward_[state] : 0;
641
 
        }
642
 
 
643
 
        DoubleReal HiddenMarkovModel::getBackwardVariable_(HMMState* state)
644
 
        {
645
 
                return backward_.find(state) != backward_.end() ? backward_[state] : 0;
646
 
        }
647
 
 
648
 
        void HiddenMarkovModel::evaluate()
649
 
        {
650
 
                for (Map<HMMState*, Map<HMMState*, DoubleReal> >::const_iterator it1 = count_trans_.begin(); it1 != count_trans_.end(); ++it1)
651
 
                {
652
 
#ifdef EVALUATE_DEBUG
653
 
                        cerr <<  it1->first->getName() << endl;
654
 
#endif
655
 
                        DoubleReal sum(0);
656
 
                        for (Map<HMMState*, DoubleReal>::const_iterator it2 = it1->second.begin(); it2 != it1->second.end(); ++it2)
657
 
                        {
658
 
                                if (count_trans_.find(it1->first) != count_trans_.end() && 
659
 
                                                count_trans_[it1->first].find(it2->first) != count_trans_[it1->first].end())
660
 
                                {
661
 
                                        sum += count_trans_[it1->first][it2->first];
662
 
#ifdef EVALUATE_DEBUG
663
 
                                        cerr << it1->first->getName() << " " << it2->first->getName() << " ";
664
 
                                        
665
 
                                        //<< count_trans_[it1->first][it2->first] << endl;
666
 
                                        for (vector<DoubleReal>::const_iterator it = train_count_trans_all_[it1->first][it2->first].begin(); it != train_count_trans_all_[it1->first][it2->first].end(); ++it)
667
 
                                        {
668
 
                                                cerr << *it << " ";
669
 
                                        }
670
 
                                        vector<DoubleReal> data = train_count_trans_all_[it1->first][it2->first];
671
 
                                        std::sort(data.begin(),data.end());
672
 
                      DoubleReal mean = gsl_stats_mean(&data.front(),1,data.size());
673
 
                      DoubleReal variance = gsl_stats_variance_m(&data.front(),1,data.size(),mean);
674
 
                                        cerr << "mean=" << mean << ", variance=" << variance << endl;
675
 
#endif
676
 
                                }
677
 
                        }
678
 
 
679
 
                        if (sum != 0)
680
 
                        {
681
 
                                for (Map<HMMState*, DoubleReal>::const_iterator it2 = it1->second.begin(); it2 != it1->second.end(); ++it2)
682
 
                                {
683
 
                                        if (count_trans_.find(it1->first) != count_trans_.end() && 
684
 
                                                        count_trans_[it1->first].find(it2->first) != count_trans_[it1->first].end())
685
 
                                        {
686
 
                                                trans_[it1->first][it2->first] = count_trans_[it1->first][it2->first] / sum;
687
 
                                        }
688
 
                                }
689
 
                        }
690
 
                }
691
 
        }
692
 
 
693
 
        void HiddenMarkovModel::setInitialTransitionProbability(const String& state, DoubleReal prob)
694
 
        {
695
 
                OPENMS_PRECONDITION(name_to_state_.find(state) != name_to_state_.end(), String("HiddenMarkovModel::setInitialTransitionProbability(" + state + ", " + String(prob) + "), no suchstate!").c_str());
696
 
                //cerr << state << " " << prob << endl;
697
 
                init_prob_[name_to_state_[state]] = prob;
698
 
        }
699
 
 
700
 
        void HiddenMarkovModel::clearInitialTransitionProbabilities()
701
 
        {
702
 
                init_prob_.clear();
703
 
        }
704
 
 
705
 
        void HiddenMarkovModel::setTrainingEmissionProbability(const String& state, DoubleReal prob)
706
 
        {
707
 
#ifdef SIMPLE_DEBUG2
708
 
                cerr << "setTrainingEmissionProbability(" << state << "(" << name_to_state_[state] << "), " << prob << ")" << endl;
709
 
#endif
710
 
 
711
 
                OPENMS_PRECONDITION(name_to_state_.find(state) != name_to_state_.end(), String("HiddenMarkovModel::setTrainingEmissionProbability(" + state + ", " + String(prob) + "), no such state!").c_str());
712
 
                train_emission_prob_[name_to_state_[state]] = prob;
713
 
        }
714
 
 
715
 
        void HiddenMarkovModel::clearTrainingEmissionProbabilities()
716
 
        {
717
 
                train_emission_prob_.clear();
718
 
        }
719
 
 
720
 
        void HiddenMarkovModel::enableTransition(const String& s1, const String& s2)
721
 
        {
722
 
#ifdef SIMPLE_DEBUG2
723
 
                cerr << "enableTransition: '" << s1 << "' -> '" << s2 << "'" << endl;
724
 
#endif
725
 
                enableTransition_(name_to_state_[s1], name_to_state_[s2]);
726
 
        }
727
 
 
728
 
        void HiddenMarkovModel::enableTransition_(HMMState* s1, HMMState* s2)
729
 
        {
730
 
                s1->addSuccessorState(s2);
731
 
                s2->addPredecessorState(s1);
732
 
                enabled_trans_[s1].insert(s2);
733
 
        }
734
 
 
735
 
        void HiddenMarkovModel::disableTransition(const String& s1, const String& s2)
736
 
        {
737
 
                disableTransition_(name_to_state_[s1], name_to_state_[s2]);
738
 
        }
739
 
        
740
 
        void HiddenMarkovModel::disableTransition_(HMMState* s1, HMMState* s2)
741
 
        {
742
 
                s1->deleteSuccessorState(s2);
743
 
                s2->deletePredecessorState(s1);
744
 
                enabled_trans_[s1].erase(s2);
745
 
        }
746
 
 
747
 
        void HiddenMarkovModel::disableTransitions()
748
 
        {
749
 
                for (Map<HMMState*, set<HMMState*> >::const_iterator it = enabled_trans_.begin(); it != enabled_trans_.end(); ++it)
750
 
                {
751
 
                        for (set<HMMState*>::const_iterator it2 = it->second.begin(); it2 != it->second.end(); ++it2)
752
 
                        {
753
 
                                it->first->deleteSuccessorState(*it2);
754
 
                                (*it2)->deletePredecessorState(it->first);
755
 
                        }
756
 
                }
757
 
                enabled_trans_.clear();
758
 
        }
759
 
 
760
 
        void HiddenMarkovModel::calculateEmissionProbabilities(Map<HMMState*, DoubleReal>& emission_probs)
761
 
        {
762
 
                Map<HMMState*, DoubleReal> states = init_prob_;
763
 
 
764
 
    while ( !states.empty() )
765
 
                {
766
 
                        Map<HMMState*, DoubleReal> tmp = states;
767
 
                        for (Map<HMMState*, DoubleReal>::const_iterator it = tmp.begin(); it != tmp.end(); ++it)
768
 
                        {
769
 
                                for (set<HMMState*>::const_iterator it2 = it->first->getSuccessorStates().begin(); it2 != it->first->getSuccessorStates().end(); ++it2)
770
 
                                {
771
 
                                        //cerr << "->" << (*it2)->getName() << "=(from " << it->first->getName() << "): "; // << states[it->first] << " * " << getTransitionProbability_(it->first, *it2) << " ";
772
 
                                        if (states.find(*it2) != states.end())
773
 
                                        {
774
 
                                                //cerr << " += " << states[it->first] * getTransitionProbability_(it->first, *it2);
775
 
                                                states[*it2] += states[it->first] * getTransitionProbability_(it->first, *it2);
776
 
                                        }
777
 
                                        else
778
 
                                        {
779
 
                                                //cerr << " = " << states[it->first] * getTransitionProbability_(it->first, *it2);
780
 
                                                states[*it2] = states[it->first] * getTransitionProbability_(it->first, *it2);
781
 
                                        }
782
 
                                        if (!(*it2)->isHidden())
783
 
                                        {
784
 
                                                if (emission_probs.find(*it2) != emission_probs.end())
785
 
                                                {
786
 
                                                        //cerr << " emission: += " << states[it->first] * getTransitionProbability_(it->first, *it2);
787
 
                                                        emission_probs[*it2] += states[it->first] * getTransitionProbability_(it->first, *it2);
788
 
                                                }
789
 
                                                else
790
 
                                                {       
791
 
                                                        //cerr << " emission: += " << states[it->first] * getTransitionProbability_(it->first, *it2);
792
 
                                                        emission_probs[*it2] = states[it->first] * getTransitionProbability_(it->first, *it2);
793
 
                                                }
794
 
                                        }
795
 
                                        //cerr << endl;
796
 
                                }
797
 
                                states.erase(it->first);
798
 
                        }
799
 
                }
800
 
                
801
 
                return;
802
 
        }
803
 
 
804
 
        void HiddenMarkovModel::dump()
805
 
        {
806
 
                cerr << "dump of transitions: " << endl;
807
 
                for (Map<HMMState*, Map<HMMState*, DoubleReal> >::const_iterator it = trans_.begin(); it != trans_.end(); ++it)
808
 
                {
809
 
                        for (Map<HMMState*, DoubleReal>::const_iterator it1 = it->second.begin(); it1 != it->second.end(); ++it1)
810
 
                        {
811
 
                                cout << it->first->getName() << " -> " << it1->first->getName() << " " << it1->second << " " << training_steps_count_[it->first][it1->first] << ": ";
812
 
                                vector<DoubleReal> all_trans = train_count_trans_all_[it->first][it1->first];
813
 
                                
814
 
        if ( !all_trans.empty() )
815
 
                                {
816
 
                                        DoubleReal sum = accumulate(all_trans.begin(), all_trans.end(), 0.0);
817
 
                                        DoubleReal avg(sum/DoubleReal(all_trans.size()));
818
 
                                        DoubleReal rsd(0);
819
 
                                        for (Size i = 0; i != all_trans.size(); ++i)
820
 
                                        {
821
 
                                                cout << all_trans[i] << " ";
822
 
                                                rsd += abs(all_trans[i] - avg);
823
 
                                        }
824
 
                                        cout << "rsd=" << rsd / DoubleReal(all_trans.size()) / avg;
825
 
                                        cout << ", avg=" << avg;
826
 
                                }
827
 
                                
828
 
                                cout << endl;
829
 
                        }
830
 
                }
831
 
                
832
 
                cerr << "dump completed" << endl;
833
 
        }
834
 
 
835
 
        void HiddenMarkovModel::forwardDump()
836
 
        {
837
 
   set<HMMState*> succ;
838
 
    for (Map<HMMState*, DoubleReal>::iterator it = init_prob_.begin(); it != init_prob_.end(); ++it)
839
 
    {
840
 
      succ.insert(it->first->getSuccessorStates().begin(), it->first->getSuccessorStates().end());
841
 
 
842
 
      while ( !succ.empty() )
843
 
      {
844
 
        set<HMMState*> succ_new;
845
 
        for (set<HMMState*>::const_iterator it = succ.begin(); it != succ.end(); ++it)
846
 
        {
847
 
                                        cerr << (*it)->getName() << endl;
848
 
          succ_new.insert((*it)->getSuccessorStates().begin(), (*it)->getSuccessorStates().end());
849
 
        }
850
 
        succ = succ_new;
851
 
      }
852
 
    }
853
 
        }
854
 
 
855
 
        void HiddenMarkovModel::estimateUntrainedTransitions()
856
 
        {
857
 
                #ifdef HIDDEN_MARKOV_MODEL_DEBUG
858
 
                cerr << "Transition training stats" << endl;
859
 
                for (Map<HMMState*, Map<HMMState*, Size> >::const_iterator it = training_steps_count_.begin(); it != training_steps_count_.end(); ++it)
860
 
                {
861
 
                        for (Map<HMMState*, Size>::const_iterator it1 = it->second.begin(); it1 != it->second.end(); ++it1)
862
 
                        {
863
 
                                cerr << it->first->getName() << " " << it1->first->getName() << " " << it1->second << endl;
864
 
                        }
865
 
                }
866
 
 
867
 
                cerr << "Estimation" << endl;
868
 
                #endif
869
 
 
870
 
 
871
 
                set<const Residue*> residues(ResidueDB::getInstance()->getResidues("Natural20"));
 
422
    trans_[name_to_state_[s1]][name_to_state_[s2]] = trans_prob;
 
423
    name_to_state_[s1]->addSuccessorState(name_to_state_[s2]);
 
424
    name_to_state_[s2]->addPredecessorState(name_to_state_[s1]);
 
425
    enabled_trans_[name_to_state_[s1]].insert(name_to_state_[s2]);
 
426
    training_steps_count_[name_to_state_[s1]][name_to_state_[s2]] = 0;
 
427
  }
 
428
 
 
429
  DoubleReal HiddenMarkovModel::getTransitionProbability(const String & s1, const String & s2) const
 
430
  {
 
431
    if (name_to_state_.find(s1) == name_to_state_.end())
 
432
    {
 
433
      throw Exception::ElementNotFound(__FILE__, __LINE__, __PRETTY_FUNCTION__, s1);
 
434
    }
 
435
    HMMState * state1 = name_to_state_[s1];
 
436
    if (name_to_state_.find(s2) == name_to_state_.end())
 
437
    {
 
438
      throw Exception::ElementNotFound(__FILE__, __LINE__, __PRETTY_FUNCTION__, s2);
 
439
    }
 
440
    HMMState * state2 = name_to_state_[s2];
 
441
    return getTransitionProbability_(state1, state2);
 
442
  }
 
443
 
 
444
  DoubleReal HiddenMarkovModel::getTransitionProbability_(HMMState * s1, HMMState * s2) const
 
445
  {
 
446
    HMMState * state1 = s1;
 
447
    HMMState * state2 = s2;
 
448
 
 
449
#ifdef SIMPLE_DEBUG2
 
450
    cerr << "getTransitionProbability(" << s1->getName() << ", " << s2->getName() << ")" << endl;
 
451
#endif
 
452
 
 
453
    // check if transition is a synonym transition
 
454
    if (synonym_trans_.find(state1) != synonym_trans_.end() &&
 
455
        synonym_trans_.find(state1)->second.find(state2) != synonym_trans_.find(state1)->second.end())
 
456
    {
 
457
      HMMState * tmp = synonym_trans_.find(state1)->second.find(state2)->second.first;
 
458
      state2 = synonym_trans_.find(state1)->second.find(state2)->second.second;       //name_to_state_[p.second];
 
459
      state1 = tmp;
 
460
    }
 
461
 
 
462
#ifdef SIMPLE_DEBUG2
 
463
    cerr << "getTransitionProbability: " << state1->getName() << " " << state2->getName() << endl;
 
464
#endif
 
465
 
 
466
    // get the transition prob
 
467
    if (trans_.find(state1) != trans_.end() && trans_.find(state1)->second.find(state2) != trans_.find(state1)->second.end())
 
468
    {
 
469
      /*
 
470
          // TODO !!!! ???? path length correction
 
471
          if (state2->getName().hasSubstring("next"))
 
472
          {
 
473
              return sqrt(trans_[state1][state2]);
 
474
          }
 
475
          else
 
476
          {*/
 
477
      return trans_.find(state1)->second.find(state2)->second;
 
478
      /*
 
479
  }
 
480
  */
 
481
    }
 
482
    else
 
483
    {
 
484
      return 0;
 
485
    }
 
486
  }
 
487
 
 
488
  void HiddenMarkovModel::train()
 
489
  {
 
490
#ifdef HIDDEN_MARKOV_MODEL_DEBUG
 
491
    cerr << "HiddenMarkovModel::train()" << endl;
 
492
#endif
 
493
    trained_trans_.clear();
 
494
 
 
495
#ifdef HIDDEN_MARKOV_MODEL_DEBUG
 
496
    cerr << "calc forward part" << endl;
 
497
#endif
 
498
 
 
499
    calculateForwardPart_();
 
500
 
 
501
#ifdef HIDDEN_MARKOV_MODEL_DEBUG
 
502
    cerr << "calc backward part" << endl;
 
503
#endif
 
504
 
 
505
    calculateBackwardPart_();
 
506
 
 
507
 
 
508
#ifdef HIDDEN_MARKOV_MODEL_DEBUG
 
509
    cerr << "calc px" << endl;
 
510
#endif
 
511
 
 
512
    // calc p_x from forward part
 
513
    DoubleReal px(0);
 
514
 
 
515
    for (Map<HMMState *, DoubleReal>::const_iterator it1 = train_emission_prob_.begin(); it1 != train_emission_prob_.end(); ++it1)
 
516
    {
 
517
      for (set<HMMState *>::const_iterator it2 = it1->first->getPredecessorStates().begin(); it2 != it1->first->getPredecessorStates().end(); ++it2)
 
518
      {
 
519
        px += getForwardVariable_(*it2);
 
520
      }
 
521
    }
 
522
 
 
523
    DoubleReal num_px(0);
 
524
    if (px != 0)
 
525
    {
 
526
      num_px = 1.0 / px;
 
527
    }
 
528
 
 
529
#ifdef HIDDEN_MARKOV_MODEL_DEBUG
 
530
    cerr << "px=" << num_px << endl;
 
531
    cerr << "add contributions to count_trans" << endl;
 
532
#endif
 
533
 
 
534
    // add contributions to count_trans_
 
535
    for (set<pair<HMMState *, HMMState *> >::const_iterator it = trained_trans_.begin(); it != trained_trans_.end(); ++it)
 
536
    {
 
537
      DoubleReal tmp(0);
 
538
      tmp = num_px * getForwardVariable_(it->first) * getBackwardVariable_(it->second) * getTransitionProbability_(it->first, it->second);
 
539
      tmp += pseudo_counts_;
 
540
      HMMState * s1 = it->first;
 
541
      HMMState * s2 = it->second;
 
542
 
 
543
      if (synonym_trans_.find(s1) != synonym_trans_.end() && synonym_trans_[s1].find(s2) != synonym_trans_[s1].end())
 
544
      {
 
545
        HMMState * tmp = synonym_trans_[s1][s2].first;
 
546
        s2 = synonym_trans_[s1][s2].second;
 
547
        s1 = tmp;
 
548
      }
 
549
 
 
550
      train_count_trans_all_[s1][s2].push_back(tmp);
 
551
      if (count_trans_.find(s1) != count_trans_.end() && count_trans_[s1].find(s2) != count_trans_[s1].end())
 
552
      {
 
553
        count_trans_[s1][s2] += tmp;
 
554
      }
 
555
      else
 
556
      {
 
557
        count_trans_[s1][s2] = tmp;
 
558
      }
 
559
      training_steps_count_[s1][s2]++;
 
560
    }
 
561
  }
 
562
 
 
563
  void HiddenMarkovModel::calculateForwardPart_()
 
564
  {
 
565
    forward_.clear();
 
566
    set<HMMState *> succ;
 
567
    for (Map<HMMState *, DoubleReal>::iterator it = init_prob_.begin(); it != init_prob_.end(); ++it)
 
568
    {
 
569
      //cerr << it->first << " " << it->second << endl;
 
570
      //cerr << it->first->getName() << endl;
 
571
      forward_[it->first] = it->second;
 
572
    }
 
573
 
 
574
    for (Map<HMMState *, DoubleReal>::iterator it = init_prob_.begin(); it != init_prob_.end(); ++it)
 
575
    {
 
576
      succ.insert(it->first->getSuccessorStates().begin(), it->first->getSuccessorStates().end());
 
577
 
 
578
      while (!succ.empty())
 
579
      {
 
580
        set<HMMState *> succ_new;
 
581
        for (set<HMMState *>::const_iterator it = succ.begin(); it != succ.end(); ++it)
 
582
        {
 
583
          set<HMMState *> pre = (*it)->getPredecessorStates();
 
584
          DoubleReal sum(0);
 
585
          for (set<HMMState *>::const_iterator it2 = pre.begin(); it2 != pre.end(); ++it2)
 
586
          {
 
587
#ifdef HIDDEN_MARKOV_MODEL_DEBUG
 
588
            cerr << "\tadding from " << (*it2)->getName() << " " << getForwardVariable_(*it2) << " * " << getTransitionProbability(*it2, *it) << endl;
 
589
#endif
 
590
            sum += getForwardVariable_(*it2) * getTransitionProbability_(*it2, *it);
 
591
            trained_trans_.insert(make_pair(*it2, *it));
 
592
          }
 
593
          forward_[*it] = sum;
 
594
          succ_new.insert((*it)->getSuccessorStates().begin(), (*it)->getSuccessorStates().end());
 
595
#ifdef HIDDEN_MARKOV_MODEL_DEBUG
 
596
          cerr << "f: " << (*it)->getName() << "\t" << sum << endl;
 
597
#endif
 
598
        }
 
599
        succ = succ_new;
 
600
      }
 
601
    }
 
602
  }
 
603
 
 
604
  void HiddenMarkovModel::calculateBackwardPart_()
 
605
  {
 
606
    backward_.clear();
 
607
    set<HMMState *> pre;
 
608
    for (Map<HMMState *, DoubleReal>::iterator it = train_emission_prob_.begin(); it != train_emission_prob_.end(); ++it)
 
609
    {
 
610
      backward_[it->first] = it->second;
 
611
    }
 
612
 
 
613
    for (Map<HMMState *, DoubleReal>::iterator it = train_emission_prob_.begin(); it != train_emission_prob_.end(); ++it)
 
614
    {
 
615
#ifdef HIDDEN_MARKOV_MODEL_DEBUG
 
616
      cerr << "b:" << it->first << " " <<  it->first->getName() << " " << it->second << endl;
 
617
#endif
 
618
      pre.insert(it->first->getPredecessorStates().begin(), it->first->getPredecessorStates().end());
 
619
 
 
620
      while (!pre.empty())
 
621
      {
 
622
        set<HMMState *> pre_new;
 
623
        for (set<HMMState *>::const_iterator it = pre.begin(); it != pre.end(); ++it)
 
624
        {
 
625
          set<HMMState *> succ = (*it)->getSuccessorStates();
 
626
          DoubleReal sum(0);
 
627
          for (set<HMMState *>::const_iterator it2 = succ.begin(); it2 != succ.end(); ++it2)
 
628
          {
 
629
#ifdef HIDDEN_MARKOV_MODEL_DEBUG
 
630
            cerr << "\tadding from " << (*it2)->getName() << " " << getBackwardVariable_(*it2) << " * " << getTransitionProbability(*it, *it2) << endl;
 
631
#endif
 
632
            sum += getBackwardVariable_(*it2) * getTransitionProbability_(*it, *it2);
 
633
            trained_trans_.insert(make_pair(*it, *it2));
 
634
          }
 
635
          backward_[*it] = sum;
 
636
          pre_new.insert((*it)->getPredecessorStates().begin(), (*it)->getPredecessorStates().end());
 
637
#ifdef HIDDEN_MARKOV_MODEL_DEBUG
 
638
          cerr << "b: " << (*it)->getName() << "\t" << sum << endl;
 
639
#endif
 
640
        }
 
641
        pre = pre_new;
 
642
      }
 
643
    }
 
644
  }
 
645
 
 
646
  DoubleReal HiddenMarkovModel::getForwardVariable_(HMMState * state)
 
647
  {
 
648
    return forward_.find(state) != forward_.end() ? forward_[state] : 0;
 
649
  }
 
650
 
 
651
  DoubleReal HiddenMarkovModel::getBackwardVariable_(HMMState * state)
 
652
  {
 
653
    return backward_.find(state) != backward_.end() ? backward_[state] : 0;
 
654
  }
 
655
 
 
656
  void HiddenMarkovModel::evaluate()
 
657
  {
 
658
    for (Map<HMMState *, Map<HMMState *, DoubleReal> >::const_iterator it1 = count_trans_.begin(); it1 != count_trans_.end(); ++it1)
 
659
    {
 
660
#ifdef EVALUATE_DEBUG
 
661
      cerr <<  it1->first->getName() << endl;
 
662
#endif
 
663
      DoubleReal sum(0);
 
664
      for (Map<HMMState *, DoubleReal>::const_iterator it2 = it1->second.begin(); it2 != it1->second.end(); ++it2)
 
665
      {
 
666
        if (count_trans_.find(it1->first) != count_trans_.end() &&
 
667
            count_trans_[it1->first].find(it2->first) != count_trans_[it1->first].end())
 
668
        {
 
669
          sum += count_trans_[it1->first][it2->first];
 
670
#ifdef EVALUATE_DEBUG
 
671
          cerr << it1->first->getName() << " " << it2->first->getName() << " ";
 
672
 
 
673
          //<< count_trans_[it1->first][it2->first] << endl;
 
674
          for (vector<DoubleReal>::const_iterator it = train_count_trans_all_[it1->first][it2->first].begin(); it != train_count_trans_all_[it1->first][it2->first].end(); ++it)
 
675
          {
 
676
            cerr << *it << " ";
 
677
          }
 
678
          vector<DoubleReal> data = train_count_trans_all_[it1->first][it2->first];
 
679
          std::sort(data.begin(), data.end());
 
680
          DoubleReal mean = gsl_stats_mean(&data.front(), 1, data.size());
 
681
          DoubleReal variance = gsl_stats_variance_m(&data.front(), 1, data.size(), mean);
 
682
          cerr << "mean=" << mean << ", variance=" << variance << endl;
 
683
#endif
 
684
        }
 
685
      }
 
686
 
 
687
      if (sum != 0)
 
688
      {
 
689
        for (Map<HMMState *, DoubleReal>::const_iterator it2 = it1->second.begin(); it2 != it1->second.end(); ++it2)
 
690
        {
 
691
          if (count_trans_.find(it1->first) != count_trans_.end() &&
 
692
              count_trans_[it1->first].find(it2->first) != count_trans_[it1->first].end())
 
693
          {
 
694
            trans_[it1->first][it2->first] = count_trans_[it1->first][it2->first] / sum;
 
695
          }
 
696
        }
 
697
      }
 
698
    }
 
699
  }
 
700
 
 
701
  void HiddenMarkovModel::setInitialTransitionProbability(const String & state, DoubleReal prob)
 
702
  {
 
703
    OPENMS_PRECONDITION(name_to_state_.find(state) != name_to_state_.end(), String("HiddenMarkovModel::setInitialTransitionProbability(" + state + ", " + String(prob) + "), no suchstate!").c_str());
 
704
    //cerr << state << " " << prob << endl;
 
705
    init_prob_[name_to_state_[state]] = prob;
 
706
  }
 
707
 
 
708
  void HiddenMarkovModel::clearInitialTransitionProbabilities()
 
709
  {
 
710
    init_prob_.clear();
 
711
  }
 
712
 
 
713
  void HiddenMarkovModel::setTrainingEmissionProbability(const String & state, DoubleReal prob)
 
714
  {
 
715
#ifdef SIMPLE_DEBUG2
 
716
    cerr << "setTrainingEmissionProbability(" << state << "(" << name_to_state_[state] << "), " << prob << ")" << endl;
 
717
#endif
 
718
 
 
719
    OPENMS_PRECONDITION(name_to_state_.find(state) != name_to_state_.end(), String("HiddenMarkovModel::setTrainingEmissionProbability(" + state + ", " + String(prob) + "), no such state!").c_str());
 
720
    train_emission_prob_[name_to_state_[state]] = prob;
 
721
  }
 
722
 
 
723
  void HiddenMarkovModel::clearTrainingEmissionProbabilities()
 
724
  {
 
725
    train_emission_prob_.clear();
 
726
  }
 
727
 
 
728
  void HiddenMarkovModel::enableTransition(const String & s1, const String & s2)
 
729
  {
 
730
#ifdef SIMPLE_DEBUG2
 
731
    cerr << "enableTransition: '" << s1 << "' -> '" << s2 << "'" << endl;
 
732
#endif
 
733
    enableTransition_(name_to_state_[s1], name_to_state_[s2]);
 
734
  }
 
735
 
 
736
  void HiddenMarkovModel::enableTransition_(HMMState * s1, HMMState * s2)
 
737
  {
 
738
    s1->addSuccessorState(s2);
 
739
    s2->addPredecessorState(s1);
 
740
    enabled_trans_[s1].insert(s2);
 
741
  }
 
742
 
 
743
  void HiddenMarkovModel::disableTransition(const String & s1, const String & s2)
 
744
  {
 
745
    disableTransition_(name_to_state_[s1], name_to_state_[s2]);
 
746
  }
 
747
 
 
748
  void HiddenMarkovModel::disableTransition_(HMMState * s1, HMMState * s2)
 
749
  {
 
750
    s1->deleteSuccessorState(s2);
 
751
    s2->deletePredecessorState(s1);
 
752
    enabled_trans_[s1].erase(s2);
 
753
  }
 
754
 
 
755
  void HiddenMarkovModel::disableTransitions()
 
756
  {
 
757
    for (Map<HMMState *, set<HMMState *> >::const_iterator it = enabled_trans_.begin(); it != enabled_trans_.end(); ++it)
 
758
    {
 
759
      for (set<HMMState *>::const_iterator it2 = it->second.begin(); it2 != it->second.end(); ++it2)
 
760
      {
 
761
        it->first->deleteSuccessorState(*it2);
 
762
        (*it2)->deletePredecessorState(it->first);
 
763
      }
 
764
    }
 
765
    enabled_trans_.clear();
 
766
  }
 
767
 
 
768
  void HiddenMarkovModel::calculateEmissionProbabilities(Map<HMMState *, DoubleReal> & emission_probs)
 
769
  {
 
770
    Map<HMMState *, DoubleReal> states = init_prob_;
 
771
 
 
772
    while (!states.empty())
 
773
    {
 
774
      Map<HMMState *, DoubleReal> tmp = states;
 
775
      for (Map<HMMState *, DoubleReal>::const_iterator it = tmp.begin(); it != tmp.end(); ++it)
 
776
      {
 
777
        for (set<HMMState *>::const_iterator it2 = it->first->getSuccessorStates().begin(); it2 != it->first->getSuccessorStates().end(); ++it2)
 
778
        {
 
779
          //cerr << "->" << (*it2)->getName() << "=(from " << it->first->getName() << "): "; // << states[it->first] << " * " << getTransitionProbability_(it->first, *it2) << " ";
 
780
          if (states.find(*it2) != states.end())
 
781
          {
 
782
            //cerr << " += " << states[it->first] * getTransitionProbability_(it->first, *it2);
 
783
            states[*it2] += states[it->first] * getTransitionProbability_(it->first, *it2);
 
784
          }
 
785
          else
 
786
          {
 
787
            //cerr << " = " << states[it->first] * getTransitionProbability_(it->first, *it2);
 
788
            states[*it2] = states[it->first] * getTransitionProbability_(it->first, *it2);
 
789
          }
 
790
          if (!(*it2)->isHidden())
 
791
          {
 
792
            if (emission_probs.find(*it2) != emission_probs.end())
 
793
            {
 
794
              //cerr << " emission: += " << states[it->first] * getTransitionProbability_(it->first, *it2);
 
795
              emission_probs[*it2] += states[it->first] * getTransitionProbability_(it->first, *it2);
 
796
            }
 
797
            else
 
798
            {
 
799
              //cerr << " emission: += " << states[it->first] * getTransitionProbability_(it->first, *it2);
 
800
              emission_probs[*it2] = states[it->first] * getTransitionProbability_(it->first, *it2);
 
801
            }
 
802
          }
 
803
          //cerr << endl;
 
804
        }
 
805
        states.erase(it->first);
 
806
      }
 
807
    }
 
808
 
 
809
    return;
 
810
  }
 
811
 
 
812
  void HiddenMarkovModel::dump()
 
813
  {
 
814
    cerr << "dump of transitions: " << endl;
 
815
    for (Map<HMMState *, Map<HMMState *, DoubleReal> >::const_iterator it = trans_.begin(); it != trans_.end(); ++it)
 
816
    {
 
817
      for (Map<HMMState *, DoubleReal>::const_iterator it1 = it->second.begin(); it1 != it->second.end(); ++it1)
 
818
      {
 
819
        cout << it->first->getName() << " -> " << it1->first->getName() << " " << it1->second << " " << training_steps_count_[it->first][it1->first] << ": ";
 
820
        vector<DoubleReal> all_trans = train_count_trans_all_[it->first][it1->first];
 
821
 
 
822
        if (!all_trans.empty())
 
823
        {
 
824
          DoubleReal sum = accumulate(all_trans.begin(), all_trans.end(), 0.0);
 
825
          DoubleReal avg(sum / DoubleReal(all_trans.size()));
 
826
          DoubleReal rsd(0);
 
827
          for (Size i = 0; i != all_trans.size(); ++i)
 
828
          {
 
829
            cout << all_trans[i] << " ";
 
830
            rsd += abs(all_trans[i] - avg);
 
831
          }
 
832
          cout << "rsd=" << rsd / DoubleReal(all_trans.size()) / avg;
 
833
          cout << ", avg=" << avg;
 
834
        }
 
835
 
 
836
        cout << endl;
 
837
      }
 
838
    }
 
839
 
 
840
    cerr << "dump completed" << endl;
 
841
  }
 
842
 
 
843
  void HiddenMarkovModel::forwardDump()
 
844
  {
 
845
    set<HMMState *> succ;
 
846
    for (Map<HMMState *, DoubleReal>::iterator it = init_prob_.begin(); it != init_prob_.end(); ++it)
 
847
    {
 
848
      succ.insert(it->first->getSuccessorStates().begin(), it->first->getSuccessorStates().end());
 
849
 
 
850
      while (!succ.empty())
 
851
      {
 
852
        set<HMMState *> succ_new;
 
853
        for (set<HMMState *>::const_iterator it = succ.begin(); it != succ.end(); ++it)
 
854
        {
 
855
          cerr << (*it)->getName() << endl;
 
856
          succ_new.insert((*it)->getSuccessorStates().begin(), (*it)->getSuccessorStates().end());
 
857
        }
 
858
        succ = succ_new;
 
859
      }
 
860
    }
 
861
  }
 
862
 
 
863
  void HiddenMarkovModel::estimateUntrainedTransitions()
 
864
  {
 
865
#ifdef HIDDEN_MARKOV_MODEL_DEBUG
 
866
    cerr << "Transition training stats" << endl;
 
867
    for (Map<HMMState *, Map<HMMState *, Size> >::const_iterator it = training_steps_count_.begin(); it != training_steps_count_.end(); ++it)
 
868
    {
 
869
      for (Map<HMMState *, Size>::const_iterator it1 = it->second.begin(); it1 != it->second.end(); ++it1)
 
870
      {
 
871
        cerr << it->first->getName() << " " << it1->first->getName() << " " << it1->second << endl;
 
872
      }
 
873
    }
 
874
 
 
875
    cerr << "Estimation" << endl;
 
876
#endif
 
877
 
 
878
 
 
879
    set<const Residue *> residues(ResidueDB::getInstance()->getResidues("Natural20"));
872
880
    for (StringList::const_iterator it = var_modifications_.begin(); it != var_modifications_.end(); ++it)
873
881
    {
874
882
      residues.insert(ResidueDB::getInstance()->getModifiedResidue(*it));
875
883
    }
876
884
 
877
 
                // pathways axyz and bxyz and the first two explicitely modeled ones
878
 
                HMMState* s2 = 0;
879
 
                HMMState* end_state = name_to_state_["end"];
880
 
                StringList pathways = StringList::create("axyz,axyz1,axyz1,bxyz,bxyz1,bxyz2");  
881
 
                for (StringList::const_iterator pathway_it = pathways.begin(); pathway_it != pathways.end(); ++pathway_it)
882
 
                {
883
 
                        String pathway = *pathway_it;
884
 
                        for (set<const Residue*>::const_iterator it = residues.begin(); it != residues.end(); ++it)
885
 
                        {
886
 
                                s2 = name_to_state_[pathway];
887
 
                                for (set<const Residue*>::const_iterator jt = residues.begin(); jt != residues.end(); ++jt)
888
 
                                {
889
 
                                        AASequence first_aa, second_aa;
890
 
                                        first_aa += *it;
891
 
                                        second_aa += *jt;
892
 
                                        String aa1(first_aa.toString()), aa2(second_aa.toString());
893
 
                                        #ifdef HIDDEN_MARKOV_MODEL_DEBUG
894
 
                                        cerr << "Estimating: " << aa1 << " -> " << aa2 << " (" << training_steps_count_[name_to_state_[aa1 + aa2 + "_" + pathway]][s2] << ") :: " ;
895
 
                                        #endif
896
 
 
897
 
                                        if (training_steps_count_[name_to_state_[aa1 + aa2 + "_" + pathway]][s2] == 0)
898
 
                                        {
899
 
                                                Size count(0);
900
 
                                                DoubleReal sum(0);
901
 
                                                // "rows" of the amino acid matrix
902
 
                                                for (set<const Residue*>::const_iterator kt = residues.begin(); kt != residues.end(); ++kt)
903
 
                                                {
904
 
                                                        AASequence third_aa;
905
 
                                                        third_aa += *kt;
906
 
                                                        String aa3(third_aa.toString());
907
 
                                                        if (training_steps_count_[name_to_state_[aa1 + aa3 + "_" + pathway]][s2] != 0)
908
 
                                                        {
909
 
                                                                sum += trans_[name_to_state_[aa1 + aa3 + "_" + pathway]][s2];
910
 
                                                                count++;
911
 
                                                        }
912
 
                                                }
913
 
                                                // "columns" of the amino acid matrix
914
 
                                                for (set<const Residue*>::const_iterator kt = residues.begin(); kt != residues.end(); ++kt)
915
 
                                                {
916
 
                                                        AASequence third_aa;
917
 
                                                        third_aa += *kt;
918
 
                                                        String aa3(third_aa.toString());
919
 
 
920
 
                                                        if (training_steps_count_[name_to_state_[aa3 + aa2 + "_" + pathway]][s2] != 0)
921
 
                                                        {
922
 
                                                                sum += trans_[name_to_state_[aa3 + aa2 + "_" + pathway]][s2];
923
 
                                                                count++;
924
 
                                                        }
925
 
                                                }
926
 
 
927
 
                                                #ifdef HIDDEN_MARKOV_MODEL_DEBUG
928
 
                                                cerr << trans_[name_to_state_[aa1 + aa2 + "_" + pathway]][s2] << " to ";
929
 
                                                #endif
930
 
                                                if (count != 0)
931
 
                                                {
932
 
                                                        #ifdef HIDDEN_MARKOV_MODEL_DEBUG
933
 
                                                        cerr << "setting transitions of " << aa1 << aa2 << "_" << pathway << " -> " << pathway << " to " << sum/DoubleReal(count) << endl;
934
 
                                                        #endif
935
 
                                                        trans_[name_to_state_[aa1 + aa2 + "_" + pathway]][s2] = sum/DoubleReal(count);
936
 
                                                        trans_[name_to_state_[aa1 + aa2 + "_" + pathway]][end_state] = 1 - sum/DoubleReal(count);
937
 
                                                }
938
 
                                                #ifdef HIDDEN_MARKOV_MODEL_DEBUG
939
 
                                                cerr << sum/DoubleReal(count) << endl;
940
 
                                        }
941
 
                                        else
942
 
                                        {
943
 
                                                cerr << "not needed" << endl;
944
 
                                                #endif
945
 
                                        }
946
 
                                }
947
 
                        }
948
 
                }
949
 
 
950
 
                // sc and cr
951
 
                String sc_residues("HKRDE");
952
 
                for (String::ConstIterator it = sc_residues.begin(); it != sc_residues.end(); ++it)
953
 
                {
954
 
                        String sc_res = *it;
955
 
 
956
 
                        for (set<const Residue*>::const_iterator jt = residues.begin(); jt != residues.end(); ++jt)
957
 
                        {
 
885
    // pathways axyz and bxyz and the first two explicitly modeled ones
 
886
    HMMState * s2 = 0;
 
887
    HMMState * end_state = name_to_state_["end"];
 
888
    StringList pathways = StringList::create("axyz,axyz1,axyz1,bxyz,bxyz1,bxyz2");
 
889
    for (StringList::const_iterator pathway_it = pathways.begin(); pathway_it != pathways.end(); ++pathway_it)
 
890
    {
 
891
      String pathway = *pathway_it;
 
892
      for (set<const Residue *>::const_iterator it = residues.begin(); it != residues.end(); ++it)
 
893
      {
 
894
        s2 = name_to_state_[pathway];
 
895
        for (set<const Residue *>::const_iterator jt = residues.begin(); jt != residues.end(); ++jt)
 
896
        {
 
897
          AASequence first_aa, second_aa;
 
898
          first_aa += *it;
 
899
          second_aa += *jt;
 
900
          String aa1(first_aa.toString()), aa2(second_aa.toString());
 
901
#ifdef HIDDEN_MARKOV_MODEL_DEBUG
 
902
          cerr << "Estimating: " << aa1 << " -> " << aa2 << " (" << training_steps_count_[name_to_state_[aa1 + aa2 + "_" + pathway]][s2] << ") :: ";
 
903
#endif
 
904
 
 
905
          if (training_steps_count_[name_to_state_[aa1 + aa2 + "_" + pathway]][s2] == 0)
 
906
          {
 
907
            Size count(0);
 
908
            DoubleReal sum(0);
 
909
            // "rows" of the amino acid matrix
 
910
            for (set<const Residue *>::const_iterator kt = residues.begin(); kt != residues.end(); ++kt)
 
911
            {
 
912
              AASequence third_aa;
 
913
              third_aa += *kt;
 
914
              String aa3(third_aa.toString());
 
915
              if (training_steps_count_[name_to_state_[aa1 + aa3 + "_" + pathway]][s2] != 0)
 
916
              {
 
917
                sum += trans_[name_to_state_[aa1 + aa3 + "_" + pathway]][s2];
 
918
                count++;
 
919
              }
 
920
            }
 
921
            // "columns" of the amino acid matrix
 
922
            for (set<const Residue *>::const_iterator kt = residues.begin(); kt != residues.end(); ++kt)
 
923
            {
 
924
              AASequence third_aa;
 
925
              third_aa += *kt;
 
926
              String aa3(third_aa.toString());
 
927
 
 
928
              if (training_steps_count_[name_to_state_[aa3 + aa2 + "_" + pathway]][s2] != 0)
 
929
              {
 
930
                sum += trans_[name_to_state_[aa3 + aa2 + "_" + pathway]][s2];
 
931
                count++;
 
932
              }
 
933
            }
 
934
 
 
935
#ifdef HIDDEN_MARKOV_MODEL_DEBUG
 
936
            cerr << trans_[name_to_state_[aa1 + aa2 + "_" + pathway]][s2] << " to ";
 
937
#endif
 
938
            if (count != 0)
 
939
            {
 
940
#ifdef HIDDEN_MARKOV_MODEL_DEBUG
 
941
              cerr << "setting transitions of " << aa1 << aa2 << "_" << pathway << " -> " << pathway << " to " << sum / DoubleReal(count) << endl;
 
942
#endif
 
943
              trans_[name_to_state_[aa1 + aa2 + "_" + pathway]][s2] = sum / DoubleReal(count);
 
944
              trans_[name_to_state_[aa1 + aa2 + "_" + pathway]][end_state] = 1 - sum / DoubleReal(count);
 
945
            }
 
946
#ifdef HIDDEN_MARKOV_MODEL_DEBUG
 
947
            cerr << sum / DoubleReal(count) << endl;
 
948
          }
 
949
          else
 
950
          {
 
951
            cerr << "not needed" << endl;
 
952
#endif
 
953
          }
 
954
        }
 
955
      }
 
956
    }
 
957
 
 
958
    // sc and cr
 
959
    String sc_residues("HKRDE");
 
960
    for (String::ConstIterator it = sc_residues.begin(); it != sc_residues.end(); ++it)
 
961
    {
 
962
      String sc_res = *it;
 
963
 
 
964
      for (set<const Residue *>::const_iterator jt = residues.begin(); jt != residues.end(); ++jt)
 
965
      {
958
966
        AASequence second_aa;
959
 
                                second_aa += *jt;
960
 
                                String aa2(second_aa.toString());
961
 
                        
962
 
                                s2 = name_to_state_[sc_res];
963
 
                                if (training_steps_count_[name_to_state_[aa2 + "_" + sc_res]][s2] == 0)
964
 
                                {
965
 
                                        Size count(0);
966
 
                                        DoubleReal sum(0);
967
 
                                        for (set<const Residue*>::const_iterator kt = residues.begin(); kt != residues.end(); ++kt)
968
 
                                        {
969
 
                                                AASequence third_aa;
970
 
                                                third_aa += *kt;
971
 
                                                String aa3(third_aa.toString());
972
 
 
973
 
                                                HMMState* s1 = name_to_state_[aa3 + "_" + sc_res];
974
 
                                                if (training_steps_count_[s1][s2] != 0)
975
 
                                                {
976
 
                                                        sum += trans_[s1][s2];
977
 
                                                        count++;
978
 
                                                }
979
 
                                        }
980
 
 
981
 
                                        if (count != 0)
982
 
                                        {
983
 
                                                trans_[name_to_state_[aa2 + "_" + sc_res]][s2] = sum/DoubleReal(count);
984
 
                                                trans_[name_to_state_[aa2 + "_" + sc_res]][end_state] = 1 - sum/DoubleReal(count);
985
 
                                        }
986
 
                                }
987
 
                        }
988
 
                }
989
 
 
990
 
                StringList bk_pathways = StringList::create("bk-1,bk-2");
991
 
 
992
 
                for (StringList::const_iterator pathway_it = bk_pathways.begin(); pathway_it != bk_pathways.end(); ++pathway_it)
993
 
                {
994
 
                        String pathway = *pathway_it;
995
 
                        s2 = name_to_state_[pathway];
996
 
                        for (set<const Residue*>::const_iterator it = residues.begin(); it != residues.end(); ++it)
997
 
                        {
998
 
                                AASequence first_aa;
999
 
                                first_aa += *it;
1000
 
                                String aa1(first_aa.toString());
1001
 
                                if (training_steps_count_[name_to_state_[aa1 + "_" + pathway]][s2] == 0)
1002
 
                                {
1003
 
                                        Size count(0);
1004
 
                                        DoubleReal sum(0);
1005
 
                                        for (set<const Residue*>::const_iterator jt = residues.begin(); jt != residues.end(); ++jt)
1006
 
                                        {
1007
 
                                                AASequence second_aa;
1008
 
                                                second_aa += *jt;
1009
 
                                                String aa2(second_aa.toString());
1010
 
                                                HMMState* s1 = name_to_state_[aa2 + "_" + pathway];
1011
 
                                                if (training_steps_count_[s1][s2] != 0)
1012
 
                                                {
1013
 
                                                        sum += trans_[s1][s2];
1014
 
                                                        count++;
1015
 
                                                }
1016
 
                                                //cerr << "Estimating transition of '" << aa1 << pathway << "' -> '" << pathway << "' to " << sum/(DoubleReal)count << endl;
1017
 
                                                if (count != 0)
1018
 
                                                {
1019
 
                                                        trans_[name_to_state_[aa1 + "_" + pathway]][s2] = sum/(DoubleReal)count;
1020
 
                                                        trans_[name_to_state_[aa1 + "_" + pathway]][end_state] = 1 - sum/(DoubleReal)count;
1021
 
                                                }
1022
 
                                        }
1023
 
                                }
1024
 
                        }
1025
 
                }
1026
 
 
1027
 
        }
1028
 
 
1029
 
        void HiddenMarkovModel::setPseudoCounts(DoubleReal pseudo_counts)
1030
 
        {
1031
 
                pseudo_counts_ = pseudo_counts;
1032
 
                return;
1033
 
        }
1034
 
 
1035
 
        DoubleReal HiddenMarkovModel::getPseudoCounts() const
1036
 
        {
1037
 
                return pseudo_counts_;
1038
 
        }
1039
 
 
1040
 
        void HiddenMarkovModel::setVariableModifications(const StringList& modifications)
1041
 
        {
1042
 
                var_modifications_ = modifications;
1043
 
        }
1044
 
 
1045
 
        void HiddenMarkovModel::copy_(const HiddenMarkovModel& source)
1046
 
        {
1047
 
                Map<HMMState*, HMMState*> old_to_new;
1048
 
    for (set<HMMState*>::const_iterator it = source.states_.begin(); it != source.states_.end(); ++it)
1049
 
    {
1050
 
      HMMState* s = new HMMState(**it);
 
967
        second_aa += *jt;
 
968
        String aa2(second_aa.toString());
 
969
 
 
970
        s2 = name_to_state_[sc_res];
 
971
        if (training_steps_count_[name_to_state_[aa2 + "_" + sc_res]][s2] == 0)
 
972
        {
 
973
          Size count(0);
 
974
          DoubleReal sum(0);
 
975
          for (set<const Residue *>::const_iterator kt = residues.begin(); kt != residues.end(); ++kt)
 
976
          {
 
977
            AASequence third_aa;
 
978
            third_aa += *kt;
 
979
            String aa3(third_aa.toString());
 
980
 
 
981
            HMMState * s1 = name_to_state_[aa3 + "_" + sc_res];
 
982
            if (training_steps_count_[s1][s2] != 0)
 
983
            {
 
984
              sum += trans_[s1][s2];
 
985
              count++;
 
986
            }
 
987
          }
 
988
 
 
989
          if (count != 0)
 
990
          {
 
991
            trans_[name_to_state_[aa2 + "_" + sc_res]][s2] = sum / DoubleReal(count);
 
992
            trans_[name_to_state_[aa2 + "_" + sc_res]][end_state] = 1 - sum / DoubleReal(count);
 
993
          }
 
994
        }
 
995
      }
 
996
    }
 
997
 
 
998
    StringList bk_pathways = StringList::create("bk-1,bk-2");
 
999
 
 
1000
    for (StringList::const_iterator pathway_it = bk_pathways.begin(); pathway_it != bk_pathways.end(); ++pathway_it)
 
1001
    {
 
1002
      String pathway = *pathway_it;
 
1003
      s2 = name_to_state_[pathway];
 
1004
      for (set<const Residue *>::const_iterator it = residues.begin(); it != residues.end(); ++it)
 
1005
      {
 
1006
        AASequence first_aa;
 
1007
        first_aa += *it;
 
1008
        String aa1(first_aa.toString());
 
1009
        if (training_steps_count_[name_to_state_[aa1 + "_" + pathway]][s2] == 0)
 
1010
        {
 
1011
          Size count(0);
 
1012
          DoubleReal sum(0);
 
1013
          for (set<const Residue *>::const_iterator jt = residues.begin(); jt != residues.end(); ++jt)
 
1014
          {
 
1015
            AASequence second_aa;
 
1016
            second_aa += *jt;
 
1017
            String aa2(second_aa.toString());
 
1018
            HMMState * s1 = name_to_state_[aa2 + "_" + pathway];
 
1019
            if (training_steps_count_[s1][s2] != 0)
 
1020
            {
 
1021
              sum += trans_[s1][s2];
 
1022
              count++;
 
1023
            }
 
1024
            //cerr << "Estimating transition of '" << aa1 << pathway << "' -> '" << pathway << "' to " << sum/(DoubleReal)count << endl;
 
1025
            if (count != 0)
 
1026
            {
 
1027
              trans_[name_to_state_[aa1 + "_" + pathway]][s2] = sum / (DoubleReal)count;
 
1028
              trans_[name_to_state_[aa1 + "_" + pathway]][end_state] = 1 - sum / (DoubleReal)count;
 
1029
            }
 
1030
          }
 
1031
        }
 
1032
      }
 
1033
    }
 
1034
 
 
1035
  }
 
1036
 
 
1037
  void HiddenMarkovModel::setPseudoCounts(DoubleReal pseudo_counts)
 
1038
  {
 
1039
    pseudo_counts_ = pseudo_counts;
 
1040
    return;
 
1041
  }
 
1042
 
 
1043
  DoubleReal HiddenMarkovModel::getPseudoCounts() const
 
1044
  {
 
1045
    return pseudo_counts_;
 
1046
  }
 
1047
 
 
1048
  void HiddenMarkovModel::setVariableModifications(const StringList & modifications)
 
1049
  {
 
1050
    var_modifications_ = modifications;
 
1051
  }
 
1052
 
 
1053
  void HiddenMarkovModel::copy_(const HiddenMarkovModel & source)
 
1054
  {
 
1055
    Map<HMMState *, HMMState *> old_to_new;
 
1056
    for (set<HMMState *>::const_iterator it = source.states_.begin(); it != source.states_.end(); ++it)
 
1057
    {
 
1058
      HMMState * s = new HMMState(**it);
1051
1059
      states_.insert(s);
1052
1060
      name_to_state_[s->getName()] = s;
1053
1061
      old_to_new[*it] = s;
1054
1062
    }
1055
 
                
1056
 
                // trans_
1057
 
                for (Map<HMMState*, Map<HMMState*, DoubleReal> >::const_iterator it1 = source.trans_.begin(); it1 != source.trans_.end(); ++it1)
1058
 
                {
1059
 
                        for (Map<HMMState*, DoubleReal>::const_iterator it2 = it1->second.begin(); it2 != it1->second.end(); ++it2)
1060
 
                        {
1061
 
                                trans_[old_to_new[it1->first]][old_to_new[it2->first]] = it2->second;
1062
 
                        }
1063
 
                }
1064
 
 
1065
 
                // count_trans_
1066
 
    for (Map<HMMState*, Map<HMMState*, DoubleReal> >::const_iterator it1 = source.count_trans_.begin(); it1 != source.count_trans_.end(); ++it1)
1067
 
    {
1068
 
      for (Map<HMMState*, DoubleReal>::const_iterator it2 = it1->second.begin(); it2 != it1->second.end(); ++it2)
 
1063
 
 
1064
    // trans_
 
1065
    for (Map<HMMState *, Map<HMMState *, DoubleReal> >::const_iterator it1 = source.trans_.begin(); it1 != source.trans_.end(); ++it1)
 
1066
    {
 
1067
      for (Map<HMMState *, DoubleReal>::const_iterator it2 = it1->second.begin(); it2 != it1->second.end(); ++it2)
 
1068
      {
 
1069
        trans_[old_to_new[it1->first]][old_to_new[it2->first]] = it2->second;
 
1070
      }
 
1071
    }
 
1072
 
 
1073
    // count_trans_
 
1074
    for (Map<HMMState *, Map<HMMState *, DoubleReal> >::const_iterator it1 = source.count_trans_.begin(); it1 != source.count_trans_.end(); ++it1)
 
1075
    {
 
1076
      for (Map<HMMState *, DoubleReal>::const_iterator it2 = it1->second.begin(); it2 != it1->second.end(); ++it2)
1069
1077
      {
1070
1078
        count_trans_[old_to_new[it1->first]][old_to_new[it2->first]] = it2->second;
1071
1079
      }
1072
1080
    }
1073
1081
 
1074
 
                for (Map<HMMState*, Map<HMMState*, std::vector<DoubleReal> > >::const_iterator it1 = source.train_count_trans_all_.begin(); it1 != source.train_count_trans_all_.end(); ++it1)
1075
 
                {
1076
 
                        for (Map<HMMState*, vector<DoubleReal> >::const_iterator it2 = it1->second.begin(); it2 != it1->second.end(); ++it2)
 
1082
    for (Map<HMMState *, Map<HMMState *, std::vector<DoubleReal> > >::const_iterator it1 = source.train_count_trans_all_.begin(); it1 != source.train_count_trans_all_.end(); ++it1)
 
1083
    {
 
1084
      for (Map<HMMState *, vector<DoubleReal> >::const_iterator it2 = it1->second.begin(); it2 != it1->second.end(); ++it2)
1077
1085
      {
1078
1086
        train_count_trans_all_[old_to_new[it1->first]][old_to_new[it2->first]] = it2->second;
1079
1087
      }
1080
1088
 
1081
 
                }
1082
 
 
1083
 
 
1084
 
                for (Map<HMMState*, Map<HMMState*, Size> >::const_iterator it1 = source.training_steps_count_.begin(); it1 != source.training_steps_count_.end(); ++it1)
 
1089
    }
 
1090
 
 
1091
 
 
1092
    for (Map<HMMState *, Map<HMMState *, Size> >::const_iterator it1 = source.training_steps_count_.begin(); it1 != source.training_steps_count_.end(); ++it1)
1085
1093
    {
1086
 
      for (Map<HMMState*, Size>::const_iterator it2 = it1->second.begin(); it2 != it1->second.end(); ++it2)
 
1094
      for (Map<HMMState *, Size>::const_iterator it2 = it1->second.begin(); it2 != it1->second.end(); ++it2)
1087
1095
      {
1088
1096
        training_steps_count_[old_to_new[it1->first]][old_to_new[it2->first]] = it2->second;
1089
1097
      }
1090
1098
    }
1091
1099
 
1092
 
                // forward and backward are just temporary objects
1093
 
 
1094
 
                for (Map<HMMState*, DoubleReal>::const_iterator it = source.train_emission_prob_.begin(); it != source.train_emission_prob_.end(); ++it)
1095
 
                {
1096
 
                        train_emission_prob_[old_to_new[it->first]] = it->second;
1097
 
                }
1098
 
 
1099
 
                for (Map<HMMState*, DoubleReal>::const_iterator it = source.init_prob_.begin(); it != source.init_prob_.end(); ++it)
1100
 
                {
1101
 
                        init_prob_[old_to_new[it->first]] = it->second;
1102
 
                }
1103
 
 
1104
 
                for (std::set<std::pair<HMMState*, HMMState*> >::const_iterator it = source.trained_trans_.begin(); it != source.trained_trans_.end(); ++it)
1105
 
                {
1106
 
                        trained_trans_.insert(make_pair(old_to_new[it->first], old_to_new[it->second]));
1107
 
                }
1108
 
 
1109
 
                synonym_trans_names_ = source.synonym_trans_names_;
1110
 
                pseudo_counts_ = source.pseudo_counts_;
1111
 
                var_modifications_ = source.var_modifications_;
 
1100
    // forward and backward are just temporary objects
 
1101
 
 
1102
    for (Map<HMMState *, DoubleReal>::const_iterator it = source.train_emission_prob_.begin(); it != source.train_emission_prob_.end(); ++it)
 
1103
    {
 
1104
      train_emission_prob_[old_to_new[it->first]] = it->second;
 
1105
    }
 
1106
 
 
1107
    for (Map<HMMState *, DoubleReal>::const_iterator it = source.init_prob_.begin(); it != source.init_prob_.end(); ++it)
 
1108
    {
 
1109
      init_prob_[old_to_new[it->first]] = it->second;
 
1110
    }
 
1111
 
 
1112
    for (std::set<std::pair<HMMState *, HMMState *> >::const_iterator it = source.trained_trans_.begin(); it != source.trained_trans_.end(); ++it)
 
1113
    {
 
1114
      trained_trans_.insert(make_pair(old_to_new[it->first], old_to_new[it->second]));
 
1115
    }
 
1116
 
 
1117
    synonym_trans_names_ = source.synonym_trans_names_;
 
1118
    pseudo_counts_ = source.pseudo_counts_;
 
1119
    var_modifications_ = source.var_modifications_;
1112
1120
 
1113
1121
 
1114
1122
    for (Map<String, Map<String, std::pair<String, String> > >::const_iterator it = synonym_trans_names_.begin(); it != synonym_trans_names_.end(); ++it)
1116
1124
      for (Map<String, pair<String, String> >::const_iterator it2 = it->second.begin(); it2 != it->second.end(); ++it2)
1117
1125
      {
1118
1126
        synonym_trans_[name_to_state_[it->first]][name_to_state_[it2->first]] =
1119
 
            make_pair(name_to_state_[it2->second.first], name_to_state_[it2->second.second]);
1120
 
      }
1121
 
    }
1122
 
 
1123
 
 
1124
 
                for (Map<HMMState*, set<HMMState*> >::const_iterator it1 = source.enabled_trans_.begin(); it1 != source.enabled_trans_.end(); ++it1)
1125
 
                {
1126
 
                        for (set<HMMState*>::const_iterator it2 = it1->second.begin(); it2 != it1->second.end(); ++it2)
1127
 
                        {
1128
 
                                enabled_trans_[old_to_new[it1->first]].insert(old_to_new[*it2]);
1129
 
                        }
1130
 
                }
1131
 
 
1132
 
        }
 
1127
          make_pair(name_to_state_[it2->second.first], name_to_state_[it2->second.second]);
 
1128
      }
 
1129
    }
 
1130
 
 
1131
 
 
1132
    for (Map<HMMState *, set<HMMState *> >::const_iterator it1 = source.enabled_trans_.begin(); it1 != source.enabled_trans_.end(); ++it1)
 
1133
    {
 
1134
      for (set<HMMState *>::const_iterator it2 = it1->second.begin(); it2 != it1->second.end(); ++it2)
 
1135
      {
 
1136
        enabled_trans_[old_to_new[it1->first]].insert(old_to_new[*it2]);
 
1137
      }
 
1138
    }
 
1139
 
 
1140
  }
 
1141
 
1133
1142
}
1134
 
 
1135