typedef struct {
  char *ALF;
  int nullchar;
  int *alfinds;
  int junkl;
  int junka;
  int nomorph;
  int nograd;
  int nalf;
  int Nwords;
  int infrloop;
  double infrtol;
  double gradstep;
  int learnloop;
  int Nlocal;
  int Ntopic;
  /*int favins;*/
  int maxsenlen;
  int maxstrlen;
} auxparamstype;


#define NAUX 30


#define Ut(t) (((int)(_(U,0,t)))-1)
#define ui(t) (u[t])
#define vj(t) (v[t])
#define UVARR(uv) mxGetCell(RAWVOC,uv)
#define UVLEN(uv) mxGetN(UVARR(uv))
#define E10 (e1==NULLCHAR)
#define E20 (e2==NULLCHAR)
#define E1u (e1==ui(i-1))
#define E2v (e2==vj(j-1))
#define PMPuv PMP(ui(i-1),vj(j-1))
#define PMPu0 PMP(ui(i-1),NULLCHAR)
#define PMP0v PMP(NULLCHAR,vj(j-1))
#define PMP00 PMP(NULLCHAR,NULLCHAR)
#define RHO0 RHO(NULLCHAR)
#define RHOu RHO(ui(i-1))
#define DRHO(x,y,z) DRhoDnu(x,y,z,AUXPARAMS)

#define WEQ(j,t) ( ABS(_(Uniq,0,j)-_(U,0,t)) < EPS)
/* P(rimitive) M(orphological) P(robability) */
#define PMP(a,b) _3(pmp,AUXPARAMS.alfinds[a-1],AUXPARAMS.alfinds[b-1],ell)
#define RHO(a) _(rho,AUXPARAMS.alfinds[a-1],ell)
/*#define JUNKELL 0*/
#define JUNKELL (AUXPARAMS.junkl)
#define exists_junkell (JUNKELL>=0)
/* #define NULLCHAR ('@')*/
#define NULLCHAR (AUXPARAMS.nullchar)
#define JUNKTOP (AUXPARAMS.junka)
#define NOMORPH (AUXPARAMS.nomorph)
#define INFRLOOP (AUXPARAMS.infrloop)
#define ITOL (AUXPARAMS.infrtol)
#define MAXSTRLEN (AUXPARAMS.maxstrlen)
#define MAXSENLEN (AUXPARAMS.maxsenlen)
#define NTOPIC (AUXPARAMS.Ntopic)
#define NLOCAL (AUXPARAMS.Nlocal)
#define NVOCAB (AUXPARAMS.Nwords)
#define NALF (AUXPARAMS.nalf)
#define GETFLDR(X,s) (*mxGetPr(mxGetField(X,0,s)))


