EMAN2
Defines | Functions
steepest.cpp File Reference
#include <cstdio>
#include <cmath>
#include "emdata.h"
Include dependency graph for steepest.cpp:

Go to the source code of this file.

Defines

#define MACHEPS   1e-15

Functions

void Utilit1 (double *D, double *dd, int l)
void Utilit2 (double *X, double *X1, double *Y, double *D, double *dd, double xk, int l, float(*my_func)(EMData *, EMData *, EMData *, float, float, float), EMData *image, EMData *refim, EMData *mask)
void Derivatives (double *X, double *D, double *Y, double *dd, double xk, int l, float(*my_func)(EMData *, EMData *, EMData *, float, float, float), EMData *image, EMData *refim, EMData *mask)
void Steepda (double *X, double xk, double e, int l, int m, int *n, float(*my_func)(EMData *, EMData *, EMData *, float, float, float), EMData *image, EMData *refim, EMData *mask)
void Utilit2_G (double *X, double *X1, double *Y, double *D, double *dd, double xk, int l, float(*my_func)(EMData *, EMData *, EMData *, Util::KaiserBessel &, float, float, float), EMData *image, EMData *refim, EMData *mask, Util::KaiserBessel &kb)
void Derivatives_G (double *X, double *D, double *Y, double *dd, double xk, int l, float(*my_func)(EMData *, EMData *, EMData *, Util::KaiserBessel &, float, float, float), EMData *image, EMData *refim, EMData *mask, Util::KaiserBessel &kb)
void Steepda_G (double *X, double xk, double e, int l, int m, int *n, float(*my_func)(EMData *, EMData *, EMData *, Util::KaiserBessel &, float, float, float), EMData *image, EMData *refim, EMData *mask, Util::KaiserBessel &kb)

Define Documentation

#define MACHEPS   1e-15

Definition at line 44 of file steepest.cpp.

Referenced by Steepda(), and Steepda_G().


Function Documentation

void Derivatives ( double *  X,
double *  D,
double *  Y,
double *  dd,
double  xk,
int  l,
float(*)(EMData *, EMData *, EMData *, float, float, float)  my_func,
EMData image,
EMData refim,
EMData mask 
)

Definition at line 79 of file steepest.cpp.

References b, and Utilit1().

Referenced by Steepda().

                         {
    double a,b,yy;
    int i;
    for (i=1; i<l+1; i++) {
      // Save X(i)
      a=X[i];
      // Find increment
      b=D[i]*xk/(2.0*(*dd));
      // Move increment in X(i)
      X[i]=X[i]+b;
      // Obtain yy
      yy=(*my_func)(image, refim, mask, (float)X[1], (float)X[2], (float)X[3]);
      // Guard against divide by zero near maximum
      if (b==0) b=1e-12;
      // Update D(i)
      D[i]=(yy-Y[3])/b;
      // Guard against locked up derivative
      if (D[i]==0) D[i]=1e-5;
      // Restore X(i) and yy
      X[i]=a; yy=Y[3];
    }
    // Obtain dd
    Utilit1(D, dd, l);
  }
void Derivatives_G ( double *  X,
double *  D,
double *  Y,
double *  dd,
double  xk,
int  l,
float(*)(EMData *, EMData *, EMData *, Util::KaiserBessel &, float, float, float)  my_func,
EMData image,
EMData refim,
EMData mask,
Util::KaiserBessel kb 
)

Definition at line 223 of file steepest.cpp.

References b, and Utilit1().

Referenced by Steepda_G().

                                               {
    double a,b,yy;
    int i;
    for (i=1; i<l+1; i++) {
      // Save X(i)
      a=X[i];
      // Find increment
      b=D[i]*xk/(2.0*(*dd));
      // Move increment in X(i)
      X[i]=X[i]+b;
      // Obtain yy
      yy=(*my_func)(image, refim, mask, kb, (float)X[1], (float)X[2], (float)X[3]);
      // Guard against divide by zero near maximum
      if (b==0) b=1e-12;
      // Update D(i)
      D[i]=(yy-Y[3])/b;
      // Guard against locked up derivative
      if (D[i]==0) D[i]=1e-5;
      // Restore X(i) and yy
      X[i]=a; yy=Y[3];
    }
    // Obtain dd
    Utilit1(D, dd, l);
  }
