[BACK]Return to drivers_for_dynamic_lifting.adb CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Dynlift

Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Dynlift/drivers_for_dynamic_lifting.adb, Revision 1.1.1.1

1.1       maekawa     1: with integer_io;                         use integer_io;
                      2: with Communications_with_User;           use Communications_with_User;
                      3: with Timing_Package;                     use Timing_Package;
                      4: with Numbers_io;                         use Numbers_io;
                      5: with Standard_Integer_Vectors;           use Standard_Integer_Vectors;
                      6: with Standard_Integer_Vectors_io;        use Standard_Integer_Vectors_io;
                      7: with Standard_Floating_Vectors;
                      8: with Standard_Complex_Poly_Systems_io;   use Standard_Complex_Poly_Systems_io;
                      9: with Lists_of_Integer_Vectors;           use Lists_of_Integer_Vectors;
                     10: with Lists_of_Integer_Vectors_io;        use Lists_of_Integer_Vectors_io;
                     11: with Lists_of_Vectors_Utilities;         use Lists_of_Vectors_Utilities;
                     12: with Arrays_of_Integer_Vector_Lists;     use Arrays_of_Integer_Vector_Lists;
                     13: with Arrays_of_Integer_Vector_Lists_io;  use Arrays_of_Integer_Vector_Lists_io;
                     14: with Standard_Complex_Solutions_io;      use Standard_Complex_Solutions_io;
                     15: with Power_Lists;                        use Power_Lists;
                     16: with Drivers_for_Vertex_Points;          use Drivers_for_Vertex_Points;
                     17: with Integer_Lifting_Utilities;          use Integer_Lifting_Utilities;
                     18: with Integer_Mixed_Subdivisions;         use Integer_Mixed_Subdivisions;
                     19: with Integer_Mixed_Subdivisions_io;      use Integer_Mixed_Subdivisions_io;
                     20: with Mixed_Volume_Computation;           use Mixed_Volume_Computation;
                     21: with Simplices,Triangulations;           use Simplices,Triangulations;
                     22: with Triangulations_io;                  use Triangulations_io;
                     23: with Dynamic_Triangulations;             use Dynamic_Triangulations;
                     24: with Cayley_Trick;                       use Cayley_Trick;
                     25: with Driver_for_Minkowski_Polynomials;
                     26: with Flatten_Mixed_Subdivisions;         use Flatten_Mixed_Subdivisions;
                     27: with Triangulations_and_Subdivisions;    use Triangulations_and_Subdivisions;
                     28: with Dynamic_Mixed_Subdivisions;         use Dynamic_Mixed_Subdivisions;
                     29: with Dynamic_Polyhedral_Continuation;    use Dynamic_Polyhedral_Continuation;
                     30: with Driver_for_Polyhedral_Continuation;
                     31: with Pruning_Statistics;
                     32:
                     33: package body Drivers_for_Dynamic_Lifting is
                     34:
                     35:   procedure Dynamic_Lifting_Info is
                     36:
                     37:     i : array(1..6) of string(1..65);
                     38:
                     39:   begin
                     40:     i(1):="  Dynamic  lifting  can  be  used  to   compute   mixed   volumes";
                     41:     i(2):="incrementally,  i.e.:  by  adding  the  points  repeatedly to the";
                     42:     i(3):="already constructed subdivision.  This method  works  efficiently";
                     43:     i(4):="when  all  Newton polytopes are (almost) equal.  The Cayley trick";
                     44:     i(5):="is implemented by means of dynamic lifting.  This trick  computes";
                     45:     i(6):="all cells in a mixed subdivision.                                ";
                     46:     for k in i'range loop
                     47:       put_line(i(k));
                     48:     end loop;
                     49:   end Dynamic_Lifting_Info;
                     50:
                     51:   procedure Write_Type_of_Mixture
                     52:               ( file : in file_type; mix,per : in Vector ) is
                     53:
                     54:   -- DESCRIPTION :
                     55:   --   Writes the information about the type of mixture on file
                     56:   --   and the permutations in the support.
                     57:
                     58:   begin
                     59:     new_line(file);
                     60:     put(file,"TYPE OF MIXTURE : ");
                     61:     put(file,"#supports : "); put(file,mix'last,1);
                     62:     put(file,"  occurrences : "); put(file,mix);
                     63:     new_line(file);
                     64:     put(file,"  permutation : "); put(file,per);
                     65:     new_line(file);
                     66:   end Write_Type_of_Mixture;
                     67:
                     68:   procedure Driver_for_Dynamic_Mixed_Volume_Computation
                     69:                 ( file : in file_type; p : in Poly_Sys; byebye : in boolean;
                     70:                   q : out Poly_Sys; qsols : out Solution_List;
                     71:                   mv : out natural ) is
                     72:
                     73:     welcome : constant string := "Mixed-Volume Computation by Dynamic Lifting";
                     74:
                     75:   -- GLOBAL VARIABLES :
                     76:
                     77:     supports : Array_of_Lists(p'range);
                     78:     n : natural := p'last;
                     79:     timer : timing_widget;
                     80:     r,max : natural;
                     81:     mix,perms : Link_to_Vector;
                     82:     ans : character;
                     83:     permp,qq : Poly_Sys(p'range);
                     84:     qqsols : Solution_List;
                     85:     subfile,solsft,qft : file_type;
                     86:     vol : natural := 0;
                     87:     mixsub : Mixed_Subdivision;
                     88:
                     89:   -- GLOBAL SWITCHES :
                     90:
                     91:     verpts     : boolean;  -- if the set of vertex points is computed
                     92:     order      : boolean;  -- process points in fixed instead of random order
                     93:     inter      : boolean;  -- if interior points are possible
                     94:     conmv      : boolean;  -- if checks on zero contributions have to be made
                     95:     caytrick   : boolean;  -- if the Cayley trick has to be applied
                     96:     reportnew  : boolean;  -- if the new cells have to be reported
                     97:     reportflat : boolean;  -- if before flattening, reporting has to be done
                     98:     subonfile  : boolean;  -- put the subdivision on separate file
                     99:     tosolve    : boolean;  -- if the system needs to be solved
                    100:     contrep    : boolean;  -- if intermediate output during continuation
                    101:     ranstart   : boolean;  -- if random coefficient start system
                    102:
                    103:     minkpoly   : natural;  -- 0 : no; 1 : only poly, > 1 : all subdivisions
                    104:
                    105:   -- EXTRACT THE ADDITIONAL POINTS :
                    106:
                    107:     function Is_In_Lifted ( pt : Link_to_Vector; lifted : List )
                    108:                           return boolean is
                    109:
                    110:     -- DESCRIPTION :
                    111:     --   Returns true if the point is in the lifted list.
                    112:
                    113:       tmp : List := lifted;
                    114:       lpt : Link_to_Vector;
                    115:
                    116:     begin
                    117:       while not Is_Null(tmp) loop
                    118:         lpt := Head_Of(tmp);
                    119:         if pt(pt'range) = lpt(pt'range)
                    120:          then return true;
                    121:          else tmp := Tail_Of(tmp);
                    122:         end if;
                    123:       end loop;
                    124:       return false;
                    125:     end Is_In_Lifted;
                    126:
                    127:     function Difference ( supp,liftsupp : in List ) return List is
                    128:
                    129:       res,res_last : List;
                    130:       tmp : List := supp;
                    131:       pt : Link_to_Vector;
                    132:
                    133:     begin
                    134:       tmp := supp;
                    135:       while not Is_Null(tmp) loop
                    136:         pt := Head_Of(tmp);
                    137:         if not Is_In_Lifted(pt,liftsupp)
                    138:          then Append(res,res_last,pt.all);
                    139:         end if;
                    140:         tmp := Tail_Of(tmp);
                    141:       end loop;
                    142:       return res;
                    143:     end Difference;
                    144:
                    145:     function Difference ( supp,liftsupp : in Array_of_Lists )
                    146:                         return Array_of_Lists is
                    147:
                    148:     -- DESCRIPTION :
                    149:     --   Returns a tuple of point lists, made of points in supp
                    150:     --   that do not belong to the corresponding lifted supports.
                    151:
                    152:       res : Array_of_Lists(supp'range);
                    153:
                    154:     begin
                    155:       for i in supp'range loop
                    156:         res(i) := Difference(supp(i),liftsupp(i));
                    157:       end loop;
                    158:       return res;
                    159:     end Difference;
                    160:
                    161:   -- DETERMINING THE ORDER OF PROCESSING THE POINTS :
                    162:
                    163:     function Determine_Order ( l : List ) return List is
                    164:
                    165:     -- DESCRIPTION :
                    166:     --   Interactive ordering of the points in the list.
                    167:     --   This function displays all points and asks the user for a position.
                    168:
                    169:       len : constant natural := Length_Of(l);
                    170:       pos : vector(1..len);
                    171:       res : List;
                    172:
                    173:       function Read_New_Positions
                    174:                   ( l : List; length : natural ) return vector is
                    175:
                    176:       -- DESCRIPTION :
                    177:       --   Lists all points in the lists and prompts for a new position.
                    178:       --   Returns the position vector.
                    179:
                    180:         newpos : vector(1..length);
                    181:         tmp : List := l;
                    182:         cnt : natural := 0;
                    183:
                    184:       begin
                    185:         put("There are "); put(length,1); put_line(" points to order.");
                    186:         put_line("Give for each separate point its new position :");
                    187:         while not Is_Null(tmp) loop
                    188:           cnt := cnt + 1;
                    189:           loop
                    190:             put(Head_Of(tmp)); put(" : ");
                    191:             Read_Natural(newpos(cnt));
                    192:             exit when (newpos(cnt) >= 1) and (newpos(cnt) <= length);
                    193:             put("New position out of range 1.."); put(length,1);
                    194:             put_line(".  Please try again.");
                    195:           end loop;
                    196:           tmp := Tail_Of(tmp);
                    197:         end loop;
                    198:         return newpos;
                    199:       end Read_New_Positions;
                    200:
                    201:       function Get ( l : List; pos : natural ) return Link_to_Vector is
                    202:
                    203:       -- DESCRIPTION :
                    204:       --   Returns the point on the indicated position in the list l.
                    205:
                    206:         tmp : List := l;
                    207:         res : Link_to_Vector;
                    208:
                    209:       begin
                    210:         if not Is_Null(l)
                    211:          then for i in 1..(pos-1) loop
                    212:                 tmp := Tail_Of(tmp);
                    213:                 exit when Is_Null(tmp);
                    214:               end loop;
                    215:               if not Is_Null(tmp)
                    216:                then res := Head_Of(tmp);
                    217:               end if;
                    218:         end if;
                    219:         return res;
                    220:       end Get;
                    221:
                    222:       function Sort ( l : in List; pos : in vector ) return List is
                    223:
                    224:       -- DESCRIPTION :
                    225:       --   Sorts the given list according to the given position vector:
                    226:       --   pos(i) determines the new position of the ith point in the list.
                    227:       --   If the returning list is empty, then the position vector was
                    228:       --   not a permutation.
                    229:
                    230:         empty,res,res_last : List;
                    231:         index : natural;
                    232:
                    233:       begin
                    234:         for i in pos'range loop   -- search index : pos(index) = i
                    235:         index := 0;
                    236:           for j in pos'range loop
                    237:             if pos(j) = i
                    238:              then index := j;
                    239:             end if;
                    240:             exit when (index /= 0);
                    241:           end loop;
                    242:           exit when (index = 0);
                    243:           Append(res,res_last,get(l,index));  -- append the vector
                    244:         end loop;
                    245:         if index = 0
                    246:          then return empty;
                    247:          else return res;
                    248:         end if;
                    249:       end Sort;
                    250:
                    251:     begin
                    252:       if Is_Null(l)
                    253:        then return l;
                    254:        else loop
                    255:               pos := Read_New_Positions(l,len);
                    256:               res := Sort(l,pos);
                    257:               exit when not Is_Null(res);
                    258:               put_line("The given position vector was not a permutation.");
                    259:               put_line("Please try again...");
                    260:             end loop;
                    261:             return res;
                    262:       end if;
                    263:     end Determine_Order;
                    264:
                    265:     procedure Determine_Processing_Order
                    266:                  ( supports : in out Array_of_Lists; mix : in Link_to_Vector;
                    267:                    fixed : out boolean ) is
                    268:
                    269:       choice : character;
                    270:       cnt : natural;
                    271:
                    272:     begin
                    273:       new_line;
                    274:       put_line("MENU for the Order of the points to add : ");
                    275:       put_line("  1. fixed order, given by the monomial ordering");
                    276:       put_line("  2. random order, generated by the algorithm");
                    277:       put_line("  3. interactively defined by you");
                    278:       put("Type 1,2, or 3 : "); Ask_Alternative(choice,"123");
                    279:       case choice is
                    280:         when '1'    => fixed := true;
                    281:         when '2'    => fixed := false;
                    282:         when others => fixed := true;
                    283:                        cnt := supports'first;
                    284:                        for i in mix'range loop
                    285:                          supports(cnt) := Determine_Order(supports(cnt));
                    286:                          cnt := cnt + mix(i);
                    287:                        end loop;
                    288:       end case;
                    289:     end Determine_Processing_Order;
                    290:
                    291:   -- INSTANTIATIONS OF THE GENERICS :
                    292:
                    293:     procedure Report_New_Simplices
                    294:                      ( t : in Triangulation; point : in Vector ) is
                    295:
                    296:     -- DESCRIPTION :
                    297:     --   Writes the new simplices on file and computes their volume.
                    298:
                    299:       v : natural;
                    300:
                    301:     begin
                    302:       new_line(file);
                    303:       put(file,"The new simplices by adding "); put(file,point);
                    304:       put_line(file," : ");
                    305:       put(file,n,t,v);
                    306:       put(file," with volume addition : ");
                    307:       put(file,vol,1); put(file," + "); put(file,v,1);
                    308:       vol := vol + v; put(file," = "); put(file,vol,1); put_line(file,".");
                    309:     end Report_New_Simplices;
                    310:     procedure R_Dynamic_Lifting is
                    311:       new Dynamic_Triangulations.Dynamic_Lifting_with_New(Report_New_Simplices);
                    312:
                    313:     procedure Collect_Flattening ( t : in Triangulation; l : List ) is
                    314:
                    315:     -- DESCRIPTION :
                    316:     --   Updates the subdivision mixsub with the flattened cells.
                    317:     --   The triangulation on entry contains the whole triangulation,
                    318:     --   not just the new cells.
                    319:
                    320:       cells : Mixed_Subdivision;
                    321:
                    322:     begin
                    323:       if Is_Null(mixsub)
                    324:        then cells := Deep_Create(n,t);
                    325:        else cells := Non_Flat_Deep_Create(n,t);
                    326:             Construct(Head_Of(mixsub),cells);
                    327:       end if;
                    328:       Flatten(cells);
                    329:       mixsub := cells;
                    330:     end Collect_Flattening;
                    331:
                    332:     procedure Report_Flattening
                    333:                     ( t : in Triangulation; l : in List ) is
                    334:
                    335:     -- DESCRIPTION :
                    336:     --   Writes the list of lifted points and the triangulation on file
                    337:     --   and updates the mixed subdivision.
                    338:
                    339:     begin
                    340:       new_line(file);
                    341:       put_line(file,"The list of lifted points before flattening : ");
                    342:       put(file,l);
                    343:       new_line(file);
                    344:       put_line(file,"The triangulation before flattening : ");
                    345:       put(file,n,t,vol);
                    346:       put(file," with volume "); put(file,vol,1); put_line(file,".");
                    347:       Collect_Flattening(t,l);
                    348:     end Report_Flattening;
                    349:     procedure C_Dynamic_Lifting is
                    350:       new Dynamic_Triangulations.Dynamic_Lifting_with_Flat(Collect_Flattening);
                    351:     procedure F_Dynamic_Lifting is
                    352:       new Dynamic_Triangulations.Dynamic_Lifting_with_Flat(Report_Flattening);
                    353:     procedure FR_Dynamic_Lifting is
                    354:       new Dynamic_Triangulations.Dynamic_Lifting_with_Flat_and_New
                    355:            ( Before_Flattening     => Report_Flattening,
                    356:              Process_New_Simplices => Report_New_Simplices);
                    357:
                    358:     procedure Report_New_Cells
                    359:                     ( mixsub : in out Mixed_Subdivision;
                    360:                       i : in natural; point : in Vector ) is
                    361:
                    362:     -- DESCRIPTION :
                    363:     --   Writes the new mixed cells on file and computes the mixed volume.
                    364:
                    365:       v : natural;
                    366:
                    367:     begin
                    368:       if not Is_Null(mixsub)
                    369:        then
                    370:          new_line(file);
                    371:          put(file,"The new mixed cells by adding "); put(file,point);
                    372:          new_line(file);
                    373:          put(file," to the "); put(file,i,1); put_line(file,"th component : ");
                    374:          put(file,n,mix.all,mixsub,v);
                    375:          put(file," with volume addition : ");
                    376:          put(file,vol,1); put(file," + "); put(file,v,1);
                    377:          vol := vol + v; put(file," = "); put(file,vol,1); new_line(file);
                    378:       end if;
                    379:     end Report_New_Cells;
                    380:     procedure R_Dynamic_Cayley is
                    381:       new Cayley_Trick.Dynamic_Cayley_with_New(Report_New_Cells);
                    382:     procedure Rt_Dynamic_Cayley is
                    383:       new Cayley_Trick.Dynamic_Cayley_with_Newt(Report_New_Cells);
                    384:     procedure R_Dynamic_Lifting is
                    385:       new Dynamic_Mixed_Subdivisions.Dynamic_Lifting_with_New(Report_New_Cells);
                    386:
                    387:     procedure Report_Flattening
                    388:                  ( mixsub : in out Mixed_Subdivision;
                    389:                    lifted : in Array_of_Lists ) is
                    390:
                    391:     -- DESCRIPTION :
                    392:     --   Writes the list of lifted points and the subdivision on file.
                    393:
                    394:     begin
                    395:       new_line(file);
                    396:       put_line(file,"The list of lifted points before flattening : ");
                    397:       for i in lifted'range loop
                    398:         put(file," points of "); put(file,i,1);
                    399:         put_line(file,"th component : ");
                    400:         put(file,lifted(i));
                    401:       end loop;
                    402:       new_line(file);
                    403:       put_line(file,"The mixed subdivision before flattening : ");
                    404:       put(file,n,mix.all,mixsub,vol);
                    405:       put(file," with volume "); put(file,vol,1); put_line(file,".");
                    406:     end Report_Flattening;
                    407:     procedure F_Dynamic_Cayley is
                    408:       new Cayley_Trick.Dynamic_Cayley_with_Flat(Report_Flattening);
                    409:     procedure Ft_Dynamic_Cayley is
                    410:       new Cayley_Trick.Dynamic_Cayley_with_Flatt(Report_Flattening);
                    411:     procedure FR_Dynamic_Cayley is
                    412:       new Cayley_Trick.Dynamic_Cayley_with_Flat_and_New
                    413:             (Before_Flattening => Report_Flattening,
                    414:              Process_New_Cells => Report_New_Cells);
                    415:     procedure FRt_Dynamic_Cayley is
                    416:       new Cayley_Trick.Dynamic_Cayley_with_Flat_and_Newt
                    417:             (Before_Flattening => Report_Flattening,
                    418:              Process_New_Cells => Report_New_Cells);
                    419:
                    420:     procedure Report_Flattening
                    421:                  ( mixsub : in out Mixed_Subdivision;
                    422:                    fs : in Face_Structures ) is
                    423:
                    424:     -- DESCRIPTION :
                    425:     --   Writes the list of lifted points and the subdivision on file.
                    426:
                    427:     begin
                    428:       new_line(file);
                    429:       put_line(file,"The lists of lifted points before flattening : ");
                    430:       for i in fs'range loop
                    431:         put(file," points of "); put(file,i,1);
                    432:         put_line(file,"th component : ");
                    433:         put(file,fs(i).l);
                    434:       end loop;
                    435:       new_line(file);
                    436:       put_line(file,"The mixed subdivision before flattening : ");
                    437:       put(file,n,mix.all,mixsub,vol);
                    438:       put(file," with volume "); put(file,vol,1); put_line(file,".");
                    439:     end Report_Flattening;
                    440:     procedure F_Dynamic_Lifting is
                    441:       new Dynamic_Mixed_Subdivisions.Dynamic_Lifting_with_Flat
                    442:             (Report_Flattening);
                    443:     procedure FR_Dynamic_Lifting is
                    444:       new Dynamic_Mixed_Subdivisions.Dynamic_Lifting_with_Flat_and_New
                    445:             (Before_Flattening => Report_Flattening,
                    446:              Process_New_Cells => Report_New_Cells);
                    447:
                    448:   -- MAIN CONSTRUCTORS :
                    449:
                    450:     procedure Compute_Triangulation is
                    451:
                    452:     -- DESCRIPTION :
                    453:     --   Application of the dynamic lifting algorithm
                    454:     --   to compute a triangulation of one polytope.
                    455:
                    456:       t : Triangulation;
                    457:       support,lifted,lifted_last : List;
                    458:       arlifted : Array_of_Lists(mix'range);
                    459:
                    460:     begin
                    461:       support := supports(supports'first);
                    462:       if verpts
                    463:        then Vertex_Points(file,support);
                    464:       end if;
                    465:       new_line(file);
                    466:       put_line(file,"CREATION OF THE TRIANGULATION :");
                    467:       new_line(file);
                    468:       tstart(timer);
                    469:       if reportnew
                    470:        then
                    471:          if reportflat
                    472:           then FR_Dynamic_Lifting(support,order,inter,max,lifted,lifted_last,t);
                    473:           else R_Dynamic_Lifting(support,order,inter,max,lifted,lifted_last,t);
                    474:          end if;
                    475:        elsif reportflat
                    476:            then
                    477:              F_Dynamic_Lifting(support,order,inter,max,lifted,lifted_last,t);
                    478:            elsif subonfile
                    479:                then
                    480:                  C_Dynamic_Lifting
                    481:                                 (support,order,inter,max,lifted,lifted_last,t);
                    482:                else
                    483:                  Dynamic_Lifting(support,order,inter,max,lifted,lifted_last,t);
                    484:       end if;
                    485:       tstop(timer);
                    486:       new_line(file);
                    487:       print_times(file,timer,"computing the triangulation");
                    488:       new_line(file);
                    489:       put_line(file,"THE LIFTED SUPPORTS :"); new_line(file);
                    490:       put(file,lifted);
                    491:       new_line(file);
                    492:       put_line(file,"THE TRIANGULATION :"); new_line(file);
                    493:       tstart(timer);
                    494:       put(file,n,t,vol);
                    495:       tstop(timer);
                    496:       new_line(file);
                    497:       put(file,"The volume : "); put(file,vol,1); new_line(file);
                    498:       new_line(file);
                    499:       print_times(file,timer,"computing the volume");
                    500:       if subonfile
                    501:        then if Is_Null(mixsub)
                    502:              then put(subfile,n,1); new_line(subfile);
                    503:                   put(subfile,1,1); new_line(subfile); -- type of mixture
                    504:                   put(subfile,n,t);
                    505:              else declare
                    506:                     lastcells : Mixed_Subdivision := Non_Flat_Deep_Create(n,t);
                    507:                   begin
                    508:                     Construct(Head_Of(mixsub),lastcells);
                    509:                     mixsub := lastcells;
                    510:                     put(subfile,n,mix.all,mixsub);
                    511:                   end;
                    512:             end if;
                    513:             Close(subfile);
                    514:       end if;
                    515:       mv := vol;
                    516:     end Compute_Triangulation;
                    517:
                    518:     procedure Compute_Cayley_Triangulation is
                    519:
                    520:     -- DESCRIPTION :
                    521:     --   Application of the dynamic lifting algorithm to compute a mixed
                    522:     --   subdivision of a tuple of polytopes by means of the Cayley trick.
                    523:
                    524:       supp,lifted : Array_of_Lists(1..r);
                    525:       t : Triangulation;
                    526:       numtri,mr : natural;
                    527:       newperms : Link_to_Vector;
                    528:
                    529:     begin
                    530:       if verpts
                    531:        then Vertex_Points(file,mix,supports);
                    532:             Clear(mix);
                    533:             Compute_Mixture(supports,mix,newperms);
                    534:             Write_Type_of_Mixture(file,mix.all,newperms.all);
                    535:       end if;
                    536:       mr := mix'last;
                    537:       supp(1..mr) := Typed_Lists(mix.all,supports);
                    538:       new_line(file);
                    539:       put_line(file,"CREATION OF THE MIXED SUBDIVISION :");
                    540:       new_line(file);
                    541:       tstart(timer);
                    542:       if reportnew
                    543:        then
                    544:          if reportflat
                    545:           then
                    546:             if minkpoly > 0
                    547:              then FRt_Dynamic_Cayley
                    548:                     (n,mix.all,supp(1..mr),order,inter,max,lifted(1..mr),t);
                    549:              else FR_Dynamic_Cayley
                    550:                     (n,mix.all,supp(1..mr),order,inter,max,
                    551:                      lifted(1..mr),mixsub,numtri);
                    552:             end if;
                    553:           else
                    554:             if minkpoly > 0
                    555:              then Rt_Dynamic_Cayley
                    556:                     (n,mix.all,supp(1..mr),order,inter,max,lifted(1..mr),t);
                    557:              else R_Dynamic_Cayley
                    558:                     (n,mix.all,supp(1..mr),order,inter,max,lifted(1..mr),
                    559:                      mixsub,numtri);
                    560:             end if;
                    561:          end if;
                    562:        elsif reportflat
                    563:            then
                    564:              if minkpoly > 0
                    565:               then Ft_Dynamic_Cayley
                    566:                      (n,mix.all,supp(1..mr),order,inter,max,lifted(1..mr),t);
                    567:               else F_Dynamic_Cayley
                    568:                      (n,mix.all,supp(1..mr),order,inter,max,
                    569:                       lifted(1..mr),mixsub,numtri);
                    570:              end if;
                    571:            else
                    572:              if minkpoly > 0
                    573:               then Dynamic_Cayley
                    574:                      (n,mix.all,supp(1..mr),order,inter,max,lifted(1..mr),t);
                    575:               else Dynamic_Cayley
                    576:                      (n,mix.all,supp(1..mr),order,inter,max,
                    577:                       lifted(1..mr),mixsub,numtri);
                    578:              end if;
                    579:       end if;
                    580:       tstop(timer);
                    581:       new_line(file);
                    582:       print_times(file,timer,"Computing the mixed subdivision");
                    583:       new_line(file);
                    584:       put_line(file,"THE LIFTED SUPPORTS :");
                    585:       new_line(file);
                    586:       put(file,lifted);
                    587:       if minkpoly > 0
                    588:        then declare
                    589:               alltri : boolean := (minkpoly > 1);
                    590:             begin
                    591:               Driver_for_Minkowski_Polynomials(file,n,mix.all,t,alltri,mixsub);
                    592:               numtri := Length_Of(t);
                    593:             end;
                    594:       end if;
                    595:       new_line(file);
                    596:       put_line(file,"THE MIXED SUBDIVISION :");
                    597:       new_line(file);
                    598:       tstart(timer);
                    599:       put(file,n,mix.all,mixsub,vol);
                    600:       tstop(timer);
                    601:       new_line(file);
                    602:       put(file,"The mixed volume equals : "); put(file,vol,1);
                    603:       new_line(file);
                    604:       put(file,"Number of cells in auxiliary triangulation : ");
                    605:       put(file,numtri,1); new_line(file);
                    606:       new_line(file);
                    607:       print_times(file,timer,"Computing the mixed volume");
                    608:       if subonfile
                    609:        then put(subfile,n,mix.all,mixsub);
                    610:             Close(subfile);
                    611:       end if;
                    612:       mv := vol;
                    613:     end Compute_Cayley_Triangulation;
                    614:
                    615:     procedure Report_Results
                    616:                 ( file : in file_type; n : in natural; mix : in Link_to_Vector;
                    617:                   mixsub : in out Mixed_Subdivision;
                    618:                   fs : in Face_Structures ) is
                    619:     begin
                    620:       new_line(file);
                    621:       put_line(file,"THE LIFTED SUPPORTS :");
                    622:       new_line(file);
                    623:       for i in fs'range loop
                    624:         put(file,fs(i).l);  new_line(file);
                    625:       end loop;
                    626:       put_line(file,"THE MIXED SUBDIVISION :");
                    627:       new_line(file);
                    628:       tstart(timer);
                    629:       if r = 1
                    630:        then put(file,n,fs(fs'first).t,vol);
                    631:        else put(file,n,mix.all,mixsub,vol);
                    632:       end if;
                    633:       tstop(timer);
                    634:       new_line(file);
                    635:       put(file,"The mixed volume equals : "); put(file,vol,1); new_line(file);
                    636:       new_line(file);
                    637:       print_times(file,timer,"Computing the mixed volume");
                    638:       if subonfile
                    639:        then put(subfile,n,mix.all,mixsub);
                    640:             Close(subfile);
                    641:       end if;
                    642:       mv := vol;
                    643:     end Report_Results;
                    644:
                    645:     procedure Compute_Mixed_Subdivision is
                    646:
                    647:     -- DESCRIPTION :
                    648:     --   Application of the dynamic lifting algorithm
                    649:     --   to compute a mixed subdivision of a tuple of polytopes.
                    650:
                    651:       supp,lifted : Array_of_Lists(1..r);
                    652:       fs : Face_Structures(1..r);
                    653:       nbsucc,nbfail : Standard_Floating_Vectors.Vector(1..r) := (1..r => 0.0);
                    654:       mr : natural;
                    655:       newperms : Link_to_Vector;
                    656:
                    657:     begin
                    658:       if verpts
                    659:        then Vertex_Points(file,mix,supports);
                    660:             Clear(mix);
                    661:             Compute_Mixture(supports,mix,newperms);
                    662:             Write_Type_of_Mixture(file,mix.all,newperms.all);
                    663:       end if;
                    664:       mr := mix'last;
                    665:       supp(1..mr) := Typed_Lists(mix.all,supports);
                    666:       new_line(file);
                    667:       put_line(file,"CREATION OF THE MIXED SUBDIVISION :");
                    668:       new_line(file);
                    669:       tstart(timer);
                    670:       if reportnew
                    671:        then
                    672:          if reportflat
                    673:           then FR_Dynamic_Lifting
                    674:                  (n,mix.all,supp(1..mr),order,inter,conmv,max,mixsub,
                    675:                   fs(1..mr),nbsucc(1..mr),nbfail(1..mr));
                    676:           else R_Dynamic_Lifting
                    677:                  (n,mix.all,supp(1..mr),order,inter,conmv,max,mixsub,
                    678:                   fs(1..mr),nbsucc(1..mr),nbfail(1..mr));
                    679:          end if;
                    680:        elsif reportflat
                    681:            then F_Dynamic_Lifting
                    682:                  (n,mix.all,supp(1..mr),order,inter,conmv,max,mixsub,
                    683:                   fs(1..mr),nbsucc(1..mr),nbfail(1..mr));
                    684:            else Dynamic_Lifting
                    685:                  (n,mix.all,supp(1..mr),order,inter,conmv,max,mixsub,
                    686:                   fs(1..mr),nbsucc(1..mr),nbfail(1..mr));
                    687:       end if;
                    688:       tstop(timer);
                    689:       Pruning_Statistics(file,nbsucc(1..mr),nbfail(1..mr));
                    690:       new_line(file);
                    691:       print_times(file,timer,"Computing the mixed subdivision");
                    692:       Report_Results(file,n,mix,mixsub,fs(1..mr));
                    693:     end Compute_Mixed_Subdivision;
                    694:
                    695:     procedure Solve_Coefficient_System is
                    696:
                    697:     -- DESCRIPTION :
                    698:     --   Application of the dynamic lifting algorithm
                    699:     --   to compute a mixed subdivision of a tuple of polytopes and
                    700:     --   to solve a start system, with randomized coefficients.
                    701:
                    702:       supp : Array_of_Lists(1..r);
                    703:       fs : Face_Structures(1..r);
                    704:       lifted : Array_of_Lists(1..r);
                    705:       numtri : natural := 0;
                    706:       lif,lif_last : List;
                    707:       nbsucc,nbfail : Standard_Floating_Vectors.Vector(1..r) := (1..r => 0.0);
                    708:       mr : natural;
                    709:       newperms : Link_to_Vector;
                    710:
                    711:     begin
                    712:       if verpts
                    713:        then Vertex_Points(file,mix,supports);
                    714:             if r > 1
                    715:              then Clear(mix);
                    716:                   Compute_Mixture(supports,mix,newperms);
                    717:                   Write_Type_of_Mixture(file,mix.all,newperms.all);
                    718:                   qq := Permute(qq,newperms);
                    719:                   for i in supports'range loop
                    720:                     qq(i) := Select_Terms(qq(i),supports(i));
                    721:                   end loop;
                    722:             end if;
                    723:       end if;
                    724:       mr := mix'last;
                    725:       supp(1..mr) := Typed_Lists(mix.all,supports);
                    726:       new_line(file);
                    727:       put_line(file,"SOLVING THE RANDOM COEFFICIENT SYSTEM :");
                    728:       new_line(file);
                    729:       tstart(timer);
                    730:       if mix'last = mix'first
                    731:        then
                    732:          Dynamic_Unmixed_Solve
                    733:                (file,n,supp(supp'first),order,inter,max,fs(fs'first).l,
                    734:                 fs(fs'first).last,fs(fs'first).t,qq,qqsols);
                    735:        else
                    736:          if caytrick
                    737:           then
                    738:             Dynamic_Cayley_Solve(file,n,mix.all,supp(1..mr),order,inter,max,
                    739:                                  lifted(1..mr),mixsub,numtri,qq,qqsols);
                    740:             for i in 1..mr loop
                    741:               fs(i).l := lifted(i);
                    742:             end loop;
                    743:           else
                    744:             Dynamic_Mixed_Solve
                    745:                      (file,n,mix.all,supp(1..mr),order,inter,conmv,max,mixsub,
                    746:                       fs(1..mr),nbsucc(1..mr),nbfail(1..mr),qq,qqsols);
                    747:          end if;
                    748:       end if;
                    749:       tstop(timer);
                    750:       if mix'last > mix'first and not caytrick
                    751:        then Pruning_Statistics(file,nbsucc(1..mr),nbfail(1..mr));
                    752:       end if;
                    753:       new_line(file);
                    754:       print_times(file,timer,"Computing the solution list");
                    755:       Report_Results(file,n,mix,mixsub,fs(1..mr));
                    756:       q := qq; qsols := qqsols;
                    757:       if not ranstart
                    758:        then put(solsft,qqsols);
                    759:             Close(solsft);
                    760:       end if;
                    761:       if ranstart
                    762:        then new_line(qft); put_line(qft,"THE SOLUTIONS :"); new_line(qft);
                    763:             put(qft,Length_Of(qqsols),n,qqsols);
                    764:             Close(qft);
                    765:       end if;
                    766:     end Solve_Coefficient_System;
                    767:
                    768:   begin
                    769:     new_line; put_line(welcome);
                    770:    -- READING GENERAL INPUT INFORMATION :
                    771:     supports := Create(p);
                    772:     new_line;
                    773:     put("Do you want to enforce a type of mixture ? (y/n) ");
                    774:     Ask_Yes_or_No(ans);
                    775:     if ans /= 'y'
                    776:      then Compute_Mixture(supports,mix,perms); r := mix'last;
                    777:      else put("Give number of different supports : "); Read_Natural(r);
                    778:           put("Give vector of occurrences : "); get(r,mix);
                    779:           perms := new Vector(1..n);
                    780:           for i in perms'range loop
                    781:             perms(i) := i;
                    782:           end loop;
                    783:     end if;
                    784:     Write_Type_of_Mixture(file,mix.all,perms.all);
                    785:    -- DETERMINE THE GLOBAL SWITCHES :
                    786:     put("Do you first want to extract the vertex points ? (y/n) ");
                    787:     Ask_Yes_or_No(ans);
                    788:     verpts := (ans = 'y');
                    789:     inter := not verpts;
                    790:     put("Do you have a maximum lifting value ? (y/n) ");
                    791:     Ask_Yes_or_No(ans);
                    792:     if ans = 'y'
                    793:      then put("  Give the maximum lifting value : ");
                    794:           Read_Positive(max);
                    795:      else max := 0;
                    796:     end if;
                    797:     Determine_Processing_Order(supports,mix,order);
                    798:     if (r > 1)
                    799:      then new_line;
                    800:           put_line("MENU for Cayley trick : ");
                    801:           put_line("  0. No Cayley trick, pruning for mixed cells.");
                    802:           put_line("  1. Cayley trick : auxiliary triangulation.");
                    803:           put_line("  2. Cayley trick with Minkowski-polynomial.");
                    804:           put_line("  3. Cayley trick with all subdivisions.");
                    805:           put("Type 0,1,2, or 3 : ");
                    806:           Ask_Alternative(ans,"0123");
                    807:           caytrick := not (ans = '0');
                    808:           case ans is
                    809:             when '2' => minkpoly := 1;
                    810:             when '3' => minkpoly := 2;
                    811:             when others => minkpoly := 0;
                    812:           end case;
                    813:           if not caytrick
                    814:            then put("Do you want online checks on zero contributions ? (y/n) ");
                    815:                 Ask_Yes_or_No(ans);
                    816:                 conmv := (ans = 'y');
                    817:            else conmv := false;
                    818:           end if;
                    819:      else caytrick := false; conmv := false;
                    820:     end if;
                    821:     put("Do you want to have the subdivision on separate file ? (y/n) ");
                    822:     Ask_Yes_or_No(ans);
                    823:     if ans = 'y'
                    824:      then subonfile := true;
                    825:           put_line("Reading the name of the file.");
                    826:           Read_Name_and_Create_File(subfile);
                    827:      else subonfile := false;
                    828:     end if;
                    829:     new_line;
                    830:     put("Are the cells to be written on file, during computation ? (y/n) ");
                    831:     Ask_Yes_or_No(ans);
                    832:     reportnew := (ans = 'y');
                    833:     put("Are the cells to be written on file, before flattening ? (y/n) ");
                    834:     Ask_Yes_or_No(ans);
                    835:     reportflat := (ans = 'y');
                    836:     permp := Permute(p,perms);
                    837:     Driver_for_Polyhedral_Continuation
                    838:       (file,permp,0,byebye,qq,qft,solsft,tosolve,ranstart,contrep);
                    839:    -- HANDLING THE UNMIXED AND THE MIXED CASE SEPARATELY :
                    840:     if not tosolve
                    841:      then if r = 1
                    842:            then Compute_Triangulation;
                    843:            else if caytrick
                    844:                  then Compute_Cayley_Triangulation;
                    845:                  else Compute_Mixed_Subdivision;
                    846:                 end if;
                    847:           end if;
                    848:      else Solve_Coefficient_System;
                    849:     end if;
                    850:   end Driver_for_Dynamic_Mixed_Volume_Computation;
                    851:
                    852: end Drivers_for_Dynamic_Lifting;

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>