MATLAB_Codes_for_Dimensionality_Reduction.rar_dimention reductio

/* * * mexCCACollectdata.c * prepare data for cca.m to form SDP * * by feisha@cis.upenn.edu */ #include "mex.h" #include "matrix.h" #include <stdlib.h> #include <float.h> #include <string.h> #include <math.h> /* the computation engine */ void collectdata(double *x, double *y, int* edgerow, int *edgecol, int *relative, int *nv, int *vidx, int D, int d, int n, int ks, double *a, double *b, double *g); /* auxiliary functions */ void sanity_check(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]); /* the gateway */ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { int n, ks=1, D, d; double *x, *y, *a, *b, *g; int *neighbors; /* sanity_check(nlhs, plhs, nrhs, prhs); */ n = mxGetN(prhs[0]); D = mxGetM(prhs[0]); d = mxGetM(prhs[1]); /*printf("%d data points, reducing from dimension %d ---> dimension %d, using \ %d nearest neighbors\n", n, D, d, ks); */ plhs[0] = mxCreateDoubleMatrix(d*d, d*d, mxREAL); plhs[1] = mxCreateDoubleMatrix(d*d, n, mxREAL); plhs[2] = mxCreateDoubleMatrix(n,1, mxREAL); collectdata( (double*)mxGetData(prhs[0]), (double*)mxGetData(prhs[1]), (int*)mxGetData(prhs[2]), (int*)mxGetData(prhs[3]), (int*)mxGetData(prhs[4]), (int*)mxGetData(prhs[5]), (int*)mxGetData(prhs[6]), D, d, n, ks, (double*)mxGetData(plhs[0]), (double*)mxGetData(plhs[1]), (double*)mxGetData(plhs[2])); return; } /* check input arguments */ void sanity_check(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { int n1, n2, n3; /* Check for proper number of input and output arguments. */ if (nrhs != 4) { mexErrMsgTxt("Four input argument required."); } if (nlhs != 3) { mexErrMsgTxt("Three output arguments required."); } /* Check data type of input argument. */ /*if (!(mxIsDouble(prhs[0])) || !(mxIsDouble(prhs[1])) || !(mxIsDouble(prhs[2]))) { mexErrMsgTxt("Input array must be of type double."); }*/ /* Check dimensions */ if(mxGetNumberOfDimensions(prhs[0])!=2 || mxGetNumberOfDimensions(prhs[1])!=2 || mxGetNumberOfDimensions(prhs[2])!=2) { mexErrMsgTxt("Input arrays must be two-dimension arrays"); } /* matching dimensions */ n1 = mxGetN(prhs[0]); n2 = mxGetN(prhs[1]); n3 = mxGetN(prhs[2]); if (n1 != n2 || n1 !=n3) { mexErrMsgTxt("Input arrays must have same number of columns!"); } return; } /* computing */ double vecdot(double* v1, int n) { /* v1'*v2 */ double x = 0; int i; for(i = 0; i < n; i++) { x += v1[i]*v1[i]; } return x; } void vecsub(double* v1,double* v2, double* v3, int n) { /* v3=v1-v2 */ int incx=1; double none = -1; /* dcopy_(&n, v1, &incx,v3, &incx); daxpy_(&n, &none, v2, &incx, v3, &incx);*/ int i; for(i=0;i<n;i++) v3[i]=v1[i]-v2[i]; } void vecout(double* v1, double* v3, int d) { /* v3 = v1*v2' */ double temp,one = 1.0; int inca = 1; int i,j,c; /* dger_(&d , &d, &one, v1, &inca,v1, &inca, v3, &d); */ c=0; for(i=0;i<d;i++){ temp=v1[i]; for(j=0;j<d;j++){ v3[c]=temp*v1[j]; c++; } } } void makevecsym(double *v, int d) { /* assume V is a dxd matrix, return (V+V')/2 */ int row, col, i; double temp; for(col=0; col < d; col++) { for(row=col; row < d; row++) { temp = (v[row*d+col] + v[col*d+row])/2; v[row*d+col] = v[col*d+row] = temp; } } } void updateA(double *a, double *v, double alpha, int d) { /* update A */ int inca = 1; char *uplo= "U"; double temp; int i,j,c; /* dspr_(uplo, &d , &alpha, v, &inca, a);*/ if(alpha==1.0){ /*optimize the common case */ c=0; for(i=0;i<d;i++){ temp=v[i]; for(j=0;j<=i;j++) { a[c]+=temp*v[j]; c++; } } } else{ c=0; for(i=0;i<d;i++){ temp=v[i]*alpha; for(j=0;j<=i;j++) { a[c]+=temp*v[j]; c++; } } } } void updateB(double * b, double alpha, double *v, int d) { /* update B(:,i) */ int i; int incx = 1; /* daxpy_(&d, &alpha, v, &incx, b, &incx);*/ for(i=0;i<d;i++) b[i]+=alpha*v[i]; } /* transform upper triangular matrix stored in A as full matrix */ void recoverA(int D, double *A, double *tempA) { int l = 0, idx, jdx,s, incx=1, incy=1; int n = D*(D+1)/2; for(jdx=0; jdx <D; jdx++) { for(idx=0; idx <=jdx; idx++) { A[jdx*D+idx] = tempA[l]; l++; } } for(jdx=0; jdx <D; jdx++) for(idx=jdx; idx <D; idx++) { A[jdx*D+idx] = A[idx*D+jdx]; } } void collectdata(double *x, double *y, int* edgesrow, int *edgescol, int *relative, int *nv, int *vidx, int D, int d, int n, int ks, double *a, double *b, double *g) { double *diffx, *diffy, *yij, *nid; double gij; int i, j, k, iEdge, iVertex,itemp, ii; int nn_start, nn_end, nn; double *diffxjk, *diffyjk; double gjk, *yjk, *tempA; int relflg = *relative; /* clear data area */ memset(a, 0, sizeof(double)*d*d*d*d); memset(b, 0, sizeof(double)*d*d*n); memset(g, 0, sizeof(double)*n); /* temporary working space */ diffx = (double*)calloc(D, sizeof(double)); diffy = (double*)calloc(d, sizeof(double)); yij = (double*)calloc(d*d, sizeof(double)); tempA = (double*)calloc((d*d)*(d*d+1)/2, sizeof(double)); if(!diffx || !diffy || !yij || !tempA) { mexErrMsgTxt("Out of memory..cannot allocate working space.."); return; } iEdge = 0; iVertex = 0; if(relflg==0) { for(i=0; i <n ; i++) { /* figure out the nearest neighbors */ nn_start = edgescol[i]; nn_end = edgescol[i+1]-1; if (nn_end >= nn_start) { for(nn=nn_start; nn<=nn_end; nn++) { j = edgesrow[nn]; /* compute diffx and diffy */ vecsub(x+i*D, x+j*D, diffx, D); gij = vecdot(diffx, D); vecsub(y+i*d, y+j*d, diffy, d); vecout(diffy, yij, d); /* update A, b */ updateA(tempA, yij, 1.0*nv[iEdge],d*d); for(itemp=0; itemp < nv[iEdge]; itemp++) { ii = vidx[iVertex]-1; updateB(b+ii*d*d, gij, yij, d*d); g[ii] = g[ii]+gij*gij; iVertex++; } iEdge++; } } } printf("Iedge :%d, ivertex: %d\n", iEdge, iVertex); } else { /* reproducing the code above to some degree so that we don't have to do if statement inside cache-intensive loop */ for(i=0; i <n ; i++) { /* figure out the nearest neighbors */ nn_start = edgescol[i]; nn_end = edgescol[i+1]-1; if (nn_end >= nn_start) { for(nn=nn_start; nn<=nn_end; nn++) { j = edgesrow[nn]; /* compute diffx and diffy */ vecsub(x+i*D, x+j*D, diffx, D); gij = vecdot(diffx, D); vecsub(y+i*d, y+j*d, diffy, d); vecout(diffy, yij, d); /* update A, b */ updateA(tempA, yij, 1*nv[iEdge]/(gij*gij), d*d); for(itemp=0; itemp < nv[iEdge]; itemp++) { ii = vidx[iVertex]-1; updateB(b+ii*d*d, 1/gij, yij, d*d); g[ii] = g[ii]+1; iVertex++;