sparse_rank_one_mex.c

关于有直接稀疏PCA的方法

// MEX wrapper for the sparse_rank_one function in sparsesvd.c
// function [X,U,u,F,k] = sparse_rank_one_mex(A,rho,tol,MaxIter,Info,[F],[X0],[k0]) 

// Last Modified: A. d'Aspremont, Laurent El Ghaoui, Ronny Luss, November 2006.
// http://www.carva.org/alexandre.daspremont


#include "sparsesvd.h" 

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
	double *Amat, *Bmat, *Rmat, *Xmatout, *Fmatout, *Umat, *uvec, *iter;
	// The following are the extra parameters that are returned to Matlab for analyzing results
	double *dualitygap_alliter, *cputime_alliter, *perceigs_alliter;
	double rho,tol,MaxIter,info,algo,numeigs,addeigs,checkgap,perceigs,check_for_more_eigs;
	int n,m,p;
	
	// Inputs:
	Amat = mxGetPr(prhs[0]);
	n = mxGetN(prhs[0]);
	Bmat = mxGetPr(prhs[1]);
	Rmat= mxGetPr(prhs[2]);
	m = mxGetM(prhs[2]);
	rho = mxGetScalar(prhs[3]);
	tol = mxGetScalar(prhs[4]);
	MaxIter = mxGetScalar(prhs[5]);
	info = mxGetScalar(prhs[6]);
	algo = mxGetScalar(prhs[7]); // 1 for full eig decomp, 3 for partial eig decomp, 5 for Gen. Eig. w/ full eig.
	numeigs = mxGetScalar(prhs[8]);
	addeigs = mxGetScalar(prhs[9]);
	checkgap = mxGetScalar(prhs[10]);
	perceigs = mxGetScalar(prhs[11]);
	check_for_more_eigs = mxGetScalar(prhs[12]);
	// Outputs:
	if (algo==5) {p=m;} else {p=n;}
	plhs[0] = mxCreateDoubleMatrix(p, p, mxREAL);
	Xmatout = mxGetPr(plhs[0]);
	plhs[1] = mxCreateDoubleMatrix(p, p, mxREAL);
	Umat = mxGetPr(plhs[1]);
	plhs[2] = mxCreateDoubleMatrix(p, 1, mxREAL);
	uvec = mxGetPr(plhs[2]);
	plhs[3] = mxCreateDoubleMatrix(p, p, mxREAL);
	Fmatout = mxGetPr(plhs[3]);
	plhs[4] = mxCreateDoubleMatrix(1, 1, mxREAL);
	iter=mxGetPr(plhs[4]);
	plhs[5] = mxCreateDoubleMatrix((int)MaxIter, 1, mxREAL);
	dualitygap_alliter=mxGetPr(plhs[5]);
	plhs[6] = mxCreateDoubleMatrix((int)MaxIter, 1, mxREAL);
	cputime_alliter=mxGetPr(plhs[6]);
	plhs[7] = mxCreateDoubleMatrix((int)MaxIter, 1, mxREAL);
	perceigs_alliter=mxGetPr(plhs[7]);
	
	switch ((int)algo){
		case 1: {sparse_rank_one(Amat,n,rho,tol,(int)(MaxIter),Xmatout,Umat,uvec,Fmatout,iter,(int)(info),(int)(checkgap),dualitygap_alliter,cputime_alliter);break;}
		case 3: {sparse_rank_one_partialeig(Amat,n,rho,tol,(int)(MaxIter),Xmatout,Umat,uvec,Fmatout,iter,(int)(info),(int)(numeigs),(int)(addeigs),(int)(checkgap),perceigs,(int)(check_for_more_eigs),dualitygap_alliter,cputime_alliter,perceigs_alliter);break;}
		case 5:	{sparse_geneig(Amat,Bmat,Rmat,n,m,rho,tol,(int)(MaxIter),Xmatout,Umat,uvec,Fmatout,iter,(int)(info),(int)(checkgap),dualitygap_alliter,cputime_alliter);break;}	  
		case 9: {sparse_rank_one_partialeig_matlab(Amat,n,rho,tol,(int)(MaxIter),Xmatout,Umat,uvec,Fmatout,iter,(int)(info),(int)(numeigs),(int)(addeigs),(int)(checkgap),perceigs,(int)(check_for_more_eigs),dualitygap_alliter,cputime_alliter,perceigs_alliter);break;}
		default:{mexPrintf("DSPCA: algo option not supported \n");mexEvalString("drawnow;");}
	} 
}