void setauxparams(auxparamstype *AP,const mxArray *parampr)
{
  const mxArray *ALFAr,*ALFindAr;
  int i,nalf,*alfinds,NCI;
  char *ALF;
      ALFAr    =       mxGetField(parampr,0,"ALF0");
      ALFindAr =       mxGetField(parampr,0,"alfinds");
      NCI      =(int)GETFLDR(parampr,"NCI")-1;     /*index of NULLCHAR in ALFINDS */
  AP->junkl    =(int)GETFLDR(parampr,"JUNKELL")-1; /*junk state index */
  AP->junka    =(int)GETFLDR(parampr,"JUNKTOP")-1; /*junk topic index */
  AP->nomorph  =(int)GETFLDR(parampr,"nomorph");   /*don't do any morphology */
  AP->nograd   =(int)GETFLDR(parampr,"nograd");   /*don't do the morph. gradient */
  AP->infrloop =(int)GETFLDR(parampr,"infrloop");  /*how long to iterate for qa and rv? */
  AP->learnloop=(int)GETFLDR(parampr,"learnloop"); /*how many times to iterate updates */
  AP->Nlocal   =(int)GETFLDR(parampr,"Nlocal");    /*number of local states */
  AP->Ntopic   =(int)GETFLDR(parampr,"Ntopic");    /*number of topics */
  AP->Nwords   =(int)GETFLDR(parampr,"Nwords");    /*number of raw vocab words */
  /*AP->favins   =(int)GETFLDR(parampr,"favins");    /*favor insertions over deletions */
  AP->infrtol  =     GETFLDR(parampr,"infrtol");   /*precision for inference */
  AP->gradstep =     GETFLDR(parampr,"gradstep");  /*gradient step size */
  AP->maxsenlen=(int)GETFLDR(parampr,"maxsenlen"); /*maximum doc length */
  AP->maxstrlen=(int)GETFLDR(parampr,"maxstrlen"); /*maximum string length */

  /* alphabet-handling bureaucracy */
  nalf = mxGetN(ALFAr);
  ALF  = mxCalloc(nalf+1,sizeof(mxChar));
  i = mxGetString(ALFAr,ALF,nalf*sizeof(mxChar)+1);
  alfinds = mxCalloc(mxGetM(ALFindAr),sizeof(int));
  for(i=0;i<mxGetM(ALFindAr);i++)
    alfinds[i] = (int)(mxGetPr(ALFindAr)[i])-1;

  AP->ALF = ALF;
  AP->alfinds = alfinds;
  AP->nullchar = ALF[NCI];
  AP->nalf = nalf;
}

void freeauxparams(auxparamstype AUXPARAMS)
{
  mxFree(AUXPARAMS.ALF);
  mxFree(AUXPARAMS.alfinds);
}



void normpmp(TENSOR3 pmp,auxparamstype AUXPARAMS)
{
  int ell,i,j;
  double s;
  for(ell=0;ell<NLOCAL;ell++) {
    s = 0.0;
    for(i=0;i<NALF;i++)
      for(j=0;j<NALF;j++)
	s += _3(pmp,i,j,ell);
    for(i=0;i<NALF;i++)
      for(j=0;j<NALF;j++)
	_3(pmp,i,j,ell) /= s;
  }
}


void init_rho(rho,pmp,AUXPARAMS)
     TENSOR3 pmp;
     TENSOR2 rho;
     auxparamstype AUXPARAMS;
{
  int ei,ej,ell;
  char a,b;
  for(ei=0;ei<NALF;ei++) {
    a = AUXPARAMS.ALF[ei];
    for(ell=0;ell<NLOCAL;ell++) {
      RHO(a) = 0.0;
      for(ej=0;ej<NALF;ej++) {
	b = AUXPARAMS.ALF[ej];
	RHO(a) += PMP(a,b) * (b!=NULLCHAR);
      }
    }
  }
}


void jointstrprob(psi,u,v,ell,pmp,AUXPARAMS)
     TENSOR2 psi;
     char *u,*v;
     int ell;
     TENSOR3 pmp;
     auxparamstype AUXPARAMS;
{
  int i,j;
  _(psi,0,0) = PMP00;
  for(i=0;i<d1(psi);i++)
    for(j=0;j<d2(psi);j++) {
      if(i||j) _(psi,i,j) = 0.0;
      if( i  ) _(psi,i,j) += PMPu0 * _(psi,i-1  ,j);
      if(   j) _(psi,i,j) += PMP0v * _(psi,i  ,j-1);
      if(i&&j) _(psi,i,j) += PMPuv * _(psi,i-1,j-1);
    }
}


