File Coverage

HMM.xs

Criterion	Covered	Total	%
statement	35	53	66.0
branch	114	456	25.0
condition			n/a
subroutine			n/a
pod			n/a
total	149	509	29.2

#line 453 "lib/PDL/PP.pm"
/*
 * THIS FILE WAS GENERATED BY PDL::PP from HMM.pd! Do not modify!
 */

#define PDL_FREE_CODE(trans, destroy, comp_free_code, ntpriv_free_code) \
    if (destroy) { \
	comp_free_code \
    } \
    if ((trans)->dims_redone) { \
	ntpriv_free_code \
    }

#include "EXTERN.h"
#include "perl.h"
#include "XSUB.h"
#include "pdl.h"
#include "pdlcore.h"
#define PDL PDL_HMM
extern Core* PDL; /* Structure hold core C functions */
#line 23 "HMM.xs"

#include

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


/* 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)

pdl_error pdl_run_logzero(pdl  *a);
pdl_error pdl_run_logadd(pdl  *a,pdl  *b,pdl  *c);
pdl_error pdl_run_logdiff(pdl  *a,pdl  *b,pdl  *c);
pdl_error pdl_run_logsumover(pdl  *a,pdl  *b);
pdl_error pdl_run_hmmfw(pdl  *a,pdl  *b,pdl  *pi,pdl  *o,pdl  *alpha);
pdl_error pdl_run_hmmfwq(pdl  *a,pdl  *b,pdl  *pi,pdl  *o,pdl  *oq,pdl  *alphaq);
pdl_error pdl_run_hmmbw(pdl  *a,pdl  *b,pdl  *omega,pdl  *o,pdl  *beta);
pdl_error pdl_run_hmmbwq(pdl  *a,pdl  *b,pdl  *omega,pdl  *o,pdl  *oq,pdl  *betaq);
pdl_error pdl_run_hmmexpect(pdl  *a,pdl  *b,pdl  *pi,pdl  *omega,pdl  *o,pdl  *alpha,pdl  *beta,pdl  *ea,pdl  *eb,pdl  *epi,pdl  *eomega);
pdl_error pdl_run_hmmexpectq(pdl  *a,pdl  *b,pdl  *pi,pdl  *omega,pdl  *o,pdl  *oq,pdl  *alphaq,pdl  *betaq,pdl  *ea,pdl  *eb,pdl  *epi,pdl  *eomega);
pdl_error pdl_run_hmmviterbi(pdl  *a,pdl  *b,pdl  *pi,pdl  *o,pdl  *delta,pdl  *psi);
pdl_error pdl_run_hmmviterbiq(pdl  *a,pdl  *b,pdl  *pi,pdl  *o,pdl  *oq,pdl  *deltaq,pdl  *psiq);
pdl_error pdl_run_hmmpath(pdl  *psi,pdl  *qfinal,pdl  *path);
pdl_error pdl_run_hmmpathq(pdl  *oq,pdl  *psiq,pdl  *qfinalq,pdl  *path);

#line 449 "lib/PDL/PP.pm"
#include "pdlperl.h"
#line 93 "HMM.xs"

Core* PDL = NULL; /* Structure hold core C functions */

MODULE = PDL::HMM PACKAGE = PDL PREFIX=pdl_run_

PROTOTYPES: DISABLE


void
pdl_run_logzero(a=PDL_XS_PERLINIT_initsv(a_SV))
  SV *a_SV = items > 0 ? ST(0) : NULL;
 PREINIT:
  PDL_XS_PREAMBLE((items == 1) ? 0 : 1);
 INPUT:
  pdl *a
 PPCODE:
  PDL->barf_if_error(pdl_run_logzero(a));
  PDL_XS_RETURN(ST(0) = a_SV)

void
pdl_run_logadd(a, b, c=c)
  SV *c_SV = items > 2 ? ST(2) : NULL;
 PREINIT:
  PDL_XS_PREAMBLE((items == 3) ? 0 : 1);
 INPUT:
  pdl *a
  pdl *b
  pdl *c=NO_INIT
 PPCODE:
  PDL_XS_INPLACE(a, c)
  PDL->barf_if_error(pdl_run_logadd(a,b,c));
  PDL_XS_RETURN(ST(0) = c_SV)

void
pdl_run_logdiff(a, b, c=c)
  SV *c_SV = items > 2 ? ST(2) : NULL;
 PREINIT:
  PDL_XS_PREAMBLE((items == 3) ? 0 : 1);
 INPUT:
  pdl *a
  pdl *b
  pdl *c=NO_INIT
 PPCODE:
  PDL_XS_INPLACE(a, c)
  PDL->barf_if_error(pdl_run_logdiff(a,b,c));
  PDL_XS_RETURN(ST(0) = c_SV)

