File: [local] / OpenXM_contrib / PHC / Ada / Schubert / driver_for_quantum_pieri.adb (download)
Revision 1.1.1.1 (vendor branch), Sun Oct 29 17:45:32 2000 UTC (23 years, 8 months ago) by maekawa
Branch: PHC, MAIN
CVS Tags: v2, maekawa-ipv6, RELEASE_1_2_3, RELEASE_1_2_2_KNOPPIX_b, RELEASE_1_2_2_KNOPPIX, RELEASE_1_2_2, RELEASE_1_2_1, HEAD Changes since 1.1: +0 -0
lines
Import the second public release of PHCpack.
OKed by Jan Verschelde.
|
with text_io,integer_io; use text_io,integer_io;
with Communications_with_User; use Communications_with_User;
with Timing_Package; use Timing_Package;
with Standard_Floating_Numbers; use Standard_Floating_Numbers;
with Standard_Complex_Numbers; use Standard_Complex_Numbers;
with Standard_Floating_Vectors;
with Standard_Natural_Vectors;
with Standard_Complex_Vectors;
with Standard_Complex_Norms_Equals; use Standard_Complex_Norms_Equals;
with Standard_Natural_Matrices;
with Standard_Natural_Matrices_io; use Standard_Natural_Matrices_io;
with Standard_Complex_Matrices; use Standard_Complex_Matrices;
with Standard_Random_Vectors; use Standard_Random_Vectors;
with Standard_Complex_VecMats; use Standard_Complex_VecMats;
with Symbol_Table; use Symbol_Table;
with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
with Standard_Complex_Poly_Matrices;
with Standard_Complex_Poly_Matrices_io; use Standard_Complex_Poly_Matrices_io;
with Standard_Complex_Poly_Systems; use Standard_Complex_Poly_Systems;
with Standard_Complex_Poly_Systems_io; use Standard_Complex_Poly_Systems_io;
with Standard_Complex_Solutions; use Standard_Complex_Solutions;
with Homotopy;
with Drivers_for_Poly_Continuation; use Drivers_for_Poly_Continuation;
with Increment_and_Fix_Continuation; use Increment_and_Fix_Continuation;
with Standard_Root_Refiners; use Standard_Root_Refiners;
with Brackets,Brackets_io; use Brackets,Brackets_io;
with Localization_Posets; use Localization_Posets;
with Localization_Posets_io; use Localization_Posets_io;
with Curves_into_Grassmannian; use Curves_into_Grassmannian;
with Curves_into_Grassmannian_io; use Curves_into_Grassmannian_io;
with Deformation_Posets; use Deformation_Posets;
with Determinantal_Systems; use Determinantal_Systems;
with Plane_Representations; use Plane_Representations;
with Drivers_for_Input_Planes; use Drivers_for_Input_planes;
procedure Driver_for_Quantum_Pieri
( file : in file_type; n,d,q : in natural ) is
function Solution_Plane
( top,bottom : Bracket;
locmap : Standard_Natural_Matrices.Matrix;
mat : Standard_Complex_Matrices.Matrix )
return Solution is
-- DESCRIPTION :
-- Returns the representation of the solution plane as a vector.
solloc : constant Standard_Complex_Vectors.Vector
:= Column_Vector_Rep(locmap,Localize(locmap,mat));
sol : Solution(solloc'length);
begin
sol.m := 1;
sol.t := Create(0.0);
sol.res := 0.0;
sol.err := 0.0;
sol.rco := 0.0;
sol.v := solloc;
return sol;
end Solution_Plane;
function Solution_Planes ( top,bottom : Bracket;
locmap : Standard_Natural_Matrices.Matrix;
vm : VecMat ) return Solution_List is
-- DESCRIPTION :
-- Returns the representation of the vector of planes as a solution list.
res,res_last : Solution_List;
begin
for i in vm'range loop
Append(res,res_last,Solution_Plane(top,bottom,locmap,vm(i).all));
end loop;
return res;
end Solution_Planes;
function Create_Polynomial_System
( top,bottom : Bracket;
locmap : Standard_Natural_Matrices.Matrix;
xpm : Standard_Complex_Poly_Matrices.Matrix;
svals : Standard_Complex_Vectors.Vector;
planes : VecMat ) return Poly_Sys is
-- DESCRIPTION :
-- Returns the polynomial system that collects the intersection
-- conditions for meeting the given m-planes at the specified s-values.
-- The system is localized according to the given localization map.
res,wrksys : Poly_Sys(svals'range);
begin
for i in svals'range loop
declare
eva : Standard_Complex_Poly_Matrices.Matrix(xpm'range(1),xpm'range(2))
:= Elim(xpm,svals(i),Create(1.0));
wrk : constant Poly_Sys := Polynomial_Equations(planes(i).all,eva);
begin
wrksys(i) := wrk(wrk'first);
Standard_Complex_Poly_Matrices.Clear(eva);
end;
end loop;
res := Column_Localize(top,bottom,locmap,wrksys);
Clear(wrksys);
return res;
end Create_Polynomial_System;
procedure Refine_Roots ( file : in file_type;
p : in Poly_Sys; sols : in out Solution_List ) is
-- DESCRIPTION :
-- Calls the root refiner, for square polynomial systems only.
epsxa : constant double_float := 10.0**(-8);
epsfa : constant double_float := 10.0**(-8);
tolsing : constant double_float := 10.0**(-8);
max : constant natural := 3;
numit : natural := 0;
begin
Reporting_Root_Refiner(file,p,sols,epsxa,epsfa,tolsing,numit,max,false);
end Refine_Roots;
procedure Solve_Target_System
( file : in file_type; start,target : in Poly_Sys;
sols : in out Solution_List; report : in boolean ) is
-- DESCRIPTION :
-- Calls the standard continuation routines to solve a specific
-- target system, starting at the solutions of the start system.
-- REQUIRED : not Is_Null(sols).
timer : Timing_Widget;
n : constant natural := target'last;
a : Standard_Complex_Vectors.Vector(1..n) := Random_Vector(1,n);
b : Standard_Complex_Vectors.Vector(1..n) := Random_Vector(1,n);
begin
tstart(timer);
Homotopy.Create(target,start,1,a,b,true); -- linear cheater
Set_Continuation_Parameter(sols,Create(0.0));
declare
procedure Sil_Cont is
new Silent_Continue(Max_Norm,
Homotopy.Eval,Homotopy.Diff,Homotopy.Diff);
procedure Rep_Cont is
new Reporting_Continue(Max_Norm,
Homotopy.Eval,Homotopy.Diff,Homotopy.Diff);
begin
if report
then Rep_Cont(file,sols,false,Create(1.0));
else Sil_Cont(sols,false,Create(1.0));
end if;
end;
tstop(timer);
new_line(file);
print_times(file,timer,"Cheater's homotopy to target system");
Refine_Roots(file,target,sols);
end Solve_Target_System;
procedure Solve_Hypersurface_Target_System
( file : in file_type; m,p,q : in natural;
start_svals,target_svals
: in Standard_Complex_Vectors.Vector;
start_planes,target_planes : in VecMat;
index_poset : in Array_of_Array_of_Nodes;
deform_poset : in Array_of_Array_of_VecMats;
report : in boolean ) is
-- DESCRIPTION :
-- This procedure tests the output of the deformation poset,
-- creating polynomial representations of the intersection conditions
-- and solution lists representing the solution planes.
-- ON ENTRY :
-- file to write intermediate output on;
-- m dimension of the input planes;
-- p dimension of the solution planes;
-- q degree of the maps;
-- start_svals interpolation points at the start;
-- target_svals interpolation points at the target;
-- start_planes input m-planes in general position;
-- target_planes specific input m-planes;
-- index_poset indexed localization poset;
-- deform_poset poset with the solution p-planes;
-- report switch for intermediate output during continuation.
dim : constant natural := m*p+q*(m+p);
top : constant Bracket := index_poset(dim)(1).top;
bot : constant Bracket := index_poset(dim)(1).bottom;
xpm : Standard_Complex_Poly_Matrices.Matrix(1..m+p,1..p)
:= Symbolic_Create(m,p,q,top,bot);
solplanes : constant VecMat := deform_poset(dim)(1).all;
locmap : Standard_Natural_Matrices.Matrix(1..(m+p)*(q+1),1..p)
:= Standard_Coordinate_Frame(m,p,q,top,bot,solplanes(1).all);
locsys : Poly_Sys(start_planes'range)
:= Create_Polynomial_System
(top,bot,locmap,xpm,start_svals,start_planes);
target : Poly_Sys(target_planes'range)
:= Create_Polynomial_System
(top,bot,locmap,xpm,target_svals,target_planes);
sols : Solution_List := Solution_Planes(top,bot,locmap,solplanes);
begin
One_Set_up_Symbol_Table(m,p,q,top,bot);
new_line(file);
put(file,"The "); put(file,q,1); put(file,"-map of "); put(file,p,1);
put_line(file,"-planes representation : ");
put(file,xpm);
put_line(file,"with as localization map :"); put(file,locmap);
new_line(file);
Reduce_Symbols(top,bot,locmap);
put_line(file,"THE GENERIC SYSTEM : ");
put_line(file,locsys);
new_line(file);
Refine_Roots(file,locsys,sols);
new_line(file);
put_line(file,"THE TARGET SYSTEM : ");
put_line(file,target);
new_line(file);
Solve_Target_System(file,locsys,target,sols,report);
end Solve_Hypersurface_Target_System;
procedure Set_Parameters ( file : in file_type; report : out boolean ) is
-- DESCRIPTION :
-- Interactive determination of the continuation and output parameters.
oc : natural;
begin
new_line;
Driver_for_Continuation_Parameters(file);
new_line;
Driver_for_Process_io(file,oc);
report := not (oc = 0);
new_line;
put_line("No more input expected. See output file for results...");
new_line;
new_line(file);
end Set_Parameters;
procedure Write_Poset_Times
( file : in file_type; timer : in Timing_Widget;
npaths : in Standard_Natural_Vectors.Vector;
timings : in Duration_Array ) is
-- DESCRIPTION :
-- Writes a overview of #paths and timings spent during deformations
-- along the poset structure.
begin
new_line(file);
put_line(file,"--------------------------------------");
put_line(file,"| TIMING INFORMATION OVERVIEW |");
put_line(file,"--------------------------------------");
put_line(file,"| n | #paths | user cpu time |");
put_line(file,"--------------------------------------");
for i in npaths'range loop
if npaths(i) /= 0
then put(file,"|"); put(file,i,4); put(file," |");
put(file,npaths(i),7); put(file," | ");
print_hms(file,timings(i)); put(file," |"); new_line(file);
end if;
end loop;
put_line(file,"--------------------------------------");
put(file,"| total |");
put(file,Standard_Natural_Vectors.Sum(npaths),7);
put(file," | ");
print_hms(file,Elapsed_User_Time(timer));
put(file," |"); new_line(file);
put_line(file,"--------------------------------------");
end Write_Poset_Times;
procedure Solve_Deformation_Poset
( file : in file_type; m,p,q : in natural;
index_poset : in out Array_of_Array_of_Nodes ) is
-- DESCRIPTION :
-- Writes the symbolic form of the maps.
lnd : Link_to_Node;
deform_poset : Array_of_Array_of_VecMats(index_poset'range)
:= Create(index_poset);
dim : constant natural := (m*p)+q*(m+p);
input : VecMat(1..dim) := Random_Complex_Planes(m,p,q);
svals : Standard_Complex_Vectors.Vector(1..dim) := Random_Vector(1,dim);
target_planes : VecMat(1..dim);
target_svals : Standard_Floating_Vectors.Vector(1..dim);
comp_target_svals : Standard_Complex_Vectors.Vector(1..dim);
root : Node := index_poset(dim)(1).all;
ans : character;
report,outlog : boolean;
timer : Timing_Widget;
npaths : Standard_Natural_Vectors.Vector(1..dim) := (1..dim => 0);
timings : Duration_Array(1..dim) := (1..dim => 0.0);
begin
new_line;
put("Do you want to have the homotopies on file ? (y/n) ");
Ask_Yes_or_No(ans);
outlog := (ans = 'y');
Driver_for_Input_Planes(file,m,p,q,target_svals,target_planes);
Set_Parameters(file,report);
tstart(timer);
Solve(file,m+p,q,deform_poset,root,input,svals,report,outlog,
npaths,timings);
tstop(timer);
new_line(file);
print_times(file,timer,"Solving along the deformation poset");
Write_Poset_Times(file,timer,npaths,timings);
for i in target_svals'range loop
comp_target_svals(i) := Create(target_svals(i));
end loop;
Solve_Hypersurface_Target_System
(file,m,p,q,svals,comp_target_svals,input,target_planes,
index_poset,deform_poset,report);
end Solve_Deformation_Poset;
procedure Create_Hypersurface_Localization_Poset
( file : in file_type;
lnkroot : in Link_to_Node; m,p,q : in natural ) is
-- DESCRIPTION :
-- Creates the posets and outputs them to the screen and on file.
-- Calls the solver afterwards.
nq : constant natural := m*p + q*(m+p);
level_poset : Array_of_Nodes(0..nq);
index_poset : Array_of_Array_of_Nodes(0..nq);
nbp : natural;
begin
level_poset := Create_Leveled_Poset(lnkroot);
Count_Roots(level_poset);
index_poset := Create_Indexed_Poset(level_poset);
put(index_poset);
put(file,index_poset);
put_line("The size of the poset : "); put_roco(index_poset);
put_line(file,"The size of the poset : "); put_roco(file,index_poset);
nbp := Root_Count_Sum(level_poset);
put("The number of paths : "); put(nbp,1); new_line;
put(file,"The number of paths : "); put(file,nbp,1); new_line(file);
Solve_Deformation_Poset(file,m,p,q,index_poset);
end Create_Hypersurface_Localization_Poset;
procedure Create_General_Localization_Poset
( file : in file_type; lnkroot : in Link_to_Node;
m,p,q : in natural; codim : in Bracket ) is
-- DESCRIPTION :
-- Creates the posets and outputs them to the screen and on file.
nq : constant natural := m*p + q*(m+p);
level_poset : Array_of_Nodes(0..nq);
index_poset : Array_of_Array_of_Nodes(0..nq);
nbp : natural;
begin
level_poset := Create_Leveled_Poset(lnkroot);
Count_Roots(level_poset);
index_poset := Create_Indexed_Poset(level_poset);
put(index_poset);
put(file,index_poset);
put_line("The size of the poset : "); put_roco(index_poset);
put_line(file,"The size of the poset : "); put_roco(file,index_poset);
nbp := Root_Count_Sum(level_poset);
put("The number of paths : "); put(nbp,1); new_line;
put(file,"The number of paths : "); put(file,nbp,1); new_line(file);
-- Solve_Deformation_Poset(file,m,p,q,index_poset);
end Create_General_Localization_Poset;
procedure Create_Top_Hypersurface_Poset
( file : in file_type; m,p,q : in natural ) is
-- DESCRIPTION :
-- Creates the poset by incrementing only top pivots.
timer : Timing_Widget;
root : Node(p) := Trivial_Root(m,p,q);
lnkroot : Link_to_Node := new Node'(root);
begin
tstart(timer);
Q_Top_Create(lnkroot,root.bottom(p),m+p);
put_line("The poset created from the top : ");
put_line(file,"The poset created from the top : ");
Create_Hypersurface_Localization_Poset(file,lnkroot,m,p,q);
tstop(timer);
new_line(file);
print_times(file,timer,"Total time for Quantum Pieri Homotopy Algorithm");
end Create_Top_Hypersurface_Poset;
procedure Create_Bottom_Hypersurface_Poset
( file : in file_type; m,p,q : in natural ) is
-- DESCRIPTION :
-- Creates the poset by decrementing only bottom pivots.
timer : Timing_Widget;
root : Node(p) := Trivial_Root(m,p,q);
lnkroot : Link_to_Node := new Node'(root);
begin
tstart(timer);
Q_Bottom_Create(lnkroot,m+p);
put_line("The poset created from the bottom : ");
put_line(file,"The poset created from the bottom : ");
Create_Hypersurface_Localization_Poset(file,lnkroot,m,p,q);
tstop(timer);
new_line(file);
print_times(file,timer,"Total time for Quantum Pieri Homotopy Algorithm");
end Create_Bottom_Hypersurface_Poset;
procedure Create_Mixed_Hypersurface_Poset
( file : in file_type; m,p,q : in natural ) is
-- DESCRIPTION :
-- Creates the poset by incrementing top and decrementing bottom pivots.
timer : Timing_Widget;
root : Node(p) := Trivial_Root(m,p,q);
lnkroot : Link_to_Node := new Node'(root);
begin
tstart(timer);
Q_Top_Bottom_Create(lnkroot,root.bottom(p),m+p);
put_line("The poset created in a mixed fashion : ");
put_line(file,"The poset created in a mixed fashion :");
Create_Hypersurface_Localization_Poset(file,lnkroot,m,p,q);
tstop(timer);
new_line(file);
print_times(file,timer,"Total time for Quantum Pieri Homotopy Algorithm");
end Create_Mixed_Hypersurface_Poset;
procedure Create_Top_General_Poset
( file : in file_type; m,p,q : in natural ) is
-- DESCRIPTION :
-- Creates the poset by incrementing top pivots.
timer : Timing_Widget;
root : Node(p) := Trivial_Root(m,p,q);
lnkroot : Link_to_Node := new Node'(root);
codim : constant Bracket := Read_Codimensions(m,p,q);
begin
tstart(timer);
Q_Top_Create(lnkroot,codim,root.bottom(p),m+p);
put_line("The poset created from the top : ");
put_line(file,"The poset created from the top :");
Create_General_Localization_Poset(file,lnkroot,m,p,q,codim);
tstop(timer);
new_line(file);
print_times(file,timer,"Total time for Quantum Pieri Homotopy Algorithm");
end Create_Top_General_Poset;
procedure Create_Bottom_General_Poset
( file : in file_type; m,p,q : in natural ) is
-- DESCRIPTION :
-- Creates the poset by decrementing bottom pivots.
timer : Timing_Widget;
root : Node(p) := Trivial_Root(m,p,q);
lnkroot : Link_to_Node := new Node'(root);
codim : constant Bracket := Read_Codimensions(m,p,q);
begin
tstart(timer);
Q_Bottom_Create(lnkroot,codim,m+p);
put_line("The poset created from the bottom : ");
put_line(file,"The poset created from the bottom :");
Create_General_Localization_Poset(file,lnkroot,m,p,q,codim);
tstop(timer);
new_line(file);
print_times(file,timer,"Total time for Quantum Pieri Homotopy Algorithm");
end Create_Bottom_General_Poset;
procedure Create_Mixed_General_Poset
( file : in file_type; m,p,q : in natural ) is
-- DESCRIPTION :
-- Creates the poset by incrementing top and decrementing bottom pivots.
timer : Timing_Widget;
root : Node(p) := Trivial_Root(m,p,q);
lnkroot : Link_to_Node := new Node'(root);
codim : constant Bracket := Read_Codimensions(m,p,q);
begin
tstart(timer);
Q_Top_Bottom_Create(lnkroot,codim,root.bottom(p),m+p);
put_line("The poset created from the bottom : ");
put_line(file,"The poset created from the bottom :");
Create_General_Localization_Poset(file,lnkroot,m,p,q,codim);
tstop(timer);
new_line(file);
print_times(file,timer,"Total time for Quantum Pieri Homotopy Algorithm");
end Create_Mixed_General_Poset;
procedure Main is
p : constant natural := d;
m : constant natural := n-d;
ans : character;
begin
new_line;
put_line("MENU for interpolating maps of fixed degree in Grassmannian.");
put_line(" 1. k_i = 1 consistently incrementing top pivots.");
put_line(" 2. consistently decrementing bottom pivots.");
put_line(" 3. mixed top-bottom sequence for poset creation.");
put_line(" 4. k_i >= 1 consistently incrementing top pivots.");
put_line(" 5. consistently incrementing bottom pivots.");
put_line(" 6. mixed top-bottom sequence for poset creation.");
put("Type 1, 2, 3, 4, 5, or 6 to choose : ");
Ask_Alternative(ans,"123456");
new_line;
case ans is
when '1' => Create_Top_Hypersurface_Poset(file,m,p,q);
when '2' => Create_Bottom_Hypersurface_Poset(file,m,p,q);
when '3' => Create_Mixed_Hypersurface_Poset(file,m,p,q);
when '4' => Create_Top_General_Poset(file,m,p,q);
when '5' => Create_Bottom_General_Poset(file,m,p,q);
when '6' => Create_Mixed_General_Poset(file,m,p,q);
when others => put_line("Option not recognized. Please try again.");
end case;
end Main;
begin
Main;
end Driver_for_Quantum_Pieri;