void margstrprob(psm,psi1,A,B,u,ell,pmp,rho,AUXPARAMS,I_A)
     /* B = inv(I-A) */
     TENSOR2 psm,psi1,A,B,rho;
     char *u;
     int ell;
     TENSOR3 pmp;
     auxparamstype AUXPARAMS;
     TENSOR2 I_A; /* I-A */
{
  int i,j;
  /* compute psm */
  _(psm,0,0) = PMP00;
  for(i=1;i<d1(psm);i++)
    _(psm,i,0) = PMPu0 * _(psm,i-1,0);
  /* compute A */
  for(i=0;i<d1(A);i++) {
    _(A,i,i) = RHO0;
    if(i) _(A,i,i-1) = RHOu + RHO0 * PMPu0;
    for(j=0;j<i-1;j++)
      _(A,i,j) = _(A,i-1,j) * PMPu0;
    /* zero out rest of A */
    for(j=i+1;j<d2(A);j++)
      _(A,i,j) = 0.0;
  }
  /* compute (I-A) */
  for(i=0;i<d1(A);i++)
    for(j=0;j<d1(A);j++)
      _(I_A,i,j) = 1.0*(i==j) - _(A,i,j);
  /* B = inv(I_A); */
  inv_ld(I_A,B);
  /* psi1 = B*psm */
  multTENSOR2(B,psm,psi1);
}



double DRhoDnu(t,e1,e2,AUXPARAMS)
     char t,e1,e2;
     auxparamstype AUXPARAMS;
{
  return(1.0*(t==e1)*(e2!=NULLCHAR));
}

void DpsmDnu(dpsm,u,e1,e2,psm,ell,pmp,AUXPARAMS)
     TENSOR2 dpsm,psm;
     char *u,e1,e2;
     int ell;
     TENSOR3 pmp;
     auxparamstype AUXPARAMS;
{
  int i;
  _(dpsm,0,0) = 1.0 * E10 * E20;
  for(i=1;i<d1(dpsm);i++)
    _(dpsm,i,0) = _(psm,i-1,0)*E1u*E20 + _(dpsm,i-1,0) * PMPu0;
}


void DpsiDnu(dpsi,psi,u,v,e1,e2,ell,pmp,AUXPARAMS)
     TENSOR2 dpsi,psi;
     char *u,*v,e1,e2;
     int ell;
     TENSOR3 pmp;
     auxparamstype AUXPARAMS;
{
  int i,j;
  _(dpsi,0,0) = 1.0 * E10 * E20;
  for(i=0;i<d1(dpsi);i++)
    for(j=0;j<d2(dpsi);j++) {
      if(i||j) _(dpsi,i,j) = 0.0;
      if( i  ) _(dpsi,i,j) += _(psi,i-1,j  ) * E1u * E20 +  _(dpsi,i-1,j  ) * PMPu0;
      if(   j) _(dpsi,i,j) += _(psi,i  ,j-1) * E10 * E2v +  _(dpsi,i  ,j-1) * PMP0v;
      if(i&&j) _(dpsi,i,j) += _(psi,i-1,j-1) * E1u * E2v +  _(dpsi,i-1,j-1) * PMPuv;
    }
}


void DABDnu(dA,dB,u,e1,e2,A,B,ell,pmp,rho,AUXPARAMS)
     TENSOR2 A,B,dA,dB,rho;
     char *u,e1,e2;
     int ell;
     TENSOR3 pmp;
     auxparamstype AUXPARAMS;
{
  int i,j,k,l;
  /* compute dA */
  for(i=0;i<d1(dA);i++) {
    _(dA,i,i) = DRHO(NULLCHAR,e1,e2);
    if(i) _(dA,i,i-1) = DRHO(ui(i-1),e1,e2) + DRHO(NULLCHAR,e1,e2) * PMPu0 + RHO0 * E1u * E20;
    for(j=0;j<i-1;j++)
      _(dA,i,j) = _(dA,i-1,j) * PMPu0 + _(A,i-1,j) * E1u * E20;
    /* zero out rest of dA */
    for(j=i+1;j<d2(dA);j++)
      _(dA,i,j) = 0.0;
  }
  /* compute dB = B * dA * B */
  for(i=0;i<d1(dB);i++)
    for(j=0;j<d1(dB);j++) {
      _(dB,i,j) = 0.0;
      for(k=0;k<d1(dB);k++)
	for(l=0;l<d1(dB);l++)
	  _(dB,i,j) += _(B,i,k) * _(B,l,j) * _(dA,k,l);
    }
}



