base_variant.py

# base_variant.py is based on modification of the base.py file from the hmmlearn package

from __future__ import print_function

import string
import sys
import os
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

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 a posteriori 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, ), optional
		Parameters of the Dirichlet prior distribution for
		:attr:`startprob_`.

	transmat_prior : array, shape (n_components, n_components), optional
		Parameters of the Dirichlet prior distribution for each row
		of the transition probabilities :attr:`transmat_`.

	algorithm : string, optional
		Decoder algorithm. Must be one of "viterbi" or "map".
		Defaults to "viterbi".

	random_state: RandomState or an int seed, optional
		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, run_id=0,
				 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
		self.run_id = run_id

	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    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_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    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_")
		self._check()

		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.
		"""
		print("Initilization...")
		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)
				#if iter%3==0:
				#    filename = "log1a/stats_iter_%d_%d"%(self.run_id,iter)
				#    with open(filename,'a') as f_handle:
				#        np.savetxt(f_handle,posteriors, fmt='%.4f', delimiter='\t')
			
			# filename = "log1/stats_iter_%d_%d"%(self.run_id,iter)
			# np.savetxt(filename, posteriors, fmt='%.4f', delimiter='\t')

			# XXX must be before convergence check, because otherwise
			#     there won't be any updates for the case ``n_iter=1``.
			#print("Maximization...")
			#print(iter)
			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

			log_xi_sum = np.full((n_components, n_components), -np.inf)
			_hmmc._compute_log_xi_sum(n_samples, n_components, fwdlattice,
			                          log_mask_zero(self.transmat_),
			                          bwdlattice, framelogprob,
			                          log_xi_sum)
			stats['trans'] += np.exp(log_xi_sum)

	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_)
			#print 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)
			#print("updating transmat")
			#print(self.transmat_)