EMAN2
steepest.cpp
Go to the documentation of this file.
00001 /******************************************************
00002 * Program to demonstrate the use of multi-dimensional *
00003 *     Steepest Descent Optimization subroutine        *
00004 * --------------------------------------------------- *
00005 *  Reference: BASIC Scientific Subroutines, Vol. II   *
00006 *  By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1].  *
00007 *                                                     *
00008 *                  C++ version by J-P Moreau, Paris.  *
00009 * --------------------------------------------------- *
00010 * Example:   Find a local maximum of                  *
00011 *            F(x,y,z) = sin(x)+2*cos(y)-sin(z)        *
00012 *                                                     *
00013 * SAMPLE RUN:                                         *
00014 *                                                     *
00015 * How many dimensions: 3                              *
00016 *                                                     *
00017 * Convergence criterion: .000000001                   *
00018 *                                                     *
00019 * Maximum number of iterations: 50                    *
00020 *                                                     *
00021 * Starting constant: 1                                *
00022 *                                                     *
00023 * Input the starting point:                           *
00024 *                                                     *
00025 *     X[1] = 1                                        *
00026 *     X[2] = 1                                        *
00027 *     X[3] = 1                                        *
00028 *                                                     *
00029 * The results are:                                    *
00030 *                                                     *
00031 *     X(1) = 1.5707973                                *
00032 *     X(2) = -0.0000170                               *
00033 *     X(3) = -4.7123856                               *
00034 *                                                     *
00035 * Local maximum = 4.0000000                           *
00036 *                                                     *
00037 * The number of steps was: 33                         *
00038 *                                                     *
00039 ******************************************************/
00040 #include <cstdio>
00041 #include <cmath>
00042 #include "emdata.h"
00043 
00044 #define  MACHEPS  1e-15
00045 
00046 using namespace EMAN;
00047 
00048 
00049 /*******************************************
00050   Function subroutine                     */   
00051 /*double Eval1(double *X) {
00052   return sin(X[1])+2.0*cos(X[2])-sin(X[3]);
00053 }*/
00054 /*******************************************/
00055 
00056   // Functions called by Steepda()
00057   void Utilit1(double *D, double *dd, int l) {
00058     int i;
00059     // Find the magnitude of the gradient
00060     *dd=0.0;
00061     for (i=1; i<l+1; i++) *dd += D[i]*D[i];
00062     *dd=sqrt(*dd);
00063   }
00064 
00065   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,
00066   EMData *refim, EMData *mask) {
00067         int i;
00068     // Update the X[i] 
00069     for (i=1; i<l+1; i++) {
00070       // Save old values
00071       X1[i]=X[i];
00072       X[i] += xk*D[i]/(*dd);
00073     }
00074     Y[3]=(*my_func)(image, refim, mask, (float)X[1], (float)X[2], (float)X[3]);
00075   }
00076 
00077   // Find approximations of partial derivatives D(i)
00078   // by finite differences
00079   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
00080   *refim, EMData *mask)  {
00081     double a,b,yy;
00082     int i;
00083     for (i=1; i<l+1; i++) {
00084       // Save X(i)
00085       a=X[i];
00086       // Find increment
00087       b=D[i]*xk/(2.0*(*dd));
00088       // Move increment in X(i)
00089       X[i]=X[i]+b;
00090       // Obtain yy
00091       yy=(*my_func)(image, refim, mask, (float)X[1], (float)X[2], (float)X[3]);
00092       // Guard against divide by zero near maximum
00093       if (b==0) b=1e-12;
00094       // Update D(i)
00095       D[i]=(yy-Y[3])/b;
00096       // Guard against locked up derivative
00097       if (D[i]==0) D[i]=1e-5;
00098       // Restore X(i) and yy
00099       X[i]=a; yy=Y[3];
00100     }
00101     // Obtain dd
00102     Utilit1(D, dd, l);
00103   }
00104 
00105 
00106 /**************************************************
00107 *    Steepest descent optimization subroutine     *
00108 * ----------------------------------------------- *
00109 * This routine finds the local maximum or minimum *
00110 * of an L-dimensional function using the method   *
00111 * of steepest decent, or the gradient.            *
00112 * The function must be available in subroutine    *
00113 * Eval(). In this version, finite differences are *
00114 * used to calculate the L partial derivatives.    *
00115 * ----------------------------------------------- *
00116 * INPUTS:                                         *
00117 *   l - The dimension of function to study        *
00118 *   e - The convergence criteria                  *
00119 *   m - The maximum number of iterations          *
00120 *   xk - A starting constant                      *
00121 *   X(i) - Initial values of variables            *
00122 * OUTPUTS:                                        *
00123 *   X(i) - The locally optimum set                *
00124 *   Eval - The local maximum found                *
00125 *   n - The number of iterations performed,       *
00126 **************************************************/
00127   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
00128   *mask)  {
00129   // Labels: e50,e51,e100,e200
00130   int i;
00131   double dd;
00132   double  D[4], Y[4];
00133   double  X1[11];
00134   
00135   *n=0;
00136   //The routine needs three values of Y to get started
00137   //Generate starting D(i) values
00138   //These are not even good guesses and slow the program a little
00139   dd=1.0;
00140   D[1]=1.0/sqrt((double)l);
00141   for (i=2; i<l+1; i++)  D[i]=D[i-1];
00142   // Start initial probe
00143   for (i=1; i<l+1; i++) {
00144     // Obtain yy and D[i]
00145     Y[i]=(*my_func)(image, refim, mask, (float)X[1], (float)X[2], (float)X[3]);
00146     // Update X[i]
00147     Utilit1(D, &dd, l);
00148     Utilit2(X, X1, Y, D, &dd, xk, l, my_func, image, refim, mask);
00149   }
00150   // We now have a history to base the subsequent search on
00151   // Accelerate search if approach is monotonic 
00152 e50: if (fabs(Y[2]-Y[1])<MACHEPS) goto e51;
00153   if ((Y[3]-Y[2])/(Y[2]-Y[1])>0.0) xk=xk*1.2;
00154   // Decelerate if heading the wrong way
00155 e51: if (Y[3]<Y[2]) xk=xk/2.0;
00156   // Update the Y[i] if value has decreased
00157   if (Y[3]>Y[2]) goto e100;
00158   // Restore the X[i]
00159   for (i=1; i<l+1; i++) {
00160     X[i]=X1[i];
00161   }
00162   goto e200;
00163 e100: Y[1]=Y[2]; Y[2]=Y[3];
00164   // Obtain new values
00165 e200: Y[3]=(*my_func)(image, refim, mask, (float)X[1], (float)X[2], (float)X[3]);
00166   Derivatives(X, D, Y, &dd, xk, l, my_func, image, refim, mask); // Get D(i)
00167   //if dd=0 then the precision limit of the computer has been reached
00168   if (dd==0) return;
00169   // Update X[i]
00170   Utilit2(X, X1, Y, D, &dd, xk, l, my_func, image, refim, mask);
00171   // Check for maximum iterations and convergence
00172   (*n)++;
00173   //printf("Step %3d: X[0]=%12.6f  X[1]=%12.6f  X[2]=%12.6f\n",*n,X[1],X[2],X[3]);
00174   if (*n>=m) return;
00175   if (fabs(Y[3]-Y[2])<e) return;
00176   // Try another iteration
00177   goto e50;
00178 } // Steepds()
00179 
00180 /*
00181 int main() {
00182 
00183   double X[11];
00184   int i,l,m,n;
00185   double xk, e;
00186   double  (*func)(double *) = Eval1;
00187         
00188   printf("\n How many dimensions: "); scanf("%d",&l);
00189   printf("\n Convergence criterion: "); scanf("%lf",&e);
00190   printf("\n Maximum number of iterations: "); scanf("%d",&m);
00191   printf("\n Starting constant: "); scanf("%lf",&xk);
00192   printf("\n Input the starting point:\n\n");
00193   for (i=1; i<l+1; i++) {
00194     printf("     X(%d) = ",i); scanf("%lf",&X[i]);
00195   }
00196 
00197   Steepda(X, xk, e, l, m, &n, func);   // Call steepest descent optimization subroutine
00198 
00199   printf("\n\n The results are:\n\n");
00200   for (i=1; i<l+1; i++) {
00201     printf("     X(%d) = %1.7f\n",i,X[i]);   
00202   }
00203   printf("\n Local maximum found = %10.7f\n",(*func)(X));
00204   printf("\n The number of iterations was %d\n\n",n);
00205 }*/
00206 
00207 // End of file Steepda.cpp
00208 
00209   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,
00210   EMData *refim, EMData *mask, Util::KaiserBessel& kb) {
00211         int i;
00212     // Update the X[i] 
00213     for (i=1; i<l+1; i++) {
00214       // Save old values
00215       X1[i]=X[i];
00216       X[i] += xk*D[i]/(*dd);
00217     }
00218     Y[3]=(*my_func)(image, refim, mask, kb, (float)X[1], (float)X[2], (float)X[3]);
00219   }
00220 
00221   // Find approximations of partial derivatives D(i)
00222   // by finite differences
00223   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
00224   *refim, EMData *mask, Util::KaiserBessel& kb)  {
00225     double a,b,yy;
00226     int i;
00227     for (i=1; i<l+1; i++) {
00228       // Save X(i)
00229       a=X[i];
00230       // Find increment
00231       b=D[i]*xk/(2.0*(*dd));
00232       // Move increment in X(i)
00233       X[i]=X[i]+b;
00234       // Obtain yy
00235       yy=(*my_func)(image, refim, mask, kb, (float)X[1], (float)X[2], (float)X[3]);
00236       // Guard against divide by zero near maximum
00237       if (b==0) b=1e-12;
00238       // Update D(i)
00239       D[i]=(yy-Y[3])/b;
00240       // Guard against locked up derivative
00241       if (D[i]==0) D[i]=1e-5;
00242       // Restore X(i) and yy
00243       X[i]=a; yy=Y[3];
00244     }
00245     // Obtain dd
00246     Utilit1(D, dd, l);
00247   }
00248 
00249 
00250 
00251   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
00252   *mask, Util::KaiserBessel& kb)  {
00253   // Labels: e50,e51,e100,e200
00254   int i;
00255   double dd;
00256   double  D[4], Y[4];
00257   double  X1[11];
00258   
00259   *n=0;
00260   //The routine needs three values of Y to get started
00261   //Generate starting D(i) values
00262   //These are not even good guesses and slow the program a little
00263   dd=1.0;
00264   D[1]=1.0/sqrt((double)l);
00265   for (i=2; i<l+1; i++)  D[i]=D[i-1];
00266   // Start initial probe
00267   for (i=1; i<l+1; i++) {
00268     // Obtain yy and D[i]
00269     Y[i]=(*my_func)(image, refim, mask, kb, (float)X[1], (float)X[2], (float)X[3]);
00270     // Update X[i]
00271     Utilit1(D, &dd, l);
00272     Utilit2_G(X, X1, Y, D, &dd, xk, l, my_func, image, refim, mask, kb);
00273   }
00274   // We now have a history to base the subsequent search on
00275   // Accelerate search if approach is monotonic 
00276 e50: if (fabs(Y[2]-Y[1])<MACHEPS) goto e51;
00277   if ((Y[3]-Y[2])/(Y[2]-Y[1])>0.0) xk=xk*1.2;
00278   // Decelerate if heading the wrong way
00279 e51: if (Y[3]<Y[2]) xk=xk/2.0;
00280   // Update the Y[i] if value has decreased
00281   if (Y[3]>Y[2]) goto e100;
00282   // Restore the X[i]
00283   for (i=1; i<l+1; i++) {
00284     X[i]=X1[i];
00285   }
00286   goto e200;
00287 e100: Y[1]=Y[2]; Y[2]=Y[3];
00288   // Obtain new values
00289 e200: Y[3]=(*my_func)(image, refim, mask, kb, (float)X[1], (float)X[2], (float)X[3]);
00290   Derivatives_G(X, D, Y, &dd, xk, l, my_func, image, refim, mask, kb); // Get D(i)
00291   //if dd=0 then the precision limit of the computer has been reached
00292   if (dd==0) return;
00293   // Update X[i]
00294   Utilit2_G(X, X1, Y, D, &dd, xk, l, my_func, image, refim, mask, kb);
00295   // Check for maximum iterations and convergence
00296   (*n)++;
00297   //printf("Step %3d: X[0]=%12.6f  X[1]=%12.6f  X[2]=%12.6f\n",*n,X[1],X[2],X[3]);
00298   if (*n>=m) return;
00299   if (fabs(Y[3]-Y[2])<e) return;
00300   // Try another iteration
00301   goto e50;
00302 } // Steepds()