double dPidNu(ustr,vstr,ell,a,b,pmp,rho,AUXPARAMS,psm,psj,ps1,A,B,dpsm,dpsj,dA,dB)
     TENSOR2 rho;
     char a,b;
     int ell;
     TENSOR3 pmp;
     auxparamstype AUXPARAMS;
     TENSOR2 psm,psj,ps1,A,B,dpsm,dpsj,dA,dB;
     char *ustr,*vstr;
{
  int j;
  double pjo,pma,dpjo,dpma;
  pjo = _(psj,d1(psj)-1,d2(psj)-1); /* joint prob */
  pma = _(ps1,d1(ps1)-1,0); /* marg prob */
  DpsiDnu(dpsj,psj,vstr,ustr,a,b,ell,pmp,AUXPARAMS);
  dpjo = _(dpsj,d1(dpsj)-1,d2(dpsj)-1); /* diff joint prob */
  DpsmDnu(dpsm,vstr,a,b,psm,ell,pmp,AUXPARAMS);
  DABDnu(dA,dB,vstr,a,b,A,B,ell,pmp,rho,AUXPARAMS);
  dpma = 0.0; /* diff marg prob */
  for(j=0;j<d1(psm);j++)
    dpma += _(dB,d1(dB)-1,j)*_(psm,j,0) + _(B,d1(B)-1,j)*_(dpsm,j,0);
  return((pma*dpjo - pjo*dpma)/(REALMIN+pma)/(REALMIN+pma)); /* diff quotient rule */
}


void buildPWORD(PWORD,PCOND,mu,mu0,AUXPARAMS)
     /* PWORD is P(w|l,a) */
     TENSOR3 PWORD,PCOND;
     TENSOR2 mu,mu0;
     auxparamstype AUXPARAMS;
{
  int w,v,a,ell;
  for(w=0;w<NVOCAB;w++)
    for(ell=0;ell<NLOCAL;ell++)
      for(a=0;a<NTOPIC;a++) 
	if(ell==JUNKELL)
	  _3(PWORD,w,ell,a) = _(mu0,w,0);
	else {
	  _3(PWORD,w,ell,a) = 0.0;
	  for(v=0;v<NVOCAB;v++)
	    _3(PWORD,w,ell,a) += _(mu,v,a) * _3(PCOND,w,v,ell);
	}
}


