/* This file contains a function R_mstep() called by R wraps m.step() in
   "R/fcn_m.step.r", and this function calls the relative functions mstep()
   in "src/emcluster.c" using in
   Dr. Maitra's ten-clusters programs.

   Writen: Wei-Chen Chen on 2008/10/19.
*/

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


/* This function is coded in "src/emcluster.c".
   Input:
     **X: double[n, p], data matrix of n*p.
     n: int[1], number of observations.
     p: int[1], number of dimersions.
     k: int[1], number of classes.		# nclass
     *pi: double[nclass], proportions of classes.
     **Mu: double[nclass, p], means of MVNs.
     **LTSigma: double[nclass, p * (p + 1) / 2], lower triangular sigma
                matrices.
     **Gamma: double[n, p], posterios matrix of n*p. 
   Output:
     *pi, **Mu, **LTSigma.
*/
void mstep(double **X,int n,int p,int k,double *pi,double **Mu,
           double **LTSigma,double **Gamma);

/* Allocate a pointer array with double precision. See "src/R_tool.c". */
double** allocate_double_array(int n);


/* This function calls mstep() in "src/emcluster.c" and is called by
   m.step() using .Call() in "R/fcn_m_step.r".
   Input:
     R_x: SEXP[R_n * R_p], data matrix of R_n*R_p.
     R_n: SEXP[1], number of observations.
     R_p: SEXP[1], number of dimersions.
     R_nclass: SEXP[1], number of classes.		# k
     R_Gamma: SEXP[R_n, R_p], posterios matrix of R_n*R_p. 
   Output:
     ret: a list contains
       pi: SEXP[R_nclass], proportions of classes.
       Mu: SEXP[R_nclass, R_p], means of MVNs.
       LTSigma: SEXP[R_nclass, R_p * (R_p + 1) / 2], lower triangular sigma
                matrices.
*/
SEXP R_mstep(SEXP R_x, SEXP R_n, SEXP R_p, SEXP R_nclass, SEXP R_Gamma){
  /* Declare variables for calling C. */
  double **C_Gamma, **C_x, *C_pi, **C_Mu, **C_LTSigma;
  int *C_n, *C_p, *C_nclass;

  /* Declare variables for R's returning. */
  SEXP pi, Mu, LTSigma, ret, ret_names;

  /* Declare variables for processing. */
  double *tmp_1, *tmp_2;
  int i, p_LTSigma;
  char *names[3] = {"pi", "Mu", "LTSigma"};

  /* Set initial values. */
  C_n = INTEGER(R_n);
  C_p = INTEGER(R_p);
  C_nclass = INTEGER(R_nclass);
  p_LTSigma = *C_p * (*C_p + 1) / 2;

  /* Allocate and protate storages. */
  PROTECT(pi = allocVector(REALSXP, *C_nclass));
  PROTECT(Mu = allocVector(REALSXP, *C_nclass * *C_p));
  PROTECT(LTSigma = allocVector(REALSXP, *C_nclass * p_LTSigma));
  PROTECT(ret = allocVector(VECSXP, 3));
  PROTECT(ret_names = allocVector(STRSXP, 3));

  i = 0;
  SET_VECTOR_ELT(ret, i++, pi);
  SET_VECTOR_ELT(ret, i++, Mu);
  SET_VECTOR_ELT(ret, i++, LTSigma);
  
  for(i = 0; i < 3; i++){
    SET_STRING_ELT(ret_names, i, mkChar(names[i])); 
  }
  setAttrib(ret, R_NamesSymbol, ret_names);

  /* Assign data. */
  C_Gamma = allocate_double_array(*C_n);
  C_x = allocate_double_array(*C_n);
  C_Mu = allocate_double_array(*C_nclass);
  C_LTSigma = allocate_double_array(*C_nclass);

  tmp_1 = REAL(R_Gamma);
  tmp_2 = REAL(R_x);
  for(i = 0; i < *C_n; i++){
    C_Gamma[i] = tmp_1;
    C_x[i] = tmp_2;
    tmp_1 += *C_nclass;
    tmp_2 += *C_p;
  }

  tmp_1 = REAL(Mu);
  tmp_2 = REAL(LTSigma);
  for(i = 0; i < *C_nclass; i++){
    C_Mu[i] = tmp_1;
    C_LTSigma[i] = tmp_2;
    tmp_1 += *C_p;
    tmp_2 += p_LTSigma;
  }

  C_pi = REAL(pi);

  /* Compute. */
  mstep(C_x, *C_n, *C_p, *C_nclass, C_pi, C_Mu, C_LTSigma, C_Gamma);

  /* Free memory and release protectation. */
  free(C_Gamma);
  free(C_x);
  free(C_Mu);
  free(C_LTSigma);
  UNPROTECT(5);

  return(ret);
} /* End of R_mstep(). */