ÎÒÓÐ 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 schedule double 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 step param_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_ij Q = 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 points nblock = 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 dp for (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 score p_model->setParameters(param_trial); cost_trial = p_costWild->getCost(); cost_trial += p_costLHY->getCost(); cost_trial += p_costTOC1->getCost(); cost_trial += p_costAPRR->getCost(); // Metropolis delta = 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' statistics ebar += cost_current; evar += cost_current*cost_current; if (cost_current < emin) emin = cost_current; // Per time step log fprintf(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 covariance for (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 covariance for (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 / M S[i*nactive + j] = 6.0*chi_s*rm*S[i*nactive + j]; } } // Reset the step distribution matrix for the next block i = do_cholesky(nactive, S, Q); j = test_cholesky(nactive, S, Q); printf("Cholesky %d %d\n", i, j); // Block statistics ebar = 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;。