direcmin.c

统计模式识别算法包

   xlow = x1 ;
   xhigh = x3 ;

   fbest = fsecbest = fthirdbest = current_err ;


/*
   Main loop.  For safety we impose a limit on iterations.
*/

   for (iter=0 ; iter<itmax ; iter++) {

      xmid = 0.5 * (xlow + xhigh) ;
      tol1 = tol * (fabs ( xbest ) + eps) ;
      tol2 = 2. * tol1 ;

/***************************************************************************
      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 ;
            }
         }
**************************************************************************/


/*
   The following convergence test simultaneously makes sure xhigh and
   xlow are close relative to tol2, and that xbest is near the midpoint.
*/
      if (fabs ( xbest - xmid )  <=  (tol2 - 0.5 * (xhigh - xlow)))
         break ;

      if (fabs ( prevdist )  >  tol1) {  // If we moved far enough try parabolic fit
         t1 = (xbest - xsecbest) * (fbest - fthirdbest) ; // Temps for the
         t2 = (xbest - xthirdbest) * (fbest - fsecbest) ; // parabolic estimate
         numer = (xbest - xthirdbest) * t2  -  (xbest - xsecbest) * t1 ;
         denom = 2. * (t1 - t2) ;  // Estimate will be numer / denom
         testdist = prevdist ;     // Will soon verify interval is shrinking
         prevdist = step ;         // Save for next iteration
         if (denom != 0.0)         // Avoid dividing by zero
            step = numer / denom ; // This is the parabolic estimate to min
         else 
            step = 1.e30 ;         // Assures failure of next test

         if ((fabs ( step ) < fabs ( 0.5 * testdist ))// If shrinking
          && (step + xbest > xlow)             // and within known bounds
          && (step + xbest < xhigh)) {         // then we can use the
            xrecent = xbest + step ;           // parabolic estimate
            if ((xrecent - xlow  <  tol2)  ||  // If we are very close
                (xhigh - xrecent <  tol2)) {   // to known bounds
               if (xbest < xmid)               // then stabilize
                  step = tol1 ;
               else
                  step = -tol1 ;
               }
            }
         else {  // Parabolic estimate poor, so use golden section
            prevdist = (xbest >= xmid)  ?  xlow - xbest  :  xhigh - xbest ;          // Poor so use    
            step = .3819660 * prevdist ;
            }
         }
      else { // prevdist did not exceed tol1: we did not move far enough
             // to justify a parabolic fit.  Use golden section.
         prevdist = (xbest >= xmid)  ?  xlow - xbest  :  xhigh - xbest ;
         step = .3819660 * prevdist ;
         }

      if (fabs (step)  >=  tol1)     // In order to numerically justify
         xrecent = xbest + step ;    // another trial we must move a
      else {                         // decent distance.
         if (step > 0.)
            xrecent = xbest + tol1 ;
         else
            xrecent = xbest - tol1 ;
         }

/*
   At long last we have a trial point 'xrecent'.  Evaluate the function.
*/

      step_out ( xrecent , direc , base ) ;
      frecent = trial_error ( tptr ) ;

      if (frecent <= fbest) {    // If we improved...
         if (xrecent >= xbest)   // Shrink the (xlow,xhigh) interval by
            xlow = xbest ;       // replacing the appropriate endpoint
         else
            xhigh = xbest ;
         xthirdbest = xsecbest ; // Update x and f values for best,
         xsecbest = xbest ;      // second and third best
         xbest = xrecent ;
         fthirdbest = fsecbest ;
         fsecbest = fbest ;
         fbest = frecent ;
         }

      else {                  // We did not improve
         if (xrecent < xbest) // Shrink the (xlow,xhigh) interval by
            xlow = xrecent ;  // replacing the appropriate endpoint
         else
            xhigh = xrecent ;

         if ((frecent <= fsecbest)   // If we at least beat the second best
          || (xsecbest == xbest)) {  // or we had a duplication
            xthirdbest = xsecbest ;  // we can update the second and third
            xsecbest = xrecent ;     // best, though not the best.
            fthirdbest = fsecbest ;  // Recall that we started iters with
            fsecbest = frecent ;     // best, sec and third all equal.
            }
         else if ((frecent <= fthirdbest) // Oh well.  Maybe at least we can
           || (xthirdbest == xbest)       // beat the third best or rid
           || (xthirdbest == xsecbest)) { // ourselves of a duplication
            xthirdbest = xrecent ;        // (which is how we start the
            fthirdbest = frecent ;        // iterations)
            }
         }
      }
 
   step_out ( xbest , direc , base );//Leave coefs at min
   update_dir ( xbest , direc ) ;// Make it be the actual distance moved

   if (user_quit)
      return -fbest ;
   else 
      return fbest ;
}


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

   Local routine to preserve coefs in 'base'

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