void
pdl_run_logsumover(a, b=PDL_XS_PERLINIT_initsv(b_SV))
  SV *b_SV = items > 1 ? ST(1) : NULL;
 PREINIT:
  PDL_XS_PREAMBLE((items == 2) ? 0 : 1);
 INPUT:
  pdl *a
  pdl *b
 PPCODE:
  PDL->barf_if_error(pdl_run_logsumover(a,b));
  PDL_XS_RETURN(ST(0) = b_SV)

void
pdl_run_hmmfw(a, b, pi, o, alpha=PDL_XS_PERLINIT_initsv(alpha_SV))
  SV *alpha_SV = items > 4 ? ST(4) : NULL;
 PREINIT:
  PDL_XS_PREAMBLE((items == 5) ? 0 : 1);
 INPUT:
  pdl *a
  pdl *b
  pdl *pi
  pdl *o
  pdl *alpha
 PPCODE:
  PDL->barf_if_error(pdl_run_hmmfw(a,b,pi,o,alpha));
  PDL_XS_RETURN(ST(0) = alpha_SV)

void
pdl_run_hmmfwq(a, b, pi, o, oq, alphaq=PDL_XS_PERLINIT_initsv(alphaq_SV))
  SV *alphaq_SV = items > 5 ? ST(5) : NULL;
 PREINIT:
  PDL_XS_PREAMBLE((items == 6) ? 0 : 1);
 INPUT:
  pdl *a
  pdl *b
  pdl *pi
  pdl *o
  pdl *oq
  pdl *alphaq
 PPCODE:
  PDL->barf_if_error(pdl_run_hmmfwq(a,b,pi,o,oq,alphaq));
  PDL_XS_RETURN(ST(0) = alphaq_SV)

void
pdl_run_hmmbw(a, b, omega, o, beta=PDL_XS_PERLINIT_initsv(beta_SV))
  SV *beta_SV = items > 4 ? ST(4) : NULL;
 PREINIT:
  PDL_XS_PREAMBLE((items == 5) ? 0 : 1);
 INPUT:
  pdl *a
  pdl *b
  pdl *omega
  pdl *o
  pdl *beta
 PPCODE:
  PDL->barf_if_error(pdl_run_hmmbw(a,b,omega,o,beta));
  PDL_XS_RETURN(ST(0) = beta_SV)

void
pdl_run_hmmbwq(a, b, omega, o, oq, betaq=PDL_XS_PERLINIT_initsv(betaq_SV))
  SV *betaq_SV = items > 5 ? ST(5) : NULL;
 PREINIT:
  PDL_XS_PREAMBLE((items == 6) ? 0 : 1);
 INPUT:
  pdl *a
  pdl *b
  pdl *omega
  pdl *o
  pdl *oq
  pdl *betaq
 PPCODE:
  PDL->barf_if_error(pdl_run_hmmbwq(a,b,omega,o,oq,betaq));
  PDL_XS_RETURN(ST(0) = betaq_SV)

void
pdl_run_hmmexpect(a, b, pi, omega, o, alpha, beta, ea=PDL_XS_PERLINIT_initsv(ea_SV), eb=PDL_XS_PERLINIT_initsv(eb_SV), epi=PDL_XS_PERLINIT_initsv(epi_SV), eomega=PDL_XS_PERLINIT_initsv(eomega_SV))
  SV *ea_SV = items > 7 ? ST(7) : NULL;
  SV *eb_SV = items > 8 ? ST(8) : NULL;
  SV *epi_SV = items > 9 ? ST(9) : NULL;
  SV *eomega_SV = items > 10 ? ST(10) : NULL;
 PREINIT:
  PDL_XS_PREAMBLE((items == 11) ? 0 : 4);
 INPUT:
  pdl *a
  pdl *b
  pdl *pi
  pdl *omega
  pdl *o
  pdl *alpha
  pdl *beta
  pdl *ea
  pdl *eb
  pdl *epi
  pdl *eomega
 PPCODE:
  PDL->barf_if_error(pdl_run_hmmexpect(a,b,pi,omega,o,alpha,beta,ea,eb,epi,eomega));
  PDL_XS_RETURN(ST(0) = ea_SV;ST(1) = eb_SV;ST(2) = epi_SV;ST(3) = eomega_SV)

