1 #ifndef MATRIX_H
   2 #define MATRIX_H
   3 
   4 #ifdef X
   5 #undef X
   6 #endif
   7 
   8 #include <stdio.h>
   9 
  10 /* matrices and vectors are the same data type, vector will only
  11    have one column (by default) or one row.
  12 */
  13 typedef struct
  14 {
  15   short   type ;
  16   int     rows ;
  17   int     cols ;
  18   float **rptr;    /* pointer to an array of rows */
  19   float *data;     /* pointer to base of data */
  20   FILE *mmapfile;
  21 } MATRIX, VECTOR ;
  22 
  23 typedef struct
  24 {
  25   float  real ;
  26   float  imag ;
  27 } COMPLEX_FLOAT, *CPTR ;
  28 
  29 #define MATRIX_CELT(m,r,c)      (((COMPLEX_FLOAT **)m->rptr)[r]+c)
  30 #define MATRIX_RELT(m,r,c)      (m->rptr[r]+c)
  31 #define MATRIX_ELT(m,r,c)       (m->type == MATRIX_REAL ? \
$^
                                  *MATRIX_RELT(m,r,c) : \
  32                                   *MATRIX_CELT(m,r,c))
  33 #define MATRIX_PTR(m,r,c)       (m->type == MATRIX_REAL ? \
