content: from __future__ import print_function
  
  import string
  import sys
  from collections import deque
  
  import numpy as np
  from scipy.misc import logsumexp
  from sklearn.base import BaseEstimator, _pprint
  from sklearn.utils import check_array, check_random_state
  from sklearn.utils.validation import check_is_fitted
  
  from . import _hmmc
  from .utils import normalize, log_normalize, iter_from_X_lengths, log_mask_zero
  
  
  #: Supported decoder algorithms.
  DECODER_ALGORITHMS = frozenset(("viterbi", "map"))
  
  
  class ConvergenceMonitor(object):
      """Monitors and reports convergence to :data:`sys.stderr`.
  
      Parameters
      ----------
      tol : double
          Convergence threshold. EM has converged either if the maximum
          number of iterations is reached or the log probability
          improvement between the two consecutive iterations is less
          than threshold.
  
      n_iter : int
          Maximum number of iterations to perform.
  
      verbose : bool
          If ``True`` then per-iteration convergence reports are printed,
          otherwise the monitor is mute.
  
      Attributes
      ----------
      history : deque
          The log probability of the data for the last two training
          iterations. If the values are not strictly increasing, the
          model did not converge.
  
      iter : int
          Number of iterations performed while training the model.
      """
      _template = "{iter:>10d} {logprob:>16.4f} {delta:>+16.4f}"
  
      def __init__(self, tol, n_iter, verbose):
          self.tol = tol
          self.n_iter = n_iter
          self.verbose = verbose
          self.history = deque(maxlen=2)
          self.iter = 0
  
      def __repr__(self):
          class_name = self.__class__.__name__
          params = dict(vars(self), history=list(self.history))
          return "{0}({1})".format(
              class_name, _pprint(params, offset=len(class_name)))
  
      def report(self, logprob):
          """Reports convergence to :data:`sys.stderr`.
  
          The output consists of three columns: iteration number, log
          probability of the data at the current iteration and convergence
          rate.  At the first iteration convergence rate is unknown and
          is thus denoted by NaN.
  
          Parameters
          ----------
          logprob : float
              The log probability of the data as computed by EM algorithm
              in the current iteration.
          """
          if self.verbose:
              delta = logprob - self.history[-1] if self.history else np.nan
              message = self._template.format(
                  iter=self.iter + 1, logprob=logprob, delta=delta)
              print(message, file=sys.stderr)
  
          self.history.append(logprob)
          self.iter += 1
  
      @property
      def converged(self):
          """``True`` if the EM algorithm converged and ``False`` otherwise."""
          # XXX we might want to check that ``logprob`` is non-decreasing.
          return (self.iter == self.n_iter or
                  (len(self.history) == 2 and
                   self.history[1] - self.history[0] < self.tol))
  
  
  class _BaseHMM(BaseEstimator):
      """Base class for Hidden Markov Models.
  
      This class allows for easy evaluation of, sampling from, and
      maximum-likelihood estimation of the parameters of a HMM.
  
      See the instance documentation for details specific to a
      particular object.
  
      Parameters
      ----------
      n_components : int
          Number of states in the model.
  
      startprob_prior : array, shape (n_components, )
          Initial state occupation prior distribution.
  
      transmat_prior : array, shape (n_components, n_components)
          Matrix of prior transition probabilities between states.
  
      algorithm : string
          Decoder algorithm. Must be one of "viterbi" or "map".
          Defaults to "viterbi".
  
      random_state: RandomState or an int seed
          A random number generator instance.
  
      n_iter : int, optional
          Maximum number of iterations to perform.
  
      tol : float, optional
          Convergence threshold. EM will stop if the gain in log-likelihood
          is below this value.
  
      verbose : bool, optional
          When ``True`` per-iteration convergence reports are printed
          to :data:`sys.stderr`. You can diagnose convergence via the
          :attr:`monitor_` attribute.
  
      params : string, optional
          Controls which parameters are updated in the training
          process.  Can contain any combination of 's' for startprob,
          't' for transmat, and other characters for subclass-specific
          emission parameters. Defaults to all parameters.
  
      init_params : string, optional
          Controls which parameters are initialized prior to
          training.  Can contain any combination of 's' for
          startprob, 't' for transmat, and other characters for
          subclass-specific emission parameters. Defaults to all
          parameters.
  
      Attributes
      ----------
      monitor\_ : ConvergenceMonitor
          Monitor object used to check the convergence of EM.
  
      startprob\_ : array, shape (n_components, )
          Initial state occupation distribution.
  
      transmat\_ : array, shape (n_components, n_components)
          Matrix of transition probabilities between states.
      """
      def __init__(self, n_components=1,
                   startprob_prior=1.0, transmat_prior=1.0,
                   algorithm="viterbi", random_state=None,
                   n_iter=10, tol=1e-2, verbose=False,
                   params=string.ascii_letters,
                   init_params=string.ascii_letters):
          self.n_components = n_components
          self.params = params
          self.init_params = init_params
          self.startprob_prior = startprob_prior
          self.transmat_prior = transmat_prior
          self.algorithm = algorithm
          self.random_state = random_state
          self.n_iter = n_iter
          self.tol = tol
          self.verbose = verbose
  
      def score_samples(self, X, lengths=None):
          """Compute the log probability under the model and compute posteriors.
  
          Parameters
          ----------
          X : array-like, shape (n_samples, n_features)
              Feature matrix of individual samples.
  
          lengths : array-like of integers, shape (n_sequences, ), optional
              Lengths of the individual sequences in ``X``. The sum of
              these should be ``n_samples``.
  
          Returns
          -------
          logprob : float
              Log likelihood of ``X``.
  
          posteriors : array, shape (n_samples, n_components)
              State-membership probabilities for each sample in ``X``.
  
          See Also
          --------
          score : Compute the log probability under the model.
          decode : Find most likely state sequence corresponding to ``X``.
          """
          check_is_fitted(self, "startprob_")
          self._check()
  
          X = check_array(X)
          n_samples = X.shape[0]
          logprob = 0
          posteriors = np.zeros((n_samples, self.n_components))
          for i, j in iter_from_X_lengths(X, lengths):
              framelogprob = self._compute_log_likelihood(X[i:j])
              logprobij, fwdlattice = self._do_forward_pass(framelogprob)
              logprob += logprobij
  
              bwdlattice = self._do_backward_pass(framelogprob)
              posteriors[i:j] = self._compute_posteriors(fwdlattice, bwdlattice)
          return logprob, posteriors
  
      def score(self, X, lengths=None):
          """Compute the log probability under the model.
  
          Parameters
          ----------
          X : array-like, shape (n_samples, n_features)
              Feature matrix of individual samples.
  
          lengths : array-like of integers, shape (n_sequences, ), optional
              Lengths of the individual sequences in ``X``. The sum of
              these should be ``n_samples``.
  
          Returns
          -------
          logprob : float
              Log likelihood of ``X``.
  
          See Also
          --------
          score_samples : Compute the log probability under the model and
              posteriors.
          decode : Find most likely state sequence corresponding to ``X``.
          """
          check_is_fitted(self, "startprob_")
          self._check()
  
          X = check_array(X)
          # XXX we can unroll forward pass for speed and memory efficiency.
          logprob = 0
          for i, j in iter_from_X_lengths(X, lengths):
              framelogprob = self._compute_log_likelihood(X[i:j])
              logprobij, _fwdlattice = self._do_forward_pass(framelogprob)
              logprob += logprobij
          return logprob
  
      def _decode_viterbi(self, X):
          framelogprob = self._compute_log_likelihood(X)
          return self._do_viterbi_pass(framelogprob)
  
      def _decode_map(self, X):
          _, posteriors = self.score_samples(X)
          logprob = np.max(posteriors, axis=1).sum()
          state_sequence = np.argmax(posteriors, axis=1)
          return logprob, state_sequence
  
      def decode(self, X, lengths=None, algorithm=None):
          """Find most likely state sequence corresponding to ``X``.
  
          Parameters
          ----------
          X : array-like, shape (n_samples, n_features)
              Feature matrix of individual samples.
  
          lengths : array-like of integers, shape (n_sequences, ), optional
              Lengths of the individual sequences in ``X``. The sum of
              these should be ``n_samples``.
  
          algorithm : string
              Decoder algorithm. Must be one of "viterbi" or "map".
              If not given, :attr:`decoder` is used.
  
          Returns
          -------
          logprob : float
              Log probability of the produced state sequence.
  
          state_sequence : array, shape (n_samples, )
              Labels for each sample from ``X`` obtained via a given
              decoder ``algorithm``.
  
          See Also
          --------
          score_samples : Compute the log probability under the model and
              posteriors.
          score : Compute the log probability under the model.
          """
          check_is_fitted(self, "startprob_")
          self._check()
  
          algorithm = algorithm or self.algorithm
          if algorithm not in DECODER_ALGORITHMS:
              raise ValueError("Unknown decoder {0!r}".format(algorithm))
  
          decoder = {
              "viterbi": self._decode_viterbi,
              "map": self._decode_map
          }[algorithm]
  
          X = check_array(X)
          n_samples = X.shape[0]
          logprob = 0
          state_sequence = np.empty(n_samples, dtype=int)
          for i, j in iter_from_X_lengths(X, lengths):
              # XXX decoder works on a single sample at a time!
              logprobij, state_sequenceij = decoder(X[i:j])
              logprob += logprobij
              state_sequence[i:j] = state_sequenceij
  
          return logprob, state_sequence
  
      def predict(self, X, lengths=None):
          """Find most likely state sequence corresponding to ``X``.
  
          Parameters
          ----------
          X : array-like, shape (n_samples, n_features)
              Feature matrix of individual samples.
  
          lengths : array-like of integers, shape (n_sequences, ), optional
              Lengths of the individual sequences in ``X``. The sum of
              these should be ``n_samples``.
  
          Returns
          -------
          state_sequence : array, shape (n_samples, )
              Labels for each sample from ``X``.
          """
          _, state_sequence = self.decode(X, lengths)
          return state_sequence
  
      def predict_proba(self, X, lengths=None):
          """Compute the posterior probability for each state in the model.
  
          X : array-like, shape (n_samples, n_features)
              Feature matrix of individual samples.
  
          lengths : array-like of integers, shape (n_sequences, ), optional
              Lengths of the individual sequences in ``X``. The sum of
              these should be ``n_samples``.
  
          Returns
          -------
          posteriors : array, shape (n_samples, n_components)
              State-membership probabilities for each sample from ``X``.
          """
          _, posteriors = self.score_samples(X, lengths)
          return posteriors
  
      def sample(self, n_samples=1, random_state=None):
          """Generate random samples from the model.
  
          Parameters
          ----------
          n_samples : int
              Number of samples to generate.
  
          random_state : RandomState or an int seed
              A random number generator instance. If ``None``, the object's
              ``random_state`` is used.
  
          Returns
          -------
          X : array, shape (n_samples, n_features)
              Feature matrix.
  
          state_sequence : array, shape (n_samples, )
              State sequence produced by the model.
          """
          check_is_fitted(self, "startprob_")
  
          if random_state is None:
              random_state = self.random_state
          random_state = check_random_state(random_state)
  
          startprob_cdf = np.cumsum(self.startprob_)
          transmat_cdf = np.cumsum(self.transmat_, axis=1)
  
          currstate = (startprob_cdf > random_state.rand()).argmax()
          state_sequence = [currstate]
          X = [self._generate_sample_from_state(
              currstate, random_state=random_state)]
  
          for t in range(n_samples - 1):
              currstate = (transmat_cdf[currstate] > random_state.rand()) \
                  .argmax()
              state_sequence.append(currstate)
              X.append(self._generate_sample_from_state(
                  currstate, random_state=random_state))
  
          return np.atleast_2d(X), np.array(state_sequence, dtype=int)
  
      def fit(self, X, lengths=None):
          """Estimate model parameters.
  
          An initialization step is performed before entering the
          EM algorithm. If you want to avoid this step for a subset of
          the parameters, pass proper ``init_params`` keyword argument
          to estimator's constructor.
  
          Parameters
          ----------
          X : array-like, shape (n_samples, n_features)
              Feature matrix of individual samples.
  
          lengths : array-like of integers, shape (n_sequences, )
              Lengths of the individual sequences in ``X``. The sum of
              these should be ``n_samples``.
  
          Returns
          -------
          self : object
              Returns self.
          """
          X = check_array(X)
          self._init(X, lengths=lengths)
          self._check()
  
          self.monitor_ = ConvergenceMonitor(self.tol, self.n_iter, self.verbose)
          for iter in range(self.n_iter):
              stats = self._initialize_sufficient_statistics()
              curr_logprob = 0
              for i, j in iter_from_X_lengths(X, lengths):
                  framelogprob = self._compute_log_likelihood(X[i:j])
                  logprob, fwdlattice = self._do_forward_pass(framelogprob)
                  curr_logprob += logprob
                  bwdlattice = self._do_backward_pass(framelogprob)
                  posteriors = self._compute_posteriors(fwdlattice, bwdlattice)
                  self._accumulate_sufficient_statistics(
                      stats, X[i:j], framelogprob, posteriors, fwdlattice,
                      bwdlattice)
  
              # XXX must be before convergence check, because otherwise
              #     there won't be any updates for the case ``n_iter=1``.
              self._do_mstep(stats)
  
              self.monitor_.report(curr_logprob)
              if self.monitor_.converged:
                  break
  
          return self
  
      def _do_viterbi_pass(self, framelogprob):
          n_samples, n_components = framelogprob.shape
          state_sequence, logprob = _hmmc._viterbi(
              n_samples, n_components, log_mask_zero(self.startprob_),
              log_mask_zero(self.transmat_), framelogprob)
          return logprob, state_sequence
  
      def _do_forward_pass(self, framelogprob):
          n_samples, n_components = framelogprob.shape
          fwdlattice = np.zeros((n_samples, n_components))
          _hmmc._forward(n_samples, n_components,
                         log_mask_zero(self.startprob_),
                         log_mask_zero(self.transmat_),
                         framelogprob, fwdlattice)
          return logsumexp(fwdlattice[-1]), fwdlattice
  
      def _do_backward_pass(self, framelogprob):
          n_samples, n_components = framelogprob.shape
          bwdlattice = np.zeros((n_samples, n_components))
          _hmmc._backward(n_samples, n_components,
                          log_mask_zero(self.startprob_),
                          log_mask_zero(self.transmat_),
                          framelogprob, bwdlattice)
          return bwdlattice
  
      def _compute_posteriors(self, fwdlattice, bwdlattice):
          # gamma is guaranteed to be correctly normalized by logprob at
          # all frames, unless we do approximate inference using pruning.
          # So, we will normalize each frame explicitly in case we
          # pruned too aggressively.
          log_gamma = fwdlattice + bwdlattice
          log_normalize(log_gamma, axis=1)
          return np.exp(log_gamma)
  
      def _init(self, X, lengths):
          """Initializes model parameters prior to fitting.
  
          Parameters
          ----------
          X : array-like, shape (n_samples, n_features)
              Feature matrix of individual samples.
  
          lengths : array-like of integers, shape (n_sequences, )
              Lengths of the individual sequences in ``X``. The sum of
              these should be ``n_samples``.
          """
          init = 1. / self.n_components
          if 's' in self.init_params or not hasattr(self, "startprob_"):
              self.startprob_ = np.full(self.n_components, init)
          if 't' in self.init_params or not hasattr(self, "transmat_"):
              self.transmat_ = np.full((self.n_components, self.n_components),
                                       init)
  
      def _check(self):
          """Validates model parameters prior to fitting.
  
          Raises
          ------
  
          ValueError
              If any of the parameters are invalid, e.g. if :attr:`startprob_`
              don't sum to 1.
          """
          self.startprob_ = np.asarray(self.startprob_)
          if len(self.startprob_) != self.n_components:
              raise ValueError("startprob_ must have length n_components")
          if not np.allclose(self.startprob_.sum(), 1.0):
              raise ValueError("startprob_ must sum to 1.0 (got {0:.4f})"
                               .format(self.startprob_.sum()))
  
          self.transmat_ = np.asarray(self.transmat_)
          if self.transmat_.shape != (self.n_components, self.n_components):
              raise ValueError(
                  "transmat_ must have shape (n_components, n_components)")
          if not np.allclose(self.transmat_.sum(axis=1), 1.0):
              raise ValueError("rows of transmat_ must sum to 1.0 (got {0})"
                               .format(self.transmat_.sum(axis=1)))
  
      def _compute_log_likelihood(self, X):
          """Computes per-component log probability under the model.
  
          Parameters
          ----------
          X : array-like, shape (n_samples, n_features)
              Feature matrix of individual samples.
  
          Returns
          -------
          logprob : array, shape (n_samples, n_components)
              Log probability of each sample in ``X`` for each of the
              model states.
          """
  
      def _generate_sample_from_state(self, state, random_state=None):
          """Generates a random sample from a given component.
  
          Parameters
          ----------
          state : int
              Index of the component to condition on.
  
          random_state: RandomState or an int seed
              A random number generator instance. If ``None``, the object's
              ``random_state`` is used.
  
          Returns
          -------
          X : array, shape (n_features, )
              A random sample from the emission distribution corresponding
              to a given component.
          """
  
      # Methods used by self.fit()
  
      def _initialize_sufficient_statistics(self):
          """Initializes sufficient statistics required for M-step.
  
          The method is *pure*, meaning that it doesn't change the state of
          the instance.  For extensibility computed statistics are stored
          in a dictionary.
  
          Returns
          -------
          nobs : int
              Number of samples in the data.
  
          start : array, shape (n_components, )
              An array where the i-th element corresponds to the posterior
              probability of the first sample being generated by the i-th
              state.
  
          trans : array, shape (n_components, n_components)
              An array where the (i, j)-th element corresponds to the
              posterior probability of transitioning between the i-th to j-th
              states.
          """
          stats = {'nobs': 0,
                   'start': np.zeros(self.n_components),
                   'trans': np.zeros((self.n_components, self.n_components))}
          return stats
  
      def _accumulate_sufficient_statistics(self, stats, X, framelogprob,
                                            posteriors, fwdlattice, bwdlattice):
          """Updates sufficient statistics from a given sample.
  
          Parameters
          ----------
          stats : dict
              Sufficient statistics as returned by
              :meth:`~base._BaseHMM._initialize_sufficient_statistics`.
  
          X : array, shape (n_samples, n_features)
              Sample sequence.
  
          framelogprob : array, shape (n_samples, n_components)
              Log-probabilities of each sample under each of the model states.
  
          posteriors : array, shape (n_samples, n_components)
              Posterior probabilities of each sample being generated by each
              of the model states.
  
          fwdlattice, bwdlattice : array, shape (n_samples, n_components)
              Log-forward and log-backward probabilities.
          """
          stats['nobs'] += 1
          if 's' in self.params:
              stats['start'] += posteriors[0]
          if 't' in self.params:
              n_samples, n_components = framelogprob.shape
              # when the sample is of length 1, it contains no transitions
              # so there is no reason to update our trans. matrix estimate
              if n_samples <= 1:
                  return
  
              lneta = np.zeros((n_samples - 1, n_components, n_components))
              _hmmc._compute_lneta(n_samples, n_components, fwdlattice,
                                   log_mask_zero(self.transmat_),
                                   bwdlattice, framelogprob, lneta)
              stats['trans'] += np.exp(logsumexp(lneta, axis=0))
  
      def _do_mstep(self, stats):
          """Performs the M-step of EM algorithm.
  
          Parameters
          ----------
          stats : dict
              Sufficient statistics updated from all available samples.
          """
          # The ``np.where`` calls guard against updating forbidden states
          # or transitions in e.g. a left-right HMM.
          if 's' in self.params:
              startprob_ = self.startprob_prior - 1.0 + stats['start']
              self.startprob_ = np.where(self.startprob_ == 0.0,
                                         self.startprob_, startprob_)
              normalize(self.startprob_)
          if 't' in self.params:
              transmat_ = self.transmat_prior - 1.0 + stats['trans']
              self.transmat_ = np.where(self.transmat_ == 0.0,
                                        self.transmat_, transmat_)
              normalize(self.transmat_, axis=1)
  