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