$^
                                  MATRIX_RELT(m,r,c) : \
  34                                   (float *)MATRIX_CELT(m,r,c))
  35 
  36 #define MATRIX_CELT_REAL(m,r,c)  (MATRIX_CELT(m,r,c)->real)
  37 #define MATRIX_CELT_IMAG(m,r,c)  (MATRIX_CELT(m,r,c)->imag)
  38 
  39 #define MATRIX_REAL        1
  40 #define MATRIX_COMPLEX     2
  41 #define MATRIX_SYM    0
  42 #define MATRIX_UPPER  1
  43 #define MATRIX_LOWER  2
  44 
  45 MATRIX  *MatrixReshape(MATRIX *m_src, MATRIX *m_dst, int rows, int cols) ;
  46 int     MatrixCheck(MATRIX *m) ;
  47 MATRIX  *MatrixInverse(MATRIX *mIn, MATRIX *mOut) ;
  48 MATRIX  *MatrixPseudoInverse(MATRIX *m, MATRIX *m_pseudo_inv) ;
  49 MATRIX  *MatrixSVDPseudoInverse(MATRIX *m, MATRIX *m_pseudo_inv) ;
  50 MATRIX  *MatrixRightPseudoInverse(MATRIX *m, MATRIX *m_pseudo_inv) ;
  51 #define MatrixLeftPseudoInverse MatrixPseudoInverse
  52 MATRIX  *MatrixAlloc(int rows, int cols, int type) ;
  53 int     MatrixFree(MATRIX **pmat) ;
  54 MATRIX  *MatrixMultiply(MATRIX *m1, MATRIX *m2, MATRIX *m3) ;
  55 MATRIX  *MatrixCopy(MATRIX *mIn, MATRIX *mOut) ;
  56 int     MatrixWriteTxt(char *fname, MATRIX *mat) ;
  57 MATRIX  *MatrixReadTxt(char *fname, MATRIX *mat) ;
  58 MATRIX  *MatrixRead(char *fname) ;
  59 int     MatrixWrite(MATRIX *mIn, char *fname, char *name) ;
  60 MATRIX  *MatrixIdentity(int n, MATRIX *mI) ;
  61 int     MatrixPrint(FILE *fp, MATRIX *mat) ;
  62 int     MatrixPrintOneLine(FILE *fp, MATRIX *mat) ;
  63 int     MatrixPrintTranspose(FILE *fp, MATRIX *mat) ;
  64 MATRIX  *MatrixTranspose(MATRIX *mIn, MATRIX *mOut) ;
  65 MATRIX  *MatrixAdd(MATRIX *m1, MATRIX *m2, MATRIX *mOut) ;
  66 MATRIX  *MatrixSubtract(MATRIX *m1, MATRIX *m2, MATRIX *mOut) ;
  67 MATRIX  *MatrixScalarMul(MATRIX *mIn, float val, MATRIX *mOut) ;
  68 MATRIX  *MatrixClear(MATRIX *mat) ;
  69 MATRIX  *MatrixSquareElts(MATRIX *mIn, MATRIX *mOut) ;
  70 MATRIX  *MatrixSignedSquareElts(MATRIX *mIn, MATRIX *mOut) ;
  71 MATRIX  *MatrixSqrtElts(MATRIX *mIn, MATRIX *mOut) ;
  72 MATRIX  *MatrixDiag(MATRIX *mDiag, MATRIX *mOut) ;
  73 MATRIX  *MatrixMakeDiagonal(MATRIX *mSrc, MATRIX *mDst) ;
  74 MATRIX  *MatrixCopyRegion(MATRIX *mSrc, MATRIX *mDst, int start_row,
  75                           int start_col, int rows, int cols,
  76                           int dest_row, int dest_col) ;
  77 MATRIX  *MatrixCopyRealRegion(MATRIX *mSrc, MATRIX *mDst,int start_row,
  78                               int start_col, int rows, int cols,
  79                               int dest_row, int dest_col) ;
  80 MATRIX  *MatrixCopyImagRegion(MATRIX *mSrc, MATRIX *mDst, int start_row,
  81                               int start_col, int rows, int cols,
  82                               int dest_row, int dest_col) ;
  83 MATRIX  *MatrixSetRegion(MATRIX *mSrc, MATRIX *mDst, int start_row,
  84                          int start_col, int rows, int cols, float val);
  85 MATRIX *MatrixRealToComplex(MATRIX *mReal, MATRIX *mImag, MATRIX *mOut);
  86 MATRIX *MatrixRegularize(MATRIX *mIn, MATRIX *mOut) ;
  87 int    MatrixSingular(MATRIX *m) ;
  88 MATRIX *MatrixToeplitz(VECTOR *v, MATRIX *T, int Type);
  89 
  90 /* determinants and eigenvectors */
  91 float  MatrixDeterminant(MATRIX *m) ;
  92 MATRIX *MatrixEigenSystem(MATRIX *m, float *evalues, MATRIX *m_dst) ;
  93 MATRIX *MatrixSVD(MATRIX *mA, VECTOR *v_z, MATRIX *mV) ;
  94 MATRIX *MatrixSVDInverse(MATRIX *m, MATRIX *m_inverse) ;
  95 float  MatrixNSConditionNumber(MATRIX *m);
  96 float  MatrixConditionNumber(MATRIX *m) ;
  97 float  MatrixSVDEigenValues(MATRIX *m, float *evalues) ;
  98 MATRIX *MatrixFactorSqrSVD(MATRIX *M, int Invert, MATRIX *D);
  99 
 100 /* statistical stuff */
 101 MATRIX *MatrixCovariance(MATRIX *mInputs, MATRIX *mCov, VECTOR *mMeans) ;
 102 MATRIX *MatrixUpdateCovariance(MATRIX *mInputs, MATRIX *mCov, VECTOR *mMeans);
 103 MATRIX *MatrixUpdateMeans(MATRIX *mInputs, VECTOR *mMeans, VECTOR *mNobs) ;
 104 MATRIX *MatrixFinalMeans(VECTOR *mMeans, VECTOR *mNobs) ;
 105 MATRIX *MatrixFinalCovariance(MATRIX *mInputs, MATRIX *mCov, VECTOR *mNobs);
 106 
 107 /* misc. I/O functions */
 108 int    MatrixAsciiWriteInto(FILE *fp, MATRIX *m) ;
 109 MATRIX *MatrixAsciiReadFrom(FILE *fp, MATRIX *m) ;
 110 int    MatrixAsciiWrite(char *fname, MATRIX *m) ;
 111 MATRIX *MatrixAsciiRead(char *fname, MATRIX *m) ;
 112 
 113 #define VectorAlloc(n, type)       MatrixAlloc(n, 1, type)
 114 #define RVectorAlloc(n, type)      MatrixAlloc(1, n, type)
 115 #define VectorFree(pm)             MatrixFree(pm)
 116 #define VectorAdd(v1, v2, v3)      MatrixAdd(v1, v2, v3)
 117 #define VectorSubtract(v1,v2,v3)   MatrixSubtract(v1, v2, v3)
 118 #define VectorScalarMul(v1,val,v2) MatrixScalarMul(v1, val, v2)
 119 #define VectorCopy(v1, v2)         MatrixCopy(v1, v2)
 120 #define VectorClear(v)             MatrixClear(v)
 121 #define VectorTranspose(vsrc,vdst) MatrixTranspose(vsrc, vdst)
 122 #define VectorAsciiWriteInto       MatrixAsciiWriteInto
 123 #define VectorAsciiReadFrom        MatrixAsciiReadFrom
 124 
 125 #define VECTOR_ELT(v,i)            ((v)->rptr[i][1])
 126 #define RVECTOR_ELT(v,i)            ((v)->rptr[1][i])
 127 #define VECTOR3_LOAD(v,x,y,z)    (VECTOR_ELT(v,1)=x, VECTOR_ELT(v,2)=y, \
