Return to homotopy.adb CVS log

Up to [local] / OpenXM_contrib / PHC / Ada / Homotopy

Import the second public release of PHCpack.

OKed by Jan Verschelde.

with unchecked_deallocation;
with Standard_Floating_Numbers;          use Standard_Floating_Numbers;
with Standard_Natural_Vectors;
with Standard_Complex_Vectors;           use Standard_Complex_Vectors;
with Standard_Complex_Matrices;          use Standard_Complex_Matrices;
with Standard_Complex_Polynomials;       use Standard_Complex_Polynomials;
with Standard_Complex_Poly_Functions;    use Standard_Complex_Poly_Functions;
with Standard_Complex_Poly_SysFun;       use Standard_Complex_Poly_SysFun;
with Standard_Complex_Jaco_Matrices;     use Standard_Complex_Jaco_Matrices;

package body Homotopy is

-- INTERNAL DATA -- to perform the numeric routines as fast as possible.

  type homtype is (nat,art);  -- natural or artificial parameter homotopy

  type homdata ( ht : homtype; n,n1 : natural ) is record

    p : Poly_Sys(1..n);                      -- target system
    pe : Eval_Poly_Sys(1..n);                -- evaluable form of target
    dh : Jaco_Mat(1..n,1..n1);               -- Jacobian matrix of homotopy
    dhe : Eval_Jaco_Mat(1..n,1..n1);         -- evaluable form of dh

    case ht is
      when nat =>
        i : integer;                         -- which variable is parameter
      when art =>
        q,h : Poly_Sys(1..n);                -- start system and homotopy
        qe,he : Eval_Poly_Sys(1..n);         -- evaluable form of q and h
        dpe,dqe : Eval_Jaco_Mat(1..n,1..n);  -- evaluable Jacobians
        k : positive;                        -- relaxation power
        gamma,beta : Vector(1..n);           -- accessibility constants
        linear : boolean;                    -- linear in t if true
      end case;

  end record;
 
  type homtp is access homdata;
  hom : homtp;