void buildPCOND(PCOND,DPCOND,RAWVOC,pmp,rho,AUXPARAMS,ustr,vstr,psm,psj,ps1,A,B,dpsm,dpsj,dA,dB,I_A)
     TENSOR3 PCOND,pmp;
     TENSOR5 DPCOND;
     TENSOR2 rho;
     const mxArray *RAWVOC;
     auxparamstype AUXPARAMS;
     TENSOR2 psm,psj,ps1,A,B,dpsm,dpsj,dA,dB,I_A;
     char *ustr,*vstr;
{
  int u,v,ell,ei,ej;
  char a,b;
  double pma,pjo;

  for(ell=0;ell<NLOCAL;ell++)
    if((ell==JUNKELL)||NOMORPH)
      for(u=0;u<NVOCAB;u++)
	for(v=0;v<NVOCAB;v++)
	  _3(PCOND,u,v,ell) = EPS*(u!=v) + (1.0-EPS*(NVOCAB-1))*(u==v);
    else {
      for(v=0;v<NVOCAB;v++) {
	mxGetString(UVARR(v),vstr,UVLEN(v)*sizeof(mxChar)+1);
	setdim2(psm,UVLEN(v)+1,1);
	setdim2(ps1,UVLEN(v)+1,1);
	setdim2(A,UVLEN(v)+1,UVLEN(v)+1);
	setdim2(I_A,UVLEN(v)+1,UVLEN(v)+1);
	setdim2(B,UVLEN(v)+1,UVLEN(v)+1);
	setdim2(dpsm,UVLEN(v)+1,1);
	setdim2(dA,UVLEN(v)+1,UVLEN(v)+1);
	setdim2(dB,UVLEN(v)+1,UVLEN(v)+1);
	margstrprob(psm,ps1,A,B,vstr,ell,pmp,rho,AUXPARAMS,I_A);
	pma = _(ps1,d1(ps1)-1,0); /* marg prob */
	for(u=0;u<NVOCAB;u++) {
	  mxGetString(UVARR(u),ustr,UVLEN(u)*sizeof(mxChar)+1);
	  setdim2(psj,UVLEN(v)+1,UVLEN(u)+1);
	  setdim2(dpsj,UVLEN(v)+1,UVLEN(u)+1);
	  jointstrprob(psj,vstr,ustr,ell,pmp,AUXPARAMS);
	  pjo = _(psj,d1(psj)-1,d2(psj)-1); /* joint prob */
	  _3(PCOND,u,v,ell) = pjo / (pma+REALMIN);
	  if(!AUXPARAMS.nograd)
	    for(ei=0;ei<NALF;ei++)
	      for(ej=0;ej<NALF;ej++) {
		a = AUXPARAMS.ALF[ei];
		b = AUXPARAMS.ALF[ej];
		_5(DPCOND,u,v,ell,ei,ej) = dPidNu(ustr,vstr,ell,a,b,pmp,rho,AUXPARAMS,psm,psj,ps1,A,B,dpsm,dpsj,dA,dB);
	      }
	}
      }
    }
}


double val_se1(mxArray *A,int i,int j)
/* returns the value (im/ex-plicit) of A(i,j) */
/* if an entry does not exist, its value is LOGREALMIN */
{
  double *xpr;
  xpr = get_se(A,i,j);
  if(xpr==NULL)
    return(LOGREALMIN);   /* entry does not exist */
  else
    return(*xpr);
}

void markovBq(U,qa,PWORD,lB,AUXPARAMS)
     TENSOR2 U,qa;
     TENSOR3 PWORD;
     mxArray *lB;
     auxparamstype AUXPARAMS;
{
  double s;
  int Uj,j,a,ell;
  for(j=0;j<d2(U);j++) {
    Uj = Ut(j);
    for(ell=0;ell<NLOCAL;ell++) {
      s = 0.0;
      for(a=0;a<NTOPIC;a++)
	s += _(qa,a,j)*LOG(_3(PWORD,Uj,ell,a));
      set_se(lB,Uj,ell,s);
    }
  }
}


void markovBr(U,qa,rv,PCOND,mu,mu0,lB,AUXPARAMS)
     TENSOR2 U,qa,rv,mu,mu0;
     TENSOR3 PCOND;
     mxArray *lB;
     auxparamstype AUXPARAMS;
{
  double sa,sv;
  int Uj,j,a,ell,v;

  for(j=0;j<d2(U);j++) {
    Uj = Ut(j);
    for(ell=0;ell<NLOCAL;ell++) {
      sa = 0.0;
      for(a=0;a<NTOPIC;a++) {
	sv = 0.0;
	for(v=0;v<NVOCAB;v++)
	  sv += _(rv,v,j)*LOG((_(mu,v,a)*(ell!=JUNKELL)+_(mu0,Uj,0)*(ell==JUNKELL))*_3(PCOND,Uj,v,ell));
	sa += _(qa,a,j)*sv;
      }
      set_se(lB,Uj,ell,sa);
    }
  }
}




void compute_lksi(lksi,lalf,lbet,lA,lB,U,lfp)
     TENSOR2 U,lA,lalf,lbet;
     TENSOR3 lksi;
     double lfp;
     mxArray *lB;
{
  int i,j,t;
  FOR3(i,j,t,d1(lA),d1(lA),d2(U)-1)
    _3(lksi,j,i,t) = _(lalf,i,t)+_(lA,j,i)+val_se1(lB,Ut(t+1),j)+_(lbet,j,t+1)-lfp;
}








