direcmin.c

统计模式识别算法包

/******************************************************************************/
/*                                                                            */
/*  direcmin - Minimize along a direction                                     */
/*                                                                            */
/*  Normally this returns the mean square error, which will be 0-1.           */
/*  If the user interrupted, it returns the negative mean square error.       */
/*                                                                            */
/*  This is a two step process.  First we find three points whose center has  */
/*  the smallest function value (we bound the minimum).                       */
/*  Then we use Brent's algorithm to refine the interval.                     */
/*                                                                            */
/*  We leave 'coefs' set at the point that produced the minimum and return    */
/*  the error function at that point.  We change the direction 'dir' to       */
/*  be the actual distance moved.                                             */
/*                                                                            */
/* Copyright (c) 1993 by Academic Press, Inc.                                 */
/*                                                                            */
/* All rights reserved.  Permission is hereby granted, until further notice,  */
/* to make copies of this diskette, which are not for resale, provided these  */
/* copies are made from this master diskette only, and provided that the      */
/* following copyright notice appears on the diskette label:                  */
/* (c) 1993 by Academic Press, Inc.                                           */
/*                                                                            */
/* Except as previously stated, no part of the computer program embodied in   */
/* this diskette may be reproduced or transmitted in any form or by any means,*/
/* electronic or mechanical, including input into storage in any information  */
/* system for resale, without permission in writing from the publisher.       */
/*                                                                            */
/* Produced in the United States of America.                                  */
/*                                                                            */
/* ISBN 0-12-479041-0                                                         */
/*                                                                            */
/******************************************************************************/

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <ctype.h>
#include <stdlib.h>
#include "const.h"       // System and limitation constants, typedefs, structs
#include "classes.h"     // Includes all class headers
#include "funcdefs.h"    // Function prototypes