-- GENERAL AUXILIARIES :

  function Mul ( v : Vector; p : Poly_Sys ) return Poly_Sys is

  -- DESCRIPTION :
  --   Multiplies the ith polynomial with the ith entry in v.

    res : Poly_Sys(p'range);

  begin
    for i in p'range loop
      res(i) := v(i)*p(i);
    end loop;
    return res;
  end Mul;

  function Mul ( m : Matrix; v : Vector ) return Matrix is

  -- DESCRIPTION :
  --   Multiplies the ith row of m with the ith entry in v.

    res : Matrix(m'range(1),m'range(2));

  begin
    for i in m'range(1) loop
      for j in m'range(2) loop
        res(i,j) := v(i)*m(i,j);
      end loop;
    end loop;
    return res;
  end Mul;

-- AUXILIARIES TO THE CONSTRUCTORS :

  function Linear_Start_Factor
             ( n,k : in natural; a : in Complex_Number ) return Poly is

  -- DESCRIPTION :
  --   Returns a*(1-t)^k, with t as the (n+1)-th variable in the polynomial.

    res,tmp : Poly;
    t : Term;

  begin
    t.cf := a;
    t.dg := new Standard_Natural_Vectors.Vector'(1..n+1 => 0);
    res := Create(t);                                           -- res = a
    t.cf := Create(1.0);
    tmp := Create(t);                                           -- tmp = 1
    t.dg(n+1) := 1;
    Sub(tmp,t);                                                 -- tmp = 1-t
    Standard_Natural_Vectors.Clear
      (Standard_Natural_Vectors.Link_to_Vector(t.dg));
    for i in 1..k loop
      Mul(res,tmp);
    end loop;
    Clear(tmp);
    return res;
  end Linear_Start_Factor;

  function Nonlinear_Start_Factor
             ( n,k : in natural; a : in Complex_Number ) return Poly is

  -- DESCRIPTION :
  --   Returns (1 - (t - t*(1-t)*a))^k, with t as the (n+1)-th variable
  --   in the polynomial.

    res,tmp : Poly;
    t : Term;

  begin
    t.cf := Create(1.0);
    t.dg := new Standard_Natural_Vectors.Vector'(1..n+1 => 0);
    res := Create(t);                                         -- res = 1
    tmp := Create(t);                                         -- tmp = 1
    t.dg(n+1) := 1;
    Sub(tmp,t);                                           -- tmp = 1-t
    Mul(tmp,t);                                           -- tmp = t*(1-t)
    Mul(tmp,-a);                                          -- tmp = -a*t*(1-t)
    Add(tmp,t);                                       -- tmp = t - a*t*(1-t)
    Sub(res,tmp);                                 -- res = 1 - t + a*t*(1-t)
    if k > 1
     then Copy(res,tmp);
          for i in 1..(k-1) loop
            Mul(res,tmp);
          end loop;
    end if;
    Standard_Natural_Vectors.Clear
      (Standard_Natural_Vectors.Link_to_Vector(t.dg));
    Clear(tmp);
    return res;
  end Nonlinear_Start_Factor;

  function Linear_Target_Factor ( n,k : in natural ) return Poly is

  -- DESCRIPTION :
  --   Returns t^k, with t as the (n+1)-th variable in the polynomial.

    res : Poly;
    t : Term;

  begin
    t.cf := Create(1.0);
    t.dg := new Standard_Natural_Vectors.Vector'(1..n+1 => 0);
    t.dg(n+1) := k;
    res := Create(t);
    Standard_Natural_Vectors.Clear
      (Standard_Natural_Vectors.Link_to_Vector(t.dg));
    return res;
  end Linear_Target_Factor;

  function Linear_Target_Factor
             ( n,k : in natural; a : in Complex_Number ) return Poly is

  -- DESCRIPTION :
  --   Returns a*t^k, with t as the (n+1)-th variable in the polynomial.

    res : Poly;
    t : Term;

  begin
    t.cf := a;
    t.dg := new Standard_Natural_Vectors.Vector'(1..n+1 => 0);
    t.dg(n+1) := k;
    res := Create(t);
    Standard_Natural_Vectors.Clear
      (Standard_Natural_Vectors.Link_to_Vector(t.dg));
    return res;
  end Linear_Target_Factor;

  function Nonlinear_Target_Factor
             ( n,k : in natural; a : in Complex_Number ) return Poly is

  -- DESCRIPTION :
  --   Returns (t - t*(1-t)*a)^k, with t as the (n+1)-th variable in
  --   the polynomial.

    res,tmp : Poly;
    t : Term;

  begin
    t.cf := Create(1.0);
    t.dg := new Standard_Natural_Vectors.Vector'(1..n+1 => 0);
    tmp := Create(t);                                         -- tmp = 1
    t.dg(n+1) := 1;
    res := Create(t);                                         -- res = t
    Sub(tmp,t);                                           -- tmp = 1-t
    Mul(tmp,t);                                           -- tmp = t*(1-t)
    Mul(tmp,a);                                           -- tmp = a*t*(1-t)
    Sub(res,tmp);                                     -- res = t - a*t*(1-t)
    if k > 1
     then Copy(res,tmp);
          for i in 1..(k-1) loop
            Mul(res,tmp);
          end loop;
    end if;
    Standard_Natural_Vectors.Clear
      (Standard_Natural_Vectors.Link_to_Vector(t.dg));
    Clear(tmp);
    return res;
  end Nonlinear_Target_Factor;

  function Plus_one_Unknown ( p : in Poly ) return Poly is

  -- DESCRIPTION :
  --   The returning polynomial has place for one additional unknown.

    res : Poly;

    procedure Plus_Unknown_In_Term ( t : in out Term; c : out boolean ) is

      n : constant natural := t.dg'length;
      temp : Standard_Natural_Vectors.Vector(1..(n+1));

    begin
      temp(1..n) := t.dg.all;
      temp(n+1) := 0;
      Standard_Natural_Vectors.Clear
        (Standard_Natural_Vectors.Link_to_Vector(t.dg));
      t.dg := new Standard_Natural_Vectors.Vector'(temp);
      c := true;
    end Plus_Unknown_In_Term;
    procedure Plus_Unknown_In_Terms is 
       new Changing_Iterator (process => Plus_Unknown_In_Term);

  begin
    Copy(p,res);
    Plus_Unknown_in_Terms(res);
    return res;
  end Plus_one_Unknown;

  function Linear_Homotopy
             ( p,q : in Poly_Sys; k : in positive; a : in Complex_Number )
             return Poly_Sys is

    h : Poly_Sys(p'range);
    tempp,tempq,q_fac,p_fac : Poly;
    n : natural := p'length;

  begin
    q_fac := Linear_Start_Factor(n,k,a);
    p_fac := Linear_Target_Factor(n,k);
    for i in h'range loop
      tempq := Plus_one_Unknown(q(i));
      tempp := Plus_one_Unknown(p(i));
      Mul(tempq,q_fac);
      Mul(tempp,p_fac);
      h(i) := tempq + tempp;
      Clear(tempq); Clear(tempp);
    end loop;
    Clear(q_fac); Clear(p_fac);
    return h;
  end Linear_Homotopy;

  function Linear_Homotopy
             ( p,q : in Poly_Sys; k : in positive; a : in Vector )
             return Poly_Sys is

    h : Poly_Sys(p'range);
    tempp,tempq,q_fac,p_fac : Poly;
    n : natural := p'length;

  begin
    for i in h'range loop
      q_fac := Linear_Start_Factor(n,k,a(i));
      p_fac := Linear_Target_Factor(n,k);
      tempq := Plus_one_Unknown(q(i));
      tempp := Plus_one_Unknown(p(i));
      Mul(tempq,q_fac);
      Mul(tempp,p_fac);
      h(i) := tempq + tempp;
      Clear(tempq); Clear(tempp);
      Clear(q_fac); Clear(p_fac);
    end loop;
    return h;
  end Linear_Homotopy;

  function Linear_Homotopy
             ( p,q : in Poly_Sys; k : in positive; a,b : in Vector )
             return Poly_Sys is

    h : Poly_Sys(p'range);
    tempp,tempq,q_fac,p_fac : Poly;
    n : natural := p'length;

  begin
    for i in h'range loop
      q_fac := Linear_Start_Factor(n,k,a(i));
      p_fac := Linear_Target_Factor(n,k,b(i));
      tempq := Plus_one_Unknown(q(i));
      tempp := Plus_one_Unknown(p(i));
      Mul(tempq,q_fac);
      Mul(tempp,p_fac);
      h(i) := tempq + tempp;
      Clear(tempq); Clear(tempp);
      Clear(q_fac); Clear(p_fac);
    end loop;
    return h;
  end Linear_Homotopy;

  function Nonlinear_Homotopy
             ( p,q : in Poly_Sys; k : in positive; a,b : in Vector )
             return Poly_Sys is

    h : Poly_Sys(p'range);
    tempp,tempq,q_fac,p_fac : Poly;
    n : natural := p'length;

  begin
    for i in h'range loop
      q_fac := Nonlinear_Start_Factor(n,k,a(i));
      p_fac := Nonlinear_Target_Factor(n,k,b(i));
      tempq := Plus_one_Unknown(q(i));
      tempp := Plus_one_Unknown(p(i));
      Mul(tempq,q_fac);
      Mul(tempp,p_fac);
      h(i) := tempq + tempp;
      Clear(tempq); Clear(tempp);
      Clear(q_fac); Clear(p_fac);
    end loop;
    return h;
  end Nonlinear_Homotopy;

  procedure Create ( p,q : in Poly_Sys; k : in positive;
                     a : in Complex_Number ) is

    n : constant natural := p'length;
    dp, dq : Jaco_Mat(1..n,1..n);
    ho : homdata(art,n,n+1);

  begin
    Copy(p,ho.p); Copy(q,ho.q);
    ho.h := Linear_Homotopy(p,q,k,a);
    ho.pe := Create(ho.p);
    ho.qe := Create(ho.q);
    ho.he := Create(ho.h);
    dp := Create(ho.p);
    dq := Create(ho.q);
    ho.dh := Create(ho.h);
    ho.dpe := Create(dp);
    ho.dqe := Create(dq);
    ho.dhe := Create(ho.dh);
    Clear(dp); Clear(dq);
    ho.k := k;
    for i in 1..n loop
      ho.gamma(i) := a;
    end loop;
    for i in 1..n loop
      ho.beta(i) := Create(1.0);
    end loop;
    ho.linear := true;
    hom := new homdata'(ho);
  end Create; 

  procedure Create ( p,q : in Poly_Sys; k : in positive; a : in Vector ) is

    n : constant natural := p'length;
    dp, dq : Jaco_Mat(1..n,1..n);
    ho : homdata(art,n,n+1);

  begin
    Copy(p,ho.p); Copy(q,ho.q);
    ho.h := Linear_Homotopy(p,q,k,a);
    ho.pe := Create(ho.p);
    ho.qe := Create(ho.q);
    ho.he := Create(ho.h);
    dp := Create(ho.p);
    dq := Create(ho.q);
    ho.dh := Create(ho.h);
    ho.dpe := Create(dp);
    ho.dqe := Create(dq);
    ho.dhe := Create(ho.dh);
    Clear(dp); Clear(dq);
    ho.k := k;
    ho.gamma := a;
    for i in 1..n loop
      ho.beta(i) := Create(1.0);
    end loop;
    ho.linear := true;
    hom := new homdata'(ho);
  end Create;

  procedure Create ( p,q : in Poly_Sys; k : in positive; a,b : in Vector;
                     linear : in boolean ) is

    n : constant natural := p'length;
    dp, dq : Jaco_Mat(1..n,1..n);
    ho : homdata(art,n,n+1);

  begin
    Copy(p,ho.p); Copy(q,ho.q);
    ho.linear := linear;
    if linear
     then ho.h := Linear_Homotopy(p,q,k,a,b);
     else ho.h := Nonlinear_Homotopy(p,q,k,a,b);      
    end if;
    ho.pe := Create(ho.p);
    ho.qe := Create(ho.q);
    ho.he := Create(ho.h);
    dp := Create(ho.p);
    dq := Create(ho.q);
    ho.dh := Create(ho.h);
    ho.dpe := Create(dp);
    ho.dqe := Create(dq);
    ho.dhe := Create(ho.dh);
    Clear(dp); Clear(dq);
    ho.k := k;
    ho.gamma := a;
    ho.beta := b;
    hom := new homdata'(ho);
  end Create;

  procedure Create ( p : in Poly_Sys; k : in integer ) is

    n : constant natural := p'last-p'first+1; 
    ho : homdata(nat,n,n+1);

  begin
    Copy(p,ho.p);
    ho.pe := Create(ho.p);
    ho.dh := Create(ho.p);
    ho.dhe := Create(ho.dh);
    ho.i := k;
    hom := new homdata'(ho);
  end Create;

-- SELECTOR :

  function Homotopy_System return Poly_Sys is

    ho : homdata renames hom.all;

  begin
    case ho.ht is
      when nat => return ho.p;
      when art => return ho.h;
    end case;
  end Homotopy_System;

-- SYMBOLIC ROUTINES :

  function Eval ( t : Complex_Number ) return Poly_Sys is

    ho : homdata renames hom.all;

  begin
    case ho.ht is
      when nat =>  -- t = x(ho.i)
        return Eval(ho.p,t,ho.i);
      when art =>  -- compute : a * ((1 - t)^k) * q + (t^k) * p
        declare
          p_factor,q_factor : Vector(1..ho.n);
          res,tmp : Poly_Sys(1..ho.n);
        begin
          if ho.linear 
           then
             if AbsVal(t) = 0.0
              then return Mul(ho.gamma,ho.q);
              elsif abs(REAL_PART(t) - 1.0 ) + 1.0 = 1.0 
                   and then abs(IMAG_PART(t)) + 1.0 = 1.0
                  then return Mul(ho.beta,ho.p);
                  else for i in 1..ho.n loop
                         q_factor(i) := ho.gamma(i);
                         p_factor(i) := ho.beta(i);
                         for i in 1..ho.k loop
                           q_factor(i) := (Create(1.0)-t) * q_factor(i);
                           p_factor(i) := t * p_factor(i);
                         end loop;
                       end loop;
                       res := Mul(p_factor,ho.p); 
                       tmp := Mul(q_factor,ho.q);
                       Add(res,tmp);
                       Clear(tmp);
                       return res;
             end if;
           else return res;  -- still to do !
          end if;
        end;
    end case;
  end Eval;

  function Diff ( t : Complex_Number ) return Poly_Sys is

    ho : homdata renames hom.all;

  begin
    case ho.ht is
      when nat =>  -- t = x(ho.i)
        return Diff(ho.p,ho.i);
      when art =>  -- compute  - a*k*(1 - t)^(k-1)*q + b*k*t^(k-1)*p
        declare
          q_factor,p_factor : Vector(1..ho.n);
          tmp : Poly_Sys(1..ho.n);
          res : Poly_Sys(1..ho.n);
        begin
          if ho.linear
           then
             for i in 1..ho.n loop
               q_factor(i) := (-ho.gamma(i)) * Create(double_float(ho.k));
               p_factor(i) := ho.beta(i) * Create(double_float(ho.k));
             end loop;
             if AbsVal(t) = 0.0
              then if ho.k = 1
                    then res := Mul(p_factor,ho.p);
                         tmp := Mul(q_factor,ho.q);
                         Add(res,tmp);
                         Clear(tmp);
                         return res;
                    else return Mul(q_factor,ho.q);
                   end if;
              elsif abs( REAL_PART(t) - 1.0 ) + 1.0 = 1.0
                   and then abs(IMAG_PART(t)) + 1.0 = 1.0
                  then return Create(double_float(ho.k)) * ho.p;
                  else for i in 1..(ho.k-1) loop
                         q_factor := (Create(1.0)-t) * q_factor;
                         p_factor := t * p_factor;
                       end loop;
                       res := Mul(p_factor,ho.p);
                       tmp := Mul(q_factor,ho.q);
                       Add(res,tmp);
                       Clear(tmp);
                       return res;
             end if;
           else return res;   -- still left to do !!!
          end if;
        end;
    end case;
  end Diff;

-- NUMERIC ROUTINES :

  function Eval ( x : Vector; t : Complex_Number ) return Vector is

    ho : homdata renames hom.all;

  begin
    case ho.ht is
      when nat =>
        declare
          y : Vector(x'first..x'last+1);
        begin
          y(1..ho.i-1) := x(1..ho.i-1);
          y(ho.i) := t;
          y(ho.i+1..y'last) := x(ho.i..x'last);
          return Eval(ho.pe,y);
        end;
      when art =>
        if AbsVal(t) + 1.0 = 1.0
         then if ho.linear
               then return ho.gamma * Eval(ho.qe,x);
               else return Eval(ho.qe,x);
              end if;
         elsif abs( REAL_PART(t) - 1.0 ) + 1.0 = 1.0
              and then abs(IMAG_PART(t)) + 1.0 = 1.0
             then if ho.linear
                   then return ho.beta * Eval(ho.pe,x);
                   else return Eval(ho.pe,x);
                  end if;
             else declare
                    y : Vector(x'first..x'last+1);
                  begin
                    y(x'range) := x;
                    y(y'last) := t;
                    return Eval(ho.he,y);
                  end;
        end if;
    end case;
  end Eval;

  function Diff ( x : Vector; t : Complex_Number ) return Vector is

    n : natural := x'length;

  begin
    case hom.ht is
      when nat => return Diff(x,t,hom.i);
      when art => return Diff(x,t,n+1);
    end case;
  end Diff;
 
  function Diff ( x : Vector; t : Complex_Number ) return Matrix is

    ho : homdata renames hom.all;
    n : natural renames ho.n;

  begin
    case ho.ht is
      when nat =>
        declare
          m : Matrix(1..n,1..n);
          y : Vector(1..n+1);
        begin
          y(1..ho.i-1) := x(1..ho.i-1);
          y(ho.i) := t;
          y(ho.i+1..n+1) := x(ho.i..n);
          for i in 1..n loop
            for j in 1..n loop
              m(i,j) := Eval(ho.dhe(i,j),y);
            end loop;
          end loop;
          return m;
        end;
      when art =>
        if AbsVal(t) + 1.0 = 1.0
         then declare
                m : Matrix(1..n,1..n) := Eval(ho.dqe,x);
              begin
                if ho.linear
                 then return Mul(m,ho.gamma);
                 else return m;
                end if;
              end;
         elsif abs( REAL_PART(t) - 1.0 ) + 1.0 = 1.0
              and then abs(IMAG_PART(t)) + 1.0 = 1.0
             then declare
                    m : Matrix(1..n,1..n) := Eval(ho.dpe,x);
                  begin
                    if ho.linear
                     then return Mul(m,ho.beta);
                     else return m;
                    end if;
                  end;
             else declare
                    m : Matrix(1..n,1..n);
                    y : Vector(1..n+1);
                  begin
                    y(1..n) := x;
                    y(n+1) := t;
                    for i in 1..n loop
                      for j in 1..n loop
                        m(i,j) := Eval(ho.dhe(i,j),y);
                      end loop;
                    end loop;
                    return m;
                  end;
        end if;
    end case;
  end Diff;

  function Diff ( x : Vector; t : Complex_Number; k : natural ) return Vector is

    ho : homdata renames hom.all;
    n : natural renames ho.n;
    y : Vector(1..n+1);
    res : Vector(1..n);

  begin
    case ho.ht is
      when nat => y(1..ho.k-1) := x(1..ho.i-1);
                  y(ho.i) := t;
                  y(ho.i+1..n+1) := x(ho.i..n);
      when art => y(1..n) := x;
                  y(n+1) := t;
    end case;
    for i in 1..n loop
      res(i) := Eval(ho.dhe(i,k),y);
    end loop;
    return res;
  end Diff;

  function Diff ( x : Vector; t : Complex_Number; k : natural ) return Matrix is

    ho : homdata renames hom.all;
    n : natural renames ho.n;
    y : Vector(1..n+1);
    res : Matrix(1..n,1..n);

  begin
    case ho.ht is
      when nat => y(1..ho.i-1) := x(1..ho.i-1);
                  y(ho.i) := t;
                  y(ho.i+1..n+1) := x(ho.i..n); 
      when art => y(1..n) := x;
                  y(n+1) := t;
    end case;
    for j in 1..(k-1) loop
      for i in 1..n loop
        res(i,j) := Eval(ho.dhe(i,j),y);
      end loop;
    end loop;
    for j in (k+1)..(n+1) loop
      for i in 1..n loop
        res(i,j-1) := Eval(ho.dhe(i,j),y);
      end loop;
    end loop;
    return res;
  end Diff;

-- DESTRUCTOR :

  procedure free is new unchecked_deallocation (homdata,homtp);

  procedure Clear is
  begin
    if hom /= null
     then declare
            ho : homdata renames hom.all;
          begin
            Clear(ho.p);  Clear(ho.pe);
            Clear(ho.dh); Clear(ho.dhe);
            case ho.ht is
              when nat => null;
              when art =>
                Clear(ho.q);   Clear(ho.qe);
                Clear(ho.h);   Clear(ho.he);
                Clear(ho.dpe); Clear(ho.dqe);
            end case;
          end;
          free(hom);
    end if;
  end Clear;

end Homotopy;