pcaL1_1.2.1.tar.gz_L1-2_PCA L1_PCA-L1_PCA_L1_R1-PCA

#include "Clp_C_Interface.h" #include <stdlib.h> #include <math.h> #include "type.h" int solveL1PCAStar (ENTITYINFOptr entityinfo, SOLVERINFOptr solverinfo, PROBLEMINFOptr probleminfo); int initialSVD (ENTITYINFOptr entityinfo, PROBLEMINFOptr probleminfo); int setupCLP (SOLVERINFOptr solverinfo); int loadClpProblem (ENTITYINFOptr entityinfo, SOLVERINFOptr solverinfo, PROBLEMINFOptr probleminfo); int changeBounds (SOLVERINFOptr solverinfo, PROBLEMINFOptr probleminfo, int l); int optimize (SOLVERINFOptr solverinfo, PROBLEMINFOptr probleminfo); int getBeta (SOLVERINFOptr solverinfo, PROBLEMINFOptr probleminfo); int getProjectedPoints (ENTITYINFOptr entityinfo, PROBLEMINFOptr probleminfo); int dgesvd (char jobu, char jobvt, int m, int n, double *A, int lda, double *S, double *VT, int ldu, double *Umat, int ldvt, double *work, int lwork); /* SVD */ void dgemm (char transa, char transb, int m, int n, int k, double alpha, double *A, int lda, double *B, int ldb, double beta, double *C, int ldc); /* multiply A *B = C */ void dgemv (char trans, int m, int n, double alpha, double *A, int lda, double *x, int incx, double beta, double *y, int incy); /* multiply Ax = y */ int solveL1PCAStar (ENTITYINFOptr entityinfo, SOLVERINFOptr solverinfo, PROBLEMINFOptr probleminfo) { int status = probleminfo->status; int l = probleminfo->l; double minobjective = probleminfo->minobjective; int numprojdim = probleminfo->numprojdim; int numfactors = probleminfo->numfactors; int numentities_n = entityinfo->numentities_n; int numattributes_m = entityinfo->numattributes_m; numprojdim = 0; status = setupCLP (solverinfo); /*initialize CLP */ if (status) { REprintf ("Error setting up CLP\n"); return 1; } /* if numattributes_m > numentities_n, project into numentities_n space, then set numattributes_m = numentities_n and continue */ if (numattributes_m > numentities_n) { status = initialSVD (entityinfo, probleminfo); if (status) { REprintf ("Error with initial SVD\n"); return 1; } } for (probleminfo->projdim = numfactors - 1; probleminfo->projdim > 0; --probleminfo->projdim) { /* the main loop */ minobjective = (OBJ_INIT); /* objective of the best l1 regression */ Clp_setDualObjectiveLimit(solverinfo->model, minobjective); status = loadClpProblem(entityinfo, solverinfo, probleminfo); /*set up columns/variables */ if (status) { REprintf ("Error with initial load of problem\n"); return 1; } REprintf ("%d ", probleminfo->projdim); if ((VERBOSITY) >= 1) { REprintf ("projdim %d\n", probleminfo->projdim); } for (l = 0; l <= probleminfo->projdim; ++l) { /* solve l1 regression for each direction, lth attribute is response */ if ((VERBOSITY) >= 1) { REprintf ("l %d minobjective %f\n", l, minobjective); } status = changeBounds (solverinfo, probleminfo, l); /* set beta_l = -1, others unbounded */ if (status) { REprintf ("Error changing beta bounds\n"); return 1; } status = optimize (solverinfo, probleminfo); /* solve with Clp, get objective value */ if (status) { REprintf ("Error solving l1 regression subproblem\n"); return 1; } if (minobjective > probleminfo->objective) { /* determine if objective is best so far */ probleminfo->bestdir[probleminfo->projdim] = l; minobjective = probleminfo->objective; probleminfo->minobjective = minobjective; if ((VERBOSITY) >= 4) { REprintf ("new bestdir %d new best objective %f\n", l, probleminfo->objective); } status = getBeta(solverinfo, probleminfo); /* store coefficients of regression */ if (status) { REprintf ("Error storing beta's\n"); return 1; } Clp_setDualObjectiveLimit(solverinfo->model, minobjective); } if (minobjective < (EPSILON)) { /* if error is 0, go to next dimension */ l = probleminfo->projdim + 1; } } /* end loop on l1 regressions */ if ((VERBOSITY) >= 1) { REprintf ("avg error %f\n", minobjective/((double) entityinfo->numentities_n)); } if (minobjective/((double) entityinfo->numentities_n) < 1.0) { ++numprojdim; } if ((VERBOSITY) >= 4) { REprintf ("bestdir %d best objective %f\n", probleminfo->bestdir[probleminfo->projdim], minobjective); } status = getProjectedPoints (entityinfo, probleminfo);/* get projected points */ if (status) { REprintf ("Error getting projected points, or done\n"); return 1; } } return 0; } /*end solveL1PCAStar */ int setupCLP (SOLVERINFOptr solverinfo) { solverinfo->model=Clp_newModel(); Clp_setLogLevel(solverinfo->model, 0); return 0; } /* end setupCLP */ int initialSVD (ENTITYINFOptr entityinfo, PROBLEMINFOptr probleminfo) { int i = probleminfo->i; int j = probleminfo->j; int status = probleminfo->status; double *work = probleminfo->work; int lwork = probleminfo->lwork; double *S = probleminfo->S; double *VT = probleminfo->VT; int numattributes_m = entityinfo->numattributes_m; int numentities_n = entityinfo->numentities_n; double *points_XT = entityinfo->points_XT; if ((VERBOSITY) >= 7) { REprintf ("points\n"); for (i = 0; i < numattributes_m; ++i) { for (j = 0; j < numentities_n; ++j) { REprintf ("%f ", points_XT[numattributes_m * j + i]); } REprintf ("\n"); } } status = dgesvd ( 'A', 'A', numattributes_m, numentities_n, points_XT, numattributes_m, S, probleminfo->preVj, numattributes_m, VT, numentities_n, work, lwork); /* get A_q = preVj, points gets destroyed here; The columns of preVj define the new subspace */ if (status) { /* preV^jT is the first projdim columns of U, derived by doing an SVD on points_XT (which has dimension numattributes_m X numentities_n). It's like doing PCA on the transpose of the data matrix; U*S has the scores */ REprintf ("Error getting SVD, error %d\n", status); return 1; } for (i = 0; i < numentities_n; ++i) { for (j = 0; j < numentities_n; ++j) { entityinfo->points_XT[numentities_n * j + i] = probleminfo->VT[numentities_n * j + i] * S[i]; } } return 0; } /* end initial SVD */ int loadClpProblem (ENTITYINFOptr entityinfo, SOLVERINFOptr solverinfo, PROBLEMINFOptr probleminfo) { int numcols = probleminfo->numcols; int i = probleminfo->i; int j = probleminfo->j; double *obj = probleminfo->obj; char **colname = probleminfo->colname; int rcnt = probleminfo->rcnt; int nzcnt = probleminfo->nzcnt; double *rhs = probleminfo->rhs; int *matind = probleminfo->matind; double *matval = probleminfo->matval; int *matbeg = probleminfo->matbeg; int projdim = probleminfo->projdim; double *points_XT = entityinfo->points_XT; int numentities_n = entityinfo->numentities_n; rcnt = numentities_n; for (i = 0; i < numentities_n; ++i) { rhs[i] = 0.0; } nzcnt = 0; numcols = 0; for (j = 0; j < (projdim + 1); ++j) { matbeg[numcols] = nzcnt; probleminfo->betaind[j] = numcols; obj[numcols] = 0.0; probleminfo->lb[numcols]=-(DBL_MAX); probleminfo->ub[numcols]=DBL_MAX; sprintf (colname[numcols], "beta_%d", j); for (i = 0; i < numentities_n; ++i) { if (points_XT[i*(projdim + 1) + j]!=0.0) { matind[nzcnt]=i; matval[nzcnt]=points_XT[i * (projdim + 1) + j]; ++nzcnt; } } ++numcols; } for (i = 0; i < numentities_n; ++i) { matbeg[numcols] = nzcnt; probleminfo->eplusind[i] = numcols; obj[numcols] = 1.0; probleminfo->lb[numcols] = 0.0; probleminfo->ub[numcols]=(DBL_MAX); sprintf (colname[numcols], "eplus_%d", i); matind[nzcnt] = i; matval[nzcnt] = 1.0; ++nzcnt; ++numcols; } for (i = 0; i < numenti