invtools.c

su 的源代码库

	    }
	    return(fnew);
	 }
      }
   }
   if (fOld < fnew)
   {
      restore(lambdaOld, search);
      fnew = fOld;
   }
   return(fnew);
}
/*  Function dot()                                              */
/*                                                              */
/*  Implements dot products for two vectors                     */
/*                                                              */
/*  Input parameters:                                           */
/*  a......................float vector                         */
/*  b......................float vector                         */
/*  n......................length of the vectors a and b        */
/*                                                              */
/*  Output parameters:                                          */
/*  dotP...................dot product                          */
/*                                                              */
/*                                              Wences Gouveia  */
/*                                              September 1995  */
float dot(float *a, float *b, int n)
{
   /* declaration of variables */ 
   int i;                        /* counter */
   float dotP;                   /* dot product */

   dotP = 0;
   for (i = 0; i < n; i++)
   {
      dotP += a[i] * b[i];
   }
   
   return(dotP);
}
/*  Function update()                                           */
/*                                                              */
/*  Update the model based on the search direction and          */
/*  step length                                                 */
/*                                                              */
/*  Input parameters:                                           */
/*  search.................search direction                     */
/*  step...................step length                          */
/*                                                              */
/*  Output parameters:                                          */
/*  alpha, beta or density.updated model depending on which one */
/*                         is active                            */
/*                         global                               */
/*  int....................returns 1 if the relative change in  */
/*                         all the parameters is less than      */
/*                         MIN_CHANGE                           */
/*                                              Wences Gouveia  */
/*                                              September 1995  */
int update(float step, float *search)
{
   /* declaration of variables */
   int i, k;                      /* counters */
   int flag;                      /* =1 => change in all the parameter is */
                                  /* smaller than MIN_CHANGE */
   float relChange=0;             /* relative change in parameter */
   float *a, *b, *r;              /* used in computing relative change */
   FILE *fp;                      /* report file */
   
   /* allocating memory */
   if (vpFrechet) a = alloc1float(limRange);
   if (vsFrechet) b = alloc1float(limRange);
   if (rhoFrechet) r = alloc1float(limRange);
   
   for (flag = 0, k = lim[0], i = 0; i < limRange; i++, k++)
   {
      if (vpFrechet)
      {
	 a[i] = alpha[k];
	 alpha[k] = alpha0[k] + step * search[i];
	 
	 /* implementing a hard constraint */
	 if (alpha[k] < (1 - BOUND) * alphaMean[k])
	    alpha[k] = (1 - BOUND) * alphaMean[k];
	 else if (alpha[k] > (1 + BOUND) * alphaMean[k])  
	    alpha[k] = (1 + BOUND) * alphaMean[k];

	 relChange = ABS(alpha[k] - a[i]) / a[i];
	 if (relChange < MIN_CHANGE) flag++;
      }
      
      if (vsFrechet)
      {
	 b[i] = beta[k];
	 beta[k] = beta0[k] + step * search[i + limRange];

	 /* implementing a hard constraint */
	 if (beta[k] < (1 - BOUND) * betaMean[k])
	    beta[k] = (1 - BOUND) * betaMean[k];
	 else if (beta[k] > (1 + BOUND) * betaMean[k])  
	    beta[k] = (1 + BOUND) * betaMean[k];

	 relChange = ABS(beta[k] - b[i]) / b[i];
	 if (relChange < MIN_CHANGE) flag++;
      }
      
      if (rhoFrechet)
      {
	 r[i] = rho[k];
	 rho[k] = rho0[k] + step * search[i + 2 * limRange];

	 /* implementing a hard constraint */
	 if (rho[k] < (1 - BOUND) * rhoMean[k])
	    rho[k] = (1 - BOUND) * rhoMean[k];
	 else if (rho[k] > (1 + BOUND) * rhoMean[k])  
	    rho[k] = (1 + BOUND) * rhoMean[k];

	 relChange = ABS(rho[k] - r[i]) / r[i];
	 if (relChange < MIN_CHANGE) flag++;
      }
   }

   if (flag == numberPar * limRange) 
   {
      flag = 1;
      /* reseting back */
      for (k = lim[0], i = 0; i < limRange; i++, k++)
      {
	 if (vpFrechet) alpha[k] = a[i];
	 if (vsFrechet) beta[k] = b[i];
	 if (rhoFrechet) rho[k] = r[i];
      }
   }
   else
   {
      flag = 0;
   }

   /* reporting */
   fp = fopen("report", "a");
   fprintf(fp,"---------------------------------------------------------\n");
   fprintf(fp,"Objective function evaluations: %d\n", modCount);
   fprintf(fp, "Step length: %f\n", step);
   for (i = 0; i <= info->nL; i++)
      fprintf(fp, "alpha[%d] = %6.2f beta[%d] = %6.2f rho[%d] = %6.2f\n",
	      i, alpha[i], i, beta[i], i, rho[i]);
   if (flag)
   {
      fprintf(fp, 
	      "Change in all parameters smaller than minimum! Quitting LS.\n");
   }
   fprintf(fp,"---------------------------------------------------------\n");
   fclose(fp);
   
   /* freeing memory */
   if (vpFrechet) free1float(a);
   if (vsFrechet) free1float(b);
   if (rhoFrechet) free1float(r);

   return(flag);
}
/*  Function restore()                                          */
/*                                                              */
/*  Restore model in line search if necessary                   */
/*                                                              */
/*  Input parameters:                                           */
/*  search.................search direction                     */
/*  step...................step length                          */
/*                                                              */
/*  Output parameters:                                          */
/*  alpha, beta or density.restored model depending             */
/*                         global                               */
/*                                              Wences Gouveia  */
/*                                              September 1995  */
void restore(float step, float *search)
{
   /* declaration of variables */
   int i, k;                      /* counters */
   FILE *fp;                      /* to report */

   for (k = lim[0], i = 0; i < limRange; i++, k++)
   {
      if (vpFrechet)
      {
	 alpha[k] = alpha0[k] + step * search[i];
	 /* implementing a hard constraint */
	 if (alpha[k] < (1 - BOUND) * alphaMean[k])
	    alpha[k] = (1 - BOUND) * alphaMean[k];
	 else if (alpha[k] > (1 + BOUND) * alphaMean[k])  
	    alpha[k] = (1 + BOUND) * alphaMean[k];
      }
      
      if (vsFrechet)
      {
	 beta[k] = beta0[k] + step * search[i + limRange];
	 /* implementing a hard constraint */
	 if (beta[k] < (1 - BOUND) * betaMean[k])
	    beta[k] = (1 - BOUND) * betaMean[k];
	 else if (beta[k] > (1 + BOUND) * betaMean[k])  
	    beta[k] = (1 + BOUND) * betaMean[k];
      }
      
      if (rhoFrechet)
      {
	 rho[k] = rho0[k] + step * search[i + 2 * limRange];
	 /* implementing a hard constraint */
	 if (rho[k] < (1 - BOUND) * rhoMean[k])
	    rho[k] = (1 - BOUND) * rhoMean[k];
	 else if (rho[k] > (1 + BOUND) * rhoMean[k])  
	    rho[k] = (1 + BOUND) * rhoMean[k];
      }
   }

   /* reporting */
   fp = fopen("report", "a");
   fprintf(fp,"---------------------------------------------------------\n");
   fprintf(fp,"Restoring model in line search\n");
   for (i = 0; i <= info->nL; i++)
      fprintf(fp, "alpha[%d] = %6.2f beta[%d] = %6.2f rho[%d] = %6.2f\n",
	      i, alpha[i], i, beta[i], i, rho[i]);
   fclose(fp);
}
/*  Function newSearch()                                        */
/*                                                              */
/*  Computation of the new search direction based on the        */
/*  previous (2) gradients and search direction using the       */
/*  Polak-Ribiere procedure                                     */
/*                                                              */
/*  Input parameters:                                           */
/*  grad0..................gradient at iteration i - 1          */
/*  grad1..................gradient at iteration i              */
/*  search.................direction at iteration i             */
/*                                                              */
/*  Output parameters:                                          */
/*  search.................new search direction                 */
/*                                              Wences Gouveia  */
/*                                              September 1995  */
float newSearch(float *grad0, float *grad1, float *search)
{
   /* declaration of variables */
   FILE *fp;                   /* report results */
   int i;                      /* counter */
   float descent;              /* grad1 dotted into search */
   float den, num;             /* used in computing the step length */
                               /* for search update */
   float beta;                 /* step length for search update */

   /* computing step length for search update */
   if (modCount > 1)
   {
      den = dot(grad0, grad0, numberPar * limRange);

      for (num = 0, i = 0; i < numberPar * limRange; i++)
      {
	 num += grad1[i] * (grad1[i] - grad0[i]);
      }
      
      beta = num / den;
   }
   else
   {
      beta = 0;
   }
   
   for (i = 0; i < numberPar * limRange; i++)
   {
      search[i] = -grad1[i] + beta * search[i];
   }
   
   /* checking descent direction */
   for (descent = 0, i = 0; i < numberPar * limRange; i++)
   {
      descent += grad1[i] * search[i];
   }
   
   if (descent > 0)
   {
      /* reporting reset */
      fp = fopen("report", "a");
      fprintf(fp,"-----------------------\n");
      fprintf(fp, 
	      "At iteration [OF:%d][GR:%d] search was reset to gradient:\n",
	      modCount, gradCount);
      fprintf(fp,"-----------------------\n");
      fclose(fp);
      
      /* reset search to gradient direction */
      for (i = 0; i < numberPar * limRange; i++)
      {
	 search[i] = -grad1[i];
      }
      
      /* computing descent direction */
      for (descent = 0, i = 0; i < numberPar * limRange; i++)
      {
	 descent += grad1[i] * search[i];
      }
   }

   /* reporting */
   fp = fopen("report", "a");
   fprintf(fp,"-----------------------\n");
   fprintf(fp, 
	   "At iteration [OF:%d][GR:%d] the slope is: %f:\n",
	   modCount, gradCount, descent);
   fprintf(fp,"-----------------------\n");
   fclose(fp);
   
   /* returning slope */
   return(descent);
}
/*  Function stop()                                             */
/*                                                              */
/*  Checking the stopping criterions used to verify             */
/*  convergence                                                 */
/*                                                              */
/*  Input parameters:                                           */
/*  grad0..................gradient at iteration i - 1          */
/*  grad1..................gradient at iteration i              */
/*  of0....................objective function at iteration i - 1*/
/*  of1....................objective function at iteration i    */
/*  maxIter................maximum iterations at conjugate      */    
/*                         gradient                             */
/*                                                              */
/*  Output parameters:                                          */
/*  notDone................=1 stop optimization                 */
/*                                              Wences Gouveia  */
/*                                              September 1995  */
int stop(float *grad0, float *grad1, float oF0, float oF1, int maxIter)
{
   /* declaration of variables */
   FILE *fp;                       /* to report */
   float ofChange, gradChange;     /* relative changes on objective */
                                   /* function and gradient */
   int i;                          /* counter */
   int flag;                       /* work-not-done flag */
   
   for (auxm1 = 0, auxm2 = 0, i = 0; i < limRange; i++)
   {
      auxm1 += grad0[i] * grad0[i];
      auxm2 += grad1[i] * grad1[i];
   }
   auxm1 = sqrt(auxm1); auxm2 = sqrt(auxm2);
   
   gradChange = ABS((auxm2 - auxm1) / auxm2);
   ofChange = ABS((oF1 - oF0) / oF1);
 
   if (gradChange < RELATIVE_CHANGE || ofChange < RELATIVE_CHANGE ||
       ABS(auxm2) < MIN_GRAD || gradCount > maxIter)
   {
      flag = 0;  /* convergence */
   }
   else
   {
      flag = 1;
   }
   
   /* reporting */
   fp = fopen("report", "a");
   fprintf(fp,"---------------------------------------------------------\n");
   fprintf(fp, "Number of conjugate gradient iterations: %d\n", gradCount);
   fprintf(fp, "Relative change on gradient module: %f\n", gradChange);
   fprintf(fp, "Relative change on objective function module: %f\n", 
	        ofChange);
   fprintf(fp, "Status: %d - (0) work done (1) work not done\n", flag);
   fprintf(fp,"---------------------------------------------------------\n");
   fclose(fp);
   return(flag);
}
/*  Function normalize()                                        */
/*                                                              */
/*  Normalize a vector in place                                 */
/*                                                              */
/*  Input parameters:                                           */
/*  vector.................input vector                         */
/*  int n..................length of input vector               */
/*                                                              */
/*  Output parameters:                                          */
/*  vector.................normalized vector                    */
/*                                              Wences Gouveia  */
/*                                              September 1995  */
void normalize(float *vector, int n)
{
   /* declaration of variables */
   int i;                         /* counter */
   float mod;                     /* vector module */

   for (mod = 0, i = 0; i < n; i++)
      mod += vector[i] * vector[i];
   mod = sqrt(mod);
   
   for (i = 0; i < n; i++) 
      vector[i] /= mod;
}