#include <gsl/gsl_rng.h>

#define MSG(str) printf(str);

#define Ut(t) (((int)(_(U,0,t)))-1)
/*#define ui(t) (((int)(_(u,0,t)))-1)
  #define vj(t) (((int)(_(v,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 RAND gsl_rng_uniform(RNG)
#define NVOC mxGetM(RAWVOC)
#define NALF d1(pmp)

#define WEQ(j,t) ( ABS(_(Uniq,0,j)-_(U,0,t)) < EPS)
/* P(rimitive) M(orphological) P(robability) */
/*#define PMP(a,b) pmpfunc((a),(b),(ell),(pmp))*/
/*#define RHO(a) rhofunc(a,ell,rho)*/
/*#define PMP(a,b) pmpfunc(alfinds[a-1],alfinds[b-1],(ell),(pmp))
  #define RHO(a) rhofunc(alfinds[a-1],ell,rho)*/
#define PMP(a,b) _3(pmp,alfinds[a-1],alfinds[b-1],ell)
#define RHO(a) _(rho,alfinds[a-1],ell)
#define NULLCHAR ('@')
#define JUNKELL 0


gsl_rng *init_rand(unsigned long int randseed)
{
  const gsl_rng_type *T;
  gsl_rng *RNG;
  gsl_rng_env_setup();
  T = gsl_rng_default;
  RNG = gsl_rng_alloc (T);
  gsl_rng_set(RNG,randseed);
  return RNG;
}


void randTENSOR2(TENSOR2 x,gsl_rng *RNG)
{
  int i,j;
  FOR2(i,j,d1(x),d2(x)) _(x,i,j) = RAND;
}

void randTENSOR3(TENSOR3 x,gsl_rng *RNG)
{
  int i,j,k;
  FOR3(i,j,k,d1(x),d2(x),d3(x)) _3(x,i,j,k) = RAND;
}


double pmpfunc(a,b,ell,pmp,rho)
     int ell;
     int a,b;
     TENSOR3 pmp;
{
  return(_3(pmp,a,b,ell));
}


void normpmp(TENSOR3 pmp)
{
  int ell,i,j;
  double s;
  for(ell=0;ell<d3(pmp);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;
  }
}


double rhofunc(a,ell,rho)
     int ell;
     int a;
     TENSOR2 rho;
{
  return(_(rho,a,ell));
}


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


void jointstrprob(psi,u,v,ell,pmp,alfinds)
     /* assumes psi has been allocated to the right size */
     TENSOR2 psi;
     char *u,*v;
     int ell;
     TENSOR3 pmp;
     int *alfinds;
{
  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,alfinds)
     /* assumes psi,psi1,A,B have been allocated to the right size */
     /* B = inv(I-A) */
     TENSOR2 psm,psi1,A,B,rho;
     char *u;
     int ell;
     TENSOR3 pmp;
     int *alfinds;
{
  int i,j;
  TENSOR2 I_A; /* I-A */
  mxArray *IAar;

  IAar = allocTENSOR2(&I_A,d1(B),d2(B)); /* must free later */

  /* 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);
  mxDestroyArray(IAar);
}


double condstrprob(u,v,ell,pmp,rho,RAWVOC,alfinds)
     TENSOR2 rho;
     int u,v;
     int ell;
     TENSOR3 pmp;
     const mxArray *RAWVOC;
     int *alfinds;
{
  TENSOR2 psm,psj,ps1,A,B;
  char *ustr,*vstr;
  int b;
  double prob;
  mxArray *psmAr,*psjAr,*ps1Ar,*BAr,*AAr;
  ustr = mxCalloc(UVLEN(u)+1,sizeof(mxChar));
  vstr = mxCalloc(UVLEN(v)+1,sizeof(mxChar));
  /* int mxGetString(const mxArray *array_ptr, char *buf, int buflen); */
  /* really ought to do some error trapping here, but won't bother */
  b = mxGetString(UVARR(u),ustr,UVLEN(u)*sizeof(mxChar)+1);
  b = mxGetString(UVARR(v),vstr,UVLEN(v)*sizeof(mxChar)+1);
  if(ell==JUNKELL)
    prob = REALMIN + streq(ustr,vstr,UVLEN(u),UVLEN(v));
  else {
    psjAr = allocTENSOR2(&psj,UVLEN(v)+1,UVLEN(u)+1);
    psmAr = allocTENSOR2(&psm,UVLEN(v)+1,1);
    ps1Ar = allocTENSOR2(&ps1,UVLEN(v)+1,1);
    AAr = allocTENSOR2(&A,UVLEN(v)+1,UVLEN(v)+1);
    BAr = allocTENSOR2(&B,UVLEN(v)+1,UVLEN(v)+1);
    jointstrprob(psj,vstr,ustr,ell,pmp,alfinds);
    margstrprob(psm,ps1,A,B,vstr,ell,pmp,rho,alfinds);
    prob = _(psj,d1(psj)-1,d2(psj)-1) / (_(ps1,d1(ps1)-1,0)+REALMIN);
    mxDestroyArray(psjAr);
    mxDestroyArray(psmAr);
    mxDestroyArray(ps1Ar);
    mxDestroyArray(AAr);
    mxDestroyArray(BAr);
  }
  mxFree(ustr);
  mxFree(vstr);
  return prob;
}


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