void
pdl_run_hmmexpectq(a, b, pi, omega, o, oq, alphaq, betaq, ea=PDL_XS_PERLINIT_initsv(ea_SV), eb=PDL_XS_PERLINIT_initsv(eb_SV), epi=PDL_XS_PERLINIT_initsv(epi_SV), eomega=PDL_XS_PERLINIT_initsv(eomega_SV))
  SV *ea_SV = items > 8 ? ST(8) : NULL;
  SV *eb_SV = items > 9 ? ST(9) : NULL;
  SV *epi_SV = items > 10 ? ST(10) : NULL;
  SV *eomega_SV = items > 11 ? ST(11) : NULL;
 PREINIT:
  PDL_XS_PREAMBLE((items == 12) ? 0 : 4);
 INPUT:
  pdl *a
  pdl *b
  pdl *pi
  pdl *omega
  pdl *o
  pdl *oq
  pdl *alphaq
  pdl *betaq
  pdl *ea
  pdl *eb
  pdl *epi
  pdl *eomega
 PPCODE:
  PDL->barf_if_error(pdl_run_hmmexpectq(a,b,pi,omega,o,oq,alphaq,betaq,ea,eb,epi,eomega));
  PDL_XS_RETURN(ST(0) = ea_SV;ST(1) = eb_SV;ST(2) = epi_SV;ST(3) = eomega_SV)

void
pdl_run_hmmviterbi(a, b, pi, o, delta=PDL_XS_PERLINIT_initsv(delta_SV), psi=PDL_XS_PERLINIT_initsv(psi_SV))
  SV *delta_SV = items > 4 ? ST(4) : NULL;
  SV *psi_SV = items > 5 ? ST(5) : NULL;
 PREINIT:
  PDL_XS_PREAMBLE((items == 6) ? 0 : 2);
 INPUT:
  pdl *a
  pdl *b
  pdl *pi
  pdl *o
  pdl *delta
  pdl *psi
 PPCODE:
  PDL->barf_if_error(pdl_run_hmmviterbi(a,b,pi,o,delta,psi));
  PDL_XS_RETURN(ST(0) = delta_SV;ST(1) = psi_SV)

void
pdl_run_hmmviterbiq(a, b, pi, o, oq, deltaq=PDL_XS_PERLINIT_initsv(deltaq_SV), psiq=PDL_XS_PERLINIT_initsv(psiq_SV))
  SV *deltaq_SV = items > 5 ? ST(5) : NULL;
  SV *psiq_SV = items > 6 ? ST(6) : NULL;
 PREINIT:
  PDL_XS_PREAMBLE((items == 7) ? 0 : 2);
 INPUT:
  pdl *a
  pdl *b
  pdl *pi
  pdl *o
  pdl *oq
  pdl *deltaq
  pdl *psiq
 PPCODE:
  PDL->barf_if_error(pdl_run_hmmviterbiq(a,b,pi,o,oq,deltaq,psiq));
  PDL_XS_RETURN(ST(0) = deltaq_SV;ST(1) = psiq_SV)

void
pdl_run_hmmpath(psi, qfinal, path=PDL_XS_PERLINIT_initsv(path_SV))
  SV *path_SV = items > 2 ? ST(2) : NULL;
 PREINIT:
  PDL_XS_PREAMBLE((items == 3) ? 0 : 1);
 INPUT:
  pdl *psi
  pdl *qfinal
  pdl *path
 PPCODE:
  PDL->barf_if_error(pdl_run_hmmpath(psi,qfinal,path));
  PDL_XS_RETURN(ST(0) = path_SV)

void
pdl_run_hmmpathq(oq, psiq, qfinalq, path=PDL_XS_PERLINIT_initsv(path_SV))
  SV *path_SV = items > 3 ? ST(3) : NULL;
 PREINIT:
  PDL_XS_PREAMBLE((items == 4) ? 0 : 1);
 INPUT:
  pdl *oq
  pdl *psiq
  pdl *qfinalq
  pdl *path
 PPCODE:
  PDL->barf_if_error(pdl_run_hmmpathq(oq,psiq,qfinalq,path));
  PDL_XS_RETURN(ST(0) = path_SV)


#line 483 "lib/PDL/PP.pm"
BOOT:
   /* Get pointer to structure of core shared C routines */
   /* make sure PDL::Core is loaded */
#line 331 "HMM.xs"
   perl_require_pv ("PDL/Core.pm"); /* make sure PDL::Core is loaded */
#ifndef aTHX_
#define aTHX_
#endif
   if (SvTRUE (ERRSV)) Perl_croak(aTHX_ "%s",SvPV_nolen (ERRSV));
   SV* CoreSV = perl_get_sv("PDL::SHARE",FALSE); /* var with core structure */
   if (!CoreSV)
     Perl_croak(aTHX_ "We require the PDL::Core module, which was not found");
   if (!(PDL = INT2PTR(Core*,SvIV( CoreSV )))) /* Core* value */
     Perl_croak(aTHX_ "Got NULL pointer for PDL");
   if (PDL->Version != PDL_CORE_VERSION)
     Perl_croak(aTHX_ "[PDL->Version: %ld PDL_CORE_VERSION: %ld XS_VERSION: %s] PDL::HMM needs to be recompiled against the newly installed PDL", (long int)PDL->Version, (long int)PDL_CORE_VERSION, XS_VERSION);