PageRenderTime 72ms CodeModel.GetById 24ms RepoModel.GetById 1ms app.codeStats 0ms

/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

  1. using System;
  2. using System.Collections.Generic;
  3. using System.Text;
  4. using System.IO;
  5. namespace CddSharp
  6. {
  7. public partial class Cdd
  8. {
  9. /* cddio.c: Basic Input and Output Procedures for cddlib
  10. written by Komei Fukuda, fukuda@ifor.math.ethz.ch
  11. Version 0.94, Aug. 4, 2005
  12. */
  13. /* cddlib : C-library of the double description method for
  14. computing all vertices and extreme rays of the polyhedron
  15. P= {x : b - A x >= 0}.
  16. Please read COPYING (GNU General Public Licence) and
  17. the manual cddlibman.tex for detail.
  18. */
  19. #if nemkell
  20. void SetInputFile(FILE **f,DataFileType inputfile,ErrorType *Error)
  21. {
  22. int opened=0,stop,quit=0;
  23. int i,dotpos=0,trial=0;
  24. char ch;
  25. char *tempname;
  26. *Error=NoError;
  27. while (!opened && !quit) {
  28. fprintf(stderr,"\n>> Input file: ");
  29. scanf("%s",inputfile);
  30. ch=getchar();
  31. stop=false;
  32. for (i=0; i<filenamelen && !stop; i++){
  33. ch=inputfile[i];
  34. switch (ch) {
  35. case '.':
  36. dotpos=i+1;
  37. break;
  38. case ';': case ' ': case '\0': case '\n': case '\t':
  39. stop=true;
  40. tempname=(char*)calloc(filenamelen,sizeof(ch));
  41. strncpy(tempname,inputfile,i);
  42. strcpy(inputfile,tempname);
  43. free(tempname);
  44. break;
  45. }
  46. }
  47. if ( ( *f = fopen(inputfile,"r") )!= null) {
  48. fprintf(stderr,"input file %s is open\n",inputfile);
  49. opened=1;
  50. *Error=NoError;
  51. }
  52. else{
  53. fprintf(stderr,"The file %s not found\n",inputfile);
  54. trial++;
  55. if (trial>5) {
  56. *Error=IFileNotFound;
  57. quit=1;
  58. }
  59. }
  60. }
  61. }
  62. void SetWriteFileName(DataFileType inputfile, DataFileType outfile, char cflag, RepresentationType rep)
  63. {
  64. char *extension;
  65. DataFileType ifilehead="";
  66. int i,dotpos;
  67. switch (cflag) {
  68. case 'o':
  69. switch (rep) {
  70. case Generator:
  71. extension=".ine"; break; /* output file for ine data */
  72. case Inequality:
  73. extension=".ext"; break; /* output file for ext data */
  74. default:
  75. extension=".xxx";break;
  76. }
  77. break;
  78. case 'a': /* decide for output adjacence */
  79. if (rep==Inequality)
  80. extension=".ead"; /* adjacency file for ext data */
  81. else
  82. extension=".iad"; /* adjacency file for ine data */
  83. break;
  84. case 'i': /* decide for output incidence */
  85. if (rep==Inequality)
  86. extension=".ecd"; /* ext incidence file */
  87. else
  88. extension=".icd"; /* ine incidence file */
  89. break;
  90. case 'n': /* decide for input incidence */
  91. if (rep==Inequality)
  92. extension=".icd"; /* ine incidence file */
  93. else
  94. extension=".ecd"; /* ext incidence file */
  95. break;
  96. case 'j': /* decide for input adjacency */
  97. if (rep==Inequality)
  98. extension=".iad"; /* ine adjacency file */
  99. else
  100. extension=".ead"; /* ext adjacency file */
  101. break;
  102. case 'l':
  103. extension=".ddl";break; /* log file */
  104. case 'd':
  105. extension=".dex";break; /* decomposition output */
  106. case 'p':
  107. extension="sub.ine";break; /* preprojection sub inequality file */
  108. case 'v':
  109. extension=".solved";break; /* verify_input file */
  110. case 's':
  111. extension=".lps";break; /* LP solution file */
  112. default:
  113. extension=".xxx";break;
  114. }
  115. dotpos=-1;
  116. for (i=0; i< strlen(inputfile); i++){
  117. if (inputfile[i]=='.') dotpos=i;
  118. }
  119. if (dotpos>1) strncpy(ifilehead, inputfile, dotpos);
  120. else strcpy(ifilehead,inputfile);
  121. if (strlen(inputfile)<=0) strcpy(ifilehead,"tempcdd");
  122. strcpy(outfile,ifilehead);
  123. strcat(outfile,extension);
  124. if (strcmp(inputfile, outfile)==0) {
  125. strcpy(outfile,inputfile);
  126. strcat(outfile,extension);
  127. }
  128. /* fprintf(stderr,"outfile name = %s\n",outfile); */
  129. }
  130. NumberType GetNumberType(char *line)
  131. {
  132. NumberType nt;
  133. if (strncmp(line, "integer", 7)==0) {
  134. nt = Integer;
  135. }
  136. else if (strncmp(line, "rational", 8)==0) {
  137. nt = Rational;
  138. }
  139. else if (strncmp(line, "real", 4)==0) {
  140. nt = Real;
  141. }
  142. else {
  143. nt=Unknown;
  144. }
  145. return nt;
  146. }
  147. void ProcessCommandLine(FILE *f, Matrix M, char *line)
  148. {
  149. char newline[linelenmax];
  150. int j;
  151. double value;
  152. init(value);
  153. if (strncmp(line, "hull", 4)==0) {
  154. M.representation = Generator;
  155. }
  156. if (strncmp(line, "debug", 5)==0) {
  157. debug = true;
  158. //#ifdef GMPRATIONAL
  159. ddf_debug = ddf_TRUE;
  160. //#endif
  161. }
  162. if (strncmp(line, "partial_enum", 12)==0 ||
  163. strncmp(line, "equality", 8)==0 ||
  164. strncmp(line, "linearity", 9)==0 ) {
  165. fgets(newline,linelenmax,f);
  166. SetLinearity(M,newline);
  167. }
  168. if (strncmp(line, "maximize", 8)==0 ||
  169. strncmp(line, "minimize", 8)==0) {
  170. if (strncmp(line, "maximize", 8)==0) M.objective=LPmax;
  171. else M.objective=LPmin;
  172. for (j = 1; j <= M.colsize; j++) {
  173. if (M.numbtype==Real) {
  174. //#if !defined(GMPRATIONAL)
  175. double rvalue;
  176. fscanf(f, "%lf", &rvalue);
  177. set_d(value, rvalue);
  178. //#endif
  179. } else {
  180. fread_rational_value (f, value);
  181. }
  182. set(M.rowvec[j - 1],value);
  183. if (debug) {fprintf(stderr,"cost(%5ld) =",j); WriteNumber(stderr,value);}
  184. } /*of j*/
  185. }
  186. clear(value);
  187. }
  188. #endif
  189. public bool AppendMatrix2Poly(Polyhedra poly, Matrix M)
  190. {
  191. bool success = false;
  192. Matrix Mpoly, Mnew = null;
  193. ErrorType err;
  194. if (poly != null && poly.m >= 0 && poly.d >= 0 &&
  195. poly.d == M.colsize && M.rowsize > 0)
  196. {
  197. Mpoly = CopyInput(poly);
  198. Mnew = AppendMatrix(Mpoly, M);
  199. poly = DDMatrix2Poly(Mnew, out err);
  200. if (err == ErrorType.NoError) success = true;
  201. }
  202. return success;
  203. }
  204. public Matrix MatrixCopy(Matrix M)
  205. {
  206. Matrix Mcopy = null;
  207. int m;
  208. int d;
  209. m = M.rowsize;
  210. d = M.colsize;
  211. if (m >= 0 && d >= 0)
  212. {
  213. Mcopy = CreateMatrix(m, d);
  214. CopyAmatrix(Mcopy.matrix, M.matrix, m, d);
  215. CopyArow(Mcopy.rowvec, M.rowvec, d);
  216. M.linset.CopyTo(Mcopy.linset);
  217. Mcopy.numbtype = M.numbtype;
  218. Mcopy.representation = M.representation;
  219. Mcopy.objective = M.objective;
  220. }
  221. return Mcopy;
  222. }
  223. public Matrix CopyMatrix(Matrix M)
  224. {
  225. return MatrixCopy(M);
  226. }
  227. public Matrix MatrixNormalizedCopy(Matrix M)
  228. {
  229. Matrix Mcopy = null;
  230. int m;
  231. int d;
  232. m = M.rowsize;
  233. d = M.colsize;
  234. if (m >= 0 && d >= 0)
  235. {
  236. Mcopy = CreateMatrix(m, d);
  237. CopyNormalizedAmatrix(Mcopy.matrix, M.matrix, m, d);
  238. CopyArow(Mcopy.rowvec, M.rowvec, d);
  239. M.linset.CopyTo(Mcopy.linset);
  240. Mcopy.numbtype = M.numbtype;
  241. Mcopy.representation = M.representation;
  242. Mcopy.objective = M.objective;
  243. }
  244. return Mcopy;
  245. }
  246. public Matrix MatrixAppend(Matrix M1, Matrix M2)
  247. {
  248. Matrix M = null;
  249. int i, m, m1, m2;
  250. int j, d, d1, d2;
  251. m1 = M1.rowsize;
  252. d1 = M1.colsize;
  253. m2 = M2.rowsize;
  254. d2 = M2.colsize;
  255. m = m1 + m2;
  256. d = d1;
  257. if (d1 >= 0 && d1 == d2 && m1 >= 0 && m2 >= 0)
  258. {
  259. M = CreateMatrix(m, d);
  260. CopyAmatrix(M.matrix, M1.matrix, m1, d);
  261. CopyArow(M.rowvec, M1.rowvec, d);
  262. for (i = 0; i < m1; i++)
  263. {
  264. if (M1.linset.Get(i + 1)) M.linset.Set(i + 1, true);
  265. }
  266. for (i = 0; i < m2; i++)
  267. {
  268. for (j = 0; j < d; j++)
  269. M.matrix[m1 + i][j] = M2.matrix[i][j];
  270. /* append the second matrix */
  271. if (M2.linset.Get(i + 1)) M.linset.Set(m1 + i + 1, true);
  272. }
  273. M.numbtype = M1.numbtype;
  274. }
  275. return M;
  276. }
  277. public Matrix MatrixNormalizedSortedCopy(Matrix M, out int[] newpos) /* 094 */
  278. {
  279. /* Sort the rows of double[][] lexicographically, and return a link to this sorted copy.
  280. The vector newpos is allocated, where newpos[i] returns the new row index
  281. of the original row i (i=1,...,M.rowsize). */
  282. Matrix Mcopy = null, Mnorm = null;
  283. int m, i;
  284. int d;
  285. int[] roworder;
  286. /* if (newpos!=null) free(newpos); */
  287. m = M.rowsize;
  288. d = M.colsize;
  289. roworder = new int[m + 1];
  290. newpos = new int[m + 1];
  291. if (m >= 0 && d >= 0)
  292. {
  293. Mnorm = MatrixNormalizedCopy(M);
  294. Mcopy = CreateMatrix(m, d);
  295. for (i = 1; i <= m; i++) roworder[i] = i;
  296. RandomPermutation(roworder, m, 123);
  297. QuickSort(roworder, 1, m, Mnorm.matrix, d);
  298. PermuteCopyAmatrix(Mcopy.matrix, Mnorm.matrix, m, d, roworder);
  299. CopyArow(Mcopy.rowvec, M.rowvec, d);
  300. for (i = 1; i <= m; i++)
  301. {
  302. if (M.linset.Get(roworder[i])) Mcopy.linset.Set(i, true);
  303. newpos[roworder[i]] = i;
  304. }
  305. Mcopy.numbtype = M.numbtype;
  306. Mcopy.representation = M.representation;
  307. Mcopy.objective = M.objective;
  308. }
  309. return Mcopy;
  310. }
  311. public Matrix MatrixUniqueCopy(Matrix M, out int[] newpos)
  312. {
  313. /* Remove row duplicates, and return a link to this sorted copy.
  314. Linearity rows have priority over the other rows.
  315. It is better to call this after sorting with MatrixNormalizedSortedCopy.
  316. The vector newpos is allocated, where *newpos[i] returns the new row index
  317. of the original row i (i=1,...,M.rowsize). *newpos[i] is negative if the original
  318. row is dominated by -*newpos[i] and eliminated in the new copy.
  319. */
  320. Matrix Mcopy = null;
  321. int m, i, uniqrows;
  322. BitArray preferredrows;
  323. int d;
  324. int[] roworder;
  325. /* if (newpos!=null) free(newpos); */
  326. m = M.rowsize;
  327. d = M.colsize;
  328. preferredrows = M.linset;
  329. roworder = new int[m + 1];
  330. if (m >= 0 && d >= 0)
  331. {
  332. for (i = 1; i <= m; i++) roworder[i] = i;
  333. UniqueRows(roworder, 1, m, M.matrix, d, preferredrows, out uniqrows);
  334. Mcopy = CreateMatrix(uniqrows, d);
  335. PermutePartialCopyAmatrix(Mcopy.matrix, M.matrix, m, d, roworder, 1, m);
  336. CopyArow(Mcopy.rowvec, M.rowvec, d);
  337. for (i = 1; i <= m; i++)
  338. {
  339. if (roworder[i] > 0 && M.linset.Get(i)) Mcopy.linset.Set(roworder[i], true);
  340. }
  341. Mcopy.numbtype = M.numbtype;
  342. Mcopy.representation = M.representation;
  343. Mcopy.objective = M.objective;
  344. }
  345. newpos = roworder;
  346. return Mcopy;
  347. }
  348. public Matrix MatrixNormalizedSortedUniqueCopy(Matrix M, out int[] newpos) /* 094 */
  349. {
  350. /* Sort and remove row duplicates, and return a link to this sorted copy.
  351. Linearity rows have priority over the other rows.
  352. It is better to call this after sorting with MatrixNormalizedSortedCopy.
  353. The vector newpos is allocated, where *newpos[i] returns the new row index
  354. of the original row i (i=1,...,M.rowsize). *newpos[i] is negative if the original
  355. row is dominated by -*newpos[i] and eliminated in the new copy.
  356. */
  357. Matrix M1 = null, M2 = null;
  358. int m, i;
  359. int d;
  360. int[] newpos1 = null, newpos1r = null, newpos2 = null;
  361. /* if (newpos!=null) free(newpos); */
  362. m = M.rowsize;
  363. d = M.colsize;
  364. newpos = new int[m + 1];
  365. newpos1r = new int[m + 1];
  366. if (m >= 0 && d >= 0)
  367. {
  368. M1 = MatrixNormalizedSortedCopy(M, out newpos1);
  369. for (i = 1; i <= m; i++) newpos1r[newpos1[i]] = i; /* reverse of newpos1 */
  370. M2 = MatrixUniqueCopy(M1, out newpos2);
  371. M2.linset.SetAll(false); //emptyset
  372. for (i = 1; i <= m; i++)
  373. {
  374. if (newpos2[newpos1[i]] > 0)
  375. {
  376. //printf("newpos1[%ld]=%ld, newpos2[newpos1[%ld]]=%ld\n",i,newpos1[i], i,newpos2[newpos1[i]]);
  377. if (M.linset.Get(i)) M2.linset.Set(newpos2[newpos1[i]], true);
  378. newpos[i] = newpos2[newpos1[i]];
  379. }
  380. else
  381. {
  382. newpos[i] = -newpos1r[-newpos2[newpos1[i]]];
  383. }
  384. }
  385. }
  386. return M2;
  387. }
  388. public Matrix MatrixSortedUniqueCopy(Matrix M, out int[] newpos) /* 094 */
  389. {
  390. /* Same as MatrixNormalizedSortedUniqueCopy except that it returns a unnormalized origial data
  391. with original ordering.
  392. */
  393. Matrix M1 = null, M2 = null;
  394. int m, i, k, ii;
  395. int d;
  396. int[] newpos1 = null, newpos1r = null, newpos2 = null;
  397. /* if (newpos!=null) free(newpos); */
  398. m = M.rowsize;
  399. d = M.colsize;
  400. newpos = new int[m + 1];
  401. newpos1r = new int[m + 1];
  402. if (m >= 0 && d >= 0)
  403. {
  404. M1 = MatrixNormalizedSortedCopy(M, out newpos1);
  405. for (i = 1; i <= m; i++) newpos1r[newpos1[i]] = i; /* reverse of newpos1 */
  406. M2 = MatrixUniqueCopy(M1, out newpos2);
  407. M2.linset.SetAll(false);
  408. for (i = 1; i <= m; i++)
  409. {
  410. if (newpos2[newpos1[i]] > 0)
  411. {
  412. if (M.linset.Get(i)) M2.linset.Set(newpos2[newpos1[i]], true);
  413. newpos[i] = newpos2[newpos1[i]];
  414. }
  415. else
  416. {
  417. newpos[i] = -newpos1r[-newpos2[newpos1[i]]];
  418. }
  419. }
  420. ii = 0;
  421. M2.linset.SetAll(false);
  422. for (i = 1; i <= m; i++)
  423. {
  424. k = newpos[i];
  425. if (k > 0)
  426. {
  427. ii += 1;
  428. newpos[i] = ii;
  429. CopyArow(M2.matrix[ii - 1], M.matrix[i - 1], d);
  430. if (M.linset.Get(i)) M2.linset.Set(ii, true);
  431. }
  432. }
  433. }
  434. return M2;
  435. }
  436. public Matrix AppendMatrix(Matrix M1, Matrix M2)
  437. {
  438. return MatrixAppend(M1, M2);
  439. }
  440. public bool MatrixAppendTo(ref Matrix M1, Matrix M2)
  441. {
  442. Matrix M = null;
  443. int i, m, m1, m2;
  444. int j, d, d1, d2;
  445. bool success = false;
  446. m1 = M1.rowsize;
  447. d1 = M1.colsize;
  448. m2 = M2.rowsize;
  449. d2 = M2.colsize;
  450. m = m1 + m2;
  451. d = d1;
  452. if (d1 >= 0 && d1 == d2 && m1 >= 0 && m2 >= 0)
  453. {
  454. M = CreateMatrix(m, d);
  455. CopyAmatrix(M.matrix, M1.matrix, m1, d);
  456. CopyArow(M.rowvec, M1.rowvec, d);
  457. for (i = 0; i < m1; i++)
  458. {
  459. if (M1.linset.Get(i + 1)) M.linset.Set(i + 1, true);
  460. }
  461. for (i = 0; i < m2; i++)
  462. {
  463. for (j = 0; j < d; j++)
  464. M.matrix[m1 + i][j] = M2.matrix[i][j];
  465. /* append the second matrix */
  466. if (M2.linset.Get(i + 1)) M.linset.Set(m1 + i + 1, true);
  467. }
  468. M.numbtype = M1.numbtype;
  469. M1 = M;
  470. success = true;
  471. }
  472. return success;
  473. }
  474. public bool MatrixRowRemove(Matrix M, int r) /* 092 */
  475. {
  476. int i, m;
  477. int d;
  478. bool success = false;
  479. m = M.rowsize;
  480. d = M.colsize;
  481. if (r >= 1 && r <= m)
  482. {
  483. M.rowsize = m - 1;
  484. M.linset.Set(r, false);
  485. /* slide the row headers */
  486. for (i = r; i < m; i++)
  487. {
  488. M.matrix[i - 1] = M.matrix[i];
  489. if (M.linset.Get(i + 1))
  490. {
  491. M.linset.Set(i + 1, false);
  492. M.linset.Set(i, true);
  493. }
  494. }
  495. success = true;
  496. }
  497. return success;
  498. }
  499. public bool MatrixRowRemove2(Matrix M, int r, out int[] newpos) /* 094 */
  500. {
  501. int i, m;
  502. int d;
  503. bool success = false;
  504. int[] roworder;
  505. newpos = null;
  506. m = M.rowsize;
  507. d = M.colsize;
  508. if (r >= 1 && r <= m)
  509. {
  510. roworder = new int[m + 1];
  511. M.rowsize = m - 1;
  512. M.linset.Set(r, false);
  513. /* slide the row headers */
  514. for (i = 1; i < r; i++) roworder[i] = i;
  515. roworder[r] = 0; /* meaning it is removed */
  516. for (i = r; i < m; i++)
  517. {
  518. M.matrix[i - 1] = M.matrix[i];
  519. roworder[i + 1] = i;
  520. if (M.linset.Get(i + 1))
  521. {
  522. M.linset.Set(i + 1, false);
  523. M.linset.Set(i, true);
  524. }
  525. }
  526. success = true;
  527. }
  528. return success;
  529. }
  530. public Matrix MatrixSubmatrix(Matrix M, BitArray delset) /* 092 */
  531. {
  532. Matrix Msub = null;
  533. int i, isub = 1, m, msub;
  534. int d;
  535. m = M.rowsize;
  536. d = M.colsize;
  537. msub = m;
  538. if (m >= 0 && d >= 0)
  539. {
  540. for (i = 1; i <= m; i++)
  541. {
  542. if (delset.Get(i)) msub -= 1;
  543. }
  544. Msub = CreateMatrix(msub, d);
  545. for (i = 1; i <= m; i++)
  546. {
  547. if (!delset.Get(i))
  548. {
  549. CopyArow(Msub.matrix[isub - 1], M.matrix[i - 1], d);
  550. if (M.linset.Get(i))
  551. {
  552. Msub.linset.Set(isub, true);
  553. }
  554. isub++;
  555. }
  556. }
  557. CopyArow(Msub.rowvec, M.rowvec, d);
  558. Msub.numbtype = M.numbtype;
  559. Msub.representation = M.representation;
  560. Msub.objective = M.objective;
  561. }
  562. return Msub;
  563. }
  564. public Matrix MatrixSubmatrix2(Matrix M, BitArray delset, out int[] newpos) /* 092 */
  565. { /* returns a pointer to a new matrix which is a submatrix of M with rows in delset
  566. removed. *newpos[i] returns the position of the original row i in the new matrix.
  567. It is -1 if and only if it is deleted.
  568. */
  569. Matrix Msub = null;
  570. int i, isub = 1, m, msub;
  571. int d;
  572. int[] roworder;
  573. newpos = null;
  574. m = M.rowsize;
  575. d = M.colsize;
  576. msub = m;
  577. if (m >= 0 && d >= 0)
  578. {
  579. roworder = new int[m + 1];
  580. for (i = 1; i <= m; i++)
  581. {
  582. if (delset.Get(i)) msub -= 1;
  583. }
  584. Msub = CreateMatrix(msub, d);
  585. for (i = 1; i <= m; i++)
  586. {
  587. if (delset.Get(i))
  588. {
  589. roworder[i] = 0; /* zero means the row i is removed */
  590. }
  591. else
  592. {
  593. CopyArow(Msub.matrix[isub - 1], M.matrix[i - 1], d);
  594. if (M.linset.Get(i))
  595. {
  596. Msub.linset.Set(isub, true);
  597. }
  598. roworder[i] = isub;
  599. isub++;
  600. }
  601. }
  602. newpos = roworder;
  603. CopyArow(Msub.rowvec, M.rowvec, d);
  604. Msub.numbtype = M.numbtype;
  605. Msub.representation = M.representation;
  606. Msub.objective = M.objective;
  607. }
  608. return Msub;
  609. }
  610. public Matrix MatrixSubmatrix2L(Matrix M, BitArray delset, out int[] newpos) /* 094 */
  611. { /* This is same as MatrixSubmatrix2 except that the linearity rows will be shifted up
  612. so that they are at the top of the matrix.
  613. */
  614. Matrix Msub = null;
  615. int i, iL, iI, m, msub;
  616. int d;
  617. int[] roworder;
  618. newpos = null;
  619. m = M.rowsize;
  620. d = M.colsize;
  621. msub = m;
  622. if (m >= 0 && d >= 0)
  623. {
  624. roworder = new int[m + 1];
  625. for (i = 1; i <= m; i++)
  626. {
  627. if (delset.Get(i)) msub -= 1;
  628. }
  629. Msub = CreateMatrix(msub, d);
  630. iL = 1; iI = M.linset.Count + 1; /* starting positions */
  631. for (i = 1; i <= m; i++)
  632. {
  633. if (delset.Get(i))
  634. {
  635. roworder[i] = 0; /* zero means the row i is removed */
  636. }
  637. else
  638. {
  639. if (M.linset.Get(i))
  640. {
  641. CopyArow(Msub.matrix[iL - 1], M.matrix[i - 1], d);
  642. Msub.linset.Set(i, false);
  643. Msub.linset.Set(iL, true);
  644. roworder[i] = iL;
  645. iL += 1;
  646. }
  647. else
  648. {
  649. CopyArow(Msub.matrix[iI - 1], M.matrix[i - 1], d);
  650. roworder[i] = iI;
  651. iI += 1;
  652. }
  653. }
  654. }
  655. newpos = roworder;
  656. CopyArow(Msub.rowvec, M.rowvec, d);
  657. Msub.numbtype = M.numbtype;
  658. Msub.representation = M.representation;
  659. Msub.objective = M.objective;
  660. }
  661. return Msub;
  662. }
  663. public int MatrixRowsRemove(ref Matrix M, BitArray delset) /* 094 */
  664. {
  665. Matrix Msub = null;
  666. int success;
  667. Msub = MatrixSubmatrix(M, delset);
  668. M = Msub;
  669. success = 1;
  670. return success;
  671. }
  672. public int MatrixRowsRemove2(ref Matrix M, BitArray delset, out int[] newpos) /* 094 */
  673. {
  674. Matrix Msub = null;
  675. int success;
  676. Msub = MatrixSubmatrix2(M, delset, out newpos);
  677. M = Msub;
  678. success = 1;
  679. return success;
  680. }
  681. public int MatrixShiftupLinearity(ref Matrix M, out int[] newpos) /* 094 */
  682. {
  683. Matrix Msub = null;
  684. int success;
  685. BitArray delset;
  686. delset = new BitArray(M.rowsize); /* emptyset */
  687. Msub = MatrixSubmatrix2L(M, delset, out newpos);
  688. M = Msub;
  689. success = 1;
  690. return success;
  691. }
  692. public Polyhedra CreatePolyhedraData(int m, int d)
  693. {
  694. int i;
  695. Polyhedra poly = null;
  696. poly = new Polyhedra();
  697. poly.child = null; /* this links the homogenized cone data */
  698. poly.m = m;
  699. poly.d = d;
  700. poly.n = -1; /* the size of output is not known */
  701. poly.m_alloc = m + 2; /* the allocated row size of matrix A */
  702. poly.d_alloc = d; /* the allocated col size of matrix A */
  703. poly.numbtype = NumberType.Real;
  704. InitializeAmatrix(poly.m_alloc, poly.d_alloc, out poly.A);
  705. InitializeArow(d, out poly.c); /* cost vector */
  706. poly.representation = RepresentationType.Inequality;
  707. poly.homogeneous = false;
  708. poly.EqualityIndex = new int[m + 2];
  709. /* size increased to m+2 in 092b because it is used by the child cone,
  710. This is a bug fix suggested by Thao Dang. */
  711. /* ith component is 1 if it is equality, -1 if it is strict inequality, 0 otherwise. */
  712. for (i = 0; i <= m + 1; i++) poly.EqualityIndex[i] = 0;
  713. poly.IsEmpty = -1; /* initially set to -1, neither true nor false, meaning unknown */
  714. poly.NondegAssumed = false;
  715. poly.InitBasisAtBottom = false;
  716. poly.RestrictedEnumeration = false;
  717. poly.RelaxedEnumeration = false;
  718. poly.AincGenerated = false; /* Ainc is a set array to store the input incidence. */
  719. return poly;
  720. }
  721. public bool InitializeConeData(int m, int d, out Cone cone)
  722. {
  723. bool success = true;
  724. int j;
  725. cone = new Cone();
  726. /* INPUT: A given representation of a cone: inequality */
  727. cone.m = m;
  728. cone.d = d;
  729. cone.m_alloc = m + 2; /* allocated row size of matrix A */
  730. cone.d_alloc = d; /* allocated col size of matrix A, B and Bsave */
  731. cone.numbtype = NumberType.Real;
  732. cone.parent = null;
  733. /* CONTROL: variables to control computation */
  734. cone.Iteration = 0;
  735. cone.HalfspaceOrder = RowOrderType.LexMin;
  736. cone.ArtificialRay = null;
  737. cone.FirstRay = null;
  738. cone.LastRay = null; /* The second description: Generator */
  739. cone.PosHead = null;
  740. cone.ZeroHead = null;
  741. cone.NegHead = null;
  742. cone.PosLast = null;
  743. cone.ZeroLast = null;
  744. cone.NegLast = null;
  745. cone.RecomputeRowOrder = true;
  746. cone.PreOrderedRun = false;
  747. cone.GroundSet = new BitArray(cone.m_alloc);
  748. cone.EqualitySet = new BitArray(cone.m_alloc);
  749. cone.NonequalitySet = new BitArray(cone.m_alloc);
  750. cone.AddedHalfspaces = new BitArray(cone.m_alloc);
  751. cone.WeaklyAddedHalfspaces = new BitArray(cone.m_alloc);
  752. cone.InitialHalfspaces = new BitArray(cone.m_alloc);
  753. cone.RayCount = 0;
  754. cone.FeasibleRayCount = 0;
  755. cone.WeaklyFeasibleRayCount = 0;
  756. cone.TotalRayCount = 0;
  757. cone.ZeroRayCount = 0;
  758. cone.EdgeCount = 0;
  759. cone.TotalEdgeCount = 0;
  760. cone.count_int = 0;
  761. cone.count_int_good = 0;
  762. cone.count_int_bad = 0;
  763. cone.rseed = 1; /* random seed for random row permutation */
  764. InitializeBmatrix(cone.d_alloc, out cone.B);
  765. InitializeBmatrix(cone.d_alloc, out cone.Bsave);
  766. InitializeAmatrix(cone.m_alloc, cone.d_alloc, out cone.A);
  767. cone.Edges = new Adjacency[cone.m_alloc];
  768. cone.InitialRayIndex = new int[d + 1];
  769. cone.OrderVector = new int[cone.m_alloc + 1];
  770. cone.newcol = new int[cone.d + 1];
  771. for (j = 0; j <= cone.d; j++) cone.newcol[j] = j; /* identity map, initially */
  772. cone.LinearityDim = -2; /* -2 if it is not computed */
  773. cone.ColReduced = false;
  774. cone.d_orig = d;
  775. /* STATES: variables to represent current state. */
  776. /*cone.Error;
  777. cone.CompStatus;
  778. cone.starttime;
  779. cone.endtime;
  780. */
  781. return success;
  782. }
  783. public Cone ConeDataLoad(Polyhedra poly)
  784. {
  785. Cone cone = null;
  786. int d, j;
  787. int m, i;
  788. m = poly.m;
  789. d = poly.d;
  790. if (!(poly.homogeneous) && poly.representation == RepresentationType.Inequality)
  791. {
  792. m = poly.m + 1;
  793. }
  794. poly.m1 = m;
  795. InitializeConeData(m, d, out cone);
  796. cone.representation = poly.representation;
  797. /* Points to the original polyhedra data, and reversely */
  798. cone.parent = poly;
  799. poly.child = cone;
  800. for (i = 1; i <= poly.m; i++)
  801. for (j = 1; j <= cone.d; j++)
  802. cone.A[i - 1][j - 1] = poly.A[i - 1][j - 1];
  803. if (poly.representation == RepresentationType.Inequality && !(poly.homogeneous))
  804. {
  805. cone.A[m - 1][0] = 1.0;
  806. for (j = 2; j <= d; j++) cone.A[m - 1][j - 1] = 0.0;
  807. }
  808. return cone;
  809. }
  810. #if nemkell
  811. void SetLinearity(Matrix M, char *line)
  812. {
  813. int i=0;
  814. int eqsize,var;
  815. char *next;
  816. const char ct[]=", "; /* allows separators "," and " ". */
  817. next=strtok(line,ct);
  818. eqsize=atol(next);
  819. while (i < eqsize && (next=strtok(null,ct))!=null) {
  820. var=atol(next);
  821. set_addelem(M.linset,var); i++;
  822. }
  823. if (i!=eqsize) {
  824. fprintf(stderr,"* Warning: there are inconsistencies in linearity setting.\n");
  825. }
  826. return;
  827. }
  828. Matrix PolyFile2Matrix (FILE *f, ErrorType *Error)
  829. {
  830. Matrix M=null;
  831. int m_input,i;
  832. int d_input,j;
  833. RepresentationType rep=Inequality;
  834. double value;
  835. bool found=false, newformat=false, successful=false, linearity=false;
  836. char command[linelenmax], comsave[linelenmax], numbtype[wordlenmax];
  837. NumberType NT;
  838. //#if !defined(GMPRATIONAL)
  839. double rvalue;
  840. //#endif
  841. init(value);
  842. (*Error)=NoError;
  843. while (!found)
  844. {
  845. if (fscanf(f,"%s",command)==EOF) {
  846. (*Error)=ImproperInputFormat;
  847. goto _L99;
  848. }
  849. else {
  850. if (strncmp(command, "V-representation", 16)==0) {
  851. rep=Generator; newformat=true;
  852. }
  853. if (strncmp(command, "H-representation", 16)==0){
  854. rep=Inequality; newformat=true;
  855. }
  856. if (strncmp(command, "partial_enum", 12)==0 ||
  857. strncmp(command, "equality", 8)==0 ||
  858. strncmp(command, "linearity", 9)==0 ) {
  859. linearity=true;
  860. fgets(comsave,linelenmax,f);
  861. }
  862. if (strncmp(command, "begin", 5)==0) found=true;
  863. }
  864. }
  865. fscanf(f, "%ld %ld %s", &m_input, &d_input, numbtype);
  866. fprintf(stderr,"size = %ld x %ld\nNumber Type = %s\n", m_input, d_input, numbtype);
  867. NT=GetNumberType(numbtype);
  868. if (NT==Unknown) {
  869. (*Error)=ImproperInputFormat;
  870. goto _L99;
  871. }
  872. M=CreateMatrix(m_input, d_input);
  873. M.representation=rep;
  874. M.numbtype=NT;
  875. for (i = 1; i <= m_input; i++) {
  876. for (j = 1; j <= d_input; j++) {
  877. if (NT==Real) {
  878. //#if defined GMPRATIONAL
  879. *Error=NoRealNumberSupport;
  880. goto _L99;
  881. //#else
  882. fscanf(f, "%lf", &rvalue);
  883. set_d(value, rvalue);
  884. //#endif
  885. } else {
  886. fread_rational_value (f, value);
  887. }
  888. set(M.matrix[i-1][j - 1],value);
  889. if (debug) {fprintf(stderr,"a(%3ld,%5ld) = ",i,j); WriteNumber(stderr,value);}
  890. } /*of j*/
  891. } /*of i*/
  892. if (fscanf(f,"%s",command)==EOF) {
  893. (*Error)=ImproperInputFormat;
  894. goto _L99;
  895. }
  896. else if (strncmp(command, "end", 3)!=0) {
  897. if (debug) fprintf(stderr,"'end' missing or illegal extra data: %s\n",command);
  898. (*Error)=ImproperInputFormat;
  899. goto _L99;
  900. }
  901. successful=true;
  902. if (linearity) {
  903. SetLinearity(M,comsave);
  904. }
  905. while (!feof(f)) {
  906. fscanf(f,"%s", command);
  907. ProcessCommandLine(f, M, command);
  908. fgets(command,linelenmax,f); /* skip the CR/LF */
  909. }
  910. _L99: ;
  911. clear(value);
  912. /* if (f!=null) fclose(f); */
  913. return M;
  914. }
  915. #endif
  916. public Polyhedra DDMatrix2Poly(Matrix M, out ErrorType err)
  917. {
  918. int i;
  919. int j;
  920. Polyhedra poly = null;
  921. err = ErrorType.NoError;
  922. if (M.rowsize < 0 || M.colsize < 0)
  923. {
  924. err = ErrorType.NegativeMatrixSize;
  925. return null;
  926. }
  927. poly = CreatePolyhedraData(M.rowsize, M.colsize);
  928. poly.representation = M.representation;
  929. poly.homogeneous = true;
  930. for (i = 1; i <= M.rowsize; i++)
  931. {
  932. if (M.linset.Get(i))
  933. {
  934. poly.EqualityIndex[i] = 1;
  935. }
  936. for (j = 1; j <= M.colsize; j++)
  937. {
  938. poly.A[i - 1][j - 1] = M.matrix[i - 1][j - 1];
  939. if (j == 1 && M.matrix[i - 1][j - 1] != 0.0) poly.homogeneous = false;
  940. } /*of j*/
  941. } /*of i*/
  942. DoubleDescription(poly, out err);
  943. return poly;
  944. }
  945. public Polyhedra DDMatrix2Poly2(Matrix M, RowOrderType horder, out ErrorType err)
  946. {
  947. int i;
  948. int j;
  949. Polyhedra poly = null;
  950. err = ErrorType.NoError;
  951. if (M.rowsize < 0 || M.colsize < 0)
  952. {
  953. err = ErrorType.NegativeMatrixSize;
  954. return null;
  955. }
  956. poly = CreatePolyhedraData(M.rowsize, M.colsize);
  957. poly.representation = M.representation;
  958. poly.homogeneous = true;
  959. for (i = 1; i <= M.rowsize; i++)
  960. {
  961. if (M.linset.Get(i))
  962. {
  963. poly.EqualityIndex[i] = 1;
  964. }
  965. for (j = 1; j <= M.colsize; j++)
  966. {
  967. poly.A[i - 1][j - 1] = M.matrix[i - 1][j - 1];
  968. if (j == 1 && M.matrix[i - 1][j - 1] != 0.0) poly.homogeneous = false;
  969. } /*of j*/
  970. } /*of i*/
  971. DoubleDescription2(poly, horder, out err);
  972. return poly;
  973. }
  974. void MatrixIntegerFilter(Matrix M)
  975. { /* setting an almost integer to the integer. */
  976. int i;
  977. int j;
  978. double x;
  979. x = 0.0;
  980. for (i = 0; i < M.rowsize; i++)
  981. for (j = 0; j < M.colsize; j++)
  982. {
  983. SnapToInteger(out x, M.matrix[i][j]);
  984. M.matrix[i][j] = x;
  985. }
  986. }
  987. void CopyRay(double[] a, int d_origsize, Ray RR,
  988. RepresentationType rep, int[] reducedcol)
  989. {
  990. int j, j1;
  991. double b;
  992. b = 0.0;
  993. for (j = 1; j <= d_origsize; j++)
  994. {
  995. j1 = reducedcol[j];
  996. if (j1 > 0)
  997. {
  998. a[j - 1] = RR.RRay[j1 - 1];
  999. /* the original column j is mapped to j1, and thus
  1000. copy the corresponding component */
  1001. }
  1002. else
  1003. {
  1004. a[j - 1] = 0.0;
  1005. /* original column is redundant and removed for computation */
  1006. }
  1007. }
  1008. b = a[0];
  1009. if (rep == RepresentationType.Generator && b != 0)
  1010. {
  1011. a[0] = 1.0;
  1012. for (j = 2; j <= d_origsize; j++)
  1013. a[j - 1] /= b; /* normalization for generators */
  1014. }
  1015. }
  1016. void WriteRay(TextWriter f, int d_origsize, Ray RR, RepresentationType rep, int[] reducedcol)
  1017. {
  1018. int j;
  1019. double[] a;
  1020. InitializeArow(d_origsize + 1, out a);
  1021. CopyRay(a, d_origsize, RR, rep, reducedcol);
  1022. for (j = 0; j < d_origsize; j++) WriteNumber(f, a[j]);
  1023. f.WriteLine();
  1024. }
  1025. void WriteArow(TextWriter f, double[] a, int d)
  1026. {
  1027. int j;
  1028. for (j = 0; j < d; j++) WriteNumber(f, a[j]);
  1029. f.WriteLine();
  1030. }
  1031. void WriteAmatrix(TextWriter f, double[][] A, int rowmax, int colmax)
  1032. {
  1033. int i, j;
  1034. if (A == null)
  1035. {
  1036. f.WriteLine("WriteAmatrix: The requested matrix is empty");
  1037. return;
  1038. }
  1039. f.WriteLine("begin");
  1040. #if GMPRATIONAL
  1041. fprintf(f, " %ld %ld rational\n",rowmax, colmax);
  1042. #else
  1043. f.WriteLine(" {0} {1} real", rowmax, colmax);
  1044. #endif
  1045. for (i = 1; i <= rowmax; i++)
  1046. {
  1047. for (j = 1; j <= colmax; j++)
  1048. {
  1049. WriteNumber(f, A[i - 1][j - 1]);
  1050. }
  1051. f.WriteLine();
  1052. }
  1053. f.WriteLine("end");
  1054. }
  1055. void WriteBmatrix(TextWriter f, int d_size, double[][] B)
  1056. {
  1057. int j1, j2;
  1058. if (B == null)
  1059. {
  1060. f.WriteLine("WriteBmatrix: The requested matrix is empty");
  1061. return;
  1062. }
  1063. for (j1 = 0; j1 < d_size; j1++)
  1064. {
  1065. for (j2 = 0; j2 < d_size; j2++)
  1066. {
  1067. WriteNumber(f, B[j1][j2]);
  1068. } /*of j2*/
  1069. f.WriteLine();
  1070. } /*of j1*/
  1071. f.WriteLine();
  1072. }
  1073. #if nemkell
  1074. void WriteSetFamily(TextWriter f, SetFamily F)
  1075. {
  1076. int i;
  1077. if (F==null){
  1078. f.WriteLine("WriteSetFamily: The requested family is empty");
  1079. return;
  1080. }
  1081. f.WriteLine("begin");
  1082. f.WriteLine(" {0} {1}", F.famsize, F.setsize);
  1083. for (i=0; i<F.famsize; i++) {
  1084. f.Write(" {0} {1} : ", i+1, F.set[i].Count);
  1085. //*****set_fwrite(f, F.set[i]);
  1086. }
  1087. f.WriteLine("end");
  1088. }
  1089. void WriteSetFamilyCompressed(FILE *f, SetFamily F)
  1090. {
  1091. int i,card;
  1092. if (F==null){
  1093. fprintf(f, "WriteSetFamily: The requested family is empty\n");
  1094. goto _L99;
  1095. }
  1096. fprintf(f,"begin\n");
  1097. fprintf(f," %ld %ld\n", F.famsize, F.setsize);
  1098. for (i=0; i<F.famsize; i++) {
  1099. card=set_card(F.set[i]);
  1100. if (F.setsize - card >= card){
  1101. fprintf(f, " %ld %ld : ", i+1, card);
  1102. set_fwrite(f, F.set[i]);
  1103. } else {
  1104. fprintf(f, " %ld %ld : ", i+1, -card);
  1105. set_fwrite_compl(f, F.set[i]);
  1106. }
  1107. }
  1108. fprintf(f,"end\n");
  1109. _L99:;
  1110. }
  1111. #endif
  1112. void WriteMatrix(TextWriter f, Matrix M)
  1113. {
  1114. int i, linsize;
  1115. if (M == null)
  1116. {
  1117. f.WriteLine("WriteMatrix: The requested matrix is empty");
  1118. return;
  1119. }
  1120. switch (M.representation)
  1121. {
  1122. case RepresentationType.Inequality:
  1123. f.WriteLine("H-representation"); break;
  1124. case RepresentationType.Generator:
  1125. f.WriteLine("V-representation"); break;
  1126. case RepresentationType.Unspecified:
  1127. break;
  1128. }
  1129. linsize = M.linset.Count;
  1130. if (linsize > 0)
  1131. {
  1132. f.Write("linearity {0} ", linsize);
  1133. for (i = 1; i <= M.rowsize; i++)
  1134. if (M.linset.Get(i)) f.Write(" {0}", i);
  1135. f.WriteLine();
  1136. }
  1137. WriteAmatrix(f, M.matrix, M.rowsize, M.colsize);
  1138. if (M.objective != LPObjectiveType.LPnone)
  1139. {
  1140. if (M.objective == LPObjectiveType.LPmax)
  1141. f.WriteLine("maximize");
  1142. else
  1143. f.WriteLine("minimize");
  1144. WriteArow(f, M.rowvec, M.colsize);
  1145. }
  1146. }
  1147. void WriteLP(TextWriter f, LP lp)
  1148. {
  1149. if (lp == null)
  1150. {
  1151. f.WriteLine("WriteLP: The requested lp is empty");
  1152. return;
  1153. }
  1154. f.WriteLine("H-representation");
  1155. WriteAmatrix(f, lp.A, (lp.m) - 1, lp.d);
  1156. if (lp.objective != LPObjectiveType.LPnone)
  1157. {
  1158. if (lp.objective == LPObjectiveType.LPmax)
  1159. f.WriteLine("maximize");
  1160. else
  1161. f.WriteLine("minimize");
  1162. WriteArow(f, lp.A[lp.objrow - 1], lp.d);
  1163. }
  1164. }
  1165. void SnapToInteger(out double y, double x)
  1166. {
  1167. /* this is broken. It does nothing. */
  1168. y = x;
  1169. }
  1170. void WriteReal(TextWriter f, double x)
  1171. {
  1172. /*
  1173. int ix1,ix2,ix;
  1174. double ax;
  1175. ax=x;
  1176. ix1= (int) (Math.Abs(ax) * 10000.0 + 0.5);
  1177. ix2= (int) (Math.Abs(ax) + 0.5);
  1178. ix2= ix2*10000;
  1179. if ( ix1 == ix2) {
  1180. if (x > 0) {
  1181. ix = (int)(ax + 0.5);
  1182. } else {
  1183. ix = (int)(-ax + 0.5);
  1184. ix = -ix;
  1185. }
  1186. f.Write(" {0}", ix); //%2ld
  1187. } else
  1188. f.Write(" {0}.9E",ax); //% .9E
  1189. * */
  1190. f.Write(x);
  1191. }
  1192. void WriteNumber(TextWriter f, double x)
  1193. {
  1194. #if GMPRATIONAL
  1195. mpz_t zn,zd;
  1196. mpz_init(zn); mpz_init(zd);
  1197. mpq_canonicalize(x);
  1198. mpq_get_num(zn,x);
  1199. mpq_get_den(zd,x);
  1200. fprintf(f," ");
  1201. if (mpz_sgn(zn)==0){
  1202. fprintf(f,"0");
  1203. } else if (mpz_cmp_ui(zd,1U)==0){
  1204. mpz_out_str(f,10,zn);
  1205. } else {
  1206. mpz_out_str(f,10,zn);fprintf(f,"/");mpz_out_str(f,10,zd);
  1207. }
  1208. mpz_clear(zn); mpz_clear(zd);
  1209. #else
  1210. WriteReal(f, x);
  1211. #endif
  1212. }
  1213. #if nemkell
  1214. void WriteIncidence(FILE *f, Polyhedra poly)
  1215. {
  1216. SetFamily I;
  1217. switch (poly.representation) {
  1218. case Inequality:
  1219. fprintf(f, "ecd_file: Incidence of generators and inequalities\n");
  1220. break;
  1221. case Generator:
  1222. fprintf(f, "icd_file: Incidence of inequalities and generators\n");
  1223. break;
  1224. default:
  1225. break;
  1226. }
  1227. I=CopyIncidence(poly);
  1228. WriteSetFamilyCompressed(f,I);
  1229. FreeSetFamily(I);
  1230. }
  1231. void WriteAdjacency(FILE *f, Polyhedra poly)
  1232. {
  1233. SetFamily A;
  1234. switch (poly.representation) {
  1235. case Inequality:
  1236. fprintf(f, "ead_file: Adjacency of generators\n");
  1237. break;
  1238. case Generator:
  1239. fprintf(f, "iad_file: Adjacency of inequalities\n");
  1240. break;
  1241. default:
  1242. break;
  1243. }
  1244. A=CopyAdjacency(poly);
  1245. WriteSetFamilyCompressed(f,A);
  1246. FreeSetFamily(A);
  1247. }
  1248. #endif
  1249. void ComputeAinc(Polyhedra poly)
  1250. {
  1251. /* This generates the input incidence array poly.Ainc, and
  1252. two sets: poly.Ared, poly.Adom.
  1253. */
  1254. int k;
  1255. int i, m1;
  1256. int j;
  1257. bool redundant;
  1258. Matrix M = null;
  1259. double sum, temp;
  1260. sum = 0.0; temp = 0.0;
  1261. if (poly.AincGenerated == true) return;
  1262. M = CopyOutput(poly);
  1263. poly.n = M.rowsize;
  1264. m1 = poly.m1;
  1265. /* this number is same as poly.m, except when
  1266. poly is given by nonhomogeneous inequalty:
  1267. !(poly.homogeneous) && poly.representation==Inequality,
  1268. it is poly.m+1. See ConeDataLoad.
  1269. */
  1270. poly.Ainc = new BitArray[m1];
  1271. for (i = 1; i <= m1; i++) poly.Ainc[i - 1] = new BitArray(poly.n);
  1272. poly.Ared = new BitArray(m1);
  1273. poly.Adom = new BitArray(m1);
  1274. for (k = 1; k <= poly.n; k++)
  1275. {
  1276. for (i = 1; i <= poly.m; i++)
  1277. {
  1278. sum = 0.0;
  1279. for (j = 1; j <= poly.d; j++)
  1280. {
  1281. temp = poly.A[i - 1][j - 1] * M.matrix[k - 1][j - 1];
  1282. sum += temp;
  1283. }
  1284. if (sum == 0.0)
  1285. {
  1286. poly.Ainc[i - 1].Set(k, true);
  1287. }
  1288. }
  1289. if (!(poly.homogeneous) && poly.representation == RepresentationType.Inequality)
  1290. {
  1291. if (M.matrix[k - 1][0] == 0.0)
  1292. {
  1293. poly.Ainc[m1 - 1].Set(k, true); /* added infinity inequality (1,0,0,...,0) */
  1294. }
  1295. }
  1296. }
  1297. for (i = 1; i <= m1; i++)
  1298. {
  1299. if (poly.Ainc[i - 1].Count == M.rowsize)
  1300. {
  1301. poly.Adom.Set(i, true);
  1302. }
  1303. }
  1304. for (i = m1; i >= 1; i--)
  1305. {
  1306. if (poly.Ainc[i - 1].Count == 0)
  1307. {
  1308. redundant = true;
  1309. poly.Ared.Set(i, true);
  1310. }
  1311. else
  1312. {
  1313. redundant = false;
  1314. for (k = 1; k <= m1; k++)
  1315. {
  1316. if (k != i && !poly.Ared.Get(k) && !poly.Adom.Get(k) &&
  1317. poly.Ainc[i - 1].IsSubsetOf(poly.Ainc[k - 1]))
  1318. { //****
  1319. if (!redundant)
  1320. {
  1321. redundant = true;
  1322. }
  1323. poly.Ared.Set(i, true);
  1324. }
  1325. }
  1326. }
  1327. }
  1328. poly.AincGenerated = true;
  1329. }
  1330. bool InputAdjacentQ(Polyhedra poly,
  1331. int i1, int i2)
  1332. /* Before calling this function, RedundantSet must be
  1333. a set of row indices whose removal results in a minimal
  1334. nonredundant system to represent the input polyhedron,
  1335. DominantSet must be the set of row indices which are
  1336. active at every extreme points/rays.
  1337. */
  1338. {
  1339. bool adj = true;
  1340. int i;
  1341. BitArray common;
  1342. int lastn = 0;
  1343. if (poly.AincGenerated == false) ComputeAinc(poly);
  1344. if (lastn != poly.n)
  1345. {
  1346. //if (lastn >0) set_free(common);
  1347. common = new BitArray(poly.n);
  1348. lastn = poly.n;
  1349. }
  1350. if (poly.Ared.Get(i1) || poly.Ared.Get(i2))
  1351. {
  1352. adj = false;
  1353. return adj;
  1354. }
  1355. if (poly.Adom.Get(i1) || poly.Adom.Get(i2))
  1356. {
  1357. // dominant inequality is considered adjacencent to all others.
  1358. adj = true;
  1359. return adj;
  1360. }
  1361. common = poly.Ainc[i1 - 1].And(poly.Ainc[i2 - 1]); //**** intersection
  1362. i = 0;
  1363. while (i < poly.m1 && adj == true)
  1364. {
  1365. i++;
  1366. if (i != i1 && i != i2 && !poly.Ared.Get(i) &&
  1367. !poly.Adom.Get(i) && common.IsSubsetOf(poly.Ainc[i - 1]))
  1368. {
  1369. adj = false;
  1370. }
  1371. }
  1372. return adj;
  1373. }
  1374. #if nemkell
  1375. void WriteInputIncidence(FILE *f, Polyhedra poly)
  1376. {
  1377. SetFamily I;
  1378. if (poly.AincGenerated==false) ComputeAinc(poly);
  1379. switch (poly.representation) {
  1380. case Inequality:
  1381. fprintf(f,"icd_file: Incidence of inequalities and generators\n");
  1382. break;
  1383. case Generator:
  1384. fprintf(f,"ecd_file: Incidence of generators and inequalities\n");
  1385. break;
  1386. default:
  1387. break;
  1388. }
  1389. I=CopyInputIncidence(poly);
  1390. WriteSetFamilyCompressed(f,I);
  1391. FreeSetFamily(I);
  1392. }
  1393. void WriteInputAdjacency(FILE *f, Polyhedra poly)
  1394. {
  1395. SetFamily A;
  1396. if (poly.AincGenerated==false){
  1397. ComputeAinc(poly);
  1398. }
  1399. switch (poly.representation) {
  1400. case Inequality:
  1401. fprintf(f, "iad_file: Adjacency of inequalities\n");
  1402. break;
  1403. case Generator:
  1404. fprintf(f, "ead_file: Adjacency of generators\n");
  1405. break;
  1406. default:
  1407. break;
  1408. }
  1409. A=CopyInputAdjacency(poly);
  1410. WriteSetFamilyCompressed(f,A);
  1411. FreeSetFamily(A);
  1412. }
  1413. void WriteProgramDescription(FILE *f)
  1414. {
  1415. fprintf(f, "* cddlib: a double description library:%s\n", DDVERSION);
  1416. fprintf(f, "* compiled for %s arithmetic.\n", ARITHMETIC);
  1417. fprintf(f,"* %s\n",CO

Large files files are truncated, but you can click here to view the full file