#include "mex.h"


void attributesSimilarity5(mxLogical* nodeAtts, 
        mwSize numNodes, mwSize numAttrs, double threshold, 
        double* idxNodesToKeep, mwSize numNodesToKeep,
        mwIndex* ir_out, mwIndex* jc_out, double* outSimAtts)
{

	double  total=0, common=0, affinity=0;
    int     nnz = 0;
	int     prevCol = -1;
    int     currCol, currRow, i, j, k;
    int     a;
    bool    a1, a2;
    
	
// 	printf("attributesSimilarity5: nnz = %d\n", nnz_in);
	printf("start attributesSimilarity5 ... \n");
    
    for (j=0; j<numNodesToKeep; j++){
        currCol = (int)idxNodesToKeep[j]-1;
        for (k=prevCol+1; k<=currCol; k++){
            jc_out[k] = nnz;
        }

        for (i=0; i<numNodesToKeep; i++) {
            currRow = (int)idxNodesToKeep[i]-1;
            
            affinity = 0;
            
           if (currCol == currRow){
                affinity = 1;
           } else {
                total    = 0;
                common   = 0;
                for (a=0; a<numAttrs; a++){
                    a1 = (bool)nodeAtts[currCol*numAttrs+a];
                    a2 = (bool)nodeAtts[currRow*numAttrs+a];
                    if (a1 || a2) {total++;}
                    if (a1 && a2) {common++;}
                }
                    
                if (total > 0) affinity=common/total;
            }
            
            if (affinity > threshold) {
//                 printf("attributesSimilarity5 [j=%d,i=%d][col=%d,row=%d]\n",j,i,currCol,currRow);
                outSimAtts[nnz] = affinity;
                ir_out[nnz] = currRow;
                nnz++;
                if (nnz % 10000000 == 0)
                    printf("attributesSimilarity5 nnz=%d M\n", nnz/1000000);
             }
        }
        prevCol = currCol;
    }

    for (k=prevCol+1; k<=numNodes; k++){
        jc_out[k] = nnz;
    }
    printf("end attributesSimilarity5 => nnz=%d\n", nnz);
    
//     mxDestroyArray(va1_P);
//     mxDestroyArray(va2_P);
//     mxDestroyArray(vaC_P);
//     mxDestroyArray(vaC_P);
    
    
}


/* The gateway function 
//this function accepts four parameters:
//1. an attributes matrix 
//2. an array with idxNodesToKeep
//3. min threshold
*/
void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[])
{
    /* variable declarations here */
    mwSize  numNodes, numAttrs;
    mxLogical* nodeAtts;
    double* idxNodesToKeep;
    double  threshold;
    mwSize  numNodesToKeep, estimateNNZ;
    double  percentSparse = 0.25; //0.33;

    
    // ouput
    mwIndex*    ir_out;
    mwIndex*    jc_out;
    double*     outSimAtts;
    
    printf("nlhs = %d, nrhs = %d\n", nlhs, nrhs);
    
    /* Check for the correct number of outputs */
    if(nlhs != 1){
        mexErrMsgIdAndTxt( "MATLAB:attributesSimilarity:wrongrhs",
                "Wrong number of output arguments.");
    }
    
    /* Check for the correct number of inputs */
    if(nrhs != 4){
        mexErrMsgIdAndTxt( "MATLAB:attributesSimilarity:wrongrhs",
                "Wrong number of input arguments.");
    }
    
    nodeAtts = mxGetLogicals(prhs[0]);    /* pointer to first input  */
    numAttrs = mxGetM(prhs[0]);     // rows
    numNodes = mxGetN(prhs[0]);     // cols
    
    idxNodesToKeep = mxGetPr(prhs[1]);  /* pointer to second input matrix */
    numNodesToKeep = mxGetM(prhs[1]);
    
    threshold = mxGetScalar(prhs[2]);  /* pointer to thrid input matrix */
    percentSparse = mxGetScalar(prhs[3]);  /* pointer to fourth input matrix */
    
    // *** estimateNNZ as part of population 
    estimateNNZ =  numNodesToKeep*numNodesToKeep*percentSparse; 

    printf("nodeAtts dimensions: numNodes=%d, numAttrs=%d\n", numNodes, numAttrs);
    printf("params: threshold = %f, percentSparse = %f, numNodesToKeep= %d, estNNZ = %d\n", threshold, percentSparse, numNodesToKeep, estimateNNZ);
    

    // prepare output: nxn
    // allocate memeory according to estimateNNZ
    plhs[0] = mxCreateSparse(numNodes, numNodes, estimateNNZ, mxREAL); 
    outSimAtts = mxGetPr(plhs[0]);
    ir_out = mxGetIr(plhs[0]);
    jc_out = mxGetJc(plhs[0]);
    
    // calculate the affinity
    attributesSimilarity5(nodeAtts, numNodes, numAttrs, threshold, idxNodesToKeep, numNodesToKeep, ir_out, jc_out, outSimAtts);


}