Index: /liacs/dbdm/dbdm_4/UnfairCasino.py
===================================================================
--- /liacs/dbdm/dbdm_4/UnfairCasino.py	(revision 60)
+++ /liacs/dbdm/dbdm_4/UnfairCasino.py	(revision 60)
@@ -0,0 +1,9 @@
+from ghmm import *
+# File: UnfairCasino.py example from the GHMM tutorial
+#
+alph = IntegerRange(1,7)
+train_seq = EmissionSequence(alph, [1, 6, 2, 5, 5, 4, 5, 1, 2, 1, 3, 6, 6, 3, 2, 1, 4, 4, 1, 1, 4, 2, 1, 1, 6, 3, 3, 2, 1, 4, 4, 3, 3, 5, 3, 3, 3, 3, 4, 3, 1, 5, 4, 1, 4, 5, 1, 1, 3, 4, 3, 5, 5, 1, 5, 2, 1, 5, 3, 6, 3, 6, 5, 6, 5, 3, 3, 1, 2, 6, 3, 3, 2, 2, 5, 4, 1, 5, 6, 3, 3, 5, 1, 5, 2, 3, 1, 1, 1, 5, 6, 4, 5, 5, 1, 6, 2, 6, 5, 3, 1, 1, 3, 3, 1, 1, 2, 5, 2, 3, 2, 4, 1, 5, 5, 5, 4, 6, 5, 6, 3, 1, 6, 1, 5, 4, 3, 1, 3, 1, 2, 6, 3, 5, 2, 1, 3, 6, 4, 4, 4, 4, 3, 5, 1, 6, 4, 4, 3, 5, 1, 5, 5, 5, 5, 6, 5, 1, 6, 1, 1, 4, 1, 4, 2, 6, 6, 2, 5, 4, 5, 5, 4, 3, 3, 6, 5, 4, 1, 5, 3, 3, 3, 2, 5, 6, 2, 3, 3, 5, 1, 3, 4, 5, 1, 6, 3, 2, 4, 5, 2, 2, 5, 1, 4, 5, 1, 5, 6, 5, 3, 4, 6, 1, 3, 2, 4, 6, 6, 3, 1, 6, 5, 2, 5, 4, 4, 4, 4, 2, 6, 3, 3, 2, 1, 1, 6, 5, 3, 3, 3, 4, 5, 6, 1, 6, 2, 3, 1, 4, 6, 5, 3, 1, 2, 4, 6, 2, 6, 2, 1, 2, 6, 5, 6, 1, 4, 4, 1, 5, 5, 5, 3, 4, 4, 5, 2, 2, 5, 1, 2, 1, 3, 2, 3, 3, 3, 5, 3, 5, 2, 3, 5, 5, 5, 2, 6, 2, 5, 6, 4, 5, 4, 4, 3, 3, 6, 2, 3, 2, 2, 6, 1, 1, 1, 4, 1, 2, 6, 1, 4, 5, 2, 6, 2, 6, 6, 3, 2, 4, 2, 2, 2, 3, 3, 4, 4, 1, 1, 1, 5, 3, 4, 3, 5, 3, 3, 3, 4, 6, 1, 1, 6, 4, 2, 4, 6, 4, 1, 5, 1, 4, 5, 5, 4, 6, 1, 3, 5, 1, 5, 3, 1, 6, 6, 3, 3, 2, 4, 4, 3, 1, 2, 5, 3, 5, 4, 5, 3, 2, 2, 6, 4, 4, 2, 6, 4, 6, 1, 4, 6, 3, 3, 4, 6, 3, 2, 3, 5, 3, 5, 4, 2, 2, 2, 1, 4, 1, 5, 1, 1, 1, 6, 6, 3, 4, 4, 5, 1, 2, 4, 3, 3, 1, 1, 5, 6, 1, 3, 4, 6, 2, 4, 5, 2, 1, 1, 2, 6, 5, 4, 3, 5, 4, 5, 2, 5, 1, 1, 5, 4, 3, 6, 5, 3, 6, 4, 2, 3, 5, 1, 4, 6, 1, 2, 6, 1, 3, 5, 2, 1, 2, 6, 1, 2, 5, 4, 4, 5, 5, 1, 4, 3, 1, 6, 5, 4, 3, 4, 2, 4, 2, 6, 6, 3, 3, 6, 1, 5, 6, 4, 1, 3, 2, 2, 1, 3, 4, 6, 3, 5, 5, 1, 5, 5, 2, 5, 3, 3, 1, 3, 4, 2, 1, 1, 2, 1, 5, 1, 5, 4, 2, 6, 3, 5, 6, 6, 2, 2, 5, 6, 1, 5, 3, 1, 2, 1, 2, 3, 4, 5, 5, 4, 3, 3, 3, 4, 6, 2, 5, 2, 5, 5, 3, 1, 4, 3, 6, 2, 5, 6, 3, 4, 2, 6, 2, 5, 1, 3, 1, 2, 2, 1, 4, 1, 6, 1, 4, 5, 5, 5, 2, 3, 2, 5, 1, 6, 2, 5, 5, 1, 4, 4, 6, 1, 3, 1, 1, 4, 3, 1, 2, 1, 1, 1, 4, 2, 2, 2, 1, 2, 3, 5, 5, 4, 4, 4, 6, 3, 6, 4, 3, 5, 2, 4, 4, 5, 6, 2, 6, 5, 5, 3, 1, 5, 4, 4, 1, 5, 3, 3, 1, 2, 6, 3, 6, 6, 6, 6, 1, 1, 3, 2, 5, 6, 4, 3, 4, 6, 5, 2, 2, 1, 2, 5, 5, 6, 1, 5, 2, 5, 6, 2, 4, 5, 4, 2, 5, 5, 4, 3, 4, 1, 6, 1, 1, 6, 2, 2, 1, 3, 4, 3, 1, 1, 6, 1, 1, 1, 4, 4, 4, 1, 5, 6, 4, 4, 2, 1, 4, 6, 6, 2, 3, 6, 2, 4, 3, 4, 3, 3, 3, 4, 1, 2, 4, 3, 3, 6, 2, 4, 1, 1, 2, 5, 4, 6, 3, 1, 2, 1, 3, 2, 1, 3, 6, 2, 2, 2, 5, 1, 1, 6, 6, 6, 1, 4, 2, 6, 2, 3, 3, 2, 4, 1, 2, 5, 2, 6, 3, 5, 6, 2, 1, 3, 6, 5, 1, 4, 1, 4, 6, 4, 4, 4, 2, 2, 3, 6, 4, 6, 5, 4, 2, 2, 5, 5, 3, 4, 2, 3, 3, 6, 6, 3, 1, 5, 5, 5, 6, 4, 2, 6, 5, 4, 5, 6, 3, 1, 1, 5, 4, 5, 5, 3, 5, 3, 2, 2, 1, 3, 3, 6, 5, 3, 5, 1, 2, 5, 2, 4, 1, 5, 2, 1, 2, 5, 3, 1, 3, 4, 4, 1, 5, 1, 3, 1, 2, 6, 3, 6, 1, 3, 5, 2, 6, 6, 5, 2, 3, 6, 3, 3, 4, 5, 1, 6, 4, 2, 4, 1, 4, 3, 3, 4, 2, 3, 1, 6, 2, 4, 3, 5, 2, 2, 4, 5, 1, 3, 2, 2, 1, 6, 3, 2, 1, 2, 1, 5, 4, 2, 1, 4, 5, 5, 3, 4, 4, 6, 5, 3, 6, 5, 3, 6, 5, 6, 3, 3, 5, 3, 3, 2, 4, 2, 6, 3, 2, 6, 5, 5, 5, 3, 3, 1, 2, 5, 1, 3, 6, 3, 3, 5, 1, 4, 4, 4, 2, 3, 4, 2, 6, 1, 6, 1, 4, 3, 2, 1, 3, 5, 4, 6, 5, 6, 6, 6, 4, 2, 4, 3, 5, 1, 1, 1, 5, 4, 6, 2, 2, 4, 3, 2, 1, 6, 2, 1, 6, 1, 6, 5, 1, 6, 4, 3, 6, 5, 4, 4, 1, 3, 4, 3, 1, 6, 4, 4, 2, 1, 2, 5, 5, 1, 4, 3, 2, 3, 2, 1, 4, 6, 3, 6, 6, 1, 5, 3, 6, 1, 3, 6, 6, 5, 3, 2, 4, 6, 2, 4, 4, 4, 2, 4, 4, 3, 6, 5, 3, 4, 1, 3, 3, 1, 6, 1, 4, 1, 6, 3, 5, 2, 1, 5, 4, 2, 4, 6, 1, 4, 3, 1, 1, 3, 2, 5, 5, 1, 1, 5, 3, 4, 3, 1, 3, 2, 5, 6, 2, 1, 3, 2, 6, 3, 6, 4, 4, 3, 5, 3, 2, 5, 2, 3, 2, 4, 2, 6, 2, 4, 6, 5, 5, 2, 3, 5, 6, 4, 3, 1, 3, 3, 2, 2, 2, 3, 6, 1, 6, 3, 1, 6, 3, 1, 3, 1, 1, 1, 4, 3, 1, 5, 4, 6, 2, 1, 6, 2, 2, 2, 1, 5, 5, 1, 2, 5, 5, 2, 5, 2, 4, 5, 1, 5, 5, 6, 6, 5, 3, 5, 6, 2, 5, 5, 5, 1, 2, 2, 4, 2, 4, 5, 5, 1, 4, 1, 5, 3, 5, 1, 4, 2, 1, 2, 2, 2, 4, 4, 4, 4, 2, 1, 1, 4, 5, 4, 1, 2, 2, 3, 5, 6, 4, 1, 1, 1, 3, 6, 1, 2, 4, 3, 2, 3, 2, 3, 6, 6, 3, 4, 4, 4, 6, 6, 2, 6, 6, 3, 2, 3, 5, 1, 2, 4, 1, 3, 5, 5, 1, 2, 5, 1, 5, 6, 2, 6, 2, 1, 1, 4, 4, 2, 1, 5, 2, 3, 4, 4, 2, 5, 2, 5, 5, 5, 3, 1, 4, 6, 6, 5, 5, 1, 3, 3, 6, 5, 6, 2, 1, 1, 1, 5, 4, 3, 1, 1, 2, 1, 3, 1, 6, 1, 5, 1, 2, 6, 1, 2, 2, 2, 5, 1, 2, 6, 5, 2, 2, 1, 3, 6, 2, 6, 1, 1, 6, 3, 5, 2, 6, 3, 1, 1, 4, 3, 5, 2, 3, 2, 4, 5, 1, 2, 5, 5, 3, 1, 4, 4, 5, 5, 5, 3, 2, 1, 3, 1, 3, 2, 2, 5, 6, 5, 2, 2, 6, 4, 1, 4, 2, 3, 5, 2, 5, 5, 6, 4, 6, 4, 3, 5, 4, 3, 5, 2, 6, 5, 5, 5, 1, 4, 6, 5, 3, 5, 5, 4, 5, 4, 4, 3, 1, 6, 4, 1, 5, 5, 6, 3, 2, 4, 3, 4, 4, 4, 5, 3, 3, 4, 2, 5, 3, 6, 1, 4, 5, 5, 6, 1, 4, 1, 1, 2, 2, 1, 3, 5, 1, 1, 2, 1, 2, 1, 1, 2, 3, 1, 3, 5, 5, 6, 4, 3, 5, 6, 1, 3, 4, 5, 6, 3, 1, 5, 4, 2, 3, 2, 6, 4, 6, 6, 6, 4, 1, 5, 1, 6, 4, 6, 6, 1, 3, 1, 4, 6, 2, 1, 4, 2, 4, 4, 4, 6, 6, 1, 6, 3, 6, 5, 3, 5, 3, 5, 6, 6, 1, 6, 2, 2, 1, 6, 1, 1, 1, 6, 3, 3, 6, 5, 1, 1, 4, 4, 2, 6, 5, 3, 4, 6, 5, 4, 2, 2, 1, 5, 3, 5, 5, 5, 4, 4, 6, 4, 4, 3, 1, 1, 6, 2, 3, 3, 3, 6, 5, 2, 1, 3, 1, 4, 4, 2, 3, 1, 2, 6, 1, 6, 4, 4, 3, 1, 4, 6, 1, 5, 4, 5, 3, 2, 2, 5, 5, 3, 6, 1, 2, 5, 5, 5, 4, 4, 4, 5, 1, 2, 4, 1, 6, 4, 3, 4, 4, 2, 4, 3, 4, 3, 5, 2, 3, 2, 3, 6, 6, 5, 2, 5, 5, 1, 4, 3, 5, 6, 5, 5, 1, 5, 3, 1, 3, 2, 2, 5, 6, 2, 6, 6, 4, 5, 2, 3, 1, 5, 5, 2, 2, 6, 5, 2, 5, 5, 3, 4, 6, 2, 1, 5, 6, 3, 2, 3, 2, 6, 6, 1, 4, 6, 5, 6, 6, 6, 5, 6, 6, 3, 6, 6, 4, 5, 5, 4, 6, 5, 2, 6, 6, 3, 3, 3, 3, 6, 5, 5, 6, 1, 3, 6, 1, 1, 1, 3, 4, 4, 4, 1, 4, 6, 6, 1, 5, 4, 1, 3, 2, 4, 5, 3, 6, 2, 3, 3, 3, 4, 5, 3, 2, 4, 4, 2, 3, 6, 4, 6, 6, 1, 2, 1, 2, 3, 1, 6, 5, 4, 4, 4, 5, 6, 6, 1, 5, 4, 1, 4, 3, 2, 1, 3, 2, 3, 1, 6, 5, 2, 3, 6, 6, 1, 1, 5, 2, 2, 1, 3, 2, 5, 2, 6, 5, 4, 5, 2, 2, 1, 6, 3, 1, 5, 1, 6, 4, 6, 3, 3, 1, 1, 4, 4, 6, 4, 2, 1, 6, 2, 4, 3, 4, 4, 5, 5, 2, 1, 4, 3, 5, 6, 6, 6, 2, 5, 6, 2, 3, 2, 4, 5, 6, 3, 3, 5, 4, 1, 1, 1, 2, 5, 6, 1, 4, 5, 3, 4, 2, 1, 6, 6, 4, 2, 6, 4, 3, 1, 2, 4, 2, 4, 6, 5, 6, 5, 6, 4, 4, 1, 2, 4, 6, 6, 1, 5, 2, 3, 2, 1, 5, 4, 1, 4, 4, 3, 6, 4, 6, 1, 6, 6, 2, 3, 5, 4, 3, 2, 1, 4, 3, 2, 2, 4, 5, 2, 1, 6, 1, 6, 1, 5, 1, 3, 1, 1, 6, 5, 6, 6, 6, 4, 6, 3, 4, 1, 4, 1, 1, 6, 3, 1, 3, 3, 4, 5, 2, 5, 4, 3, 1, 5, 5, 6, 6, 6, 5, 3, 6, 2, 4, 5, 6, 2, 3, 2, 3, 1, 3, 6, 4, 4, 5, 5, 5, 3, 1, 1, 6, 1, 4, 1, 2, 3, 3, 2, 3, 4, 1])
+
+test_seq = EmissionSequence(alph,[4, 6, 1, 1, 5, 2, 5, 3, 4, 2, 1, 6, 1, 5, 5, 5, 6, 2, 4, 3, 4, 3, 4, 1, 3, 4, 2, 2, 3, 3, 2, 6, 6, 3, 6, 4, 1, 4, 4, 4, 6, 2, 1, 1, 2, 2, 2, 3, 5, 1, 2, 1, 4, 2, 6, 1, 6, 4, 4, 1, 1, 4, 6, 5, 1, 2, 5, 6, 3, 5, 1, 1, 2, 2, 1, 1, 5, 4, 6, 6, 3, 5, 4, 4, 4, 3, 3, 6, 6, 2, 1, 2, 1, 3, 2, 6, 2, 4, 2, 4])
+
+test_vit = []
Index: /liacs/dbdm/dbdm_4/ecoli_hmm.py
===================================================================
--- /liacs/dbdm/dbdm_4/ecoli_hmm.py	(revision 60)
+++ /liacs/dbdm/dbdm_4/ecoli_hmm.py	(revision 60)
@@ -0,0 +1,62 @@
+#!/usr/bin/env python
+
+# http://en.wikipedia.org/wiki/Stop_codon
+# http://en.wikipedia.org/wiki/Escherichia_coli
+# http://en.wikipedia.org/wiki/Open_reading_frame
+
+import ghmm
+import sys
+
+
+def ecoli_hmm(seq):
+    """Try to find genes inside e sequence using a HMM"""
+    # Model 4 bases A C G T and unknown state N
+    sigma = ghmm.Alphabet(['a', 'c', 'g', 't', 'n' ])
+    print sigma
+   
+    # XXX: Proper values, based of statistics
+    # The transition matrix A is chosen such that it reflects the statistics
+    # Part one will try to get us into a gene, part two will tell us when to
+    # get out of it
+    A = [[0.6,0.4], [0.3, 0.7]]
+
+    # The emission probabilities matrix is modeled after the statistics with a
+    # total of 1
+    etogene = [ 0.2, 0.3, 0.3, 0.1, 0.2 ]
+    eingene = [ 0.2, 0.3, 0.3, 0.2, 0.1 ]
+
+    B = [etogene,eingene]
+    
+    # Initial distribution favors outside
+    pi  = [0.9, 0.1]
+
+    m = ghmm.HMMFromMatrices(sigma,ghmm.DiscreteDistribution(sigma),A ,B, pi)
+    print "Initial HMM", m
+
+    obs_seq = m.sampleSingle(20)
+    print "Observation sequence : ", obs_seq
+    obs = map(sigma.external, obs_seq)
+    print "Observations         : ", ''.join(obs)
+
+    # XXX: Training
+    
+    # XXX: Validation
+
+    # XXX: Results
+
+
+
+if __name__ == "__main__":
+    # Load data
+    try:
+        handle = open(sys.argv[1],"rU")
+    except IndexError:
+        print "Usage %s <data_file>" % (sys.argv[0])
+        sys.exit(64)
+    except IOError:
+        print "Unable to open '%s'" % (sys.argv[1])
+        sys.exit(64)
+    
+    seq = handle.read()
+    handle.close()
+    ecoli_hmm(seq)
Index: /liacs/dbdm/dbdm_4/unfair_casino_example.py
===================================================================
--- /liacs/dbdm/dbdm_4/unfair_casino_example.py	(revision 60)
+++ /liacs/dbdm/dbdm_4/unfair_casino_example.py	(revision 60)
@@ -0,0 +1,112 @@
+#!/usr/bin/env python
+#
+# Tutorial on using GHMM with Python
+# http://ghmm.sourceforge.net/ghmm-python-tutorial.html
+# Ported to python code for lazy programmers ;-)
+#
+# In the following, we assume that you have installed GHMM including the Python
+# bindings. Download the UnfairCasino.py-file.
+#
+# You might have seen the unfair casino example (Chair Biological Sequence
+# Analysis, Durbin et. al, 1998), where a dealer in a casino occasionally
+# exchanges a fair dice with a loaded one. We do not know if and when this
+# exchange happens, we only can observe the throws. This stochastic process we
+# will model with a HMM.
+
+import ghmm
+print "For aditional help use: help('ghmm')"
+
+#Creating the first model:
+
+# There are two states in our example. Either the dice is fair (state 0; Python
+# indexes arrays like C and C++ from 0) or it is loaded (state 1). We will
+# define the transition and emission matrices explicitly.
+
+# To save us some typing (namely ghmm. before any class or function from GHMM),
+# we will import all symbols from ghmm directly.
+
+from ghmm import *
+
+# Our alphabet are the numbers on the dice {1,2,3,4,5,6}. To be consistent with
+# Python's range function the upper limit is not part of a range.
+
+sigma = IntegerRange(1,7)
+
+# The transition matrix A is chosen such that we will stay longer in the fair
+# than the loaded state. The emission probabilities are set accordingly in B ---
+# guessing equal probabilities for the fair dice and higher probabilities for
+# ones and twos etc. for the loaded dice --- and the initial state distribution
+# pi has no preference for either state.
+
+A = [[0.9, 0.1], [0.3, 0.7]]
+efair = [1.0 / 6] * 6
+print "Fair probabilities:", efair
+
+eloaded = [3.0 / 13, 3.0 / 13, 2.0 / 13, 2.0 / 13, 2.0 / 13, 1.0 / 13]
+B = [efair, eloaded]
+pi = [0.5] * 2
+m = HMMFromMatrices(sigma, DiscreteDistribution(sigma), A, B, pi)
+print "Inital HHM", m
+
+# Generating sequences
+
+# At this time GHMM still exposes some internals. For example, symbols from
+# discrete alphabets are internally represented as 0,...,|alphabet| -1 and the
+# alphabet is responsible for translating between internal and external
+# representations.
+
+# First we sample from the model, specifying the length of the sequence, and
+# then transfer the observations to the external representation.
+
+obs_seq = m.sampleSingle(20)
+print "Observation sequence : ", obs_seq
+obs = map(sigma.external, obs_seq)
+print "Observations         : ", obs
+
+# Learning from data
+
+# In a real application we would want to train the model. That is, we would
+# like to estimate the parameters of the model from data with the goal to
+# maximize the likelihood of that data given the model. This is done for HMMs
+# with the Baum-Welch algorithm which is actually an instance of the well-known
+# Expectation-Maximization procedure for missing data problems in statistics.
+# The process is iterative and hence sometime called re-estimation.
+
+# We can use m as a starting point and use some 'real' data to train it. The
+# variable train_seq is an EmissionSequence or SequenceSet object.
+
+ 
+from UnfairCasino import train_seq 
+m.baumWelch(train_seq)
+print "Trained using baumWelch and 'train_seq'"
+print m
+
+# If you want more fine-grained control over the learning procedure, you can do
+# single steps and monitor the relevant diagnostics yourself, or employ
+# meta-heuristics such as noise-injection to avoid getting stuck in local
+# maxima. [Currently for continous emissions only] Computing a Viterbi-path
+
+# A Viterbi-path is a sequence of states Q maximizing the probability
+# P[Q|observations, model]. If we compute a Viterbi-path for a unfair casino
+# sequence of throws, the state sequence tells us when the casino was fair
+# (state 0) and when it was unfair (state 1).
+
+
+from UnfairCasino import test_seq
+v = m.viterbi(test_seq)
+print "fairness of test_seq: ", v
+
+# You can investigate how successful your model is at detecting cheating. Just
+# create artificial observations and see if you can detect them.
+
+my_seq = EmissionSequence(sigma, [1] * 20 + [6] * 10 + [1] * 40)
+print "artificial sequence    : ", my_seq
+print "artificial observation : ", m.viterbi(my_seq)
+
+print '''
+#
+# XXX: Some error on closing data, save to ignore as we are finished anyways
+#
+'''
+
+# EOF
