Annotation of OpenXM_contrib/PHC/Ada/Homotopy/drivers_for_homotopy_creation.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 Standard_Floating_Numbers; use Standard_Floating_Numbers;
4: with Standard_Complex_Numbers_io; use Standard_Complex_Numbers_io;
5: with Standard_Random_Numbers; use Standard_Random_Numbers;
6: with Numbers_io; use Numbers_io;
7: with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
8: with Homotopy;
9: with Projective_Transformations; use Projective_Transformations;
10: with Homogenization; use Homogenization;
11:
12: package body Drivers_for_Homotopy_Creation is
13:
14: procedure Info_on_Power is
15:
16: i : array(1..6) of string(1..65);
17:
18: begin
19: i(1):=" The homotopy parameter k determines the power of the";
20: i(2):="continuation parameter t. Taking k>1 has as effect that at the";
21: i(3):="beginning and at the end of the continuation, changes in t do not";
22: i(4):="change the homotopy as much as with a homotopy that is linear in";
23: i(5):="t so that paths are better to follow. The default value k=2 is";
24: i(6):="recommended. ";
25: for j in i'range loop
26: put_line(i(j));
27: end loop;
28: end Info_on_Power;
29:
30: procedure Info_on_Constant is
31:
32: i : array(1..3) of string(1..65);
33:
34: begin
35: i(1):=" The homotopy parameter a ensures the regularity of the solution";
36: i(2):="paths, i.e.: by choosing a random complex number for a, all paths";
37: i(3):="are regular and do not diverge for t<1. ";
38: for j in i'range loop
39: put_line(i(j));
40: end loop;
41: end Info_on_Constant;
42:
43: procedure Info_on_Target is
44:
45: i : array(1..7) of string(1..65);
46:
47: begin
48: i(1):=" The target value for the continuation parameter t is by default";
49: i(2):="1. To create stepping stones in the continuation stage, it is";
50: i(3):="possible to let the continuation stop at t<1, for instance at t =";
51: i(4):="0.9 or even at a complex value for t. The solutions at t<1 will";
52: i(5):="serve at start solutions to take up the continuation in a later";
53: i(6):="stage. In this stage, the same homotopy parameters k and a must";
54: i(7):="be used. ";
55: for j in i'range loop
56: put_line(i(j));
57: end loop;
58: end Info_on_Target;
59:
60: procedure Info_on_Exit is
61:
62: i : constant string
63: :="By typing 0, the current settings are used in the homotopy.";
64:
65: begin
66: put_line(i);
67: end Info_on_Exit;
68:
69: procedure Info_on_Projective_Transformation is
70:
71: i : array(1..5) of string(1..65);
72:
73: begin
74: i(1):=" A projective transformation of the homotopy and start solutions";
75: i(2):="makes the equations in the polynomials homogeneous and adds an a";
76: i(3):="random hyperplane. The vectors of the start solutions are";
77: i(4):="extended with an additional unknown. For solutions at infinity,";
78: i(5):="this additional unknown equals zero. ";
79: for j in i'range loop
80: put_line(i(j));
81: end loop;
82: end Info_on_Projective_Transformation;
83:
84: procedure Read_Complex ( x : in out Complex_Number ) is
85:
86: -- DESCRIPTION :
87: -- User friendly reading of a complex number.
88:
89: re,im : double_float;
90:
91: begin
92: put(" Give a real number for real part : "); Read_Double_Float(re);
93: put(" Give a real number for imaginary part : "); Read_Double_Float(im);
94: x := Create(re,im);
95: end Read_Complex;
96:
97: procedure Default_Homotopy_Settings
98: ( k : out positive; a,t : out Complex_Number ) is
99: begin
100: k := 2;
101: a := Random1;
102: t := Create(1.0);
103: end Default_Homotopy_Settings;
104:
105: procedure Default_Homotopy_Settings
106: ( k : out positive; a,t : out Complex_Number;
107: proj : out boolean ) is
108: begin
109: Default_Homotopy_Settings(k,a,t);
110: proj := false;
111: end Default_Homotopy_Settings;
112:
113: procedure Menu_for_Homotopy_Settings
114: ( file : in file_type; k : in out positive;
115: a,t : in out Complex_Number; proj : in out boolean ) is
116:
117: ans : character;
118: choice : string(1..2) := " ";
119:
120: begin
121: new_line;
122: put_line
123: ("Homotopy is H(x,t) = a*(1-t)^k * Q(x) + t^k * P(x) = 0, t in [0,1],");
124: put_line
125: (" with Q(x) = 0 a start system, and P(x) = 0 the target system.");
126: loop
127: new_line;
128: put_line("MENU with settings for homotopy :");
129: put(" relaxation constant k : "); put(k,2); new_line;
130: put(" regularity constant a : "); put(a); new_line;
131: put(" target value for t : "); put(t); new_line;
132: put(" projective transformation : ");
133: if proj then put("yes"); else put("no"); end if; new_line;
134: put("Type k, a, t, or p to change, preceded by i for info."
135: & " Type 0 to exit : ");
136: Ask_Alternative(choice,"katp0",'i');
137: exit when choice(1) = '0';
138: case choice(1) is
139: when 'k' =>
140: put("Give a value (1,2, or 3) for the relaxation constant k : ");
141: Read_Positive(k);
142: when 'a' =>
143: put_line("Reading a complex value for the regularity constant a.");
144: loop
145: Read_Complex(a);
146: exit when AbsVal(a) /= 0.0;
147: put_line("The value 0 for a is not allowed! Try again.");
148: end loop;
149: when 't' =>
150: put_line("Reading the (complex) target value for t.");
151: Read_Complex(t);
152: when 'p' =>
153: put("Do you want a projective transformation? ");
154: Ask_Yes_or_No(ans); proj := (ans ='y');
155: when 'i' =>
156: new_line;
157: case choice(2) is
158: when 'k' => Info_on_Power;
159: when 'a' => Info_on_Constant;
160: when 't' => Info_on_Target;
161: when 'p' => Info_on_Projective_Transformation;
162: when '0' => Info_on_Exit;
163: when others => null;
164: end case;
165: when others => null;
166: end case;
167: end loop;
168: new_line(file);
169: put_line(file,"HOMOTOPY PARAMETERS :");
170: put(file," k : "); put(file,k,2); new_line(file);
171: put(file," a : "); put(file,a); new_line(file);
172: put(file," t : "); put(file,t); new_line(file);
173: if proj
174: then put_line(file," projective transformation");
175: else put_line(file," no projective transformation");
176: end if;
177: end Menu_for_Homotopy_Settings;
178:
179: procedure Menu_for_Homotopy_Settings
180: ( file : in file_type; k : in out positive;
181: a,t : in out Complex_Number ) is
182:
183: choice : string(1..2) := " ";
184:
185: begin
186: new_line;
187: put_line
188: ("Homotopy is H(x,t) = a*(1-t)^k * Q(x) + t^k * P(x) = 0, t in [0,1],");
189: put_line
190: (" with Q(x) = 0 a start system, and P(x) = 0 the target system.");
191: loop
192: new_line;
193: put_line("MENU with settings for homotopy :");
194: put(" relaxation constant k : "); put(k,2); new_line;
195: put(" regularity constant a : "); put(a); new_line;
196: put(" target value for t : "); put(t); new_line;
197: put("Type k, a, or t to change, preceded by i for info."
198: & " Type 0 to exit : ");
199: Ask_Alternative(choice,"kat0",'i');
200: exit when choice(1) = '0';
201: case choice(1) is
202: when 'k' =>
203: put("Give a value (1,2, or 3) for the relaxation constant k : ");
204: Read_Positive(k);
205: when 'a' =>
206: put_line("Reading a complex value for the regularity constant a.");
207: loop
208: Read_Complex(a);
209: exit when AbsVal(a) /= 0.0;
210: put_line("The value 0 for a is not allowed! Try again.");
211: end loop;
212: when 't' =>
213: put_line("Reading the (complex) target value for t.");
214: Read_Complex(t);
215: when 'i' =>
216: new_line;
217: case choice(2) is
218: when 'k' => Info_on_Power;
219: when 'a' => Info_on_Constant;
220: when 't' => Info_on_Target;
221: when '0' => Info_on_Exit;
222: when others => null;
223: end case;
224: when others => null;
225: end case;
226: end loop;
227: new_line(file);
228: put_line(file,"HOMOTOPY PARAMETERS :");
229: put(file," k : "); put(file,k,2); new_line(file);
230: put(file," a : "); put(file,a); new_line(file);
231: put(file," t : "); put(file,t); new_line(file);
232: end Menu_for_Homotopy_Settings;
233:
234: procedure Driver_for_Homotopy_Construction
235: ( file : in file_type; p,q : in out Poly_Sys;
236: qsols : in out Solution_List; target : out Complex_Number ) is
237:
238: k : positive;
239: a,t : Complex_Number;
240: prt : boolean;
241:
242: begin
243: Default_Homotopy_Settings(k,a,t,prt);
244: Menu_for_Homotopy_Settings(file,k,a,t,prt);
245: target := t;
246: if prt
247: then Projective_Transformation(p);
248: Projective_Transformation(q);
249: Projective_Transformation(qsols);
250: declare
251: pp,sysp : Poly_Sys(p'first..p'last+1);
252: begin
253: pp := Add_Random_Hyperplanes(p,1,true);
254: sysp := Add_Standard_Hyperplanes(q,1);
255: Homotopy.Create(pp,sysp,k,a);
256: Clear(pp); Clear(sysp);
257: end;
258: else Homotopy.Create(p,q,k,a);
259: end if;
260: end Driver_for_Homotopy_Construction;
261:
262: end Drivers_for_Homotopy_Creation;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to drivers_for_homotopy_creation.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Homotopy