ÎÒÓÐ simulated annealing with metropolies(Monte Carlo)×öµÄÒ»¸öÏîÄ¿µÄ´úÂ룬ÄãÒª¿´¿´Ã´£¿void anneal(int nparam, int nstep, int nstep_per_block, double t0,const double * param_in,double cost_in, double * params_out, double * cost_out) {int nblock;int step;int block;int nactive;int rank;int n_accepted = 0;int i, j, n;double cost_current, cost_trial;int * param_index;double * param_current;double * param_trial;double * Q;double * S;double * u;double * dp;double * A;FILE * fp_log_file;char fname[FILENAME_MAX];double temp = t0;double tempmax = temp;double ebar, evar, emin, eta, specific_heat;double delta;double chi = 0.8; // Annealing scheduledouble chi_s = 3.0; // Vanderbilt/Louie 'growth factor' double rm;double root3 = sqrt(3.0);double p = 0.02/sqrt(3.0); //max size of annealing stepparam_current = new double[nparam];param_trial = new double[nparam];cost_current = cost_in;MPI_Comm_rank(MPI_COMM_WORLD, &rank);sprintf(fname, "a_%4.4d.log", rank);fp_log_file = fopen(fname, "a");if (fp_log_file == (FILE *) NULL) errorMessage("fopen(log) failed\n");// Work out the number of active parameters, and set up the// index table of the active parameters.// Note that the complete array of parameters (param_trial) must // be used to evaluate the cost function.nactive = 0;for (n = 0; n < nparam; n++) {param_current[n] = param_in[n];param_trial[n] = param_in[n];if (P.is_active[n]) nactive++;}param_index = new int[nactive];i = 0;for (n = 0; n < nparam; n++) {if (P.is_active[n]) param_index[i++] = n;}// Initialise the step distribution matrix Q_ijQ = new double[nactive*nactive];S = new double[nactive*nactive];u = new double[nactive];dp = new double[nactive];A = new double[nactive];double * Qtmp;Qtmp = new double[nactive*nactive];for (i = 0; i < nactive; i++) {for (j = 0; j < nactive; j++) {delta = (i == j);Q[i*nactive + j] = p*delta*param_current[param_index[j]]; }}// carry out annealing pointsnblock = nstep/nstep_per_block;rm = 1.0/(double) nstep_per_block;for (block = 0; block < nblock; block++) {// Set the schedule for this block, and initialise blockwise quantities. // We also ensure the step distribution matrix is diagonal.temp = chi*temp;for (i = 0; i < nactive; i++) {A[i] = 0.0;for (j = 0; j < nactive; j++) {S[i*nactive + j] = 0.0;delta = (i == j);Q[i*nactive + j] *= delta;}}ebar = 0.0;evar = 0.0;emin = cost_current;for (i = 0; i < nactive; i++) {printf("Step: %d %g\n", i, Q[i*nactive + i]);}for (step = 0; step < nstep_per_block; step++) {// Set the random vector u, and compute the step size dpfor (i = 0; i < nactive; i++) {u[i] = root3*(r_uniform()*2.0 - 1.0);}for (i = 0; i < nactive; i++) {dp[i] = 0.0;for (j = 0; j < nactive; j++) {dp[i] += Q[i*nactive + j]*u[j];}}for (i = 0; i < nactive; i++) {n = param_index[i];param_trial[n] = param_current[n] + dp[i];if (param_trial[n] < P.min[n]) param_trial[n] = P.min[n];if (param_trial[n] > P.max[n]) param_trial[n] = P.max[n];}// calculate new cost function scorep_model->setParameters(param_trial);cost_trial = p_costWild->getCost();cost_trial += p_costLHY->getCost();cost_trial += p_costTOC1->getCost();cost_trial += p_costAPRR->getCost();// Metropolisdelta = cost_trial - cost_current;if (delta < 0.0 || r_uniform() < exp(-delta/temp)) {for (n = 0; n < nparam; n++) {param_current[n] = param_trial[n];}cost_current = cost_trial;++n_accepted;}// 'Energy' statisticsebar += cost_current;evar += cost_current*cost_current;if (cost_current < emin) emin = cost_current;// Per time step logfprintf(fp_log_file, "%6d %6d %10.4f %10.4f %10.4f %10.4f\n", block, step, temp,cost_current, cost_trial,(float) n_accepted / (float) (block*nstep_per_block + step)); // Accumulate average, measured covariancefor (i = 0; i < nactive; i++) {A[i] += param_current[param_index[i]];for (j = 0; j < nactive; j++) {S[i*nactive + j]+= param_current[param_index[i]]*param_current[param_index[j]]; }}/* Next step*/}// Set the previous block average and measured covariancefor (i = 0; i < nactive; i++) {A[i] = rm*A[i];}for (i = 0; i < nactive; i++) {for (j = 0; j < nactive; j++) {S[i*nactive + j] = rm*S[i*nactive + j] - A[i]*A[j];if (i == j) printf("Average: %d %g %g\n", i, A[i], S[i*nactive+j]); // Set the convarience for the next iteration s = 6 chi_s S / MS[i*nactive + j] = 6.0*chi_s*rm*S[i*nactive + j];}}// Reset the step distribution matrix for the next blocki = do_cholesky(nactive, S, Q);j = test_cholesky(nactive, S, Q);printf("Cholesky %d %d\n", i, j);// Block statisticsebar = rm*ebar;evar = rm*evar;specific_heat = (evar - ebar*ebar) / temp*temp;eta = (ebar - emin)/ebar;fprintf(fp_log_file, "%d %d %f %f %f %f %f %f\n", block, nstep_per_block, temp, ebar, evar, emin, specific_heat, eta);/* Next block */}*cost_out = cost_current;for (n = 0; n < nparam; n++) {params_out[n] = param_current[n];}fclose(fp_log_file);delete param_index;delete param_current; delete param_trial;delete Q;delete u;delete dp;delete S;delete A;return;}。