搜索
ÎÒÓÐsimulatedannealingwithmetropolies(MonteCarlo)×öµÄÒ»¸öÏîÄ¿µÄ´úÂ룬ÄãÒª¿´¿´Ã´£¿voidanneal(intnparam,intnstep,intnstep_per_block,doublet0,constdouble*param_in,doublecost_in,double*params_out,double*cost_out){intnblock;intstep;intblock;intnactive;intrank;intn_accepted=0;inti,j,n;doublecost_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;charfname[FILENAME_MAX];doubletemp=t0;doubletempmax=temp;doubleebar,evar,emin,eta,specific_heat;doubledelta;doublechi=0.8;//Annealingscheduledoublechi_s=3.0;//Vanderbilt/Louie'growthfactor'doublerm;doubleroot3=sqrt(3.0);doublep=0.02/sqrt(3.0);//maxsizeofannealingstepparam_current=newdouble[nparam];param_trial=newdouble[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);//Workoutthenumberofactiveparameters,andsetupthe//indextableoftheactiveparameters.//Notethatthecompletearrayofparameters(param_trial)must//beusedtoevaluatethecostfunction.nactive=0;for(n=0;nnparam;n++){param_current[n]=param_in[n];param_trial[n]=param_in[n];if(P.is_active[n])nactive++;}param_index=newint[nactive];i=0;for(n=0;nnparam;n++){if(P.is_active[n])param_index[i++]=n;}//InitialisethestepdistributionmatrixQ_ijQ=newdouble[nactive*nactive];S=newdouble[nactive*nactive];u=newdouble[nactive];dp=newdouble[nactive];A=newdouble[nactive];double*Qtmp;Qtmp=newdouble[nactive*nactive];for(i=0;inactive;i++){for(j=0;jnactive;j++){delta=(i==j);Q[i*nactive+j]=p*delta*param_current[param_index[j]];}}//carryoutannealingpointsnblock=nstep/nstep_per_block;rm=1.0/(double)nstep_per_block;for(block=0;blocknblock;block++){//Setthescheduleforthisblock,andinitialiseblockwisequantities.//Wealsoensurethestepdistributionmatrixisdiagonal.temp=chi*temp;for(i=0;inactive;i++){A[i]=0.0;for(j=0;jnactive;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;inactive;i++){printf(Step:%d%g\n,i,Q[i*nactive+i]);}for(step=0;stepnstep_per_block;step++){//Settherandomvectoru,andcomputethestepsizedpfor(i=0;inactive;i++){u[i]=root3*(r_uniform()*2.0-1.0);}for(i=0;inactive;i++){dp[i]=0.0;for(j=0;jnactive;j++){dp[i]+=Q[i*nactive+j]*u[j];}}for(i=0;inactive;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];}//calculatenewcostfunctionscorep_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(delta0.0||r_uniform()exp(-delta/temp)){for(n=0;nnparam;n++){param_current[n]=param_trial[n];}cost_current=cost_trial;++n_accepted;}//'Energy'statisticsebar+=cost_current;evar+=cost_current*cost_current;if(cost_currentemin)emin=cost_current;//Pertimesteplogfprintf(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));//Accumulateaverage,measuredcovariancefor(i=0;inactive;i++){A[i]+=param_current[param_index[i]];for(j=0;jnactive;j++){S[i*nactive+j]+=param_current[param_index[i]]*param_current[param_index[j]];}}/*Nextstep*/}//Setthepreviousblockaverageandmeasuredcovariancefor(i=0;inactive;i++){A[i]=rm*A[i];}for(i=0;inactive;i++){for(j=0;jnactive;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]);//Settheconvarienceforthenextiterations=6chi_sS/MS[i*nactive+j]=6.0*chi_s*rm*S[i*nactive+j];}}//Resetthestepdistributionmatrixforthenextblocki=do_cholesky(nactive,S,Q);j=test_cholesky(nactive,S,Q);printf(Cholesky%d%d\n,i,j);//Blockstatisticsebar=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);/*Nextblock*/}*cost_out=cost_current;for(n=0;nnparam;n++){params_out[n]=param_current[n];}fclose(fp_log_file);deleteparam_index;deleteparam_current;deleteparam_trial;deleteQ;deleteu;deletedp;deleteS;deleteA;return;}