This graph shows which files directly or indirectly include this file:
Go to the source code of this file.
Defines | |
#define | integer int |
#define | frand() ((float) rand() / (RAND_MAX+1.0)) |
#define | max(a, b) a > b ? a : b |
#define | min(a, b) a < b ? a : b |
Functions | |
float ** | matrix (int m, int n) |
float ** | submatrix (float **A, int m, int n, int newm, int newn) |
void | free_matrix (float **A) |
void | fill_rand (float **A, int m, int n) |
void | print_matrix (float **A, int m, int n) |
void | transpose_matrix (float **A, int m, int n) |
void | column_orient (float **A, int m, int n, float **B) |
void | row_sum (float **A, int m, int n, float *e) |
void | copy_matrix (float **A, float **B, int m, int n) |
float * | myvector (int n) |
float * | subvector (float *x, int n, int newn) |
void | free_vector (float *x) |
void | print_vector (float *x, int n) |
void | fill_ones (float *x, int n) |
float | sum_vector (float *x, int n) |
void | copy_vector (float *x, float *y, int n) |
void | fillveczeros (float *x, int n) |
float | norm (float *x, int n) |
void | alphaxpy (float alpha, float *x, float *y, int n) |
void | matvec_mult (float **A, float *x, int m, int n, float *b) |
void | matvec_multT (float **A, float *x, int m, int n, float *b) |
void | get_col (float **A, int m, int k, float *ak) |
void | get_row (float **A, int n, int k, float *ak) |
void | put_col (float **A, int m, float *v, int k) |
void | svd (float **A, int m, int n, float **U, float *x, float **V) |
void | tsvd (float **U, float *s, float **V, float *b, int m, int n, float omega, float *x, float *tol) |
void | gcv_tsvd (float *s, float *bhat, int m, int n, float omega, float *tol) |
void | tikhonov (float **U, float *s, float **V, float *b, int m, int n, float omega, float *x, float *alpha) |
float | GCVmin_Tik (float *s, float *bhat, int m, int n, int MaxIts, float omega) |
float | GCVfun_Tik (float alpha, float *s, float *bhat, int m, int n, float omega) |
void | generate_noise (double *E, int m, int n) |
|
Definition at line 4 of file hybr.h. Referenced by fill_rand(). |
|
|
|
|
Definition at line 822 of file utilnum.cpp.
|
|
Definition at line 690 of file utilnum.cpp. 00691 { 00692 int i,j; 00693 for (j = 0; j < n; j++){ 00694 for (i = 0; i < m; i++) 00695 { 00696 B[j][i] = A[i][j]; 00697 } 00698 } 00699 }
|
|
Definition at line 714 of file utilnum.cpp. 00715 { 00716 int i,j; 00717 for(i = 0; i < m; i++) { 00718 for (j = 0; j < n; j++) 00719 B[i][j]=A[i][j]; 00720 00721 } 00722 }
|
|
Definition at line 791 of file utilnum.cpp.
|
|
Definition at line 771 of file utilnum.cpp. 00772 { 00773 int i; 00774 for (i = 0; i < m; i++){ 00775 x[i] = (float) 1.0; 00776 } 00777 }
|
|
Definition at line 651 of file utilnum.cpp. 00652 { 00653 int i,j; 00654 for(i = 0; i < m; i++) { 00655 for (j = 0; j < n; j++) 00656 A[i][j]=frand(); 00657 00658 } 00659 }
|
|
Definition at line 799 of file utilnum.cpp. 00800 { 00801 int i; 00802 for (i=0; i<n; i++) 00803 x[i] = 0.0; 00804 }
|
|
Definition at line 644 of file utilnum.cpp. 00645 { 00646 free(A[0]); 00647 free(A); 00648 }
|
|
Definition at line 755 of file utilnum.cpp. 00756 { 00757 free(x); 00758 }
|
|
Definition at line 200 of file utilnum.cpp. 00206 : s : Vector containing singular values. 00207 * bhat : Vector containing spectral coefficients of right 00208 * hand side vector. That is, bhat = U^T*b. 00209 * 00210 * Output: tol : Scalar truncation tolerance for TSVD. 00211 */ 00212 { 00213 float * bhat_copy = myvector(m); 00214 float * rho = myvector(m-1); 00215 float * G = myvector(m-1); 00216 int i, kk=0; 00217 float ni; 00218 /* 00219 * NEED TO: Make sure singular values are sorted from 00220 * largest to smallest, and do the same for the 00221 * bhat values. In MATLAB this looks like: 00222 * [s, idx] = sort(abs(s)); 00223 * s = flipud(s); 00224 * idx = flipud(idx); 00225 * bhat = abs( bhat(idx) ); 00226 * 00227 * Initialize residual and GCV function values. 00228 */ 00229 00230 copy_vector(bhat, bhat_copy, m); 00231 for (i=0; i<m; i++) 00232 bhat_copy[i] = fabs( bhat[i] ); 00233 00234 rho[m-2] = bhat_copy[m-1]*bhat_copy[m-1]; 00235 00236 /* 00237 * Recursively compute residual and GCV function values. 00238 */ 00239 for (i = m-3; i>-1; i--){ 00240 rho[i] = rho[i+1] + bhat_copy[i+1]*bhat_copy[i+1]; 00241 } 00242 00243 // print_vector(bhat_copy, m-1); 00244 for(i=0; i< m-1; i++){ 00245 ni = (float) m - omega*(i+1); 00246 G[i] = (float) n* rho[i]/(ni*ni); 00247 } 00248 // print_vector(G,m-1); 00249 /* 00250 * Now find minimum of the discrete GCV function. 00251 */ 00252 for (i=0; i < m-2; i++){ 00253 if (G[i+1]<G[i]) 00254 kk = i+1; 00255 } 00256 00257 *tol = s[kk]; 00258 // printf("kk = %d, s[kk] = %f\n", kk, *tol); 00259 00260 free_vector(bhat_copy); 00261 free_vector(rho); 00262 free_vector(G); 00263 }
|
|
Definition at line 483 of file utilnum.cpp. 00488 : alpha: value at which to evaluate the GCV function 00489 * s : Vector containing singular or spectral values of A. 00490 * bhat : Vector containing spectral coefficients of right 00491 * hand side vector. That is, bhat = U^T*b. 00492 * 00493 * Output: G : GCV function for Tikhonov at alpha 00494 */ 00495 { 00496 float t0, G; 00497 float * hold = myvector(m-n); 00498 float * tt= myvector(n); 00499 float * t1= myvector(n); 00500 float * t2= myvector(n); 00501 float * t3= myvector(n); 00502 int i; 00503 00504 for (i=0; i< m-n; i++) 00505 hold[i] = fabs(bhat[n])*fabs(bhat[n]); 00506 t0 = sum_vector(hold,m-n); 00507 00508 for (i = 0; i<n; i++){ 00509 tt[i] = (float)1.0 /(s[i]*s[i] +alpha*alpha); 00510 t1[i] = alpha*alpha * tt[i]; 00511 t2[i] = fabs(bhat[i] * t1[i]) * fabs(bhat[i] * t1[i]); 00512 t3[i] = t1[i] + ((float)1.0 - omega)* (s[i]*s[i]) * tt[i]; 00513 } 00514 G = (float) n*(sum_vector(t2,n) + t0) / ((sum_vector(t3,n) + (float)m - (float)n)* (sum_vector(t3,n) + (float)m - (float)n)); 00515 00516 //printf(" svec = %f, to = %f, Num = %f; Dem = %f \n", sum_vector(t2,n), t0, (float) n*(sum_vector(t2,n) + t0), ((sum_vector(t3,n) + (float)m - (float)n)* (sum_vector(t3,n) + (float)m - (float)n))); 00517 00518 free_vector(hold); 00519 free_vector(tt); 00520 free_vector(t1); 00521 free_vector(t2); 00522 free_vector(t3); 00523 00524 return G; 00525 }
|
|
Definition at line 329 of file utilnum.cpp. 00336 : s : Vector containing singular or spectral values of A. 00337 * bhat : Vector containing spectral coefficients of right 00338 * hand side vector. That is, bhat = U^T*b. 00339 * MaxIts : Maximum number of iterations to find the minimum 00340 * 00341 * Output: xf : minimum of the GCV function for Tikhonov. 00342 */ 00343 { 00344 float ax, bx, seps, c, a, b, v, w, sf, d, e, x, fx, fv, fw, xm, tol, tol1, tol2, r, q, p, si, fu, fval, xf; 00345 int funccount, iter, maxfun, gs; 00346 00347 /*float tt; 00348 00349 tt = GCVfun_Tik(0.01, s, bhat, m,n, omega); 00350 printf("alpha = 0.01; GCV(a) = %f\n", tt); 00351 tt = GCVfun_Tik(0.59, s, bhat, m,n, omega); 00352 printf("alpha = 0.59; GCV(a) = %f\n", tt); 00353 tt = GCVfun_Tik(.999, s, bhat, m,n, omega); 00354 printf("alpha = .999; GCV(a) = %f\n", tt); 00355 tt = GCVfun_Tik(1.2, s, bhat, m,n, omega); 00356 printf("alpha = 1.2; GCV(a) = %f\n", tt); 00357 */ 00358 00359 ax = (float) 0.0; 00360 // bx = (float)1.0; 00361 bx = s[0]; 00362 00363 /* Initialize some values */ 00364 seps = sqrt(FLT_EPSILON); 00365 c = (float)0.5*((float)3.0 - sqrt((float)5.0)); 00366 a = ax; 00367 b = bx; 00368 v = a + c*(b-a); 00369 w = v; 00370 xf = v; 00371 d = (float)0.0; 00372 e = (float)0.0; 00373 x = xf; 00374 00375 fx = GCVfun_Tik(x, s, bhat, m,n, omega); 00376 00377 fv = fx; 00378 fw = fx; 00379 xm = (float)0.5*(a+b); 00380 tol = (float).0001; 00381 tol1 = seps*fabs(xf) + tol/(float)3.0; 00382 tol2 = (float)2.0*tol1; 00383 funccount = 0; 00384 iter = 0; 00385 maxfun = 500; 00386 00387 for (iter = 1; iter< MaxIts; iter++){ 00388 gs = 1; 00389 /* Is a parabolic fit possible*/ 00390 if (fabs(e) > tol1){ 00391 gs = 0; 00392 r = (xf-w)*(fx-fv); 00393 q = (xf-v)*(fx-fw); 00394 p = (xf-v)*q-(xf-w)*r; 00395 q = (float)2.0*(q-r); 00396 if (q > (float)0.0) 00397 p = -p; 00398 q = fabs(q); 00399 r = e; 00400 e = d; 00401 /*Is the parabola acceptable*/ 00402 if ( (fabs(p)<fabs((float)0.5*q*r)) && (p>q*(a-xf)) && (p<q*(b-xf)) ){ 00403 /* Yes, parabolic interpolation step*/ 00404 d = p/q; 00405 x = xf+d; 00406 if (((x-a) < tol2) || ((b-x) < tol2)){ 00407 if(xm-xf <0) 00408 si = -(float)1.0; 00409 else 00410 si = (float)1.0; 00411 00412 d = tol1*si; 00413 } 00414 } 00415 else{ 00416 /* Not acceptable, must do a golden section step*/ 00417 gs = 1; 00418 } 00419 } 00420 if (gs){ 00421 /* A golden-section step is required*/ 00422 if (xf >= xm) 00423 e = a-xf; 00424 else 00425 e = b-xf; 00426 00427 d = c*e; 00428 } 00429 /* The function must not be evaluated too close to xf*/ 00430 if(d <0) 00431 si = -(float)1.0; 00432 else 00433 si = (float)1.0; 00434 00435 x = xf + si * max( fabs(d), tol1 ); 00436 fu = GCVfun_Tik( x, s, bhat, m,n, omega); 00437 funccount = funccount + 1; 00438 if (fu <= fx){ 00439 if (x >= xf) 00440 a = xf; 00441 else 00442 b = xf; 00443 00444 v = w; 00445 fv = fw; 00446 w = xf; 00447 fw = fx; 00448 xf = x; 00449 fx = fu; 00450 } 00451 else{ 00452 if (x < xf) 00453 a = x; 00454 else 00455 b = x; 00456 00457 if ( (fu <= fw) || (w == xf) ){ 00458 v = w; 00459 fv = fw; 00460 w = x; 00461 fw = fu; 00462 } 00463 else if ((fu <= fv) || (v == xf) || (v == w)){ 00464 v = x; 00465 fv = fu; 00466 } 00467 } 00468 xm = (float)0.5*(a+b); 00469 tol1 = seps*fabs(xf) + tol/(float)3.0; 00470 tol2 = (float)2.0*tol1; 00471 if (funccount >= maxfun || fabs(xf-xm) <= (tol2 - (float)0.5*(b-a))){ 00472 fval = fx; 00473 return xf; 00474 } 00475 } 00476 fval = fx; 00477 return xf; 00478 }
|
|
|
|
Definition at line 862 of file utilnum.cpp. 00863 { 00864 int i; 00865 for (i = 0; i < m; i++) 00866 ak[i] = A[i][k]; 00867 00868 }
|
|
Definition at line 873 of file utilnum.cpp. 00874 { 00875 int i; 00876 for (i = 0; i < n; i++) 00877 ak[i] = A[k][i]; 00878 00879 }
|
|
Definition at line 606 of file utilnum.cpp. Referenced by EMAN::Util::assign_groups(), EMAN::Transform::at(), EMAN::Transform::get_matrix3_row(), EMAN::Transform::get_rotation(), EMAN::Icosahedral2Sym::get_sym(), EMAN::Transform::operator[](), EMAN::Transform::printme(), EMAN::AreaProcessor::process_inplace(), EMAN::ZeroConstantProcessor::process_pixel(), EMAN::Transform::set(), and EMAN::Transform::transform(). 00607 { 00608 int i; 00609 float **A=(float **) calloc(m, sizeof(float*)); 00610 if (A == NULL) 00611 { 00612 fprintf(stderr, "Allocation failure in matrix()"); 00613 exit (1); 00614 } 00615 A[0]=(float *) calloc(m * n, sizeof(float)); 00616 if (A[0] == NULL) 00617 { 00618 fprintf(stderr, "Allocation failure in matrix()"); 00619 exit (1); 00620 } 00621 for(i=1; i < m; i++) 00622 A[i]=A[0] + i * n; 00623 return A; 00624 }
|
|
Definition at line 832 of file utilnum.cpp. 00833 { 00834 int i,j; 00835 for(i=0; i< m; i++){ 00836 b[i] = 0.0; 00837 for(j=0; j<n; j++){ 00838 b[i] = b[i] + A[i][j]*x[j]; 00839 } 00840 } 00841 00842 //cblas_dgemv(CblasRowMajor, CblasNoTrans, m, n, 1.0, *A, n, x, 1, 0.0, b, 1); 00843 }
|
|
Definition at line 846 of file utilnum.cpp. 00847 { 00848 int i,j; 00849 for(i=0; i< n; i++){ 00850 b[i] = 0.0; 00851 for(j=0; j<m; j++){ 00852 b[i] = b[i] + A[j][i]*x[j]; 00853 } 00854 } 00855 00856 //cblas_dgemv(CblasRowMajor, CblasTrans, m, n, 1.0, *A, m, x, 1, 0.0, b, 1); 00857 }
|
|
Definition at line 728 of file utilnum.cpp. 00729 { 00730 float *x=(float *) calloc(n, sizeof(float)); 00731 if (x == NULL) 00732 { 00733 fprintf(stderr, "Allocation failure in matrix()"); 00734 exit (1); 00735 } 00736 return x; 00737 }
|
|
|
Definition at line 662 of file utilnum.cpp. 00663 { 00664 int i,j; 00665 for (i=0; i< m; i++) 00666 { 00667 for (j=0; j< n; j++) 00668 printf ("%.15f ", A[i][j]); 00669 putchar ('\n'); 00670 } 00671 }
|
|
Definition at line 761 of file utilnum.cpp. 00762 { 00763 int i; 00764 for (i=0; i< n; i++) 00765 { 00766 printf ("%.15f \n ", x[i]); 00767 } 00768 }
|
|
Definition at line 882 of file utilnum.cpp. 00883 { 00884 int i; 00885 for (i = 0; i < m; i++) 00886 A[i][k] = v[i]; 00887 00888 }
|
|
Definition at line 702 of file utilnum.cpp. Referenced by EMAN::NormalizeRowProcessor::process_inplace(). 00703 { 00704 int i,j; 00705 for (i = 0; i < m; i++){ 00706 for (j = 0; j < n; j++) 00707 { 00708 e[i] = e[i]+ A[i][j]; 00709 } 00710 } 00711 }
|
|
Definition at line 627 of file utilnum.cpp. 00628 { 00629 int i,j; 00630 float **B = matrix(newm, newn); 00631 if (newm > m || newn > n){ 00632 printf("error in submatrix\n"); 00633 00634 return B; 00635 } 00636 for(i=0;i<newm;i++) 00637 for(j=0;j<newn;j++) 00638 B[i][j]=A[i][j]; 00639 00640 return B; 00641 }
|
|
Definition at line 740 of file utilnum.cpp. 00741 { 00742 int i; 00743 float *y = myvector(newn); 00744 if (newn > n){ 00745 printf("error in subvector"); 00746 return y; 00747 } 00748 for(i=0;i<newn;i++) 00749 y[i]=x[i]; 00750 00751 return y; 00752 }
|
|
Definition at line 780 of file utilnum.cpp. 00781 { 00782 int i; 00783 float sum = 0.0; 00784 for (i = 0; i < n; i++){ 00785 sum = sum + x[i]; 00786 } 00787 return sum; 00788 }
|
|
Definition at line 65 of file utilnum.cpp. 00067 : 00068 * 00069 * A = USV^T 00070 * 00071 * NOTE: We assume m >= n 00072 * 00073 * Input: A : m-by-n matrix 00074 * m, n : row and column dimensions of A 00075 * 00076 * Output: U : Orthogonal matrix containing left singular vectors. 00077 * s : vector containing singular values of A. 00078 * : That is, s = diag(S) 00079 * V : Orthogonal matrix containing right singular vectors. 00080 */ 00081 00082 { 00083 char jobu = 'A'; 00084 char jobv = 'A'; 00085 00086 int i; 00087 00088 // chao something funny about min/max here... hardwire lwork to 100 for now 00089 // integer lwork = max(3*min(m,n) + max(m,n), 5*min(m,n)); 00090 integer lwork = 100; 00091 integer lda=m; 00092 integer ldu=m; 00093 integer ldv=n; 00094 integer info; 00095 integer mm = m; 00096 integer nn = n; 00097 00098 float ** A_column = matrix(n,m); 00099 float * work; 00100 00101 work = myvector((int) lwork); 00102 00103 /* We have to column orient the matrix going in because dgesvd expects the 00104 * matrix to be column oriented and in C, the matrix is row oriented. */ 00105 column_orient(A, m, n, A_column); 00106 00107 /* The LAPACK routines destroy the input matrix, so use a 00108 * copy of it. 00109 * 00110 * Check to see if we are using single or float precision, 00111 * and use appropriate LAPACK routine to compute SVD. 00112 */ 00113 00114 sgesvd_(&jobu, &jobv, &mm, &nn, *A_column, &lda, s, *U, &ldu, *V, &ldv, 00115 work, &lwork, &info); 00116 00117 //sgesvd_(jobu, jobv, mm, nn, *A_column, lda, s, *U, ldu, *V, ldv, &info); 00118 00119 /* We have to transpose the U matrix to make it row oriented again. 00120 * We don't worry about V because we would have had to transpose it anyways. */ 00121 transpose_matrix(U, m, m); 00122 00123 free_vector(work); 00124 free_matrix(A_column); 00125 00126 00127 }
|
|
Definition at line 268 of file utilnum.cpp. 00272 : We assume m >= n 00273 * 00274 * Input: U : Orthogonal matrix containing left singular vectors. 00275 * s : vector containing singular values of A. 00276 * : That is, s = diag(S) 00277 * V : Orthogonal matrix containing right singular vectors. 00278 * b : Right hand side vector 00279 * m, n : row and column dimensions of A. 00280 * alpha : Regularization parameter. If alpha < 0, then 00281 * GCV is used to choose it. 00282 * omega : Omega parameter 00283 * 00284 * Output: x : Tikhonov solution 00285 * alpha : Regularization parameter 00286 */ 00287 { 00288 int i, MaxIts = 500; 00289 float * bhat = myvector(m); 00290 float * xhat = myvector(n); 00291 float * D = myvector(n); 00292 00293 matvec_multT(U, b, m, m, bhat); 00294 00295 /*printf("b:\n"); 00296 print_vector(b, m); 00297 printf("m = %d, n=%d, omega=%f\n", m,n,omega); 00298 printf("bhat:\n"); 00299 print_vector(bhat, m); 00300 printf("s:\n"); 00301 print_vector(s, n); 00302 */ 00303 00304 00305 /* 00306 * Use GCV to choose alpha. 00307 */ 00308 00309 if (*alpha < (float)0.0) 00310 *alpha = GCVmin_Tik(s, bhat, m, n, MaxIts, omega); 00311 00312 /* 00313 * Now we need to compute the the filtered Tikhonov solution 00314 */ 00315 for (i=0; i<n;i++){ 00316 D[i] = fabs(s[i])*fabs(s[i]) + (*alpha)*(*alpha); 00317 xhat[i] = (s[i] * bhat[i])/D[i]; 00318 } 00319 00320 matvec_mult(V, xhat, n, n, x); 00321 00322 free_vector(bhat); 00323 free_vector(xhat); 00324 free_vector(D); 00325 }
|
|
Definition at line 674 of file utilnum.cpp. 00675 { 00676 float save; 00677 int i,j; 00678 for (i = 0; i < m; i++){ 00679 for (j = i+1; j < n; j++) 00680 { 00681 save = A[i][j]; 00682 A[i][j] = A[j][i]; 00683 A[j][i] = save; 00684 } 00685 } 00686 }
|
|
Definition at line 130 of file utilnum.cpp. 00134 : We assume m >= n 00135 * 00136 * Input: U : Orthogonal matrix containing left singular vectors. 00137 * s : vector containing singular values of A. 00138 * : That is, s = diag(S) 00139 * V : Orthogonal matrix containing right singular vectors. 00140 * b : Right hand side vector 00141 * m, n : row and column dimensions of A. 00142 * tol : Truncation tolerance. If tol < 0, then 00143 * GCV is used to choose it. 00144 * 00145 * Output: x : TSVD solution 00146 * tol : Truncation tolerance. 00147 */ 00148 { 00149 int i, k; 00150 float * sinv = myvector(n); 00151 float * bhat = myvector(m); 00152 00153 matvec_multT(U, b, m, m, bhat); 00154 00155 //printf("m=%d, n=%d\n", m,n); 00156 /* 00157 * Use GCV to choose tol. 00158 */ 00159 k = 0; 00160 if (*tol < (float) 0.0) 00161 gcv_tsvd(s, bhat, m, n, omega, tol); 00162 00163 /* 00164 * Now we need to compute the the filtered inverse of the 00165 * singular values. That is, 1/s(i) if s(i) >= tol, and 00166 * 0 otherwise. I'm not sure this is the best way to do this 00167 * in Fortran -- it would be better to not have to check each 00168 * time through the loop. But, I do not want to assume the 00169 * the singular values are ordered -- they may not be in 00170 * some situations. Maybe a pre-sort might be helpful??? 00171 */ 00172 00173 for (i=0; i<n; i++){ 00174 if (fabs(s[i]) >= *tol){ 00175 sinv[i] = (float) 1.0 / s[i]; 00176 k++; 00177 } 00178 } 00179 *tol = k; 00180 //printf("truncindex = %f\n", *tol); 00181 //print_vector(s,n); 00182 00183 /* 00184 * Now the TSVD solution is easy to compute. 00185 */ 00186 for (i=0; i<n; i++) 00187 bhat[i] = sinv[i]*bhat[i]; 00188 00189 matvec_mult(V, bhat, n, n, x); 00190 00191 00192 free_vector(sinv); 00193 free_vector(bhat); 00194 }
|