asamin.c

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

		{
		  mxGetString(prhs[1], cmd, cmdlen);
		} else {
		  *cmd = (char) 0;
		}
	    }
	}
      if ((nrhs == 1) || (strcmp ("rand_seed",cmd) == 0))
	{
	  option_matched = TRUE;
	  if (nrhs > 2)
	    {
	      set_longint_to_option (prhs[2], rand_seed);
	    }
	  else
	    {
	      mexPrintf ("  rand_seed = %ld\n",*rand_seed);
	    }
	}
      if ((nrhs == 1) || (strcmp ("test_in_cost_func",cmd) == 0))
	{
	  option_matched = TRUE;
	  if (nrhs > 2)
	    {
	      set_int_to_option (prhs[2], test_in_cost_func);
	    }
	  else
	    {
	      mexPrintf ("  test_in_cost_func = %ld\n",*test_in_cost_func);
	    }
	}
      if ((nrhs == 1) || (strcmp ("use_rejected_cost",cmd) == 0))
	{
	  option_matched = TRUE;
	  if (nrhs > 2)
	    {
	      set_int_to_option (prhs[2], use_rejected_cost);
	    }
	  else
	    {
	      mexPrintf ("  use_rejected_cost = %ld\n",*use_rejected_cost);
	    }
	}
      if ((nrhs == 1) || (strcmp ("asa_out_file",cmd) == 0))
	{
	  option_matched = TRUE;
	  if (nrhs > 2)
	    {
	      set_string_to_option (prhs[2], USER_ASA_OPTIONS->Asa_Out_File);
	    }
	  else
	    {
	      mexPrintf ("  asa_out_file = '%s'\n",USER_ASA_OPTIONS->Asa_Out_File);
	    }
	}
      if ((nrhs == 1) || (strcmp ("limit_acceptances",cmd) == 0))
	{
	  option_matched = TRUE;
	  if (nrhs > 2)
	    {
	      set_longint_to_option (prhs[2], &USER_ASA_OPTIONS->Limit_Acceptances);
	    }
	  else
	    {
	      mexPrintf ("  limit_acceptances = %ld\n",
			 USER_ASA_OPTIONS->Limit_Acceptances);
	    }
	}
      if ((nrhs == 1) || (strcmp ("limit_generated",cmd) == 0))
	{
	  option_matched = TRUE;
	  if (nrhs > 2)
	    {
	      set_longint_to_option (prhs[2], &USER_ASA_OPTIONS->Limit_Generated);
	    }
	  else
	    {
	      mexPrintf ("  limit_generated = %ld\n",
			 USER_ASA_OPTIONS->Limit_Generated);
	    }
	}
      if ((nrhs == 1) || (strcmp ("limit_invalid",cmd) == 0))
	{
	  option_matched = TRUE;
	  if (nrhs > 2)
	    {
	      set_int_to_option (prhs[2],
				 &USER_ASA_OPTIONS->Limit_Invalid_Generated_States);
	    }
	  else
	    {
	      mexPrintf ("  limit_invalid = %d\n",
			 USER_ASA_OPTIONS->Limit_Invalid_Generated_States);
	    }
	}
      if ((nrhs == 1) || (strcmp ("accepted_to_generated_ratio",cmd) == 0))
	{
	  option_matched = TRUE;
	  if (nrhs > 2)
	    {
	      set_real_to_option (prhs[2],
				  &USER_ASA_OPTIONS->Accepted_To_Generated_Ratio);
	    }
	  else
	    {

	      mexPrintf ("  accepted_to_generated_ratio = %g\n",
			 USER_ASA_OPTIONS->Accepted_To_Generated_Ratio);
	    }
	}
      if ((nrhs == 1) ||  (strcmp ("cost_precision",cmd) == 0 ))
	{
	  option_matched = TRUE;
	  if (nrhs > 2)
	    {
	      set_real_to_option (prhs[2], &USER_ASA_OPTIONS->Cost_Precision);
	    }
	  else
	    {
	      mexPrintf ("  cost_precision = %g\n",
			 USER_ASA_OPTIONS->Cost_Precision);
	    }
	}
      if ((nrhs == 1) ||  (strcmp ("maximum_cost_repeat",cmd) == 0 ))
	{
	  option_matched = TRUE;
	  if (nrhs > 2)
	    {
	      set_int_to_option (prhs[2], &USER_ASA_OPTIONS->Maximum_Cost_Repeat);
	    }
	  else
	    {
	      mexPrintf ("  maximum_cost_repeat = %d\n",
			 USER_ASA_OPTIONS->Maximum_Cost_Repeat);
	    }
	}
      if ((nrhs == 1) ||  (strcmp ("number_cost_samples",cmd) == 0 ))
	{
	  option_matched = TRUE;
	  if (nrhs > 2)
	    {
	      set_int_to_option (prhs[2], &USER_ASA_OPTIONS->Number_Cost_Samples);
	    }
	  else
	    {
	      mexPrintf ("  number_cost_samples = %d\n",
			 USER_ASA_OPTIONS->Number_Cost_Samples);
	    }
	}
      if ((nrhs == 1) ||  (strcmp ("temperature_ratio_scale",cmd) == 0 ))
	{
	  option_matched = TRUE;
	  if (nrhs > 2)
	    {
	      set_real_to_option (prhs[2],
				  &USER_ASA_OPTIONS->Temperature_Ratio_Scale);
	    }
	  else
	    {
	      mexPrintf ("  temperature_ratio_scale = %g\n",
			 USER_ASA_OPTIONS->Temperature_Ratio_Scale);
	    }
	}
      if ((nrhs == 1) ||  (strcmp ("cost_parameter_scale_ratio",cmd) == 0 ))
	{
	  option_matched = TRUE;
	  if (nrhs > 2)
	    {
	      set_real_to_option (prhs[2],
				  &USER_ASA_OPTIONS->Cost_Parameter_Scale_Ratio);
	    }
	  else
	    {
	      mexPrintf ("  cost_parameter_scale_ratio = %g\n",
			 USER_ASA_OPTIONS->Cost_Parameter_Scale_Ratio);
	    }
	}
      if ((nrhs == 1) ||  (strcmp ("temperature_anneal_scale",cmd) == 0 ))
	{
	  option_matched = TRUE;
	  if (nrhs > 2)
	    {
	      set_real_to_option (prhs[2],
				  &USER_ASA_OPTIONS->Temperature_Anneal_Scale);
	    }
	  else
	    {
	      mexPrintf ("  temperature_anneal_scale = %g\n",
			 USER_ASA_OPTIONS->Temperature_Anneal_Scale);
	    }
	}
      if ((nrhs == 1) ||  (strcmp ("include_integer_parameters",cmd) == 0 ))
	{
	  option_matched = TRUE;
	  if (nrhs > 2)
	    {
	      set_int_to_option (prhs[2],
				 &USER_ASA_OPTIONS->Include_Integer_Parameters);
	    }
	  else
	    {
	      mexPrintf ("  include_integer_parameters = %d\n",
			 USER_ASA_OPTIONS->Include_Integer_Parameters);
	    }
	}
      if ((nrhs == 1) ||  (strcmp ("user_initial_parameters",cmd) == 0 ))
	{
	  option_matched = TRUE;
	  if (nrhs > 2)
	    {
	      set_int_to_option (prhs[2],
				 &USER_ASA_OPTIONS->User_Initial_Parameters);
	    }
	  else
	    {
	      mexPrintf ("  user_initial_parameters = %d\n",
			 USER_ASA_OPTIONS->User_Initial_Parameters);
	    }
	}
      if ((nrhs == 1) ||  (strcmp ("sequential_parameters",cmd) == 0 ))
	{
	  option_matched = TRUE;
	  if (nrhs > 2)
	    {
	      set_longint_to_option (prhs[2],
				     &USER_ASA_OPTIONS->Sequential_Parameters);
	    }
	  else
	    {
	      mexPrintf ("  sequential_parameters = %d\n",
			 USER_ASA_OPTIONS->Sequential_Parameters);
	    }
	}
      if ((nrhs == 1) ||  (strcmp ("initial_parameter_temperature",cmd) == 0 ))
	{
	  option_matched = TRUE;
	  if (nrhs > 2)
	    {
	      set_real_to_option (prhs[2],
				  &USER_ASA_OPTIONS->Initial_Parameter_Temperature);
	    }
	  else
	    {
	      mexPrintf ("  initial_parameter_temperature = %g\n",
			 USER_ASA_OPTIONS->Initial_Parameter_Temperature);
	    }
	}
      if ((nrhs == 1) ||  (strcmp ("acceptance_frequency_modulus",cmd) == 0 ))
	{
	  option_matched = TRUE;
	  if (nrhs > 2)
	    {
	      set_int_to_option (prhs[2],
				 &USER_ASA_OPTIONS->Acceptance_Frequency_Modulus);
	    }
	  else
	    {
	      mexPrintf ("  acceptance_frequency_modulus = %d\n",
			 USER_ASA_OPTIONS->Acceptance_Frequency_Modulus);
	    }
	}
      if ((nrhs == 1) ||  (strcmp ("generated_frequency_modulus",cmd) == 0 ))
	{
	  option_matched = TRUE;
	  if (nrhs > 2)
	    {
	      set_int_to_option (prhs[2],
				 &USER_ASA_OPTIONS->Generated_Frequency_Modulus);
	    }
	  else
	    {
	      mexPrintf ("  generated_frequency_modulus = %d\n",
			 USER_ASA_OPTIONS->Generated_Frequency_Modulus);
	    }
	}
      if ((nrhs == 1) ||  (strcmp ("reanneal_cost",cmd) == 0 ))
	{
	  option_matched = TRUE;
	  if (nrhs > 2)
	    {
	      set_int_to_option (prhs[2],
				 &USER_ASA_OPTIONS->Reanneal_Cost);
	    }
	  else
	    {
	      mexPrintf ("  reanneal_cost = %d\n",
			 USER_ASA_OPTIONS->Reanneal_Cost);
	    }
	}
      if ((nrhs == 1) ||  (strcmp ("reanneal_parameters",cmd) == 0 ))
	{
	  option_matched = TRUE;
	  if (nrhs > 2)
	    {
	      set_int_to_option (prhs[2],
				 &USER_ASA_OPTIONS->Reanneal_Parameters);
	    }
	  else
	    {
	      mexPrintf ("  reanneal_parameters = %d\n",
			 USER_ASA_OPTIONS->Reanneal_Parameters);
	    }
	}
      if ((nrhs == 1) ||  (strcmp ("delta_x",cmd) == 0 ))
	{
	  option_matched = TRUE;
	  if (nrhs > 2)
	    {
	      set_real_to_option (prhs[2],
				  &USER_ASA_OPTIONS->Delta_X);
	    }
	  else
	    {
	      mexPrintf ("  delta_x = %g\n",
			 USER_ASA_OPTIONS->Delta_X);
	    }
	}
      if (!option_matched) mexErrMsgTxt ("Error: The option name is invalid.");
      return;
    }
  if ( strcmp ("minimize",cmd) == 0 )
    {
      nlhs = 3;  /* The minimize command returns three matrices. */
      if ( nrhs < 6 )
	{
	  mexErrMsgTxt ("Error: At least 5 operands are expected for 'minimize'.");
	}
      
      if (!mxIsChar (prhs[1]))
	{
	  mexErrMsgTxt ("Error: The first operand must be the name of the cost function");
	}
      if ((cmdlen = mxGetN(prhs[1]) + 1) > MAXLEN_COST_FUNC_NAME)
	{
	  mexErrMsgTxt ("Error: The name of the cost function is too long.");
	}
      mxGetString(prhs[1], cost_func_name, cmdlen);
      
      if ((mxGetN (prhs[2]) != 1) || (mxGetN (prhs[3]) != 1) ||
	  (mxGetN (prhs[4]) != 1) || (mxGetN (prhs[5]) != 1) ||
	  (mxGetM (prhs[3]) != mxGetM (prhs[2]) ) ||
	  (mxGetM (prhs[4]) != mxGetM (prhs[2]) ) ||
	  (mxGetM (prhs[5]) != mxGetM (prhs[2]) ) ||
	  (!mxIsDouble (prhs[2])) || (mxIsComplex (prhs[2])) ||
	  (!mxIsDouble (prhs[3])) || (mxIsComplex (prhs[3])) ||
	  (!mxIsDouble (prhs[4])) || (mxIsComplex (prhs[4])) ||
	  (!mxIsDouble (prhs[5])) || (mxIsComplex (prhs[5])))
	{
	  mexErrMsgTxt ("Error: The second through fifth operands must be columns vectors of the same size.");
	}
      plhs[0] = mxCreateDoubleMatrix (1, 1, mxREAL);
      minimum_cost_value = mxGetPr (plhs[0]);

      if ((parameter_dimension =
	   (ALLOC_INT *) mxCalloc (1, sizeof (ALLOC_INT))) == NULL)
	{
	  mexErrMsgTxt ("Internal Error: parameter_dimension cannot be allocated.");
	}
      *parameter_dimension = mxGetM(prhs[2]);

      /* Allocate and initialize parameter initial values;
	 the parameter final values will be stored here later */

      plhs[1] = mxCreateDoubleMatrix (*parameter_dimension, 1, mxREAL);
      cost_parameters = mxGetPr (plhs[1]);
      memcpy (cost_parameters, mxGetPr (prhs[2]),
	      *parameter_dimension * sizeof (double));
	  
      /* Allocate and initialize parameter lower bounds */
      parameter_lower_bound = mxGetPr (prhs[3]);
      
      /* Allocate and initialize parameter upper bounds */
      parameter_upper_bound = mxGetPr (prhs[4]);

      /* Allocate and initialize the parameter types, real or integer */
      if ((parameter_int_real =
	   (int *) mxCalloc (*parameter_dimension, sizeof (int))
	   ) == NULL)
	{
	  mexErrMsgTxt ("Internal Error: parameter_int_real cannot be allocated.");
	}
      ptr = mxGetPr (prhs[5]);
      for (i=0 ; i < *parameter_dimension ; i++) {
	parameter_int_real[i] = *ptr;
	ptr++;
      }

      /* Allocate space for parameter cost_tangents -
	 used for reannealing */
      plhs[2] =	mxCreateDoubleMatrix (*parameter_dimension, 1, mxREAL);
      cost_tangents = mxGetPr (plhs[2]);

      /* Allocate space for parameter cost_curvatures/covariance */
      plhs[3] =
	mxCreateDoubleMatrix (*parameter_dimension,
			      *parameter_dimension,
			      mxREAL);
      cost_curvature = mxGetPr (plhs[3]);

      /* Allocate space for status (exit_code and cost_flag),
       which will be set after returning from asa()*/
      
      plhs[4] =
	mxCreateDoubleMatrix (2, 1, mxREAL);

      /* Allocate the space for the RHS variables of the matlab cost
	 functions (global pointer variables)*/

      if (*test_in_cost_func == 0)
	{
	  cost_function = cost_function_without_test;
	  matlab_cost_func_nrhs = nrhs - 5;
	}
      else
	{
	  cost_function = cost_function_with_test;
	  matlab_cost_func_nrhs = nrhs - 3;
	}
      if ((matlab_cost_func_prhs =
	   mxCalloc (matlab_cost_func_nrhs, sizeof (mxArray *)))
	  == NULL)
	{
	  mexErrMsgTxt ("Internal Error: matlab_cost_func_prhs cannot be allocated.");
	}

      matlab_cost_func_prhs[0] =
	mxCreateDoubleMatrix (*parameter_dimension, 1, mxREAL);
      if (*test_in_cost_func == 0) {
	for (i = 0; i < (nrhs - 6) ; i++) {
	  matlab_cost_func_prhs[1 + i] = prhs[6 + i];
	}
	if ((critical_cost_value =
	     (double *) mxCalloc (1, sizeof (double)))
	    == NULL)
	  {
	    mexErrMsgTxt ("Internal Error: critical_cost_value cannot be allocated.");
	  }
	if ((user_acceptance_flag =
	     (double *) mxCalloc (1, sizeof (double)))
	    == NULL)
	  {
	    mexErrMsgTxt ("Internal Error: user_acceptance_flag cannot be allocated.");
	  }
      } else {
	matlab_cost_func_prhs[1] = mxCreateDoubleMatrix (1, 1, mxREAL);
	critical_cost_value = mxGetPr (matlab_cost_func_prhs[1]);
	matlab_cost_func_prhs[2] = mxCreateDoubleMatrix (1, 1, mxREAL);
	user_acceptance_flag = mxGetPr(matlab_cost_func_prhs[2]);
	for (i = 0; i < (nrhs - 6) ; i++) {
	  matlab_cost_func_prhs[3 + i] = prhs[6 + i];
	}
      }

      resettable_randflt (rand_seed, 1); /* initialize random
					    number generator */

      *minimum_cost_value =
	asa (cost_function,
	     randflt,
	     rand_seed,
	     cost_parameters,
	     parameter_lower_bound,
	     parameter_upper_bound,
	     cost_tangents,
	     cost_curvature,
	     parameter_dimension,
	     parameter_int_real,
	     cost_flag,
	     exit_code,
	     USER_ASA_OPTIONS);

      *(mxGetPr (plhs[4])) = *exit_code;
      *(mxGetPr (plhs[4]) + 1) = *cost_flag;
	
      if (*exit_code == -1)
	{
	  mexErrMsgTxt ("Error: calloc in asa () failed.");
	}
      
      /* Free local memories */
      mxDestroyArray (matlab_cost_func_prhs[0]);
      if (*test_in_cost_func != 0) {
	mxDestroyArray (matlab_cost_func_prhs[1]);
	mxDestroyArray (matlab_cost_func_prhs[2]);
      }
      if (*test_in_cost_func == 0) {
	mxFree (critical_cost_value);
	mxFree (user_acceptance_flag);
      }
      mxFree (matlab_cost_func_prhs);
      mxFree (parameter_int_real);
      mxFree (parameter_dimension);
      return;
    }
  mexErrMsgTxt ("Error: No such a command.");
}
  
static void
exit_function ()
{
  mxFree (USER_ASA_OPTIONS->Asa_Out_File);
  mxFree (USER_ASA_OPTIONS);
}