#include <R.h>
#include <Rmath.h>
#include <Rinternals.h>
#include <stdio.h>

/* simultaneously calculate m different products of KDEs, 1 for each component,
   as in equation (8) of Benaglia et al for a fixed value of \ell.  
   If r is the number of coordinates in block \ell, then each
   final value is the r-fold product of KDEs */
void KDErepeatedbw(
     int *nn, /* Sample size */
     int *mm, /* Number of components */
     int *rr, /* size of current block */
     double *x, /* data:  vector of length nn*rr */
     double *hh, /* m-vector of bandwidths (compare to KDErepeated) */
     double *z, /* nn*mm vector of normalized posteriors (or indicators in
                   stochastic case), normalized by "column" */
     double *f  /* nxm matrix of KDE products */
) {
  int n=*nn, i, ii;
  int mn = *mm*n, rn=*rr*n, jn, kn, kkn;
  double sum1, sum2, tmp, xik;
  double const2 = 0.39894228040143267794/((double)(*rr)); /* .3989...=1/(sqrt(2*pi)) */
  double const1;

  for(jn=0; jn<mn; jn+=n, hh++) { /* jn is component index times n */
    /* at each iteration, *hh changes to the current bandwidth value */
    const1 = -0.5 / (*hh * *hh);
    for(i=0; i<n; i++) {
      f[i + jn] = 1.0;
      for(kn=0; kn<rn; kn+=n) { /* kn is coordinate index within block times n */
        sum1 = 0.0;
        xik = x[i + kn];
        for(ii=0; ii<n; ii++) { /* ii is the i index of equation (8) */
          sum2 = 0.0;
          for(kkn=0; kkn<rn; kkn+=n) { /* kk is the k index of equation (8) */
            tmp = xik - x[ii + kkn];
            sum2 +=  exp(tmp * tmp * const1); /* Using normal kernel */
          }
          sum1 += z[ii + jn] * sum2;
        }
        f[i + jn] *= sum1 * const2 / *hh;
      }
    }
  }
}