bvr_delta_new.c

best routing protocol

/*
 * [delta] = bvr_delta(beacon_vectors, C_k)
 * 
 *
 * Inputs:
 * Outputs:

 */

#include <mex.h>
#include <string.h>

void bvr_delta(double * delta,
			   double* bv_src, double* bv_dest, double *C_k,
			   int n_src, int n_dest, int k, int r, int op)
{
		int i,j, index, src, dest, offset, idx;
		double delta_plus;
		double delta_minus, diff;
		double *C_k_dest;
		
		/*mexPrintf("n=%d,r=%d,k=%d\n",n,r,k);
		mexPrintf("beacon(1,1)%f,(2,1)%f\n",beacon_vectors[0],beacon_vectors[1]);
		*/
		if (op==1) /* default */{
				for (dest=0; dest<n_dest ; dest++){
						C_k_dest = C_k + dest;
						
						for (src=0,offset=dest*n_dest; src<n_src; src++){
								delta_plus =0.0;
								delta_minus = 0.0;
								for (i=0; i<k; i++){
										idx = (int)(C_k_dest[i*n_dest]-1);
										
										diff = bv_src[src+idx*n_src]-bv_dest[dest+idx*n_dest];
										if (diff>0.)
												delta_plus += diff;
										else
												delta_minus -= diff;
										delta[offset+i*n_dest*n_src]=10*delta_plus+delta_minus;
										/*
										  mexPrintf("src,%d,dest %d, i%d,%f(offset %d)\n",src,dest,i, delta[offset+i*n*n],offset+i*n*n);
										*/
								}
								offset+=1;
						}
				}
		}
		else{/* sqr L2 norm*/
				for (dest=0; dest<n_dest ; dest++){
						C_k_dest = C_k + dest;
						for (src=0,offset=dest*n_dest; src<n_src; src++, offset++){
								delta_plus =0.0;
								for (i=0; i<k; i++){
										idx = (int)(C_k_dest[i*n_dest]-1);
										
										diff = bv_src[src+idx*n_src]-bv_dest[dest+idx*n_dest];
										delta_plus += diff*diff;
										
										/*
										  mexPrintf("src,%d,dest %d, i%d,%f(offset %d)\n",src,dest,i, delta[offset+i*n*n],offset+i*n*n);
										*/
								}
								delta[offset] = delta_plus;
						}
				}
				
		}
		
		return;
}


void mexFunction(int nlhs, mxArray** plhs,
                 int nrhs, const mxArray** prhs)
{
    int n_src, n_dest, n_dest1, r1, k, r,op;
	int dims[3];
	
    if (nrhs != 4)
        mxErrMsgTxt("Invalid number of arguments");
	
	/*bv_src, n_src x r*/
	n_src = mxGetM(prhs[0]);
    r = mxGetN(prhs[0]);
	/*bv_dest, n_dest x r*/
	n_dest = mxGetM(prhs[1]);
    r1 = mxGetN(prhs[1]);
	if (r1!=r)
			mxErrMsgTxt("Beacons numbers don't match");
	/*C_k,  n_dest x k*/
	n_dest1 = mxGetM(prhs[2]);
	k = mxGetN(prhs[2]);
	if (n_dest!=n_dest1)
			mxErrMsgTxt("C_k and bv_dest don't match");
	/*op: how to calculate distance function 0:standard 1: sqr L2 norm*/
    if (!mxIsNumeric(prhs[3]) || mxGetM(prhs[3]) * mxGetN(prhs[3]) != 1)
        mxErrMsgTxt("Op must be a scalar");
    op = (int) (*mxGetPr(prhs[3]));
	if ( (op<0) || (op>1) )
			mxErrMsgTxt("invalid Op: only 0 or 1 supported");
	
	/*output: n_src x n_dest x k*/
	dims[0]=n_src;
	dims[1]=n_dest;
	if (op==1)
			dims[2]=k;
	else
			dims[2]=1;
	
	plhs[0] = mxCreateNumericArray(3, dims, mxDOUBLE_CLASS, mxREAL);

	bvr_delta(mxGetPr(plhs[0]),
			  mxGetPr(prhs[0]), mxGetPr(prhs[1]), mxGetPr(prhs[2]), n_src,n_dest, k,r, op);
	
    /*dijkstra(mxGetPr(plhs[0]), mxGetPr(plhs[1]), mxGetPr(plhs[2]),
             mxGetJc(prhs[0]), mxGetIr(prhs[0]), mxGetPr(prhs[0]),
             s, n);*/
}