void DpsmDnu(dpsm,u,e1,e2,psm,ell,pmp,alfinds)
     /* assumes dpsm has been allocated to the right size */
     TENSOR2 dpsm,psm;
     char *u,e1,e2;
     int ell;
     TENSOR3 pmp;
     int *alfinds;
{
  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,0) * PMPu0;
}


void DpsiDnu(dpsi,psi,u,v,e1,e2,ell,pmp,alfinds)
     /* assumes dpsi has been allocated to the right size */
     TENSOR2 dpsi,psi;
     char *u,*v,e1,e2;
     int ell;
     TENSOR3 pmp;
     int *alfinds;
{
  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,alfinds)
     TENSOR2 A,B,dA,dB,rho;
     char *u,e1,e2;
     int ell;
     TENSOR3 pmp;
     int *alfinds;
{
  int i,j,k,l;
  /* compute dA */
  for(i=0;i<d1(dA);i++) {
    _(dA,i,i) = DRhoDnu(NULLCHAR,e1,e2);
    if(i) _(dA,i,i-1) = DRhoDnu(ui(i-1),e1,e2) + DRhoDnu(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 */
  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(u,v,ell,e1,e2,pmp,rho,RAWVOC,alfinds)
     TENSOR2 rho;
     int u,v;
     char e1,e2;
     int ell;
     TENSOR3 pmp;
     const mxArray *RAWVOC;
     int *alfinds;
{
  TENSOR2 psm,psj,ps1,A,B,dpsm,dpsj,dA,dB;
  char *ustr,*vstr;
  int j;
  double pjo,pma,dpjo,dpma;
  mxArray *psmAr,*psjAr,*ps1Ar,*BAr,*AAr,*dpsmAr,*dpsjAr,*dAAr,*dBAr;
  if(ell==JUNKELL)
    return(REALMIN);
  else {
    /* memory allocation bureaucracy */
    ustr = mxCalloc(UVLEN(u)+1,sizeof(mxChar));
    vstr = mxCalloc(UVLEN(v)+1,sizeof(mxChar));
    j = mxGetString(UVARR(u),ustr,UVLEN(u)*sizeof(mxChar)+1);
    j = mxGetString(UVARR(v),vstr,UVLEN(v)*sizeof(mxChar)+1);
    psjAr = allocTENSOR2(&psj,UVLEN(v)+1,UVLEN(u)+1);
    dpsjAr = allocTENSOR2(&dpsj,UVLEN(v)+1,UVLEN(u)+1);
    psmAr = allocTENSOR2(&psm,UVLEN(v)+1,1);
    dpsmAr = allocTENSOR2(&dpsm,UVLEN(v)+1,1);
    ps1Ar = allocTENSOR2(&ps1,UVLEN(v)+1,1);
    AAr = allocTENSOR2(&A,UVLEN(v)+1,UVLEN(v)+1);
    BAr = allocTENSOR2(&B,UVLEN(v)+1,UVLEN(v)+1);
    dAAr = allocTENSOR2(&dA,UVLEN(v)+1,UVLEN(v)+1);
    dBAr = allocTENSOR2(&dB,UVLEN(v)+1,UVLEN(v)+1);
    /* actual calculations */
    jointstrprob(psj,vstr,ustr,ell,pmp,alfinds);
    pjo = _(psj,d1(psj)-1,d2(psj)-1); /* joint prob */
    margstrprob(psm,ps1,A,B,vstr,ell,pmp,rho,alfinds);
    pma = _(ps1,d1(ps1)-1,0); /* marg prob */
    DpsiDnu(dpsj,psj,vstr,ustr,e1,e2,ell,pmp,alfinds);
    dpjo = _(dpsj,d1(dpsj)-1,d2(dpsj)-1); /* diff joint prob */
    DpsmDnu(dpsm,vstr,e1,e2,psm,ell,pmp,alfinds);
    DABDnu(dA,dB,vstr,e1,e2,A,B,ell,pmp,rho,alfinds);
    dpma = 0.0; /* diff marg prob */
    for(j=0;j<d1(psm);j++)
      dpma += _(dB,d1(dB)-1,j)*_(psm,j,0) + _(B,d1(dB)-1,j)*_(dpsm,j,0);
    /* free memory */
    mxDestroyArray(psjAr);
    mxDestroyArray(dpsjAr);
    mxDestroyArray(psmAr);
    mxDestroyArray(dpsmAr);
    mxDestroyArray(ps1Ar);
    mxDestroyArray(AAr);
    mxDestroyArray(BAr);
    mxDestroyArray(dAAr);
    mxDestroyArray(dBAr);
    mxFree(ustr);
    mxFree(vstr);
    return((pma*dpjo - pjo*dpma)/(REALMIN+pma)/(REALMIN+pma)); /* diff quotient rule */
  }
}


void buildPCOND(PCOND,RAWVOC,pmp,rho,alfinds)
     TENSOR3 PCOND,pmp;
     TENSOR2 rho;
     const mxArray *RAWVOC;
     int *alfinds;
{
  int w,v,ell;
  for(w=0;w<d1(PCOND);w++)
    for(v=0;v<d2(PCOND);v++)
      for(ell=0;ell<d3(PCOND);ell++)
	_3(PCOND,w,v,ell) = condstrprob(w,v,ell,pmp,rho,RAWVOC,alfinds);
}


void buildPWORD(PWORD,PCOND,RAWVOC,mu,pmp,rho,alfinds)
     TENSOR3 PWORD,PCOND,pmp;
     TENSOR2 mu,rho;
     const mxArray *RAWVOC;
     int *alfinds;
{
  int w,v,a,ell;
  for(w=0;w<d1(PWORD);w++)
    for(ell=0;ell<d2(PWORD);ell++)
      for(a=0;a<d3(PWORD);a++) {
	_3(PWORD,w,ell,a) = 0.0;
	for(v=0;v<d1(PWORD);v++)
	  _3(PWORD,w,ell,a) += _(mu,v,a) * _3(PCOND,w,v,ell);
      }
}


void buildDPCOND(DPCOND,RAWVOC,pmp,rho,ALF,alfinds)
     TENSOR3 pmp;
     TENSOR2 rho;
     TENSOR5 DPCOND;
     const mxArray *RAWVOC;
     char *ALF;
     int *alfinds;
{
  int w,v,ell,ei,ej;
  char a,b;
  for(w=0;w<d1(DPCOND);w++)
    for(v=0;v<d2(DPCOND);v++)
      for(ell=0;(ell<d3(DPCOND))&&(ell!=JUNKELL);ell++)
	for(ei=0;ei<NALF;ei++)
	  for(ej=0;ej<NALF;ej++) {
	    a = ALF[ei];
	    b = ALF[ej];
	    _5(DPCOND,w,v,ell,ei,ej) = dPidNu(w,v,ell,a,b,pmp,rho,RAWVOC,alfinds);
      }
}


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)
     TENSOR2 U,qa;
     TENSOR3 PWORD;
     mxArray *lB;
{
  double s;
  int Uj,j,a,ell;
  for(j=0;j<d2(U);j++) {
    Uj = Ut(j);
    for(ell=0;ell<mxGetN(lB);ell++) {
      s = 0.0;
      for(a=0;a<d1(qa);a++)
	s += _(qa,a,j)*LOG(_3(PWORD,Uj,ell,a));
      set_se(lB,Uj,ell,s);
    }
  }
}


void markovBr(U,qa,rv,PCOND,mu,lB)
     TENSOR2 U,qa,rv,mu;
     TENSOR3 PCOND;
     mxArray *lB;
{
  double sa,sv;
  int Uj,j,a,ell,v;
  for(j=0;j<d2(U);j++) {
    Uj = Ut(j);
    for(ell=0;ell<mxGetN(lB);ell++) {
      sa = 0.0;
      for(a=0;a<d1(qa);a++) {
	sv = 0.0;
	for(v=0;v<mxGetM(lB);v++)
	  sv += _(rv,v,j)*LOG(_(mu,a,v)*_3(PCOND,Uj,v,ell));
	sa += _(qa,a,j)*sv;
      }
      set_se(lB,Uj,ell,sa);
    }
  }
}



double logdirint(alpha,qa,n)
     TENSOR2 alpha,qa;
     int n;
{
  double s=0.0,sa=0.0;
  int a;
  for(a=0;a<d1(qa);a++) {
    sa += _(alpha,a,0);
    s += LGAM(_(alpha,a,0)+sumrow(qa,a)) - LGAM(_(alpha,a,0));
  }
  return(s + LGAM(sa) - LGAM(sa+n));
}


void lbloglik(ll,lfp,U,Uniq,Ucnt,qa,lA,lB,lC,alpha,lalf)
     TENSOR2 U,Uniq,Ucnt,qa,lA,lC,alpha,lalf;
     mxArray *lB;
     double *lfp,*ll;
{
  int a,j;
  *ll = 0.0; /*_(ll,0,0)=0.0;*/
  for(a=0;a<d1(qa);a++)
    for(j=0;j<d2(Uniq);j++)
      *ll -= _(qa,a,j)*_(Ucnt,0,j)*LOG(_(qa,a,j));
  *lfp = logsumexp_col(lalf,d2(lalf)-1);
  *ll += *lfp;
  *ll += logdirint(alpha,qa,d2(U));
}


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;
}


void alphabeta_lg(lalf,lbet,lA,lB,lC,U,tmp)
     TENSOR2 U,lA,lC,lalf,lbet,tmp;
     mxArray *lB;
{
  int i,j,t;

  for(i=0;i<d1(lalf);i++) _(lalf,i,0) = val_se1(lB,Ut(0),i) + _(lC,i,0);
  for(t=1;t<d2(U);t++) {
    for(i=0;i<d1(lA);i++) {
      for(j=0;j<d1(lA);j++)
	_(tmp,j,0) = _(lalf,j,t-1) + _(lA,i,j) + val_se1(lB,Ut(t),i);
      _(lalf,i,t) = logsumexp_col(tmp,0);
    }
  }

  for(i=0;i<d1(lbet);i++) _(lbet,i,d2(U)-1) = 0.0;
  for(t=d2(U)-1;t>0;t--) {
    for(i=0;i<d1(lA);i++) {
      for(j=0;j<d1(lA);j++)
	_(tmp,j,0) = _(lbet,j,t) + _(lA,j,i) + val_se1(lB,Ut(t),j);
      _(lbet,i,t-1) = logsumexp_col(tmp,0);
    }
  }

}


void update_rvrl(rv1,rl1,rv0,rl0,qa,U,Uniq,Ucnt,lA,lC,PCOND,mu,lB,tmp,lalf,lbet)
     TENSOR2 U,Uniq,Ucnt,lA,lC,mu,rv1,rl1,rv0,rl0,qa;
     TENSOR2 lalf,lbet,tmp;
     TENSOR3 PCOND;
     mxArray *lB;
{
  int v,a,j,ell,t,Uj;
  double tt,s,x,ga;

  markovBr(Uniq,qa,rv0,PCOND,mu,lB);
  alphabeta_lg(lalf,lbet,lA,lB,lC,U,tmp); 
  for(ell=0;ell<d1(rl1);ell++)
    for(j=0;j<d2(rl1);j++)
      _(rl1,ell,j) = EXP( _(lalf,ell,j) + _(lbet,ell,j));
  pnormcol(rl1,rl1);

  for(v=0;v<mxGetM(lB);v++) {
    for(j=0;j<d2(Uniq);j++) {
      Uj = ((int)(_(Uniq,0,j)))-1;
      tt = _(Ucnt,0,j);
      s=0.0;
      for(t=0;t<d2(U);t++)
	if(Uj==Ut(t)) {
	  x = 1.0;
	  for(ell=0;ell<d1(rl1);ell++)
	    for(a=0;a<d1(qa);a++)
	      x *= POW((_(mu,v,a) * _3(PCOND,Uj,v,ell)),(_(qa,a,j)*_(rl1,ell,t)/tt));
	  s += x;
	}
      _(rv1,v,j) = s;
    }
  }
  pnormcol(rv1,rv1);
}



void update_qaql(qa1,ql1,qa0,ql0,U,Uniq,Ucnt,lA,lC,PWORD,alpha,lB,tmp,lalf,lbet)
     TENSOR2 U,Uniq,Ucnt,lA,lC,alpha,qa1,ql1,qa0,ql0;
     TENSOR2 lalf,lbet,tmp;
     TENSOR3 PWORD;
     mxArray *lB;/*,*qa1,*ql1,*qa0,*ql0;*/
{
  int a,j,ell,t;
  double tt,s,x,ga;

  markovBq(Uniq,qa0,PWORD,lB);
  alphabeta_lg(lalf,lbet,lA,lB,lC,U,tmp); 
  for(ell=0;ell<d1(ql1);ell++)
    for(j=0;j<d2(ql1);j++)
      _(ql1,ell,j) = EXP( _(lalf,ell,j) + _(lbet,ell,j));
  pnormcol(ql1,ql1);


  for(a=0;a<d1(qa1);a++) {
    ga = _(alpha,a,0) + sumrow(qa0,a);
    for(j=0;j<d2(Uniq);j++) {
      tt = _(Ucnt,0,j);
      s=0.0;
      for(t=0;t<d2(U);t++)
	if(WEQ(j,t)) {
	  x=1.0;
	  for(ell=0;ell<d1(ql1);ell++)
	    x *= POW(_3(PWORD,j,ell,a),_(ql1,ell,t)/tt);
	  s += x;
	}
      _(qa1,a,j) = EXP(PSI(ga))*s;
    }
  }
  pnormcol(qa1,qa1);
}



int update_rloop(rv1,rl1,rv0,rl0,qa,U,Uniq,Ucnt,lA,lC,PCOND,mu,lB,tmp,lalf,lbet,tol,infrloop)
     TENSOR2 U,Uniq,Ucnt,*rv1,*rl1,*rv0,*rl0,qa,lA,lC,mu;
     TENSOR2 lalf,lbet,tmp;
     TENSOR3 PCOND;
     mxArray *lB;
     double tol;
     int infrloop;
{
  int i,good;
  good = 0;
  for(i=0;(i<infrloop)&&(!good);i++) {
    swapTENSOR2(rv0,rv1);
    swapTENSOR2(rl0,rl1);
    update_rvrl(*rv1,*rl1,*rv0,*rl0,qa,U,Uniq,Ucnt,lA,lC,PCOND,mu,lB,tmp,lalf,lbet);
    good = (l1diff2(*rl1,*rl0)<tol);
  }
  return good;
}



int update_qloop(qa1,ql1,qa0,ql0,U,Uniq,Ucnt,lA,lC,PWORD,alpha,lB,tmp,lalf,lbet,tol,infrloop)
     TENSOR2 U,Uniq,Ucnt,*qa1,*ql1,*qa0,*ql0,lA,lC,alpha;
     TENSOR2 lalf,lbet,tmp;
     TENSOR3 PWORD;
     mxArray *lB;
     double tol;
     int infrloop;
{
  int i,good;
  good = 0;
  for(i=0;(i<infrloop)&&(!good);i++) {
    swapTENSOR2(qa0,qa1);
    swapTENSOR2(ql0,ql1);
    update_qaql(*qa1,*ql1,*qa0,*ql0,U,Uniq,Ucnt,lA,lC,PWORD,alpha,lB,tmp,lalf,lbet);
    good = (l1diff2(*qa1,*qa0)<tol);
  }
  return good;
}


void inference_q(qa,qa0,ql,ql0,U,Uniq,Ucnt,lA,lC,PWORD,alpha,lB,tmp,lalf,lbet,tol,infrloop,RNG)
     TENSOR2 U,Uniq,Ucnt,lA,lC,alpha,*qa,*ql,*qa0,*ql0;
     TENSOR2 lalf,lbet,tmp;
     TENSOR3 PWORD;
     mxArray *lB;
     double tol;
     int infrloop;
     gsl_rng *RNG;
{
  randTENSOR2(*qa,RNG); pnormcol(*qa,*qa);
  update_qloop(qa,ql,qa0,ql0,U,Uniq,Ucnt,lA,lC,PWORD,alpha,lB,tmp,lalf,lbet,tol,infrloop);
}


void inference_r(rv,rv0,rl,rl0,qa,U,Uniq,Ucnt,lA,lC,PCOND,mu,lB,tmp,lalf,lbet,tol,infrloop,RNG)
     TENSOR2 U,Uniq,Ucnt,lA,lC,mu,*rv,*rl,*rv0,*rl0,qa;
     TENSOR2 lalf,lbet,tmp;
     TENSOR3 PCOND;
     mxArray *lB;
     double tol;
     int infrloop;
     gsl_rng *RNG;
{
  randTENSOR2(*rv,RNG); pnormcol(*rv,*rv);
  update_rloop(rv,rl,rv0,rl0,qa,U,Uniq,Ucnt,lA,lC,PCOND,mu,lB,tmp,lalf,lbet,tol,infrloop);
}


void init_alpha(alpha)
     TENSOR2 alpha;
{
  int a;
  for(a=0;a<d1(alpha);a++)
    _(alpha,a,0) = 1.0;
}

void init_eta(eta)
     TENSOR2 eta;
{
  _(eta,0,0) = 0.01;
}


void init_nu(pmp,RNG,ALF,alfinds)
     /* trivial morph ed probs */
     int *alfinds;
     gsl_rng *RNG;
     TENSOR3 pmp;
     char *ALF;
{
  int ei,ej,ell;
  char a,b;
  for(ei=0;ei<NALF;ei++) {
    a = ALF[ei];
    for(ej=0;ej<NALF;ej++) {
      b = ALF[ej];
      for(ell=0;ell<d3(pmp);ell++)
	PMP(a,b) = 1.0*(a==b) + 1e-4*RAND;
    }
  }
  PMP00 = 10.0;
  normpmp(pmp);
}

double zerofunc()
{
     return 0.0;
}

void update_musum(musum,Uniq,Ucnt,qa,rv)
     TENSOR2 musum,Uniq,Ucnt,qa,rv;
{
  int v,a,j;
  for(v=0;v<d1(musum);v++)
    for(a=0;a<d2(musum);a++)
      for(j=0;j<d2(Uniq);j++)
	_(musum,v,a) += _(Ucnt,0,j) * _(rv,v,j) * _(qa,a,j);
}

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)
     TENSOR3 gradsum,PCOND;
     TENSOR2 rv,rl,U,Uniq;
     TENSOR5 DPCOND;
{
  int e1,e2,ell,v,t,j;
  for(e1=0;e1<d1(gradsum);e1++)
    for(e2=0;e2<d2(gradsum);e2++)
      for(ell=0;(ell<d3(gradsum))&&(ell!=JUNKELL);ell++)
	for(v=0;v<d1(PCOND);v++)
	  for(t=0;t<d2(U);t++)
	    for(j=0;j<d2(Uniq);j++)
	      if(WEQ(j,t))
		_3(gradsum,e1,e2,ell) += _(rl,ell,t) * _(rv,v,j) * _5(DPCOND,Ut(t),v,ell,e1,e2) / 
		  (_3(PCOND,Ut(t),v,ell) + REALMIN);
}

