docs: def _forward(X, s=1.1, gamma=1., k=5):
      """
      Forward dynamic algorithm for burstness automaton HMM, from `Kleinberg
      (2002) <http://www.cs.cornell.edu/home/kleinber/bhs.pdf>`_.
  
      Parameters
      ----------
      X : list
          A series of time-gaps between events.
      s : float
          (default: 1.1) Scaling parameter ( > 1.)that controls graininess of
          burst detection. Lower values make the model more sensitive.
      gamma : float
          (default: 1.0) Parameter that controls the 'cost' of higher burst
          states. Higher values make it more 'difficult' to achieve a higher
          burst state.
      k : int
          (default: 5) Number of states. Higher values increase computational
          cost of the algorithm. A maximum of 25 is suggested by the literature.
  
      Returns
      -------
      states : list
          Optimal state sequence.
      """
      X = list(X)
  
      def alpha(i):
          return (n/T)*(s**i)
  
      def tau(i, j):
          if j > i:
              return (j-i)*gamma*log(n)
          return 0.
  
      def f(j, x):
          return alpha(j) * exp(-1. * alpha(j) * x)
  
      def C(j, t):
          if j == 0 and t == 0:
              return 0.
          elif t == 0:
              return float("inf")
  
          C_tau = min([C_values[l][t-1] + tau(l, j) for l in xrange(k)])
          return (-1. * log(f(j,X[t]))) + C_tau
  
      T = sum(X)
      n = len(X)
  
      # C() requires default (0) values, so we construct the "array" in advance.
      C_values = [[0 for t in xrange(len(X))] for j in xrange(k)]
      for j in xrange(k):
          for t in xrange(len(X)):
              C_values[j][t] = C(j,t)
  
      # Find the optimal state sequence.
      states = [argmin([c[t] for c in C_values]) for t in xrange(n)]
      return states