File Coverage

HMM.pd

Criterion	Covered	Total	%
statement	21	21	100.0
branch	8	16	50.0
condition			n/a
subroutine	2	2	100.0
pod	2	2	100.0
total	33	41	80.4

#-*- Mode: CPerl -*-

##======================================================================
## Header Administrivia
##======================================================================

our $VERSION = '0.06010';
pp_setversion($VERSION);

##-- accomodate typo-fix in PDL-2.008 API
eval "require PDL";
use version;
our $pdl_version       = version->parse($PDL::VERSION);
our $propagate_badflag = $pdl_version >= '2.008' ? "propagate_badflag" : "propogate_badflag";

##-- floating-point types
## PDL v2.082 doesn't like LD here (so we use 'E' instead)
##  - t/03_baum.t crashes at line 115 with `PP INTERNAL ERROR in logadd: unhandled datatype(11)`
my $FLOAT_TYPES = [qw(F D E)];

##------------------------------------------------------
## pm additions
pp_addpm({At=>'Top'},<<'EOPM');
=pod

=head1 NAME

PDL::HMM - Hidden Markov Model utilities in PDL

=head1 SYNOPSIS

 use PDL::HMM;

 ##-----------------------------------------------------
 ## Dimensions

 $N = $number_of_states;
 $M = $number_of_symbols;
 $T = $length_of_input;

 $A = $maximum_ambiguity;

 ##-----------------------------------------------------
 ## Parameters

 $af     = log(random($N,$N));
 $bf     = log(random($N,$M));
 $pif    = log(random($N));
 $omegaf = log(random($N));

 @theta = ($a,$b,$pi,$omega) = hmmmaximize($af,$bf,$pif,$omegaf);

 $o  = long(rint($M*random($T)));

 maximum_n_ind(dice_axis($a->logsumover+$pi+$b, 1,$o),
               ($oq=zeroes(long,$A,$T))); ##-- for constrained variants

 ##-----------------------------------------------------
 ## Log arithmetic

 $log0 = logzero;
 $logz = logadd(log($x),log($y));
 $logz = logdiff(log($x),log($y));
 $logz = logsumover(log($x));

 ##-----------------------------------------------------
 ## Sequence Probability

 $alpha  = hmmfw ($a,$b,$pi,    $o     ); ##-- forward (full)
 $alphaq = hmmfwq($a,$b,$pi,    $o, $oq); ##-- forward (constrained)

 $beta   = hmmbw ($a,$b,$omega,  $o    ); ##-- backward (full)
 $betaq  = hmmbwq($a,$b,$omega,  $o,$oq); ##-- backward (constrained)

 ##-----------------------------------------------------
 ## Parameter Estimation

 @expect = ($ea,$eb,$epi,$eomega) = hmmexpect0(@theta);    ##-- initialize

 hmmexpect (@theta, $o,     $alpha, $beta,  $ea,$eb,$epi); ##-- expect (full)
 hmmexpectq(@theta, $o,$oq, $alphaq,$betaq, $ea,$eb,$epi); ##-- expect (constrained)

 ($a,$b,$pi,$omega) = hmmmaximize($ea,$eb,$epi,$eomega);   ##-- maximize

 ##-----------------------------------------------------
 ## Sequence Analysis

 ($delta,$psi)   = hmmviterbi ($a,$b,$pi, $o);     ##-- trellis (full)
 ($deltaq,$psiq) = hmmviterbiq($a,$b,$pi, $o,$oq); ##-- trellis (constrained)

 $paths  = hmmpath (     $psi,  sequence($N));     ##-- backtrace (full)
 $pathsq = hmmpathq($oq, $psiq, sequence($A));     ##-- backtrace (constrained)

=cut

EOPM
## /pm additions
##------------------------------------------------------

##------------------------------------------------------
## Exports: None
#pp_export_nothing();

##------------------------------------------------------
## Includes / defines
pp_addhdr(<<'EOH');

#include

/*#define DEBUG_ALPHA*/
/*#define DEBUG_BETA*/
/*#define DEBUG_VITERBI*/

EOH


##======================================================================
## C Utilities
##======================================================================

##----------------------------------------------------------------------
## Log addition
pp_addhdr(<<'EOH');

/* logadd(x,y) = log(exp(x)+exp(y))
 *   + Code from Manning & Schütze (1997), Sec. 9.4, page 337
 *
 * LOG_BIG = log(1E31)
 */