void update_pmp(pmp,gradsum,gradstep)
     TENSOR3 pmp,gradsum;
{
  int e1,e2,ell,v,t,j;
  for(e1=0;e1<d1(gradsum);e1++)
    for(e2=0;e2<d2(gradsum);e2++)
      for(ell=0;(ell<d3(gradsum))&&(ell!=JUNKELL);ell++)
	_3(pmp,e1,e2,ell) = MAX(REALMIN,_3(pmp,e1,e2,ell)+gradstep*_3(gradsum,e1,e2,ell));
  normpmp(pmp);
}

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


void lda_learn(randseed,true_loglik,infrloop,infrtol,gradloop,learnloop,gradstep,RAWVOC,CORPUS,UNIQ,UCNTS,alpha,mu,eta,b,c,pmp,ALF,alfinds)
     unsigned long int randseed;
     double true_loglik,gradstep,infrtol;
     int infrloop,gradloop,learnloop,*alfinds;
     const mxArray *RAWVOC,*CORPUS,*UNIQ,*UCNTS;
     TENSOR2 alpha,mu,eta,b,c;
     TENSOR3 pmp;
     char *ALF;
{
  /* input parameters
     -randseed
     -ground truth loglik
     _Nlocal (note: Nlocal implicitly defined by dims of b)
     _Ntopic (note: Ntopic implicitly defined by dims of alpha)
     -infrloop (how much to loop for inference)
     -infrtol (to what precision to carry out inference)
     -gradloop (how many gradient steps to take -- 1 for now)
     -learnloop (how much to loop for learning)
     -gradstep (gradient step size)
     -alphabet

     data
     -RAWVOC
     -CORPUS (numeric)
     -UNIQ (numeric)
     -UCNTS

     output parameters
     -alpha (will hold fixed at 1 for now)
     -eta (will hold fixed at small ~ 1/1000 for now)
     -mu (key parameter to estimate)
     -b (will estimate)
     -c (will estimate)
     -nu (key parameter to estimate -- pmp)
  */

  int d,i,j;
  double ll,LL,lfp;
  gsl_rng *RNG;
  mxArray *PWORDar,*PCONDar,*DPCONDar,*lB;
  mxArray *OUTAUX[11],*INNAUX[9];
  TENSOR2 U,Uniq,Ucnt,rho,lA,lC,musum,bsum,csum,lalf,lbet,tmp,qa,qa0,qrl,qrl0,rv,rv0;
  TENSOR3 PWORD,PCOND,gradsum,lksi;
  TENSOR5 DPCOND;

  /* allocate auxilliary memory */
  OUTAUX[ 0] = allocTENSOR3(&PWORD,NVOC,d3(pmp),d2(mu));
  OUTAUX[ 1] = allocTENSOR3(&PCOND,NVOC,NVOC,d3(pmp));
  OUTAUX[ 2] = allocTENSOR5(&DPCOND,NVOC,NVOC,d3(pmp),d1(pmp),d2(pmp));
  OUTAUX[ 3] = allocTENSOR2(&rho,d1(pmp),d3(pmp));
  OUTAUX[ 4] = allocTENSOR2(&lA,d1(b),d2(b));
  OUTAUX[ 5] = allocTENSOR2(&lC,d1(c),d2(c));
  OUTAUX[ 6] = allocTENSOR2(&musum,d1(mu),d2(mu));
  OUTAUX[ 7] = allocTENSOR2(&bsum,d1(b),d2(b));
  OUTAUX[ 8] = allocTENSOR2(&csum,d1(c),d2(c));
  OUTAUX[ 9] = allocTENSOR3(&gradsum,NALF,NALF,d3(pmp));
  OUTAUX[10] = allocTENSOR2(&tmp,d1(b),1);

  RNG = init_rand(randseed);   /* init random number generator (RNG) */
  init_alpha(alpha);
  init_eta(eta);
  /* init mu */
  randTENSOR2(mu,RNG); pnormcol(mu,mu);
  /* init b,c */
  randTENSOR2(b,RNG); pnormcol(b,b);
  randTENSOR2(c,RNG); pnormcol(c,c);
  init_nu(pmp,RNG,ALF,alfinds);

  for(i=0;i<learnloop;i++) {
    init_rho(rho,pmp,ALF,alfinds);
    LL = 0.0;
    printf("starting iteration %d of %d \n",i+1,learnloop);
    applyfuncTENSOR2(lA,b,&LOG);
    applyfuncTENSOR2(lC,c,&LOG);
    applyfuncTENSOR2(musum,musum,&zerofunc); /*init_musum;*/
    applyfuncTENSOR2(bsum,bsum,&zerofunc); /*init_bsum;*/
    applyfuncTENSOR2(csum,csum,&zerofunc); /*init_csum;*/
    applyfuncTENSOR3(gradsum,gradsum,&zerofunc); /*init_gradsum;*/
    /*printf("building PCOND... ");*/
    buildPCOND(PCOND,RAWVOC,pmp,rho,alfinds);
    /*printf("building PWORD...");*/
    buildPWORD(PWORD,PCOND,RAWVOC,mu,pmp,rho,alfinds);
    /*printf("building DPCOND...\n");*/
    buildDPCOND(DPCOND,RAWVOC,pmp,rho,ALF,alfinds);
    /*printf("processing corpus: \n");*/
    for(d=0;d<mxGetM(CORPUS);d++) {
      /*printf("load doc %d of %d; ",d+1,mxGetM(CORPUS));*/
      makeT2ad(U,    mxGetCell(CORPUS,d));
      makeT2ad(Uniq, mxGetCell(UNIQ,d));
      makeT2ad(Ucnt, mxGetCell(UCNTS,d));
      /* allocate local auxilliary variables */
      lB = mxCreateSparse(NVOC,d1(b),d2(U)*d1(b),mxREAL);
      INNAUX[0] = allocTENSOR2(&lalf,d1(b),d2(U));
      INNAUX[1] = allocTENSOR2(&lbet,d1(b),d2(U));
      INNAUX[2] = allocTENSOR2(&qa,d2(mu),d2(Uniq));
      INNAUX[3] = allocTENSOR2(&qa0,d2(mu),d2(Uniq));
      INNAUX[4] = allocTENSOR2(&rv,NVOC,d2(Uniq));
      INNAUX[5] = allocTENSOR2(&rv0,NVOC,d2(Uniq));
      INNAUX[6] = allocTENSOR2(&qrl,d1(b),d2(U));
      INNAUX[7] = allocTENSOR2(&qrl0,d1(b),d2(U));
      INNAUX[8] = allocTENSOR3(&lksi,d1(b),d1(b),d2(U)-1);
      /* perform inference on dth document */
      inference_q(&qa,&qa0,&qrl,&qrl0,U,Uniq,Ucnt,lA,lC,PWORD,alpha,lB,tmp,lalf,lbet,infrtol,infrloop,RNG);
      /*
      inference_r(&rv,&rv0,&qrl,&qrl0,qa,U,Uniq,Ucnt,lA,lC,PCOND,mu   ,lB,tmp,lalf,lbet,infrtol,infrloop,RNG);
      */
      lbloglik(&ll,&lfp,U,Uniq,Ucnt,qa,lA,lB,lC,alpha,lalf);
      LL += ll;
      compute_lksi(lksi,lalf,lbet,lA,lB,U,lfp);
      /* gather statistics */
      update_musum(musum,Uniq,Ucnt,qa,rv);
      update_bcsum(bsum,csum,lksi,qrl);
      update_gradsum(gradsum,rv,qrl,PCOND,DPCOND);
      /* free local auxilliary variables */
      mxDestroyArray(lB);
      for(j=0;j<9;j++)
	mxDestroyArray(INNAUX[j]);
    }
    /* update_alpha; */
    /* update_eta; */
    /* update_mu(mu,musum,eta); */
    update_bc(b,bsum,c,csum);
    /* update_pmp(pmp,gradsum,gradstep); */
    printf("loglik: %1.5f; \n",LL);
  }

  /* free auxilliary memory */ 
  /*printf("free aux \n "); */
  for(j=0;j<11;j++)
    mxDestroyArray(OUTAUX[j]);
  
  gsl_rng_free(RNG); /* free random number generator (RNG) */
}