Logo Search packages:      
Sourcecode: objcryst-fox version File versions  Download package

void ObjCryst::MonteCarloObj::RunParallelTempering ( long &  nbSteps,
const bool  silent = false,
const REAL  finalcost = 0,
const REAL  maxTime = -1 
)

For internal use only.

Do a single Parallel Tempering run. This is called by Optimize(...) and MultiRunOptimize(), which must also prepare the optimization (PrepareRefParList(), etc..).

Definition at line 881 of file GlobalOptimObj.cpp.

References ObjCryst::RefinableObj::ClearParamSet(), ObjCryst::RefinableObj::CreateParamSet(), ObjCryst::OptimizationObj::DisplayReport(), ObjCryst::RefObjOpt::GetChoice(), ObjCryst::OptimizationObj::GetLogLikelihood(), ObjCryst::OptimizationObj::GetName(), ObjCryst::RefinableObj::GetNbPar(), ObjCryst::RefinableObj::GetPar(), ObjCryst::RefinableObj::GetParamSet(), mAnnealingScheduleMutation, mAnnealingScheduleTemp, ObjCryst::OptimizationObj::mBestCost, ObjCryst::OptimizationObj::mBestParSavedSetIndex, ObjCryst::OptimizationObj::mContext, mCurrentCost, ObjCryst::OptimizationObj::mLastOptimTime, mMutationAmplitude, mMutationAmplitudeGamma, mMutationAmplitudeMax, mMutationAmplitudeMin, ObjCryst::OptimizationObj::mNbTrial, ObjCryst::OptimizationObj::mRecursiveRefinedObjList, ObjCryst::OptimizationObj::mRefParList, ObjCryst::OptimizationObj::mStopAfterCycle, mTemperature, mTemperatureGamma, mTemperatureMax, mTemperatureMin, ObjCryst::OptimizationObj::mvContextObjStats, ObjCryst::OptimizationObj::mvObjWeight, ObjCryst::OptimizationObj::mXMLAutoSave, NewConfiguration(), Chronometer::print(), ObjCryst::RefinablePar::Print(), ObjCryst::RefinableObj::RestoreParamSet(), ObjCryst::RefinableObj::SaveParamSet(), Chronometer::seconds(), Chronometer::start(), ObjCryst::OptimizationObj::TagNewBestConfig(), ObjCryst::OptimizationObj::UpdateDisplay(), and ObjCryst::XMLCrystFileSaveGlobal().

Referenced by MultiRunOptimize(), and Optimize().

