Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

hybr.h File Reference

This graph shows which files directly or indirectly include this file:

Included by dependency graph

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)


Define Documentation

 
#define frand  )     ((float) rand() / (RAND_MAX+1.0))
 

Definition at line 4 of file hybr.h.

Referenced by fill_rand().

#define integer   int
 

Definition at line 3 of file hybr.h.

Referenced by ieeeck_(), ilaenv_(), isamax_(), lsame_(), pow_ii(), pow_ri(), saxpy_(), sbdsqr_(), scopy_(), sdot_(), sgebd2_(), sgebrd_(), sgelq2_(), sgelqf_(), sgemm_(), sgemv_(), sgeqr2_(), sgeqrf_(), sger_(), sgesvd_(), slabrd_(), slacpy_(), slaed0_(), slaed1_(), slaed2_(), slaed3_(), slaed4_(), slaed6_(), slaed7_(), slaed8_(), slaed9_(), slaeda_(), slaev2_(), slamc1_(), slamc2_(), slamc4_(), slamc5_(), slamch_(), slamrg_(), slange_(), slanst_(), slansy_(), slarf_(), slarfb_(), slarfg_(), slarft_(), slartg_(), slascl_(), slaset_(), slasq1_(), slasq2_(), slasq3_(), slasq4_(), slasq5_(), slasq6_(), slasr_(), slasrt_(), slassq_(), slasv2_(), slatrd_(), snrm2_(), sorg2l_(), sorg2r_(), sorgbr_(), sorgl2_(), sorglq_(), sorgql_(), sorgqr_(), sorgtr_(), sorm2r_(), sormbr_(), sorml2_(), sormlq_(), sormqr_(), srot_(), sscal_(), sstedc_(), ssteqr_(), ssterf_(), sstevd_(), sswap_(), ssyev_(), ssymv_(), ssyr2_(), ssyr2k_(), ssytd2_(), ssytrd_(), strmm_(), strmv_(), and svd().

#define max a,
 )     a > b ? a : b
 

Definition at line 8 of file hybr.h.

Referenced by EMAN::ScaleAligner::align(), EMAN::ScaleAlignerABS::align_using_base(), EMAN::EMData::calc_fourier_shell_correlation(), EMAN::EMData::calc_max_location(), EMAN::TomoFscCmp::cmp(), EMAN::Util::constrained_helix(), EMAN::Util::constrained_helix_test(), EMAN::Processor::EMFourierFilterFunc(), EMAN::Util::find_max(), EMAN::Util::find_min_and_max(), EMAN::nnSSNR_ctfReconstructor::finish(), EMAN::nn4_ctf_rectReconstructor::finish(), EMAN::nn4_ctfReconstructor::finish(), EMAN::nnSSNR_Reconstructor::finish(), EMAN::nn4_rectReconstructor::finish(), EMAN::nn4Reconstructor::finish(), GCVmin_Tik(), EMAN::EMUtil::getRenderMinMax(), EMAN::Util::group_proj_by_phitheta(), EMAN::Util::image_mutation(), EMAN::Util::k_means_cont_table_(), main(), max2d(), max3d(), EMAN::EMData::max_3D_pixel_error(), mono(), parabld(), EMAN::HistogramBin::process_inplace(), EMAN::ClampingProcessor::process_inplace(), EMAN::NormalizeByMassProcessor::process_inplace(), EMAN::EMData::rotavg(), EMAN::MarchingCubes::setCmapMinMax(), EMAN::MrcIO::update_stat(), EMAN::EMData::update_stat(), and wustl_mm::GraySkeletonCPP::VolumeSkeletonizer::VoxelOr().

#define min a,
 )     a < b ? a : b
 

Definition at line 9 of file hybr.h.

Referenced by EMAN::ScaleAligner::align(), EMAN::ScaleAlignerABS::align_using_base(), apmq(), EMAN::EMData::calc_min_location(), EMAN::Util::constrained_helix(), EMAN::Util::constrained_helix_test(), EMAN::Util::find_min_and_max(), EMAN::BoxSVDClassifier::getMapping(), EMAN::EMUtil::getRenderMinMax(), EMAN::Util::histc(), EMAN::Util::image_mutation(), mono(), parabld(), EMAN::HistogramBin::process_inplace(), EMAN::ClampingProcessor::process_inplace(), EMAN::TestImageCylinder::process_inplace(), EMAN::TestImageAxes::process_inplace(), EMAN::NormalizeByMassProcessor::process_inplace(), recons3d_HyBR_mpi_Cart(), EMAN::EMData::rotavg(), EMAN::EMData::rotavg_i(), EMAN::MarchingCubes::setCmapMinMax(), EMAN::MrcIO::update_stat(), and EMAN::EMData::update_stat().


