/SqlSpatialIndexing/Lib/Cdd/IO.cs
C# | 2168 lines | 1751 code | 272 blank | 145 comment | 377 complexity | 7cf76b086103d7ab21ee0643618c47d2 MD5 | raw file
Large files files are truncated, but you can click here to view the full file
- using System;
- using System.Collections.Generic;
- using System.Text;
- using System.IO;
-
- namespace CddSharp
- {
- public partial class Cdd
- {
- /* cddio.c: Basic Input and Output Procedures for cddlib
- written by Komei Fukuda, fukuda@ifor.math.ethz.ch
- Version 0.94, Aug. 4, 2005
- */
-
- /* cddlib : C-library of the double description method for
- computing all vertices and extreme rays of the polyhedron
- P= {x : b - A x >= 0}.
- Please read COPYING (GNU General Public Licence) and
- the manual cddlibman.tex for detail.
- */
-
- #if nemkell
- void SetInputFile(FILE **f,DataFileType inputfile,ErrorType *Error)
- {
- int opened=0,stop,quit=0;
- int i,dotpos=0,trial=0;
- char ch;
- char *tempname;
-
-
- *Error=NoError;
- while (!opened && !quit) {
- fprintf(stderr,"\n>> Input file: ");
- scanf("%s",inputfile);
- ch=getchar();
- stop=false;
- for (i=0; i<filenamelen && !stop; i++){
- ch=inputfile[i];
- switch (ch) {
- case '.':
- dotpos=i+1;
- break;
- case ';': case ' ': case '\0': case '\n': case '\t':
- stop=true;
- tempname=(char*)calloc(filenamelen,sizeof(ch));
- strncpy(tempname,inputfile,i);
- strcpy(inputfile,tempname);
- free(tempname);
- break;
- }
- }
- if ( ( *f = fopen(inputfile,"r") )!= null) {
- fprintf(stderr,"input file %s is open\n",inputfile);
- opened=1;
- *Error=NoError;
- }
- else{
- fprintf(stderr,"The file %s not found\n",inputfile);
- trial++;
- if (trial>5) {
- *Error=IFileNotFound;
- quit=1;
- }
- }
- }
- }
-
- void SetWriteFileName(DataFileType inputfile, DataFileType outfile, char cflag, RepresentationType rep)
- {
- char *extension;
- DataFileType ifilehead="";
- int i,dotpos;
-
- switch (cflag) {
- case 'o':
- switch (rep) {
- case Generator:
- extension=".ine"; break; /* output file for ine data */
- case Inequality:
- extension=".ext"; break; /* output file for ext data */
- default:
- extension=".xxx";break;
- }
- break;
-
- case 'a': /* decide for output adjacence */
- if (rep==Inequality)
- extension=".ead"; /* adjacency file for ext data */
- else
- extension=".iad"; /* adjacency file for ine data */
- break;
- case 'i': /* decide for output incidence */
- if (rep==Inequality)
- extension=".ecd"; /* ext incidence file */
- else
- extension=".icd"; /* ine incidence file */
- break;
- case 'n': /* decide for input incidence */
- if (rep==Inequality)
- extension=".icd"; /* ine incidence file */
- else
- extension=".ecd"; /* ext incidence file */
- break;
- case 'j': /* decide for input adjacency */
- if (rep==Inequality)
- extension=".iad"; /* ine adjacency file */
- else
- extension=".ead"; /* ext adjacency file */
- break;
- case 'l':
- extension=".ddl";break; /* log file */
- case 'd':
- extension=".dex";break; /* decomposition output */
- case 'p':
- extension="sub.ine";break; /* preprojection sub inequality file */
- case 'v':
- extension=".solved";break; /* verify_input file */
- case 's':
- extension=".lps";break; /* LP solution file */
- default:
- extension=".xxx";break;
- }
- dotpos=-1;
- for (i=0; i< strlen(inputfile); i++){
- if (inputfile[i]=='.') dotpos=i;
- }
- if (dotpos>1) strncpy(ifilehead, inputfile, dotpos);
- else strcpy(ifilehead,inputfile);
- if (strlen(inputfile)<=0) strcpy(ifilehead,"tempcdd");
- strcpy(outfile,ifilehead);
- strcat(outfile,extension);
- if (strcmp(inputfile, outfile)==0) {
- strcpy(outfile,inputfile);
- strcat(outfile,extension);
- }
- /* fprintf(stderr,"outfile name = %s\n",outfile); */
- }
-
-
- NumberType GetNumberType(char *line)
- {
- NumberType nt;
-
- if (strncmp(line, "integer", 7)==0) {
- nt = Integer;
- }
- else if (strncmp(line, "rational", 8)==0) {
- nt = Rational;
- }
- else if (strncmp(line, "real", 4)==0) {
- nt = Real;
- }
- else {
- nt=Unknown;
- }
- return nt;
- }
-
- void ProcessCommandLine(FILE *f, Matrix M, char *line)
- {
- char newline[linelenmax];
- int j;
- double value;
-
- init(value);
- if (strncmp(line, "hull", 4)==0) {
- M.representation = Generator;
- }
- if (strncmp(line, "debug", 5)==0) {
- debug = true;
- //#ifdef GMPRATIONAL
- ddf_debug = ddf_TRUE;
- //#endif
- }
- if (strncmp(line, "partial_enum", 12)==0 ||
- strncmp(line, "equality", 8)==0 ||
- strncmp(line, "linearity", 9)==0 ) {
- fgets(newline,linelenmax,f);
- SetLinearity(M,newline);
- }
- if (strncmp(line, "maximize", 8)==0 ||
- strncmp(line, "minimize", 8)==0) {
- if (strncmp(line, "maximize", 8)==0) M.objective=LPmax;
- else M.objective=LPmin;
- for (j = 1; j <= M.colsize; j++) {
- if (M.numbtype==Real) {
- //#if !defined(GMPRATIONAL)
- double rvalue;
- fscanf(f, "%lf", &rvalue);
- set_d(value, rvalue);
- //#endif
- } else {
- fread_rational_value (f, value);
- }
- set(M.rowvec[j - 1],value);
- if (debug) {fprintf(stderr,"cost(%5ld) =",j); WriteNumber(stderr,value);}
- } /*of j*/
- }
- clear(value);
- }
- #endif
-
- public bool AppendMatrix2Poly(Polyhedra poly, Matrix M)
- {
- bool success = false;
- Matrix Mpoly, Mnew = null;
- ErrorType err;
-
- if (poly != null && poly.m >= 0 && poly.d >= 0 &&
- poly.d == M.colsize && M.rowsize > 0)
- {
- Mpoly = CopyInput(poly);
- Mnew = AppendMatrix(Mpoly, M);
-
- poly = DDMatrix2Poly(Mnew, out err);
-
- if (err == ErrorType.NoError) success = true;
- }
- return success;
- }
-
- public Matrix MatrixCopy(Matrix M)
- {
- Matrix Mcopy = null;
- int m;
- int d;
-
- m = M.rowsize;
- d = M.colsize;
- if (m >= 0 && d >= 0)
- {
- Mcopy = CreateMatrix(m, d);
- CopyAmatrix(Mcopy.matrix, M.matrix, m, d);
- CopyArow(Mcopy.rowvec, M.rowvec, d);
- M.linset.CopyTo(Mcopy.linset);
- Mcopy.numbtype = M.numbtype;
- Mcopy.representation = M.representation;
- Mcopy.objective = M.objective;
- }
- return Mcopy;
- }
-
- public Matrix CopyMatrix(Matrix M)
- {
- return MatrixCopy(M);
- }
-
- public Matrix MatrixNormalizedCopy(Matrix M)
- {
- Matrix Mcopy = null;
- int m;
- int d;
-
- m = M.rowsize;
- d = M.colsize;
- if (m >= 0 && d >= 0)
- {
- Mcopy = CreateMatrix(m, d);
- CopyNormalizedAmatrix(Mcopy.matrix, M.matrix, m, d);
- CopyArow(Mcopy.rowvec, M.rowvec, d);
- M.linset.CopyTo(Mcopy.linset);
- Mcopy.numbtype = M.numbtype;
- Mcopy.representation = M.representation;
- Mcopy.objective = M.objective;
- }
- return Mcopy;
- }
-
-
- public Matrix MatrixAppend(Matrix M1, Matrix M2)
- {
- Matrix M = null;
- int i, m, m1, m2;
- int j, d, d1, d2;
-
- m1 = M1.rowsize;
- d1 = M1.colsize;
- m2 = M2.rowsize;
- d2 = M2.colsize;
-
- m = m1 + m2;
- d = d1;
-
- if (d1 >= 0 && d1 == d2 && m1 >= 0 && m2 >= 0)
- {
- M = CreateMatrix(m, d);
- CopyAmatrix(M.matrix, M1.matrix, m1, d);
- CopyArow(M.rowvec, M1.rowvec, d);
- for (i = 0; i < m1; i++)
- {
- if (M1.linset.Get(i + 1)) M.linset.Set(i + 1, true);
- }
- for (i = 0; i < m2; i++)
- {
- for (j = 0; j < d; j++)
- M.matrix[m1 + i][j] = M2.matrix[i][j];
- /* append the second matrix */
- if (M2.linset.Get(i + 1)) M.linset.Set(m1 + i + 1, true);
- }
- M.numbtype = M1.numbtype;
- }
- return M;
- }
-
- public Matrix MatrixNormalizedSortedCopy(Matrix M, out int[] newpos) /* 094 */
- {
- /* Sort the rows of double[][] lexicographically, and return a link to this sorted copy.
- The vector newpos is allocated, where newpos[i] returns the new row index
- of the original row i (i=1,...,M.rowsize). */
- Matrix Mcopy = null, Mnorm = null;
- int m, i;
- int d;
- int[] roworder;
-
- /* if (newpos!=null) free(newpos); */
- m = M.rowsize;
- d = M.colsize;
- roworder = new int[m + 1];
- newpos = new int[m + 1];
- if (m >= 0 && d >= 0)
- {
- Mnorm = MatrixNormalizedCopy(M);
- Mcopy = CreateMatrix(m, d);
- for (i = 1; i <= m; i++) roworder[i] = i;
-
- RandomPermutation(roworder, m, 123);
- QuickSort(roworder, 1, m, Mnorm.matrix, d);
-
- PermuteCopyAmatrix(Mcopy.matrix, Mnorm.matrix, m, d, roworder);
- CopyArow(Mcopy.rowvec, M.rowvec, d);
- for (i = 1; i <= m; i++)
- {
- if (M.linset.Get(roworder[i])) Mcopy.linset.Set(i, true);
- newpos[roworder[i]] = i;
- }
- Mcopy.numbtype = M.numbtype;
- Mcopy.representation = M.representation;
- Mcopy.objective = M.objective;
-
- }
-
- return Mcopy;
- }
-
- public Matrix MatrixUniqueCopy(Matrix M, out int[] newpos)
- {
- /* Remove row duplicates, and return a link to this sorted copy.
- Linearity rows have priority over the other rows.
- It is better to call this after sorting with MatrixNormalizedSortedCopy.
- The vector newpos is allocated, where *newpos[i] returns the new row index
- of the original row i (i=1,...,M.rowsize). *newpos[i] is negative if the original
- row is dominated by -*newpos[i] and eliminated in the new copy.
- */
- Matrix Mcopy = null;
- int m, i, uniqrows;
- BitArray preferredrows;
- int d;
- int[] roworder;
-
- /* if (newpos!=null) free(newpos); */
- m = M.rowsize;
- d = M.colsize;
- preferredrows = M.linset;
- roworder = new int[m + 1];
- if (m >= 0 && d >= 0)
- {
- for (i = 1; i <= m; i++) roworder[i] = i;
- UniqueRows(roworder, 1, m, M.matrix, d, preferredrows, out uniqrows);
-
- Mcopy = CreateMatrix(uniqrows, d);
- PermutePartialCopyAmatrix(Mcopy.matrix, M.matrix, m, d, roworder, 1, m);
- CopyArow(Mcopy.rowvec, M.rowvec, d);
- for (i = 1; i <= m; i++)
- {
- if (roworder[i] > 0 && M.linset.Get(i)) Mcopy.linset.Set(roworder[i], true);
- }
- Mcopy.numbtype = M.numbtype;
- Mcopy.representation = M.representation;
- Mcopy.objective = M.objective;
- }
- newpos = roworder;
- return Mcopy;
- }
-
-
- public Matrix MatrixNormalizedSortedUniqueCopy(Matrix M, out int[] newpos) /* 094 */
- {
- /* Sort and remove row duplicates, and return a link to this sorted copy.
- Linearity rows have priority over the other rows.
- It is better to call this after sorting with MatrixNormalizedSortedCopy.
- The vector newpos is allocated, where *newpos[i] returns the new row index
- of the original row i (i=1,...,M.rowsize). *newpos[i] is negative if the original
- row is dominated by -*newpos[i] and eliminated in the new copy.
- */
- Matrix M1 = null, M2 = null;
- int m, i;
- int d;
- int[] newpos1 = null, newpos1r = null, newpos2 = null;
-
- /* if (newpos!=null) free(newpos); */
- m = M.rowsize;
- d = M.colsize;
- newpos = new int[m + 1];
- newpos1r = new int[m + 1];
- if (m >= 0 && d >= 0)
- {
- M1 = MatrixNormalizedSortedCopy(M, out newpos1);
- for (i = 1; i <= m; i++) newpos1r[newpos1[i]] = i; /* reverse of newpos1 */
- M2 = MatrixUniqueCopy(M1, out newpos2);
- M2.linset.SetAll(false); //emptyset
- for (i = 1; i <= m; i++)
- {
- if (newpos2[newpos1[i]] > 0)
- {
- //printf("newpos1[%ld]=%ld, newpos2[newpos1[%ld]]=%ld\n",i,newpos1[i], i,newpos2[newpos1[i]]);
- if (M.linset.Get(i)) M2.linset.Set(newpos2[newpos1[i]], true);
- newpos[i] = newpos2[newpos1[i]];
- }
- else
- {
- newpos[i] = -newpos1r[-newpos2[newpos1[i]]];
- }
- }
-
- }
-
- return M2;
- }
-
- public Matrix MatrixSortedUniqueCopy(Matrix M, out int[] newpos) /* 094 */
- {
- /* Same as MatrixNormalizedSortedUniqueCopy except that it returns a unnormalized origial data
- with original ordering.
- */
- Matrix M1 = null, M2 = null;
- int m, i, k, ii;
- int d;
- int[] newpos1 = null, newpos1r = null, newpos2 = null;
-
- /* if (newpos!=null) free(newpos); */
- m = M.rowsize;
- d = M.colsize;
- newpos = new int[m + 1];
- newpos1r = new int[m + 1];
- if (m >= 0 && d >= 0)
- {
- M1 = MatrixNormalizedSortedCopy(M, out newpos1);
- for (i = 1; i <= m; i++) newpos1r[newpos1[i]] = i; /* reverse of newpos1 */
- M2 = MatrixUniqueCopy(M1, out newpos2);
-
- M2.linset.SetAll(false);
- for (i = 1; i <= m; i++)
- {
- if (newpos2[newpos1[i]] > 0)
- {
- if (M.linset.Get(i)) M2.linset.Set(newpos2[newpos1[i]], true);
- newpos[i] = newpos2[newpos1[i]];
- }
- else
- {
- newpos[i] = -newpos1r[-newpos2[newpos1[i]]];
- }
- }
-
- ii = 0;
- M2.linset.SetAll(false);
- for (i = 1; i <= m; i++)
- {
- k = newpos[i];
- if (k > 0)
- {
- ii += 1;
- newpos[i] = ii;
- CopyArow(M2.matrix[ii - 1], M.matrix[i - 1], d);
- if (M.linset.Get(i)) M2.linset.Set(ii, true);
- }
- }
-
- }
-
- return M2;
- }
-
- public Matrix AppendMatrix(Matrix M1, Matrix M2)
- {
- return MatrixAppend(M1, M2);
- }
-
- public bool MatrixAppendTo(ref Matrix M1, Matrix M2)
- {
- Matrix M = null;
- int i, m, m1, m2;
- int j, d, d1, d2;
- bool success = false;
-
- m1 = M1.rowsize;
- d1 = M1.colsize;
- m2 = M2.rowsize;
- d2 = M2.colsize;
-
- m = m1 + m2;
- d = d1;
-
- if (d1 >= 0 && d1 == d2 && m1 >= 0 && m2 >= 0)
- {
- M = CreateMatrix(m, d);
- CopyAmatrix(M.matrix, M1.matrix, m1, d);
- CopyArow(M.rowvec, M1.rowvec, d);
- for (i = 0; i < m1; i++)
- {
- if (M1.linset.Get(i + 1)) M.linset.Set(i + 1, true);
- }
- for (i = 0; i < m2; i++)
- {
- for (j = 0; j < d; j++)
- M.matrix[m1 + i][j] = M2.matrix[i][j];
- /* append the second matrix */
- if (M2.linset.Get(i + 1)) M.linset.Set(m1 + i + 1, true);
- }
- M.numbtype = M1.numbtype;
-
- M1 = M;
- success = true;
- }
- return success;
- }
-
- public bool MatrixRowRemove(Matrix M, int r) /* 092 */
- {
- int i, m;
- int d;
- bool success = false;
-
- m = M.rowsize;
- d = M.colsize;
-
- if (r >= 1 && r <= m)
- {
- M.rowsize = m - 1;
-
- M.linset.Set(r, false);
- /* slide the row headers */
- for (i = r; i < m; i++)
- {
- M.matrix[i - 1] = M.matrix[i];
- if (M.linset.Get(i + 1))
- {
- M.linset.Set(i + 1, false);
- M.linset.Set(i, true);
- }
- }
- success = true;
- }
- return success;
- }
-
- public bool MatrixRowRemove2(Matrix M, int r, out int[] newpos) /* 094 */
- {
- int i, m;
- int d;
- bool success = false;
- int[] roworder;
-
- newpos = null;
-
- m = M.rowsize;
- d = M.colsize;
-
- if (r >= 1 && r <= m)
- {
- roworder = new int[m + 1];
- M.rowsize = m - 1;
-
- M.linset.Set(r, false);
- /* slide the row headers */
- for (i = 1; i < r; i++) roworder[i] = i;
- roworder[r] = 0; /* meaning it is removed */
- for (i = r; i < m; i++)
- {
- M.matrix[i - 1] = M.matrix[i];
- roworder[i + 1] = i;
- if (M.linset.Get(i + 1))
- {
- M.linset.Set(i + 1, false);
- M.linset.Set(i, true);
- }
- }
- success = true;
- }
- return success;
- }
-
- public Matrix MatrixSubmatrix(Matrix M, BitArray delset) /* 092 */
- {
- Matrix Msub = null;
- int i, isub = 1, m, msub;
- int d;
-
- m = M.rowsize;
- d = M.colsize;
- msub = m;
- if (m >= 0 && d >= 0)
- {
- for (i = 1; i <= m; i++)
- {
- if (delset.Get(i)) msub -= 1;
- }
- Msub = CreateMatrix(msub, d);
- for (i = 1; i <= m; i++)
- {
- if (!delset.Get(i))
- {
- CopyArow(Msub.matrix[isub - 1], M.matrix[i - 1], d);
- if (M.linset.Get(i))
- {
- Msub.linset.Set(isub, true);
- }
- isub++;
- }
- }
- CopyArow(Msub.rowvec, M.rowvec, d);
- Msub.numbtype = M.numbtype;
- Msub.representation = M.representation;
- Msub.objective = M.objective;
- }
- return Msub;
- }
-
- public Matrix MatrixSubmatrix2(Matrix M, BitArray delset, out int[] newpos) /* 092 */
- { /* returns a pointer to a new matrix which is a submatrix of M with rows in delset
- removed. *newpos[i] returns the position of the original row i in the new matrix.
- It is -1 if and only if it is deleted.
- */
-
- Matrix Msub = null;
- int i, isub = 1, m, msub;
- int d;
- int[] roworder;
-
- newpos = null;
-
- m = M.rowsize;
- d = M.colsize;
- msub = m;
- if (m >= 0 && d >= 0)
- {
- roworder = new int[m + 1];
- for (i = 1; i <= m; i++)
- {
- if (delset.Get(i)) msub -= 1;
- }
- Msub = CreateMatrix(msub, d);
- for (i = 1; i <= m; i++)
- {
- if (delset.Get(i))
- {
- roworder[i] = 0; /* zero means the row i is removed */
- }
- else
- {
- CopyArow(Msub.matrix[isub - 1], M.matrix[i - 1], d);
- if (M.linset.Get(i))
- {
- Msub.linset.Set(isub, true);
- }
- roworder[i] = isub;
- isub++;
- }
- }
- newpos = roworder;
- CopyArow(Msub.rowvec, M.rowvec, d);
- Msub.numbtype = M.numbtype;
- Msub.representation = M.representation;
- Msub.objective = M.objective;
- }
- return Msub;
- }
-
-
- public Matrix MatrixSubmatrix2L(Matrix M, BitArray delset, out int[] newpos) /* 094 */
- { /* This is same as MatrixSubmatrix2 except that the linearity rows will be shifted up
- so that they are at the top of the matrix.
- */
- Matrix Msub = null;
- int i, iL, iI, m, msub;
- int d;
- int[] roworder;
-
- newpos = null;
-
- m = M.rowsize;
- d = M.colsize;
- msub = m;
- if (m >= 0 && d >= 0)
- {
- roworder = new int[m + 1];
- for (i = 1; i <= m; i++)
- {
- if (delset.Get(i)) msub -= 1;
- }
- Msub = CreateMatrix(msub, d);
- iL = 1; iI = M.linset.Count + 1; /* starting positions */
- for (i = 1; i <= m; i++)
- {
- if (delset.Get(i))
- {
- roworder[i] = 0; /* zero means the row i is removed */
- }
- else
- {
- if (M.linset.Get(i))
- {
- CopyArow(Msub.matrix[iL - 1], M.matrix[i - 1], d);
- Msub.linset.Set(i, false);
- Msub.linset.Set(iL, true);
- roworder[i] = iL;
- iL += 1;
- }
- else
- {
- CopyArow(Msub.matrix[iI - 1], M.matrix[i - 1], d);
- roworder[i] = iI;
- iI += 1;
- }
- }
- }
- newpos = roworder;
- CopyArow(Msub.rowvec, M.rowvec, d);
- Msub.numbtype = M.numbtype;
- Msub.representation = M.representation;
- Msub.objective = M.objective;
- }
- return Msub;
- }
-
- public int MatrixRowsRemove(ref Matrix M, BitArray delset) /* 094 */
- {
- Matrix Msub = null;
- int success;
-
- Msub = MatrixSubmatrix(M, delset);
-
- M = Msub;
- success = 1;
- return success;
- }
-
- public int MatrixRowsRemove2(ref Matrix M, BitArray delset, out int[] newpos) /* 094 */
- {
- Matrix Msub = null;
- int success;
-
- Msub = MatrixSubmatrix2(M, delset, out newpos);
-
- M = Msub;
- success = 1;
- return success;
- }
-
- public int MatrixShiftupLinearity(ref Matrix M, out int[] newpos) /* 094 */
- {
- Matrix Msub = null;
- int success;
- BitArray delset;
-
- delset = new BitArray(M.rowsize); /* emptyset */
- Msub = MatrixSubmatrix2L(M, delset, out newpos);
-
- M = Msub;
-
- success = 1;
- return success;
- }
-
- public Polyhedra CreatePolyhedraData(int m, int d)
- {
- int i;
- Polyhedra poly = null;
-
- poly = new Polyhedra();
- poly.child = null; /* this links the homogenized cone data */
- poly.m = m;
- poly.d = d;
- poly.n = -1; /* the size of output is not known */
- poly.m_alloc = m + 2; /* the allocated row size of matrix A */
- poly.d_alloc = d; /* the allocated col size of matrix A */
- poly.numbtype = NumberType.Real;
- InitializeAmatrix(poly.m_alloc, poly.d_alloc, out poly.A);
- InitializeArow(d, out poly.c); /* cost vector */
- poly.representation = RepresentationType.Inequality;
- poly.homogeneous = false;
-
- poly.EqualityIndex = new int[m + 2];
- /* size increased to m+2 in 092b because it is used by the child cone,
- This is a bug fix suggested by Thao Dang. */
- /* ith component is 1 if it is equality, -1 if it is strict inequality, 0 otherwise. */
- for (i = 0; i <= m + 1; i++) poly.EqualityIndex[i] = 0;
-
- poly.IsEmpty = -1; /* initially set to -1, neither true nor false, meaning unknown */
-
- poly.NondegAssumed = false;
- poly.InitBasisAtBottom = false;
- poly.RestrictedEnumeration = false;
- poly.RelaxedEnumeration = false;
-
- poly.AincGenerated = false; /* Ainc is a set array to store the input incidence. */
-
- return poly;
- }
-
- public bool InitializeConeData(int m, int d, out Cone cone)
- {
- bool success = true;
- int j;
-
- cone = new Cone();
-
- /* INPUT: A given representation of a cone: inequality */
- cone.m = m;
- cone.d = d;
- cone.m_alloc = m + 2; /* allocated row size of matrix A */
- cone.d_alloc = d; /* allocated col size of matrix A, B and Bsave */
- cone.numbtype = NumberType.Real;
- cone.parent = null;
-
- /* CONTROL: variables to control computation */
- cone.Iteration = 0;
-
- cone.HalfspaceOrder = RowOrderType.LexMin;
-
- cone.ArtificialRay = null;
- cone.FirstRay = null;
- cone.LastRay = null; /* The second description: Generator */
- cone.PosHead = null;
- cone.ZeroHead = null;
- cone.NegHead = null;
- cone.PosLast = null;
- cone.ZeroLast = null;
- cone.NegLast = null;
- cone.RecomputeRowOrder = true;
- cone.PreOrderedRun = false;
- cone.GroundSet = new BitArray(cone.m_alloc);
- cone.EqualitySet = new BitArray(cone.m_alloc);
- cone.NonequalitySet = new BitArray(cone.m_alloc);
- cone.AddedHalfspaces = new BitArray(cone.m_alloc);
- cone.WeaklyAddedHalfspaces = new BitArray(cone.m_alloc);
- cone.InitialHalfspaces = new BitArray(cone.m_alloc);
- cone.RayCount = 0;
- cone.FeasibleRayCount = 0;
- cone.WeaklyFeasibleRayCount = 0;
- cone.TotalRayCount = 0;
- cone.ZeroRayCount = 0;
- cone.EdgeCount = 0;
- cone.TotalEdgeCount = 0;
- cone.count_int = 0;
- cone.count_int_good = 0;
- cone.count_int_bad = 0;
- cone.rseed = 1; /* random seed for random row permutation */
-
- InitializeBmatrix(cone.d_alloc, out cone.B);
- InitializeBmatrix(cone.d_alloc, out cone.Bsave);
- InitializeAmatrix(cone.m_alloc, cone.d_alloc, out cone.A);
-
- cone.Edges = new Adjacency[cone.m_alloc];
- cone.InitialRayIndex = new int[d + 1];
- cone.OrderVector = new int[cone.m_alloc + 1];
-
- cone.newcol = new int[cone.d + 1];
- for (j = 0; j <= cone.d; j++) cone.newcol[j] = j; /* identity map, initially */
- cone.LinearityDim = -2; /* -2 if it is not computed */
- cone.ColReduced = false;
- cone.d_orig = d;
-
- /* STATES: variables to represent current state. */
- /*cone.Error;
- cone.CompStatus;
- cone.starttime;
- cone.endtime;
- */
-
- return success;
- }
-
- public Cone ConeDataLoad(Polyhedra poly)
- {
- Cone cone = null;
- int d, j;
- int m, i;
-
- m = poly.m;
- d = poly.d;
- if (!(poly.homogeneous) && poly.representation == RepresentationType.Inequality)
- {
- m = poly.m + 1;
- }
- poly.m1 = m;
-
- InitializeConeData(m, d, out cone);
- cone.representation = poly.representation;
-
- /* Points to the original polyhedra data, and reversely */
- cone.parent = poly;
- poly.child = cone;
-
- for (i = 1; i <= poly.m; i++)
- for (j = 1; j <= cone.d; j++)
- cone.A[i - 1][j - 1] = poly.A[i - 1][j - 1];
-
- if (poly.representation == RepresentationType.Inequality && !(poly.homogeneous))
- {
- cone.A[m - 1][0] = 1.0;
- for (j = 2; j <= d; j++) cone.A[m - 1][j - 1] = 0.0;
- }
-
- return cone;
- }
-
- #if nemkell
- void SetLinearity(Matrix M, char *line)
- {
- int i=0;
- int eqsize,var;
- char *next;
- const char ct[]=", "; /* allows separators "," and " ". */
-
- next=strtok(line,ct);
- eqsize=atol(next);
- while (i < eqsize && (next=strtok(null,ct))!=null) {
- var=atol(next);
- set_addelem(M.linset,var); i++;
- }
- if (i!=eqsize) {
- fprintf(stderr,"* Warning: there are inconsistencies in linearity setting.\n");
- }
- return;
- }
-
- Matrix PolyFile2Matrix (FILE *f, ErrorType *Error)
- {
- Matrix M=null;
- int m_input,i;
- int d_input,j;
- RepresentationType rep=Inequality;
- double value;
- bool found=false, newformat=false, successful=false, linearity=false;
- char command[linelenmax], comsave[linelenmax], numbtype[wordlenmax];
- NumberType NT;
- //#if !defined(GMPRATIONAL)
- double rvalue;
- //#endif
-
- init(value);
- (*Error)=NoError;
- while (!found)
- {
- if (fscanf(f,"%s",command)==EOF) {
- (*Error)=ImproperInputFormat;
- goto _L99;
- }
- else {
- if (strncmp(command, "V-representation", 16)==0) {
- rep=Generator; newformat=true;
- }
- if (strncmp(command, "H-representation", 16)==0){
- rep=Inequality; newformat=true;
- }
- if (strncmp(command, "partial_enum", 12)==0 ||
- strncmp(command, "equality", 8)==0 ||
- strncmp(command, "linearity", 9)==0 ) {
- linearity=true;
- fgets(comsave,linelenmax,f);
- }
- if (strncmp(command, "begin", 5)==0) found=true;
- }
- }
- fscanf(f, "%ld %ld %s", &m_input, &d_input, numbtype);
- fprintf(stderr,"size = %ld x %ld\nNumber Type = %s\n", m_input, d_input, numbtype);
- NT=GetNumberType(numbtype);
- if (NT==Unknown) {
- (*Error)=ImproperInputFormat;
- goto _L99;
- }
- M=CreateMatrix(m_input, d_input);
- M.representation=rep;
- M.numbtype=NT;
-
- for (i = 1; i <= m_input; i++) {
- for (j = 1; j <= d_input; j++) {
- if (NT==Real) {
- //#if defined GMPRATIONAL
- *Error=NoRealNumberSupport;
- goto _L99;
- //#else
- fscanf(f, "%lf", &rvalue);
- set_d(value, rvalue);
- //#endif
- } else {
- fread_rational_value (f, value);
- }
- set(M.matrix[i-1][j - 1],value);
- if (debug) {fprintf(stderr,"a(%3ld,%5ld) = ",i,j); WriteNumber(stderr,value);}
- } /*of j*/
- } /*of i*/
- if (fscanf(f,"%s",command)==EOF) {
- (*Error)=ImproperInputFormat;
- goto _L99;
- }
- else if (strncmp(command, "end", 3)!=0) {
- if (debug) fprintf(stderr,"'end' missing or illegal extra data: %s\n",command);
- (*Error)=ImproperInputFormat;
- goto _L99;
- }
-
- successful=true;
- if (linearity) {
- SetLinearity(M,comsave);
- }
- while (!feof(f)) {
- fscanf(f,"%s", command);
- ProcessCommandLine(f, M, command);
- fgets(command,linelenmax,f); /* skip the CR/LF */
- }
-
- _L99: ;
- clear(value);
- /* if (f!=null) fclose(f); */
- return M;
- }
-
- #endif
-
- public Polyhedra DDMatrix2Poly(Matrix M, out ErrorType err)
- {
- int i;
- int j;
- Polyhedra poly = null;
-
- err = ErrorType.NoError;
- if (M.rowsize < 0 || M.colsize < 0)
- {
- err = ErrorType.NegativeMatrixSize;
- return null;
- }
- poly = CreatePolyhedraData(M.rowsize, M.colsize);
- poly.representation = M.representation;
- poly.homogeneous = true;
-
- for (i = 1; i <= M.rowsize; i++)
- {
- if (M.linset.Get(i))
- {
- poly.EqualityIndex[i] = 1;
- }
- for (j = 1; j <= M.colsize; j++)
- {
- poly.A[i - 1][j - 1] = M.matrix[i - 1][j - 1];
- if (j == 1 && M.matrix[i - 1][j - 1] != 0.0) poly.homogeneous = false;
- } /*of j*/
- } /*of i*/
- DoubleDescription(poly, out err);
-
- return poly;
- }
-
- public Polyhedra DDMatrix2Poly2(Matrix M, RowOrderType horder, out ErrorType err)
- {
- int i;
- int j;
- Polyhedra poly = null;
-
- err = ErrorType.NoError;
- if (M.rowsize < 0 || M.colsize < 0)
- {
- err = ErrorType.NegativeMatrixSize;
- return null;
- }
- poly = CreatePolyhedraData(M.rowsize, M.colsize);
- poly.representation = M.representation;
- poly.homogeneous = true;
-
- for (i = 1; i <= M.rowsize; i++)
- {
- if (M.linset.Get(i))
- {
- poly.EqualityIndex[i] = 1;
- }
- for (j = 1; j <= M.colsize; j++)
- {
- poly.A[i - 1][j - 1] = M.matrix[i - 1][j - 1];
- if (j == 1 && M.matrix[i - 1][j - 1] != 0.0) poly.homogeneous = false;
- } /*of j*/
- } /*of i*/
- DoubleDescription2(poly, horder, out err);
-
- return poly;
- }
-
- void MatrixIntegerFilter(Matrix M)
- { /* setting an almost integer to the integer. */
- int i;
- int j;
- double x;
-
- x = 0.0;
- for (i = 0; i < M.rowsize; i++)
- for (j = 0; j < M.colsize; j++)
- {
- SnapToInteger(out x, M.matrix[i][j]);
- M.matrix[i][j] = x;
- }
-
- }
-
- void CopyRay(double[] a, int d_origsize, Ray RR,
- RepresentationType rep, int[] reducedcol)
- {
- int j, j1;
- double b;
-
- b = 0.0;
- for (j = 1; j <= d_origsize; j++)
- {
- j1 = reducedcol[j];
- if (j1 > 0)
- {
- a[j - 1] = RR.RRay[j1 - 1];
- /* the original column j is mapped to j1, and thus
- copy the corresponding component */
- }
- else
- {
- a[j - 1] = 0.0;
- /* original column is redundant and removed for computation */
- }
- }
-
- b = a[0];
- if (rep == RepresentationType.Generator && b != 0)
- {
- a[0] = 1.0;
- for (j = 2; j <= d_origsize; j++)
- a[j - 1] /= b; /* normalization for generators */
- }
- }
-
-
- void WriteRay(TextWriter f, int d_origsize, Ray RR, RepresentationType rep, int[] reducedcol)
- {
- int j;
- double[] a;
-
-
- InitializeArow(d_origsize + 1, out a);
-
- CopyRay(a, d_origsize, RR, rep, reducedcol);
- for (j = 0; j < d_origsize; j++) WriteNumber(f, a[j]);
- f.WriteLine();
- }
-
- void WriteArow(TextWriter f, double[] a, int d)
- {
- int j;
-
- for (j = 0; j < d; j++) WriteNumber(f, a[j]);
- f.WriteLine();
- }
-
- void WriteAmatrix(TextWriter f, double[][] A, int rowmax, int colmax)
- {
- int i, j;
-
- if (A == null)
- {
- f.WriteLine("WriteAmatrix: The requested matrix is empty");
- return;
- }
- f.WriteLine("begin");
- #if GMPRATIONAL
- fprintf(f, " %ld %ld rational\n",rowmax, colmax);
- #else
- f.WriteLine(" {0} {1} real", rowmax, colmax);
- #endif
- for (i = 1; i <= rowmax; i++)
- {
- for (j = 1; j <= colmax; j++)
- {
- WriteNumber(f, A[i - 1][j - 1]);
- }
- f.WriteLine();
- }
- f.WriteLine("end");
-
- }
-
- void WriteBmatrix(TextWriter f, int d_size, double[][] B)
- {
- int j1, j2;
-
- if (B == null)
- {
- f.WriteLine("WriteBmatrix: The requested matrix is empty");
- return;
- }
- for (j1 = 0; j1 < d_size; j1++)
- {
- for (j2 = 0; j2 < d_size; j2++)
- {
- WriteNumber(f, B[j1][j2]);
- } /*of j2*/
- f.WriteLine();
- } /*of j1*/
- f.WriteLine();
-
- }
- #if nemkell
- void WriteSetFamily(TextWriter f, SetFamily F)
- {
- int i;
-
- if (F==null){
- f.WriteLine("WriteSetFamily: The requested family is empty");
- return;
- }
- f.WriteLine("begin");
- f.WriteLine(" {0} {1}", F.famsize, F.setsize);
- for (i=0; i<F.famsize; i++) {
- f.Write(" {0} {1} : ", i+1, F.set[i].Count);
- //*****set_fwrite(f, F.set[i]);
- }
- f.WriteLine("end");
-
- }
-
- void WriteSetFamilyCompressed(FILE *f, SetFamily F)
- {
- int i,card;
-
- if (F==null){
- fprintf(f, "WriteSetFamily: The requested family is empty\n");
- goto _L99;
- }
- fprintf(f,"begin\n");
- fprintf(f," %ld %ld\n", F.famsize, F.setsize);
- for (i=0; i<F.famsize; i++) {
- card=set_card(F.set[i]);
- if (F.setsize - card >= card){
- fprintf(f, " %ld %ld : ", i+1, card);
- set_fwrite(f, F.set[i]);
- } else {
- fprintf(f, " %ld %ld : ", i+1, -card);
- set_fwrite_compl(f, F.set[i]);
- }
- }
- fprintf(f,"end\n");
- _L99:;
- }
- #endif
- void WriteMatrix(TextWriter f, Matrix M)
- {
- int i, linsize;
-
- if (M == null)
- {
- f.WriteLine("WriteMatrix: The requested matrix is empty");
- return;
- }
- switch (M.representation)
- {
- case RepresentationType.Inequality:
- f.WriteLine("H-representation"); break;
- case RepresentationType.Generator:
- f.WriteLine("V-representation"); break;
- case RepresentationType.Unspecified:
- break;
- }
- linsize = M.linset.Count;
- if (linsize > 0)
- {
- f.Write("linearity {0} ", linsize);
- for (i = 1; i <= M.rowsize; i++)
- if (M.linset.Get(i)) f.Write(" {0}", i);
- f.WriteLine();
- }
- WriteAmatrix(f, M.matrix, M.rowsize, M.colsize);
- if (M.objective != LPObjectiveType.LPnone)
- {
- if (M.objective == LPObjectiveType.LPmax)
- f.WriteLine("maximize");
- else
- f.WriteLine("minimize");
- WriteArow(f, M.rowvec, M.colsize);
- }
-
- }
-
- void WriteLP(TextWriter f, LP lp)
- {
- if (lp == null)
- {
- f.WriteLine("WriteLP: The requested lp is empty");
- return;
- }
- f.WriteLine("H-representation");
-
- WriteAmatrix(f, lp.A, (lp.m) - 1, lp.d);
- if (lp.objective != LPObjectiveType.LPnone)
- {
- if (lp.objective == LPObjectiveType.LPmax)
- f.WriteLine("maximize");
- else
- f.WriteLine("minimize");
- WriteArow(f, lp.A[lp.objrow - 1], lp.d);
- }
-
- }
-
-
- void SnapToInteger(out double y, double x)
- {
- /* this is broken. It does nothing. */
- y = x;
- }
-
-
- void WriteReal(TextWriter f, double x)
- {
- /*
- int ix1,ix2,ix;
- double ax;
-
- ax=x;
- ix1= (int) (Math.Abs(ax) * 10000.0 + 0.5);
- ix2= (int) (Math.Abs(ax) + 0.5);
- ix2= ix2*10000;
- if ( ix1 == ix2) {
- if (x > 0) {
- ix = (int)(ax + 0.5);
- } else {
- ix = (int)(-ax + 0.5);
- ix = -ix;
- }
- f.Write(" {0}", ix); //%2ld
- } else
- f.Write(" {0}.9E",ax); //% .9E
- * */
-
- f.Write(x);
- }
-
- void WriteNumber(TextWriter f, double x)
- {
- #if GMPRATIONAL
- mpz_t zn,zd;
-
- mpz_init(zn); mpz_init(zd);
- mpq_canonicalize(x);
- mpq_get_num(zn,x);
- mpq_get_den(zd,x);
- fprintf(f," ");
- if (mpz_sgn(zn)==0){
- fprintf(f,"0");
- } else if (mpz_cmp_ui(zd,1U)==0){
- mpz_out_str(f,10,zn);
- } else {
- mpz_out_str(f,10,zn);fprintf(f,"/");mpz_out_str(f,10,zd);
- }
- mpz_clear(zn); mpz_clear(zd);
- #else
- WriteReal(f, x);
- #endif
- }
-
- #if nemkell
- void WriteIncidence(FILE *f, Polyhedra poly)
- {
- SetFamily I;
-
- switch (poly.representation) {
- case Inequality:
- fprintf(f, "ecd_file: Incidence of generators and inequalities\n");
- break;
- case Generator:
- fprintf(f, "icd_file: Incidence of inequalities and generators\n");
- break;
-
- default:
- break;
- }
- I=CopyIncidence(poly);
- WriteSetFamilyCompressed(f,I);
- FreeSetFamily(I);
- }
-
- void WriteAdjacency(FILE *f, Polyhedra poly)
- {
- SetFamily A;
-
- switch (poly.representation) {
- case Inequality:
- fprintf(f, "ead_file: Adjacency of generators\n");
- break;
- case Generator:
- fprintf(f, "iad_file: Adjacency of inequalities\n");
- break;
-
- default:
- break;
- }
- A=CopyAdjacency(poly);
- WriteSetFamilyCompressed(f,A);
- FreeSetFamily(A);
- }
- #endif
-
- void ComputeAinc(Polyhedra poly)
- {
- /* This generates the input incidence array poly.Ainc, and
- two sets: poly.Ared, poly.Adom.
- */
- int k;
- int i, m1;
- int j;
- bool redundant;
- Matrix M = null;
- double sum, temp;
-
- sum = 0.0; temp = 0.0;
- if (poly.AincGenerated == true) return;
-
- M = CopyOutput(poly);
- poly.n = M.rowsize;
- m1 = poly.m1;
- /* this number is same as poly.m, except when
- poly is given by nonhomogeneous inequalty:
- !(poly.homogeneous) && poly.representation==Inequality,
- it is poly.m+1. See ConeDataLoad.
- */
- poly.Ainc = new BitArray[m1];
- for (i = 1; i <= m1; i++) poly.Ainc[i - 1] = new BitArray(poly.n);
- poly.Ared = new BitArray(m1);
- poly.Adom = new BitArray(m1);
-
- for (k = 1; k <= poly.n; k++)
- {
- for (i = 1; i <= poly.m; i++)
- {
- sum = 0.0;
- for (j = 1; j <= poly.d; j++)
- {
- temp = poly.A[i - 1][j - 1] * M.matrix[k - 1][j - 1];
- sum += temp;
- }
- if (sum == 0.0)
- {
- poly.Ainc[i - 1].Set(k, true);
- }
- }
- if (!(poly.homogeneous) && poly.representation == RepresentationType.Inequality)
- {
- if (M.matrix[k - 1][0] == 0.0)
- {
- poly.Ainc[m1 - 1].Set(k, true); /* added infinity inequality (1,0,0,...,0) */
- }
- }
- }
-
- for (i = 1; i <= m1; i++)
- {
- if (poly.Ainc[i - 1].Count == M.rowsize)
- {
- poly.Adom.Set(i, true);
- }
- }
- for (i = m1; i >= 1; i--)
- {
- if (poly.Ainc[i - 1].Count == 0)
- {
- redundant = true;
- poly.Ared.Set(i, true);
- }
- else
- {
- redundant = false;
- for (k = 1; k <= m1; k++)
- {
- if (k != i && !poly.Ared.Get(k) && !poly.Adom.Get(k) &&
- poly.Ainc[i - 1].IsSubsetOf(poly.Ainc[k - 1]))
- { //****
- if (!redundant)
- {
- redundant = true;
- }
- poly.Ared.Set(i, true);
- }
- }
- }
- }
-
- poly.AincGenerated = true;
-
-
- }
-
- bool InputAdjacentQ(Polyhedra poly,
- int i1, int i2)
- /* Before calling this function, RedundantSet must be
- a set of row indices whose removal results in a minimal
- nonredundant system to represent the input polyhedron,
- DominantSet must be the set of row indices which are
- active at every extreme points/rays.
- */
- {
- bool adj = true;
- int i;
- BitArray common;
- int lastn = 0;
-
- if (poly.AincGenerated == false) ComputeAinc(poly);
- if (lastn != poly.n)
- {
- //if (lastn >0) set_free(common);
- common = new BitArray(poly.n);
- lastn = poly.n;
- }
- if (poly.Ared.Get(i1) || poly.Ared.Get(i2))
- {
- adj = false;
- return adj;
- }
- if (poly.Adom.Get(i1) || poly.Adom.Get(i2))
- {
- // dominant inequality is considered adjacencent to all others.
- adj = true;
- return adj;
- }
- common = poly.Ainc[i1 - 1].And(poly.Ainc[i2 - 1]); //**** intersection
- i = 0;
- while (i < poly.m1 && adj == true)
- {
- i++;
- if (i != i1 && i != i2 && !poly.Ared.Get(i) &&
- !poly.Adom.Get(i) && common.IsSubsetOf(poly.Ainc[i - 1]))
- {
- adj = false;
- }
- }
-
- return adj;
- }
-
- #if nemkell
- void WriteInputIncidence(FILE *f, Polyhedra poly)
- {
- SetFamily I;
-
- if (poly.AincGenerated==false) ComputeAinc(poly);
- switch (poly.representation) {
- case Inequality:
- fprintf(f,"icd_file: Incidence of inequalities and generators\n");
- break;
-
- case Generator:
- fprintf(f,"ecd_file: Incidence of generators and inequalities\n");
- break;
-
- default:
- break;
- }
-
- I=CopyInputIncidence(poly);
- WriteSetFamilyCompressed(f,I);
- FreeSetFamily(I);
- }
-
- void WriteInputAdjacency(FILE *f, Polyhedra poly)
- {
- SetFamily A;
-
- if (poly.AincGenerated==false){
- ComputeAinc(poly);
- }
- switch (poly.representation) {
- case Inequality:
- fprintf(f, "iad_file: Adjacency of inequalities\n");
- break;
-
- case Generator:
- fprintf(f, "ead_file: Adjacency of generators\n");
- break;
-
- default:
- break;
- }
- A=CopyInputAdjacency(poly);
- WriteSetFamilyCompressed(f,A);
- FreeSetFamily(A);
- }
-
-
- void WriteProgramDescription(FILE *f)
- {
- fprintf(f, "* cddlib: a double description library:%s\n", DDVERSION);
- fprintf(f, "* compiled for %s arithmetic.\n", ARITHMETIC);
- fprintf(f,"* %s\n",CO…
Large files files are truncated, but you can click here to view the full file