cddmex.c

一个matlab的将军模型

		a = mxGetPr(tmpa);
	  	for (i = 0 ; i < (int)(M->rowsize); i++) {
	  		b[i] = dd_get_d(M->matrix[i][0]);
    			for (j = 0; j < (int)(M->colsize) - 1; j++) {
      				a[i + j * (int)(M->rowsize)] = - dd_get_d(M->matrix[i][j + 1]);
      			}
      		}
      		mxSetField(P, 0, "A", tmpa);
      		mxSetField(P, 0, "B", tmpb);
      		return P;
	}
	return 0;
}

void
hull(int nlhs, mxArray * plhs[], int nrhs, const mxArray * prhs[])
{
	dd_PolyhedraPtr P;
	dd_ErrorType err;
	dd_MatrixPtr H,V;
	
	if (nrhs  == 1 && nlhs == 1 && mxIsStruct(prhs[0])) {
		V = FT_get_V_MatrixPtr(prhs[0]);		
		dd_set_global_constants();  /* First, this must be called. */

		P = dd_DDMatrix2Poly(V, &err); /* compute the second representation */
		if (err == dd_NoError) {
			H = dd_CopyInequalities(P);
			plhs[0] = FT_set_H_MatrixPtr(H);
			dd_FreeMatrix(H);
		} else {
    			dd_WriteErrorMessages(stdout,err);
    			mexErrMsgTxt("CDD returned an error, see above(!) for details");
  		}
		dd_FreeMatrix(V);
  		dd_FreePolyhedra(P);
  		return;
	} else {
		mexErrMsgTxt("hull expects a V input struct and produces an H output struct");
	}
}

void
v_hull_extreme(int nlhs, mxArray * plhs[], int nrhs, const mxArray * prhs[])
{
	dd_PolyhedraPtr P,P1;
	dd_ErrorType err;
	dd_MatrixPtr H,V,V1;
	
	if (nrhs  == 1 && nlhs == 1 && mxIsStruct(prhs[0])) {
		V = FT_get_V_MatrixPtr(prhs[0]);		
		dd_set_global_constants();  /* First, this must be called. */

		P = dd_DDMatrix2Poly(V, &err); /* compute the second representation */
		if (err == dd_NoError) {
			H = dd_CopyInequalities(P);
			P1 = dd_DDMatrix2Poly(H, &err); /* compute the second representation */
			if (err == dd_NoError) {
				V1 = dd_CopyGenerators(P1);
				plhs[0] = FT_set_V_MatrixPtr(V1);
				dd_FreeMatrix(V1);
			} else {
    				dd_WriteErrorMessages(stdout,err);
    				mexErrMsgTxt("CDD returned an error, see above(!) for details");    			
  			}
			dd_FreePolyhedra(P1);
			dd_FreeMatrix(H);
		} else {
    			dd_WriteErrorMessages(stdout,err);
    			mexErrMsgTxt("CDD returned an error, see above(!) for details");
  		}
		dd_FreeMatrix(V);
  		dd_FreePolyhedra(P);
  		return;
	} else {
		mexErrMsgTxt("v_hull_extreme expects a V input struct and produces a V output struct");
	}
}

void
extreme(int nlhs, mxArray * plhs[], int nrhs, const mxArray * prhs[])
{
	dd_PolyhedraPtr P;
	dd_ErrorType err;
	dd_MatrixPtr H,V;
	
	if (nrhs  == 1 && nlhs == 1 && mxIsStruct(prhs[0])) {
		dd_set_global_constants();  /* First, this must be called. */
		H = FT_get_H_MatrixPtr(prhs[0]);		

		P = dd_DDMatrix2Poly(H, &err); /* compute the second representation */
		if (err == dd_NoError) {
			V = dd_CopyGenerators(P);
			plhs[0] = FT_set_V_MatrixPtr(V);
			dd_FreeMatrix(V);
		} else {
    			dd_WriteErrorMessages(stdout,err);
    			mexErrMsgTxt("CDD returned an error, see above(!) for details");    			
  		}
		dd_FreeMatrix(H);
  		dd_FreePolyhedra(P);
  		return;
	} else {
		mexErrMsgTxt("extreme expects an H input struct and produces a V output struct");
	}
}