#define LOG_BIG    71.3801378828154
#define LOG_ZERO  -1E+38
#define LOG_ONE    0
#define LOG_NONE   1
static inline double logadd1(double x, double y) {
  if      (y-x > LOG_BIG) return y;
  else if (x-y > LOG_BIG) return x;
  /*else        return min(x,y) + log(exp(x-min(x,y)) + exp(y-min(x,y))); */
  else if (x
  else          return y        + log(exp(x-y)        + 1);
}
static inline double logadd0(double x, double y) {
  return log(exp(x)+exp(y));
}

/* logdiff(x,y) = log(exp(x)-exp(y))
 *   + adapted from above
 *   + always returns positive (i.e. symmetric difference)
 */
static inline double logdiff1(double x, double y) {
  if      (y-x > LOG_BIG) { return y; }
  else if (x-y > LOG_BIG) { return x; }
  /*else        { return max(x,y) + log(exp(max(x,y)-max(x,y)) - exp(min(x,y)-max(x,y))); } */
  /*                   = max(x,y) + log( 1 - exp(min(x,y)-max(x,y))); } */
  else if (x>y) { return x        + log( 1 - exp(y-x)); }
  else          { return y        + log( 1 - exp(x-y)); }
}
static inline double logdiff0(double x, double y) {
  return log(x>y ? (exp(x)-exp(y)) : (exp(y)-exp(x)));
}

/*
#define logadd(x,y) logadd0(x,y)
#define logdiff(x,y) logdiff0(x,y)
*/

#define logadd(x,y) logadd1(x,y)
#define logdiff(x,y) logdiff1(x,y)

EOH


##======================================================================
## PDL::PP Wrappers
##======================================================================

##======================================================================
## Basic Utilities
pp_addpm(<<'EOPM');
=pod

=head1 Log Arithmetic

=cut
EOPM

##------------------------------------------------------
## logzero(): near approximation of log(0)
pp_def('logzero',
       Pars => '[o]a()',
       GenericTypes => $FLOAT_TYPES,
       Code => 'broadcastloop %{ $a() = LOG_ZERO; %} $PDLSTATESETGOOD(a);',
       HandleBad=>1,
       Doc => 'Approximates $a() = log(0), avoids nan.',
       BadDoc => 'logzero() handles bad values.  The state of the output PDL is always good.',
      );


##------------------------------------------------------
## log addition: logadd(a,b) = log(exp(a)+exp(b))
pp_def('logadd',
       Pars => 'a(); b(); [o]c()',
       GenericTypes => $FLOAT_TYPES,
       Inplace=>['a'], ##-- can run inplace on a()
       Code => '$c() = logadd($a(),$b());',
       Doc => 'Computes $c() = log(exp($a()) + exp($b())), should be more stable.',
      );

##------------------------------------------------------
## log subtraction: logdiff(a,b) = log(exp(max(a,b))-exp(min(a,b)))
pp_def('logdiff',
       Pars => 'a(); b(); [o]c()',
       GenericTypes => $FLOAT_TYPES,
       Inplace=>['a'], ##-- can run inplace on a()
       Code => '$c() = logdiff($a(),$b());',
       Doc => 'Computes log symmetric difference c = log(exp(max(a,b)) - exp(min(a,b))), may be more stable.',
      );


##------------------------------------------------------
## log sum: logsumover(a) = log(sumover(exp(a)))
pp_def('logsumover',
       Pars => 'a(n); [o]b()',
       GenericTypes => $FLOAT_TYPES,
       Code => (join(" ",
		     'double sum=LOG_ZERO;',
		     'loop (n) %{ sum = logadd($a(),sum); %}',
		     '$b() = sum;')),
       Doc => 'Computes $b() = log(sumover(exp($a()))), should be more stable.',
      );


##======================================================================
## Sequence Probability
pp_addpm(<<'EOPM');
=pod

=head1 Sequence Probability

=cut
EOPM


##------------------------------------------------------
## Forward probability: hmmfw(A,B,pi, O, [o]alpha)
pp_def
('hmmfw',
 Pars => 'a(N,N); b(N,M); pi(N);  o(T);  [o]alpha(N,T)',
 GenericTypes => $FLOAT_TYPES,
 Code =>
('
 /*-- Initialize: t==0 --*/
 int i,j,t, o_tp1 = $o(T=>0);
 loop (N) %{
   $alpha(T=>0) = $pi() + $b(M=>o_tp1);

   #ifdef DEBUG_ALPHA
   printf("INIT: j=%u,t=0,o=%d:  pi(j=%u)=%.2e  b(j=%u,o=%d)=%.2e  alpha(j=%u,t=0)=%.2e\n",
	  n,o_tp1,
	  n,      exp($pi()),
	  n,o_tp1 exp($b(M=>o_tp1)),
	  n,      exp($alpha(T=>0))
	 );
   #endif
 %}

 #ifdef DEBUG_ALPHA
 printf("\n\n");
 #endif

 /*-- Loop: time t>0 --*/
 for (t=0; t < $SIZE(T)-1; t++) {
   o_tp1 = $o(T=>t+1);

   /*-- Loop: state_(t+1)==j --*/
   for (j=0; j<$SIZE(N); j++) {
     $GENERIC(alpha) alpha_j_tp1 = ($GENERIC(alpha))LOG_ZERO;


     /*-- Loop: state_t==i --*/
     for (i=0; i<$SIZE(N); i++) {
       alpha_j_tp1 = logadd( $alpha(N=>i,T=>t) + $a(N0=>i,N1=>j),  alpha_j_tp1 );

       #ifdef DEBUG_ALPHA
        printf("i=%u,j=%u,t=%u,o=%d:  alpha(i=%u,t=%u)=%.2e  a(i=%u,j=%u)=%.2e  b(j=%u,o=%d)=%.2e  prod=%.2e  sum=%.2e\n",
	       i,j,t,o_tp1,
	       i,t,     exp($alpha(N=>i,T=>t)),
	       i,j,     exp($a(N0=>i,N1=>j)),
	       j,o_tp1, exp($b(N=>j,M=>o_tp1)),
	       exp( $alpha(N=>i,T=>t) + $a(N0=>i,N1=>j) ), exp(alpha_j_tp1));
       #endif
     }

     /*-- Storage: alpha(time=t+1, state=j) --*/
     $alpha(N=>j,T=>t+1) = alpha_j_tp1 + $b(N=>j,M=>o_tp1);

     #ifdef DEBUG_ALPHA
      printf("----> alpha(j=%u,t=%u)=%.2E\n", j,t+1, exp($alpha(N=>j,T=>t+1)));
     #endif
   }
   #ifdef DEBUG_ALPHA
    printf("\n\n");
   #endif
 }
'),

Doc=>
  ('Compute forward probability (alpha) matrix
for input $o given model parameters
@theta = ($a, $b, $pi, $omega).

Output (pseudocode) for all 0<=i

 $alpha(i,t) = log P( $o(0:t), q(t)==i | @theta )

Note that the final-state probability vector $omega() is neither
passed to this function nor used in the computation, but
can be used to compute the final sequence probability for $o as:

  log P( $o | @theta ) = logsumover( $omega() + $alpha(:,t-1) )

'),

);

pp_addpm('*hmmalpha = \&hmmfw;');
pp_add_exported('','hmmalpha');


##------------------------------------------------------
## Forward probability (constrained): hmmfwq(A,B,pi, O,Q, [o]alphaq)
pp_def
('hmmfwq',
 Pars => 'a(N,N); b(N,M); pi(N);  o(T); oq(Q,T);  [o]alphaq(Q,T)',
 GenericTypes => $FLOAT_TYPES,
 Code =>
('
 /*-- Initialize: t==0 --*/
 int i,j,t, o_tp1 = $o(T=>0);
 int qi,qj;

 for (qi=0; qi < $SIZE(Q); qi++) {
   j = $oq(Q=>qi,T=>0);
   $alphaq(Q=>qi,T=>0) = (j>=0 ? $pi(N=>j) + $b(N=>j,M=>o_tp1) : ($GENERIC(alphaq))LOG_ZERO);
 }

 #ifdef DEBUG_ALPHA
 printf("\n\n");
 #endif

 /*-- Loop: time t>0 --*/
 for (t=0; t < $SIZE(T)-1; t++) {
   o_tp1 = $o(T=>t+1);

   /*-- Loop: q_(t+1)=qj : state_(t+1)=oq(qj,t+1)=j --*/
   for (qj=0; qj < $SIZE(Q); qj++) {
     $GENERIC(alphaq) alpha_j_tp1 = ($GENERIC(alphaq))LOG_ZERO;
     j = $oq(Q=>qj,T=>t+1);

     /*-- Loop: q_(t)=qi : state_(t)=oq(qi,t)=i --*/
     for (qi=0; j>=0 && qi < $SIZE(Q); qi++) {
       i = $oq(Q=>qi,T=>t);
       if (i < 0) break;

       alpha_j_tp1 = logadd( $alphaq(Q=>qi,T=>t) + $a(N0=>i,N1=>j),  alpha_j_tp1 );

       #ifdef DEBUG_ALPHA
        printf("qi=%u,i=%d,qj=%u,j=%d,t=%u,o=%d:  alphaq(qi=%u,t=%u)=%.2e  a(i=%d,j=%d)=%.2e  b(j=%d,o=%d)=%.2e  prod=%.2e  sum=%.2e\n",
	       qi,i, qj,j, t,o_tp1,
	       i,t,     exp($alphaq(Q=>qi,T=>t)),
	       i,j,     ((i>=0 && j>=0) ? exp($a(N0=>i,N1=>j)) : 0),
	       j,o_tp1, ((i>=0 && j>=0) ? exp($b(N=>j,M=>o_tp1)) : 0),
               ((i>=0 && j>=0) ? exp( $alphaq(Q=>qi,T=>t) + $a(N0=>i,N1=>j) ) : 0),
               exp(alpha_j_tp1));
       #endif
     }
     /*-- End Loop: q_(t)=qi : state_(t)=oq(qi,t)=i --*/

     /*-- Storage: alphaq(time=t+1, stateIndex=qj) --*/
     if (j>=0) {
       $alphaq(Q=>qj,T=>t+1) = alpha_j_tp1 + $b(N=>j,M=>o_tp1);
     } else {
       $alphaq(Q=>qj,T=>t+1) = ($GENERIC(alphaq))LOG_ZERO;
     }

     #ifdef DEBUG_ALPHA
      printf("----> alphaq(qj=%u [j=%d], t=%u)=%.2E\n", qj,j,t+1, exp($alphaq(Q=>qj,T=>t+1)));
     #endif
   }
   /*-- End Loop: q_(t+1)=qj : state_(t+1)=oq(qj,t+1)=j --*/

   #ifdef DEBUG_ALPHA
    printf("\n\n");
   #endif
 }
 /*-- End Loop: time t>0 --*/
'),

Doc=>
  ('Compute constrained forward probability (alphaq) matrix
for input $o given model parameters
@theta = ($a, $b, $pi, $omega),
considering only the initial
non-negative state indices in $oq(:,t) for each observation $o(t).

Output (pseudocode) for all 0<=qi

 $alphaq(qi,t) = log P( $o(0:t), q(t)==$oq(qi,t) | @theta )

Note that the final-state probability vector $omega() is neither
passed to this function nor used in the computation, but
can be used to compute the final sequence probability for $o as:

  log P( $o | @theta ) = logsumover( $alphaq(:,t-1) + $omega($oqTi) )

where:

  $oqTi = $oq(:,t-1)->where($oq(:,t-1)>=0)

'),

);

pp_addpm('*hmmalphaq = \&hmmfwq;');
pp_add_exported('','hmmalphaq');


##------------------------------------------------------
## Backward probability: hmmbw(A,B,omega, O, [o]beta)
pp_def
#@l=
('hmmbw',
 GenericTypes => $FLOAT_TYPES,
 Pars => 'a(N,N); b(N,M); omega(N); o(T); [o]beta(N,T)',
 Code =>
('
 int i,j,t = $SIZE(T)-1;

 /*-- Initialize: time t==T --*/
 loop(N) %{ $beta(T=>t) = $omega(); %}

 /*-- Loop: time t < T --*/
 for (t--; t >= 0; t--) {
   int o_tp1 = $o(T=>t+1);

   /*-- Loop: t
   for (i=0; i<$SIZE(N); i++) {
     $GENERIC(beta) beta_i_t = ($GENERIC(beta))LOG_ZERO;

     /*-- Loop: t
     for (j=0; j<$SIZE(N); j++) {
       beta_i_t = logadd( $a(N0=>i,N1=>j) + $b(N=>j,M=>o_tp1) + $beta(N=>j,T=>t+1) ,  beta_i_t );

       #ifdef DEBUG_BETA
       printf("i=%u,j=%u,t=%u,o=%d:  a(i=%u,j=%u)=%.2e  b(j=%u,o=%u)=%.2e  beta(j=%u,t+1=%u)=%.2e   prod=%.2e  sum=%.2e\n",
	      i,j,t,o_tp1,
	      i,j,    exp($a(N0=>i,N1=>j)),
	      j,o_t,  exp($b(N=>j,M=>o_t)),
	      j,t+1,  exp($beta(N=>j,T=>t+1)),
	      exp($a(N0=>i,N1=>j)+$b(N=>j,M=>o_t)+$beta(N=>j,T=>t+1)),  exp(beta_i_t));
       #endif
     }

     /*-- t
     $beta(N=>i,T=>t) = beta_i_t;

     #ifdef DEBUG_BETA
     printf("\n");
     #endif
   }
   #ifdef DEBUG_BETA
   printf("\n\n");
   #endif
 }
'
),

Doc=>
  ('Compute backward probability (beta) matrix
for input $o given model parameters
@theta = ($a, $b, $pi, $omega).

Output (pseudocode) for all 0<=i

 $beta(i,t) = log P( $o(t+1:T-1) | q(t)==i, @theta )

Note that the initial-state probability vector $pi() is neither
passed to this function nor used in the computation, but
can be used to compute the final sequence probability for $o as:

  log P( $o | @theta ) = logsumover( $pi() + $b(:,$o(0)) + $beta(:,0) )

'),

);


pp_addpm('*hmmbeta = \&hmmbw;');
pp_add_exported('','hmmbeta');


##------------------------------------------------------
## Backward probability (constrained): hmmbwq(A,B,omega, O,Q, [o]beta)
pp_def
#@l=
('hmmbwq',
 GenericTypes => $FLOAT_TYPES,
 Pars => 'a(N,N); b(N,M); omega(N); o(T); oq(Q,T); [o]betaq(Q,T)',
 Code =>
('
 int i,j,t = $SIZE(T)-1;
 int qi, qj;

 /*-- Initialize: time t==T --*/
 for (qi=0; qi < $SIZE(Q); qi++) {
   i = $oq(Q=>qi,T=>t);
   $betaq(Q=>qi,T=>t) = (i>=0 ? $omega(N=>i) : ($GENERIC(betaq))LOG_ZERO);
 }

 /*-- Loop: time t < T --*/
 for (t--; t >= 0; t--) {
   int o_tp1 = $o(T=>t+1);

   /*-- Loop: t
   for (qi=0; qi<$SIZE(Q); qi++) {
     $GENERIC(betaq) beta_i_t = ($GENERIC(betaq))LOG_ZERO;
     i = $oq(Q=>qi,T=>t);

     /*-- Loop: t
     for (qj=0; i>=0 && qj<$SIZE(Q); qj++) {
       j = $oq(Q=>qj,T=>t+1);
       if (j < 0) break;

       beta_i_t = logadd( $a(N0=>i,N1=>j) + $b(N=>j,M=>o_tp1) + $betaq(Q=>qj,T=>t+1) ,  beta_i_t );
     }

     /*-- t
     $betaq(Q=>qi,T=>t) = beta_i_t;
   }
   /*-- End Loop: t
 }
 /*-- End Loop: time t < T --*/
'
),

Doc=>
  ('Compute constrained backward probability (betaq) matrix
for input $o given model parameters
@theta = ($a, $b, $pi, $omega),
considering only the initial non-negative state indices in $oq(:,t) for
each observation $o(t).

Output (pseudocode) for all 0<=qi

 $betaq(qi,t) = log P( $o(t+1:T-1) | q(t)==$oq(qi,t), @theta )

Note that the initial-state probability vector $pi() is neither
passed to this function nor used in the computation, but
can be used to compute the final sequence probability for $o as:

  log P( $o | @theta ) = logsumover( $betaq(:,0) + $pi($oq0i) + $b($oq0i,$o(0)) )

where:

  $oq0i = $oq(:,0)->where( $oq(:,0) >= 0 );

'),

);


pp_addpm('*hmmbetaq = \&hmmbwq;');
pp_add_exported('','hmmbetaq');


##======================================================================
## Parameter Estimation

##------------------------------------------------------
## Parameter Estimation: Initialize
pp_addpm(<<'EOPM');
=pod

=head1 Parameter Estimation

=head2 hmmexpect0

=for sig

  Signature: (a(N,N); b(N,M); pi(N); omega(N); [o]ea(N,N); [o]eb(N,M); [o]epi(N); [o]eomega(N))

Initializes expectation matrices $ea(), $eb() and $epi() to logzero().
For use with hmmexpect().

=cut

sub hmmexpect0 {
  my ($a,$b,$pi,$omega, $ea,$eb,$epi,$eomega) = @_;

  $ea  = zeroes($a->type,  $a->dims)  if (!defined($ea));
  $eb  = zeroes($b->type,  $b->dims)  if (!defined($eb));
  $epi = zeroes($pi->type, $pi->dims) if (!defined($epi));
  $eomega = zeroes($omega->type, $omega->dims) if (!defined($eomega));

  $ea  .= PDL::logzero();
  $eb  .= PDL::logzero();
  $epi .= PDL::logzero();
  $eomega .= PDL::logzero();

  return ($ea,$eb,$epi,$eomega);
}


EOPM

pp_add_exported('', 'hmmexpect0');



##------------------------------------------------------
## Parameter Estimation: Expect
pp_def
('hmmexpect',
 GenericTypes => $FLOAT_TYPES,
 Pars => join(" ",
              qw(a(N,N);
                 b(N,M);
                 pi(N);
                 omega(N);
                 o(T);
                 alpha(N,T);
                 beta(N,T);),
              qw([o]ea(N,N);
                 [o]eb(N,M);
                 [o]epi(N);
                 [o]eomega(N);)),
 Code =>
('
 int i,j,t;
 int o_tp1, o_t;
 double p_o = LOG_ZERO;
 double gamma_it;
 double xi_ijt;

 /*-- Initialize: t==(T-1): P(o|@theta) --*/
 t = $SIZE(T)-1;
 loop (N) %{ p_o = logadd(p_o, $omega() + $alpha(T=>t)); %}


 /*-- Initialize: t==(T-1): Iterate: state_t==i: get gamma(i,t) --*/
 o_t = $o(T=>t);
 for (i=0; i<$SIZE(N); i++) {
   gamma_it         = $alpha(N=>i,T=>t) + $beta(N=>i,T=>t) - p_o;
   $eb(N=>i,M=>o_t) = logadd($eb(N=>i,M=>o_t), gamma_it);
   $eomega(N=>i)    = logadd($eomega(N=>i)   , gamma_it);
 }

 /*-- Main: Iterate: T-1 > t >= 0 --*/
 for (t--; t>=0; t--) {
   o_tp1 = o_t;
   o_t   = $o(T=>t);

   /*-- Main: Iterate: state_t == i --*/
   for (i=0; i<$SIZE(N); i++) {
     gamma_it = $alpha(N=>i,T=>t) + $beta(N=>i,T=>t) - p_o;

     /*-- Main: Iterate: state_(t+1) == j --*/
     for (j=0; j<$SIZE(N); j++) {
       xi_ijt = $alpha(N=>i,T=>t) + $a(N0=>i,N1=>j) + $b(N=>j,M=>o_tp1) + $beta(N=>j,T=>t+1) - p_o;

       $ea(N0=>i,N1=>j) = logadd(xi_ijt, $ea(N0=>i,N1=>j));
     }

     /*-- Main: Update: pi --*/
     if (t==0) $epi(N=>i) = logadd(gamma_it, $epi(N=>i));

     /*-- Main: Update: b --*/
     $eb(N=>i,M=>o_t)     = logadd(gamma_it, $eb(N=>i,M=>o_t));
   }
 }
'),

 Doc =>
  ('Compute partial Baum-Welch re-estimation of the model @theta = ($a, $b, $pi, $omega)
for the observation sequence $o() with forward- and backward-probability
matrices $alpha(), $beta().  Result is recorded as log pseudo-frequencies
in the expectation matrices $ea(), $eb(), $epi(), and $eomega(), which are required parameters,
and should have been initialized (e.g. by L()) before calling this function.

Can safely be called sequentially for incremental reestimation.
'),
);


##------------------------------------------------------
## Parameter Estimation: Expect (constrained)
pp_def
('hmmexpectq',
 GenericTypes => $FLOAT_TYPES,
 Pars => join(" ",
              qw(a(N,N);
                 b(N,M);
                 pi(N);
                 omega(N);),
              qw(o(T);
                 oq(Q,T);),
              qw(alphaq(Q,T);
                 betaq(Q,T);),
              qw([o]ea(N,N);
                 [o]eb(N,M);
                 [o]epi(N);
                 [o]eomega(N);)),
 Code =>
('
 int i,j,t, qi,qj;
 int o_tp1, o_t;
 double p_o = LOG_ZERO;
 double gamma_it;
 double xi_ijt;

 /*-- Initialize: t==(T-1): P(o|@theta) --*/
 t = $SIZE(T)-1;
 for (qi=0; qi < $SIZE(Q); qi++) {
   i = $oq(Q=>qi,T=>t);
   if (i < 0) break;

   p_o = logadd(p_o, $omega(N=>i) + $alphaq(Q=>qi,T=>t));
 }

 /*-- Initialize: t==(T-1): Iterate: q_(t)=qi: state_t=oq(qi,t)=i: get gamma(i,t) --*/
 o_t = $o(T=>t);
 for (qi=0; qi < $SIZE(Q); qi++) {
   i                = $oq(Q=>qi,T=>t);
   if (i < 0) break;
   gamma_it         = $alphaq(Q=>qi,T=>t) + $betaq(Q=>qi,T=>t) - p_o;
   $eb(N=>i,M=>o_t) = logadd($eb(N=>i,M=>o_t), gamma_it);
   $eomega(N=>i)    = logadd($eomega(N=>i)   , gamma_it);
 }

 /*-- Loop: T-1 > t >= 0 --*/
 for (t--; t>=0; t--) {
   o_tp1 = o_t;
   o_t   = $o(T=>t);

   /*-- Loop: q_(t)=qi: state_(t)=oq(qi,t)=i --*/
   for (qi=0; qi<$SIZE(Q); qi++) {
     i        = $oq(Q=>qi,T=>t);
     if (i < 0) break;
     gamma_it = $alphaq(Q=>qi,T=>t) + $betaq(Q=>qi,T=>t) - p_o;

     /*-- Loop: q_(t+1)=qj: state_(t+1)=oq(qj,t+1)=j --*/
     for (qj=0; qj<$SIZE(Q); qj++) {
       j      = $oq(Q=>qj,T=>t+1);
       if (j < 0) break;

       xi_ijt = $alphaq(Q=>qi,T=>t) + $a(N0=>i,N1=>j) + $b(N=>j,M=>o_tp1) + $betaq(Q=>qj,T=>t+1) - p_o;

       $ea(N0=>i,N1=>j) = logadd(xi_ijt, $ea(N0=>i,N1=>j));
     }

     /*-- Update: pi --*/
     if (t==0) $epi(N=>i) = logadd(gamma_it, $epi(N=>i));

     /*-- Update: b --*/
     $eb(N=>i,M=>o_t)     = logadd(gamma_it, $eb(N=>i,M=>o_t));
   }
   /*-- End Loop: q_(t)=qi: state_(t)=oq(qi,t)=i --*/
 }
 /*-- End Loop: T-1 > t >= 0 --*/
'),

 Doc =>
  ('Compute constrained partial Baum-Welch re-estimation of the model @theta = ($a, $b, $pi, $omega)
for the observation sequence $o(),
with constrained forward- and backward-probability
matrices $alphaq(), $betaq(),
considering only the initial non-negative state
indices in $oq(:,t) for observation $o(t).
Result is recorded as log pseudo-frequencies
in the expectation matrices $ea(), $eb(), $epi(), and $eomega(), which are required parameters,
and should have been initialized (e.g. by L()) before calling this function.

Can safely be called sequentially for incremental reestimation.
'),
);


##------------------------------------------------------
## Parameter Estimation: Maximization
pp_addpm(<<'EOPM');
=pod

=head2 hmmmaximize

=for sig

  Signature: (Ea(N,N); Eb(N,M); Epi(N); Eomega(N); [o]ahat(N,N); [o]bhat(N,M); [o]pihat(N); [o]omegahat(N));

Maximizes expectation values from $Ea(), $Eb(), $Epi(), and $Eomega()
to log-probability matrices $ahat(), $bhat(), $pihat(), and $omegahat().
Can also be used to compile a maximum-likelihood model
from log-frequency matrices.

=cut

sub hmmmaximize {
  my ($ea,$eb,$epi,$eomega, $ahat,$bhat,$pihat,$omegahat) = @_;

  $ahat  = zeroes($ea->type,  $ea->dims)  if (!defined($ahat));
  $bhat  = zeroes($eb->type,  $eb->dims)  if (!defined($bhat));
  $pihat = zeroes($epi->type, $epi->dims) if (!defined($pihat));
  $omegahat = zeroes($eomega->type, $eomega->dims) if (!defined($omegahat));

  my $easumover = $ea->xchg(0,1)->logsumover->inplace->logadd($eomega);

  $ahat  .= $ea  - $easumover;
  $bhat  .= $eb  - $eb->xchg(0,1)->logsumover;
  $pihat .= $epi - $epi->logsumover;
  $omegahat .= $eomega - $easumover;

  return ($ahat,$bhat,$pihat,$omegahat);
}

EOPM

pp_add_exported('', 'hmmmaximize');


##======================================================================
## Sequence Analysis
pp_addpm(<<'EOPM');
=pod

=head1 Sequence Analysis

=cut
EOPM


##--------------------------------------------------------------
## Sequence Analysis: Viterbi
pp_def
('hmmviterbi',
 Pars => join(" ",
              qw(a(N,N);
                 b(N,M);
                 pi(N);),
              #qw(omega(N);),
              qw(o(T);
                 [o]delta(N,T);),
              'int [o]psi(N,T)'),
 GenericTypes => $FLOAT_TYPES,
 Code =>
('
 int i,j, t, o_t;
 double delta_jt, delta_tmp;
 int psi_jt;

 /*-- Initialize: t==0: Loop: state_0==N --*/
 o_t = $o(T=>0);
 loop (N) %{
   $delta(T=>0) = $pi() + $b(M=>o_t);
   $psi  (T=>0) = 0;
#ifdef DEBUG_VITERBI
   printf("t=0,j=%d,o_t=%d: delta(t=0,j=%d)=%.2e  psi(t=0,j=%d)=%.0g  b(j=%d,o=%d)=%.2e\n",
           N,o_t,
           N, exp($delta(T=>0)),
           N, $psi(T=>0),
           N,o_t, exp($b(M=>o_t)));
#endif
 %}

#ifdef DEBUG_VITERBI
   printf("\n");
#endif

 /*-- Main: t>0: Loop: time==t --*/
 for (t=1; t<$SIZE(T); t++) {
   o_t = $o(T=>t);

   /*-- Main: t>0: Loop: state_t==j --*/
   for (j=0; j<$SIZE(N); j++) {
     psi_jt   = 0;
     delta_jt = $delta(N=>0,T=>t-1) + $a(N0=>0,N1=>j);

     /*-- Main: t>0: Loop: state_(t-1)==i --*/
     for (i=1; i<$SIZE(N); i++) {
       delta_tmp = $delta(N=>i,T=>t-1) + $a(N0=>i,N1=>j);

       if (delta_tmp > delta_jt) {
	 delta_jt = delta_tmp;
	 psi_jt   = i;
#ifdef DEBUG_VITERBI
         printf("+");
#endif
       }

#ifdef DEBUG_VITERBI
       printf("t=%d,i=%d,j=%d,o_t=%d:  deltaX(i=%d,t=%d)=%.2e  psi(j=%d,t=%d)=%.0g  delta(j=%d,t=%d)=%.2e\n",
              t,i,j,o_t,
              i,t,  exp(delta_tmp),
              j,t,  psi_jt,
              j,t,  exp(delta_jt));
#endif
     }

     /*-- Main: t>0: Store data for state,time=(j,t) --*/
     $delta(N=>j,T=>t) = delta_jt + $b(N=>j,M=>o_t);
     $psi  (N=>j,T=>t) = psi_jt;

#ifdef DEBUG_VITERBI
     printf("\n---> t=%d: b(j=%d,o=%d)=%.2e  delta(j=%d,t=%d)=%.2e  psi(j=%d,t=%d)=%.0g\n\n",
            t, j,o_t, exp($b(N=>j,M=>o_t)),
            j,t,      exp($delta(N=>j,T=>t)),
            j,t,      $psi(N=>j,T=>t));
#endif
   }
#ifdef DEBUG_VITERBI
   printf("\n");
#endif
 }
'),
  Doc =>
  ('Computes Viterbi algorithm trellises $delta() and $psi() for the
observation sequence $o() given the model parameters @theta = ($a,$b,$pi,$omega).

Outputs:

Probability matrix $delta(): log probability of best path to state $j at time $t:

 $delta(j,t) = max_{q(0:t)} log P( $o(0:t), q(0:t-1), $q(t)==j | @theta )

Path backtrace matrix $psi(): best predecessor for state $j at time $t:

 $psi(j,t) = arg_{q(t-1)} max_{q(0:t)} P( $o(0:t), q(0:t-1), $q(t)==j | @theta )

Note that if you are using termination probabilities $omega(),
then in order to find the most likely final state, you need to
compute the contribution of $omega() yourself, which is easy
to do:

 $best_final_q = maximum_ind($delta->slice(",-1") + $omega);

'),
);


##--------------------------------------------------------------
## Sequence Analysis: Viterbi (constrained)
pp_def
('hmmviterbiq',
 GenericTypes => $FLOAT_TYPES,
 Pars => join(" ",
              qw(a(N,N);
                 b(N,M);
                 pi(N);),
              qw(o(T);
                 oq(Q,T);),
              '[o]deltaq(Q,T);',
              'indx [o]psiq(Q,T)',
             ),
 Code =>
('
 int qi,qj, i,j, t, o_t;
 double deltaq_jt, deltaq_tmp;
 int psiq_jt;

 /*-- Initialize: t=0: Loop: q_(0)=qi: state_(0)=oq(qi,0)=i --*/
 o_t = $o(T=>0);
 for (qi=0; qi<$SIZE(Q); qi++) {
   i = $oq(Q=>qi,T=>0);
   $psiq(Q=>qi,T=>0) = 0;
   $deltaq(Q=>qi,T=>0) = (i>=0 ? ($pi(N=>i)+$b(N=>i,M=>o_t)) : ($GENERIC(deltaq))LOG_ZERO);
 }

 /*-- Loop: t>0: Loop: time==t --*/
 for (t=1; t<$SIZE(T); t++) {
   o_t = $o(T=>t);

   /*-- Loop: t>0: q_(t)=qj : state_(t)=oq(qj,t)=j --*/
   for (qj=0; qj<$SIZE(Q); qj++) {
     j         = $oq(Q=>qj,T=>t);
     i         = $oq(Q=>0, T=>t-1);
     psiq_jt   = 0;

     if (j >= 0 && i >=0) {
       deltaq_jt = $deltaq(Q=>0,T=>t-1) + $a(N0=>i,N1=>j);
     } else {
       deltaq_jt = $deltaq(Q=>0,T=>t-1) + LOG_ZERO;
     }

     /*-- Loop: t>0: q_(t-1)=qi : state_(t-1)=oq(qi,t)=i --*/
     for (qi=1; qi<$SIZE(Q); qi++) {
       i         = $oq(Q=>qi,T=>t-1);
       if (j < 0 || i < 0) break;

       deltaq_tmp = $deltaq(Q=>qi,T=>t-1) + $a(N0=>i,N1=>j);

       if (deltaq_tmp > deltaq_jt) {
	 deltaq_jt = deltaq_tmp;
	 psiq_jt   = qi;
       }

     }
     /*-- End Loop: t>0: q_(t-1)=qi : state_(t-1)=oq(qi,t)=i --*/

     /*-- Main: t>0: Store data for stateIndex,time=(qj,t) --*/
     $deltaq(Q=>qj,T=>t) = deltaq_jt + (j>=0 ? $b(N=>j,M=>o_t) : LOG_ZERO);
     $psiq  (Q=>qj,T=>t) = psiq_jt;

   }
   /*-- Loop: t>0: q_(t)=qj : state_(t)=oq(qj,t)=j --*/

 }
 /*-- End Loop: t>0: Loop: time==t --*/
'
),
  Doc =>
  ('Computes constrained Viterbi algorithm trellises $deltaq() and $psiq() for the
observation sequence $o() given the model parameters @theta = ($a,$b,$pi,$omega),
considering only the initial non-negative state indices $oq(:,t) for each
observarion $o(t).

Outputs:

Constrained probability matrix $deltaq(): log probability of best path to state $oq(j,t) at time $t:

 $deltaq(j,t) = max_{j(0:t)} log P( $o(0:t), q(0:t-1)==$oq(:,j(0:t-1)), q(t)==$oq(j,t) | @theta )

Constrained path backtrace matrix $psiq(): best predecessor index for state $oq(j,t) at time $t:

 $psiq(j,t) = arg_{j(t-1)} max_{j(0:t)} P( $o(0:t), q(0:t-1)=$oq(:,j(0:t-1)), q(t)==$oq(j,t) | @theta )

Note that if you are using termination probabilities $omega(),
then in order to find the most likely final state, you need to
compute the contribution of $omega() yourself, which is quite easy
to do:

 $best_final_j = maximum_ind($deltaq->slice(",-1") + $omega->index($oq->slice(",(-1)")))

'),
);


##--------------------------------------------------------------
## Sequence Analysis: Backtrace
pp_def
('hmmpath',
 Pars => q(psi(N,T); indx qfinal(); indx [o]path(T)),
 GenericTypes => $FLOAT_TYPES,
 Code =>
('
 /*-- Initialize: t==T-1: state_(t)==final() --*/
 int t = $SIZE(T)-1;
 $path(T=>t) = $qfinal();

 /*-- Main: T-1 > t >= 0: Loop: time==t --*/
 for (t--; t>=0; t--) {
   int q_tp1   = $path(T=>t+1);
   $path(T=>t) = $psi (T=>t+1,N=>q_tp1);
 }
'),
  Doc =>
  ('Computes best-path backtrace $path() for the final state $qfinal()
from completed Viterbi trellis $psi().

Outputs:

Path backtrace $path(): state (in best sequence) at time $t:

 $path(t) = arg_{q(t)} max_{q(0:T-1)} log P( $o(), q(0:T-2), $q(T-1)==$qfinal() | @theta )

This even threads over multiple final states, if specified,
so you can align paths to their final states just by calling:

 $bestpaths = hmmpath($psi, sequence($N));

Note that $path(T-1) == $qfinal(): yes, this is redundant,
but also tends to be quite convenient.

'),
);

##--------------------------------------------------------------
## Sequence Analysis: Backtrace (constrained)
pp_def
('hmmpathq',
 Pars => q(indx oq(Q,T); psiq(Q,T); indx qfinalq(); indx [o]path(T)),
 GenericTypes => $FLOAT_TYPES,
 Code =>
('
 /*-- Initialize: t==T-1: state_(t)==final() --*/
 int t = $SIZE(T)-1;
 $path(T=>t) = $qfinalq();

 /*-- Get index backtrace --*/
 for (t--; t>=0; t--) {
   int qi_tp1  = $path(T=>t+1);
   $path(T=>t) = $psiq(T=>t+1,Q=>qi_tp1);
 }

 /*-- Convert indices to state ids --*/
 loop (T) %{
   int qi  = $path();
   $path() = $oq(Q=>qi);
 %}
'),
  Doc =>
  ('Computes constrained best-path backtrace $path() for the final state index $qfinalq()
from completed constrained Viterbi trellis $psiq().

Outputs:

Path backtrace $path(): state (in best sequence) at time $t:

 $path(t) = arg_{q(t)} max_{q(0:T-1)} log P( $o(), q(0:T-2), $q(T-1)==$oq($qfinalq(),T-1) | @theta )

This is really just a convenience method for dealing with constrained
lookup -- the same thing can be accomplished using hmmpath() and
some PDL index magic.

'),
);



##======================================================================
## Footer Administrivia
##======================================================================

##------------------------------------------------------
## pm additions
pp_addpm(<<'EOPM');


##---------------------------------------------------------------------
=pod

=head1 COMMON PARAMETERS

HMMs are specified by parameters $a(N,N), $b(N,M), $pi(N), and $omega(N);
input sequences are represented as vectors $o(T) of integer values in the range [0..M-1],
where the following notational conventions are used:

=over 4


=item States:

The model has $N states, denoted $q,
0 <= $q < $N.


=item Alphabet:

The input- (aka "observation-") alphabet of the model has $M elements,
denoted $o(t), 0 <= $o(t) < $M.


=item Time indices:

Time indices are denoted $t,
1 <= $t < $T.


=item Input Sequences:

Input- (aka "observation-") sequences are represented as vectors of
of length $T whose component values are in the range [0..M-1],
i.e. alphabet indices.



=item Initial Probabilities:

The vector $pi(N) gives the (log) initial state probability distribution:

 $pi(i) = log P( $q(0)==i )



=item Final Probabilities:

The vector $omega(N) gives the (log) final state probability distribution:

 $omega(i) = log P( $q($T)==i )

This parameter is a nonstandard extension.
You can simulate the behavior of more traditional definitions
(such as that given in Rabiner (1989)) by setting:

 $omega = zeroes($N);

wherever it is required.



=item Arc Probabilities:

The matrix $a(N,N) gives the (log) conditional state-transition probability distribution:

 $a(i,j) = log P( $q(t+1)==j | $q(t)==i )



=item Emission Probabilities:

The matrix $b(N,M) gives the (log) conditional symbol emission probability:

 $b(j,o) = log P( $o(t)==o | $q(t)==j )



=back

=cut

##---------------------------------------------------------------------
=pod

=head1 ACKNOWLEDGEMENTS

Perl by Larry Wall.

PDL by Karl Glazebrook, Tuomas J. Lukka, Christian Soeller, and others.

Implementation based largely on the formulae in:
L. E. Rabiner, "A tutorial on Hidden Markov Models and selected
applications in speech recognition," Proceedings of the IEEE 77:2,
Februrary, 1989, pages 257--286.

=cut

##----------------------------------------------------------------------
=pod

=head1 KNOWN BUGS

Probably many.

=cut


##---------------------------------------------------------------------
=pod

=head1 AUTHOR

Bryan Jurish [e-mail]

=head2 Copyright Policy

Copyright (C) 2005-2023 by Bryan Jurish. All rights reserved.

This package is free software, and entirely without warranty.
You may redistribute it and/or modify it under the same terms
as Perl itself.

=head1 SEE ALSO

perl(1), PDL(3perl).

=cut

EOPM


# Always make sure that you finish your PP declarations with
# pp_done
pp_done();
##----------------------------------------------------------------------