/*=================================================================
* syntax: SPMX = SPARSIFY(MX, THRES)
*
* SPARSIFY  - sparsify the input matrix, i.e. ignore the values 
* 			of the matrix which	are below a threshold
* 			
* 			Input: 	- MX: m-by-n matrix (sparse or full)
* 					- THRES: threshold value (double)
*
* 			Output: - SPMX: m-by-n sparse matrix only with values 
* 					whose absolut value is above the given threshold
*
* 			Written by Mirko Visontai (10/24/2003)
*=================================================================*/


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

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    /* Declare variable */
    int i,m,n,nzmax,newnnz,col,processed,passed;
	int starting_row_index, current_row_index, stopping_row_index;
    double *in_pr,*in_pi,*out_pr,*out_pi;
    int *in_ir,*in_jc,*out_ir,*out_jc;
   	double thres;

    /* Check for proper number of input and output arguments */    
    if ((nlhs != 1) || (nrhs != 2)){
		mexErrMsgTxt("usage: SPMX = SPARSIFY(MX, THRES).");
    } 
	/* if matrix is complex threshold the norm of the numbers */
	if (mxIsComplex(prhs[0])){
		/* Check data type of input argument  */
		if (mxIsSparse(prhs[0])){

			/* read input */
			in_pr  	= mxGetPr(prhs[0]);
			in_pi  	= mxGetPi(prhs[0]);
			in_ir	= mxGetIr(prhs[0]);
			in_jc	= mxGetJc(prhs[0]);
			nzmax 	= mxGetNzmax(prhs[0]);
			m		= mxGetM(prhs[0]);
			n		= mxGetN(prhs[0]);
			thres 	= mxGetScalar(prhs[1]);

			/* Count new nonzeros */
			newnnz=0;
			for(i=0; i<nzmax; i++){
				if (sqrt(in_pr[i]*in_pr[i] + in_pi[i]*in_pi[i])>thres) {newnnz++;}
			}

			if (newnnz>0){
				/* create output */
				plhs[0] 	= mxCreateSparse(m,n,newnnz,mxCOMPLEX);
				if (plhs[0]==NULL)
					mexErrMsgTxt("Could not allocate enough memory!\n");
				out_pr 		= mxGetPr(plhs[0]);
				out_pi 		= mxGetPr(plhs[0]);
				out_ir 		= mxGetIr(plhs[0]);
				out_jc 		= mxGetJc(plhs[0]);
				passed		= 0;
				out_jc[0]	= 0;
				for (col=0; col<n; col++){
					starting_row_index = in_jc[col];
					stopping_row_index = in_jc[col+1];
					out_jc[col+1] = out_jc[col];
					if (starting_row_index == stopping_row_index)
						continue;
					else {
						for (current_row_index = starting_row_index;
							current_row_index < stopping_row_index;
							current_row_index++)  {
								if (sqrt(in_pr[current_row_index]*in_pr[current_row_index] + 
									     in_pi[current_row_index]*in_pi[current_row_index] ) > thres){

									out_pr[passed]=in_pr[current_row_index];	
									out_pi[passed]=in_pi[current_row_index];	
									out_ir[passed]=in_ir[current_row_index];	
									out_jc[col+1] = out_jc[col+1]+1;
									passed++;
								}
						}
					}
				}
			}
			else{
				plhs[0] = mxCreateSparse(m,n,0,mxCOMPLEX);
			}
		}
		else{ /* for full matrices */
			/* read input */
			in_pr  	= mxGetPr(prhs[0]);
			in_pi  	= mxGetPr(prhs[0]);
			m		= mxGetM(prhs[0]);
			n		= mxGetN(prhs[0]);
			thres 	= mxGetScalar(prhs[1]);

			/* Count new nonzeros */
			newnnz=0;
			for(i=0; i<m*n; i++){
				if (sqrt(in_pr[i]*in_pr[i] + in_pi[i]*in_pi[i])>thres) {newnnz++;}
			}

			if (newnnz>0){
				/* create output */
				plhs[0]	 	= mxCreateSparse(m,n,newnnz,mxCOMPLEX);
				if (plhs[0]==NULL)
					mexErrMsgTxt("Could not allocate enough memory!\n");
				out_pr 		= mxGetPr(plhs[0]);
				out_pi 		= mxGetPi(plhs[0]);
				out_ir 		= mxGetIr(plhs[0]);
				out_jc 		= mxGetJc(plhs[0]);
				passed		= 0;
				out_jc[0] 	= 0;

				for (col=0; col<n; col++){
					out_jc[col+1] = out_jc[col];
					for (current_row_index=0; current_row_index<m; current_row_index++){
							if (sqrt(in_pr[current_row_index+m*col]*in_pr[current_row_index+m*col] +
									 in_pi[current_row_index+m*col]*in_pi[current_row_index+m*col]) > thres){
								
								out_pr[passed]=in_pr[current_row_index+m*col];	
								out_ir[passed]=current_row_index;	
								out_jc[col+1] = out_jc[col+1]+1;
								passed++;
							}
					}
				}
			}
			else{
				plhs[0] = mxCreateSparse(m,n,0,mxCOMPLEX);
			}
		}
	}
	else { 
    	/* Check data type of input argument  */
    	if (mxIsSparse(prhs[0])){

			/* read input */
   			in_pr  	= mxGetPr(prhs[0]);
   			in_ir	= mxGetIr(prhs[0]);
    		in_jc	= mxGetJc(prhs[0]);
	   		nzmax 	= mxGetNzmax(prhs[0]);
			n		= mxGetN(prhs[0]);
			m		= mxGetM(prhs[0]);
			thres 	= mxGetScalar(prhs[1]);

	  		/* Count new nonzeros */
			newnnz=0;
			for(i=0; i<nzmax; i++){
				if ((fabs(in_pr[i]))>thres) {newnnz++;}
			}

			if (newnnz>0){
				/* create output */
	   	 		plhs[0] 	= mxCreateSparse(m,n,newnnz,mxREAL);
				if (plhs[0]==NULL)
					mexErrMsgTxt("Could not allocate enough memory!\n");
   	 			out_pr 		= mxGetPr(plhs[0]);
    			out_ir 		= mxGetIr(plhs[0]);
    			out_jc 		= mxGetJc(plhs[0]);
				passed		= 0;
				out_jc[0]	= 0;
				for (col=0; col<n; col++){
					starting_row_index = in_jc[col];
					stopping_row_index = in_jc[col+1];
					out_jc[col+1] = out_jc[col];
					if (starting_row_index == stopping_row_index)
						continue;
					else {
						for (current_row_index = starting_row_index;
							current_row_index < stopping_row_index;
							current_row_index++)  {
								if (fabs(in_pr[current_row_index])>thres){
									out_pr[passed]=in_pr[current_row_index];	
									out_ir[passed]=in_ir[current_row_index];	
									out_jc[col+1] = out_jc[col+1]+1;
									passed++;
								}
						}
					}
				}
			}
			else{
    			plhs[0] = mxCreateSparse(m,n,0,mxREAL);
			}
		}
		else{ /* for full matrices */
			/* read input */
   			in_pr  	= mxGetPr(prhs[0]);
			n		= mxGetN(prhs[0]);
			m		= mxGetM(prhs[0]);
			thres 	= mxGetScalar(prhs[1]);

	  		/* Count new nonzeros */
			newnnz=0;
			for(i=0; i<m*n; i++){
				if ((fabs(in_pr[i]))>thres) {newnnz++;}
			}

			if (newnnz>0){
				/* create output */
	   	 		plhs[0]	 	= mxCreateSparse(m,n,newnnz,mxREAL);
				if (plhs[0]==NULL)
					mexErrMsgTxt("Could not allocate enough memory!\n");
   	 			out_pr 		= mxGetPr(plhs[0]);
    			out_ir 		= mxGetIr(plhs[0]);
    			out_jc 		= mxGetJc(plhs[0]);
				passed		= 0;
    			out_jc[0] 	= 0;

				for (col=0; col<n; col++){
					out_jc[col+1] = out_jc[col];
					for (current_row_index=0; current_row_index<m; current_row_index++){
							if (fabs(in_pr[current_row_index+m*col])>thres){
								out_pr[passed]=in_pr[current_row_index+m*col];	
								out_ir[passed]=current_row_index;	
								out_jc[col+1] = out_jc[col+1]+1;
								passed++;
							}
					}
				}
			}
			else{
    			plhs[0] = mxCreateSparse(m,n,0,mxREAL);
			}
		}
	}
}