home *** CD-ROM | disk | FTP | other *** search
-
-
- /*
- YAMP - Yet Another Matrix Program
- Version: 1.6
- Author: Mark Von Tress, Ph.D.
- Date: 01/11/93
-
- Copyright(c) Mark Von Tress 1993
- 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.
-
- */
-
- #include <d:dos.h>
- #include <d:stdlib.h>
- #include <d:stdio.h>
- #include <d:math.h>
- #include <d:string.h>
- #include <d:io.h>
- #include <d:fcntl.h>
- #include <d:errno.h>
- #include <d:types.h>
- #include <d:stat.h>
- #include <d:iostream.h>
-
-
- //#define IN_RAM // for in ram version
-
- #define STACKLENGTH 16352
-
- #define NDECS 10
- #define ZERO (1.0E-8)
-
- #ifndef SIGNATURE
- #define SIGNATURE 0x5a5a
- #endif
-
-
- /////////////////////////////////////////// string type class
-
- #define MAXCHARS 80
-
-
- class strtype {
- private:
- char name[MAXCHARS];
- public:
- strtype( strtype &str );
- strtype( char *str );
- strtype( void );
- strtype operator+(strtype str);
- strtype operator+(const char *str);
- strtype operator=(strtype str);
- strtype operator=(const char *str);
- void Showname( void ) { printf("%s",name) ;}
- char *StringAddress( void ) { return name; }
- };
-
-
- ////////////////////////////////////////////////////////////////
-
-
-
- ////////////////////////////// inram version
- #ifdef IN_RAM
-
- #ifndef HUGE_POINTER
-
- // define HUGE_PTR for MSC7 to be __huge. The value __far
- // works for MSC7.01 which is a patch of version 7.
-
- #ifdef _MSC_VER
- #define HUGE_POINTER __far
- #endif
- // define huge for BC3.1
- #ifdef __BCPLUSPLUS__
- #define HUGE_POINTER far
- #endif
-
- // need to define this for other compilers for huge keyword
- #ifndef HUGE_POINTER
- #define HUGE_POINTER
- #endif
-
- #endif
-
- #include <alloc.h>
-
- // for compatability with large model.
- #ifndef VCLOSE
- #define VCLOSE 1
- #define vclose()
- #endif
-
- typedef struct vdoub{
- // use huge pointer to avoid memory leak in far pointer delete
- // in MSC7. It doesn't release big chunks correctly. huge
- // pointers work ok. 06/07/92
- int nrefs;
- double HUGE_POINTER *mm;
- } vdoub;
-
- class VMatrix
- {
-
- protected:
- long curvecind;
- vdoub *v;
- double HUGE_POINTER *mcheck;
- void SetupVectors( int rr=1, int cc=1 );
- void PurgeVectors( void );
- void Replace( VMatrix &ROp );
- void NewReference(VMatrix &ROp );
- long mindex( int i, int j){ return ( (i-1)*c+(j-1) ); }
- strtype name;
- int signature;
-
- public:
-
- double Nrerror( const char *errormsg );
- void Garbage( const char *errormsg);
- void Garbage( void ){ Garbage(" ");}
- int r,c;
-
- double m( int i, int j ) { // get a value
- #ifndef NO_CHECKING
- if( i<1 || j<1 || i>r || j>c ){
- cerr << "index out of range\n";
- exit(1);
- }
- #endif
- curvecind = mindex(i,j);
- return v->mm[ curvecind ];
- }
- VMatrix& M( int i, int j ) { // load current buffer
- curvecind = mindex(i,j);
- return *this;
- }
-
- VMatrix( void ); // constructors and destructors
- VMatrix( const char *str, int i, int j);
- VMatrix( VMatrix& ROp );
-
- ~VMatrix( void );
-
- double operator = ( double d ){ // write in a value
- this->v->mm[curvecind] = d;
- return d;
- }
- void operator= (VMatrix &ROp);
-
-
- void Nameit( const char *str ) { name = str; }
- void Nameit( VMatrix &mat ){ name = mat.name; }
- void Nameit( strtype newname){ name = newname; }
- strtype Getname( VMatrix &mat){ return mat.name; }
- void Showname( void ){ name.Showname(); }
- char *StringAddress( void ) { return name.StringAddress(); }
-
- void LoadMat(void);
- void DisplayMat(void);
- void InfoMat( void );
- void Writeb (char *fid, VMatrix &mat);
-
- friend VMatrix& operator+ (VMatrix &LOp, VMatrix &ROp);
- friend VMatrix& operator+ (double scalar, VMatrix &ROp);
- friend VMatrix& operator+ (VMatrix &ROp, double scalar);
- friend VMatrix& operator- (VMatrix &LOp, VMatrix &ROp);
- friend VMatrix& operator- (double scalar, VMatrix &ROp);
- friend VMatrix& operator- (VMatrix &ROp, double scalar);
- friend VMatrix& operator- (VMatrix &ROp);
- friend VMatrix& operator* (VMatrix &LOp, VMatrix &ROp);
- friend VMatrix& operator* (double scalar, VMatrix &ROp);
- friend VMatrix& operator* (VMatrix &ROp, double scalar);
- friend VMatrix& operator% (VMatrix &LOp, VMatrix &ROp );
- friend VMatrix& operator/ (VMatrix &LOp, VMatrix &ROp);
- friend VMatrix& operator/ (VMatrix &ROp, double scalar);
- friend VMatrix& operator/ (double scalar, VMatrix &ROp);
-
- };
-
-
- //////////////// matrix stack
-
- class MStack: public VMatrix {
-
- private:
-
- int stackloc, level;
- MStack *next;
- VMatrix& Pop( void );
-
- public:
-
- MStack( void );
- void Push( VMatrix &ROp );
-
- void Inclevel();
- void Declevel();
- VMatrix &DecReturn( void ){
- Declevel();
- next->level = level;
- return * ((VMatrix *) this->next);
- }
- VMatrix &ReturnMat( void ){
- next->level = level;
- return * ((VMatrix *) this->next);
- }
-
- void PrintStack( void );
- void Cleanstack( void );
-
- };
-
-
- #else
-
- ///////////////////////// large memory model
-
-
-
- //////////////////////////////////////// large memory model
-
-
- typedef struct vmem{
- unsigned signature;
- unsigned dirty;
- unsigned ele_size;
- unsigned ele_per_blk;
- long num_ele;
- unsigned cblock;
- char *buf;
- } vmem;
-
- typedef struct hdr {
- int active;
- int nrefs;
- unsigned nblocks;
- unsigned offset;
- union { struct hdr *h;
- struct vmem *v;
- } next;
- } hdr;
-
-
- int vopen( void );
- int vclose (void);
- hdr *vmalloc ( long num_ele, unsigned int ele_size);
- int vfree ( hdr *p);
-
- // the following 3 functions must be public if you want to
- // use the vmalloc system by itself. They are not needed globally in
- // the virtual matrix package.
-
- //void *vread (hdr *p, long index);
- //void *vwrite (hdr *p, long index, char *thiss);
- //long vele(hdr *p);
-
- class vdoub{
-
- protected:
-
- unsigned signature;// garbage flag;
- hdr *hdr; // virtual memory handle
- double *cur_ele; // current element
- long cur_index; // current index
- // static int nobj;
-
- friend double val( vdoub &v);
-
- public:
-
- int Garbage( void ){ return !( signature == SIGNATURE ); }
-
- vdoub ( long array_size );
- vdoub ( void );
- vdoub ( vdoub &src );
- ~vdoub( void );
-
- double v( long index ); // get value
- vdoub& operator[]( long index ); // store value
- double operator = ( double d );
- double operator = ( vdoub &v );
- };
-
- ///////////////////////////////////////// virtual matrix object
-
- class VMatrix : private vdoub
- {
-
- protected:
- unsigned signature;
- long mindex( int i, int j){
- return (long) (((long)(i-1))*((long)c)+((long)(j-1))) ; }
- long curvecind; // current index for vector
- strtype name;
- void SetupVectors( int rr=1, int cc=1 );
- void PurgeVectors( void );
- void Replace( VMatrix &ROp );
- void NewReference( VMatrix &ROp );
-
- public:
-
- double Nrerror( const char *errormsg );
- void Garbage( const char *errormsg);
- void Garbage( void ){ Garbage(" ");}
- int r,c;
-
- double m( int i, int j ) { // get a value
- #ifndef NO_CHECKING
- if( i<1 || j<1 || i>r || j>c ){
- cerr << "index out of range\n";
- exit(1);
- }
- #endif
- curvecind = mindex(i,j);
- return v( curvecind );
- }
- VMatrix& M( int i, int j ) { // load current buffer
- curvecind = mindex(i,j);
- vdoub::operator[](curvecind);
- return *this;
- }
-
- VMatrix( void ); // constructors and destructors
- VMatrix( const char *str, int i, int j);
- VMatrix( VMatrix& ROp);
-
- ~VMatrix( void );
-
- double operator = ( double d ){ // write in a value
- vdoub::operator[](curvecind) = d;
- return d;
- }
- void operator= (VMatrix &ROp);
-
-
- void Nameit( const char *str ) { name = str; }
- void Nameit( VMatrix &mat ){ name = mat.name; }
- void Nameit( strtype newname){ name = newname; }
- strtype Getname( VMatrix &mat){ return mat.name; }
- void Showname( void ){ name.Showname(); }
- char *StringAddress( void ) { return name.StringAddress(); }
-
- void LoadMat(void);
- void DisplayMat(void);
- void InfoMat( void );
- void Writeb (char *fid, VMatrix &mat);
-
- friend VMatrix& operator+ (VMatrix &LOp, VMatrix &ROp);
- friend VMatrix& operator+ (double scalar, VMatrix &ROp);
- friend VMatrix& operator+ (VMatrix &ROp, double scalar);
- friend VMatrix& operator- (VMatrix &LOp, VMatrix &ROp);
- friend VMatrix& operator- (double scalar, VMatrix &ROp);
- friend VMatrix& operator- (VMatrix &ROp, double scalar);
- friend VMatrix& operator- (VMatrix &ROp);
- friend VMatrix& operator* (VMatrix &LOp, VMatrix &ROp);
- friend VMatrix& operator* (double scalar, VMatrix &ROp);
- friend VMatrix& operator* (VMatrix &ROp, double scalar);
- friend VMatrix& operator% (VMatrix &LOp, VMatrix &ROp );
- friend VMatrix& operator/ (VMatrix &LOp, VMatrix &ROp);
- friend VMatrix& operator/ (VMatrix &ROp, double scalar);
- friend VMatrix& operator/ (double scalar, VMatrix &ROp);
- };
-
-
- //////////////// matrix stack
-
- class MStack: public VMatrix {
-
- private:
-
- int stackloc, level;
- MStack *next;
- VMatrix& Pop( void );
-
- public:
-
- MStack( void );
- void Push( VMatrix &ROp );
-
- void Inclevel();
- void Declevel();
- VMatrix &DecReturn( void ){
- Declevel();
- next->level = level;
- return * ((VMatrix *) this->next);
- }
- VMatrix &ReturnMat( void ){
- next->level = level;
- return * ((VMatrix *) this->next);
- }
-
- void PrintStack( void );
- void Cleanstack( void );
-
- };
-
- ///////////////////////////////////////////// end of large memory model
-
- #endif
-
- extern MStack *Dispatch;
-
- typedef enum boolean { FFALSE, TTRUE } boolean;
-
- ///////////////////////////// display control
-
- int Setwid( int w ); // set display width
- int Setdec( int d ); // set decimals
- int Getwid( void ); // report display width
- int Getdec( void ); // report decimals
-
- ////////////////////////// io
-
- VMatrix& Reada( char *fid); // read ascii file
- VMatrix& Readb( char *fid); // read binary file
- void Writea( char *fid, VMatrix &mat); // write ascii file
- //void VMatrix::Writeb( char *fid, VMatrix &mat) // write binary file
-
- /////////////////////////// sorting, subsetting, and conatenation
- // sort order == 0 => sort col c and exchange rows
- // order != 0 => sort row c and exchange cols
- //
- // submat r,c is bottom right corner, lr,lc is upper left corner
- //
- // Ch horozontal concatenation
- // Cv vertical concatentation
- //
-
- VMatrix& MSort( VMatrix &a, int col, int order=0);
- VMatrix& Submat( VMatrix &a, int r, int c, int lr=1, int lc=1);
- VMatrix& Ch( VMatrix &a, VMatrix &b );
- VMatrix& Cv( VMatrix &a, VMatrix &b );
-
- /////////////////////////// special matrices
- //
- // Tran transpose of matrix
- // Mexp elementwise exponential
- // Mabs elementwise absolute value
- // Mlog elementwise log
- // Msqrt elementwise sqrt
- // Trace trace of square matrix
- // Ident identity matrix
- // Helm Helmert matrix
- // Fill constant matrix
- // Index build an index matrix: if lower>upper step downward
- //
-
- VMatrix& Tran(VMatrix &ROp );
- VMatrix& Mexp(VMatrix &ROp );
- VMatrix& Mabs(VMatrix &ROp );
- VMatrix& Mlog(VMatrix &ROp );
- VMatrix& Msqrt(VMatrix &ROp );
- double Trace(VMatrix &ROp );
- VMatrix& Ident(int n );
- VMatrix& Helm( int n );
- VMatrix& Fill(int r, int c, double d );
- VMatrix& Index( int lower, int upper, int step = 1);
-
- /////////////////////////// inversion and decomposition
- //
- // Sweep sweep on diagonal elements k1 to k2
- // Inv matrix inverse using sweep(1,ROp.c,ROp)
- // Kron Kroniker (direct) product
- // Det Determinant
- // Cholesky upper triangular S where ROp = S'S, for symmetric ROp
- // QR ROp( rxc, r>c ) = Q(rxc) R(cxc), and Q'Q = Ident(c)
- // makeq = TTRUE => make Q, it is garbage otherwise
- // Svd ROp = USTran(V),
- // where U'U = I, V'V=I, S= Column of singular values
- // makeu = TTrue => Compute U, it is garbage otherwise
- // makev = TTrue => Compute V, it is garbage otherwise
- // Ginv See Svd, Ginv(ROp) = V Inv(Diag(S+)) Tran(U)
- // Fft Almost fast Fourier transform, can be slow if
- // n > 1000 or prime factors of n > 19, does not require
- // n to be a power of 2
- // INZEE = TTRUE => compute fft
- // INZEE = FFALSE=> compute ifft
- //
-
- VMatrix& Sweep( int k1, int k2, VMatrix& ROp); /* sweep out a row */
- VMatrix& Inv( VMatrix &ROp ); /*matrix inversion using sweep */
- VMatrix& Kron( VMatrix &a, VMatrix &b );
- double Det( VMatrix &ROp );
- VMatrix& Cholesky(VMatrix& ROp);
- void QR( VMatrix& ROp, VMatrix& Q, VMatrix& R, boolean makeq = TTRUE);
- int Svd( VMatrix &A, VMatrix &U, VMatrix &S, VMatrix &V,
- boolean makeu = TTRUE, boolean makev = TTRUE);
- VMatrix& Ginv( VMatrix& ROp );
- VMatrix& Fft( VMatrix &ROp, boolean INZEE = TTRUE);
-
-
- ////////////////////////// reshaping functions
-
- VMatrix& Vec( VMatrix& ROp ); // send matrix to vector
- // stack columns into an (r*c)x1
- VMatrix& Vecdiag( VMatrix& ROp ); // extract diagonal into a vector
- VMatrix& Diag( VMatrix& ROp ); // create a diagonal matrix from
- // square matrix or column vector
- VMatrix& Shape( VMatrix& ROp, int nrows = 1 ); // reshape into nrows
-
- ////////////////////////// summation functions
-
- typedef enum margins { ALL, ROWS, COLUMNS } margins; // summation margins
-
- VMatrix& Sum( VMatrix& ROp, margins marg = ALL ); // sums over margins
- VMatrix& Sumsq( VMatrix& ROp, margins marg = ALL ); // sum of squares
- VMatrix& Cusum( VMatrix& ROp); // accumulate along rows
-
-
- /////////////////////////// max and mins
- //
- // if marg = ALL return 3x1 [ maxr, maxc, max element ]
- // if marg = ROWS return 3xc Tran([ maxr, col=1..c, max element ])
- // if marg = COLUMNS return rx3 [ row=1..r, maxc, max element]
- //
- // same structure for mins
- //
-
- VMatrix& Mmax( VMatrix& ROp, margins marg=ALL ); // matrix max
- VMatrix& Mmin( VMatrix& ROp, margins marg=ALL ); // matrix min
-
- ////////////////////////// row and column operations
- //
- // Crow c*r1
- // Srow swap r1 with r2
- // Lrow r2 + c*r1
- // Ccol c*c1
- // Scol swap c1 with c2
- // Lcol c2 + c*c1
- //
- // note default values for row and col are out of range, so they
- // must be supplied. This is for error checking.
- // default constants perform no changes.
- // !!!!!!!!!!!!! these functions operate directly on matrix
- //////////////////////////
-
- void Crow( VMatrix& ROp, int row=0, double c=1.0);
- void Srow( VMatrix &ROp, int row1 = 0, int row2 = 0);
- void Lrow( VMatrix &ROp, int row1 = 0, int row2 = 0, double c=0.0);
- void Ccol( VMatrix& ROp, int col=0, double c=1.0);
- void Scol( VMatrix &ROp, int col1 = 0, int col2 = 0);
- void Lcol( VMatrix &ROp, int col1 = 0, int col2 = 0, double c=0.0);
-
- //////////////////////////////////////////////////////////////
- //////////////////////////// XYplot objects ////////////////
- //////////////////////////////////////////////////////////////
-
-
- #ifdef VIRTGRAF
-
- #if !(__BCPLUSPLUS__ || __TCPLUSPLUS__)
- #error GMatrix class requires Borland BGI functions
- #endif
-
-
- #ifdef __TINY__
- #error GMatrix class will not work in the tiny model.
- #endif
-
- #include <conio.h>
- #include <stdarg.h>
- #include <graphics.h>
-
- // the Axis class is called from the GMatrix Class. It
- // is not needed otherwise, and will not be commented
- class Axis
- {
- protected:
- // uses fortran translations of Algorithm AS 168 of Applied Statistics
- double valmin, step;
- double *vals, fmn, fmx, offset;
- int nvals, irprin, maxpr, ifact;
- VMatrix *Vals;
- int nplt, mpv, ir, ifault , iv;
- double dmax( double x, double y);
- void scale(double fmn, double fmx, int nplt, int mpv,
- double *valmin, double *step, int *nvals,int *ir, int *ifault);
- void axis(double valmin,double step,int nvals, int maxpr,int ir,
- int *irprin, double *offset, int *ifact, double *vals, int iv,
- int * ifault);
- void getvals(double fmn,double fmx,int n, double *vals,
- double *offset,int *ifact,int *nvals,int *irprin);
- struct linesettingstype lineinfo;
-
- public:
- Axis( VMatrix &grf_vec, int col, int pixels, int pxlspplot);
- ~Axis( void );
-
- virtual void Show( void )
- {
- restorecrtmode();
- Dispatch->Nrerror("AXIS: trying to show an undefined axis object");
- }
- strtype GetAxisLabel( strtype &name );
- double Rescale( double x );
-
- };
- class YAxis:public Axis
- {
- private:
- double deltay;
- int h, w, xlen, ylen, bmarg, umarg, lmarg, rmarg;
-
- public:
- double miny, maxy;
-
- void Show( boolean gridon );
- YAxis( VMatrix &grf, int xxlen, int yylen, int ww, int hh,
- int bbmarg, int uumarg, int llmarg, int rrmarg);
- };
-
- class XAxis:public Axis
- {
- private:
- double deltax;
- int w, h, xlen, ylen, bmarg, umarg, lmarg, rmarg;
-
- public:
- double minx, maxx;
-
- void Show( boolean gridon );
- XAxis( VMatrix &grf, int xxlen, int yylen, int ww, int hh,
- int bbmarg, int uumarg, int llmarg, int rrmarg);
- };
-
-
- /////////////////////////////// graph matrix class
-
- class GMatrix
- {
- private:
- // Borland graphics information. bgi files must be available
- int GraphDriver;
- int GraphMode;
- int MaxX, MaxY;
- int MaxColors;
- int ErrorCode;
- struct linesettingstype lineinfo;
- struct viewporttype viewinfo;
-
-
- // reference lines. can have at most 20 reference lines
- double *hrefs;
- double *vrefs;
- int nvrefs, nhrefs;
-
- // output a line, and pause function ESC stops program
- int gprintf(int *xloc, int *yloc, char *fmt, ... );
- void Pause();
-
- // graph is stored in a VMatrix: [graphid, symbol, x, y] 4 cols
- int graphnum; // number of graphs in grf;
- VMatrix *grf;
-
- // plot a point on the screen
- void plotpoint( int ix, int iy, int symbol);
-
- public:
-
- // constructors and destructors
-
- GMatrix( void );
- GMatrix( VMatrix &ROp, int symbol= -' '); // ROp should be rx2 matrix
- // column 1 is x, col 2 is y
- ~GMatrix( void );
-
- // display a graph
- void Show( void );
-
- // add a new vector to plot
- void AddVec( VMatrix &ROp, int symbol = -' ');
-
- // annotation options
-
- strtype *title, *title2, *xname, *yname, *PathToDriver;
- boolean RectangleOn, YGridOn, XGridOn;
- void Href( double href );
- void Vref( double vref );
- int ClearHref( void ){ int t=nhrefs; nhrefs = -1; return t;}
- int ClearVref( void ){ int t=nvrefs; nvrefs = -1; return t;}
-
- // save the grf matrix
- void SaveGraph( char *fid ) { Writea( fid, *grf ); }
- };
-
- #endif
-