void
adj_extreme(int nlhs, mxArray * plhs[], int nrhs, const mxArray * prhs[])
{
	dd_PolyhedraPtr P;
	dd_ErrorType err;
	dd_MatrixPtr H,V;
	dd_SetFamilyPtr A;
	
	if (nrhs  == 1 && nlhs == 2 && mxIsStruct(prhs[0])) {
		dd_set_global_constants();  /* First, this must be called. */
		H = FT_get_H_MatrixPtr(prhs[0]);		

		P = dd_DDMatrix2Poly(H, &err); /* compute the second representation */
		if (err == dd_NoError) {
			V = dd_CopyGenerators(P);
			A = dd_CopyAdjacency(P);
			plhs[0] = FT_set_V_MatrixPtr(V);
			plhs[1] = FT_set_SetFamilyPtr(A);
			dd_FreeMatrix(V);
			dd_FreeSetFamily(A);
		} else {
    			dd_WriteErrorMessages(stdout,err);
    			mexErrMsgTxt("CDD returned an error, see above(!) for details");    			
  		}
		dd_FreeMatrix(H);
  		dd_FreePolyhedra(P);
  		return;
	} else {
		mexErrMsgTxt("adj_extreme expects an H input struct and produces a V output struct and the adjacency struct");
	}
}

void
reduce_h(int nlhs, mxArray * plhs[], int nrhs, const mxArray * prhs[])
{
	dd_ErrorType err;
	dd_MatrixPtr H,H1;
	dd_rowset red;
	
	if (nrhs  == 1 && nlhs >= 1 && nlhs <= 2 && mxIsStruct(prhs[0])) {
		dd_set_global_constants();  /* First, this must be called. */
		H = FT_get_H_MatrixPtr(prhs[0]);		
		red = dd_RedundantRows(H, &err); /* find redundant rows */
		if (err == dd_NoError) {
			/* remove the red rows */
			H1 = dd_MatrixSubmatrix(H, red);
			plhs[0] = FT_set_H_MatrixPtr(H1);
			dd_FreeMatrix(H1);
			if (nlhs == 2) {
				plhs[1] = FT_set_Set(red);
			}
		} else {
    			dd_WriteErrorMessages(stdout,err);
    			mexErrMsgTxt("CDD returned an error, see above(!) for details");
  		}
		dd_FreeMatrix(H);
		set_free(red);
  		return;
	} else {
		mexErrMsgTxt("reduce_h expects an H input struct and produces a H output struct and an optional vector of removed rows");
	}
}

void
reduce_v(int nlhs, mxArray * plhs[], int nrhs, const mxArray * prhs[])
{
	dd_ErrorType err;
	dd_MatrixPtr V,V1;
	dd_rowset red;
	
	if (nrhs  == 1 && nlhs >= 1 && nlhs <= 2 && mxIsStruct(prhs[0])) {
		dd_set_global_constants();  /* First, this must be called. */
		V = FT_get_V_MatrixPtr(prhs[0]);		
		red = dd_RedundantRows(V, &err); /* find redundant rows */
		if (err == dd_NoError) {
			/* remove the red rows */
			V1 = dd_MatrixSubmatrix(V, red);
			plhs[0] = FT_set_V_MatrixPtr(V1);
			dd_FreeMatrix(V1);
			if (nlhs == 2) {
				plhs[1] = FT_set_Set(red);
			}
		} else {
    			dd_WriteErrorMessages(stdout,err);
    			mexErrMsgTxt("CDD returned an error, see above(!) for details");
  		}
		dd_FreeMatrix(V);
		set_free(red);
  		return;
	} else {
		mexErrMsgTxt("reduce_v expects an V input struct and produces a V output struct and an optional vector of removed vertices");
	}
}

