cddmex.c

CheckMate is a MATLAB-based tool for modeling, simulating and investigating properties of hybrid dyn

  		}
		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;
	dd_SetFamilyPtr GI, GA;

	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);
						
			GI=dd_CopyInputIncidence(P);
			GA=dd_CopyAdjacency(P);

			plhs[0] = FT_set_V_MatrixPtr(V);
			plhs[1] = ZH_set_Vlist(GI,GA);
            dd_FreeSetFamily(GA);
            dd_FreeSetFamily(GI);
			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
solve_lp_DS(int nlhs, mxArray * plhs[], int nrhs, const mxArray * prhs[])
{
  /* uses Dual Simplex method */
	/* 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_DualSimplex;
	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
find_interior_DS(int nlhs, mxArray * plhs[], int nrhs, const mxArray * prhs[])
{
  /* uses Dual Simplex method */
	/* 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_DualSimplex;
	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 
file_ine(int nlhs, mxArray * plhs[], int nrhs, const mxArray * prhs[])
/* Ine file input, V output, similar to extreme */
{
    dd_PolyhedraPtr poly;
    dd_MatrixPtr M;
    dd_ErrorType err;
    char *inputfile;
    FILE *reading=NULL;
    dd_MatrixPtr A, G;
    dd_SetFamilyPtr GI,GA;
    int buflen, status;
    
    dd_set_global_constants();  /* First, this must be called. */
    if (nrhs  == 1 && nlhs <=2 && mxIsChar(prhs[0])) {
        /*  dd_SetInputFile(&reading,inputfile, &err); */
        buflen = mxGetN(prhs[0]) + 1;
        inputfile= mxCalloc(buflen, sizeof(char));
        status = mxGetString(prhs[0], inputfile, buflen);
        if ( (reading = fopen(inputfile,"r") )== NULL) {
            mxErrMsgTxt("Input file not found\n");
            return;
        }
        printf(" Input file opened. \n");    
        M=dd_PolyFile2Matrix(reading, &err);
        
        if (err==dd_NoError) {
            poly=dd_DDMatrix2Poly(M, &err); /* compute the second representation */
            if (err!=dd_NoError) {
                dd_WriteErrorMessages(stdout,err);
                mxErrMsgTxt("CDD internal error\n");
                return;
            }
            A=dd_CopyInequalities(poly);
            G=dd_CopyGenerators(poly);
            GI=dd_CopyInputIncidence(poly);
            GA=dd_CopyAdjacency(poly);
            plhs[0] = FT_set_V_MatrixPtr(G);
            plhs[1] = ZH_set_Vlist(GI,GA);
            dd_FreePolyhedra(poly);
            dd_FreeMatrix(M);
            return;
        }
    }
    else {
        mexErrMsgTxt("file-ine expects an file input");
    }
    return;
}

void implicit_linear(int nlhs, mxArray * plhs[], int nrhs, const mxArray * prhs[])
{
    dd_MatrixPtr H;
    dd_rowset rows;
    dd_ErrorType err;
    if (nrhs  == 1 && nlhs <= 2 && mxIsStruct(prhs[0])) {
        H = FT_get_H_MatrixPtr(prhs[0]);		
        dd_set_global_constants();  /* First, this must be called. */
        rows = dd_ImplicitLinearityRows(H, &err);

        if (err == dd_NoError) {
            int total, element, i;
            double* pr;

            total = set_card(rows);
            plhs[0] = mxCreateDoubleMatrix(total, 1, mxREAL);
            pr = mxGetPr(plhs[0]);
            
            for (i=1, element=0; i<=rows[0]; i++){
                if (set_member(i, rows)) {
                    pr[element++] = (double)i;
                }
            }
        } else {
            dd_WriteErrorMessages(stdout,err);
            mexErrMsgTxt("CDD returned an error, see above(!) for details");
        }
        dd_FreeMatrix(H);
        set_free(rows);
        return;
    } else {
        mexErrMsgTxt("hull expects a V input struct and produces an H output struct");
    }
    
}

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,"solve_lp_DS\0") == 0) {
            solve_lp_DS(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,"find_interior_DS\0") == 0) {
            find_interior_DS(nlhs, plhs, nrhs -1, prhs + 1);
            return;
        }
        if (strcmp(input_buf,"version\0") == 0) {
            printf("Version %s\n", CDDMEX_VERSION);
            return;
        }
        if (strcmp(input_buf,"file-ine\0") ==0){
            file_ine(nlhs, plhs, nrhs -1, prhs +1);
            return;
        }
        if (strcmp(input_buf,"implicit_linear\0") ==0){
            implicit_linear(nlhs, plhs, nrhs -1, prhs +1);
            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);
	}
}