{
   //Keep a copy of the total number of steps, and decrement nbStep
   const long nbSteps=nbStep;
   unsigned int accept;// 1 if last trial was accepted? 2 if new best config ? else 0
   mNbTrial=0;
   // time (in seconds) when last autoSave was made (if enabled)
      unsigned long secondsWhenAutoSave=0;

   if(!silent) cout << "Starting Parallel Tempering Optimization"<<endl;
   //Total number of parallel refinements,each is a 'World'. The most stable
   // world must be i=nbWorld-1, and the most changing World (high mutation,
   // high temperature) is i=0.
      const long nbWorld=30;
      CrystVector_long worldSwapIndex(nbWorld);
      for(int i=0;i<nbWorld;++i) worldSwapIndex(i)=i;
   // Number of successive trials for each World. At the end of these trials
   // a swap is tried with the upper World (eg i-1). This number effectvely sets
   // the rate of swapping.
      const int nbTryPerWorld=10;
   // Initialize the costs
      mCurrentCost=this->GetLogLikelihood();
      REAL runBestCost=mCurrentCost;
      CrystVector_REAL currentCost(nbWorld);
      currentCost=mCurrentCost;
   // Init the different temperatures
      CrystVector_REAL simAnnealTemp(nbWorld);
      for(int i=0;i<nbWorld;i++)
      {
         switch(mAnnealingScheduleTemp.GetChoice())
         {
            case ANNEALING_BOLTZMANN:
               simAnnealTemp(i)=
                  mTemperatureMin*log((REAL)nbWorld)/log((REAL)(i+2));break;
            case ANNEALING_CAUCHY:
               simAnnealTemp(i)=mTemperatureMin*nbWorld/(i+1);break;
            //case ANNEALING_QUENCHING:
            case ANNEALING_EXPONENTIAL:
               simAnnealTemp(i)=mTemperatureMax
                              *pow(mTemperatureMin/mTemperatureMax,
                                    i/(REAL)(nbWorld-1));break;
            case ANNEALING_GAMMA:
               simAnnealTemp(i)=mTemperatureMax+(mTemperatureMin-mTemperatureMax)
                              *pow(i/(REAL)(nbWorld-1),mTemperatureGamma);break;
            case ANNEALING_SMART:
               simAnnealTemp(i)=mCurrentCost/(100.+(REAL)i/(REAL)nbWorld*900.);break;
            default:
               simAnnealTemp(i)=mCurrentCost/(100.+(REAL)i/(REAL)nbWorld*900.);break;
         }
      }
   //Init the different mutation rate parameters
      CrystVector_REAL mutationAmplitude(nbWorld);
      for(int i=0;i<nbWorld;i++)
      {
         switch(mAnnealingScheduleMutation.GetChoice())
         {
            case ANNEALING_BOLTZMANN:
               mutationAmplitude(i)=
                  mMutationAmplitudeMin*log((REAL)(nbWorld-1))/log((REAL)(i+2));
               break;
            case ANNEALING_CAUCHY:
               mutationAmplitude(i)=mMutationAmplitudeMin*(REAL)(nbWorld-1)/(i+1);break;
            //case ANNEALING_QUENCHING:
            case ANNEALING_EXPONENTIAL:
               mutationAmplitude(i)=mMutationAmplitudeMax
                              *pow(mMutationAmplitudeMin/mMutationAmplitudeMax,
                                    i/(REAL)(nbWorld-1));break;
            case ANNEALING_GAMMA:
               mutationAmplitude(i)=mMutationAmplitudeMax+(mMutationAmplitudeMin-mMutationAmplitudeMax)
                              *pow(i/(REAL)(nbWorld-1),mMutationAmplitudeGamma);break;
            case ANNEALING_SMART:
               mutationAmplitude(i)=sqrt(mMutationAmplitudeMin*mMutationAmplitudeMax);break;
            default:
               mutationAmplitude(i)=sqrt(mMutationAmplitudeMin*mMutationAmplitudeMax);break;
         }
      }
   // Init the parameter sets for each World
   // All Worlds start from the same (current) configuration.
      CrystVector_long worldCurrentSetIndex(nbWorld);
      for(int i=nbWorld-1;i>=0;i--)
      {
         if((i!=(nbWorld-1))&&(i%2==0))
            for(int j=0;j<mRecursiveRefinedObjList.GetNb();j++)
               mRecursiveRefinedObjList.GetObj(j).RandomizeConfiguration();
         worldCurrentSetIndex(i)=mRefParList.CreateParamSet();
         mRefParList.RestoreParamSet(worldCurrentSetIndex(nbWorld-1));
      }
      //mNbTrial=nbSteps;;
      const long lastParSavedSetIndex=mRefParList.CreateParamSet("MonteCarloObj:Last parameters (PT)");
      const long runBestIndex=mRefParList.CreateParamSet("Best parameters for current run (PT)");
      CrystVector_REAL swapPar;
   //Keep track of how many trials are accepted for each World
      CrystVector_long worldNbAcceptedMoves(nbWorld);
      worldNbAcceptedMoves=0;
   //Do a report each... And check if mutation rate is OK (for annealing_smart)s
      const int nbTrialsReport=3000;
   // TEST : allow GENETIC mating of configurations
      //Get gene groups list :TODO: check for missing groups
         CrystVector_uint refParGeneGroupIndex(mRefParList.GetNbPar());
         unsigned int first=1;
         for(int i=0;i<mRecursiveRefinedObjList.GetNb();i++) 
            mRecursiveRefinedObjList.GetObj(i).GetGeneGroup(mRefParList,refParGeneGroupIndex,first);
         #if 0
         if(!silent) 
            for(int i=0;i<mRefParList.GetNbPar();i++)
            {
               cout << "Gene Group:"<<refParGeneGroupIndex(i)<<" :";
               mRefParList.GetPar(i).Print();
            }
         #endif
      // number of gene groups
      // to select which gene groups are exchanged in the mating
         //const unsigned int nbGeneGroup=refParGeneGroupIndex.max();
         //CrystVector_int crossoverGroupIndex(nbGeneGroup);
         //const long parSetOffspringA=mRefParList.CreateParamSet("Offspring A");
         //const long parSetOffspringB=mRefParList.CreateParamSet("Offspring B");
   // record the statistical distribution n=f(cost function) for each World
      //CrystMatrix_REAL trialsDensity(100,nbWorld+1);
      //trialsDensity=0;
      //for(int i=0;i<100;i++) trialsDensity(i,0)=i/(float)100;
   // Do we need to update the display ?
   bool needUpdateDisplay=false;
   //Do the refinement
   bool makeReport=false;
   Chronometer chrono;
   chrono.start();
   float lastUpdateDisplayTime=chrono.seconds();
   for(;mNbTrial<nbSteps;)
   {
      for(int i=0;i<nbWorld;i++)
      {
         mContext=i;
         //mRefParList.RestoreParamSet(worldCurrentSetIndex(i));
         mMutationAmplitude=mutationAmplitude(i);
         mTemperature=simAnnealTemp(i);
         for(int j=0;j<nbTryPerWorld;j++)
         {
            //mRefParList.SaveParamSet(lastParSavedSetIndex);
            mRefParList.RestoreParamSet(worldCurrentSetIndex(i));
            this->NewConfiguration();
            accept=0;
            REAL cost=this->GetLogLikelihood();
            //trialsDensity((long)(cost*100.),i+1)+=1;
            if(cost<currentCost(i))
            {
               accept=1;
               currentCost(i)=cost;
               mRefParList.SaveParamSet(worldCurrentSetIndex(i));
               if(cost<runBestCost)
               {
                  accept=2;
                  runBestCost=currentCost(i);
                  this->TagNewBestConfig();
                  needUpdateDisplay=true;
                  
                  mRefParList.SaveParamSet(runBestIndex);
                  if(runBestCost<mBestCost)
                  {
                     mBestCost=currentCost(i);
                     mRefParList.SaveParamSet(mBestParSavedSetIndex);
                     if(!silent) cout << "->Trial :" << mNbTrial 
                                   << " World="<< worldSwapIndex(i)
                                   << " Temp="<< mTemperature
                                   << " Mutation Ampl.: "<<mMutationAmplitude
                                   << " NEW OVERALL Best Cost="<<mBestCost<< endl;
                  }
                  else if(!silent) cout << "->Trial :" << mNbTrial 
                                   << " World="<< worldSwapIndex(i)
                                   << " Temp="<< mTemperature
                                   << " Mutation Ampl.: "<<mMutationAmplitude
                                   << " NEW RUN Best Cost="<<runBestCost<< endl;
                  if(!silent) this->DisplayReport();
               }
               worldNbAcceptedMoves(i)++;
            }
            else
            {
               if(log((rand()+1)/(REAL)RAND_MAX)<(-(cost-currentCost(i))/mTemperature) )
               {
                  accept=1;
                  currentCost(i)=cost;
                  mRefParList.SaveParamSet(worldCurrentSetIndex(i));
                  worldNbAcceptedMoves(i)++;
               }
            }
            if(  ((mXMLAutoSave.GetChoice()==1)&&((chrono.seconds()-secondsWhenAutoSave)>86400))
               ||((mXMLAutoSave.GetChoice()==2)&&((chrono.seconds()-secondsWhenAutoSave)>3600))
               ||((mXMLAutoSave.GetChoice()==3)&&((chrono.seconds()-secondsWhenAutoSave)> 600))
               ||((mXMLAutoSave.GetChoice()==4)&&(accept==2)) )
            {
               secondsWhenAutoSave=(unsigned long)chrono.seconds();
               string saveFileName=this->GetName();
               time_t date=time(0);
               char strDate[40];
               strftime(strDate,sizeof(strDate),"%Y-%m-%d_%H-%M-%S",localtime(&date));//%Y-%m-%dT%H:%M:%S%Z
               char costAsChar[30];
               if(accept!=2) mRefParList.RestoreParamSet(mBestParSavedSetIndex);
               sprintf(costAsChar,"-Cost-%f",this->GetLogLikelihood());
               saveFileName=saveFileName+(string)strDate+(string)costAsChar+(string)".xml";
               XMLCrystFileSaveGlobal(saveFileName);
               //if(accept!=2) mRefParList.RestoreParamSet(lastParSavedSetIndex);
            }
            //if(accept==0) mRefParList.RestoreParamSet(lastParSavedSetIndex);
            mNbTrial++;nbStep--;
            if((mNbTrial%nbTrialsReport)==0) makeReport=true;
         }//nbTryPerWorld trials
      }//For each World
      //Try swapping worlds
      for(int i=1;i<nbWorld;i++)
      {
         #if 0
         mRefParList.RestoreParamSet(worldCurrentSetIndex(i));
         mMutationAmplitude=mutationAmplitude(i);
         cout<<i<<":"<<currentCost(i)<<":"<<this->GetLogLikelihood()<<endl;
         #endif
         #if 1
         if( log((rand()+1)/(REAL)RAND_MAX)
                < (-(currentCost(i-1)-currentCost(i))/simAnnealTemp(i)))
         #else
         // Compare World (i-1) and World (i) with the same amplitude,
         // hence the same max likelihood error
         mRefParList.RestoreParamSet(worldCurrentSetIndex(i-1));
         mMutationAmplitude=mutationAmplitude(i);
         if( log((rand()+1)/(REAL)RAND_MAX)
                < (-(this->GetLogLikelihood()-currentCost(i))/simAnnealTemp(i)))
         #endif
         {  
         /*
            if(i>2)
            {
               cout <<"->Swapping Worlds :" << i <<"(cost="<<currentCost(i)<<")"
                    <<" with "<< (i-1) <<"(cost="<< currentCost(i-1)<<")"<<endl;
            }
            */
            swapPar=mRefParList.GetParamSet(worldCurrentSetIndex(i));
            mRefParList.GetParamSet(worldCurrentSetIndex(i))=
               mRefParList.GetParamSet(worldCurrentSetIndex(i-1));
            mRefParList.GetParamSet(worldCurrentSetIndex(i-1))=swapPar;
            const REAL tmp=currentCost(i);
            currentCost(i)=currentCost(i-1);
            currentCost(i-1)=tmp;
            const long tmpIndex=worldSwapIndex(i);
            worldSwapIndex(i)=worldSwapIndex(i-1);
            worldSwapIndex(i-1)=tmpIndex;
            #if 0
            // Compute correct costs in the case we use maximum likelihood
            mRefParList.RestoreParamSet(worldCurrentSetIndex(i));
            mMutationAmplitude=mutationAmplitude(i);
            currentCost(i)=this->GetLogLikelihood();

            mRefParList.RestoreParamSet(worldCurrentSetIndex(i-1));
            mMutationAmplitude=mutationAmplitude(i-1);
            currentCost(i-1)=this->GetLogLikelihood();
            #endif
         }
      }
      #if 0
      //Try mating worlds- NEW !
      TAU_PROFILE_TIMER(timer1,\
               "MonteCarloObj::Optimize (Try mating Worlds)"\
               ,"", TAU_FIELD);
      TAU_PROFILE_START(timer1);
      if( (rand()/(REAL)RAND_MAX)<.1)
      for(int k=nbWorld-1;k>nbWorld/2;k--)
         for(int i=k-nbWorld/3;i<k;i++)
         {
            #if 0
            // Random switching of gene groups
            for(unsigned int j=0;j<nbGeneGroup;j++) 
               crossoverGroupIndex(j)= (int) floor(rand()/((REAL)RAND_MAX-1)*2);
            for(int j=0;j<mRefParList.GetNbPar();j++)
            {
               if(0==crossoverGroupIndex(refParGeneGroupIndex(j)-1))
               {
                  mRefParList.GetParamSet(parSetOffspringA)(j)=
                     mRefParList.GetParamSet(worldCurrentSetIndex(i))(j);
                  mRefParList.GetParamSet(parSetOffspringB)(j)=
                     mRefParList.GetParamSet(worldCurrentSetIndex(k))(j);
               }
               else
               {
                  mRefParList.GetParamSet(parSetOffspringA)(j)=
                     mRefParList.GetParamSet(worldCurrentSetIndex(k))(j);
                  mRefParList.GetParamSet(parSetOffspringB)(j)=
                     mRefParList.GetParamSet(worldCurrentSetIndex(i))(j);
               }
            }
            #endif
            #if 1
            // Switch gene groups in two parts
            unsigned int crossoverPoint1=
               (int)(1+floor(rand()/((REAL)RAND_MAX-1)*(nbGeneGroup)));
            unsigned int crossoverPoint2=
               (int)(1+floor(rand()/((REAL)RAND_MAX-1)*(nbGeneGroup)));
            if(crossoverPoint2<crossoverPoint1)
            {
               int tmp=crossoverPoint1;
               crossoverPoint1=crossoverPoint2;
               crossoverPoint2=tmp;
            }
            if(crossoverPoint1==crossoverPoint2) crossoverPoint2+=1;
            for(int j=0;j<mRefParList.GetNbPar();j++)
            {
               if((refParGeneGroupIndex(j)>crossoverPoint1)&&refParGeneGroupIndex(j)<crossoverPoint2)
               {
                  mRefParList.GetParamSet(parSetOffspringA)(j)=
                     mRefParList.GetParamSet(worldCurrentSetIndex(i))(j);
                  mRefParList.GetParamSet(parSetOffspringB)(j)=
                     mRefParList.GetParamSet(worldCurrentSetIndex(k))(j);
               }
               else
               {
                  mRefParList.GetParamSet(parSetOffspringA)(j)=
                     mRefParList.GetParamSet(worldCurrentSetIndex(k))(j);
                  mRefParList.GetParamSet(parSetOffspringB)(j)=
                     mRefParList.GetParamSet(worldCurrentSetIndex(i))(j);
               }
            }
            #endif
            // Try both offspring
            for(int junk=0;junk<2;junk++)
            {
               if(junk==0) mRefParList.RestoreParamSet(parSetOffspringA);
               else mRefParList.RestoreParamSet(parSetOffspringB);
               REAL cost=this->GetLogLikelihood();
               //if(log((rand()+1)/(REAL)RAND_MAX)
               //    < (-(cost-currentCost(k))/simAnnealTemp(k)))
               if(cost<currentCost(k))
               {
                  // Also exchange genes for higher-temperature World ?
                  //if(junk==0) 
                  //   mRefParList.GetParamSet(worldCurrentSetIndex(i))=
                  //      mRefParList.GetParamSet(parSetOffspringB);
                  //else
                  //   mRefParList.GetParamSet(worldCurrentSetIndex(i))=
                  //      mRefParList.GetParamSet(parSetOffspringA);
                  currentCost(k)=cost;
                  mRefParList.SaveParamSet(worldCurrentSetIndex(k));
                  //worldNbAcceptedMoves(k)++;
                  if(!silent) cout << "Accepted mating :"<<k<<"(with"<<i<<")"
                       <<" (crossoverGene1="<< crossoverPoint1<<","
                       <<" crossoverGene2="<< crossoverPoint2<<")"
                       <<endl; 
                  if(cost<runBestCost)
                  {
                     runBestCost=cost;
                     this->TagNewBestConfig();
                     needUpdateDisplay=true;
                     mRefParList.SaveParamSet(runBestIndex);
                     if(cost<mBestCost)
                     {
                        mBestCost=cost;
                        mRefParList.SaveParamSet(mBestParSavedSetIndex);
                        if(!silent) cout << "->Trial :" << mNbTrial 
                          << " World="<< worldSwapIndex(k)
                          << " Temp="<< simAnnealTemp(k)
                          << " Mutation Ampl.: "<<mMutationAmplitude
                          << " NEW OVERALL Best Cost="<<mBestCost<< "(MATING !)"<<endl;
                     }
                     else if(!silent) cout << "->Trial :" << mNbTrial 
                          << " World="<< worldSwapIndex(k)
                          << " Temp="<< simAnnealTemp(k)
                          << " Mutation Ampl.: "<<mMutationAmplitude
                          << " NEW RUN Best Cost="<<runBestCost<< "(MATING !)"<<endl;
                     bestConfigNb=mNbTrial;
                     if(!silent) this->DisplayReport();
                     //for(int i=0;i<mRefinedObjList.GetNb();i++) 
                     //   mRefinedObjList.GetObj(i).Print();
                  }
                  i=k;//Don't test other Worlds
                  break;
               }
               //mNbTrial++;nbStep--;
               //if((mNbTrial%nbTrialsReport)==0) makeReport=true;
            }
         }
      TAU_PROFILE_STOP(timer1);
      #endif
      if(true==makeReport)
      {
         makeReport=false;
         worldNbAcceptedMoves*=nbWorld;
         if(!silent)
         {
            #if 0
            {// Experimental, dynamical weighting
               REAL max=0.;
               map<const RefinableObj*,REAL> ll,llvar;
               map<const RefinableObj*,LogLikelihoodStats>::iterator pos;
               for(pos=mvContextObjStats[0].begin();pos!=mvContextObjStats[0].end();++pos)
               {
                  ll   [pos->first]=0.;
                  llvar[pos->first]=0.;
               }
               for(int i=0;i<nbWorld;i++)
               {
                  for(pos=mvContextObjStats[0].begin();pos!=mvContextObjStats[0].end();++pos)
                  {
                     ll   [pos->first] += pos->second.mTotalLogLikelihood;
                     llvar[pos->first] += pos->second.mTotalLogLikelihoodDeltaSq;
                  }
               }
               for(pos=mvContextObjStats[0].begin();pos!=mvContextObjStats[0].end();++pos)
               {
                  cout << pos->first->GetName()
                       << " " << llvar[pos->first]
                       << " " << mvObjWeight[pos->first].mWeight
                       << " " << max<<endl;
                  llvar[pos->first] *= mvObjWeight[pos->first].mWeight;
                  if(llvar[pos->first]>max) max=llvar[pos->first];
               }
               map<const RefinableObj*,REAL>::iterator pos2;
               for(pos2=llvar.begin();pos2!=llvar.end();++pos2)
               {
                  const REAL d=pos2->second;
                  if(d<(max/mvObjWeight.size()/10.))
                  {
                     if(d<1) continue;
                     mvObjWeight[pos2->first].mWeight *=2;
                  }
               }
               REAL ll1=0;
               REAL llt=0;
               for(pos2=ll.begin();pos2!=ll.end();++pos2)
               {
                  llt += pos2->second;
                  ll1 += pos2->second * mvObjWeight[pos2->first].mWeight;
               }
               map<const RefinableObj*,DynamicObjWeight>::iterator posw;
               for(posw=mvObjWeight.begin();posw!=mvObjWeight.end();++posw)
               {
                  posw->second.mWeight *= llt/ll1;
               }
            }
            #endif //Experimental dynamical weighting
            #if 1 //def __DEBUG__
            for(int i=0;i<nbWorld;i++)
            {
               cout<<"   World :"<<worldSwapIndex(i)<<":";
               map<const RefinableObj*,LogLikelihoodStats>::iterator pos;
               for(pos=mvContextObjStats[i].begin();pos!=mvContextObjStats[i].end();++pos)
               {
                  cout << pos->first->GetName() 
                       << "(LLK="
                       << pos->second.mLastLogLikelihood
                       //<< "(<LLK>="
                       //<< pos->second.mTotalLogLikelihood/nbTrialsReport
                       //<< ", <delta(LLK)^2>="
                       //<< pos->second.mTotalLogLikelihoodDeltaSq/nbTrialsReport
                       << ", w="<<mvObjWeight[pos->first].mWeight
                       <<")  ";
                  pos->second.mTotalLogLikelihood=0;
                  pos->second.mTotalLogLikelihoodDeltaSq=0;
               }
               cout << endl;
            }
            #endif
            for(int i=0;i<nbWorld;i++)
            {
               //mRefParList.RestoreParamSet(worldCurrentSetIndex(i));
               cout <<"   World :" << worldSwapIndex(i)
                    <<" Temp.: " << simAnnealTemp(i)
                    <<" Mutation Ampl.: " << mutationAmplitude(i)
                    <<" Current Cost=" << currentCost(i)
                    <<" Accepting "
                    << (int)((REAL)worldNbAcceptedMoves(i)/nbTrialsReport*100)
                    <<"% moves " <<endl;
               //     <<"% moves " << mRefParList.GetPar("Pboccup").GetValue()<<endl;
            }
         }
         if(!silent) cout <<"Trial :" << mNbTrial << " Best Cost=" << runBestCost<< " ";
         if(!silent) chrono.print();
         //Change the mutation rate if necessary for each world
         if(ANNEALING_SMART==mAnnealingScheduleMutation.GetChoice())
         {
            for(int i=0;i<nbWorld;i++)
            {
               if((worldNbAcceptedMoves(i)/(REAL)nbTrialsReport)>0.30)
                  mutationAmplitude(i)*=2.;
               if((worldNbAcceptedMoves(i)/(REAL)nbTrialsReport)<0.10)
                  mutationAmplitude(i)/=2.;
               if(mutationAmplitude(i)>mMutationAmplitudeMax) 
                  mutationAmplitude(i)=mMutationAmplitudeMax;
               if(mutationAmplitude(i)<mMutationAmplitudeMin)
                  mutationAmplitude(i)=mMutationAmplitudeMin;
            }
         }
         if(ANNEALING_SMART==mAnnealingScheduleTemp.GetChoice())
         {
            for(int i=0;i<nbWorld;i++)
            {
               if((worldNbAcceptedMoves(i)/(REAL)nbTrialsReport)>0.30)
                  simAnnealTemp(i)/=1.5;
               if((worldNbAcceptedMoves(i)/(REAL)nbTrialsReport)>0.80)
                  simAnnealTemp(i)/=1.5;
               if((worldNbAcceptedMoves(i)/(REAL)nbTrialsReport)>0.95)
                  simAnnealTemp(i)/=1.5;

               if((worldNbAcceptedMoves(i)/(REAL)nbTrialsReport)<0.10)
                  simAnnealTemp(i)*=1.5;
               if((worldNbAcceptedMoves(i)/(REAL)nbTrialsReport)<0.04)
                  simAnnealTemp(i)*=1.5;
               //if((worldNbAcceptedMoves(i)/(REAL)nbTrialsReport)<0.01)
               //   simAnnealTemp(i)*=1.5;
               //cout<<"World#"<<i<<":"<<worldNbAcceptedMoves(i)<<":"<<nbTrialsReport<<endl;
               //if(simAnnealTemp(i)>mTemperatureMax) simAnnealTemp(i)=mTemperatureMax;
               //if(simAnnealTemp(i)<mTemperatureMin) simAnnealTemp(i)=mTemperatureMin;
            }
         }
         worldNbAcceptedMoves=0;
         //this->DisplayReport();

         #ifdef __WX__CRYST__
         if(0!=mpWXCrystObj) mpWXCrystObj->UpdateDisplayNbTrial();
         #endif
      }
      if( (needUpdateDisplay&&(lastUpdateDisplayTime<(chrono.seconds()-1)))||(lastUpdateDisplayTime<(chrono.seconds()-10)))
      {
         mRefParList.RestoreParamSet(runBestIndex);
         this->UpdateDisplay();
         needUpdateDisplay=false;
         lastUpdateDisplayTime=chrono.seconds();
      }
      
      #ifdef __WX__CRYST__
      mMutexStopAfterCycle.Lock();
      #endif
      if((runBestCost<finalcost) || mStopAfterCycle ||( (maxTime>0)&&(chrono.seconds()>maxTime))) 
      {
         #ifdef __WX__CRYST__
         mMutexStopAfterCycle.Unlock();
         #endif
         if(!silent) cout << endl <<endl << "Refinement Stopped:"<<mBestCost<<endl;
         break;
      }
      #ifdef __WX__CRYST__
      mMutexStopAfterCycle.Unlock();
      #endif
   }//Trials
   mLastOptimTime=chrono.seconds();
   //Restore Best values
      //mRefParList.Print();
      if(!silent) this->DisplayReport();
      mRefParList.RestoreParamSet(runBestIndex);
      //for(int i=0;i<mRefinedObjList.GetNb();i++) mRefinedObjList.GetObj(i).Print();
      mCurrentCost=this->GetLogLikelihood();
      if(!silent) cout<<"Run Best Cost:"<<mCurrentCost<<endl;
      if(!silent) chrono.print();
   //Save density of states
      //ofstream out("densityOfStates.txt");
      //out << trialsDensity<<endl;
      //out.close();
   // Clear temporary param set
      for(int i=0;i<nbWorld;i++)
      {
         mRefParList.ClearParamSet(worldCurrentSetIndex(i));
         //mvSavedParamSet.push_back(make_pair(worldCurrentSetIndex(i),currentCost(i)));
      }
      mRefParList.ClearParamSet(lastParSavedSetIndex);
      mRefParList.ClearParamSet(runBestIndex);
}


Generated by  Doxygen 1.6.0   Back to index