39 #include "EST_FMatrix.h"
40 #include "EST_system.h"
52 for(
int i=0; i<x.
n(); ++i)
55 return polynomial_fit(x,y,co_effs,weights,order);
63 cerr <<
"polynomial_fit : order must be >= 1" << endl;
68 cerr <<
"polynomial_fit : x and y must have same dimension" << endl;
72 if(weights.
n() != y.
n()){
73 cerr <<
"polynomial_fit : weights must have same dimension as x and y" << endl;
78 cerr <<
"polynomial_fit : x and y must have at least order+1 elements"
91 for(
int row=0;row<y.
n();row++)
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];
110 if(!inverse(At_A,At_A_inv,singularity))
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;
119 cerr <<
"singularity at point : " << singularity;
120 cerr <<
" = " << x[singularity] <<
"," << y[singularity];
121 cerr <<
" )" << endl;
127 co_effs = At_A_inv * At_y1;
168 cerr <<
"diagonalise: non-square matrix ";
196 for (i = I = 0; i < n; ++i, ++I)
200 for (j = J = 0; j < n; ++j, ++J)
219 s(0, i) = a.a(row, i);
230 s(i, 0) = a.a(i, col);
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);
285 static void row_swap(
int from,
int to,
EST_FMatrix &a)
304 return inverse(a,inv,singularity);
323 int r=0,this_row,all_zeros;
399 singularity = ((this_row > j) ? this_row : j);
411 return pseudo_inverse(a,inv,singularity);
422 return inverse(a,inv,singularity);
429 transpose(a,a_trans);
430 multiply(a_trans,a,atrans_a);
431 if (!inverse(atrans_a,atrans_a_inverse,singularity))
433 multiply(atrans_a_inverse,a_trans,inv);
446 cerr <<
"Tried to take determinant of non-square matrix\n";
460 for (i = 0; i < n; ++i)
462 p = (float)(i + j + 2);
464 A[i] = pow((
float)-1.0, p) * determinant(sub(a, i, j));
470 for (i = 0; i < n; ++i)
495 cerr <<
"Can't make non-square identity matrix !" << endl;
511 cerr <<
"Can't add vectors of differing lengths !" << endl;
516 for(
int i=0; i<a_len; i++ )
529 cerr <<
"Can't subtract vectors of differing lengths !" << endl;
534 for(
int i=0; i<a_len; i++ )
546 cerr <<
"Can't extract diagonal of non-square matrix !" << endl;
557 float polynomial_value(
const EST_FVector &coeffs,
const float x)
561 for(
int i=0;i<coeffs.
length();i++)
577 cerr <<
"Can't symmetrize non-square matrix !" << endl;
602 make_random_matrix(
EST_FMatrix &M,
const float scale)
606 for(
int row=0;row<M.
num_rows();row++)
609 r = scale * ((double)rand()/(double)RAND_MAX);
615 make_random_vector(
EST_FVector &V,
const float scale)
619 for(
int i=0;i<V.
length();i++)
621 r = scale * ((double)rand()/(double)RAND_MAX);
627 make_random_symmetric_matrix(
EST_FMatrix &M,
const float scale)
631 cerr <<
"Can't make non-square symmetric matrix !" << endl;
637 for(
int row=0;row<M.
num_rows();row++)
638 for(
int col=0;col<=row;col++)
640 r = scale * ((double)rand()/(double)RAND_MAX);
647 make_random_diagonal_matrix(
EST_FMatrix &M,
const float scale)
651 cerr <<
"Can't make non-square symmetric matrix !" << endl;
656 for(
int row=0;row<M.
num_rows();row++)
657 M.
a_no_check(row,row) = scale * ((double)rand()/(double)RAND_MAX);
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;
672 for(
int row=0;row<T.
num_rows();row++)
674 T.
a_no_check(row,col) = pow(t[row],(
float)col);
701 for(
int i=0;i<v1.
length();i++)
702 for(
int j=0;j<v2.
length();j++)
710 #if !defined(SYSTEM_IS_WIN32)
714 gettimeofday(&tp,&tzp);
715 seed = getpid() * (tp.tv_usec&0x7fff);
716 cerr <<
"seed: " << seed << endl;
729 gettimeofday(&tp,&tzp);
730 seed = (getpid()&0x7f) * (tp.tv_usec&0xff);
731 cerr <<
"seed48: " << seed << endl;