[BACK]Return to drivers_for_input_planes.adb CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / PHC / Ada / Schubert

Annotation of OpenXM_contrib/PHC/Ada/Schubert/drivers_for_input_planes.adb, Revision 1.1.1.1

1.1       maekawa     1: with integer_io,Numbers_io;              use integer_io,Numbers_io;
                      2: with Communications_with_User;           use Communications_with_User;
                      3: with Standard_Floating_Numbers;          use Standard_Floating_Numbers;
                      4: with Standard_Complex_Numbers;           use Standard_Complex_Numbers;
                      5: with Standard_Random_Numbers;            use Standard_Random_Numbers;
                      6: with Standard_Floating_Vectors_io;       use Standard_Floating_Vectors_io;
                      7: with Standard_Floating_Matrices;
                      8: with Standard_Floating_Matrices_io;      use Standard_Floating_Matrices_io;
                      9: with Standard_Complex_Matrices;
                     10: with Standard_Random_Matrices;           use Standard_Random_Matrices;
                     11: with Standard_Complex_VecMats_io;        use Standard_Complex_VecMats_io;
                     12: with Osculating_Planes;                  use Osculating_Planes;
                     13:
                     14: package body Drivers_for_Input_Planes is
                     15:
                     16:   function Finite ( dim : Bracket; fin_sum : natural ) return boolean is
                     17:
                     18:   -- DESCRIPTION :
                     19:   --   Returns true if the co-dimensions in dim will lead to a finite
                     20:   --   number of solutions.
                     21:
                     22:     sum : natural := 0;
                     23:
                     24:   begin
                     25:     for i in dim'range loop
                     26:       sum := sum + dim(i);
                     27:     end loop;
                     28:     if sum = fin_sum
                     29:      then return true;
                     30:      else return false;
                     31:     end if;
                     32:   end Finite;
                     33:
                     34:   function Read_Codimensions ( m,p,q : natural ) return Bracket is
                     35:
                     36:     mpq : constant natural := m*p + q*(m+p);
                     37:     codim : Bracket(1..mpq);
                     38:     n : natural;
                     39:
                     40:   begin
                     41:     loop
                     42:       put("Give number of intersection conditions : "); get(n);
                     43:       put("Give "); put(n,1); put(" co-dimensions (sum = ");
                     44:       put(mpq,1); put(") : ");
                     45:       for i in 1..n loop
                     46:         get(codim(i));
                     47:       end loop;
                     48:       for i in 1..n-1 loop
                     49:         put(codim(i),1); put(" + ");
                     50:       end loop;
                     51:       put(codim(n),1);
                     52:       if Finite(codim(1..n),mpq)
                     53:           then put(" = "); put(mpq,1); put_line("  Finite #sols.");
                     54:             exit;
                     55:           else put(" /= "); put(mpq,1);
                     56:             put_line("  Please try again.");
                     57:       end if;
                     58:     end loop;
                     59:     return codim(1..n);
                     60:   end Read_Codimensions;
                     61:
                     62:   function Convert ( mat : Standard_Floating_Matrices.Matrix )
                     63:                    return Standard_Complex_Matrices.Matrix is
                     64:
                     65:   -- DESCRIPTION :
                     66:   --   Converts a real matrix into a complex one.
                     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_Complex_Planes ( m,p : natural ) return VecMat is
                     80:
                     81:     res : VecMat(1..m*p);
                     82:     n : constant natural := m+p;
                     83:     use Standard_Complex_Matrices;
                     84:
                     85:   begin
                     86:     for i in res'range loop
                     87:       res(i) := new Matrix'(Random_Orthogonal_Matrix(n,m));
                     88:     end loop;
                     89:     return res;
                     90:   end Random_Complex_Planes;
                     91:
                     92:   function Random_Complex_Planes ( m,p : natural; k : Bracket ) return VecMat is
                     93:
                     94:     res : VecMat(k'range);
                     95:     n : constant natural := m+p;
                     96:     use Standard_Complex_Matrices;
                     97:
                     98:   begin
                     99:     for i in res'range loop
                    100:       res(i) := new Matrix'(Random_Orthogonal_Matrix(n,m+1-k(i)));
                    101:     end loop;
                    102:     return res;
                    103:   end Random_Complex_Planes;
                    104:
                    105:   function Random_Complex_Planes ( m,p,q : natural ) return VecMat is
                    106:
                    107:     res : VecMat(1..(m*p+q*(m+p)));
                    108:     n : constant natural := m+p;
                    109:     use Standard_Complex_Matrices;
                    110:
                    111:   begin
                    112:     for i in res'range loop
                    113:       res(i) := new Matrix'(Random_Orthogonal_Matrix(n,m));
                    114:     end loop;
                    115:     return res;
                    116:   end Random_Complex_Planes;
                    117:
                    118:   function Random_Real_Planes ( m,p : natural ) return VecMat is
                    119:
                    120:     res : VecMat(1..m*p);
                    121:     n : constant natural := m+p;
                    122:     fltmat : Standard_Floating_Matrices.Matrix(1..n,1..m);
                    123:     cmpmat : Standard_Complex_Matrices.Matrix(1..n,1..m);
                    124:
                    125:   begin
                    126:     for i in res'range loop
                    127:       fltmat := Random_Orthogonal_Matrix(n,m);
                    128:       cmpmat := Convert(fltmat);
                    129:       res(i) := new Standard_Complex_Matrices.Matrix'(cmpmat);
                    130:     end loop;
                    131:     return res;
                    132:   end Random_Real_Planes;
                    133:
                    134:   function Random_Real_Planes ( m,p : natural; k : Bracket ) return VecMat is
                    135:
                    136:     res : VecMat(k'range);
                    137:     n : constant natural := m+p;
                    138:
                    139:   begin
                    140:     for i in res'range loop
                    141:       declare
                    142:         fltmat : Standard_Floating_Matrices.Matrix(1..n,1..m+1-k(i));
                    143:         cmpmat : Standard_Complex_Matrices.Matrix(1..n,1..m+1-k(i));
                    144:       begin
                    145:         fltmat := Random_Orthogonal_Matrix(n,m+1-k(i));
                    146:         cmpmat := Convert(fltmat);
                    147:         res(i) := new Standard_Complex_Matrices.Matrix'(cmpmat);
                    148:       end;
                    149:     end loop;
                    150:     return res;
                    151:   end Random_Real_Planes;
                    152:
                    153:   function Random_Real_Planes ( m,p,q : natural ) return VecMat is
                    154:
                    155:     res : VecMat(1..(m*p+q*(m+p)));
                    156:     n : constant natural := m+p;
                    157:     fltmat : Standard_Floating_Matrices.Matrix(1..n,1..m);
                    158:     cmpmat : Standard_Complex_Matrices.Matrix(1..n,1..m);
                    159:
                    160:   begin
                    161:     for i in res'range loop
                    162:       fltmat := Random_Orthogonal_Matrix(n,m);
                    163:       cmpmat := Convert(fltmat);
                    164:       res(i) := new Standard_Complex_Matrices.Matrix'(cmpmat);
                    165:     end loop;
                    166:     return res;
                    167:   end Random_Real_Planes;
                    168:
                    169:   function Equidistant_Interpolation_Points ( n : natural ) return Vector is
                    170:
                    171:     res : Vector(1..n);
                    172:     s : double_float := Random;
                    173:     inc : constant double_float := 2.0/double_float(n);
                    174:
                    175:   begin
                    176:     for i in res'range loop
                    177:       res(i) := s;
                    178:       s := s+inc;
                    179:       if s >= 1.0
                    180:        then s := s - 2.0;
                    181:       end if;
                    182:     end loop;
                    183:     return res;
                    184:   end Equidistant_Interpolation_Points;
                    185:
                    186:   function Read_Interpolation_Points ( n : natural ) return Vector is
                    187:
                    188:     res : Vector(1..n);
                    189:
                    190:   begin
                    191:     new_line;
                    192:     put("Give "); put(n,1); put_line(" distinct real interpolation points : ");
                    193:     for i in res'range loop
                    194:       Read_Double_Float(res(i));
                    195:     end loop;
                    196:     return res;
                    197:   end Read_Interpolation_Points;
                    198:
                    199:   function Osculating_Input_Planes ( m,p : natural ) return VecMat is
                    200:
                    201:     s : constant Vector := Equidistant_Interpolation_Points(m*p);
                    202:
                    203:   begin
                    204:     return Osculating_Input_Planes(m,p,s);
                    205:   end Osculating_Input_Planes;
                    206:
                    207:   function Osculating_Input_Planes
                    208:              ( m,p : natural; s : Vector ) return VecMat is
                    209:
                    210:     res : VecMat(1..m*p);
                    211:     n : constant natural := m+p;
                    212:     realmat : Standard_Floating_Matrices.Matrix(1..n,1..m);
                    213:
                    214:   begin
                    215:     for i in res'range loop
                    216:       realmat := Orthogonal_Basis(n,m,s(i));
                    217:       res(i) := new Standard_Complex_Matrices.Matrix'(Convert(realmat));
                    218:     end loop;
                    219:     return res;
                    220:   end Osculating_Input_Planes;
                    221:
                    222:   function Osculating_Input_Planes
                    223:              ( m,p : natural; k : bracket ) return VecMat is
                    224:
                    225:     s : constant Vector := Equidistant_Interpolation_Points(k'length);
                    226:
                    227:   begin
                    228:     return Osculating_Input_Planes(m,p,k,s);
                    229:   end Osculating_Input_Planes;
                    230:
                    231:   function Osculating_Input_Planes
                    232:              ( m,p : natural; k : bracket; s : Vector ) return VecMat is
                    233:
                    234:     res : VecMat(k'range);
                    235:     n : constant natural := m+p;
                    236:
                    237:   begin
                    238:     for i in res'range loop
                    239:       declare
                    240:         realmat : Standard_Floating_Matrices.Matrix(1..n,1..m+1-k(i));
                    241:       begin
                    242:         realmat := Orthogonal_Basis(n,m+1-k(i),s(i));
                    243:         res(i) := new Standard_Complex_Matrices.Matrix'(Convert(realmat));
                    244:       end;
                    245:     end loop;
                    246:     return res;
                    247:   end Osculating_Input_Planes;
                    248:
                    249:   function Osculating_Input_Planes ( m,p,q : natural ) return VecMat is
                    250:
                    251:     dim : constant natural := m*p + q*(m+p);
                    252:     s : constant Vector := Equidistant_Interpolation_Points(dim);
                    253:
                    254:   begin
                    255:     return Osculating_Input_Planes(m,p,q,s);
                    256:   end Osculating_Input_Planes;
                    257:
                    258:   function Osculating_Input_Planes
                    259:              ( m,p,q : natural; s : Vector ) return VecMat is
                    260:
                    261:     res : VecMat(1..m*p+q*(m+p));
                    262:     n : constant natural := m+p;
                    263:     realmat : Standard_Floating_Matrices.Matrix(1..n,1..m);
                    264:
                    265:   begin
                    266:     for i in res'range loop
                    267:       realmat := Orthogonal_Basis(n,m,s(i));
                    268:       res(i) := new Standard_Complex_Matrices.Matrix'(Convert(realmat));
                    269:     end loop;
                    270:     return res;
                    271:   end Osculating_Input_Planes;
                    272:
                    273:   function Read_Input_Planes ( m,p : natural ) return VecMat is
                    274:
                    275:     res : VecMat(1..m*p);
                    276:     planesfile : file_type;
                    277:     n : constant natural := m+p;
                    278:     realmat : Standard_Floating_Matrices.Matrix(1..n,1..m);
                    279:
                    280:   begin
                    281:     new_line;
                    282:     put_line("Reading the name of the file with the input planes.");
                    283:     Read_Name_and_Open_File(planesfile);
                    284:     for i in res'range loop
                    285:       get(planesfile,realmat);
                    286:       realmat := Orthogonalize(realmat);
                    287:       res(i) := new Standard_Complex_Matrices.Matrix'(Convert(realmat));
                    288:     end loop;
                    289:     Close(planesfile);
                    290:     return res;
                    291:   end Read_Input_Planes;
                    292:
                    293:   function Read_Input_Planes ( m,p : natural; k : in Bracket ) return VecMat is
                    294:
                    295:     res : VecMat(k'range);
                    296:     planesfile : file_type;
                    297:     n : constant natural := m+p;
                    298:
                    299:   begin
                    300:     new_line;
                    301:     put_line("Reading the name of the file with the input planes.");
                    302:     Read_Name_and_Open_File(planesfile);
                    303:     for i in res'range loop
                    304:       declare
                    305:         realmat : Standard_Floating_Matrices.Matrix(1..n,1..m+1-k(i));
                    306:       begin
                    307:         get(planesfile,realmat);
                    308:         realmat := Orthogonalize(realmat);
                    309:         res(i) := new Standard_Complex_Matrices.Matrix'(Convert(realmat));
                    310:       end;
                    311:     end loop;
                    312:     Close(planesfile);
                    313:     return res;
                    314:   end Read_Input_Planes;
                    315:
                    316:   function Read_Input_Planes ( m,p,q : natural ) return VecMat is
                    317:
                    318:     res : VecMat(1..m*p+q*(m+p));
                    319:     planesfile : file_type;
                    320:     n : constant natural := m+p;
                    321:     realmat : Standard_Floating_Matrices.Matrix(1..n,1..m);
                    322:
                    323:   begin
                    324:     new_line;
                    325:     put_line("Reading the name of the file with the input planes.");
                    326:     Read_Name_and_Open_File(planesfile);
                    327:     for i in res'range loop
                    328:       get(planesfile,realmat);
                    329:       realmat := Orthogonalize(realmat);
                    330:       res(i) := new Standard_Complex_Matrices.Matrix'(Convert(realmat));
                    331:     end loop;
                    332:     Close(planesfile);
                    333:     return res;
                    334:   end Read_Input_Planes;
                    335:
                    336: -- MAIN INTERACTIVE DRIVERS :
                    337:
                    338:   function Select_Input_Choice return character is
                    339:
                    340:   -- DESCRIPTION :
                    341:   --   Displays the menu to obtain a choice for the kind of input.
                    342:
                    343:     res : character;
                    344:
                    345:   begin
                    346:     new_line;
                    347:     put_line("MENU for generating real input planes : ");
                    348:     put_line("  0. Random real planes at equidistant interpolation points.");
                    349:     put_line("  1. Generate input planes osculating a rational normal curve.");
                    350:     put_line("  2. Interactively give s-values for the "
                    351:                                & "osculating input planes.");
                    352:     put_line("  3. Give the name of the file with input planes.");
                    353:     put("Type 0, 1, 2, or 3 to select : "); Ask_Alternative(res,"0123");
                    354:     return res;
                    355:   end Select_Input_Choice;
                    356:
                    357:   procedure Driver_for_Input_Planes
                    358:               ( file : in file_type; m,p : in natural; planes : out VecMat ) is
                    359:
                    360:     input_choice : character := Select_Input_Choice;
                    361:     dim : constant natural := m*p;
                    362:     svals : Vector(1..dim);
                    363:
                    364:   begin
                    365:     case input_choice is
                    366:       when '0' => planes := Random_Real_Planes(m,p);
                    367:       when '1' => planes := Osculating_Input_Planes(m,p);
                    368:       when '2' => svals := Read_Interpolation_Points(dim);
                    369:                   planes := Osculating_Input_Planes(m,p,svals);
                    370:       when '3' => planes := Read_Input_Planes(m,p);
                    371:       when others => null;
                    372:     end case;
                    373:   end Driver_for_Input_Planes;
                    374:
                    375:   procedure Driver_for_Input_Planes
                    376:               ( file : in file_type; m,p : in natural; k : in Bracket;
                    377:                 planes : out VecMat ) is
                    378:
                    379:     input_choice : character := Select_Input_Choice;
                    380:     svals : Vector(k'range);
                    381:
                    382:   begin
                    383:     case input_choice is
                    384:       when '0' => planes := Random_Real_Planes(m,p,k);
                    385:       when '1' => planes := Osculating_Input_Planes(m,p,k);
                    386:       when '2' => svals := Read_Interpolation_Points(k'length);
                    387:                   planes := Osculating_Input_Planes(m,p,k,svals);
                    388:       when '3' => planes := Read_Input_Planes(m,p,k);
                    389:       when others => null;
                    390:     end case;
                    391:   end Driver_for_Input_Planes;
                    392:
                    393:   procedure Driver_for_Input_Planes
                    394:               ( file : in file_type; m,p,q : in natural;
                    395:                 s : out Vector; planes : out VecMat ) is
                    396:
                    397:     input_choice : character := Select_Input_Choice;
                    398:     dim : constant natural := m*p + q*(m+p);
                    399:     svals : Vector(1..dim);
                    400:
                    401:   begin
                    402:     case input_choice is
                    403:       when '0' => svals := Equidistant_Interpolation_Points(dim);
                    404:                   s := svals;
                    405:                   planes := Random_Real_Planes(m,p,q);
                    406:       when '1' => svals := Equidistant_Interpolation_Points(dim);
                    407:                   s := svals;
                    408:                   planes := Osculating_Input_Planes(m,p,q,svals);
                    409:       when '2' => svals := Read_Interpolation_Points(dim);
                    410:                   s := svals;
                    411:                   planes := Osculating_Input_Planes(m,p,q,svals);
                    412:       when '3' => svals := Read_Interpolation_Points(dim);
                    413:                   s := svals;
                    414:                   planes := Read_Input_Planes(m,p,q);
                    415:       when others => null;
                    416:     end case;
                    417:     put_line(file,"The interpolation points : ");
                    418:     put_line(file,svals);
                    419:     put_line(file,"The target planes : ");
                    420:     put(file,planes,3);
                    421:   end Driver_for_Input_Planes;
                    422:
                    423: end Drivers_for_Input_Planes;

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>