模拟退火算法解决路径优化 的源代码.doc

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 fnameFILENAME_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 doublenparam; param_trial = new doublenparam; 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) failedn); / 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_currentn = param_inn; param_trialn = param_inn; if (P.is_activen) nactive+; param_index = new intnactive; i = 0; for (n = 0; n nparam; n+) if (P.is_activen) param_indexi+ = n; / Initialise the step distribution matrix Q_ij Q = new doublenactive*nactive; S = new doublenactive*nactive; u = new doublenactive; dp = new doublenactive; A = new doublenactive; double * Qtmp; Qtmp = new doublenactive*nactive; for (i = 0; i nactive; i+) for (j = 0; j nactive; j+) delta = (i = j); Qi*nactive + j = p*delta*param_currentparam_indexj; / 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+) Ai = 0.0; for (j = 0; j nactive; j+) Si*nactive + j = 0.0; delta = (i = j); Qi*nactive + j *= delta; ebar = 0.0; evar = 0.0; emin = cost_current; for (i = 0; i nactive; i+) printf(Step: %d %gn, i, Qi*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+) ui = root3*(r_uniform()*2.0 - 1.0); for (i = 0; i nactive; i+) dpi = 0.0; for (j = 0; j nactive; j+) dpi += Qi*nactive + j*uj; for (i = 0; i nactive; i+) n = param_indexi; param_trialn = param_currentn + dpi; if (param_trialn P.maxn) param_trialn = P.maxn; / 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_currentn = param_trialn; 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.4fn, 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+) Ai += param_currentparam_indexi; for (j = 0; j nactive; j+) Si*nactive + j += param_currentparam_indexi*param_currentparam_indexj; /* Next step*/ / Set the previous block average and measured covariance for (i = 0; i nactive; i+) Ai = rm*Ai; for (i = 0; i nactive; i+) for (j = 0; j nactive; j+) Si*nactive + j = rm*Si*nactive + j - Ai*Aj; if (i = j) printf(Average: %d %g %gn, i, Ai, Si*nactive+j); / Set the convarience for the next iteration s = 6 chi_s S / M Si*nactive + j = 6.0*chi_s*rm*Si*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 %dn, 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 %fn, block, nstep_per_block, temp, ebar, evar, emin, specific_heat, eta); /* Next