#include "mex.h"


void weightedSum(mwIndex* ir1, mwIndex* jc1, double* aff1, 
        mwIndex* ir2, mwIndex* jc2, double* aff2, 
        mwSize n, mwSize end, double w, 
        mwIndex* ir_out, mwIndex* jc_out, double* outSimAtts)
{

 
    int nnz = 0;
	int currCol, currRow;
    int i1, r1, startIdx1, stopIdx1;
    int i2, r2, startIdx2, stopIdx2;
    double val, val1, val2;
    int MAX = 2147483647;
	
	printf("starting weightedSum ...\n");
    
	for (currCol = 0; currCol < n; currCol++) {
   		jc_out[currCol] = nnz;
        
        startIdx1   = (int)jc1[currCol];
		stopIdx1    = (int)jc1[currCol + 1];
		i1          = startIdx1;
        
        startIdx2   = (int)jc2[currCol];
		stopIdx2    = (int)jc2[currCol + 1];
		i2          = startIdx2;
        
//         printf("weightedSum: currCol = %d: Idx1=[%d,%d], Idx2=[%d,%d] ...\n",
//                 currCol,startIdx1,stopIdx1,startIdx2,stopIdx2);
        
        while (i1 < stopIdx1 || i2 < stopIdx2) {
         
            if (i1 < stopIdx1) 
                r1 = (int)ir1[i1];
            else 
                r1 = MAX;
            
            if (i2 < stopIdx2)
                r2 = (int)ir2[i2];
            else
                r2 = MAX;
            
 //           printf("weightedSum:currCol = %d, r1=%d, r2=%d ...\n",currCol, r1,r2);
            if (r1 < r2) {
                currRow =  r1;
                val1 = aff1[i1];
                val2 = 0;
                i1++;
            } else if (r2 < r1) {
                currRow =  r2;
                val1 = 0;
                val2 = aff2[i2];
                i2++;
            } else { // r1==r2
                currRow =  r1;
                val1 = aff1[i1];
                val2 = aff2[i2];
                i1++; i2++;
            }

 //           printf("weightedSum:currCol = %d, currRow = %d ...\n",currCol, currRow);
            if (currCol < end && currRow < end) {
                val = val1*(1-w)+val2*w;
            } else {
                val = val1;
            }
            if (val > 0) {
                outSimAtts[nnz] = val;
                ir_out[nnz] = currRow;
                nnz++;
                if (nnz % 1000000 == 0)
                    printf("weightedSum nnz=%d M\n", nnz/1000000);
           }
		}
	} 
    
	jc_out[n] = nnz;

    printf("end weightedSum => nnz=%d\n", nnz);
    
}


/* The gateway function 
//this function accepts two parameters:
//1. a sparse structure affinity matrix
//2. a sparse structure affinity matrix
//3. weight
//4. end index
*/
void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[])
{
    /* variable declarations here */
    mwIndex *ir1, *jc1;
    mwIndex *ir2, *jc2;
    double  *aff1, *aff2;
    mwSize  n,m,n2,m2;
    double  w;
    int     end;
    int     nnz;
    

    mwIndex*    ir_out;
    mwIndex*    jc_out;
    double*     outSimAtts;
    mwSize      estimateNNZ;
    double      percent_sparse = 1; //0.75;
    
    
    printf("weightedSum: nlhs = %d, nrhs = %d\n", nlhs, nrhs);
    
    /* Check for the correct number of outputs */
    if(nlhs != 1){
        mexErrMsgIdAndTxt( "MATLAB:weightedSum:wrongrhs",
                "Wrong number of output arguments.");
    }
    /* Check for the correct number of inputs */
    if(nrhs != 4){
        mexErrMsgIdAndTxt( "MATLAB:weightedSum:wrongrhs",
                "Wrong number of input arguments.");
    }
    
    aff1 = mxGetPr(prhs[0]);  /* pointer to first input  */
    ir1 = mxGetIr(prhs[0]);
    jc1 = mxGetJc(prhs[0]);
    
    aff2 = mxGetPr(prhs[1]);  /* pointer to second input  */
    ir2 = mxGetIr(prhs[1]);
    jc2 = mxGetJc(prhs[1]);
    
    w   = mxGetScalar(prhs[2]);         /* pointer to thrid input */
    end = (int)mxGetScalar(prhs[3]);    /* pointer to fourth input */

    printf("weightedSum: w=%f, end=%d\n", w, end);
    
    /* dimensions of input matrices */
    m   = mxGetM(prhs[0]);
    n   = mxGetN(prhs[0]);
    m2  = mxGetM(prhs[1]);
    n2  = mxGetN(prhs[1]);
    
    printf("weightedSum dimensions: m=%d, n=%d, m2=%d, n2=%d\n", m, n, m2, n2);
    if (m != m2 || n != n2) {
        mexErrMsgIdAndTxt("MATLAB:attributesSimilarity:matchdims",
            "Inner dimensions of matrices do not match.");
    }
    if (n != m) {
        mexErrMsgIdAndTxt("MATLAB:attributesSimilarity:square",
            "Function requires input matrix 1 must be square.");
    }
    
    nnz = jc1[n]+jc2[n];
    printf("weightedSum: nnz=%d, percent_sparse=%f, mul=%f\n", 
            nnz, percent_sparse, percent_sparse*nnz);
    estimateNNZ = (mwSize)(percent_sparse*nnz);
    printf("weightedSum: jc1[n]=%d, jc2[n]=%d, estimateNNZ=%d\n", 
            jc1[n], jc2[n], estimateNNZ);
    

    // prepare output: nxn
    // allocate memeory according to estimateNNZ
    plhs[0] = mxCreateSparse(n,n, estimateNNZ, mxREAL);
    outSimAtts = mxGetPr(plhs[0]);
    ir_out = mxGetIr(plhs[0]);
    jc_out = mxGetJc(plhs[0]);

    // calculate the weighted sum
    // affinity(1:e, 1:e)= affinity(1:e,1:e)*(1-attWeight)+attAffinity(1:e,1:e)*attWeight;
    weightedSum(ir1, jc1, aff1, ir2, jc2, aff2, n, end, w, ir_out, jc_out, outSimAtts);

}