asamin.c

模拟退化算法在C环境下的实现

/***********************************************************************
* MATLAB Gateway Routine for Lester Ingber's Adaptive Simulated
* Annealing (ASA)
* 
* Copyright (c) 1999-2001 Shinichi Sakata.  All Rights Reserved.
***********************************************************************/
/* $Id: asamin.c,v 1.24 2002/12/25 20:10:50 ssakata Exp ssakata $ */

#include <string.h>
#include "mex.h"
#include "asa.h"
#include "asamin.h"

/* Some working storage */
static USER_DEFINES *USER_ASA_OPTIONS;
static char cost_func_name[MAXLEN_COST_FUNC_NAME];
static int use_rejected_cost[1] = {FALSE};
static double (*cost_function) ();
static int matlab_cost_func_nrhs = 1;
static mxArray **matlab_cost_func_prhs;
static double *critical_cost_value;
static double *user_acceptance_flag;

static void
set_real_to_option (const mxArray *value,
		    double *var)
{
  if ( (mxGetM (value) == 1) && (mxGetN (value) == 1)
       && (mxIsDouble (value) ) && (!mxIsComplex (value)) )
    {
      *var = mxGetScalar (value);
    } else {
      mexErrMsgTxt ("Error: The second operand must be real.");
    }
}

static void
set_int_to_option (const mxArray *value,
		   int *var)
{
  if ( (mxGetM (value) == 1) && (mxGetN (value)==1)
       && (mxIsDouble (value)) && (!mxIsComplex (value)) )
    {
      *var = mxGetScalar (value);
    } else {
      mexErrMsgTxt ("Error: The second operand must be real.");
    }
}

static void
set_longint_to_option (const mxArray *value,
		       long int *var)
{
  if ( (mxGetM (value) == 1) && (mxGetN (value)==1)
       && (mxIsDouble (value)) && (!mxIsComplex (value)) )
    {
      *var = mxGetScalar (value);
    } else {
      mexErrMsgTxt ("Error: The second operand must be real.");
    }
}

static void
set_string_to_option (const mxArray *value,
		      char *var)
{
  if ( mxIsChar (value) )
    {
      mxGetString (value, var, mxGetN (value) + 1);
    } else {
      mexErrMsgTxt ("Error: The second operand must be a character string.");
    }
}
      
/*
  The next two functions, myrand and randflt, were copied from
  user.c in ASA.
*/

#define MULT ((LONG_INT) 25173)
#define MOD ((LONG_INT) 65536)
#define INCR ((LONG_INT) 13849)
#define FMOD ((double) 65536.0)

/***********************************************************************
* double myrand - returns random number between 0 and 1
*	This routine returns the random number generator between 0 and 1
***********************************************************************/

static double
myrand (LONG_INT * rand_seed)
{
#if TRUE			/* (change to FALSE for alternative RNG) */
  *rand_seed = (LONG_INT) ((MULT * (*rand_seed) + INCR) % MOD);
  return ((double) (*rand_seed) / FMOD);
#else
  /* See "Random Number Generators: Good Ones Are Hard To Find,"
     Park & Miller, CACM 31 (10) (October 1988) pp. 1192-1201.
     ***********************************************************
     THIS IMPLEMENTATION REQUIRES AT LEAST 32 BIT INTEGERS
     *********************************************************** */
#define _A_MULTIPLIER  16807L
#define _M_MODULUS     2147483647L	/* (2**31)-1 */
#define _Q_QUOTIENT    127773L	/* 2147483647 / 16807 */
#define _R_REMAINDER   2836L	/* 2147483647 % 16807 */
  long lo;
  long hi;
  long test;

  hi = *rand_seed / _Q_QUOTIENT;
  lo = *rand_seed % _Q_QUOTIENT;
  test = _A_MULTIPLIER * lo - _R_REMAINDER * hi;
  if (test > 0)
    {
      *rand_seed = test;
    }
  else
    {
      *rand_seed = test + _M_MODULUS;
    }
  return ((double) *rand_seed / _M_MODULUS);
#endif /* alternative RNG */
}

/***********************************************************************
* double randflt
***********************************************************************/

static double
randflt (LONG_INT *rand_seed)
{
  return (resettable_randflt (rand_seed, 0));
}

/***********************************************************************
* double resettable_randflt
***********************************************************************/
static double
resettable_randflt (LONG_INT * rand_seed, int reset)

  /* shuffles random numbers in random_array[SHUFFLE] array */