Function Documentation

void alphaxpy float  alpha,
float *  x,
float *  y,
int  n
 

Definition at line 822 of file utilnum.cpp.

00823 { 
00824  int i;
00825  for (i=0; i<n; i++)
00826    y[i] = alpha*x[i] + y[i];
00827 }

void column_orient float **  A,
int  m,
int  n,
float **  B
 

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 }

void copy_matrix float **  A,
float **  B,
int  m,
int  n
 

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 }

void copy_vector float *  x,
float *  y,
int  n
 

Definition at line 791 of file utilnum.cpp.

00792 {
00793   int i;
00794   for(i = 0; i < n; i++) 
00795       y[i]=x[i];    
00796 }

void fill_ones float *  x,
int  n
 

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 }

void fill_rand float **  A,
int  m,
int  n
 

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 }

void fillveczeros float *  x,
int  n
 

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 }

void free_matrix float **  A  ) 
 

Definition at line 644 of file utilnum.cpp.

00645 {
00646   free(A[0]);
00647   free(A);
00648 }

void free_vector float *  x  ) 
 

Definition at line 755 of file utilnum.cpp.

00756 {
00757   free(x);
00758 }

void gcv_tsvd float *  s,
float *  bhat,
int  m,
int  n,
float  omega,
float *  tol
 

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 }

float GCVfun_Tik float  alpha,
float *  s,
float *  bhat,
int  m,
int  n,
float  omega
 

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 }

float GCVmin_Tik float *  s,
float *  bhat,
int  m,
int  n,
int  MaxIts,
float  omega
 

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 }

void generate_noise double *  E,
int  m,
int  n
 

void get_col float **  A,
int  m,
int  k,
float *  ak
 

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 }

void get_row float **  A,
int  n,
int  k,
float *  ak
 

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 }

float** matrix int  m,
int  n
 

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 }

void matvec_mult float **  A,
float *  x,
int  m,
int  n,
float *  b
 

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 }

void matvec_multT float **  A,
float *  x,
int  m,
int  n,
float *  b
 

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 }

float* myvector int  n  ) 
 

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 }

float norm float *  x,
int  n
 

Definition at line 807 of file utilnum.cpp.

Referenced by EMAN::EMData::calc_center_density(), EMAN::EMData::calc_radial_dist(), EMAN::Util::cml_line_in3d(), EMAN::FRCCmp::cmp(), EMAN::PhaseCmp::cmp(), EMAN::TomoCccCmp::cmp(), EMAN::DoGFourierProcessor::create_radial_func(), EMAN::FourierWeightAverager::finish(), EMAN::FourierPixelInserter3D::FourierPixelInserter3D(), EMAN::PointArray::get_center(), EMAN::BoxSVDClassifier::go(), main(), EMAN::EMData::make_footprint(), EMAN::ReconstructorVolumeData::normalize_threed(), EMAN::PointArray::pdb2mrc_by_nfft(), EMAN::WienerFourierReconstructor::pixel_at(), EMAN::FourierReconstructor::pixel_at(), EMAN::TestImagePureGaussian::process_inplace(), EMAN::RealToFFTProcessor::process_inplace(), EMAN::PointArray::projection_by_nfft(), slange_(), slanst_(), slansy_(), snrm2_(), and EMAN::EMData::unwrap_largerR().

00808 {
00809   float x_norm = 0.0;
00810   int i;
00811   for (i=0; i<n; i++)
00812         x_norm = x_norm + x[i]*x[i];
00813 
00814   x_norm = sqrt(x_norm);
00815 
00816 //x_norm = cblas_dnrm2(n, x, 1);
00817 
00818   return x_norm;
00819 }

void print_matrix float **  A,
int  m,
int  n
 

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 }

void print_vector float *  x,
int  n
 

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 }

void put_col float **  A,
int  m,
float *  v,
int  k
 

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 }

void row_sum float **  A,
int  m,
int  n,
float *  e
 

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 }

float** submatrix float **  A,
int  m,
int  n,
int  newm,
int  newn
 

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 }

float* subvector float *  x,
int  n,
int  newn
 

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 }

float sum_vector float *  x,
int  n
 

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 }

void svd float **  A,
int  m,
int  n,
float **  U,
float *  x,
float **  V
 

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 }

void tikhonov float **  U,
float *  s,
float **  V,
float *  b,
int  m,
int  n,
float  omega,
float *  x,
float *  alpha
 

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 }

void transpose_matrix float **  A,
int  m,
int  n
 

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 }

void tsvd float **  U,
float *  s,
float **  V,
float *  b,
int  m,
int  n,
float  omega,
float *  x,
float *  tol
 

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 }


Generated on Tue Jun 11 13:46:51 2013 for EMAN2 by  doxygen 1.3.9.1