mexccacollectdata_windows.c

Fei Sha 等人编写的流形学习算法CCA的matlab代码

/*
 *
 * mexCCACollectdata.c
 *      prepare data for cca.m to form SDP
 *
 * by [e-mail]
 */
 
 #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);
    
}
void vecout(double* v1, double* v3, int d) { /* v3 = v1*v2' */
    double one = 1.0;
    int inca = 1;
    memset(v3, 0, sizeof(double)*d*d);
    dger(&d , &d, &one, v1, &inca,v1, &inca, v3, &d); 
}
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";
    dspr(uplo, &d , &alpha, v, &inca, a);
}
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);
}

/* 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++;
                    }
                    iEdge++;
                }   
            }
        }
    }
    /* make A symmetric */
    recoverA(d*d, a, tempA); 
}