{

  /* This RNG is a modified algorithm of that presented in
   * %A K. Binder
   * %A D. Stauffer
   * %T A simple introduction to Monte Carlo simulations and some
   *    specialized topics
   * %B Applications of the Monte Carlo Method in statistical physics
   * %E K. Binder
   * %I Springer-Verlag
   * %C Berlin
   * %D 1985
   * %P 1-36
   * where it is stated that such algorithms have been found to be
   * quite satisfactory in many statistical physics applications. */

  double rranf;
  unsigned kranf;
  int n;
  static int randflt_initial_flag = 0;
  LONG_INT initial_seed;
  static double random_array[SHUFFLE];	/* random variables */

  if (*rand_seed < 0)
    *rand_seed = -*rand_seed;

  if ((randflt_initial_flag == 0) || reset)
    {
      initial_seed = *rand_seed;

      for (n = 0; n < SHUFFLE; ++n)
	random_array[n] = myrand (&initial_seed);

      randflt_initial_flag = 1;

      for (n = 0; n < 1000; ++n)	/* warm up random generator */
	rranf = randflt (&initial_seed);

      rranf = randflt (rand_seed);

      return (rranf);
    }

  kranf = (unsigned) (myrand (rand_seed) * SHUFFLE) % SHUFFLE;
  rranf = *(random_array + kranf);
  *(random_array + kranf) = myrand (rand_seed);

  return (rranf);
}

static double
cost_function_without_test (double *x,
			    double *parameter_lower_bound,
			    double *parameter_upper_bound,
			    double *cost_tangents,
			    double *cost_curvature,
			    ALLOC_INT *parameter_dimension,
			    int *parameter_int_real,
			    int *cost_flag,
			    int *exit_code,
			    USER_DEFINES *USER_OPTIONS)
{
  double unif_rn;
  double cost_value;
  
  /* The cost function written in matlab should return
     the cost value */

  int nlhs=2;
  mxArray *plhs[2];

  
  /* The cost function written in matlab is called with arguments:
     x */

  /* Set matlab_cost_func_prhs[0] */
  memcpy (mxGetPr (matlab_cost_func_prhs[0]), x,
	  *parameter_dimension * sizeof (double));

  *user_acceptance_flag = USER_OPTIONS->User_Acceptance_Flag;  
  
  mexCallMATLAB (nlhs, plhs,
		 matlab_cost_func_nrhs,
		 matlab_cost_func_prhs,
		 cost_func_name);
  
  cost_value = mxGetScalar (plhs[0]);
  *cost_flag = mxGetScalar (plhs[1]);

  if( !USER_OPTIONS->User_Acceptance_Flag )
    {
      unif_rn = randflt (USER_OPTIONS->Random_Seed);
      *critical_cost_value =
	*USER_OPTIONS->Last_Cost -
	USER_OPTIONS->Cost_Temp_Curr * log (unif_rn + DBL_MIN);
      if (cost_value <= *critical_cost_value)
	{
	  USER_OPTIONS->User_Acceptance_Flag = TRUE;
	}
      else
	{
	  USER_OPTIONS->User_Acceptance_Flag = FALSE;
	}
    }

  /* If the cost value less than the past best value is returned to
    asa (), asa () updates various indices using the returned cost
    value, whether the current state is accepted or not. This feature
    may have some values in certain applications, but it also may be
    confusing.  Here, we set the *critical_cost_value to cost_value if
    the current state is rejected, so that the cost value returned to
    asa () is never less than the past best value if the current state
    is rejected. If you want to return the cost value computed by the
    user matlab function instead, turn off this feature by replacing
    TRUE in the following line with FALSE.  */

  if (!*use_rejected_cost)
    {
      if ((! *user_acceptance_flag) &&
	  *cost_flag &&
	  (! USER_OPTIONS->User_Acceptance_Flag))
	{
	  cost_value = *critical_cost_value;
	}
    }

  /* This gateway routine assumes that the user cost function can
     always determine if the current stage should be accepted; so,
     USER_OPTIONS->Cost_Acceptance_Flag is always set to TRUE. */
  USER_OPTIONS->Cost_Acceptance_Flag = TRUE;
  mxDestroyArray (plhs[0]);
  mxDestroyArray (plhs[1]);
  return (cost_value);
}