void Steepda ( double *  X,
double  xk,
double  e,
int  l,
int  m,
int *  n,
float(*)(EMData *, EMData *, EMData *, float, float, float)  my_func,
EMData image,
EMData refim,
EMData mask 
)

Definition at line 127 of file steepest.cpp.

References Derivatives(), MACHEPS, sqrt(), Utilit1(), and Utilit2().

Referenced by EMAN::Util::twoD_fine_ali_SD().

          {
  // Labels: e50,e51,e100,e200
  int i;
  double dd;
  double  D[4], Y[4];
  double  X1[11];
  
  *n=0;
  //The routine needs three values of Y to get started
  //Generate starting D(i) values
  //These are not even good guesses and slow the program a little
  dd=1.0;
  D[1]=1.0/sqrt((double)l);
  for (i=2; i<l+1; i++)  D[i]=D[i-1];
  // Start initial probe
  for (i=1; i<l+1; i++) {
    // Obtain yy and D[i]
    Y[i]=(*my_func)(image, refim, mask, (float)X[1], (float)X[2], (float)X[3]);
    // Update X[i]
    Utilit1(D, &dd, l);
    Utilit2(X, X1, Y, D, &dd, xk, l, my_func, image, refim, mask);
  }
  // We now have a history to base the subsequent search on
  // Accelerate search if approach is monotonic 
e50: if (fabs(Y[2]-Y[1])<MACHEPS) goto e51;
  if ((Y[3]-Y[2])/(Y[2]-Y[1])>0.0) xk=xk*1.2;
  // Decelerate if heading the wrong way
e51: if (Y[3]<Y[2]) xk=xk/2.0;
  // Update the Y[i] if value has decreased
  if (Y[3]>Y[2]) goto e100;
  // Restore the X[i]
  for (i=1; i<l+1; i++) {
    X[i]=X1[i];
  }
  goto e200;
e100: Y[1]=Y[2]; Y[2]=Y[3];
  // Obtain new values
e200: Y[3]=(*my_func)(image, refim, mask, (float)X[1], (float)X[2], (float)X[3]);
  Derivatives(X, D, Y, &dd, xk, l, my_func, image, refim, mask); // Get D(i)
  //if dd=0 then the precision limit of the computer has been reached
  if (dd==0) return;
  // Update X[i]
  Utilit2(X, X1, Y, D, &dd, xk, l, my_func, image, refim, mask);
  // Check for maximum iterations and convergence
  (*n)++;
  //printf("Step %3d: X[0]=%12.6f  X[1]=%12.6f  X[2]=%12.6f\n",*n,X[1],X[2],X[3]);
  if (*n>=m) return;
  if (fabs(Y[3]-Y[2])<e) return;
  // Try another iteration
  goto e50;
} // Steepds()
void Steepda_G ( double *  X,
double  xk,
double  e,
int  l,
int  m,
int *  n,
float(*)(EMData *, EMData *, EMData *, Util::KaiserBessel &, float, float, float)  my_func,
EMData image,
EMData refim,
EMData mask,
Util::KaiserBessel kb 
)

Definition at line 251 of file steepest.cpp.

References Derivatives_G(), MACHEPS, sqrt(), Utilit1(), and Utilit2_G().

