1 /****************************************************************************************/ 2 /* */ 3 /* vinci.h */ 4 /* */ 5 /****************************************************************************************/ 6 /* */ 7 /* Authors: Benno Bueeler (bueeler@ifor.math.ethz.ch) */ 8 /* and */ 9 /* Andreas Enge (enge@ifor.math.ethz.ch) */ 10 /* Institute for Operations Research */ 11 /* Swiss Federal Institute of Technology Zurich */ 12 /* Switzerland */ 13 /* */ 14 /* Last Changes: May 19, 2003 */ 15 /* */ 16 /****************************************************************************************/ 17 /* */ 18 /* global constants, types, variable and function declarations */ 19 /* */ 20 /****************************************************************************************/ 21 22 #define COPYRIGHT1 "(C) Benno Bueeler and Andreas Enge, {bueeler,enge}@ifor.math.ethz.ch" 23 #define COPYRIGHT2 " Developed in the ETHZ-EPFL project on Optimisation and Computational" 24 #define COPYRIGHT3 " Geometry with Komei Fukuda" 25 #define VERSION "1.0.5" 26 #define VERSION_DATE "July 2003" 27 28 #define T01 "The package computes the volume of a polytope whose vertices, defining" 29 #define T02 "hyperplanes and/or incidences of vertices and facettes are stored in files" 30 #define T03 "following Avis' and Fukuda's polytope format. The vertices are supposed to" 31 #define T04 "be in a file with extension '.ine', the hyperplanes in '.ext' and the" 32 #define T05 "incidences in '.icd' (see the sample files 'square.ext', 'square.icd' and" 33 #define T06 "'square.ine').\n" 34 #define T07 "Its basic call is 'vinci file' where 'file' stands for the polyhedron file" 35 #define T08 "name without extension, e. g. 'vinci square'. In this case the existing files" 36 #define T09 "and installed additional packages are analysed and according to the problem" 37 #define T10 "type an appropriate volume computation method is chosen automatically." 38 #define T11 "" 39 #define T12 "The following additional parameters are possible:" 40 #define T13 "-m followed by the label of a specific volume computation method; to see a" 41 #define T14 " list of the implemented methods, use 'vinci -m' without any label." 42 #define T15 "-s directly followed by a natural integer. The value determines for how many" 43 #define T16 " recursion levels intermediate results are stored. A higher value speeds" 44 #define T17 " up certain methods considerably while needing more storage space." 45 #define T18 "-r directly followed by an integer. The value sets the random seed used for" 46 #define T19 " determining the objective function for Lawrence's formula." 47 #define T20 "\nFor more information please consult the manual." 48 49 /****************************************************************************************/ 50 51 #include <stdio.h> 52 #include <stdlib.h> 53 #include <string.h> 54 #include <math.h> 55 #include <sys/times.h> 56 57 /****************************************************************************************/ 58 /* The following constants may be changed by a user */ 59 /****************************************************************************************/ 60 61 #define LRS_EXEC "lrs" 62 /* location of external programme with path */ 63 64 #define PIVOTING 1 65 #define MIN_PIVOT 0.5 66 #define PIVOTING_LASS 0 67 #define MIN_PIVOT_LASS 0.1 68 /* for choosing a pivoting strategy whenever this is needed, e.g. for determinant */ 69 /* computation. The last two constants are valid for Lasserre's method, the first */ 70 /* two for any other method. PIVOTING and PIVOTING_LASS can take the following */ 71 /* values: */ 72 /* 0: The first row with entry bigger than MIN_PIVOT (in its absolute value) is */ 73 /* chosen. */ 74 /* 1: partial pivoting; the row with maximal entry (absolute value) is chosen; */ 75 /* MIN_PIVOT is ignored. */ 76 /* 2: total pivoting; the row and column with maximal entry (absolute value) are */ 77 /* chosen, MIN_PIVOT is ignored. If total pivoting is not possible (e. g. for */ 78 /* solving linear equations) partial pivoting is performed. */ 79 80 #define DEFAULT_STORAGE 20 81 /* For so many recursion levels intermediate volumes are stored, starting from the */ 82 /* level with two planes fixed. The value is only active if no option -s is speci- */ 83 /* fied; it may be set to 0, for instance, by using the option -s0. */ 84 85 #define STATISTICS 86 /* If STATISTICS is defined, during volume computation, some statistical variables */ 87 /* like the number of simplices and their volume distribution are withheld. */ 88 89 #define NO_RATIONAL 90 /* If RATIONAL is defined, exact computation will be used in some places. In this */ 91 /* case, some real variables will be defined as rational. This option only works if */ 92 /* the g++-compiler is used, and its exact consequences are unknown... */ 93 94 /****************************************************************************************/ 95 /* The following constants should not be changed by a user */ 96 /****************************************************************************************/ 97 98 #define EPSILON 1e-10 99 #define INCIDENCE_EPSILON 1e-6 100 /* Our usual constant of 1e-10 results in wrong incidence structures for rv_10_14 */ 101 /* and rv_8_30. Using the cdd constant of 1e-6 everywhere produces a wrong result */ 102 /* with hot for ccp6. Therefore, we define a separate constant for the incidence */ 103 /* computation. */ 104 #define ARRAYSIZESTEP 5 105 /* the dynamic arrays are increased in steps of this size */ 106 107 #define CLOCKS_PER_SEC 100 108 /* defines the number of clock ticks per second returned by */ 109 /* times() */ 110 111 #define FALSE 0 112 #define TRUE 1 113 #define EOLN 10 /* end of line character */ 114 115 #define NONE 0 /* constants for data types in files */ 116 #define INTEGER_T 1 117 #define REAL_T 2 118 #define RATIONAL_T 3 119 120 #define RCH 1 /* constants for the volume computation methods */ 121 #define HOT 4 122 #define LAWD 9 123 #define LAWND 10 124 #define RLASS 11 125 #define LRS 12 126 127 #define KEY_VERTICES 2 /* constants for the key type actually used in the balanced */ 128 #define KEY_PLANES_VAR 3 /* tree routines */ 129 130 #ifdef STATISTICS 131 #define STAT_SMALLEST_EXP -200 132 #define STAT_BIGGEST_EXP 200 133 #endif 134 /* The simplices in the size classes from 1e(STAT_SMALLEST_EXP) to */ 135 /* 1e(STAT_BIGGEST_EXP) are counted separately for each class. Simplices too small */ 136 /* or too big are summarized in one variable. */ 137 138 #ifdef RATIONAL 139 #include <Rational.h> 140 #endif 141 142 /****************************************************************************************/ 143 /* type definitions */ 144 /****************************************************************************************/ 145 146 typedef double real; 147 typedef unsigned char boolean; 148 typedef unsigned char T_LassInt; 149 150 #ifdef RATIONAL 151 typedef Rational rational; 152 #else 153 typedef real rational; 154 #endif 155 156 struct T_Vertex 157 {real *coords; 158 /* the coordinates of a vertex */ 159 int no; 160 }; /* type of a vertex */ 161 typedef struct T_Vertex T_Vertex; 162 163 /* The sets of vertices are implemented as ordered lists */ 164 struct T_VertexSet 165 {int maxel; /* The list of vertices loe may contain elements 0 to maxel */ 166 int lastel; /* Effectively, loe contains elements from 0 to lastel; */ 167 /* an empty set is indicated by lastel == -1. */ 168 T_Vertex **loe; /* The elements of the set are pointers to vertices which are */ 169 /* stored in loe, in ascending order following the numbering */ 170 /* of the vertices */ 171 }; 172 typedef struct T_VertexSet T_VertexSet; 173 174 175 /* The sets of VertexSets are implemented as linked chain because only simple */ 176 /* operations are needed on them. */ 177 struct T_VertexSuperset 178 {struct T_VertexSuperset *next; 179 struct T_VertexSet content; 180 }; 181 typedef struct T_VertexSuperset T_VertexSuperset; 182 183 184 /* Types for storing face volumes in balanced trees. For storing and retrieving a key */ 185 /* is needed for every face; all possibilities for keys are defined by the union T_Key. */ 186 /* Only one of the keys is active at a time. T_Tree defines the tree itself. */ 187 union T_Key 188 { 189 struct {T_VertexSet set; 190 int d; 191 } vertices; 192 /* The vertices contained in the face and the dimension for which the volume */ 193 /* is stored. */ 194 struct {T_LassInt *hyperplanes, *variables;} hypervar; 195 /* A list of fixed constraints and variables onto which the face has been */ 196 /* projected, used for Lasserre's formula */ 197 }; 198 typedef union T_Key T_Key; 199 200 struct T_Tree 201 {struct T_Tree *tree_l, *tree_r; /* the left and right subtrees */ 202 int tree_b; 203 T_Key key; 204 rational vol; /* the stored volume */ 205 }; 206 typedef struct T_Tree T_Tree; 207 208 /****************************************************************************************/ 209 /* global variable declarations */ 210 /****************************************************************************************/ 211 212 extern int G_d; 213 extern int G_m; 214 extern int G_n; 215 extern real **G_Hyperplanes; 216 /* The first d components are the coordinates of a normal vector, the d+first com- */ 217 /* ponent is the right hand side. */ 218 extern boolean **G_Incidence; 219 /* Incidence structure of vertices and hyperplanes. The rows are indexed with the */ 220 /* vertex numbers, the columns with the hyperplane numbers, and the value is TRUE */ 221 /* if the incidence relation is fulfilled. This allows to easily renumber the */ 222 /* vertices while keeping the incidence structure up to date by simply exchanging */ 223 /* rows. */ 224 extern T_VertexSet G_Vertices; 225 /* set of the vertices of the polytope */ 226 extern int G_Storage; 227 /* see the annotations for DEFAULT_STORAGE */ 228 extern int G_RandomSeed; 229 extern rational G_Minus1; 230 /* dummy to have a pointer, namely &G_Minus1, to the value -1 */ 231 232 #ifdef STATISTICS 233 /* some variables for statistical purpose */ 234 extern unsigned int Stat_Count; 235 /* for counting the number of partial volumes computed */ 236 extern real Stat_Smallest, Stat_Biggest; 237 /* the volumes of the smallest and the biggest simplex */ 238 extern unsigned int Stat_CountPos [STAT_BIGGEST_EXP - STAT_SMALLEST_EXP + 3]; 239 /* Stat_CountPos [i] counts the number of positive partial volumes (usually */ 240 /* simplex volumes except for Lawrence's formula) which are contained between */ 241 /* 1e(i - 1 + STAT_SMALLEST_EXP) and 1e(i + STAT_SMALLEST_EXP). Smaller volumes */ 242 /* are counted in Stat_CountPos [0], bigger ones in the last entry of */ 243 /* Stat_CountPos. */ 244 extern unsigned int Stat_CountNeg [STAT_BIGGEST_EXP - STAT_SMALLEST_EXP + 3]; 245 /* the same for the absolute values of "negative volumes" in Lawrence's formula */ 246 extern unsigned int *Stat_CountStored, *Stat_CountRetrieved; 247 /* counts the number of volumes stored in and retrieved from the tree */ 248 extern unsigned int Stat_CountShifts; 249 /* counts the number of shifts performed in Lasserre's method */ 250 extern long int Stat_ActualMem; 251 /* the memory actually used on the heap */ 252 extern long int Stat_MaxMem; 253 /* the maximal heap memory used during execution */ 254 /* Both values do not comprise the memory reserve G_MemRes. */ 255 #endif 256 257 /****************************************************************************************/ 258 /* functions and procedures from 'vinci_memory' */ 259 /****************************************************************************************/ 260 261 void *my_malloc (long int size); 262 void *my_realloc (void *pointer, long int new_size, long int size_diff); 263 void my_free (void *pointer, long int size); 264 T_Vertex *create_vertex (); 265 void free_vertex (T_Vertex *v); 266 void create_hyperplanes (); 267 void free_hyperplanes (); 268 void create_incidence (); 269 void free_incidence (); 270 T_VertexSet *create_faces (); 271 void free_faces (T_VertexSet *face); 272 rational ***create_basis (); 273 void free_basis (rational *** basis); 274 rational *create_fact (); 275 void free_fact (rational *fact); 276 rational *create_vector (); 277 void free_vector (rational *v); 278 int *create_int_vector (int n); 279 void free_int_vector (int *v, int n); 280 rational **create_matrix (int m, int n); 281 void redim_matrix (rational ***A, int m_alt, int m_neu, int n); 282 void free_matrix (rational **A, int m, int n); 283 void create_key (T_Key *key, int key_choice); 284 void free_key (T_Key key, int key_choice); 285 #ifdef STATISTICS 286 void init_statistics (); 287 void update_statistics (rational volume); 288 void free_statistics (); 289 #endif 290 void tree_out (T_Tree **ppr , boolean *pi_balance, T_Key key, rational **volume, 291 T_Key **keyfound, int key_choice); 292 void add_hypervar (T_LassInt hyperplane, T_LassInt variable, T_Key *key); 293 void delete_hypervar (T_LassInt hyperplane, T_LassInt variable, T_Key *key); 294 295 /****************************************************************************************/ 296 /* functions and procedures from 'vinci_set' */ 297 /****************************************************************************************/ 298 299 T_VertexSet create_empty_set (void); 300 void free_set (T_VertexSet s); 301 void free_set_and_vertices (T_VertexSet s); 302 void clear_set (T_VertexSet *s); 303 T_VertexSet duplicate_set (T_VertexSet s); 304 void copy_set (T_VertexSet s1, T_VertexSet *s2); 305 void renumber_vertices (); 306 rational normalise_vertices (); 307 void print_set (FILE *f, T_VertexSet s); 308 boolean empty (T_VertexSet s); 309 boolean is_in_hyperplane (T_Vertex *e, int j); 310 boolean are_equal_sets (T_VertexSet s1, T_VertexSet s2); 311 void add_element (T_VertexSet *s, T_Vertex *e); 312 boolean delete_element (T_VertexSet *s, T_Vertex *e); 313 void intersect_with_hyperplane (T_VertexSet s, int j, T_VertexSet *inter); 314 T_VertexSuperset *create_empty_superset (void); 315 void free_superset (T_VertexSuperset **S); 316 void print_superset (FILE *f, T_VertexSuperset *S); 317 boolean is_in_superset (T_VertexSet s, T_VertexSuperset *S); 318 void add_superelement (T_VertexSuperset **S, T_VertexSet s); 319 320 /****************************************************************************************/ 321 /* functions and procedures from 'vinci_computation' */ 322 /****************************************************************************************/ 323 324 rational factorial (int n); 325 rational det_and_invert (rational **A, int rows, int columns, boolean verbose); 326 void simplex_volume (T_VertexSet S, rational *volume, boolean verbose); 327 rational add_orthonormal (int d, T_VertexSet face, rational **H, T_Vertex *vertex); 328 rational orthonormal (int d, T_VertexSet face, rational **H); 329 330 /****************************************************************************************/ 331 /* functions and procedures from 'vinci_volume' */ 332 /****************************************************************************************/ 333 334 void volume_ch_file (rational *volume, char *vertexfile, char *planesfile); 335 void volume_ortho_file (rational *volume, char *vertexfile, char *planesfile); 336 void volume_lawrence_file (rational *volume, char *vertexfile, char *planesfile); 337 void volume_lawrence_lrs_file (rational *volume, char *planesfile); 338 void volume_lrs_file (rational *volume, char *rational_volume, char *vertexfile); 339 340 /****************************************************************************************/ 341 /* functions and procedures from 'vinci_lass' */ 342 /****************************************************************************************/ 343 344 void volume_lasserre_file (rational *volume, char *planesfile); 345 346 /****************************************************************************************/ 347 /* functions and procedures from 'vinci_screen' */ 348 /****************************************************************************************/ 349 350 void print_hyperplanes (FILE *f); 351 void print_coords (FILE *f, T_Vertex *v); 352 void print_incidence (FILE *f); 353 void print_matrix (FILE *f, rational **A, int m, int n); 354 #ifdef STATISTICS 355 void print_statistics (FILE *f, int method); 356 #endif 357 358 /****************************************************************************************/ 359 /* functions and procedures from 'vinci_file' */ 360 /****************************************************************************************/ 361 362 FILE * open_read (char *filename); 363 int determine_data_type (char *data_type); 364 void sread_rational_value (char *s, rational *value); 365 void read_vertices (char *filename); 366 void read_hyperplanes (char *filename); 367 void compute_incidence (); 368 369 /****************************************************************************************/ 370