home *** CD-ROM | disk | FTP | other *** search
/ Stars of Shareware: Programmierung / SOURCE.mdf / programm / msdos / c / yamp / test.cpp < prev    next >
Encoding:
C/C++ Source or Header  |  1993-01-11  |  6.3 KB  |  225 lines
  1.  
  2. /*
  3. YAMP - Yet Another Matrix Program
  4. Version: 1.6
  5. Author: Mark Von Tress, Ph.D.
  6. Date: 01/11/93  
  7.  
  8. Copyright(c) Mark Von Tress 1993
  9. Portions of this code are (c) 1991 by Allen I. Holub and are used by
  10. permission.
  11.  
  12. DISCLAIMER: THIS PROGRAM IS PROVIDED AS IS, WITHOUT ANY
  13. WARRANTY, EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED
  14. TO FITNESS FOR A PARTICULAR PURPOSE. THE AUTHOR DISCLAIMS
  15. ALL LIABILITY FOR DIRECT OR CONSEQUENTIAL DAMAGES RESULTING
  16. FROM USE OF THIS PROGRAM.
  17.  
  18. */
  19. // make sure to define the global define variable  IN_RAM if you want
  20. // to use the medium model. This can be done in the compiler code
  21. // generation window of compiler options. It may also be defined for
  22. // the command line compiler
  23.  
  24. #include "virt.h"
  25.  
  26. /* required global declaration for the matrix stack object */
  27.  
  28. //unsigned int _stklen = STACKLENGTH;
  29.  
  30. MStack *Dispatch = new MStack;
  31.  
  32.  
  33. VMatrix& function1(VMatrix &a, VMatrix &b)
  34. {
  35.     // This function tests the freind functions and returns a value
  36.     
  37.     a.Garbage("function1");     // check a and b
  38.     b.Garbage("function1");
  39.     
  40.     Dispatch->Inclevel();     // increment push-pop level
  41.     VMatrix c("c", 1, 1);     // create a local matrix
  42.     
  43.     c = a + b + a + b;     // check a repeated matrix addition
  44.     c.DisplayMat();     // print c
  45.     c = 431.2 + Tran(a) *b + 2.134;     // check commutivity of scalar addition
  46.     c.DisplayMat();     // print c
  47.     c = -c;             // check uniary minus
  48.     c.DisplayMat();     // print c
  49.     c = a - b;             // matrix subraction
  50.     c.DisplayMat();     // print c
  51.     c = 5 - a - 5;     // check commutivity of scalar subtraction
  52.     c.DisplayMat();     // print c
  53.     c = 5*a*5;             // check commutivity of scalar multiplication
  54.     c.DisplayMat();     // print c
  55.     c = a % a;             // check elementwise multiplication
  56.     c.DisplayMat();     // print c
  57.     c = a / 1234;         // check scalar division
  58.     c.DisplayMat();     // print c
  59.     c = a / b;             // check elementwise division
  60.     c.DisplayMat();     // print c
  61.  
  62.     Dispatch->PrintStack();       // show stack before push
  63.     Dispatch->Push(c);               // push c onto stack
  64.     Dispatch->PrintStack();       // examine stack after push
  65.     return Dispatch->DecReturn(); // decrement push-pop level, and return
  66.                                   // stack top
  67. }
  68.  
  69. VMatrix &function0(void)
  70. {
  71.     // test some of the output functions and raw matrix functions
  72.  
  73.     Dispatch->Inclevel();
  74.     VMatrix d = VMatrix("d", 1, 1);
  75.     VMatrix a, c, H, hilb;
  76.  
  77.     a = Reada("catchv.dat");     // read an ascii matrix
  78.     a.DisplayMat();             // display a
  79.     Writea("junk.dat", a);         // write ascii matrix
  80.     a.Writeb("junk.bin", a);     // write binary matrix
  81.     a = Readb("junk.bin");         // read binary matrix
  82.     a.InfoMat();                 // display matrix info
  83.     a.DisplayMat();             // display a
  84.     d = Submat(a, a.r, 4, 1, 2);     // take a submatrix of a
  85.     d.DisplayMat();             // display it
  86.     c = Ch(d, d);                 // horizontal concatenation
  87.     c.DisplayMat();
  88.     c = Cv(d, d);                 // vertical contatenation
  89.     c.DisplayMat();
  90.     c = Kron(Ident(3), d);         // Kroniker's product
  91.     Setdec(1);                     // set number of decimals to print
  92.     Setwid(5);                     // set print width
  93.     c.DisplayMat();
  94.     H = Helm(4);                 // make a helmert matrix
  95.     H.DisplayMat();             // show that it is orthonormal
  96.     d = Tran(H) *H;
  97.     d.DisplayMat();
  98.     d = H*Tran(H);
  99.     d.DisplayMat();
  100.  
  101.     a = Ident(4) + Fill(4, 4, 0.5);     // redefine a
  102.     a.DisplayMat();                     // print a
  103.     c = Tran(Helm(4)) *a*H;             // diagonalize a
  104.     c.DisplayMat();                     // display c
  105.     c = Inv(a);                         // invert a
  106.     c.DisplayMat();                     // display c
  107.     VMatrix b;                             // create b
  108.     b = c*a;                             // redefine b
  109.     b.DisplayMat();                      // display b
  110.     c = function1(a, a);                 // call a function
  111.     c.DisplayMat();                     // display c
  112.     Dispatch->Push(c);
  113.     return Dispatch->DecReturn();
  114. }
  115.  
  116. VMatrix ®ression(void)     // do a multiple linear regression
  117. {
  118.     Dispatch->Inclevel();
  119.     VMatrix a, xy, reg;
  120.     a = Reada("catchv.dat");                 // read data
  121.     VMatrix b = VMatrix("b", a.r, a.c);     //simplify indexes
  122.     int N = a.r;                             // N
  123.     int p = 3;                                 // params
  124.     xy = Ch(Fill(N, 1, 1), Submat(a, N, 4, 1, 2));     // make x y
  125.     reg = Sweep(1, p, Tran(xy) *xy);         // solve regression
  126.                                             // using sweep
  127.     reg.M(p + 1, p + 1) = reg.m(p + 1, p + 1) / ((double) (N - p));
  128.     // divide to get mse
  129.  
  130.     VMatrix x = Submat(xy, N, 3, 1, 1);
  131.     VMatrix y = Submat(xy, N, 4, 1, 4);
  132.     VMatrix betahat = Inv(Tran(x) *x) *Tran(x) *y;
  133.                                            //solve regression using
  134.                                            //normal equations
  135.     betahat.DisplayMat();
  136.     reg.DisplayMat();
  137.     Dispatch->Push(reg);                    // put return mat on stack
  138.     return Dispatch->DecReturn();            // dec level & return
  139. }
  140.  
  141. void altermatrix(VMatrix &t)     // alter a matrix element
  142. {
  143.     t.M(1, 1) = 2525.0;
  144.     VMatrix temp = MSort(t, 1);
  145.     temp.DisplayMat();
  146. }
  147.  
  148. void testhuge(void)
  149. {
  150.     //VMatrix VeryBig("vb", 1000, 1000);
  151.         cout << "Inverting " << endl;
  152.         VMatrix VeryBig = Inv( Helm(1000 ) );
  153.     VeryBig.InfoMat();
  154.  
  155. }
  156.  
  157. void version1p1(void)
  158. {
  159.     //
  160.     // testreg.cpp for test of regression functions.
  161.     //
  162.     Dispatch->Inclevel();
  163.     VMatrix A = Fill(6, 6, 1.0);
  164.     VMatrix B = Fill(6, 1, 2.0);
  165.  
  166.     (Vec(A)).DisplayMat();
  167.     (Vecdiag(A)).DisplayMat();
  168.     (Diag(B)).DisplayMat();
  169.     (Shape(A, 3)).DisplayMat();
  170.     (Sum(A, ROWS)).DisplayMat();
  171.     (Sumsq(A, COLUMNS)).DisplayMat();
  172.     (Cusum(A)).DisplayMat();
  173.     (Mmin(B)).DisplayMat();
  174.     (Mmax(B, ROWS)).DisplayMat();
  175.  
  176.     VMatrix C = Ch(A, Diag(B));
  177.     Setwid(5);
  178.     Setdec(2);
  179.     C.DisplayMat();
  180.     Crow(C, 1, 0.6);
  181.     C.DisplayMat();
  182.     Srow(C, 1, 6);
  183.     C.DisplayMat();
  184.     Lrow(C, 2, 3,.5);
  185.     C.DisplayMat();
  186.     Ccol(C, 1, 0.6);
  187.     C.DisplayMat();
  188.     Scol(C, 1, 6);
  189.     C.DisplayMat();
  190.     Lcol(C, 2, 3,.5);
  191.     C.DisplayMat();
  192.     Dispatch->Cleanstack();
  193.     Dispatch->Declevel();
  194. }
  195.  
  196. main()
  197. {
  198.     VMatrix testit;
  199.  
  200.     printf("%d\n", sizeof (testit));
  201.  
  202.     testit = regression();
  203.     testit.DisplayMat();
  204.  
  205.     VMatrix b = testit;
  206.     VMatrix a = 2*testit;
  207.  
  208.     VMatrix c = a + b + a + b;
  209.     c.DisplayMat();
  210.  
  211.     VMatrix t = function0() + function0();
  212.     t.DisplayMat();
  213.  
  214.     altermatrix(testit);
  215.     testit.DisplayMat();
  216.  
  217.     version1p1();
  218.  
  219.     #ifndef IN_RAM
  220.     testhuge();
  221.     #endif
  222.  
  223.     vclose();
  224. }
  225.