void
copy_v(int nlhs, mxArray * plhs[], int nrhs, const mxArray * prhs[])
{
	dd_MatrixPtr V;
	
	if (nrhs  == 1 && nlhs == 1 && mxIsStruct(prhs[0])) {
		dd_set_global_constants();  /* First, this must be called. */
		V = FT_get_V_MatrixPtr(prhs[0]);		
		plhs[0] = FT_set_V_MatrixPtr(V);
		dd_FreeMatrix(V);
  		return;
	} else {
		mexErrMsgTxt("copy_v expects a V input struct and produces a V output struct");
	}
}

void
copy_h(int nlhs, mxArray * plhs[], int nrhs, const mxArray * prhs[])
{
	dd_MatrixPtr H;
	
	if (nrhs  == 1 && nlhs == 1 && mxIsStruct(prhs[0])) {
		dd_set_global_constants();  /* First, this must be called. */
		H = FT_get_H_MatrixPtr(prhs[0]);		
		plhs[0] = FT_set_H_MatrixPtr(H);
		dd_FreeMatrix(H);
  		return;
	} else {
		mexErrMsgTxt("copy_h expects a H input struct and produces a H output struct");
	}
}

void
solve_lp(int nlhs, mxArray * plhs[], int nrhs, const mxArray * prhs[])
{
	/* The original LP data  m x n matrix 
	   = | b   -A  |
	     | c0  c^T |,

	where the LP to be solved is to
	maximize  c^T x  +   c0
	subj. to
	     A   x  <=  b.
	*/

	dd_ErrorType error=dd_NoError;
	dd_LPSolverType solver=dd_CrissCross;  /* either DualSimplex or CrissCross */
	dd_LPPtr lp;   /* pointer to LP data structure that is not visible by user. */
  
	dd_set_global_constants(); /* First, this must be called once to use cddlib. */


	/* Input an LP using the cdd library  */
	lp = MB_get_LP_MatrixPtr(prhs[0]);

    
	/* Solve the LP by cdd LP solver. */
	dd_LPSolve(lp,solver,&error);
	if (error!=dd_NoError) dd_WriteErrorMessages(stdout, error);
	
	/* Take the solution. */
	plhs[0] = MB_set_LPsol_MatrixPtr(lp);


	/* Free allocated spaces. */
	dd_FreeLPData(lp);
}


void
find_interior(int nlhs, mxArray * plhs[], int nrhs, const mxArray * prhs[])
{
	/* We would like to find an iterior point
	   for a polyhedron in H representation 
	     A   x  <=  b.
	*/

	dd_ErrorType error=dd_NoError;
	dd_LPSolverType solver=dd_CrissCross;  /* either DualSimplex or CrissCross */
	dd_LPPtr lp, lp1;   /* pointer to LP data structure that is not visible by user. */
	dd_MatrixPtr A;
	int j;  
  
	dd_set_global_constants(); /* First, this must be called once to use cddlib. */


	/* Input an LP using the cdd library  */
	/* lp = MB_get_LP_MatrixPtr(prhs[0]); */

	if (A=FT_get_H_MatrixPtr(prhs[0])) {
		/* set objective */
		A->objective = dd_LPmin;
		for (j = 0; j < A->colsize; j++)
			dd_set_d(A->rowvec[j],0.0);
		lp=dd_Matrix2LP(A, &error);
  		dd_FreeMatrix(A);
	}else{
		mexErrMsgTxt("Error in the setting of LP matrix.");
	}

	
	lp1=dd_MakeLPforInteriorFinding(lp);
	dd_LPSolve(lp1,solver,&error);
	if (error!=dd_NoError) dd_WriteErrorMessages(stdout, error);
	    
	/* Take the solution. */
	plhs[0] = MB_set_LPsol_MatrixPtr(lp1);


	/* Free allocated spaces. */
	dd_FreeLPData(lp);
	dd_FreeLPData(lp1);
}