static double
cost_function_with_test (double *x,
			 double *parameter_lower_bound,
			 double *parameter_upper_bound,
			 double *cost_tangents,
			 double *cost_curvature,
			 ALLOC_INT *parameter_dimension,
			 int *parameter_int_real,
			 int *cost_flag,
			 int *exit_code,
			 USER_DEFINES *USER_OPTIONS)
{
  double unif_rn;
  double cost_value;
  
  /* The cost function written in matlab should return
     the cost value
     cost_flag (the constraints are valid?)
     user_acceptance_flag (the current state is accepted?) */

  int nlhs=3;
  mxArray *plhs[3];

  
  /* The cost function written in matlab is called with arguments:
     x
     critical_cost_value
     User_Acceptance_Flag (=1 if acceptancetest is not necessary)  */

  /* Set matlab_cost_func_prhs[0] */
  memcpy (mxGetPr (matlab_cost_func_prhs[0]), x,
	  *parameter_dimension * (ALLOC_INT) sizeof (double));

  /* Set matlab_cost_func_prhs[1] */

  if( !USER_OPTIONS->User_Acceptance_Flag )
    {
      unif_rn = randflt (USER_OPTIONS->Random_Seed);
      *critical_cost_value =
	*USER_OPTIONS->Last_Cost -
	USER_OPTIONS->Cost_Temp_Curr * log (unif_rn + DBL_MIN);
    }

  /* Set matlab_cost_func_prhs[2] */
  *user_acceptance_flag = USER_OPTIONS->User_Acceptance_Flag;  
       
  mexCallMATLAB (nlhs, plhs,
		 matlab_cost_func_nrhs,
		 matlab_cost_func_prhs,
		 cost_func_name);
  
  cost_value = mxGetScalar (plhs[0]);
  *cost_flag = mxGetScalar (plhs[1]);
  USER_OPTIONS->User_Acceptance_Flag = mxGetScalar (plhs[2]);

  /* If the cost value less than the past best value is returned to
    asa (), asa () updates various indices using the returned cost
    value, whether the current state is accepted or not. This feature
    may have some values in certain applications, but it also may be
    confusing.  Here, we set the *critical_cost_value to cost_value if
    the current state is rejected, so that the cost value returned to
    asa () is never less than the past best value if the current state
    is rejected. If you want to return the cost value computed by the
    user matlab function instead, turn off this feature by replacing
    TRUE in the following line with FALSE.  */

  if (!*use_rejected_cost)
    {
      if ((! *user_acceptance_flag) &&
	  *cost_flag &&
	  (! USER_OPTIONS->User_Acceptance_Flag))
	{
	  cost_value = *critical_cost_value;
	}
    }

  /* This gateway routine assumes that the user cost function can
     always determine if the current stage should be accepted; so,
     USER_OPTIONS->Cost_Acceptance_Flag is always set to TRUE. */
  USER_OPTIONS->Cost_Acceptance_Flag = TRUE;
  mxDestroyArray (plhs[0]);
  mxDestroyArray (plhs[1]);
  mxDestroyArray (plhs[2]);
  return (cost_value);
}

static void
user_acceptance_test (double current_cost,
		      ALLOC_INT * parameter_dimension,
		      USER_DEFINES * USER_OPTIONS)
{
  mexErrMsgTxt ("Internal Error: user_acceptance_test () is called.");
}

/* An exit function to clean up the working storage.
   Registered with matlab and called by matlab on exit */


