text: #!/usr/bin/env python
  
  import numpy
  import bisect
  
  Float = numpy.core.numerictypes.sctype2char(float)
  
  __author__ = "Peter Maxwell"
  __copyright__ = "Copyright 2007-2012, The Cogent Project"
  __credits__ = ["Peter Maxwell", "Gavin Huttley"]
  __license__ = "GPL"
  __version__ = "1.5.3"
  __maintainer__ = "Peter Maxwell"
  __email__ = "pm67nz@gmail.com"
  __status__ = "Production"
  
  class TransitionMatrix(object):
      """The transition matrix for a Markov process.  Just a square
      numpy array plus a list of state 'tags', eg:
      
      >>> a = numpy.array([
          [.9, .0, .1],
          [.0, .9, .1],
          [.5, .5, .0],], Float)
      >>> T = TransitionMatrix(a, ['x', 'y', 'm'])
      """
      
      def __init__(self, matrix, tags, stationary_probs=None):
          self.Matrix = numpy.array(matrix, Float)
          self.Tags = list(tags)
          self.size = len(matrix)
          assert matrix.shape == (self.size, self.size)
          assert len(tags) == self.size
          if stationary_probs is None:
              self._stationary_probs = None
          else:
              # Could recalculate, but if it is provided then faster to
              # just trust it
              assert len(stationary_probs) == self.size
              self._stationary_probs = numpy.array(
                      stationary_probs, Float)
      
      def _getStationaryProbs(self):
          if self._stationary_probs is None:
              matrix = self.Matrix
              for i in range(10):
                  matrix = numpy.core.multiarray.dot(matrix, matrix)
              self._stationary_probs = matrix[0]
          return self._stationary_probs
      StationaryProbs = property(_getStationaryProbs)
      
      def emit(self, random_series):
          """Generates an infinite sequence of states"""
          partitions = numpy.add.accumulate(self.Matrix, axis=1)
          for (state, row) in enumerate(partitions[:-1]):
              assert abs(row[-1]-1.0) < 1e-6, (state, self.Matrix[state])
          x = random_series.uniform(0.0, 1.0)
          state = bisect.bisect_left(
                  numpy.add.accumulate(self.StationaryProbs), x)
          while 1:
              yield self.Tags[state]
              x = random_series.uniform(0.0, 1.0)
              state = bisect.bisect_left(partitions[state], x)
      
      def __repr__(self):
          from cogent.util.table import Table
          labels = []
          for (i, label) in enumerate(self.Tags):
              if hasattr(label, '__len__') and not isinstance(
                      label, basestring):
                  label = ','.join(str(z) for z in label)
              # Table needs unique labels
              label = "%s (%s)" % (label, i)
              labels.append(label)
          heading = [''] + labels
          a = [[name] + list(row) for (name, row) in zip(labels, self.Matrix)]
          return str(Table(header = heading, rows = a))
      
      def withoutSilentStates(self):
          """An equivalent matrix without any of the states that have a
          false tag value"""
          N = self.size
          silent = numpy.array(
                  [(not max(tag)) for tag in self.Tags], Float)
          matrix = numpy.zeros([N, N], Float)
          for i in range(N):
              row = numpy.zeros([N], Float)
              emt = numpy.zeros([N], Float)
              row[i] = 1.0
              while max((row+emt)-emt):
                  row = numpy.dot(row, self.Matrix)
                  nul = silent * row
                  emt += row - nul
                  row = nul
              matrix[i] = emt
          keep = [i for i in range(self.size) if max(matrix[:,i])]
          matrix = numpy.take(matrix, keep, axis=0)
          matrix = numpy.take(matrix, keep, axis=1)
          tags = numpy.take(self.Tags, keep, axis=0)
          return type(self)(matrix, tags)
      
      def getLikelihoodOfSequence(self, obs, backward=False):
          """Just for testing really"""
          profile = numpy.zeros([len(obs), self.size], Float)
          for (i,a) in enumerate(obs):
              # This is suspiciously alphabet-like!
              profile[i, self.Tags.index(obs[i])] = 1.0
          return self.getLikelihoodOfProfile(profile, backward=backward)
      
      def getLikelihoodOfProfile(self, obs, backward=False):
          """Just for testing really"""
          if not backward:
              state_probs = self.StationaryProbs.copy()
              for i in range(len(obs)):
                  state_probs = numpy.dot(state_probs, self.Matrix) * obs[i]
              return sum(state_probs)
          else:
              state_probs = numpy.ones([self.size], Float)
              for i in range(len(obs)-1, -1, -1):
                  state_probs = numpy.dot(self.Matrix, state_probs * obs[i])
              return sum(state_probs * self.StationaryProbs)
      
      def getPosteriorProbs(self, obs):
          """'obs' is a sequence of state probability vectors"""
          result = numpy.zeros([len(obs), self.size], Float)
          
          # Forward
          state_probs = self.StationaryProbs.copy()
          for i in range(len(obs)):
              state_probs = numpy.dot(state_probs, self.Matrix) * obs[i]
              state_probs /= sum(state_probs)
              result[i] = state_probs
          
          # and Backward
          state_probs = numpy.ones([self.size], Float)
          for i in range(len(obs)-1, -1, -1):
              state_probs = numpy.dot(self.Matrix, state_probs)
              state_probs /= sum(state_probs)
              result[i] *= state_probs
              result[i] /= sum(result[i])
              state_probs *= obs[i]
          
          return result
      
      def nestTransitionMatricies(self, Ts, blended=None):
          """Useful for combining several X/Y/M Pair HMMs into
          one large HMM. The transition matricies 'Ts' end up
          along the diagonal blocks of the result, and the off
          diagonal values depend on the region switching probablities
          defined by 'self'"""
          
          if blended is None:
              blended = lambda a,b: (a+b)/2.0
              #blended = lambda a,b: numpy.sqrt(a*b)
              #blended = lambda a,b: b
          
          R = self.Matrix
          n = len(R)
          assert len(Ts) == n
          result = None
          for (x, a) in enumerate(self.Tags):
              a = Ts[a-1]
              if result is None:
                  tags = a.Tags
                  c = len(tags)
                  result = numpy.zeros([c*n, c*n], Float)
              else:
                  assert a.Tags == tags, (a.Tags, tags)
              a = a.Matrix
              for (y, b) in enumerate(self.Tags):
                  b = Ts[b-1].Matrix
                  block = self.Matrix[x,y] * blended(a, b)
                  result[x*c:(x+1)*c, y*c:(y+1)*c] = block
          all_tags = []
          for i in self.Tags:
              all_tags.extend([[i*e for e in tag] for tag in tags])
          return TransitionMatrix(result, all_tags)
      
  
  def SiteClassTransitionMatrix(switch, probs):
      """TM defined by stationary probabilities and a 'switch' parameter,
      switch=0.0 gives an identity Matrix, switch=1.0 gives independence,
      ie: zero-order markov process"""
      probs = numpy.asarray(probs)
      assert numpy.allclose(sum(probs), 1.0), probs
      I = numpy.identity(len(probs), Float)
      switch_probs =  (1.0 - I) * (probs * switch) + \
              I * (1.0 - (1.0 - probs) * switch)
      tags = [i+1 for i in range(len(switch_probs))]
      return TransitionMatrix(
              switch_probs, tags, stationary_probs=probs.copy())
      
  