// MBCategMean.c
// Calculate the mean vector/matrix for each category.
//
//   The algorithm for MicroBayes is described in the following paper
//
//   D. Zhang, M.T. Wells, C.D. Smart and W.E. Fry (2003). Bayesian
//       Normalization and Identification for Differential Gene
//       Expression Data.
//
//   PLASE CITE THIS PAPER AFTER YOU USE THIS SOFTWARE.
//
//	Dabao Zhang		February 24, 2003
//	Revised by Dabao Zhang
//       March 9, 2003 -- Consider genes' location effect
//
//	Copyright (c) 2003 by Dabao Zhang.

#include <math.h>
#include "mex.h"

void MBCategMean(int nrows,double inM[],double inCateg[],double mM[])
{
    int n, k;
    int idxS, idxCurr;
    double mCurr;
    
    idxS = 0;
    idxCurr = (int)inCateg[0];
    mCurr = 0;
    for(n=0; n<nrows; n++)
    {
        if( idxCurr != ((int)inCateg[n]) )
        {
            mM[idxCurr-1] = mCurr/(n-idxS);

            mCurr = inM[n];
            idxS = n;
            idxCurr = (int)inCateg[n];
        }
        else
        {
            mCurr = mCurr + inM[n];
        }

        if( n==(nrows-1) )
        {
            mM[idxCurr-1] = mCurr/(n-idxS+1);
        }
    }
}

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    double *inM,*inCateg,*mM;
    int nrows;

    // Check for proper number of arguments.
    if( nrhs!=2 )
    {
        mexErrMsgTxt("Two inputs required.");
    }
    else if( nlhs>1 )
    {
        mexErrMsgTxt("Too many output arguments.");
    }

    nrows = mxGetM(prhs[0]);
    
    // RHS: inM, inCateg
    inM = mxGetPr(prhs[0]);
    inCateg = mxGetPr(prhs[1]);
    
    plhs[0] = mxCreateDoubleMatrix(inCateg[nrows-1],1,mxREAL);
    mM = mxGetPr(plhs[0]);

    MBCategMean(nrows,inM,inCateg,mM);
}