$^
                                  VECTOR_ELT(v,3)=z) ;
 128 #define VECTOR_LOAD   VECTOR3_LOAD
 129 #define V3_LOAD       VECTOR3_LOAD
 130 
 131 double Vector3Angle(VECTOR *v1, VECTOR *v2) ;
 132 float  VectorLen(VECTOR *v) ;
 133 float  VectorAngle(VECTOR *v1, VECTOR *v2) ;
 134 float  VectorDot(VECTOR *v1, VECTOR *v2) ;
 135 float  VectorNormalizedDot(VECTOR *v1, VECTOR *v2) ;
 136 float  VectorDistance(VECTOR *v1, VECTOR *v2) ;
 137 VECTOR *MatrixColumn(MATRIX *m, VECTOR *v, int col) ;
 138 MATRIX *VectorOuterProduct(VECTOR *v1, VECTOR *v2, MATRIX *m) ;
 139 VECTOR *VectorCrossProduct(VECTOR *v1, VECTOR *v2, VECTOR *vdst) ;
 140 float  VectorTripleProduct(VECTOR *v1, VECTOR *v2, VECTOR *v3) ;
 141 VECTOR *VectorNormalize(VECTOR *vin, VECTOR *vout) ;
 142 MATRIX *MatrixNormalizeCol(MATRIX *m, MATRIX *mcnorm);
 143 
 144 /* these are macro versions that work on 3-d vectors */
 145 #define RV3_X(v)     (RVECTOR_ELT(v,1))
 146 #define RV3_Y(v)     (RVECTOR_ELT(v,2))
 147 #define RV3_Z(v)     (RVECTOR_ELT(v,3))
 148 #define V3_X(v)      (VECTOR_ELT(v,1))
 149 #define V3_Y(v)      (VECTOR_ELT(v,2))
 150 #define V3_Z(v)      (VECTOR_ELT(v,3))
 151 #define V3_LEN(v)    (sqrt(V3_X(v)*V3_X(v)+V3_Y(v)*V3_Y(v)+V3_Z(v)*V3_Z(v)))
 152 #define V3_LEN_IS_ZERO(v)    (DZERO(V3_LEN_SQ(v)))
 153 #define V3_LEN_SQ(v) (V3_X(v)*V3_X(v)+V3_Y(v)*V3_Y(v)+V3_Z(v)*V3_Z(v))
 154 #define V3_CROSS_PRODUCT(va,vb,vc) \
$^
                 V3_X(vc) = V3_Y(va)*V3_Z(vb)- V3_Z(va)*V3_Y(vb),  \
 155                  V3_Y(vc) = V3_Z(va)*V3_X(vb)- V3_X(va)*V3_Z(vb),  \
 156                  V3_Z(vc) = V3_X(va)*V3_Y(vb)- V3_Y(va)*V3_X(vb) ;
 157 #define V3_TRIPLE_PRODUCT(va,vb,vc) \
