1 /*************************************************************************/
2 /*                                                                       */
3 /*                Centre for Speech Technology Research                  */
4 /*                     University of Edinburgh, UK                       */
5 /*                         Copyright (c) 1996                            */
6 /*                        All Rights Reserved.                           */
7 /*                                                                       */
8 /*  Permission is hereby granted, free of charge, to use and distribute  */
9 /*  this software and its documentation without restriction, including   */
10 /*  without limitation the rights to use, copy, modify, merge, publish,  */
11 /*  distribute, sublicense, and/or sell copies of this work, and to      */
12 /*  permit persons to whom this work is furnished to do so, subject to   */
13 /*  the following conditions:                                            */
14 /*   1. The code must retain the above copyright notice, this list of    */
15 /*      conditions and the following disclaimer.                         */
16 /*   2. Any modifications must be clearly marked as such.                */
17 /*   3. Original authors' names are not deleted.                         */
18 /*   4. The authors' names are not used to endorse or promote products   */
19 /*      derived from this software without specific prior written        */
20 /*      permission.                                                      */
21 /*                                                                       */
22 /*  THE UNIVERSITY OF EDINBURGH AND THE CONTRIBUTORS TO THIS WORK        */
23 /*  DISCLAIM ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING      */
24 /*  ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, IN NO EVENT   */
25 /*  SHALL THE UNIVERSITY OF EDINBURGH NOR THE CONTRIBUTORS BE LIABLE     */
26 /*  FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES    */
27 /*  WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN   */
28 /*  AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,          */
29 /*  ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF       */
30 /*  THIS SOFTWARE.                                                       */
31 /*                                                                       */
32 /*************************************************************************/
33 /*                                                                       */
34 /*                     Author :  Simon King                              */
35 /*                     Date   :  February 1999                           */
36 /* --------------------------------------------------------------------- */
37 /*          Double matrix class - copied from FMatrix !                  */
38 /*                                                                       */
39 /*************************************************************************/
40 
41 #ifndef __DMatrix_H__
42 #define __DMatrix_H__
43 
44 #include "EST_TSimpleMatrix.h"
45 #include "EST_TSimpleVector.h"
46 #include "EST_FMatrix.h"
47 
48 
49 class EST_DVector;
50 
51 /** A matrix class for double precision floating point numbers.
52 EST_DMatrix x should be used instead of double **x wherever
53 possible.*/
54 class EST_DMatrix : public EST_TSimpleMatrix<double> {
55 private:
56 public:
57     /// size constructor
EST_DMatrix(int m,int n)58     EST_DMatrix(int m, int n):EST_TSimpleMatrix<double>(m, n)  {}
59     /// copy constructor
EST_DMatrix(const EST_DMatrix & a)60     EST_DMatrix(const EST_DMatrix &a):EST_TSimpleMatrix<double>(a)  {}
61 
62     static EST_String default_file_type;
63     /// CHECK  - what does this do???
64     EST_DMatrix(const EST_DMatrix &a, int b);
65     /// default constructor
EST_DMatrix()66     EST_DMatrix():EST_TSimpleMatrix<double>()  {}
67 
68     /// Save in file (ascii or binary)
69     EST_write_status save(const EST_String &filename,
70 			  const EST_String &type =
71 			        EST_DMatrix::default_file_type);
72     /// Load from file (ascii or binary as defined in file)
73     EST_read_status load(const EST_String &filename);
74     /// Save in file in est format
75     EST_write_status est_save(const EST_String &filename,
76 			      const EST_String &type);
77     /// Load from file in est format (binary/ascii defined in file itself)
78     EST_read_status est_load(const EST_String &filename);
79 
80     /// Copy 2-d array {\tt x} of size {\tt rows x cols} into matrix.
81     void copyin(double **x, int rows, int cols);
82 
83     /// Add elements of 2 same sized matrices.
84     EST_DMatrix &operator+=(const EST_DMatrix &a);
85 
86     /// Subtract elements of 2 same sized matrices.
87     EST_DMatrix &operator-=(const EST_DMatrix &a);
88 
89     /// elementwise multiply by scalar
90     EST_DMatrix &operator*=(const double f);
91 
92     /// elementwise divide by scalar
93     EST_DMatrix &operator/=(const double f);
94 
95     /// Multiply all elements of matrix by {\tt x}.
96     friend EST_DMatrix operator*(const EST_DMatrix &a, const double x);
97 
98     /// Multiply matrix by vector.
99     friend EST_DVector operator*(const EST_DMatrix &a, const EST_DVector &v);
100 
101     /// Multiply vector by matrix
102     friend EST_DVector operator*(const EST_DVector &v,const EST_DMatrix &a);
103 
104     /// Multiply matrix by matrix.
105     friend EST_DMatrix operator*(const EST_DMatrix &a, const EST_DMatrix &b);
106 };
107 
108 
109 /** A vector class for double precision floating point
110     numbers. {\tt EST_DVector x} should be used instead of
111     {\tt float *x} wherever possible.
112 */
113 class EST_DVector: public EST_TSimpleVector<double> {
114 public:
115     /// Size constructor.
EST_DVector(int n)116     EST_DVector(int n): EST_TSimpleVector<double>(n) {}
117     /// Copy constructor.
EST_DVector(const EST_DVector & a)118     EST_DVector(const EST_DVector &a): EST_TSimpleVector<double>(a) {}
119     /// Default constructor.
EST_DVector()120     EST_DVector(): EST_TSimpleVector<double>() {}
121 
122     /// elementwise multiply
123     EST_DVector &operator*=(const EST_DVector &s);
124 
125     /// elementwise add
126     EST_DVector &operator+=(const EST_DVector &s);
127 
128     /// elementwise multiply by scalar
129     EST_DVector &operator*=(const double d);
130 
131     /// elementwise divide by scalar
132     EST_DVector &operator/=(const double d);
133 
134     EST_write_status est_save(const EST_String &filename,
135 			      const EST_String &type);
136 
137     /// save vector to file <tt> filename</tt>.
138     EST_write_status save(const EST_String &filename,
139 			  const EST_String &type);
140 
141     /// load vector from file <tt> filename</tt>.
142     EST_read_status load(const EST_String &filename);
143     /// Load from file in est format (binary/ascii defined in file itself)
144     EST_read_status est_load(const EST_String &filename);
145 };
146 
147 int square(const EST_DMatrix &a);
148 /// inverse
149 int inverse(const EST_DMatrix &a, EST_DMatrix &inv);
150 int inverse(const EST_DMatrix &a, EST_DMatrix &inv, int &singularity);
151 /// pseudo inverse (for non-square matrices)
152 int pseudo_inverse(const EST_DMatrix &a, EST_DMatrix &inv);
153 int pseudo_inverse(const EST_DMatrix &a, EST_DMatrix &inv,int &singularity);
154 
155 /// some useful matrix creators
156 /// make an identity matrix of dimension n
157 void eye(EST_DMatrix &a, const int n);
158 /// make already square matrix into I without resizing
159 void eye(EST_DMatrix &a);
160 
161 /// the user should use est_seed to seed the random number generator
162 void est_seed();
163 void est_seed48();
164 /// all elements are randomised
165 void make_random_vector(EST_DVector &M, const double scale);
166 /// all elements are randomised
167 void make_random_matrix(EST_DMatrix &M, const double scale);
168 /// used for variance
169 void make_random_diagonal_matrix(EST_DMatrix &M, const double scale);
170 /// used for covariance
171 void make_random_symmetric_matrix(EST_DMatrix &M, const double scale);
172 
173 void make_poly_basis_function(EST_DMatrix &T, EST_DVector t);
174 
175 /// elementwise add
176 EST_DVector add(const EST_DVector &a,const EST_DVector &b);
177 /// elementwise subtract
178 EST_DVector subtract(const EST_DVector &a,const EST_DVector &b);
179 
180 /// enforce symmetry
181 void symmetrize(EST_DMatrix &a);
182 /// stack columns on top of each other to make a vector
183 void stack_matrix(const EST_DMatrix &M, EST_DVector &v);
184 /// inplace diagonalise
185 void inplace_diagonalise(EST_DMatrix &a);
186 
187 
188 double determinant(const EST_DMatrix &a);
189 /// not implemented ??
190 int singular(EST_DMatrix &a);
191 /// exchange rows and columns
192 void transpose(const EST_DMatrix &a,EST_DMatrix &b);
193 EST_DMatrix triangulate(const EST_DMatrix &a);
194 
195 /// extract leading diagonal as a matrix
196 EST_DMatrix diagonalise(const EST_DMatrix &a);
197 /// extract leading diagonal as a vector
198 EST_DVector diagonal(const EST_DMatrix &a);
199 /// sum of elements
200 double sum(const EST_DMatrix &a);
201 void multiply(const EST_DMatrix &a, const EST_DMatrix &b, EST_DMatrix &c);
202 int  floor_matrix(EST_DMatrix &M, const double floor);
203 
204 /// matrix product of two vectors (#rows = length of first vector, #cols = length of second vector)
205 EST_DMatrix cov_prod(const EST_DVector &v1,const EST_DVector &v2);
206 
207 EST_DMatrix operator*(const EST_DMatrix &a, const EST_DMatrix &b);
208 EST_DMatrix operator-(const EST_DMatrix &a, const EST_DMatrix &b);
209 EST_DMatrix operator+(const EST_DMatrix &a, const EST_DMatrix &b);
210 
211 EST_DVector operator-(const EST_DVector &a, const EST_DVector &b);
212 EST_DVector operator+(const EST_DVector &a, const EST_DVector &b);
213 
214 EST_DMatrix sub(const EST_DMatrix &a, int row, int col);
215 EST_DMatrix DMatrix_abs(const EST_DMatrix &a);
216 
217 EST_DMatrix row(const EST_DMatrix &a, int row);
218 EST_DMatrix column(const EST_DMatrix &a, int col);
219 
220 
221 /// least squares fit
222 bool
223 polynomial_fit(EST_DVector &x, EST_DVector &y, EST_DVector &co_effs, int order);
224 
225 /// weighted least squares fit
226 bool
227 polynomial_fit(EST_DVector &x, EST_DVector &y, EST_DVector &co_effs,
228 	       EST_DVector &weights, int order);
229 
230 double
231 polynomial_value(const EST_DVector &coeffs, const double x);
232 
233 /// vector dot product
234 double operator*(const EST_DVector &v1, const EST_DVector &v2);
235 
236 
237 #endif
238