Annotation of OpenXM_contrib/PHC/Ada/Schubert/driver_for_sagbi_homotopies.adb, Revision 1.1
1.1 ! maekawa 1: with integer_io; use integer_io;
! 2: with Timing_Package; use Timing_Package;
! 3: with Communications_with_User; use Communications_with_User;
! 4: with Numbers_io; use Numbers_io;
! 5: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
! 6: with Standard_Floating_Numbers_io; use Standard_Floating_Numbers_io;
! 7: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
! 8: with Standard_Random_Numbers; use Standard_Random_Numbers;
! 9: with Standard_Integer_Vectors;
! 10: with Standard_Natural_Matrices;
! 11: with Standard_Natural_Matrices_io; use Standard_Natural_Matrices_io;
! 12: with Standard_Complex_Vectors;
! 13: with Standard_Complex_VecVecs;
! 14: with Standard_Random_Vectors; use Standard_Random_Vectors;
! 15: with Standard_Floating_Vectors;
! 16: with Standard_Floating_Vectors_io; use Standard_Floating_Vectors_io;
! 17: with Standard_Floating_Matrices;
! 18: with Standard_Floating_Matrices_io; use Standard_Floating_Matrices_io;
! 19: with Standard_Complex_Matrices;
! 20: with Standard_Complex_Norms_Equals; use Standard_Complex_Norms_Equals;
! 21: with Standard_Random_Matrices; use Standard_Random_Matrices;
! 22: with Symbol_Table; use Symbol_Table;
! 23: with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
! 24: with Standard_Complex_Polynomials_io; use Standard_Complex_Polynomials_io;
! 25: with Standard_Complex_Poly_Functions; use Standard_Complex_Poly_Functions;
! 26: with Standard_Complex_Poly_Systems; use Standard_Complex_Poly_Systems;
! 27: with Standard_Complex_Poly_Systems_io; use Standard_Complex_Poly_Systems_io;
! 28: with Standard_Complex_Poly_SysFun; use Standard_Complex_Poly_SysFun;
! 29: with Standard_Complex_Jaco_Matrices;
! 30: with Standard_Complex_Laur_Systems; use Standard_Complex_Laur_Systems;
! 31: with Standard_Complex_Laur_Functions; use Standard_Complex_Laur_Functions;
! 32: with Standard_Complex_Laur_Systems; use Standard_Complex_Laur_Systems;
! 33: with Standard_Complex_Laur_SysFun; use Standard_Complex_Laur_SysFun;
! 34: with Standard_Complex_Laur_Jacomats;
! 35: with Exponent_Vectors; use Exponent_Vectors;
! 36: with Standard_Poly_Laur_Convertors; use Standard_Poly_Laur_Convertors;
! 37: with Matrix_Indeterminates;
! 38: with Homotopy;
! 39: with Standard_Complex_Solutions; use Standard_Complex_Solutions;
! 40: with Standard_Complex_Solutions_io; use Standard_Complex_Solutions_io;
! 41: with Lists_of_Integer_Vectors; use Lists_of_Integer_Vectors;
! 42: with Lists_of_Integer_Vectors_io; use Lists_of_Integer_Vectors_io;
! 43: with Arrays_of_Integer_Vector_Lists; use Arrays_of_Integer_Vector_Lists;
! 44: with Increment_and_Fix_Continuation; use Increment_and_Fix_Continuation;
! 45: with Standard_Root_Refiners; use Standard_Root_Refiners;
! 46: with Drivers_for_Poly_Continuation; use Drivers_for_Poly_Continuation;
! 47: with Power_Lists; use Power_Lists;
! 48: with Integer_Lifting_Utilities; use Integer_Lifting_Utilities;
! 49: with Integer_Mixed_Subdivisions; use Integer_Mixed_Subdivisions;
! 50: with Integer_Polyhedral_Continuation; use Integer_Polyhedral_Continuation;
! 51: with Triangulations; use Triangulations;
! 52: with Dynamic_Triangulations; use Dynamic_Triangulations;
! 53: with Triangulations_and_Subdivisions; use Triangulations_and_Subdivisions;
! 54: with Bracket_Expansions; use Bracket_Expansions;
! 55: with SAGBI_Homotopies; use SAGBI_Homotopies;
! 56: with Matrix_Homotopies;
! 57: with Matrix_Homotopies_io;
! 58: with Osculating_Planes; use Osculating_Planes;
! 59:
! 60: procedure Driver_for_SAGBI_Homotopies
! 61: ( file : in file_type; n,d : in natural ) is
! 62:
! 63: dim : constant natural := (n-d)*d;
! 64:
! 65: function Convert ( mat : Standard_Floating_Matrices.Matrix )
! 66: return Standard_Complex_Matrices.Matrix is
! 67:
! 68: res : Standard_Complex_Matrices.Matrix(mat'range(1),mat'range(2));
! 69:
! 70: begin
! 71: for i in mat'range(1) loop
! 72: for j in mat'range(2) loop
! 73: res(i,j) := Create(mat(i,j));
! 74: end loop;
! 75: end loop;
! 76: return res;
! 77: end Convert;
! 78:
! 79: function Random_Osculating_SAGBI_Homotopy
! 80: ( file : file_type; locmap : Standard_Natural_Matrices.Matrix )
! 81: return Poly_Sys is
! 82:
! 83: -- DESCRIPTION :
! 84: -- Generates planes that are osculating a rational normal curve.
! 85: -- The system on return is the SAGBI homotopy.
! 86:
! 87: p : Poly;
! 88: l : Poly_Sys(1..dim);
! 89: s : double_float;
! 90: inc : double_float := 2.0/(double_float(dim));
! 91: mat : Standard_Floating_Matrices.Matrix(1..n,1..n-d);
! 92:
! 93: begin
! 94: -- p := Lifted_Localized_Laplace_Expansion(n,d);
! 95: p := Lifted_Localized_Laplace_Expansion(locmap);
! 96: new_line(file);
! 97: put_line(file,"The selected s-values : ");
! 98: s := Random;
! 99: for i in l'range loop
! 100: s := s + inc;
! 101: if s >= 1.0
! 102: then s := s - 2.0;
! 103: end if; -- s lies in [-1.0,+1.0]
! 104: put(file,s); new_line(file);
! 105: mat := Orthogonal_Basis(n,n-d,s);
! 106: Matrix_Homotopies.Add_Target(i,Convert(mat));
! 107: l(i) := Intersection_Condition(mat,p);
! 108: end loop;
! 109: return l;
! 110: end Random_Osculating_SAGBI_Homotopy;
! 111:
! 112: function Given_Osculating_SAGBI_Homotopy
! 113: ( file : file_type; locmap : Standard_Natural_Matrices.Matrix )
! 114: return Poly_Sys is
! 115:
! 116: -- DESCRIPTION :
! 117: -- Reads s-values and generates planes that are osculating a
! 118: -- rational normal curve.
! 119: -- The system on return is the SAGBI homotopy.
! 120:
! 121: p : Poly;
! 122: l : Poly_Sys(1..dim);
! 123: s : Standard_Floating_Vectors.Vector(1..dim);
! 124: mat : Standard_Floating_Matrices.Matrix(1..n,1..n-d);
! 125:
! 126: begin
! 127: -- p := Lifted_Localized_Laplace_Expansion(n,d);
! 128: p := Lifted_Localized_Laplace_Expansion(locmap);
! 129: new_line;
! 130: put("Give "); put(dim,1); put_line(" distinct real values for s : ");
! 131: for i in s'range loop
! 132: Read_Double_Float(s(i));
! 133: end loop;
! 134: new_line(file);
! 135: put_line(file,"The selected s-values : ");
! 136: put_line(file,s);
! 137: for i in l'range loop
! 138: mat := Orthogonal_Basis(n,n-d,s(i));
! 139: Matrix_Homotopies.Add_Target(i,Convert(mat));
! 140: l(i) := Intersection_Condition(mat,p);
! 141: end loop;
! 142: return l;
! 143: end Given_Osculating_SAGBI_Homotopy;
! 144:
! 145: function Read_Input_SAGBI_Homotopy
! 146: ( file : file_type; locmap : Standard_Natural_Matrices.Matrix )
! 147: return Poly_Sys is
! 148:
! 149: -- DESCRIPTION :
! 150: -- Reads the input planes from file.
! 151: -- The system on return is the SAGBI homotopy.
! 152:
! 153: planesfile : file_type;
! 154: p : Poly;
! 155: l : Poly_Sys(1..dim);
! 156: mat : Standard_Floating_Matrices.Matrix(1..n,1..n-d);
! 157:
! 158: begin
! 159: -- p := Lifted_Localized_Laplace_Expansion(n,d);
! 160: p := Lifted_Localized_Laplace_Expansion(locmap);
! 161: new_line;
! 162: put_line("Reading the name of the file with the input planes.");
! 163: Read_Name_and_Open_File(planesfile);
! 164: new_line(file);
! 165: put_line(file,"The input planes : ");
! 166: for i in l'range loop
! 167: get(planesfile,mat);
! 168: new_line(file); put(file,mat);
! 169: mat := Orthogonalize(mat);
! 170: Matrix_Homotopies.Add_Target(i,Convert(mat));
! 171: l(i) := Intersection_Condition(mat,p);
! 172: end loop;
! 173: Close(planesfile);
! 174: return l;
! 175: end Read_Input_SAGBI_Homotopy;
! 176:
! 177: function Complex_Random_SAGBI_Homotopy
! 178: ( locmap : Standard_Natural_Matrices.Matrix ) return Poly_Sys is
! 179:
! 180: -- DESCRIPTION :
! 181: -- Generates a SAGBI homotopy by taking random complex (n-d)-planes.
! 182: -- This system will be the start system in the Cheater's homotopy.
! 183:
! 184: p : Poly;
! 185: l : Poly_Sys(1..dim);
! 186: mat : Standard_Complex_Matrices.Matrix(1..n,1..n-d);
! 187:
! 188: begin
! 189: -- p := Lifted_Localized_Laplace_Expansion(n,d);
! 190: p := Lifted_Localized_Laplace_Expansion(locmap);
! 191: for i in l'range loop
! 192: mat := Random_Orthogonal_Matrix(n,n-d);
! 193: Matrix_Homotopies.Add_Start(i,mat);
! 194: l(i) := Intersection_Condition(mat,p);
! 195: end loop;
! 196: return l;
! 197: end Complex_Random_SAGBI_Homotopy;
! 198:
! 199: function Start_System ( l : Poly_Sys ) return Poly_Sys is
! 200:
! 201: -- DESCRIPTION :
! 202: -- Returns the system l after evaluation for t = 0.
! 203: -- This will be the start system in the SAGBI homotopy, flat deformation.
! 204:
! 205: r : Poly_Sys(l'range);
! 206: m : constant natural := Number_of_Unknowns(l(l'first));
! 207:
! 208: begin
! 209: for i in l'range loop
! 210: r(i) := Eval(l(i),Create(0.0),m);
! 211: end loop;
! 212: return r;
! 213: end Start_System;
! 214:
! 215: function Target_System ( l : Poly_Sys ) return Poly_Sys is
! 216:
! 217: -- DESCRIPTION :
! 218: -- Returns the system l after evaluation for t = 1.
! 219: -- This will be the target system in the SAGBI homotopy, flat deformation.
! 220:
! 221: r : Poly_Sys(l'range);
! 222: m : constant natural := Number_of_Unknowns(l(l'first));
! 223:
! 224: begin
! 225: for i in l'range loop
! 226: r(i) := Eval(l(i),Create(1.0),m);
! 227: end loop;
! 228: return r;
! 229: end Target_System;
! 230:
! 231: procedure Polyhedral_Homotopy_Continuation
! 232: ( file : in file_type;
! 233: p : in Poly_Sys; sols : in out Solution_List;
! 234: t : in Triangulation; lifted : in List; rep : in boolean ) is
! 235:
! 236: -- DESCRIPTION :
! 237: -- This is the first continuation stage, the resolution of the start system
! 238: -- in the SAGBI homotopy by means of polyhedral homotopy continuation.
! 239:
! 240: use Standard_Complex_Laur_JacoMats;
! 241:
! 242: timer : Timing_Widget;
! 243: lifted_lq,lq : Laur_Sys(p'range);
! 244: mix : Standard_Integer_Vectors.Vector(1..1) := (1..1 => p'last);
! 245: lif : Array_of_Lists(1..1) := (1..1 => lifted);
! 246: h : Eval_Coeff_Laur_Sys(p'range);
! 247: c : Standard_Complex_VecVecs.VecVec(h'range);
! 248: e : Exponent_Vectors_Array(h'range);
! 249: j : Eval_Coeff_Jaco_Mat(h'range,h'first..h'last+1);
! 250: m : Mult_Factors(j'range(1),j'range(2));
! 251: mixsub : Mixed_Subdivision;
! 252:
! 253: begin
! 254: tstart(timer);
! 255: mixsub := Shallow_Create(p'last,t);
! 256: lq := Polynomial_to_Laurent_System(p);
! 257: lifted_lq := Perform_Lifting(p'last,mix,lif,lq);
! 258: h := Create(lq);
! 259: for i in c'range loop
! 260: declare
! 261: coeff_lq : constant Standard_Complex_Vectors.Vector := Coeff(lq(i));
! 262: begin
! 263: c(i) := new Standard_Complex_Vectors.Vector(coeff_lq'range);
! 264: for k in coeff_lq'range loop
! 265: c(i)(k) := coeff_lq(k);
! 266: end loop;
! 267: end;
! 268: end loop;
! 269: e := Create(lq);
! 270: Create(lq,j,m);
! 271: if rep
! 272: then Mixed_Solve(file,lifted_lq,lif,h,c,e,j,m,mix,mixsub,sols);
! 273: else Mixed_Solve(lifted_lq,lif,h,c,e,j,m,mix,mixsub,sols);
! 274: end if;
! 275: tstop(timer);
! 276: new_line(file);
! 277: print_times(file,timer,"Polyhedral Homotopy Continuation");
! 278: end Polyhedral_Homotopy_Continuation;
! 279:
! 280: procedure Solve_Start_System
! 281: ( file : in file_type;
! 282: p : in Poly_Sys; sols : in out Solution_List;
! 283: report : in boolean ) is
! 284:
! 285: -- DESCRIPTION :
! 286: -- Computes the volume of the support of p. This is the
! 287: -- combinatorial-geometric set up for the first continuation stage.
! 288:
! 289: timer : Timing_Widget;
! 290: support : List := Create(p(p'first));
! 291: lifted,lifted_last : List;
! 292: t : Triangulation;
! 293: vol : natural;
! 294:
! 295: begin
! 296: new_line(file);
! 297: put_line(file,"The support of the start system : ");
! 298: new_line(file);
! 299: put(file,support);
! 300: tstart(timer);
! 301: Dynamic_Lifting(support,false,false,0,lifted,lifted_last,t);
! 302: new_line(file);
! 303: put_line(file,"The lifted support : ");
! 304: new_line(file);
! 305: put(file,lifted);
! 306: vol := Volume(t);
! 307: tstop(timer);
! 308: new_line(file);
! 309: put(file,"The volume : "); put(file,vol,1); new_line(file);
! 310: new_line(file);
! 311: print_times(file,timer,"Triangulation and Volume Computation");
! 312: Polyhedral_Homotopy_Continuation(file,p,sols,t,lifted,report);
! 313: Clear(t);
! 314: end Solve_Start_System;
! 315:
! 316: procedure Flat_Deformation
! 317: ( file : in file_type; sh : in Poly_Sys;
! 318: sols : in out Solution_List; report : in boolean ) is
! 319:
! 320: -- DESCRIPTION :
! 321: -- Performs flat deformation as defined by the SAGBI homotopy.
! 322: -- This is the second stage in the continuation.
! 323:
! 324: use Standard_Complex_Jaco_Matrices;
! 325:
! 326: timer : Timing_Widget;
! 327: sh_eval : Eval_Poly_Sys(sh'range);
! 328: jac_mat : Jaco_Mat(sh'range,sh'first..sh'last+1);
! 329: eva_jac : Eval_Jaco_Mat(jac_mat'range(1),jac_mat'range(2));
! 330:
! 331: function Eval ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
! 332: return Standard_Complex_Vectors.Vector is
! 333:
! 334: xt : Standard_Complex_Vectors.Vector(x'first..x'last+1);
! 335:
! 336: begin
! 337: xt(x'range) := x;
! 338: xt(xt'last) := t;
! 339: return Eval(sh_eval,xt);
! 340: end Eval;
! 341:
! 342: function Diff ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
! 343: return Standard_Complex_Matrices.Matrix is
! 344:
! 345: res : Standard_Complex_Matrices.Matrix(x'range,x'range);
! 346: xt : Standard_Complex_Vectors.Vector(x'first..x'last+1);
! 347:
! 348: begin
! 349: xt(x'range) := x;
! 350: xt(xt'last) := t;
! 351: for i in res'range(1) loop
! 352: for j in res'range(2) loop
! 353: res(i,j) := Eval(eva_jac(i,j),xt);
! 354: end loop;
! 355: end loop;
! 356: return res;
! 357: end Diff;
! 358:
! 359: function Diff ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
! 360: return Standard_Complex_Vectors.Vector is
! 361:
! 362: res : Standard_Complex_Vectors.Vector(x'range);
! 363: xt : Standard_Complex_Vectors.Vector(x'first..x'last+1);
! 364:
! 365: begin
! 366: xt(x'range) := x;
! 367: xt(xt'last) := t;
! 368: for i in res'range loop
! 369: res(i) := Eval(eva_jac(i,eva_jac'last(2)),xt);
! 370: end loop;
! 371: return res;
! 372: end Diff;
! 373:
! 374: procedure Sil_Cont is new Silent_Continue(Max_Norm,Eval,Diff,Diff);
! 375: procedure Rep_Cont is new Reporting_Continue(Max_Norm,Eval,Diff,Diff);
! 376:
! 377: begin
! 378: tstart(timer);
! 379: sh_eval := Create(sh);
! 380: jac_mat := Create(sh);
! 381: eva_jac := Create(jac_mat);
! 382: Set_Continuation_Parameter(sols,Create(0.0));
! 383: if report
! 384: then Rep_Cont(file,sols,false,Create(1.0));
! 385: else Sil_Cont(sols,false,Create(1.0));
! 386: end if;
! 387: tstop(timer);
! 388: new_line(file); print_times(file,timer,"Flat Deformation");
! 389: Clear(sh_eval);
! 390: Clear(jac_mat);
! 391: Clear(eva_jac);
! 392: end Flat_Deformation;
! 393:
! 394: procedure Solve_Target_System
! 395: ( file : in file_type; start,target : in Poly_Sys;
! 396: sols : in out Solution_List; report : in boolean ) is
! 397:
! 398: -- DESCRIPTION :
! 399: -- This is the third and last continuation stage: Cheater's homotopy.
! 400: -- It is implemented in the space of polynomials.
! 401:
! 402: timer : Timing_Widget;
! 403: n : constant natural := target'last;
! 404: a : Standard_Complex_Vectors.Vector(1..n) := Random_Vector(1,n);
! 405: b : Standard_Complex_Vectors.Vector(1..n) := Random_Vector(1,n);
! 406:
! 407: begin
! 408: tstart(timer);
! 409: Homotopy.Create(target,start,2,a,b,true); -- linear cheater
! 410: Set_Continuation_Parameter(sols,Create(0.0));
! 411: declare
! 412: procedure Sil_Cont is
! 413: new Silent_Continue(Max_Norm,
! 414: Homotopy.Eval,Homotopy.Diff,Homotopy.Diff);
! 415: procedure Rep_Cont is
! 416: new Reporting_Continue(Max_Norm,
! 417: Homotopy.Eval,Homotopy.Diff,Homotopy.Diff);
! 418: begin
! 419: if report
! 420: then Rep_Cont(file,sols,false,Create(1.0));
! 421: else Sil_Cont(sols,false,Create(1.0));
! 422: end if;
! 423: end;
! 424: tstop(timer);
! 425: new_line(file);
! 426: print_times(file,timer,"Cheater's homotopy to target system");
! 427: end Solve_Target_System;
! 428:
! 429: procedure Solve_Target_System
! 430: ( file : in file_type; symtarget : in Poly;
! 431: sols : in out Solution_List; report : in boolean ) is
! 432:
! 433: -- DESCRIPTION :
! 434: -- This is the third and last continuation stage: Cheater's homotopy.
! 435: -- It uses a coefficient-matrix homotopy and takes as input the target
! 436: -- system with coefficients as brackets.
! 437: -- This is far more expensive than the linear Cheater's homotopy.
! 438:
! 439: use Standard_Complex_Jaco_Matrices;
! 440:
! 441: timer : Timing_Widget;
! 442: h : Standard_Complex_Poly_Functions.Eval_Coeff_Poly := Create(symtarget);
! 443: homsys : Poly_Sys(1..dim);
! 444: jacmat : Eval_Coeff_Jaco_Mat(1..dim,1..dim);
! 445: mulfac : Mult_Factors(1..dim,1..dim);
! 446: brkcff : Standard_Complex_Vectors.Vector := Coeff(symtarget);
! 447: syscff : Standard_Complex_VecVecs.VecVec(1..dim);
! 448:
! 449: begin
! 450: tstart(timer);
! 451: Set_Continuation_Parameter(sols,Create(0.0));
! 452: for i in homsys'range loop
! 453: Copy(symtarget,homsys(i));
! 454: end loop;
! 455: Create(homsys,jacmat,mulfac);
! 456: for i in syscff'range loop
! 457: syscff(i) := new Standard_Complex_Vectors.Vector(brkcff'range);
! 458: end loop;
! 459: declare
! 460:
! 461: function Eval ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
! 462: return Standard_Complex_Vectors.Vector is
! 463:
! 464: y : Standard_Complex_Vectors.Vector(x'range);
! 465: c : Standard_Complex_Vectors.Vector(brkcff'range);
! 466:
! 467: begin
! 468: for i in y'range loop
! 469: c := Intersection_Coefficients(Matrix_Homotopies.Eval(i,t),brkcff);
! 470: y(i) := Eval(h,c,x);
! 471: end loop;
! 472: return y;
! 473: end Eval;
! 474:
! 475: function dHt ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
! 476: return Standard_Complex_Vectors.Vector is
! 477:
! 478: y : Standard_Complex_Vectors.Vector(x'range) := x;
! 479:
! 480: begin
! 481: return y;
! 482: end dHt;
! 483:
! 484: function dHx ( x : Standard_Complex_Vectors.Vector; t : Complex_Number )
! 485: return Standard_Complex_Matrices.Matrix is
! 486: begin
! 487: for i in syscff'range loop
! 488: syscff(i).all
! 489: := Intersection_Coefficients(Matrix_Homotopies.Eval(i,t),brkcff);
! 490: end loop;
! 491: return Eval(jacmat,mulfac,syscff,x);
! 492: end dHx;
! 493:
! 494: procedure Sil_Cont is new Silent_Continue(Max_Norm,Eval,dHt,dHx);
! 495: procedure Rep_Cont is new Reporting_Continue(Max_Norm,Eval,dHt,dHx);
! 496:
! 497: begin
! 498: if report
! 499: then Rep_Cont(file,sols,false,Create(1.0));
! 500: else Sil_Cont(sols,false,Create(1.0));
! 501: end if;
! 502: end;
! 503: Standard_Complex_VecVecs.Clear(syscff);
! 504: tstop(timer);
! 505: new_line(file);
! 506: print_times(file,timer,"Cheater's homotopy to target system");
! 507: end Solve_Target_System;
! 508:
! 509: procedure Refine_Roots ( file : in file_type;
! 510: p : in Poly_Sys; sols : in out Solution_List ) is
! 511:
! 512: epsxa : constant double_float := 10.0**(-8);
! 513: epsfa : constant double_float := 10.0**(-8);
! 514: tolsing : constant double_float := 10.0**(-8);
! 515: max : constant natural := 3;
! 516: numit : natural := 0;
! 517:
! 518: begin
! 519: Reporting_Root_Refiner(file,p,sols,epsxa,epsfa,tolsing,numit,max,false);
! 520: end Refine_Roots;
! 521:
! 522: procedure Set_Parameters ( file : in file_type; report : out boolean ) is
! 523:
! 524: -- DESCRIPTION :
! 525: -- Interactive determination of the continuation and output parameters.
! 526:
! 527: oc : natural;
! 528:
! 529: begin
! 530: new_line;
! 531: Driver_for_Continuation_Parameters(file);
! 532: new_line;
! 533: Driver_for_Process_io(file,oc);
! 534: report := not (oc = 0);
! 535: new_line;
! 536: put_line("See the output file for results.");
! 537: new_line;
! 538: end Set_Parameters;
! 539:
! 540: function t_Symbol return Symbol is
! 541:
! 542: -- DESCRIPTION :
! 543: -- Returns the symbol to represent the continuation parameter t.
! 544:
! 545: res : Symbol;
! 546:
! 547: begin
! 548: res(1) := 't';
! 549: for i in 2..res'last loop
! 550: res(i) := ' ';
! 551: end loop;
! 552: return res;
! 553: end t_Symbol;
! 554:
! 555: function Lowest_Degree_Localization_Map
! 556: ( n,d : natural ) return Standard_Natural_Matrices.Matrix is
! 557:
! 558: -- DESCRIPTION :
! 559: -- Returns the lowest-degree localization map.
! 560:
! 561: res : Standard_Natural_Matrices.Matrix(1..n,1..d);
! 562:
! 563: begin
! 564: for i in 1..n-d loop
! 565: for j in 1..d loop
! 566: res(i,j) := 2;
! 567: end loop;
! 568: end loop;
! 569: for i in n-d+1..n loop
! 570: for j in 1..d loop
! 571: if i = j+n-d
! 572: then res(i,j) := 1;
! 573: else res(i,j) := 0;
! 574: end if;
! 575: end loop;
! 576: end loop;
! 577: return res;
! 578: end Lowest_Degree_Localization_Map;
! 579:
! 580: function Read_Localization_Map
! 581: ( n,d : natural ) return Standard_Natural_Matrices.Matrix is
! 582:
! 583: -- DESCRIPTION :
! 584: -- Reads a localization map from standard input.
! 585:
! 586: res : Standard_Natural_Matrices.Matrix(1..n,1..d);
! 587:
! 588: begin
! 589: put_line("Reading localization map: 0,1 for I and 2 = free element.");
! 590: put_line("The lowest-degree localization map looks like :");
! 591: put(Lowest_Degree_Localization_Map(n,d));
! 592: put_line("and the general localization map is :");
! 593: put(Localization_Map(n,d));
! 594: put("Give a "); put(n,1); put("-by-"); put(d,1);
! 595: put_line(" matrix to represent your own map :");
! 596: get(res); skip_line;
! 597: return res;
! 598: end Read_Localization_Map;
! 599:
! 600: procedure Main is
! 601:
! 602: -- DESCRIPTION :
! 603: -- There are four parts in the elaboration of the SAGBI homotopy :
! 604: -- 0. Set up of the polynomials in the SAGBI homotopy;
! 605: -- 1. Solve the start system by polyhedral continuation;
! 606: -- 2. Apply the flat deformations to solve complex instance;
! 607: -- 3. Cheater's homotopy from complex to real instance.
! 608:
! 609: timer,totaltimer : Timing_Widget;
! 610: realsagbih,compsagbih,realtarget,comptarget,starts : Poly_Sys(1..dim);
! 611: sols : Solution_List;
! 612: report : boolean;
! 613: inputchoice,localchoice : character;
! 614: locmap : Standard_Natural_Matrices.Matrix(1..n,1..d);
! 615:
! 616: begin
! 617: new_line;
! 618: put_line("MENU for a localization for the output planes : ");
! 619: put_line(" 1. Lowest-degree map, with I in lower-right corner.");
! 620: put_line(" 2. General map, able to represent any intersection.");
! 621: put_line(" 3. Give in your own localization map.");
! 622: put("Type 1, 2, or 3 to select : "); Ask_Alternative(localchoice,"123");
! 623: case localchoice is
! 624: when '1' => locmap := Lowest_Degree_Localization_Map(n,d);
! 625: when '2' => locmap := Localization_Map(n,d);
! 626: when '3' => locmap := Read_Localization_Map(n,d);
! 627: when others => null;
! 628: end case;
! 629: new_line;
! 630: put_line("MENU for constructing target system based on input planes : ");
! 631: put_line(" 1. Generate input planes osculating a rational normal curve.");
! 632: put_line(" 2. Interactively give s-values for the "
! 633: & "osculating input planes.");
! 634: put_line(" 3. Give the name of the file with input planes.");
! 635: put("Type 1, 2, or 3 to select : "); Ask_Alternative(inputchoice,"123");
! 636: tstart(totaltimer);
! 637: tstart(timer);
! 638: Matrix_Indeterminates.Initialize_Symbols(n,d);
! 639: Matrix_Indeterminates.Reduce_Symbols(locmap);
! 640: Matrix_Homotopies.Init(dim);
! 641: case inputchoice is
! 642: when '1' => realsagbih := Random_Osculating_SAGBI_Homotopy(file,locmap);
! 643: when '2' => realsagbih := Given_Osculating_SAGBI_Homotopy(file,locmap);
! 644: when '3' => realsagbih := Read_Input_SAGBI_Homotopy(file,locmap);
! 645: when others => null;
! 646: end case;
! 647: realtarget := Target_System(realsagbih);
! 648: compsagbih := Complex_Random_SAGBI_Homotopy(locmap);
! 649: comptarget := Target_System(compsagbih);
! 650: starts := Start_System(compsagbih);
! 651: Symbol_Table.Add(t_Symbol);
! 652: new_line(file);
! 653: put_line(file,"The target polynomial system with real coefficients : ");
! 654: new_line(file);
! 655: put(file,realtarget'last,realtarget);
! 656: new_line(file);
! 657: put_line(file,"The localization map : "); put(file,locmap);
! 658: new_line(file);
! 659: Matrix_Homotopies_io.Write(file);
! 660: put_line(file,"The target polynomial system with complex coefficients : ");
! 661: put_line(file,comptarget);
! 662: new_line(file);
! 663: put_line(file,"The SAGBI Homotopy as complex polynomial system : ");
! 664: new_line(file);
! 665: put_line(file,compsagbih);
! 666: new_line(file);
! 667: put_line(file,"The SAGBI Homotopy at t=0, the start system : ");
! 668: new_line(file);
! 669: put_line(file,starts);
! 670: tstop(timer);
! 671: new_line(file);
! 672: print_times(file,timer,"Setting up the polynomials in the SAGBI homotopy");
! 673: Set_Parameters(file,report);
! 674: Solve_Start_System(file,starts,sols,report);
! 675: Flat_Deformation(file,compsagbih,sols,report);
! 676: Solve_Target_System(file,comptarget,realtarget,sols,report);
! 677: -- Solve_Target_System(file,symtarget,sols,report);
! 678: -- this is the determinantal cheater's homotopy, very expensive
! 679: Refine_Roots(file,realtarget,sols);
! 680: Matrix_Indeterminates.Clear_Symbols;
! 681: tstop(totaltimer);
! 682: new_line(file);
! 683: print_times(file,totaltimer,"Total time for Solving the System");
! 684: end Main;
! 685:
! 686: begin
! 687: Main;
! 688: end Driver_for_SAGBI_Homotopies;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to driver_for_sagbi_homotopies.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Schubert