double LayerNet::direcmin (
   TrainingSet *tptr , // Training set to use
   double start_err ,  // Error (function value) at starting coefficients
   int itmax ,         // Upper limit on number of iterations allowed
   double eps ,        // Small, but greater than machine precision
   double tol ,        // Brent's tolerance (>= sqrt machine precision)
   double *base ,      // Work area (stepping out point)
   double *direc )     // Work area (stepping out direction)
{
   int key, user_quit, iter ;
   double step, x1, x2, x3, t1, t2, numer, denom, max_step ;
   double xlow, xhigh, xbest, testdist ;
   double current_err, err, previous_err, step_err ;
   double prevdist, etemp, frecent, fthirdbest, fsecbest, fbest ;
   double tol1, tol2, xrecent, xthirdbest, xsecbest, xmid;
   double  first_step = 2.5 ; // Heuristically found best

   user_quit = 0 ;

/*
   Take one step out in the gradient direction.  First preserve
   original weights for use as departure point parameterized by STEP.
*/

   preserve ( base ) ;   // Establishes a base for stepping out
   step_out ( first_step , direc , base ) ;
   err = trial_error ( tptr ) ;

/*
   If it increased, we had numerical problems computing the direction or
   the direction itself is too large a step.
   Negate the direction and use -1, 0 and 1.618 as first three steps.
   Otherwise use 0, 1 and 2.618 as first three steps.
*/

   if (err > start_err) {
      negate_dir ( direc ) ;
      x1 = -first_step ;
      x2 = 0. ;
      previous_err = err ;
      current_err = start_err ;
      }
   else {
      x1 = 0. ;
      x2 = first_step ;
      previous_err = start_err ;
      current_err = err ;
      }

/*
   At this point we have taken a single step and the function decreased.
   Take one more step in the golden ratio.
   Also keep errors lined up as 'previous_err', 'current_err' and 'err'.
   The corresponding abscissae will be x1, x2 and x3.
*/

/************************************************************************
   if (kbhit()) {          // Was a key pressed?
      key = getch () ;     // Read it if so
      while (kbhit())      // Flush key buffer in case function key
         getch () ;        // or key was held down
      if (key == 27)       // ESCape
         return (- err) ;
      }
***********************************************************************/

   x3 = x2 + 1.618034 * first_step ;
   step_out ( x3 , direc , base ) ;
   err = trial_error ( tptr ) ;

/*
   We now have three points x1, x2 and x3 with corresponding errors
   of 'previous_err', 'current_err' and 'err'.
   Endlessly loop until we bracket the minimum with the outer two.
*/

   while (err < current_err) {   // As long as we are descending...

/*********************************************************************
      if (kbhit()) {          // Was a key pressed?
         key = getch () ;     // Read it if so
         while (kbhit())      // Flush key buffer in case function key
            getch () ;        // or key was held down
         if (key == 27) {     // ESCape
            user_quit = 1 ;
            break ;
            }
         }

********************************************************************/

/*
   Try a parabolic fit to estimate the location of the minimum.
*/

      t1 = (x2 - x1) * (current_err - err) ;
      t2 = (x2 - x3) * (current_err - previous_err) ;
      denom = 2. * ( t2 - t1 ) ;
      if (fabs ( denom ) < eps) {
         if (denom > 0.)
            denom = eps ;
         else
            denom = -eps ;
         }
      step = x2 + ((x2 - x1) * t1  -  (x2 - x3) * t2) / denom ;//Here if perfect
      max_step = x2 + 200. * (x3 - x2) ; // Don't jump too far

      if ((x2 - step) * (step - x3)  >  0.) {         // It's between x2 and x3
         step_out ( step , direc , base ) ;
         step_err = trial_error ( tptr ) ;

         if (step_err < err) {   // It worked!  We found min between b and c.

            x1 = x2 ;
            x2 = step ;
            previous_err = current_err ;
            current_err = step_err ;
            goto BOUNDED ;
            }
         else if (step_err > current_err) { // Slight miscalc.  Min at x2.
            x3 = step ;
            err = step_err ;
            goto BOUNDED ;
            }
         else {             // Parabolic fit was total waste.  Use default.
            step = x3 + 1.618034 * (x3 - x2) ;
            step_out ( step , direc , base ) ;
            step_err = trial_error ( tptr ) ;
            }
         }

      else if ((x3 - step) * (step - max_step) > 0.0) { // Between x3 and lim
         step_out ( step , direc , base ) ;
         step_err = trial_error ( tptr ) ;
         if (step_err < err) {  // Decreased, so advance by golden ratio
            x2 = x3 ;
            x3 = step ;
            step = x3 + 1.618034 * (x3 - x2) ;
            current_err = err ;
            err = step_err ;
            step_out ( step , direc , base ) ;
            step_err = trial_error ( tptr ) ;
            }
         }

      else if ((step - max_step) * (max_step - x3)  >= 0.) {  // Beyond limit
         step = max_step ;
         step_out ( step , direc , base ) ;
         step_err = trial_error ( tptr ) ;
         if (step_err < err) {  // Decreased, so advance by golden ratio
            x2 = x3 ;
            x3 = step ;
            step = x3 + 1.618034 * (x3 - x2) ;
            current_err = err ;
            err = step_err ;
            step_out ( step , direc , base ) ;
            step_err = trial_error ( tptr ) ;
            }
         }

      else {  // Wild!  Reject parabolic and use golden ratio.
         step = x3 + 1.618034 * (x3 - x2) ;
         step_out ( step , direc , base ) ;
         step_err = trial_error ( tptr ) ;
         }

/*
   Shift three points and continue endless loop
*/

      x1 = x2 ;
      x2 = x3 ;
      x3 = step ;
      previous_err = current_err ;
      current_err = err ;
      err = step_err ;
      } // Endless stepping out loop

BOUNDED:
   step_out ( x2 , direc , base);//Leave coefs at min

   if (x1 > x3) {  // We may have switched direction at start.
      t1 = x1 ;    // Brent's method which follows assumes ordered parameter.
      x1 = x3 ;
      x3 = t1 ;
      }

   if (user_quit) {
      update_dir ( x2 , direc ) ;// Make it be the actual dist moved
      return -current_err ;
      }

/*
--------------------------------------------------------------------------------

   At this point we have bounded the minimum between x1 and x3.

   Go to the refinement stage.  We use Brent's algorithm.

--------------------------------------------------------------------------------
*/

/*
  Initialize prevdist, the distance moved on the previous step, to 0 so that
  the 'if (fabs ( prevdist )  >  tol1)' encountered on the first iteration
  below will fail, forcing a golden section the first time.  Also initialize
  step to 0 to avoid a zealous compiler from pointing out that it was
  referenced before being set.
*/

   prevdist = step = 0.0 ;

/*
   We always keep the minimum bracketed between xlow and xhigh.
   xbest has the min function so far (or latest if tie).
   xsecbest and xthirdbest are the second and third best.
*/

   xbest = xsecbest = xthirdbest = x2 ;