Return to dynamic_polyhedral_continuation.adb CVS log

Up to [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Dynlift

Import the second public release of PHCpack.

OKed by Jan Verschelde.

with integer_io;                         use integer_io;
with Standard_Complex_Numbers;           use Standard_Complex_Numbers;
with Standard_Complex_Laur_Polys;        use Standard_Complex_Laur_Polys;
with Standard_Complex_Laur_Systems;      use Standard_Complex_Laur_Systems;
with Standard_Poly_Laur_Convertors;      use Standard_Poly_Laur_Convertors;
with Fewnomial_System_Solvers;           use Fewnomial_System_Solvers;
with Transforming_Laurent_Systems;
with Simplices;                          use Simplices;
with Dynamic_Triangulations;             use Dynamic_Triangulations;
with Cayley_Trick;                       use Cayley_Trick;
with Integer_Lifting_Utilities;          use Integer_Lifting_Utilities;
with Unfolding_Subdivisions;             use Unfolding_Subdivisions;
with Integer_Mixed_Subdivisions_io;      use Integer_Mixed_Subdivisions_io;
with Integer_Polyhedral_Continuation;    use Integer_Polyhedral_Continuation;

package body Dynamic_Polyhedral_Continuation is

-- CAUTION : PATCH FOR FEWNOMIALS => NO SHIFT !!!

-- AUXILIAIRIES :

  procedure Flatten ( t : in out Term ) is

  -- DESCRIPTION :
  --   Flattens the Laurent term, i.e., the last exponent of t becomes zero.

  begin
    t.dg(t.dg'last) := 0;
  end Flatten;

  procedure Flatten ( p : in out Poly ) is

  -- DESCRIPTION :
  --   Flattens the Laurent polynomial,
  --   i.e., the last exponent of every monomial becomes zero.

    procedure Flatten_Term ( t : in out Term; cont : out boolean ) is
    begin
      Flatten(t); cont := true;
    end Flatten_Term;
    procedure Flatten_Terms is new Changing_Iterator(Flatten_Term);

  begin
    Flatten_Terms(p);
  end Flatten;

  procedure Flatten ( p : in out Laur_Sys ) is

  -- DESCRIPTION :
  --   Flattens the Laurent polynomial system,
  --   i.e., the last exponent of every monomial becomes zero.

  begin
    for i in p'range loop
      Flatten(p(i));
    end loop;
  end Flatten;

  function Non_Flattened_Points ( l : List ) return List is

  -- DESCRIPTION :
  --   Returns the list of points with last coordinate /= 0.
   
    tmp,res,res_last : List;
    pt : Link_to_Vector;

  begin
    tmp := l;
    while not Is_Null(tmp) loop
      pt := Head_Of(tmp);
      if pt(pt'last) /= 0
       then Append(res,res_last,pt);
      end if;
      tmp := Tail_Of(tmp);
    end loop;
    return res;
  end Non_Flattened_Points;

  function Non_Flattened_Points ( l : Array_of_Lists )
                                return Array_of_Lists is

  -- DESCRIPTION :
  --   Returns the points with last coordinate /= 0.

    res : Array_of_Lists(l'range);

  begin
    for i in l'range loop
      res(i) := Non_Flattened_Points(l(i));
    end loop;
    return res;
  end Non_Flattened_Points;

  function Non_Flattened_Points ( fs : Face_Structures )
                                return Array_of_Lists is

  -- DESCRIPTION :
  --   Returns the points with last coordinate /= 0.

    res : Array_of_Lists(fs'range);

  begin
    for i in fs'range loop
      res(i) := Non_Flattened_Points(fs(i).l);
    end loop;
    return res;
  end Non_Flattened_Points;

  function Non_Flattened_Points
               ( mix : Vector; mixsub : Mixed_Subdivision )
               return Array_of_Lists is

  -- DESCRIPTION :
  --   Returns the points with last coordinate /= 0.

    res,res_last : Array_of_Lists(mix'range);
    tmp : Mixed_Subdivision := mixsub;

  begin
    while not Is_Null(tmp) loop
      declare
        mic : Mixed_Cell := Head_Of(tmp);
      begin
        for i in mic.pts'range loop
          Deep_Concat_Diff(res(i),res_last(i),Non_Flattened_Points(mic.pts(i)));
        end loop;
      end;
      tmp := Tail_Of(tmp);
    end loop;
    return res;
  end Non_Flattened_Points;

  function Convert ( s : Simplex ) return Mixed_Cell is

    res : Mixed_Cell;

  begin
    res.nor := new vector'(Normal(s));
    res.pts := new Array_of_Lists(1..1);
    res.pts(1) := Shallow_Create(Vertices(s));
    return res;
  end Convert; 

  function Convert ( normal : in vector; s : Simplex ) return Mixed_Cell is

    res : Mixed_Cell;

  begin
    res.nor := new vector'(normal);
    res.pts := new Array_of_Lists(1..1);
    res.pts(1) := Shallow_Create(Vertices(s));
    return res;
  end Convert;

  function Non_Flattened_Cells
                 ( flatnor : Link_to_Vector; mixsub : Mixed_Subdivision )
                 return Mixed_Subdivision is

  -- DESCRIPTION :
  --   Returns the cells which are not flattened.

    tmp,res,res_last : Mixed_Subdivision;
    mic : Mixed_Cell;

  begin
    tmp := mixsub;
    while not Is_Null(tmp) loop
      mic := Head_Of(tmp);
      if mic.nor.all /= flatnor.all
       then Append(res,res_last,mic);
      end if;
      tmp := Tail_Of(tmp);
    end loop;
    return res;
  end Non_Flattened_Cells;

  function Convert ( t : Triangulation ) return Mixed_Subdivision is

    res,res_last : Mixed_Subdivision;
    tmp : Triangulation := t;

  begin
    while not Is_Null(tmp) loop
      declare
        mic : Mixed_Cell := Convert(Head_Of(tmp));
      begin
        Append(res,res_last,mic);
      end;
      tmp := Tail_Of(tmp);
    end loop;
    return res;
  end Convert;

  function Non_Flattened_Cells
                 ( flatnor : Link_to_Vector; t : Triangulation )
                 return Mixed_Subdivision is

    tmp : Triangulation;
    res,res_last : Mixed_Subdivision;
   
  begin
    tmp := t;
    while not Is_Null(tmp) loop
      declare
        s : Simplex := Head_Of(tmp);
        nor : constant vector := Normal(s);
      begin
        if nor /= flatnor.all
         then declare
                mic : Mixed_Cell := Convert(nor,s);
              begin
                Append(res,res_last,mic);
              end;
        end if;
      end;
      tmp := Tail_Of(tmp);
    end loop;
    return res;
  end Non_Flattened_Cells;

-- UTILITIES :

  function Pointer_to_Last ( l : Solution_List ) return Solution_List is

  -- DESCRIPTION :
  --   Returns a pointer to the last element of l.

  begin
    if Is_Null(l)
     then return l;
     elsif Is_Null(Tail_Of(l))
         then return l;
         else return Pointer_to_Last(Tail_Of(l));
    end if;
  end Pointer_to_Last;

  function Initialize_Polyhedral_Homotopy
                ( n :  natural; mix : Vector; fs : Face_Structures;
                  p : Laur_Sys ) return Laur_Sys is

  -- DESCRIPTION :
  --   Initializes the polyhedral homotopy w.r.t. to the already lifted
  --   points in the face structure.

    res : Laur_Sys(p'range);

  begin
    if not Is_Null(fs(fs'first).l)
     then declare
            lifted : Array_of_Lists(fs'range);
          begin
            for i in lifted'range loop
              lifted(i) := fs(i).l;
            end loop;
            res := Perform_Lifting(n,mix,lifted,p);
          end;
    end if;
    return res;
  end Initialize_Polyhedral_Homotopy;

  function Initialize_Polyhedral_Homotopy
                ( n : natural; l : List; p : Laur_Sys ) return Laur_Sys is

  -- DESCRIPTION :
  --   Initializes the polyhedral homotopy w.r.t. to the already lifted
  --   points in the list.

    res : Laur_Sys(p'range);

  begin
    if not Is_Null(l)
     then declare
            lifted : Array_of_Lists(p'range);
            mix : Vector(1..1) := (1..1 => n);
          begin
            for i in lifted'range loop
              lifted(i) := l;
            end loop;
            res := Perform_Lifting(n,mix,lifted,p);
          end;
    end if;
    return res;
  end Initialize_Polyhedral_Homotopy;

  function Select_Terms ( p : Poly; l : List ) return Poly is

  -- DESCRIPTION :
  --   Given a tuple of lifted points, selected terms of p, whose exponents
  --   occur in l, will be returned.

    res : Poly := Null_Poly;
    tmp : List := l;
    point : Link_to_Vector;

  begin
    while not Is_Null(tmp) loop
      point := Head_Of(tmp);
      declare
        d : degrees := new Vector(point'first..point'last-1);
        t : Term;
      begin
        d.all := point(point'first..point'last-1);
        t.cf := Coeff(p,d);
        t.dg := d;
        Add(res,t);
        Standard_Integer_Vectors.Clear
          (Standard_Integer_Vectors.Link_to_Vector(d));
      end;
      tmp := Tail_Of(tmp);
    end loop;
    return res;
  end Select_Terms;

  function Select_Terms ( p : Laur_Sys; mix : Vector; l : Array_of_Lists )
                        return Laur_Sys is

  -- DESCRIPTION :
  --   Given a tuple of lifted points, selected terms of p, whose exponents
  --   occur in l, will be returned.

    res : Laur_Sys(p'range);
    cnt : natural := p'first;

  begin
    for i in mix'range loop
      for j in 1..mix(i) loop
        res(cnt) := Select_Terms(p(cnt),l(i));
        cnt := cnt + 1;
      end loop;
    end loop;
    return res;
  end Select_Terms;

  procedure Update_Polyhedral_Homotopy
               ( p : in Laur_Sys; point : in Vector; i : in integer;
                 hom : in out Laur_Sys ) is

  -- DESCRIPTION :
  --   Given the lifted point of the ith support list,
  --   the ith polynomial of hom will be updated with a new term.

    d : degrees
      := new Standard_Integer_Vectors.Vector(point'first..point'last-1);
    t : Term;

  begin
    d.all := point(point'first..point'last-1);
    t.cf := Coeff(p(i),d);
    t.dg := new Standard_Integer_Vectors.Vector'(point);
    Add(hom(i),t);
    Standard_Integer_Vectors.Clear(Standard_Integer_Vectors.Link_to_Vector(d));
    Clear(t);
  end Update_Polyhedral_homotopy;

  procedure Update_Polyhedral_Homotopy
               ( p : in Laur_Sys; l : in List; i : in integer;
                 hom : in out Laur_Sys ) is

  -- DESCRIPTION :
  --   Given a list of lifted points of the ith support list,
  --   the ith polynomial of hom will be updated with new terms.

    tmp : List := l;

  begin
    while not Is_Null(tmp) loop
      Update_Polyhedral_Homotopy(p,Head_Of(tmp).all,i,hom);
      tmp := Tail_Of(tmp);
    end loop;
  end Update_Polyhedral_Homotopy;

  procedure Update_Polyhedral_Homotopy
               ( p : in Laur_Sys; lifted : in Array_of_Lists;
                 mix : in Vector; hom : in out Laur_Sys ) is

  -- DESCRIPTION :
  --   Given a lists of lifted points, the polyhedral homotopy hom
  --   will be updated with new terms.

    cnt : natural := hom'first;

  begin
    for i in mix'range loop
      for j in 1..mix(i) loop
        Update_Polyhedral_Homotopy(p,lifted(i),cnt,hom);
        cnt := cnt + 1;
      end loop;
    end loop;
  end Update_Polyhedral_Homotopy;

  procedure Update_Polyhedral_Homotopy
               ( p : in Laur_Sys; lifted : in List; hom : in out Laur_Sys ) is

  -- DESCRIPTION :
  --   Given a list of lifted points, the unmixed polyhedral homotopy hom
  --   will be updated with new terms.

  begin
    for i in hom'range loop
      Update_Polyhedral_Homotopy(p,lifted,i,hom);
    end loop;
  end Update_Polyhedral_Homotopy;

  procedure Solve_by_Unfolding
                ( file : in file_type; p,hom : in Laur_Sys; n : in natural;
                  mix : in Vector; nor : in Link_to_Vector;
                  mixsub : in Mixed_Subdivision;
                  sols : in out Solution_List ) is

  -- DESCRIPTION :
  --   All cells in the given mixed subdivision have the same normal.
  --   The solutions which correspond to these cells will be computed
  --   by means of the given polyhedral homotopy.

  -- ON ENTRY :
  --   file       to write intermediate results on;
  --   p          randomized system, without lifting;
  --   hom        the polyhedral homotopy;
  --   n          dimension before lifting;
  --   mix        type of mixture;
  --   nor        the same normal to all cells in the subdivision;
  --   mixsub     the mixed subdivision all with the same normal nor.

  -- ON RETURN :
  --   sols       the solutions of p, w.r.t. the cells in mixsub.

    work : Mixed_Subdivision;
    mv : natural;
    homsub : Laur_Sys(p'range);
    first : boolean := true;
    flatnor : Link_to_Vector;
    sols_last : Solution_List;

    procedure Process ( mic : in Mixed_Cell; newpts : in Array_of_Lists ) is

      subsys : Laur_Sys(p'range) := Select_Terms(p,mix,mic.pts.all);
      subsols : Solution_List;
      fail : boolean;

    begin
      Transforming_Laurent_Systems.Shift(subsys);  -- patch for Fewnomials !!
      Solve(subsys,subsols,fail);
      put(file,"Number of solutions of subsystem : ");
      put(file,Length_Of(subsols),1); new_line(file);
      if first
       then homsub := Perform_Lifting(n,mix,mic.pts.all,p);
            sols := subsols;
            first := false;
       else Update_Polyhedral_Homotopy(p,newpts,mix,homsub);
            Mixed_Continuation(file,homsub,mic.nor.all,subsols);
            Set_Continuation_Parameter(sols,Create(0.0));
            Mixed_Continuation(file,homsub,flatnor.all,sols);
            Flatten(homsub);
            sols_last := Pointer_to_Last(sols);
            Concat(sols,sols_last,subsols);
            Shallow_Clear(subsols);
      end if;
      Clear(subsys);
    end Process;
    procedure Solve_Subsystem is new Unfolding(Process);

  begin
    new_line(file);
    put_line(file,"****  SOLVING BY UNFOLDING  ****");
    new_line(file);
    flatnor := new vector(1..n+1);
    flatnor(1..n) := (1..n => 0);
    flatnor(n+1) := 1;
    Copy(mixsub,work);
    put(file,n,mix,work,mv);
    Solve_Subsystem(work);
    --Deep_Clear(work);
    Clear(flatnor); Clear(homsub);
    Set_Continuation_Parameter(sols,Create(0.0));
    Mixed_Continuation(file,hom,nor.all,sols);
  end Solve_by_Unfolding;

  procedure Solve_with_Unfolding
                ( file : in file_type; p,hom : in Laur_Sys; n : in natural;
                  mix : in Vector; mixsub : in Mixed_Subdivision;
                  sols : in out Solution_List ) is

  -- DESCRIPTION :
  --   Computes the new solutions corresponding to the cells in the
  --   mixed subdivision.

    sols_last : Solution_List;
    newnor : List := Different_Normals(mixsub);
    tmp : List := newnor;
    normal : Link_to_Vector;

  begin
    while not Is_Null(tmp) loop
      normal := Head_Of(tmp);
      declare
        newcells : Mixed_Subdivision := Extract(normal.all,mixsub);
        bigcell : Mixed_Subdivision;
        newsols : Solution_List;
      begin
      -- put_line(file,"THE LIST OF NEW CELLS"); put(file,n,mix,newcells);
        if Length_Of(newcells) = 1
         then Mixed_Solve(file,hom,mix,newcells,newsols);
      --   else Solve_by_Unfolding(file,p,hom,n,mix,normal,newcells,newsols);
         else bigcell := Merge_Same_Normal(newcells);
             -- put_line(file,"THE BIG CELL"); put(file,n,mix,bigcell);
              Mixed_Solve(file,hom,mix,bigcell,newsols);
             -- Deep_Clear(bigcell);
              Shallow_Clear(bigcell);
        end if;
        sols_last := Pointer_to_Last(sols);
        Concat(sols,sols_last,newsols);
        Shallow_Clear(newsols);
        Shallow_Clear(newcells);
      end;
      tmp := Tail_Of(tmp);
    end loop;
    Deep_Clear(newnor);
  end Solve_with_Unfolding;

-- DYNAMIC LIFTING FOR UNMIXED SYSTEMS :

  procedure Dynamic_Unmixed_Solve
                ( file : in file_type; n : in natural;
                  l : in List; order,inter : in boolean; maxli : in natural;
                  lifted,lifted_last : in out List; t : in out Triangulation;
                  q : in Poly_Sys; qsols : in out Solution_List ) is

    p : Laur_Sys(q'range) := Polynomial_to_Laurent_System(q);
    qt : Laur_Sys(q'range);         -- the polyhedral homotopy
    firstflat : boolean := true;    -- first time flattening
    flatnor : Link_to_Vector;       -- the flat normal
    qsols_last : Solution_List;
    mix : Vector(1..1) := (1..1 => n);

    procedure Solve_Before_Flattening
                 ( t : in Triangulation; lifted : in List ) is

    -- DESCRIPTION :
    --   Computes the new solutions, right before the subdivision
    --   is flattened.

      newpts : List;
      newcells : Mixed_Subdivision;

    begin
      if firstflat
       then newpts := lifted;
            newcells := Convert(t);
       else newcells := Non_Flattened_Cells(flatnor,t);
            newpts := Non_Flattened_Points(lifted);
      end if;
      Update_Polyhedral_Homotopy(p,newpts,qt);
      if not firstflat
       then new_line(file);
            put_line(file,"***  EXTENDING THE SOLUTIONS  ***");
            new_line(file);
            Set_Continuation_Parameter(qsols,Create(0.0));
            Mixed_Continuation(file,qt,flatnor.all,qsols);
      end if;
      Solve_with_Unfolding(file,p,qt,n,mix,newcells,qsols);
      if not firstflat
       then Shallow_Clear(newpts); Shallow_Clear(newcells);
       else firstflat := false;
      end if;
      Flatten(qt);
    end Solve_Before_Flattening;

    procedure S_Dynamic_Lifting is
      new Dynamic_Triangulations.Dynamic_Lifting_with_Flat
                 ( Before_Flattening => Solve_Before_Flattening );

  begin
    qt := Initialize_Polyhedral_Homotopy(n,lifted,p);
    flatnor := new vector(1..n+1);
    flatnor(1..n) := (1..n => 0);
    flatnor(n+1) := 1;
    S_Dynamic_Lifting(l,order,inter,maxli,lifted,lifted_last,t);
    Solve_Before_Flattening(t,lifted);
    Clear(flatnor); Clear(mix);
    Clear(qt); Clear(p);
  end Dynamic_Unmixed_Solve;

-- DYNAMIC LIFTING FOR THE CAYLEY TRICK :

  procedure Dynamic_Cayley_Solve
                ( file : in file_type; n : in natural; mix : in Vector;
                  supports : in Array_of_Lists; order,inter : in boolean;
                  maxli : in natural; lifted : in out Array_of_Lists;
                  mixsub : in out Mixed_Subdivision; numtri : out natural;
                  q : in Poly_Sys; qsols : in out Solution_List ) is

    p : Laur_Sys(q'range) := Polynomial_to_Laurent_System(q);
    qt : Laur_Sys(q'range);         -- the polyhedral homotopy
    firstflat : boolean := true;    -- first time flattening
    flatnor : Link_to_Vector;       -- the flat normal
    qsols_last : Solution_List;

    procedure Solve_Before_Flattening
                ( mixsub : in out Mixed_Subdivision;
                  lifted : in Array_of_Lists ) is

    -- DESCRIPTION :
    --   Computes the new solutions, right before the subdivision
    --   is flattened.

      newpts : Array_of_Lists(lifted'range);
      newcells : Mixed_Subdivision;

    begin
      if firstflat
       then newpts := lifted;
            newcells := mixsub;
       else newcells := Non_Flattened_Cells(flatnor,mixsub);
            newpts := Non_Flattened_Points(lifted);
      end if;
      Update_Polyhedral_Homotopy(p,newpts,mix,qt);
      if not Is_Null(qsols)
       then new_line(file);
            put_line(file,"***  EXTENDING THE SOLUTIONS  ***");
            new_line(file);
            Set_Continuation_Parameter(qsols,Create(0.0));
            Mixed_Continuation(file,qt,flatnor.all,qsols);
      end if;
      if not Is_Null(newcells)
       then
         Solve_with_Unfolding(file,p,qt,n,mix,newcells,qsols);
         if not firstflat
          then Shallow_Clear(newpts); Shallow_Clear(newcells);
          else firstflat := false;
         end if;
      end if;
      Flatten(qt);
    end Solve_Before_Flattening;

    procedure S_Dynamic_Lifting is
      new Cayley_Trick.Dynamic_Cayley_with_Flat
                ( Before_Flattening => Solve_Before_Flattening );

  begin
    flatnor := new vector(1..n+1);
    flatnor(1..n) := (1..n => 0);
    flatnor(n+1) := 1;
    S_Dynamic_Lifting(n,mix,supports,order,inter,maxli,lifted,mixsub,numtri);
    Solve_Before_Flattening(mixsub,lifted);
    Clear(flatnor);
  end Dynamic_Cayley_Solve;

-- DYNAMIC LIFTING FOR SEMI-MIXED SYSTEMS :

  procedure Dynamic_Mixed_Solve
                ( file : in file_type; n : in natural;
                  mix : in Vector; supports : in Array_of_Lists;
                  order,inter,conmv : in boolean; maxli : in natural;
                  mixsub : in out Mixed_Subdivision;
                  fs : in out Face_Structures;
                  nbsucc,nbfail : in out Standard_Floating_Vectors.Vector;
                  q : in Poly_Sys; qsols : in out Solution_List ) is

    p : Laur_Sys(q'range) := Polynomial_to_Laurent_System(q);
    qt : Laur_Sys(q'range);         -- the polyhedral homotopy
    firstflat : boolean := true;    -- first time flattening 
    flatnor : Link_to_Vector;       -- the flat normal
    qsols_last : Solution_List;

    procedure Solve_Before_Flattening
                ( mixsub : in out Mixed_Subdivision;
                  fs : in Face_Structures ) is

    -- DESCRIPTION :
    --   Computes the new solutions, right before the subdivision 
    --   is flattened.

      newpts : Array_of_Lists(fs'range);
      newcells : Mixed_Subdivision;
      newsols : Solution_List;

    begin
      if firstflat
       then for i in fs'range loop
              newpts(i) := fs(i).l;
            end loop;
            newcells := mixsub;
       else newcells := Non_Flattened_Cells(flatnor,mixsub);
            newpts := Non_Flattened_Points(fs);
      end if;
      Update_Polyhedral_Homotopy(p,newpts,mix,qt);
      if not firstflat
       then new_line(file);
            put_line(file,"***  EXTENDING THE SOLUTIONS  ***");
            new_line(file);
            Set_Continuation_Parameter(qsols,Create(0.0));
            Mixed_Continuation(file,qt,flatnor.all,qsols);
      end if;
      Mixed_Solve(file,qt,mix,newcells,newsols);
      qsols_last := Pointer_to_Last(qsols);
      Concat(qsols,qsols_last,newsols);
      Shallow_Clear(newsols);
      if not firstflat
       then Shallow_Clear(newpts); Shallow_Clear(newcells);
       else firstflat := false;
      end if;
      Flatten(qt);
    end Solve_Before_Flattening;

    procedure S_Dynamic_Lifting is
      new Dynamic_Mixed_Subdivisions.Dynamic_Lifting_with_Flat
                ( Before_Flattening => Solve_Before_Flattening );

  begin
    qt := Initialize_Polyhedral_Homotopy(n,mix,fs,p);
    flatnor := new vector(1..n+1);
    flatnor(1..n) := (1..n => 0);
    flatnor(n+1) := 1;
    S_Dynamic_Lifting(n,mix,supports,order,inter,conmv,maxli,mixsub,fs,
                      nbsucc,nbfail);
    Solve_Before_Flattening(mixsub,fs);
    Clear(flatnor);
    Clear(qt); Clear(p);
  end Dynamic_Mixed_Solve;

end Dynamic_Polyhedral_Continuation;