$^
                (V3_X(vc) * (V3_Y(va)*V3_Z(vb)- V3_Z(va)*V3_Y(vb)) + \
 158                  V3_Y(vc) * (V3_Z(va)*V3_X(vb)- V3_X(va)*V3_Z(vb)) + \
 159                  V3_Z(vc) * (V3_X(va)*V3_Y(vb)- V3_Y(va)*V3_X(vb)))
 160 #define V3_ADD(va,vb,vc) \
$^
                 V3_X(vc) = V3_X(va)+V3_X(vb), \
 161                  V3_Y(vc) = V3_Y(va)+V3_Y(vb), \
 162                  V3_Z(vc) = V3_Z(va)+V3_Z(vb) ;
 163 #define V3_SUBTRACT(va,vb,vc) \
$^
                 V3_X(vc) = V3_X(va)-V3_X(vb), \
 164                  V3_Y(vc) = V3_Y(va)-V3_Y(vb), \
 165                  V3_Z(vc) = V3_Z(va)-V3_Z(vb) ;
 166 #define V3_DOT(va,vb) (V3_X(va)*V3_X(vb)+V3_Y(va)*V3_Y(vb)+V3_Z(va)*V3_Z(vb))
 167 #define V3_SCALAR_MUL(va,s,vb)  (V3_X(vb)=V3_X(va)*s,\
$^
                                 V3_Y(vb)=V3_Y(va)*s,\
 168                                  V3_Z(vb)=V3_Z(va)*s)
 169 
 170 #define V3_NORMALIZE(va,vb)  { float len = (V3_LEN(va)) ; \
$^
                                  if (FZERO(len)) len = 1.0f ; \
 171                                   else len = 1.0f / len ; \
 172                                   V3_SCALAR_MUL(va,len,vb) ; }
 173 #define V3_CLEAR(v)          (V3_X(v) = 0, V3_Y(v) = 0, V3_Z(v) = 0)
 174 
 175 #define X_ROTATION   0
 176 #define Y_ROTATION   1
 177 #define Z_ROTATION   2
 178 
 179 MATRIX *MatrixAllocRotation(int n, float angle, int which) ;
 180 MATRIX *MatrixReallocRotation(int n, float angle, int which, MATRIX *m) ;
 181 MATRIX *MatrixAllocTranslation(int n, double *trans) ;
 182 #define MatrixClone(mat)   MatrixCopy(mat, NULL)
 183 #define VectorClone        MatrixClone
 184 
 185 float MatrixTrace(MATRIX *M);
 186 MATRIX *MatrixVertCat(MATRIX *m1, MATRIX *m2, MATRIX *mcat);
 187 MATRIX *MatrixHorCat(MATRIX *m1, MATRIX *m2, MATRIX *mcat);
 188 
 189 MATRIX *MatrixConstVal(float val, int rows, int cols, MATRIX *X);
 190 MATRIX *MatrixZero(int rows, int cols, MATRIX *X);
 191 MATRIX *MatrixSum(MATRIX *m, int dim, MATRIX *msum);
 192 MATRIX *MatrixDRand48(int rows, int cols, MATRIX *m);
 193 MATRIX *MatrixSimilarityTransform(MATRIX *m_src, MATRIX *m_mul, MATRIX *m_dst);
 194 
 195 double VectorSum(MATRIX *v);
 196 double VectorMean(MATRIX *v);
 197 double VectorVar(MATRIX *v, double *pMean);
 198 double VectorStdDev(MATRIX *v, double *pMean);
 199 double VectorRange(MATRIX *v, double *pVmin, double *pVmax);
 200 
 201 MATRIX *GaussianMatrix(int len, float std, int norm, MATRIX *G);
 202 MATRIX *GaussianVector(int len, float mean, float std, int norm, MATRIX *g);
 203 MATRIX *MatrixReorderRows(MATRIX *X, int *NewRowOrder, MATRIX *XRO);
 204 int MatrixRandPermRows(MATRIX *X);
 205 
 206 #endif