double zerofunc()
{
     return 0.0;
}

void update_musum(musum,mu0sum,U,Uniq,qa,rv,rl,AUXPARAMS)
     TENSOR2 musum,mu0sum,U,Uniq,qa,rv,rl;
     auxparamstype AUXPARAMS;
{
  int v,a,j,t,ell;

  for(t=0;t<d2(U);t++)
    for(j=0;j<d2(Uniq);j++)
      if(WEQ(j,t))
	for(v=0;v<NVOCAB;v++)
	  for(a=0;a<NTOPIC;a++)
	    for(ell=0;ell<NLOCAL;ell++)
	      if(ell!=JUNKELL)
		_(musum,v,a) += _(rl,ell,t) * _(rv,v,j) * _(qa,a,j);

  if(exists_junkell)
    for(t=0;t<d2(U);t++)
      _(mu0sum,Ut(t),0) += _(rl,JUNKELL,t);
}

void update_bcsum(bsum,csum,lksi,rl)
     TENSOR2 bsum,csum,rl;
     TENSOR3 lksi;
{
  int i,j,k;
  for(i=0;i<d1(lksi);i++)
    for(j=0;j<d2(lksi);j++)
      for(k=0;k<d3(lksi);k++)
	_(bsum,i,j) += EXP(_3(lksi,i,j,k));

  for(i=0;i<d1(csum);i++)
    _(csum,i,0) += _(rl,i,0);
}


void update_bc(b,bsum,c,csum)
     TENSOR2 b,bsum,c,csum;
{
  pnormcol(b,bsum);
  pnormcol(c,csum);
}

void update_gradsum(gradsum,rv,rl,PCOND,DPCOND,U,Uniq,AUXPARAMS)
     TENSOR3 gradsum,PCOND;
     TENSOR2 rv,rl,U,Uniq;
     TENSOR5 DPCOND;
     auxparamstype AUXPARAMS;
{
  int e1,e2,ell,v,t,j,ut;
  for(e1=0;e1<NALF;e1++)
    for(e2=0;e2<NALF;e2++)
      for(ell=0;ell<NLOCAL;ell++)
	if(ell!=JUNKELL)
	  for(v=0;v<NVOCAB;v++)
	    for(t=0;t<d2(U);t++)
	      for(j=0;j<d2(Uniq);j++)
		if(WEQ(j,t)) {
		  ut = Ut(t);
		  _3(gradsum,e1,e2,ell) += _(rl,ell,t) * _(rv,v,j) * _5(DPCOND,ut,v,ell,e1,e2) / 
		    (_3(PCOND,ut,v,ell) + REALMIN);
		}
}


void update_pmp(pmp,gradsum,AUXPARAMS)
     TENSOR3 pmp,gradsum;
     auxparamstype AUXPARAMS;
{
  int e1,e2,ell,v,t,j;
  for(e1=0;e1<NALF;e1++)
    for(e2=0;e2<NALF;e2++)
      for(ell=0;ell<NLOCAL;ell++)
	if(ell!=JUNKELL)
	  _3(pmp,e1,e2,ell) = MAX(REALMIN,_3(pmp,e1,e2,ell)+AUXPARAMS.gradstep*_3(gradsum,e1,e2,ell));
  normpmp(pmp,AUXPARAMS);
}

 
void update_mu(mu,musum,eta,AUXPARAMS)
     TENSOR2 mu,musum,eta;
     auxparamstype AUXPARAMS;
{
  int a,v;
  for(v=0;v<NVOCAB;v++)
    for(a=0;a<NTOPIC;a++)
      _(mu,v,a) = MAX(EPS,_(eta,0,0)-1.0+_(musum,v,a));
  pnormcol(mu,mu);
}