/*
YAMP - Yet Another Matrix Program
Version: 1.5
Author: Mark Von Tress, Ph.D.
Date: 08/28/92

Portions of this code are (c) 1991 by Allen I. Holub and are used by
permission.

DISCLAIMER: THIS PROGRAM IS PROVIDED AS IS, WITHOUT ANY
WARRANTY, EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED
TO FITNESS FOR A PARTICULAR PURPOSE. THE AUTHOR DISCLAIMS
ALL LIABILITY FOR DIRECT OR CONSEQUENTIAL DAMAGES RESULTING
FROM USE OF THIS PROGRAM.

*/
// make sure to define the global define variable IN_RAM if you want
// to use the medium model. This can be done in the compiler code
// generation window of compiler options. It may also be defined for
// the command line compiler

#include "virt.h"

/* required global declaration for the matrix stack object */

//unsigned int _stklen = STACKLENGTH;

MStack *Dispatch = new MStack;

VMatrix& function1(VMatrix &a, VMatrix &b)
{
// This function tests the freind functions and returns a value

a.Garbage("function1"); // check a and b
b.Garbage("function1");

Dispatch->Inclevel(); // increment push-pop level
VMatrix c("c", 1, 1); // create a local matrix

c = a + b + a + b; // check a repeated matrix addition
c.DisplayMat(); // print c
c = 431.2 + Tran(a) *b + 2.134; // check commutivity of scalar addition
c.DisplayMat(); // print c
c = -c; // check uniary minus
c.DisplayMat(); // print c
c = a - b; // matrix subraction
c.DisplayMat(); // print c
c = 5 - a - 5; // check commutivity of scalar subtraction
c.DisplayMat(); // print c
c = 5*a*5; // check commutivity of scalar multiplication
c.DisplayMat(); // print c
c = a % a; // check elementwise multiplication
c.DisplayMat(); // print c
c = a / 1234; // check scalar division
c.DisplayMat(); // print c
c = a / b; // check elementwise division
c.DisplayMat(); // print c

Dispatch->PrintStack(); // show stack before push
Dispatch->Push(c); // push c onto stack
Dispatch->PrintStack(); // examine stack after push
return Dispatch->DecReturn(); // decrement push-pop level, and return
// stack top
}

VMatrix &function0(void)
{
// test some of the output functions and raw matrix functions

Dispatch->Inclevel();
VMatrix d = VMatrix("d", 1, 1);
VMatrix a, c, H, hilb;

a.DisplayMat(); // display a
Writea("junk.dat", a); // write ascii matrix
a.Writeb("junk.bin", a); // write binary matrix
a.InfoMat(); // display matrix info
a.DisplayMat(); // display a
d = Submat(a, a.r, 4, 1, 2); // take a submatrix of a
d.DisplayMat(); // display it
c = Ch(d, d); // horizontal concatenation
c.DisplayMat();
c = Cv(d, d); // vertical contatenation
c.DisplayMat();
c = Kron(Ident(3), d); // Kroniker's product
Setdec(1); // set number of decimals to print
Setwid(5); // set print width
c.DisplayMat();
H = Helm(4); // make a helmert matrix
H.DisplayMat(); // show that it is orthonormal
d = Tran(H) *H;
d.DisplayMat();
d = H*Tran(H);
d.DisplayMat();

a = Ident(4) + Fill(4, 4, 0.5); // redefine a
a.DisplayMat(); // print a
c = Tran(Helm(4)) *a*H; // diagonalize a
c.DisplayMat(); // display c
c = Inv(a); // invert a
c.DisplayMat(); // display c
VMatrix b; // create b
b = c*a; // redefine b
b.DisplayMat(); // display b
c = function1(a, a); // call a function
c.DisplayMat(); // display c
Dispatch->Push(c);
return Dispatch->DecReturn();
}

VMatrix ®ression(void) // do a multiple linear regression
{
Dispatch->Inclevel();
VMatrix a, xy, reg;
VMatrix b = VMatrix("b", a.r, a.c); //simplify indexes
int N = a.r; // N
int p = 3; // params
xy = Ch(Fill(N, 1, 1), Submat(a, N, 4, 1, 2)); // make x y
reg = Sweep(1, p, Tran(xy) *xy); // solve regression
// using sweep
reg.M(p + 1, p + 1) = reg.m(p + 1, p + 1) / ((double) (N - p));
// divide to get mse

VMatrix x = Submat(xy, N, 3, 1, 1);
VMatrix y = Submat(xy, N, 4, 1, 4);
VMatrix betahat = Inv(Tran(x) *x) *Tran(x) *y;
//solve regression using
//normal equations
betahat.DisplayMat();
reg.DisplayMat();
Dispatch->Push(reg); // put return mat on stack
return Dispatch->DecReturn(); // dec level & return
}

void altermatrix(VMatrix &t) // alter a matrix element
{
t.M(1, 1) = 2525.0;
VMatrix temp = MSort(t, 1);
temp.DisplayMat();
}

void testhuge(void)
{
VMatrix VeryBig("vb", 300, 300);
VeryBig.InfoMat();

}

void version1p1(void)
{
//
// testreg.cpp for test of regression functions.
//
Dispatch->Inclevel();
VMatrix A = Fill(6, 6, 1.0);
VMatrix B = Fill(6, 1, 2.0);

(Vec(A)).DisplayMat();
(Vecdiag(A)).DisplayMat();
(Diag(B)).DisplayMat();
(Shape(A, 3)).DisplayMat();
(Sum(A, ROWS)).DisplayMat();
(Sumsq(A, COLUMNS)).DisplayMat();
(Cusum(A)).DisplayMat();
(Mmin(B)).DisplayMat();
(Mmax(B, ROWS)).DisplayMat();

VMatrix C = Ch(A, Diag(B));
Setwid(5);
Setdec(2);
C.DisplayMat();
Crow(C, 1, 0.6);
C.DisplayMat();
Srow(C, 1, 6);
C.DisplayMat();
Lrow(C, 2, 3,.5);
C.DisplayMat();
Ccol(C, 1, 0.6);
C.DisplayMat();
Scol(C, 1, 6);
C.DisplayMat();
Lcol(C, 2, 3,.5);
C.DisplayMat();
Dispatch->Cleanstack();
Dispatch->Declevel();
}

main()
{
VMatrix testit;

printf("%d\n", sizeof (testit));

testit = regression();
testit.DisplayMat();

VMatrix b = testit;
VMatrix a = 2*testit;

VMatrix c = a + b + a + b;
c.DisplayMat();

VMatrix t = function0() + function0();
t.DisplayMat();

altermatrix(testit);
testit.DisplayMat();

version1p1();

#ifndef IN_RAM
testhuge();
#endif

vclose();
}