Edinburgh Speech Tools 2.4-release
vec_mat_aux.cc
1/*************************************************************************/
2/* */
3/* Centre for Speech Technology Research */
4/* University of Edinburgh, UK */
5/* Copyright (c) 1995,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/* Author : Simon King */
34/* Date : April 1995 */
35/*-----------------------------------------------------------------------*/
36/* EST_FMatrix Class auxiliary functions */
37/* */
38/*=======================================================================*/
39#include "EST_FMatrix.h"
40#include "EST_system.h"
41#include <cstdlib>
42#include <climits>
43#include "EST_unix.h"
44#include "EST_math.h"
45#include <ctime>
46
47bool polynomial_fit(EST_FVector &x, EST_FVector &y,
48 EST_FVector &co_effs, int order)
49{
50 EST_FVector weights;
51 weights.resize(x.n());
52 for(int i=0; i<x.n(); ++i)
53 weights[i] = 1.0;
54
55 return polynomial_fit(x,y,co_effs,weights,order);
56}
57
58bool polynomial_fit(EST_FVector &x, EST_FVector &y, EST_FVector &co_effs,
59 EST_FVector &weights, int order)
60{
61
62 if(order <= 0){
63 cerr << "polynomial_fit : order must be >= 1" << endl;
64 return false;
65 }
66
67 if(x.n() != y.n()){
68 cerr << "polynomial_fit : x and y must have same dimension" << endl;
69 return false;
70 }
71
72 if(weights.n() != y.n()){
73 cerr << "polynomial_fit : weights must have same dimension as x and y" << endl;
74 return false;
75 }
76
77 if(x.n() <= order){
78 cerr << "polynomial_fit : x and y must have at least order+1 elements"
79 << endl;
80 return false;
81 }
82
83
84 // matrix of basis function values
86 A.resize(x.n(),order+1);
87
88 EST_FVector y1;
89 y1.resize(y.n());
90
91 for(int row=0;row<y.n();row++)
92 {
93 y1[row] = y[row] * weights[row];
94 for(int col=0;col<=order;col++){
95 A(row,col) = pow(x[row],(float)col) * weights[row];
96
97 }
98 }
99
100 // could call pseudo_inverse, but save a bit by doing
101 // it here since we need transpose(A) anyway
102
103 EST_FMatrix At, At_A, At_A_inv;
104 int singularity=-2;
105
106 transpose(A,At);
107 multiply(At,A,At_A);
108
109 // error check
110 if(!inverse(At_A,At_A_inv,singularity))
111 {
112 cerr << "polynomial_fit : inverse failed (";
113 if(singularity == -2)
114 cerr << "unspecified reason)" << endl;
115 else if(singularity == -1)
116 cerr << "non-square !!)" << endl;
117 else
118 {
119 cerr << "singularity at point : " << singularity;
120 cerr << " = " << x[singularity] << "," << y[singularity];
121 cerr << " )" << endl;
122 }
123 return false;
124 }
125
126 EST_FVector At_y1 = At * y1;
127 co_effs = At_A_inv * At_y1;
128 return true;
129
130}
131
132float matrix_max(const EST_FMatrix &a)
133{
134 int i, j;
135 float v = INT_MIN;
136
137 for (i = 0; i < a.num_rows(); ++i)
138 for (j = 0; j < a.num_columns(); ++j)
139 if (a.a_no_check(i, j) > v)
140 v = a.a_no_check(i, j);
141
142 return v;
143}
144
145int square(const EST_FMatrix &a)
146{
147 return a.num_rows() == a.num_columns();
148}
149// add all elements in matrix.
150float sum(const EST_FMatrix &a)
151{
152 int i, j;
153 float t = 0.0;
154 for (i = 0; i < a.num_rows(); ++i)
155 for (j = 0; j < a.num_columns(); ++j)
156 t += a.a_no_check(i, j);
157 return t;
158}
159
160// set all elements not on the diagonal to zero.
161EST_FMatrix diagonalise(const EST_FMatrix &a)
162{
163 int i;
164 EST_FMatrix b(a, 0); // initialise and fill b with zeros
165
166 if (a.num_rows() != a.num_columns())
167 {
168 cerr << "diagonalise: non-square matrix ";
169 return b;
170 }
171
172 for (i = 0; i < a.num_rows(); ++i)
173 b(i, i) = a.a_no_check(i, i);
174
175 return b;
176}
177
178// set all elements not on the diagonal to zero.
179void inplace_diagonalise(EST_FMatrix &a)
180{
181 // NB - will work on non-square matrices without warning
182 int i,j;
183
184 for (i = 0; i < a.num_rows(); ++i)
185 for (j = 0; j < a.num_columns(); ++j)
186 if(i != j)
187 a.a_no_check(i, j) = 0;
188}
189
190EST_FMatrix sub(const EST_FMatrix &a, int row, int col)
191{
192 int i, j, I, J;
193 int n = a.num_rows() - 1;
194 EST_FMatrix s(n, n);
195
196 for (i = I = 0; i < n; ++i, ++I)
197 {
198 if (I == row)
199 ++I;
200 for (j = J = 0; j < n; ++j, ++J)
201 {
202 if (J == col)
203 ++J;
204 s(i, j) = a.a(I, J);
205 }
206 }
207
208 // cout << "sub: row " << row << " col " << col << "\n" << s;
209
210 return s;
211}
212
213EST_FMatrix row(const EST_FMatrix &a, int row)
214{
215 EST_FMatrix s(1, a.num_columns());
216 int i;
217
218 for (i = 0; i < a.num_columns(); ++i)
219 s(0, i) = a.a(row, i);
220
221 return s;
222}
223
224EST_FMatrix column(const EST_FMatrix &a, int col)
225{
226 EST_FMatrix s(a.num_rows(), 1);
227 int i;
228
229 for (i = 0; i < a.num_rows(); ++i)
230 s(i, 0) = a.a(i, col);
231
232 return s;
233}
234
235EST_FMatrix triangulate(const EST_FMatrix &a)
236{
237 EST_FMatrix b(a, 0);
238 int i, j;
239
240 for (i = 0; i < a.num_rows(); ++i)
241 for (j = i; j < a.num_rows(); ++j)
242 b(j, i) = a.a(j, i);
243
244 return b;
245}
246
247void transpose(const EST_FMatrix &a,EST_FMatrix &b)
248{
249 int i, j;
250 b.resize(a.num_columns(), a.num_rows());
251
252 for (i = 0; i < b.num_rows(); ++i)
253 for (j = 0; j < b.num_columns(); ++j)
254 b.a_no_check(i, j) = a.a_no_check(j, i);
255}
256
257EST_FMatrix backwards(EST_FMatrix &a)
258{
259 int i, j, n;
260 n = a.num_columns();
261 EST_FMatrix t(n, n);
262
263 for (i = 0; i < n; ++i)
264 for (j = 0; j < n; ++j)
265 t(n - i - 1, n - j - 1) = a.a(i, j);
266
267 return t;
268}
269
270
271// changed name from abs as it causes on at least on POSIX machine
272// where int abs(int) is a macro
273EST_FMatrix fmatrix_abs(const EST_FMatrix &a)
274{
275 int i, j;
276 EST_FMatrix b(a, 0);
277
278 for (i = 0; i < a.num_rows(); ++i)
279 for (j = 0; j < a.num_columns(); ++j)
280 b.a_no_check(i, j) = fabs(a.a_no_check(i, j));
281
282 return b;
283}
284
285static void row_swap(int from, int to, EST_FMatrix &a)
286{
287 int i;
288 float f;
289
290 if (from == to)
291 return;
292
293 for (i=0; i < a.num_columns(); i++)
294 {
295 f = a.a_no_check(to,i);
296 a.a_no_check(to,i) = a.a_no_check(from,i);
297 a.a_no_check(from,i) = f;
298 }
299}
300
301int inverse(const EST_FMatrix &a,EST_FMatrix &inv)
302{
303 int singularity=0;
304 return inverse(a,inv,singularity);
305}
306
307int inverse(const EST_FMatrix &a,EST_FMatrix &inv,int &singularity)
308{
309
310 // Used to use a function written by Richard Tobin (in C) but
311 // we no longer need C functionality any more. This algorithm
312 // follows that in "Introduction to Algorithms", Cormen, Leiserson
313 // and Rivest (p759) and the term Gauss-Jordon is used for some part,
314 // As well as looking back at Richard's.
315 // This also keeps a record of which rows are which from the original
316 // so that it can return which column actually has the singularity
317 // in it if it fails to find an inverse.
318 int i, j, k;
319 int n = a.num_rows();
320 EST_FMatrix b = a; // going to destructively manipulate b to get inv
321 EST_FMatrix pos; // the original position
322 float biggest,s;
323 int r=0,this_row,all_zeros;
324
325 singularity = -1;
326 if (a.num_rows() != a.num_columns())
327 return FALSE;
328
329 // Make the inverse the identity matrix.
330 inv.resize(n,n);
331 pos.resize(n,1);
332 for (i=0; i<n; i++)
333 for (j=0; j<n; j++)
334 inv.a_no_check(i,j) = 0.0;
335 for (i=0; i<n; i++)
336 {
337 inv.a_no_check(i,i) = 1.0;
338 pos.a_no_check(i,0) = (float)i;
339 }
340
341 // Manipulate b to make it into the identity matrix, while performing
342 // the same manipulation on inv. Once b becomes the identity inv will
343 // be the inverse (unless theres a singularity)
344
345 for (i=0; i<n; i++)
346 {
347 // Find the absolute largest val in this col as the next to
348 // manipulate.
349 biggest = 0.0;
350 r = 0;
351 for (j=i; j<n; j++)
352 {
353 if (fabs(b.a_no_check(j,i)) > biggest)
354 {
355 r = j;
356 biggest = fabs(b.a_no_check(j,i));
357 }
358 }
359
360 if (biggest == 0.0) // oops found a singularity
361 {
362 singularity = (int)pos.a_no_check(i,0);
363 return FALSE;
364 }
365
366 // Swap current with biggest
367 this_row = (int)pos.a_no_check(i,0); // in case we need this number
368 row_swap(r,i,b);
369 row_swap(r,i,inv);
370 row_swap(r,i,pos);
371
372 // Make b(i,i) = 1
373 s = b(i,i);
374 for (k=0; k<n; k++)
375 {
376 b.a_no_check(i,k) /= s;
377 inv.a_no_check(i,k) /= s;
378 }
379
380 // make rest in col(i) 0
381 for (j=0; j<n; j++)
382 {
383 if (j==i) continue;
384 s = b.a_no_check(j,i);
385 all_zeros = TRUE;
386 for (k=0; k<n; k++)
387 {
388 b.a_no_check(j,k) -= b.a_no_check(i,k) * s;
389 if (b.a_no_check(j,k) != 0)
390 all_zeros = FALSE;
391 inv.a_no_check(j,k) -= inv.a_no_check(i,k) * s;
392 }
393 if (all_zeros)
394 {
395 // printf("singularity between (internal) columns %d %d\n",
396 // this_row,j);
397 // always identify greater row so linear regression
398 // can preserve intercept in column 0
399 singularity = ((this_row > j) ? this_row : j);
400 return FALSE;
401 }
402 }
403 }
404
405 return TRUE;
406}
407
408int pseudo_inverse(const EST_FMatrix &a, EST_FMatrix &inv)
409{
410 int singularity=0;
411 return pseudo_inverse(a,inv,singularity);
412}
413
414int pseudo_inverse(const EST_FMatrix &a, EST_FMatrix &inv,int &singularity)
415{
416 // for non-square matrices
417 // useful for solving linear eqns
418 // (e.g. polynomial fitting)
419
420 // is it square ?
421 if( a.num_rows() == a.num_columns() )
422 return inverse(a,inv,singularity);
423
424 if( a.num_rows() < a.num_columns() )
425 return FALSE;
426
427 EST_FMatrix a_trans,atrans_a,atrans_a_inverse;
428
429 transpose(a,a_trans);
430 multiply(a_trans,a,atrans_a);
431 if (!inverse(atrans_a,atrans_a_inverse,singularity))
432 return FALSE;
433 multiply(atrans_a_inverse,a_trans,inv);
434
435 return TRUE;
436}
437
438
439float determinant(const EST_FMatrix &a)
440{
441 int i, j;
442 int n = a.num_rows();
443 float det;
444 if (!square(a))
445 {
446 cerr << "Tried to take determinant of non-square matrix\n";
447 return 0.0;
448 }
449
450 EST_FVector A(n);
451
452 if (n == 2) // special case of 2x2 determinant
453 return (a.a_no_check(0,0) * a.a_no_check(1,1)) -
454 (a.a_no_check(0,1) * a.a_no_check(1,0));
455
456 float p;
457
458 // create cofactor matrix
459 j = 1;
460 for (i = 0; i < n; ++i)
461 {
462 p = (float)(i + j + 2); // because i & j should start at 1
463 // cout << "power " <<p << endl;
464 A[i] = pow((float)-1.0, p) * determinant(sub(a, i, j));
465 }
466 // cout << "cofactor " << A;
467
468 // sum confactor and original matrix
469 det = 0.0;
470 for (i = 0; i < n; ++i)
471 det += a.a_no_check(i, j) * A[i];
472
473 return det;
474}
475
476void eye(EST_FMatrix &a, const int n)
477{
478 int i,j;
479 a.resize(n,n);
480 for (i=0; i<n; i++)
481 {
482 for (j=0; j<n; j++)
483 a.a_no_check(i,j) = 0.0;
484
485 a.a_no_check(i,i) = 1.0;
486 }
487}
488
489void eye(EST_FMatrix &a)
490{
491 int i,n;
492 n=a.num_rows();
493 if(n != a.num_columns())
494 {
495 cerr << "Can't make non-square identity matrix !" << endl;
496 return;
497 }
498
499 a.fill(0.0);
500 for (i=0; i<n; i++)
501 a.a_no_check(i,i) = 1.0;
502}
503
504EST_FVector add(const EST_FVector &a,const EST_FVector &b)
505{
506 // a + b
507 int a_len = a.length();
508 EST_FVector ans( a_len );
509
510 if(a_len != b.length()){
511 cerr << "Can't add vectors of differing lengths !" << endl;
512 ans.resize(0);
513 return ans;
514 };
515
516 for( int i=0; i<a_len; i++ )
517 ans.a_no_check(i) = a.a_no_check(i) + b.a_no_check(i);
518
519 return ans;
520}
521
522EST_FVector subtract(const EST_FVector &a,const EST_FVector &b)
523{
524 // a - b
525 int a_len = a.length();
526 EST_FVector ans( a_len );
527
528 if(a_len != b.length()){
529 cerr << "Can't subtract vectors of differing lengths !" << endl;
530 ans.resize(0);
531 return ans;
532 };
533
534 for( int i=0; i<a_len; i++ )
535 ans.a_no_check(i) = a.a_no_check(i) - b.a_no_check(i);
536
537 return ans;
538}
539
540EST_FVector diagonal(const EST_FMatrix &a)
541{
542
543 EST_FVector ans;
544 if(a.num_rows() != a.num_columns())
545 {
546 cerr << "Can't extract diagonal of non-square matrix !" << endl;
547 return ans;
548 }
549 int i;
550 ans.resize(a.num_rows());
551 for(i=0;i<a.num_rows();i++)
552 ans.a_no_check(i) = a.a_no_check(i,i);
553
554 return ans;
555}
556
557float polynomial_value(const EST_FVector &coeffs, const float x)
558{
559 float y=0;
560
561 for(int i=0;i<coeffs.length();i++)
562 y += coeffs.a_no_check(i) * pow(x,(float)i);
563
564 return y;
565}
566
567void symmetrize(EST_FMatrix &a)
568{
569 // uses include enforcing symmetry
570 // of covariance matrices (to compensate
571 // for rounding errors)
572
573 float f;
574
575 if(a.num_columns() != a.num_rows())
576 {
577 cerr << "Can't symmetrize non-square matrix !" << endl;
578 return;
579 }
580
581 // no need to look at entries on the diagonal !
582 for(int i=0;i<a.num_rows();i++)
583 for(int j=i+1;j<a.num_columns();j++)
584 {
585 f = 0.5 * (a.a_no_check(i,j) + a.a_no_check(j,i));
586 a.a_no_check(i,j) = a.a_no_check(j,i) = f;
587 }
588}
589
590void
591stack_matrix(const EST_FMatrix &M, EST_FVector &v)
592{
593 v.resize(M.num_rows() * M.num_columns());
594 int i,j,k=0;
595 for(i=0;i<M.num_rows();i++)
596 for(j=0;j<M.num_columns();j++)
597 v.a_no_check(k++) = M(i,j);
598}
599
600
601void
602make_random_matrix(EST_FMatrix &M, const float scale)
603{
604
605 float r;
606 for(int row=0;row<M.num_rows();row++)
607 for(int col=0;col<M.num_columns();col++)
608 {
609 r = scale * ((double)rand()/(double)RAND_MAX);
610 M.a_no_check(row,col) = r;
611 }
612}
613
614void
615make_random_vector(EST_FVector &V, const float scale)
616{
617
618 float r;
619 for(int i=0;i<V.length();i++)
620 {
621 r = scale * ((double)rand()/(double)RAND_MAX);
622 V.a_no_check(i) = r;
623 }
624}
625
626void
627make_random_symmetric_matrix(EST_FMatrix &M, const float scale)
628{
629 if(M.num_rows() != M.num_columns())
630 {
631 cerr << "Can't make non-square symmetric matrix !" << endl;
632 return;
633 }
634
635 float r;
636
637 for(int row=0;row<M.num_rows();row++)
638 for(int col=0;col<=row;col++)
639 {
640 r = scale * ((double)rand()/(double)RAND_MAX);
641 M.a_no_check(row,col) = r;
642 M.a_no_check(col,row) = r;
643 }
644}
645
646void
647make_random_diagonal_matrix(EST_FMatrix &M, const float scale)
648{
649 if(M.num_rows() != M.num_columns())
650 {
651 cerr << "Can't make non-square symmetric matrix !" << endl;
652 return;
653 }
654
655 M.fill(0.0);
656 for(int row=0;row<M.num_rows();row++)
657 M.a_no_check(row,row) = scale * ((double)rand()/(double)RAND_MAX);
658
659
660}
661
662void
663make_poly_basis_function(EST_FMatrix &T, EST_FVector t)
664{
665 if(t.length() != T.num_rows())
666 {
667 cerr << "Can't make polynomial basis function : dimension mismatch !" << endl;
668 cerr << "t.length()=" << t.length();
669 cerr << " T.num_rows()=" << T.num_rows() << endl;
670 return;
671 }
672 for(int row=0;row<T.num_rows();row++)
673 for(int col=0;col<T.num_columns();col++)
674 T.a_no_check(row,col) = pow(t[row],(float)col);
675
676}
677
678int
679floor_matrix(EST_FMatrix &M, const float floor)
680{
681 int i,j,k=0;
682 for(i=0;i<M.num_rows();i++)
683 for(j=0;j<M.num_columns();j++)
684 if(M.a_no_check(i,j) < floor)
685 {
686 M.a_no_check(i,j) = floor;
687 k++;
688 }
689 return k;
690}
691
693cov_prod(const EST_FVector &v1,const EST_FVector &v2)
694{
695 // form matrix of vector product, e.g. for covariance
696 // treat first arg as a column vector and second as a row vector
697
698 EST_FMatrix m;
699 m.resize(v1.length(),v2.length());
700
701 for(int i=0;i<v1.length();i++)
702 for(int j=0;j<v2.length();j++)
703 m.a_no_check(i,j) = v1.a_no_check(i) * v2.a_no_check(j);
704
705 return m;
706}
707
708void est_seed()
709{
710#if !defined(SYSTEM_IS_WIN32)
711 unsigned int seed;
712 struct timeval tp;
713 struct timezone tzp;
714 gettimeofday(&tp,&tzp);
715 seed = getpid() * (tp.tv_usec&0x7fff);
716 cerr << "seed: " << seed << endl;
717 srand(seed);
718#else
719 srand(time(NULL));
720#endif
721}
722
723#if 0
724void est_seed48()
725{
726 unsigned short seed;
727 struct timeval tp;
728 struct timezone tzp;
729 gettimeofday(&tp,&tzp);
730 seed = (getpid()&0x7f) * (tp.tv_usec&0xff);
731 cerr << "seed48: " << seed << endl;
732 seed48(&seed);
733}
734#endif
int num_columns() const
return number of columns
Definition: EST_TMatrix.h:181
void fill(const T &v)
fill matrix with value v
Definition: EST_TMatrix.cc:314
INLINE const T & a_no_check(int row, int col) const
const access with no bounds check, care recommend
Definition: EST_TMatrix.h:184
int num_rows() const
return number of rows
Definition: EST_TMatrix.h:179
void resize(int rows, int cols, int set=1)
resize matrix
void resize(int n, int set=1)
resize vector
INLINE int length() const
number of items in vector.
Definition: EST_TVector.h:252
INLINE int n() const
number of items in vector.
Definition: EST_TVector.h:254
INLINE const T & a_no_check(int n) const
read-only const access operator: without bounds checking
Definition: EST_TVector.h:257