Referenced by EMAN::Util::twoD_fine_ali_SD_G().

                                {
  // Labels: e50,e51,e100,e200
  int i;
  double dd;
  double  D[4], Y[4];
  double  X1[11];
  
  *n=0;
  //The routine needs three values of Y to get started
  //Generate starting D(i) values
  //These are not even good guesses and slow the program a little
  dd=1.0;
  D[1]=1.0/sqrt((double)l);
  for (i=2; i<l+1; i++)  D[i]=D[i-1];
  // Start initial probe
  for (i=1; i<l+1; i++) {
    // Obtain yy and D[i]
    Y[i]=(*my_func)(image, refim, mask, kb, (float)X[1], (float)X[2], (float)X[3]);
    // Update X[i]
    Utilit1(D, &dd, l);
    Utilit2_G(X, X1, Y, D, &dd, xk, l, my_func, image, refim, mask, kb);
  }
  // We now have a history to base the subsequent search on
  // Accelerate search if approach is monotonic 
e50: if (fabs(Y[2]-Y[1])<MACHEPS) goto e51;
  if ((Y[3]-Y[2])/(Y[2]-Y[1])>0.0) xk=xk*1.2;
  // Decelerate if heading the wrong way
e51: if (Y[3]<Y[2]) xk=xk/2.0;
  // Update the Y[i] if value has decreased
  if (Y[3]>Y[2]) goto e100;
  // Restore the X[i]
  for (i=1; i<l+1; i++) {
    X[i]=X1[i];
  }
  goto e200;
e100: Y[1]=Y[2]; Y[2]=Y[3];
  // Obtain new values
e200: Y[3]=(*my_func)(image, refim, mask, kb, (float)X[1], (float)X[2], (float)X[3]);
  Derivatives_G(X, D, Y, &dd, xk, l, my_func, image, refim, mask, kb); // Get D(i)
  //if dd=0 then the precision limit of the computer has been reached
  if (dd==0) return;
  // Update X[i]
  Utilit2_G(X, X1, Y, D, &dd, xk, l, my_func, image, refim, mask, kb);
  // Check for maximum iterations and convergence
  (*n)++;
  //printf("Step %3d: X[0]=%12.6f  X[1]=%12.6f  X[2]=%12.6f\n",*n,X[1],X[2],X[3]);
  if (*n>=m) return;
  if (fabs(Y[3]-Y[2])<e) return;
  // Try another iteration
  goto e50;
} // Steepds()
void Utilit1 ( double *  D,
double *  dd,
int  l 
)

Definition at line 57 of file steepest.cpp.

References sqrt().

Referenced by Derivatives(), Derivatives_G(), Steepda(), and Steepda_G().

                                             {
    int i;
    // Find the magnitude of the gradient
    *dd=0.0;
    for (i=1; i<l+1; i++) *dd += D[i]*D[i];
    *dd=sqrt(*dd);
  }
void Utilit2 ( double *  X,
double *  X1,
double *  Y,
double *  D,
double *  dd,
double  xk,
int  l,
float(*)(EMData *, EMData *, EMData *, float, float, float)  my_func,
EMData image,
EMData refim,
EMData mask 
)

Definition at line 65 of file steepest.cpp.

Referenced by Steepda().

                               {
        int i;
    // Update the X[i] 
    for (i=1; i<l+1; i++) {
      // Save old values
      X1[i]=X[i];
      X[i] += xk*D[i]/(*dd);
    }
    Y[3]=(*my_func)(image, refim, mask, (float)X[1], (float)X[2], (float)X[3]);
  }
void Utilit2_G ( double *  X,
double *  X1,
double *  Y,
double *  D,
double *  dd,
double  xk,
int  l,
float(*)(EMData *, EMData *, EMData *, Util::KaiserBessel &, float, float, float)  my_func,
EMData image,
EMData refim,
EMData mask,
Util::KaiserBessel kb 
)

Definition at line 209 of file steepest.cpp.

Referenced by Steepda_G().

                                                     {
        int i;
    // Update the X[i] 
    for (i=1; i<l+1; i++) {
      // Save old values
      X1[i]=X[i];
      X[i] += xk*D[i]/(*dd);
    }
    Y[3]=(*my_func)(image, refim, mask, kb, (float)X[1], (float)X[2], (float)X[3]);
  }