static void
reset_options (LONG_INT * rand_seed,
	       int * test_in_cost_func,
	       int * use_rejected_cost,
	       USER_DEFINES * USER_OPTIONS)
{
  *rand_seed = 696969;
  *test_in_cost_func = 0;
  *use_rejected_cost = FALSE;
  /* USER_OPTIONS->Limit_Acceptances = 10000; */
  USER_OPTIONS->Limit_Acceptances = 1000;
  USER_OPTIONS->Limit_Generated = 99999;
  USER_OPTIONS->Limit_Invalid_Generated_States = 1000;
  /* USER_OPTIONS->Accepted_To_Generated_Ratio = 1.0E-6; */
  USER_OPTIONS->Accepted_To_Generated_Ratio = 1.0E-4;

  USER_OPTIONS->Cost_Precision = 1.0E-18;
  USER_OPTIONS->Maximum_Cost_Repeat = 5;
  USER_OPTIONS->Number_Cost_Samples = 5;
  USER_OPTIONS->Temperature_Ratio_Scale = 1.0E-5;
  USER_OPTIONS->Cost_Parameter_Scale_Ratio = 1.0;
  USER_OPTIONS->Temperature_Anneal_Scale = 100.0;

  USER_OPTIONS->Include_Integer_Parameters = FALSE;
  USER_OPTIONS->User_Initial_Parameters = FALSE;
  USER_OPTIONS->Sequential_Parameters = -1;
  USER_OPTIONS->Initial_Parameter_Temperature = 1.0;

  USER_OPTIONS->Acceptance_Frequency_Modulus = 100;
  USER_OPTIONS->Generated_Frequency_Modulus = 10000;
  USER_OPTIONS->Reanneal_Cost = 1;
  USER_OPTIONS->Reanneal_Parameters = TRUE;

  USER_OPTIONS->Delta_X = 0.001;
  USER_OPTIONS->User_Tangents = FALSE;
  USER_OPTIONS->Curvature_0 = FALSE;
  strcpy (USER_OPTIONS->Asa_Out_File, "asa.log");
}

void
mexFunction (int nlhs, mxArray *plhs[],
	     int nrhs, const mxArray *prhs[])
{
  char cmd[MAXLEN_CMD];
  int cmdlen;
  static LONG_INT rand_seed[1] = {696969};
  static int test_in_cost_func[1] = {0};

  static short initialized = FALSE;

  int option_matched = FALSE;
  int i;
  double *ptr;

  int exit_code[1] = {0};
  double *parameter_lower_bound, *parameter_upper_bound,
    *cost_parameters, *cost_tangents, *cost_curvature;
  double *minimum_cost_value;

  /* the number of parameters to optimize */
  ALLOC_INT *parameter_dimension;

  /* pointer to array storage for parameter type flags */
  int *parameter_int_real;

  /* valid flag for cost function */
  int cost_flag[1] = {0};

  if (!initialized) {
    /* Register the Exit Function */

    if (0 != mexAtExit(exit_function)) {
      mexErrMsgTxt("Internal Error: Failed to register exit function.\n");
    }
    initialized = TRUE;
    if ((USER_ASA_OPTIONS =
	 (USER_DEFINES *) mxCalloc (1, sizeof (USER_DEFINES))) == NULL)
      {
	mexErrMsgTxt ("Internal Error: USER_DEFINES cannot be allocated.");
      }
    mexMakeMemoryPersistent (USER_ASA_OPTIONS);
    if ((USER_ASA_OPTIONS->Asa_Out_File =
	 (char *) mxCalloc (MAXLEN_ASA_OUT_FILE, sizeof (char))
	 ) == NULL)
      {
	mexErrMsgTxt ("Internal Error: USER_ASA_OPTIONS->Asa_Out_File cannot be allocated.");
      }
    mexMakeMemoryPersistent (USER_ASA_OPTIONS->Asa_Out_File);

    reset_options (rand_seed, test_in_cost_func,
		   use_rejected_cost, USER_ASA_OPTIONS);
    USER_ASA_OPTIONS->Acceptance_Test = user_acceptance_test;
  }

  if (nrhs<1) {
    mexErrMsgTxt ("Error: No command was given.");
  }

  if (!mxIsChar(prhs[0]))
    {
      mexErrMsgTxt("Error: The first argument must be a commond.");
    }
  if ((cmdlen = mxGetN(prhs[0]) + 1) <= MAXLEN_CMD)
    {
      mxGetString(prhs[0], cmd, cmdlen);
    } else {
      *cmd = (char) 0;
    }
  if ( strcmp ("reset",cmd) == 0 )
    {
      reset_options (rand_seed, test_in_cost_func,
		     use_rejected_cost, USER_ASA_OPTIONS);
      return;
    }
  if ( strcmp ("set",cmd) == 0 )
    {
      nlhs = 0; /* The set command returns nothing. */
      if (nrhs > 3)
	{
	  mexErrMsgTxt ("Error: Too many operands are given.");
	}
      if (nrhs >= 2)
	{
	  if (!mxIsChar (prhs[1]))
	    {
	      mexErrMsgTxt ("Error: The first operand must be an option name.");
	    }
	  else
	    {
	      if ((cmdlen = mxGetN(prhs[1]) + 1) <= MAXLEN_CMD)