void LayerNet::preserve ( double *base )
{
   int n ;

   if (nhid1 == 0) {                // No hidden layer
      n = nout * (nin+1) ;
      memcpy ( base , out_coefs , n * sizeof(double) ) ;
      }

   else if (nhid2 == 0) {           // One hidden layer
      n = nhid1 * (nin+1) ;
      memcpy ( base , hid1_coefs , n * sizeof(double) ) ;
      base += n ;
      n = nout * (nhid1+1) ;
      memcpy ( base , out_coefs , n * sizeof(double) ) ;
      }

   else {                           // Two hidden layers
      n = nhid1 * (nin+1) ;
      memcpy ( base , hid1_coefs , n * sizeof(double) ) ;
      base += n ;
      n = nhid2 * (nhid1+1) ;
      memcpy ( base , hid2_coefs , n * sizeof(double) ) ;
      base += n ;
      n = nout * (nhid2+1) ;
      memcpy ( base , out_coefs , n * sizeof(double) ) ;
      }
}


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

   Local routine to step out from base

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

void LayerNet::step_out ( double step , double *direc , double *base )
{
   int i, n ;

   if (nhid1 == 0) {                // No hidden layer
      n = nout * (nin+1) ;
      for (i=0 ; i<n ; i++)
         out_coefs[i] = *(base++) + *(direc++) * step ;
      }

   else if (nhid2 == 0) {           // One hidden layer
      n = nhid1 * (nin+1) ;
      for (i=0 ; i<n ; i++)
         hid1_coefs[i] = *(base++) + *(direc++) * step ;
      n = nout * (nhid1+1) ;
      for (i=0 ; i<n ; i++)
         out_coefs[i] = *(base++) + *(direc++) * step ;
      }

   else {                           // Two hidden layers
      n = nhid1 * (nin+1) ;
      for (i=0 ; i<n ; i++)
         hid1_coefs[i] = *(base++) + *(direc++) * step ;
      n = nhid2 * (nhid1+1) ;
      for (i=0 ; i<n ; i++)
         hid2_coefs[i] = *(base++) + *(direc++) * step ;
      n = nout * (nhid2+1) ;
      for (i=0 ; i<n ; i++)
         out_coefs[i] = *(base++) + *(direc++) * step ;
      }
}

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

   Local routine to make 'dir' be the actual distance moved

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

void LayerNet::update_dir ( double step , double *direc )
{
   int n ;

  if (nhid1 == 0)               // No hidden layer
      n = nout * (nin+1) ;
   else if (nhid2 == 0)          // One hidden layer
      n = nhid1 * (nin+1) + nout * (nhid1+1) ;
   else                          // Two hidden layers
      n = nhid1 * (nin+1) + nhid2 * (nhid1+1) + nout * (nhid2+1) ;

   while (n--) {
      *direc *= step ;
      ++direc ;
      }
}

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

   Local routine to negate 'dir'

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

void LayerNet::negate_dir ( double *direc )
{
   int n ;

  if (nhid1 == 0)               // No hidden layer
      n = nout * (nin+1) ;
   else if (nhid2 == 0)          // One hidden layer
      n = nhid1 * (nin+1) + nout * (nhid1+1) ;
   else                          // Two hidden layers
      n = nhid1 * (nin+1) + nhid2 * (nhid1+1) + nout * (nhid2+1) ;

   while (n--) {
      *direc = -*direc ;
      ++direc ;
      }
}