#include "mex.h"


void attributesSimilarity5a(double* 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, a1, a2, a3;
	
// 	printf("attributesSimilarity5a: nnz = %d\n", nnz_in);
	printf("start attributesSimilarity5a ... \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 {
                affinity = 0;
                total    = 0;
                common   = 0;
                for (a=0; a<numAttrs; a++){
                    a1 = (int)nodeAtts[a*numNodes+currCol];
                    a2 = (int)nodeAtts[a*numNodes+currRow];
                    a3 = a1+a2;
                    if (a3 == 1) {total++;}
                    if (a3 == 2) {total++; common++;}
//                     if (a1==1 || a2==1) {
//                        total++;
//                        if (a1==a2) common++;
//                     }
                }
                if (total > 0) affinity=common/total;
            }
            
            if (affinity > threshold) {
//                 printf("attributesSimilarity5a [j=%d,i=%d][col=%d,row=%d]\n",j,i,currCol,currRow);
                outSimAtts[nnz] = affinity;
                ir_out[nnz] = currRow;
                nnz++;
                if (nnz % 1000 == 0)
                    printf("attributesSimilarity5a nnz=%d K\n", nnz/1000);
             }
        }
        prevCol = currCol;
    }
    
    for (k=prevCol+1; k<=numNodes; k++){
        jc_out[k] = nnz;
    }
    printf("end attributesSimilarity5a => nnz=%d\n", nnz);
    
}


/* 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;
    double* nodeAtts;
    double* idxNodesToKeep;
    double  threshold;
    mwSize  numNodesToKeep, estimateNNZ;
    double  percent_sparse = 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 != 3){
        mexErrMsgIdAndTxt( "MATLAB:attributesSimilarity:wrongrhs",
                "Wrong number of input arguments.");
    }
    
    nodeAtts = mxGetPr(prhs[0]);    /* pointer to first input  */
    numNodes = mxGetM(prhs[0]);     // rows
    numAttrs = mxGetN(prhs[0]);     // cols
    
    idxNodesToKeep = mxGetPr(prhs[1]);  /* pointer to second input matrix */
    numNodesToKeep = mxGetM(prhs[1]);
    // *** estimateNNZ as third of population 
    estimateNNZ =  numNodesToKeep*numNodesToKeep*percent_sparse; 
    
    threshold = mxGetScalar(prhs[2]);  /* pointer to thrid input matrix */
    
    printf("nodeAtts dimensions: numNodes=%d, numAttrs=%d\n", numNodes, numAttrs);
    printf("params: threshold = %f, numNodesToKeep= %d, estNNZ = %d\n", threshold, 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
    attributesSimilarity5a(nodeAtts, numNodes, numAttrs, threshold, idxNodesToKeep, numNodesToKeep, ir_out, jc_out, outSimAtts);


}