void
mexFunction(int nlhs, mxArray * plhs[], int nrhs, const mxArray * prhs[])
{
	int buflen, status;
	char *input_buf;
	
	if (nrhs<1){
		mexErrMsgTxt("Input 1 must be a row vector string");
	}
		
	
	/* 1. input must be a string and row vector */
	if (mxIsChar(prhs[0]) && (mxGetM(prhs[0]) == 1)) {
	/* get the length of the input string */
		buflen = mxGetN(prhs[0]) + 1;
		/* allocate memory for input and output strings */
    		input_buf=mxCalloc(buflen, sizeof(char));
		/* copy the string data from prhs[0] into a C string input_ buf. */
    		status = mxGetString(prhs[0], input_buf, buflen);
		if (strcmp(input_buf,"hull\0") == 0) {
			hull(nlhs, plhs, nrhs -1, prhs + 1);
			return;
		}
		if (strcmp(input_buf,"extreme\0") == 0) {
			extreme(nlhs, plhs, nrhs -1, prhs + 1);
			return;
		}
		if (strcmp(input_buf,"reduce_h\0") == 0) {
			reduce_h(nlhs, plhs, nrhs -1, prhs + 1);
			return;
		}
		if (strcmp(input_buf,"reduce_v\0") == 0) {
			reduce_v(nlhs, plhs, nrhs -1, prhs + 1);
			return;
		}
		if (strcmp(input_buf,"copy_v\0") == 0) {
			copy_v(nlhs, plhs, nrhs -1, prhs + 1);
			return;
		}
		if (strcmp(input_buf,"copy_h\0") == 0) {
			copy_h(nlhs, plhs, nrhs -1, prhs + 1);
			return;
		}
		if (strcmp(input_buf,"v_hull_extreme\0") == 0) {
			v_hull_extreme(nlhs, plhs, nrhs -1, prhs + 1);
			return;
		}
		if (strcmp(input_buf,"adj_extreme\0") == 0) {
			adj_extreme(nlhs, plhs, nrhs -1, prhs + 1);
			return;
		}
		if (strcmp(input_buf,"solve_lp\0") == 0) {
			solve_lp(nlhs, plhs, nrhs -1, prhs + 1);
			return;
		}
		if (strcmp(input_buf,"find_interior\0") == 0) {
			find_interior(nlhs, plhs, nrhs -1, prhs + 1);
			return;
		}
		if (strcmp(input_buf,"version\0") == 0) {
			printf("Version %s\n", CDDMEX_VERSION);
			return;
		}
		mexErrMsgTxt("Unknown function");		
	} else {
		mexErrMsgTxt("Input 1 must be a row vector string");
	}	
}

void
adjacency(int nlhs, mxArray * plhs[], int nrhs, const mxArray * prhs[])
{
	dd_PolyhedraPtr P;
	dd_ErrorType err;
	dd_MatrixPtr V;
	dd_SetFamilyPtr A;
	
		
	if (mxIsStruct(prhs[0])) {
		V = FT_get_V_MatrixPtr(prhs[0]);		
		dd_set_global_constants();  /* First, this must be called. */

		P = dd_DDMatrix2Poly(V, &err); /* compute the second representation */
		if (err == dd_NoError) {
			A = dd_CopyInputAdjacency(P);
			plhs[0] = FT_set_SetFamilyPtr(A);
			dd_FreeSetFamily(A);
		} else {
    			dd_WriteErrorMessages(stdout,err);
    			mexErrMsgTxt("CDD returned an error, see above(!) for details");
  		}
		dd_FreeMatrix(V);
  		dd_FreePolyhedra(P);
	}
}