FOX/ObjCryst++  1.10.X (development)
GlobalOptimObj.cpp
1 /* ObjCryst++ Object-Oriented Crystallographic Library
2  (c) 2000-2002 Vincent Favre-Nicolin vincefn@users.sourceforge.net
3  2000-2001 University of Geneva (Switzerland)
4 
5  This program is free software; you can redistribute it and/or modify
6  it under the terms of the GNU General Public License as published by
7  the Free Software Foundation; either version 2 of the License, or
8  (at your option) any later version.
9 
10  This program is distributed in the hope that it will be useful,
11  but WITHOUT ANY WARRANTY; without even the implied warranty of
12  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  GNU General Public License for more details.
14 
15  You should have received a copy of the GNU General Public License
16  along with this program; if not, write to the Free Software
17  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
18 */
19 /*
20 * source file for Global Optimization Objects
21 *
22 */
23 #include <iomanip>
24 
25 #include "ObjCryst/RefinableObj/GlobalOptimObj.h"
26 #include "ObjCryst/ObjCryst/Crystal.h"
27 #include "ObjCryst/Quirks/VFNStreamFormat.h"
28 #include "ObjCryst/Quirks/VFNDebug.h"
29 #include "ObjCryst/Quirks/Chronometer.h"
30 #include "ObjCryst/ObjCryst/IO.h"
31 #include "ObjCryst/RefinableObj/LSQNumObj.h"
32 
33 #include "ObjCryst/ObjCryst/Molecule.h"
34 
35 #ifdef __WX__CRYST__
36  #include "ObjCryst/wxCryst/wxRefinableObj.h"
37  #undef GetClassName // Conflict from wxMSW headers ? (cygwin)
38 #endif
39 
40 //For some reason, with wxWindows this must be placed after wx headers (Borland c++)
41 #include <fstream>
42 #include <sstream>
43 #include <stdio.h>
44 #include "boost/format.hpp"
45 
46 namespace ObjCryst
47 {
48 void CompareWorlds(const CrystVector_long &idx,const CrystVector_long &swap, const RefinableObj &obj)
49 {
50  const long nb=swap.numElements();
51  const CrystVector_REAL *pv0=&(obj.GetParamSet(idx(swap(nb-1))));
52  for(long i=0;i<idx.numElements();++i)
53  {
54  REAL d=0.0;
55  const CrystVector_REAL *pv1=&(obj.GetParamSet(idx(swap(i))));
56  for(long j=0;j<pv0->numElements();++j) d += ((*pv0)(j)-(*pv1)(j))*((*pv0)(j)-(*pv1)(j));
57  cout<<"d("<<i<<")="<<sqrt(d)<<endl;
58  }
59  cout<<endl;
60 }
61 //#################################################################################
62 //
63 // OptimizationObj
64 //
65 //#################################################################################
66 ObjRegistry<OptimizationObj> gOptimizationObjRegistry("List of all Optimization objects");
67 
68 
69 OptimizationObj::OptimizationObj(const string name):
70 mName(name),mSaveFileName("GlobalOptim.save"),
71 mNbTrialPerRun(10000000),mNbTrial(0),mBestCost(-1),
72 mBestParSavedSetIndex(-1),
73 mContext(0),
74 mIsOptimizing(false),mStopAfterCycle(false),
75 mRefinedObjList("OptimizationObj: "+mName+" RefinableObj registry"),
76 mRecursiveRefinedObjList("OptimizationObj: "+mName+" recursive RefinableObj registry"),
77 mLastOptimTime(0)
78 {
79  VFN_DEBUG_ENTRY("OptimizationObj::OptimizationObj()",5)
80  // This must be done in a real class to avoid calling a pure virtual method
81  // if a graphical representation is automatically called upon registration.
82  // gOptimizationObjRegistry.Register(*this);
83 
84  static bool need_initRandomSeed=true;
85  if(need_initRandomSeed==true)
86  {
87  srand(time(NULL));
88  need_initRandomSeed=false;
89  }
90  // We only copy parameters, so do not delete them !
91  mRefParList.SetDeleteRefParInDestructor(false);
92  VFN_DEBUG_EXIT("OptimizationObj::OptimizationObj()",5)
93 }
94 
95 OptimizationObj::~OptimizationObj()
96 {
97  VFN_DEBUG_ENTRY("OptimizationObj::~OptimizationObj()",5)
98  gOptimizationObjRegistry.DeRegister(*this);
99  VFN_DEBUG_EXIT("OptimizationObj::~OptimizationObj()",5)
100 }
101 
102 void OptimizationObj::RandomizeStartingConfig()
103 {
104  VFN_DEBUG_ENTRY("OptimizationObj::RandomizeStartingConfig()",5)
105  this->PrepareRefParList();
106  for(int j=0;j<mRefParList.GetNbParNotFixed();j++)
107  {
108  if(true==mRefParList.GetParNotFixed(j).IsLimited())
109  {
110  const REAL min=mRefParList.GetParNotFixed(j).GetMin();
111  const REAL max=mRefParList.GetParNotFixed(j).GetMax();
112  mRefParList.GetParNotFixed(j).MutateTo(min+(max-min)*(rand()/(REAL)RAND_MAX) );
113  }
114  else if(true==mRefParList.GetParNotFixed(j).IsPeriodic())
115  mRefParList.GetParNotFixed(j).
116  Mutate(mRefParList.GetParNotFixed(j).GetPeriod()*rand()/(REAL)RAND_MAX);
117  }
118  //else cout << mRefParList.GetParNotFixed(j).Name() <<" Not limited :-(" <<endl;
119  VFN_DEBUG_EXIT("OptimizationObj::RandomizeStartingConfig()",5)
120 }
121 
122 void OptimizationObj::FixAllPar()
123 {
124  VFN_DEBUG_ENTRY("OptimizationObj::FixAllPar()",5)
125  this->BuildRecursiveRefObjList();
126  for(int i=0;i<mRecursiveRefinedObjList.GetNb();i++)
127  mRecursiveRefinedObjList.GetObj(i).FixAllPar();
128  VFN_DEBUG_EXIT("OptimizationObj::FixAllPar():End",5)
129 }
130 void OptimizationObj::SetParIsFixed(const string& parName,const bool fix)
131 {
132  this->BuildRecursiveRefObjList();
133  for(int i=0;i<mRecursiveRefinedObjList.GetNb();i++)
134  mRecursiveRefinedObjList.GetObj(i).SetParIsFixed(parName,fix);
135 }
136 void OptimizationObj::SetParIsFixed(const RefParType *type,const bool fix)
137 {
138  this->BuildRecursiveRefObjList();
139  for(int i=0;i<mRecursiveRefinedObjList.GetNb();i++)
140  mRecursiveRefinedObjList.GetObj(i).SetParIsFixed(type,fix);
141 }
142 
143 void OptimizationObj::UnFixAllPar()
144 {
145  this->BuildRecursiveRefObjList();
146  for(int i=0;i<mRecursiveRefinedObjList.GetNb();i++)
147  mRecursiveRefinedObjList.GetObj(i).UnFixAllPar();
148 }
149 
150 void OptimizationObj::SetParIsUsed(const string& parName,const bool use)
151 {
152  this->BuildRecursiveRefObjList();
153  for(int i=0;i<mRecursiveRefinedObjList.GetNb();i++)
154  mRecursiveRefinedObjList.GetObj(i).SetParIsUsed(parName,use);
155 }
156 void OptimizationObj::SetParIsUsed(const RefParType *type,const bool use)
157 {
158  this->BuildRecursiveRefObjList();
159  for(int i=0;i<mRecursiveRefinedObjList.GetNb();i++)
160  mRecursiveRefinedObjList.GetObj(i).SetParIsUsed(type,use);
161 }
162 void OptimizationObj::SetLimitsRelative(const string &parName,
163  const REAL min, const REAL max)
164 {
165  this->BuildRecursiveRefObjList();
166  for(int i=0;i<mRecursiveRefinedObjList.GetNb();i++)
167  mRecursiveRefinedObjList.GetObj(i).SetLimitsRelative(parName,min,max);
168 }
169 void OptimizationObj::SetLimitsRelative(const RefParType *type,
170  const REAL min, const REAL max)
171 {
172  this->BuildRecursiveRefObjList();
173  for(int i=0;i<mRecursiveRefinedObjList.GetNb();i++)
174  mRecursiveRefinedObjList.GetObj(i).SetLimitsRelative(type,min,max);
175 }
176 void OptimizationObj::SetLimitsAbsolute(const string &parName,
177  const REAL min, const REAL max)
178 {
179  this->BuildRecursiveRefObjList();
180  for(int i=0;i<mRecursiveRefinedObjList.GetNb();i++)
181  mRecursiveRefinedObjList.GetObj(i).SetLimitsAbsolute(parName,min,max);
182 }
183 void OptimizationObj::SetLimitsAbsolute(const RefParType *type,
184  const REAL min, const REAL max)
185 {
186  this->BuildRecursiveRefObjList();
187  for(int i=0;i<mRecursiveRefinedObjList.GetNb();i++)
188  mRecursiveRefinedObjList.GetObj(i).SetLimitsAbsolute(type,min,max);
189 }
190 
191 REAL OptimizationObj::GetLogLikelihood() const
192 {
193  TAU_PROFILE("OptimizationObj::GetLogLikelihood()","void ()",TAU_DEFAULT);
194  REAL cost =0.;
195  for(int i=0;i<mRecursiveRefinedObjList.GetNb();i++)
196  {
197  const REAL tmp=mRecursiveRefinedObjList.GetObj(i).GetLogLikelihood();
198  if(tmp!=0.)
199  {
200  LogLikelihoodStats* st=&((mvContextObjStats[mContext])
201  [&(mRecursiveRefinedObjList.GetObj(i))]);
202  st->mTotalLogLikelihood += tmp;
203  st->mTotalLogLikelihoodDeltaSq +=
204  (tmp-st->mLastLogLikelihood)*(tmp-st->mLastLogLikelihood);
205  st->mLastLogLikelihood=tmp;
206  }
207  cost += mvObjWeight[&(mRecursiveRefinedObjList.GetObj(i))].mWeight * tmp;
208  }
209  return cost;
210 }
211 void OptimizationObj::StopAfterCycle()
212 {
213  VFN_DEBUG_MESSAGE("OptimizationObj::StopAfterCycle()",5)
214  if(mIsOptimizing)
215  {
216  #ifdef __WX__CRYST__
217  wxMutexLocker lock(mMutexStopAfterCycle);
218  #endif
219  mStopAfterCycle=true;
220  }
221 }
222 
223 void OptimizationObj::DisplayReport()
224 {
225  //:TODO: ask all objects to print their own report ?
226 }
227 
228 void OptimizationObj::AddRefinableObj(RefinableObj &obj)
229 {
230  VFN_DEBUG_MESSAGE("OptimizationObj::AddRefinableObj():"<<obj.GetName(),5)
231  //in case some object has been modified, to avoid rebuilding the entire list
232  this->BuildRecursiveRefObjList();
233 
234  mRefinedObjList.Register(obj);
235  RefObjRegisterRecursive(obj,mRecursiveRefinedObjList);
236  #ifdef __WX__CRYST__
237  if(0!=this->WXGet()) this->WXGet()->AddRefinedObject(obj);
238  #endif
239 }
240 
241 RefinableObj& OptimizationObj::GetFullRefinableObj(const bool rebuild)
242 {
243  if(rebuild) this->PrepareRefParList();
244  return mRefParList;
245 }
246 
247 const string& OptimizationObj::GetName()const { return mName;}
248 void OptimizationObj::SetName(const string& name) {mName=name;}
249 
250 const string OptimizationObj::GetClassName()const { return "OptimizationObj";}
251 
252 void OptimizationObj::Print()const {this->XMLOutput(cout);}
253 
254 void OptimizationObj::RestoreBestConfiguration()
255 {
256  //:TODO: check list of refinableObj has not changed, and the list of
257  // RefPar has not changed in all sub-objects.
258  if(mBestParSavedSetIndex>0) mRefParList.RestoreParamSet(mBestParSavedSetIndex);
259 }
260 
261 bool OptimizationObj::IsOptimizing()const{return mIsOptimizing;}
262 
263 void OptimizationObj::TagNewBestConfig()
264 {
265  for(int i=0;i<mRecursiveRefinedObjList.GetNb();i++)
266  mRecursiveRefinedObjList.GetObj(i).TagNewBestConfig();
267  mMainTracker.AppendValues(mNbTrial);
268 }
269 
270 REAL OptimizationObj::GetLastOptimElapsedTime()const
271 {
272  return mLastOptimTime;
273 }
274 
275 MainTracker& OptimizationObj::GetMainTracker(){return mMainTracker;}
276 
277 const MainTracker& OptimizationObj::GetMainTracker()const{return mMainTracker;}
278 
279 RefObjOpt& OptimizationObj::GetXMLAutoSaveOption() {return mXMLAutoSave;}
280 const RefObjOpt& OptimizationObj::GetXMLAutoSaveOption()const {return mXMLAutoSave;}
281 
282 const REAL& OptimizationObj::GetBestCost()const{return mBestCost;}
283 REAL& OptimizationObj::GetBestCost(){return mBestCost;}
284 
285 void OptimizationObj::BeginOptimization(const bool allowApproximations, const bool enableRestraints)
286 {
287  for(int i=0;i<mRefinedObjList.GetNb();i++)
288  {
289  mRefinedObjList.GetObj(i).BeginOptimization(allowApproximations,enableRestraints);
290  }
291 }
292 
293 void OptimizationObj::EndOptimization()
294 {
295  for(int i=0;i<mRefinedObjList.GetNb();i++) mRefinedObjList.GetObj(i).EndOptimization();
296 }
297 
298 long& OptimizationObj::NbTrialPerRun() {return mNbTrialPerRun;}
299 
300 const long& OptimizationObj::NbTrialPerRun() const {return mNbTrialPerRun;}
301 
302 void OptimizationObj::PrepareRefParList()
303 {
304  VFN_DEBUG_ENTRY("OptimizationObj::PrepareRefParList()",6)
305 
306  this->BuildRecursiveRefObjList();
307  // As any parameter been added in the recursive list of objects ?
308  // or has any object been added/removed ?
309  RefinableObjClock clock;
310  GetRefParListClockRecursive(mRecursiveRefinedObjList,clock);
311  if( (clock>mRefParList.GetRefParListClock())
312  ||(mRecursiveRefinedObjList.GetRegistryClock()>mRefParList.GetRefParListClock()) )
313  {
314  VFN_DEBUG_MESSAGE("OptimizationObj::PrepareRefParList():Rebuild list",6)
315  mRefParList.ResetParList();
316  mRefParList.EraseAllParamSet();
317  for(int i=0;i<mRecursiveRefinedObjList.GetNb();i++)
318  mRefParList.AddPar(mRecursiveRefinedObjList.GetObj(i),true);
319  mvSavedParamSet.clear();
320  mBestParSavedSetIndex=mRefParList.CreateParamSet("Best Configuration");
321  mvSavedParamSet.push_back(make_pair(mBestParSavedSetIndex,mBestCost));
322 
323  mMainTracker.ClearTrackers();
324 
325  REAL (OptimizationObj::*fl)() const;
326  fl=&OptimizationObj::GetLogLikelihood;
327  mMainTracker.AddTracker(new TrackerObject<OptimizationObj>
328  (this->GetName()+"::Overall LogLikelihood",*this,fl));
329 
330  for(long i=0;i<mRecursiveRefinedObjList.GetNb();i++)
331  {
332  REAL (RefinableObj::*fp)() const;
333  fp=&RefinableObj::GetLogLikelihood;
334  mMainTracker.AddTracker(new TrackerObject<RefinableObj>
335  (mRecursiveRefinedObjList.GetObj(i).GetName()+"::LogLikelihood",mRecursiveRefinedObjList.GetObj(i),fp));
336 
337  if(mRecursiveRefinedObjList.GetObj(i).GetClassName()=="Crystal")
338  {
339  REAL (Crystal::*fc)() const;
340  const Crystal *pCryst=dynamic_cast<const Crystal *>(&(mRecursiveRefinedObjList.GetObj(i)));
341  fc=&Crystal::GetBumpMergeCost;
342  mMainTracker.AddTracker(new TrackerObject<Crystal>
343  (pCryst->GetName()+"::BumpMergeCost",*pCryst,fc));
344  fc=&Crystal::GetBondValenceCost;
345  mMainTracker.AddTracker(new TrackerObject<Crystal>
346  (pCryst->GetName()+"::BondValenceCost",*pCryst,fc));
347  }
348  }
349  }
350  // Prepare for refinement, ie get the list of not fixed parameters,
351  // and prepare the objects...
352  mRefParList.PrepareForRefinement();
353  for(int i=0;i<mRecursiveRefinedObjList.GetNb();i++)
354  mRecursiveRefinedObjList.GetObj(i).PrepareForRefinement();
355  VFN_DEBUG_EXIT("OptimizationObj::PrepareRefParList()",6)
356 }
357 
358 void OptimizationObj::InitOptions()
359 {
360  VFN_DEBUG_MESSAGE("OptimizationObj::InitOptions()",5)
361  static string xmlAutoSaveName;
362  static string xmlAutoSaveChoices[6];
363 
364  static bool needInitNames=true;
365  if(true==needInitNames)
366  {
367  xmlAutoSaveName="Save Best Config Regularly";
368  xmlAutoSaveChoices[0]="No";
369  xmlAutoSaveChoices[1]="Every day";
370  xmlAutoSaveChoices[2]="Every hour";
371  xmlAutoSaveChoices[3]="Every 10mn";
372  xmlAutoSaveChoices[4]="Every new best config (a lot ! Not Recommended !)";
373  xmlAutoSaveChoices[5]="Every Run (Recommended)";
374 
375  needInitNames=false;//Only once for the class
376  }
377  mXMLAutoSave.Init(6,&xmlAutoSaveName,xmlAutoSaveChoices);
378  VFN_DEBUG_MESSAGE("OptimizationObj::InitOptions():End",5)
379 }
380 
381 void OptimizationObj::UpdateDisplay()
382 {
383  Chronometer chrono;
384  #ifdef __WX__CRYST__
385  if(0!=this->WXGet()) this->WXGet()->CrystUpdate(true,true);
386  #endif
387  for(int i=0;i<mRefinedObjList.GetNb();i++)
388  mRefinedObjList.GetObj(i).UpdateDisplay();
389  mMainTracker.UpdateDisplay();
390 }
391 void OptimizationObj::BuildRecursiveRefObjList()
392 {
393  // First check if anything has changed (ie if a sub-object has been
394  // added or removed in the recursive refinable object list)
395  RefinableObjClock clock;
396  GetSubRefObjListClockRecursive(mRefinedObjList,clock);
397  if(clock>mRecursiveRefinedObjList.GetRegistryClock())
398  {
399  VFN_DEBUG_ENTRY("OptimizationObj::BuildRecursiveRefObjList()",5)
400  mRecursiveRefinedObjList.DeRegisterAll();
401  for(int i=0;i<mRefinedObjList.GetNb();i++)
402  RefObjRegisterRecursive(mRefinedObjList.GetObj(i),mRecursiveRefinedObjList);
403  VFN_DEBUG_EXIT("OptimizationObj::BuildRecursiveRefObjList()",5)
404  }
405 }
406 
407 //#################################################################################
408 //
409 // MonteCarloObj
410 //
411 //#################################################################################
412 MonteCarloObj::MonteCarloObj(const string name):
413 OptimizationObj(name),
414 mCurrentCost(-1),
415 mTemperatureMax(1e6),mTemperatureMin(.001),mTemperatureGamma(1.0),
416 mMutationAmplitudeMax(8.),mMutationAmplitudeMin(.125),mMutationAmplitudeGamma(1.0),
417 mNbTrialRetry(0),mMinCostRetry(0)
418 #ifdef __WX__CRYST__
419 ,mpWXCrystObj(0)
420 #endif
421 {
422  VFN_DEBUG_ENTRY("MonteCarloObj::MonteCarloObj()",5)
423  this->InitOptions();
424  mGlobalOptimType.SetChoice(GLOBAL_OPTIM_PARALLEL_TEMPERING);
425  mAnnealingScheduleTemp.SetChoice(ANNEALING_SMART);
426  mAnnealingScheduleMutation.SetChoice(ANNEALING_EXPONENTIAL);
427  mXMLAutoSave.SetChoice(5);//Save after each Run
428  mAutoLSQ.SetChoice(0);
429  gOptimizationObjRegistry.Register(*this);
430  VFN_DEBUG_EXIT("MonteCarloObj::MonteCarloObj()",5)
431 }
432 
433 MonteCarloObj::MonteCarloObj(const bool internalUseOnly):
434 OptimizationObj(""),
435 mCurrentCost(-1),
436 mTemperatureMax(.03),mTemperatureMin(.003),mTemperatureGamma(1.0),
437 mMutationAmplitudeMax(16.),mMutationAmplitudeMin(.125),mMutationAmplitudeGamma(1.0),
438 mNbTrialRetry(0),mMinCostRetry(0)
439 #ifdef __WX__CRYST__
440 ,mpWXCrystObj(0)
441 #endif
442 {
443  VFN_DEBUG_ENTRY("MonteCarloObj::MonteCarloObj(bool)",5)
444  this->InitOptions();
445  mGlobalOptimType.SetChoice(GLOBAL_OPTIM_PARALLEL_TEMPERING);
446  mAnnealingScheduleTemp.SetChoice(ANNEALING_SMART);
447  mAnnealingScheduleMutation.SetChoice(ANNEALING_EXPONENTIAL);
448  mXMLAutoSave.SetChoice(5);//Save after each Run
449  mAutoLSQ.SetChoice(0);
450  if(false==internalUseOnly) gOptimizationObjRegistry.Register(*this);
451  VFN_DEBUG_EXIT("MonteCarloObj::MonteCarloObj(bool)",5)
452 }
453 
454 MonteCarloObj::~MonteCarloObj()
455 {
456  VFN_DEBUG_ENTRY("MonteCarloObj::~MonteCarloObj()",5)
457  gOptimizationObjRegistry.DeRegister(*this);
458  VFN_DEBUG_EXIT ("MonteCarloObj::~MonteCarloObj()",5)
459 }
460 void MonteCarloObj::SetAlgorithmSimulAnnealing(const AnnealingSchedule scheduleTemp,
461  const REAL tMax, const REAL tMin,
462  const AnnealingSchedule scheduleMutation,
463  const REAL mutMax, const REAL mutMin,
464  const long nbTrialRetry,const REAL minCostRetry)
465 {
466  VFN_DEBUG_MESSAGE("MonteCarloObj::SetAlgorithmSimulAnnealing()",5)
467  mGlobalOptimType.SetChoice(GLOBAL_OPTIM_SIMULATED_ANNEALING);
468  mTemperatureMax=tMax;
469  mTemperatureMin=tMin;
470  mAnnealingScheduleTemp.SetChoice(scheduleTemp);
471 
472 
473  mMutationAmplitudeMax=mutMax;
474  mMutationAmplitudeMin=mutMin;
475  mAnnealingScheduleMutation.SetChoice(scheduleMutation);
476  mNbTrialRetry=nbTrialRetry;
477  mMinCostRetry=minCostRetry;
478  VFN_DEBUG_MESSAGE("MonteCarloObj::SetAlgorithmSimulAnnealing():End",3)
479 }
480 
481 void MonteCarloObj::SetAlgorithmParallTempering(const AnnealingSchedule scheduleTemp,
482  const REAL tMax, const REAL tMin,
483  const AnnealingSchedule scheduleMutation,
484  const REAL mutMax, const REAL mutMin)
485 {
486  VFN_DEBUG_MESSAGE("MonteCarloObj::SetAlgorithmParallTempering()",5)
487  mGlobalOptimType.SetChoice(GLOBAL_OPTIM_PARALLEL_TEMPERING);
488  mTemperatureMax=tMax;
489  mTemperatureMin=tMin;
490  mAnnealingScheduleTemp.SetChoice(scheduleTemp);
491 
492  mMutationAmplitudeMax=mutMax;
493  mMutationAmplitudeMin=mutMin;
494  mAnnealingScheduleMutation.SetChoice(scheduleMutation);
495  //mNbTrialRetry=nbTrialRetry;
496  //mMinCostRetry=minCostRetry;
497  VFN_DEBUG_MESSAGE("MonteCarloObj::SetAlgorithmParallTempering():End",3)
498 }
499 void MonteCarloObj::Optimize(long &nbStep,const bool silent,const REAL finalcost,
500  const REAL maxTime)
501 {
502  //:TODO: Other algorithms !
503  TAU_PROFILE("MonteCarloObj::Optimize()","void (long)",TAU_DEFAULT);
504  VFN_DEBUG_ENTRY("MonteCarloObj::Optimize()",5)
505  this->BeginOptimization(true);
506  this->PrepareRefParList();
507 
508  this->InitLSQ(false);
509 
510  mIsOptimizing=true;
511  if(mTemperatureGamma<0.1) mTemperatureGamma= 0.1;
512  if(mTemperatureGamma>10.0)mTemperatureGamma=10.0;
513  if(mMutationAmplitudeGamma<0.1) mMutationAmplitudeGamma= 0.1;
514  if(mMutationAmplitudeGamma>10.0)mMutationAmplitudeGamma=10.0;
515  // prepare all objects
516  this->TagNewBestConfig();
517  mCurrentCost=this->GetLogLikelihood();
518  mBestCost=mCurrentCost;
519  mvObjWeight.clear();
520  mMainTracker.ClearValues();
521  switch(mGlobalOptimType.GetChoice())
522  {
523  case GLOBAL_OPTIM_SIMULATED_ANNEALING:
524  {
525  this->RunSimulatedAnnealing(nbStep,silent,finalcost,maxTime);
526  break;
527  }//case GLOBAL_OPTIM_SIMULATED_ANNEALING
528  case GLOBAL_OPTIM_PARALLEL_TEMPERING:
529  {
530  this->RunParallelTempering(nbStep,silent,finalcost,maxTime);
531  break;
532  }//case GLOBAL_OPTIM_PARALLEL_TEMPERING
533  case GLOBAL_OPTIM_RANDOM_LSQ: //:TODO:
534  {
535  long cycles = 1;
536  this->RunRandomLSQMethod(cycles);
537  break;
538  }//case GLOBAL_OPTIM_GENETIC
539  }
540  mIsOptimizing=false;
541  #ifdef __WX__CRYST__
542  mMutexStopAfterCycle.Lock();
543  #endif
544  mStopAfterCycle=false;
545  #ifdef __WX__CRYST__
546  mMutexStopAfterCycle.Unlock();
547  #endif
548 
549  mRefParList.RestoreParamSet(mBestParSavedSetIndex);
550  this->EndOptimization();
551 
552  if(mSaveTrackedData.GetChoice()==1)
553  {
554  ofstream outTracker;
555  const string outTrackerName=this->GetName()+"-Tracker.dat";
556  outTracker.open(outTrackerName.c_str());
557  mMainTracker.SaveAll(outTracker);
558  outTracker.close();
559  }
560 
561  for(vector<pair<long,REAL> >::iterator pos=mvSavedParamSet.begin();pos!=mvSavedParamSet.end();++pos)
562  if(pos->first==mBestParSavedSetIndex)
563  {
564  if( (pos->second>mBestCost)
565  ||(pos->second<0))
566  {
567  pos->second=mBestCost;
568  break;
569  }
570  }
571 
572  this->UpdateDisplay();
573 
574  VFN_DEBUG_EXIT("MonteCarloObj::Optimize()",5)
575 }
576 void MonteCarloObj::MultiRunOptimize(long &nbCycle,long &nbStep,const bool silent,
577  const REAL finalcost,const REAL maxTime)
578 {
579  //:TODO: Other algorithms !
580  TAU_PROFILE("MonteCarloObj::MultiRunOptimize()","void (long)",TAU_DEFAULT);
581  VFN_DEBUG_ENTRY("MonteCarloObj::MultiRunOptimize()",5)
582  //Keep a copy of the total number of steps, and decrement nbStep
583  const long nbStep0=nbStep;
584  this->BeginOptimization(true);
585  this->PrepareRefParList();
586 
587  this->InitLSQ(false);
588 
589  mIsOptimizing=true;
590  if(mTemperatureGamma<0.1) mTemperatureGamma= 0.1;
591  if(mTemperatureGamma>10.0)mTemperatureGamma=10.0;
592  if(mMutationAmplitudeGamma<0.1) mMutationAmplitudeGamma= 0.1;
593  if(mMutationAmplitudeGamma>10.0)mMutationAmplitudeGamma=10.0;
594  // prepare all objects
595  mCurrentCost=this->GetLogLikelihood();
596  mBestCost=mCurrentCost;
597  this->TagNewBestConfig();
598  mvObjWeight.clear();
599  long nbTrialCumul=0;
600  const long nbCycle0=nbCycle;
601  Chronometer chrono;
602  while(nbCycle!=0)
603  {
604  if(!silent) cout <<"MonteCarloObj::MultiRunOptimize: Starting Run#"<<abs(nbCycle)<<endl;
605  nbStep=nbStep0;
606  for(int i=0;i<mRefinedObjList.GetNb();i++) mRefinedObjList.GetObj(i).RandomizeConfiguration();
607  mMainTracker.ClearValues();
608  chrono.start();
609  switch(mGlobalOptimType.GetChoice())
610  {
611  case GLOBAL_OPTIM_SIMULATED_ANNEALING:
612  {
613  try{this->RunSimulatedAnnealing(nbStep,silent,finalcost,maxTime);}
614  catch(...){cout<<"Unhandled exception in MonteCarloObj::MultiRunOptimize() ?"<<endl;}
615  break;
616  }
617  case GLOBAL_OPTIM_PARALLEL_TEMPERING:
618  {
619  try{this->RunParallelTempering(nbStep,silent,finalcost,maxTime);}
620  catch(...){cout<<"Unhandled exception in MonteCarloObj::MultiRunOptimize() ?"<<endl;}
621  break;
622  }
623  case GLOBAL_OPTIM_RANDOM_LSQ:
624  {
625  try{this->RunRandomLSQMethod(nbCycle);}
626  catch(...){cout<<"Unhandled exception in MonteCarloObj::RunRandomLSQMethod() ?"<<endl;}
627  //nbCycle=1;
628  break;
629  }
630  }
631  nbTrialCumul+=(nbStep0-nbStep);
632  if(finalcost>1)
633  cout<<"Finished Run #"<<nbCycle0-nbCycle<<", final cost="
634  <<this->GetLogLikelihood()<<", nbTrial="<< nbStep0-nbStep<<" ("<<chrono.seconds()
635  <<" seconds), so far <nbTrial>="<< nbTrialCumul/(nbCycle0-nbCycle+1)<<endl;
636  else
637  cout<<"Finished Run #"<<nbCycle0-nbCycle<<", final cost="
638  <<this->GetLogLikelihood()<<", nbTrial="<< nbStep0-nbStep<<" ("<<chrono.seconds()
639  <<" seconds)"<<endl;
640 
641  nbStep=nbStep0;
642  this->UpdateDisplay();
643  stringstream s;
644  s<<"Run #"<<abs(nbCycle);
645  mvSavedParamSet.push_back(make_pair(mRefParList.CreateParamSet(s.str()),mCurrentCost));
646  if(!silent) cout <<"MonteCarloObj::MultiRunOptimize: Finished Run#"
647  <<abs(nbCycle)<<", Run Best Cost:"<<mCurrentCost
648  <<", Overall Best Cost:"<<mBestCost<<endl;
649  if(mXMLAutoSave.GetChoice()==5)
650  {
651  string saveFileName=this->GetName();
652  time_t date=time(0);
653  char strDate[40];
654  strftime(strDate,sizeof(strDate),"%Y-%m-%d_%H-%M-%S",localtime(&date));//%Y-%m-%dT%H:%M:%S%Z
655  char costAsChar[30];
656  sprintf(costAsChar,"-Run#%ld-Cost-%f",abs(nbCycle),this->GetLogLikelihood());
657  saveFileName=saveFileName+(string)strDate+(string)costAsChar+(string)".xml";
658  XMLCrystFileSaveGlobal(saveFileName);
659  }
660  if(mSaveTrackedData.GetChoice()==1)
661  {
662  ofstream outTracker;
663  char runNum[40];
664  sprintf(runNum,"-Tracker-Run#%ld.dat",abs(nbCycle));
665  const string outTrackerName=this->GetName()+runNum;
666  outTracker.open(outTrackerName.c_str());
667  mMainTracker.SaveAll(outTracker);
668  outTracker.close();
669  }
670  nbCycle--;
671  #ifdef __WX__CRYST__
672  mMutexStopAfterCycle.Lock();
673  #endif
674  if(mStopAfterCycle)
675  {
676  #ifdef __WX__CRYST__
677  mMutexStopAfterCycle.Unlock();
678  #endif
679  break;
680  }
681  #ifdef __WX__CRYST__
682  mMutexStopAfterCycle.Unlock();
683  #endif
684  }
685  mIsOptimizing=false;
686 
687  #ifdef __WX__CRYST__
688  mMutexStopAfterCycle.Lock();
689  #endif
690  mStopAfterCycle=false;
691  #ifdef __WX__CRYST__
692  mMutexStopAfterCycle.Unlock();
693  #endif
694 
695  mRefParList.RestoreParamSet(mBestParSavedSetIndex);
696 
697  for(vector<pair<long,REAL> >::iterator pos=mvSavedParamSet.begin();pos!=mvSavedParamSet.end();++pos)
698  if(pos->first==mBestParSavedSetIndex)
699  {
700  if( (pos->second>mBestCost)
701  ||(pos->second<0))
702  {
703  pos->second=mBestCost;
704  break;
705  }
706  }
707 
708  this->EndOptimization();
709 
710  this->UpdateDisplay();
711 
712  if(finalcost>1)
713  cout<<endl<<"Finished all runs, number of trials to reach cost="
714  <<finalcost<<" : <nbTrial>="<<nbTrialCumul/(nbCycle0-nbCycle)<<endl;
715  VFN_DEBUG_EXIT("MonteCarloObj::MultiRunOptimize()",5)
716 }
717 
718 void MonteCarloObj::RunSimulatedAnnealing(long &nbStep,const bool silent,
719  const REAL finalcost,const REAL maxTime)
720 {
721  //Keep a copy of the total number of steps, and decrement nbStep
722  const long nbSteps=nbStep;
723  unsigned int accept;// 1 if last trial was accepted? 2 if new best config ? else 0
724  mNbTrial=0;
725  // time (in seconds) when last autoSave was made (if enabled)
726  unsigned long secondsWhenAutoSave=0;
727 
728  if(!silent) cout << "Starting Simulated Annealing Optimization for"<<nbSteps<<" trials"<<endl;
729  if(!silent) this->DisplayReport();
730  REAL runBestCost;
731  mCurrentCost=this->GetLogLikelihood();
732  runBestCost=mCurrentCost;
733  const long lastParSavedSetIndex=mRefParList.CreateParamSet("MonteCarloObj:Last parameters (SA)");
734  const long runBestIndex=mRefParList.CreateParamSet("Best parameters for current run (SA)");
735  //Report each ... cycles
736  const int nbTryReport=3000;
737  // Keep record of the number of accepted moves
738  long nbAcceptedMoves=0;//since last report
739  long nbAcceptedMovesTemp=0;//since last temperature/mutation rate change
740  // Number of tries since best configuration found
741  long nbTriesSinceBest=0;
742  // Change temperature (and mutation) every...
743  const int nbTryPerTemp=300;
744 
745  mTemperature=sqrt(mTemperatureMin*mTemperatureMax);
746  mMutationAmplitude=sqrt(mMutationAmplitudeMin*mMutationAmplitudeMax);
747 
748  // Do we need to update the display ?
749  bool needUpdateDisplay=false;
750  Chronometer chrono;
751  chrono.start();
752  for(mNbTrial=1;mNbTrial<=nbSteps;)
753  {
754  if((mNbTrial % nbTryPerTemp) == 1)
755  {
756  VFN_DEBUG_MESSAGE("-> Updating temperature and mutation amplitude.",3)
757  // Temperature & displacements amplitude
758  switch(mAnnealingScheduleTemp.GetChoice())
759  {
760  case ANNEALING_BOLTZMANN:
761  mTemperature=
762  mTemperatureMin*log((REAL)nbSteps)/log((REAL)(mNbTrial+1));break;
763  case ANNEALING_CAUCHY:
764  mTemperature=mTemperatureMin*nbSteps/mNbTrial;break;
765  //case ANNEALING_QUENCHING:
766  case ANNEALING_EXPONENTIAL:
767  mTemperature=mTemperatureMax
768  *pow(mTemperatureMin/mTemperatureMax,
769  mNbTrial/(REAL)nbSteps);break;
770  case ANNEALING_GAMMA:
771  mTemperature=mTemperatureMax+(mTemperatureMin-mTemperatureMax)
772  *pow(mNbTrial/(REAL)nbSteps,mTemperatureGamma);break;
773  case ANNEALING_SMART:
774  {
775  if((nbAcceptedMovesTemp/(REAL)nbTryPerTemp)>0.30)
776  mTemperature/=1.5;
777  if((nbAcceptedMovesTemp/(REAL)nbTryPerTemp)<0.10)
778  mTemperature*=1.5;
779  if(mTemperature>mTemperatureMax) mTemperature=mTemperatureMax;
780  if(mTemperature<mTemperatureMin) mTemperature=mTemperatureMin;
781  nbAcceptedMovesTemp=0;
782  break;
783  }
784  default: mTemperature=mTemperatureMin;break;
785  }
786  switch(mAnnealingScheduleMutation.GetChoice())
787  {
788  case ANNEALING_BOLTZMANN:
789  mMutationAmplitude=
790  mMutationAmplitudeMin*log((REAL)nbSteps)/log((REAL)(mNbTrial+1));
791  break;
792  case ANNEALING_CAUCHY:
793  mMutationAmplitude=mMutationAmplitudeMin*nbSteps/mNbTrial;break;
794  //case ANNEALING_QUENCHING:
795  case ANNEALING_EXPONENTIAL:
796  mMutationAmplitude=mMutationAmplitudeMax
797  *pow(mMutationAmplitudeMin/mMutationAmplitudeMax,
798  mNbTrial/(REAL)nbSteps);break;
799  case ANNEALING_GAMMA:
800  mMutationAmplitude=mMutationAmplitudeMax+(mMutationAmplitudeMin-mMutationAmplitudeMax)
801  *pow(mNbTrial/(REAL)nbSteps,mMutationAmplitudeGamma);break;
802  case ANNEALING_SMART:
803  if((nbAcceptedMovesTemp/(REAL)nbTryPerTemp)>0.3) mMutationAmplitude*=2.;
804  if((nbAcceptedMovesTemp/(REAL)nbTryPerTemp)<0.1) mMutationAmplitude/=2.;
805  if(mMutationAmplitude>mMutationAmplitudeMax)
806  mMutationAmplitude=mMutationAmplitudeMax;
807  if(mMutationAmplitude<mMutationAmplitudeMin)
808  mMutationAmplitude=mMutationAmplitudeMax;
809  nbAcceptedMovesTemp=0;
810  break;
811  default: mMutationAmplitude=mMutationAmplitudeMin;break;
812  }
813  }
814 
815  this->NewConfiguration();
816  accept=0;
817  REAL cost=this->GetLogLikelihood();
818  if(cost<mCurrentCost)
819  {
820  accept=1;
821  mCurrentCost=cost;
822  mRefParList.SaveParamSet(lastParSavedSetIndex);
823  if(mCurrentCost<runBestCost)
824  {
825  accept=2;
826  runBestCost=mCurrentCost;
827  this->TagNewBestConfig();
828  needUpdateDisplay=true;
829  mRefParList.SaveParamSet(runBestIndex);
830  if(runBestCost<mBestCost)
831  {
832  mBestCost=mCurrentCost;
833  mRefParList.SaveParamSet(mBestParSavedSetIndex);
834  if(!silent) cout << "Trial :" << mNbTrial
835  << " Temp="<< mTemperature
836  << " Mutation Ampl.: "<<mMutationAmplitude
837  << " NEW OVERALL Best Cost="<<runBestCost<< endl;
838  }
839  else if(!silent) cout << "Trial :" << mNbTrial
840  << " Temp="<< mTemperature
841  << " Mutation Ampl.: "<<mMutationAmplitude
842  << " NEW Run Best Cost="<<runBestCost<< endl;
843  nbTriesSinceBest=0;
844  }
845  nbAcceptedMoves++;
846  nbAcceptedMovesTemp++;
847  }
848  else
849  {
850  if( log((rand()+1)/(REAL)RAND_MAX) < (-(cost-mCurrentCost)/mTemperature) )
851  {
852  accept=1;
853  mCurrentCost=cost;
854  mRefParList.SaveParamSet(lastParSavedSetIndex);
855  nbAcceptedMoves++;
856  nbAcceptedMovesTemp++;
857  }
858  }
859  if(accept==0) mRefParList.RestoreParamSet(lastParSavedSetIndex);
860 
861  if( (mNbTrial % nbTryReport) == 0)
862  {
863  if(!silent) cout <<"Trial :" << mNbTrial << " Temp="<< mTemperature;
864  if(!silent) cout <<" Mutation Ampl.: " <<mMutationAmplitude<< " Best Cost=" << runBestCost
865  <<" Current Cost=" << mCurrentCost
866  <<" Accepting "<<(int)((REAL)nbAcceptedMoves/nbTryReport*100)
867  <<"% moves" << endl;
868  nbAcceptedMoves=0;
869  #ifdef __WX__CRYST__
870  if(0!=mpWXCrystObj) mpWXCrystObj->UpdateDisplayNbTrial();
871  #endif
872  }
873  mNbTrial++;nbStep--;
874 
875  #ifdef __WX__CRYST__
876  mMutexStopAfterCycle.Lock();
877  #endif
878  if((runBestCost<finalcost) || mStopAfterCycle ||( (maxTime>0)&&(chrono.seconds()>maxTime)))
879  {
880  #ifdef __WX__CRYST__
881  mMutexStopAfterCycle.Unlock();
882  #endif
883  if(!silent) cout << endl <<endl << "Refinement Stopped."<<endl;
884  break;
885  }
886  #ifdef __WX__CRYST__
887  mMutexStopAfterCycle.Unlock();
888  #endif
889  nbTriesSinceBest++;
890  if( ((mXMLAutoSave.GetChoice()==1)&&((chrono.seconds()-secondsWhenAutoSave)>86400))
891  ||((mXMLAutoSave.GetChoice()==2)&&((chrono.seconds()-secondsWhenAutoSave)>3600))
892  ||((mXMLAutoSave.GetChoice()==3)&&((chrono.seconds()-secondsWhenAutoSave)> 600))
893  ||((mXMLAutoSave.GetChoice()==4)&&(accept==2)) )
894  {
895  secondsWhenAutoSave=(unsigned long)chrono.seconds();
896  string saveFileName=this->GetName();
897  time_t date=time(0);
898  char strDate[40];
899  strftime(strDate,sizeof(strDate),"%Y-%m-%d_%H-%M-%S",localtime(&date));//%Y-%m-%dT%H:%M:%S%Z
900  char costAsChar[30];
901  if(accept!=2) mRefParList.RestoreParamSet(mBestParSavedSetIndex);
902  sprintf(costAsChar,"-Cost-%f",this->GetLogLikelihood());
903  saveFileName=saveFileName+(string)strDate+(string)costAsChar+(string)".xml";
904  XMLCrystFileSaveGlobal(saveFileName);
905  if(accept!=2) mRefParList.RestoreParamSet(lastParSavedSetIndex);
906  }
907  if((mNbTrial%300==0)&&needUpdateDisplay)
908  {
909  this->UpdateDisplay();
910  needUpdateDisplay=false;
911  mRefParList.RestoreParamSet(lastParSavedSetIndex);
912  }
913 
914  }
915  //cout<<"Beginning final LSQ refinement? ... ";
916  if(mAutoLSQ.GetChoice()>0)
917  {// LSQ
918  if(!silent) cout<<"Beginning final LSQ refinement"<<endl;
919  for(int i=0;i<mRefinedObjList.GetNb();i++) mRefinedObjList.GetObj(i).SetApproximationFlag(false);
920  mRefParList.RestoreParamSet(runBestIndex);
921  mCurrentCost=this->GetLogLikelihood();
922  try {mLSQ.Refine(-50,true,true,false,0.001);}
923  catch(const ObjCrystException &except){};
924  if(!silent) cout<<"LSQ cost: "<<mCurrentCost<<" -> "<<this->GetLogLikelihood()<<endl;
925 
926  // Need to go back to optimization with approximations allowed (they are not during LSQ)
927  for(int i=0;i<mRefinedObjList.GetNb();i++) mRefinedObjList.GetObj(i).SetApproximationFlag(true);
928 
929  REAL cost=this->GetLogLikelihood();
930  if(cost<mCurrentCost)
931  {
932  mCurrentCost=cost;
933  mRefParList.SaveParamSet(lastParSavedSetIndex);
934  if(mCurrentCost<runBestCost)
935  {
936  runBestCost=mCurrentCost;
937  mRefParList.SaveParamSet(runBestIndex);
938  if(runBestCost<mBestCost)
939  {
940  mBestCost=mCurrentCost;
941  mRefParList.SaveParamSet(mBestParSavedSetIndex);
942  if(!silent) cout << "LSQ : NEW OVERALL Best Cost="<<runBestCost<< endl;
943  }
944  else if(!silent) cout << " LSQ : NEW Run Best Cost="<<runBestCost<< endl;
945  }
946  }
947  if(!silent) cout<<"Finished LSQ refinement"<<endl;
948  }
949 
950 
951  mLastOptimTime=chrono.seconds();
952  //Restore Best values
953  mRefParList.RestoreParamSet(runBestIndex);
954  mRefParList.ClearParamSet(runBestIndex);
955  mRefParList.ClearParamSet(lastParSavedSetIndex);
956  mCurrentCost=this->GetLogLikelihood();
957  if(!silent) this->DisplayReport();
958  if(!silent) chrono.print();
959 }
960 /*
961 void MonteCarloObj::RunNondestructiveLSQRefinement( int nbCycle,bool useLevenbergMarquardt,
962  const bool silent, const bool callBeginEndOptimization,
963  const float minChi2var )
964 {
965  float bsigma=-1, bdelta=-1;
966  float asigma=-1, adelta=-1;
967  //set the sigma values lower - it makes the molecular model more stable for LSQ
968  for(int i=0;i<mRefinedObjList.GetNb();i++) {
969  if(mRefinedObjList.GetObj(i).GetClassName()=="Crystal") {
970  try {
971  Crystal * pCryst = dynamic_cast<Crystal *>(&(mRefinedObjList.GetObj(i)));
972  for(int s=0;s<pCryst->GetScattererRegistry().GetNb();s++) {
973  Molecule *pMol=dynamic_cast<Molecule*>(&(pCryst->GetScatt(s)));
974  if(pMol==NULL) continue; // not a Molecule
975  for(vector<MolBond*>::iterator pos = pMol->GetBondList().begin(); pos != pMol->GetBondList().end();++pos) {
976  bsigma = (*pos)->GetLengthSigma();
977  bdelta = (*pos)->GetLengthDelta();
978  (*pos)->SetLengthDelta(0.02);
979  (*pos)->SetLengthSigma(0.001);
980  }
981  for(vector<MolBondAngle*>::iterator pos=pMol->GetBondAngleList().begin();pos != pMol->GetBondAngleList().end();++pos)
982  {
983  asigma = (*pos)->GetAngleSigma();
984  adelta = (*pos)->GetAngleDelta();
985  (*pos)->SetAngleDelta(0.2*DEG2RAD);
986  (*pos)->SetAngleSigma(0.01*DEG2RAD);
987  }
988  }
989  } catch (const std::bad_cast& e) {
990 
991  }
992  }
993  }
994  for(int i=0;i<mRefinedObjList.GetNb();i++) mRefinedObjList.GetObj(i).SetApproximationFlag(false);
995  try {
996  mLSQ.Refine(nbCycle,useLevenbergMarquardt,silent,callBeginEndOptimization,minChi2var);
997  }
998  catch(const ObjCrystException &except) {
999 
1000  };
1001  for(int i=0;i<mRefinedObjList.GetNb();i++) mRefinedObjList.GetObj(i).SetApproximationFlag(true);
1002 
1003  if(bsigma<0 || bdelta<0 || asigma<0 || adelta<0) return;
1004  //restore the delta and sigma values
1005  for(int i=0;i<mRefinedObjList.GetNb();i++) {
1006  if(mRefinedObjList.GetObj(i).GetClassName()=="Crystal") {
1007  try {
1008  Crystal * pCryst = dynamic_cast<Crystal *>(&(mRefinedObjList.GetObj(i)));
1009  for(int s=0;s<pCryst->GetScattererRegistry().GetNb();s++) {
1010  Molecule *pMol=dynamic_cast<Molecule*>(&(pCryst->GetScatt(s)));
1011  if(pMol==NULL) continue; // not a Molecule
1012  for(vector<MolBond*>::iterator pos = pMol->GetBondList().begin(); pos != pMol->GetBondList().end();++pos) {
1013  (*pos)->SetLengthDelta(bdelta);
1014  (*pos)->SetLengthSigma(bsigma);
1015  }
1016  for(vector<MolBondAngle*>::iterator pos=pMol->GetBondAngleList().begin();pos != pMol->GetBondAngleList().end();++pos)
1017  {
1018  (*pos)->SetAngleDelta(adelta);
1019  (*pos)->SetAngleSigma(asigma);
1020  }
1021  }
1022  } catch (const std::bad_cast& e) {
1023 
1024  }
1025  }
1026  }
1027 }
1028 */
1029 void MonteCarloObj::RunRandomLSQMethod(long &nbCycle)
1030 {
1031  //perform random move
1032  mMutationAmplitude=mMutationAmplitudeMax;
1033  float bsigma=-1, bdelta=-1;
1034  float asigma=-1, adelta=-1;
1035 
1036  //set the delta and sigma values - low values are good for LSQ!
1037  for(int i=0;i<mRefinedObjList.GetNb();i++) {
1038  if(mRefinedObjList.GetObj(i).GetClassName()=="Crystal") {
1039  try {
1040  Crystal * pCryst = dynamic_cast<Crystal *>(&(mRefinedObjList.GetObj(i)));
1041  for(int s=0;s<pCryst->GetScattererRegistry().GetNb();s++)
1042  {
1043  Molecule *pMol=dynamic_cast<Molecule*>(&(pCryst->GetScatt(s)));
1044  if(pMol==NULL) continue; // not a Molecule
1045  for(vector<MolBond*>::iterator pos = pMol->GetBondList().begin(); pos != pMol->GetBondList().end();++pos) {
1046  bsigma = (*pos)->GetLengthSigma();
1047  bdelta = (*pos)->GetLengthDelta();
1048  (*pos)->SetLengthDelta(0.02);
1049  (*pos)->SetLengthSigma(0.001);
1050  }
1051  for(vector<MolBondAngle*>::iterator pos=pMol->GetBondAngleList().begin();pos != pMol->GetBondAngleList().end();++pos)
1052  {
1053  asigma = (*pos)->GetAngleSigma();
1054  adelta = (*pos)->GetAngleDelta();
1055  (*pos)->SetAngleDelta(0.2*DEG2RAD);
1056  (*pos)->SetAngleSigma(0.01*DEG2RAD);
1057  }
1058  }
1059  } catch (const std::bad_cast& e){
1060 
1061  }
1062  }
1063  }
1064 
1065  const long starting_point=mRefParList.CreateParamSet("MonteCarloObj:Last parameters (RANDOM-LSQ)");
1066  mRefParList.SaveParamSet(starting_point);
1067  while(nbCycle!=0) {
1068  nbCycle--;
1069  mRefParList.RestoreParamSet(starting_point);
1070  //this->NewConfiguration();
1071  for(int i=0;i<mRefinedObjList.GetNb();i++) mRefinedObjList.GetObj(i).RandomizeConfiguration();
1072  this->UpdateDisplay();
1073 
1074  //perform LSQ
1075  for(int i=0;i<mRefinedObjList.GetNb();i++) mRefinedObjList.GetObj(i).SetApproximationFlag(false);
1076  //mCurrentCost=this->GetLogLikelihood();
1077  try {
1078  mLSQ.Refine(20,true,true,false,0.001);
1079  }
1080  catch(const ObjCrystException &except) {
1081  //cout<<"Something wrong?"<<endl;
1082  };
1083  for(int i=0;i<mRefinedObjList.GetNb();i++) mRefinedObjList.GetObj(i).SetApproximationFlag(true);
1084  //cout<<"LSQ cost: "<<mCurrentCost<<" -> "<<this->GetLogLikelihood()<<endl;
1085  REAL lsq_cost=this->GetLogLikelihood();
1086  mCurrentCost = lsq_cost;
1087  //mRefParList.SaveParamSet(lsqtParSavedSetIndex);
1088  if(mCurrentCost<mBestCost)
1089  {
1090  mBestCost=mCurrentCost;
1091  mRefParList.SaveParamSet(mBestParSavedSetIndex);
1092  }
1093  this->UpdateDisplay();
1094 
1095  //save it to the file
1096  string saveFileName=this->GetName();
1097  time_t date=time(0);
1098  char strDate[40];
1099  strftime(strDate,sizeof(strDate),"%Y-%m-%d_%H-%M-%S",localtime(&date));//%Y-%m-%dT%H:%M:%S%Z
1100  char costAsChar[30];
1101  sprintf(costAsChar,"#Run%ld-Cost-%f",nbCycle, mCurrentCost);
1102  saveFileName=saveFileName+(string)strDate+(string)costAsChar+(string)".xml";
1103  XMLCrystFileSaveGlobal(saveFileName);
1104 
1105  #ifdef __WX__CRYST__
1106  mMutexStopAfterCycle.Lock();
1107  #endif
1108  if(mStopAfterCycle)
1109  {
1110  #ifdef __WX__CRYST__
1111  mMutexStopAfterCycle.Unlock();
1112  #endif
1113  break;
1114  }
1115  #ifdef __WX__CRYST__
1116  mMutexStopAfterCycle.Unlock();
1117  #endif
1118  }
1119 
1120  if(bsigma<0 || bdelta<0 || asigma<0 || adelta<0) return;
1121  //restore the delta and sigma values
1122  for(int i=0;i<mRefinedObjList.GetNb();i++) {
1123  if(mRefinedObjList.GetObj(i).GetClassName()=="Crystal") {
1124  try {
1125  Crystal * pCryst = dynamic_cast<Crystal *>(&(mRefinedObjList.GetObj(i)));
1126  for(int s=0;s<pCryst->GetScattererRegistry().GetNb();s++)
1127  {
1128  Molecule *pMol=dynamic_cast<Molecule*>(&(pCryst->GetScatt(s)));
1129  if(pMol==NULL) continue; // not a Molecule
1130  for(vector<MolBond*>::iterator pos = pMol->GetBondList().begin(); pos != pMol->GetBondList().end();++pos) {
1131  (*pos)->SetLengthDelta(bdelta);
1132  (*pos)->SetLengthSigma(bsigma);
1133  }
1134  for(vector<MolBondAngle*>::iterator pos=pMol->GetBondAngleList().begin();pos != pMol->GetBondAngleList().end();++pos)
1135  {
1136  (*pos)->SetAngleDelta(adelta);
1137  (*pos)->SetAngleSigma(asigma);
1138  }
1139  }
1140  } catch (const std::bad_cast& e){
1141 
1142  }
1143  }
1144  }
1145 }
1146 
1147 void MonteCarloObj::RunParallelTempering(long &nbStep,const bool silent,
1148  const REAL finalcost,const REAL maxTime)
1149 {
1150  TAU_PROFILE("MonteCarloObj::RunParallelTempering()","void ()",TAU_DEFAULT);
1151  TAU_PROFILE_TIMER(timer0a,"MonteCarloObj::RunParallelTempering() Begin 1","", TAU_FIELD);
1152  TAU_PROFILE_TIMER(timer0b,"MonteCarloObj::RunParallelTempering() Begin 2","", TAU_FIELD);
1153  TAU_PROFILE_TIMER(timer1,"MonteCarloObj::RunParallelTempering() New Config + LLK","", TAU_FIELD);
1154  TAU_PROFILE_TIMER(timerN,"MonteCarloObj::RunParallelTempering() Finish","", TAU_FIELD);
1155  TAU_PROFILE_START(timer0a);
1156  //Keep a copy of the total number of steps, and decrement nbStep
1157  const long nbSteps=nbStep;
1158  unsigned int accept;// 1 if last trial was accepted? 2 if new best config ? else 0
1159  mNbTrial=0;
1160  // time (in seconds) when last autoSave was made (if enabled)
1161  unsigned long secondsWhenAutoSave=0;
1162 
1163  // Periodicity of the automatic LSQ refinements (if the option is set)
1164  const unsigned int autoLSQPeriod=150000;
1165 
1166  if(!silent) cout << "Starting Parallel Tempering Optimization"<<endl;
1167  //Total number of parallel refinements,each is a 'World'. The most stable
1168  // world must be i=nbWorld-1, and the most changing World (high mutation,
1169  // high temperature) is i=0.
1170  const long nbWorld=30;
1171  CrystVector_long worldSwapIndex(nbWorld);
1172  for(int i=0;i<nbWorld;++i) worldSwapIndex(i)=i;
1173  // Number of successive trials for each World. At the end of these trials
1174  // a swap is tried with the upper World (eg i-1). This number effectvely sets
1175  // the rate of swapping.
1176  const int nbTryPerWorld=10;
1177  // Initialize the costs
1178  mCurrentCost=this->GetLogLikelihood();
1179  REAL runBestCost=mCurrentCost;
1180  CrystVector_REAL currentCost(nbWorld);
1181  currentCost=mCurrentCost;
1182  // Init the different temperatures
1183  CrystVector_REAL simAnnealTemp(nbWorld);
1184  for(int i=0;i<nbWorld;i++)
1185  {
1186  switch(mAnnealingScheduleTemp.GetChoice())
1187  {
1188  case ANNEALING_BOLTZMANN:
1189  simAnnealTemp(i)=
1190  mTemperatureMin*log((REAL)nbWorld)/log((REAL)(i+2));break;
1191  case ANNEALING_CAUCHY:
1192  simAnnealTemp(i)=mTemperatureMin*nbWorld/(i+1);break;
1193  //case ANNEALING_QUENCHING:
1194  case ANNEALING_EXPONENTIAL:
1195  simAnnealTemp(i)=mTemperatureMax
1196  *pow(mTemperatureMin/mTemperatureMax,
1197  i/(REAL)(nbWorld-1));break;
1198  case ANNEALING_GAMMA:
1199  simAnnealTemp(i)=mTemperatureMax+(mTemperatureMin-mTemperatureMax)
1200  *pow(i/(REAL)(nbWorld-1),mTemperatureGamma);break;
1201  case ANNEALING_SMART:
1202  simAnnealTemp(i)=mCurrentCost/(100.+(REAL)i/(REAL)nbWorld*900.);break;
1203  default:
1204  simAnnealTemp(i)=mCurrentCost/(100.+(REAL)i/(REAL)nbWorld*900.);break;
1205  }
1206  }
1207  //Init the different mutation rate parameters
1208  CrystVector_REAL mutationAmplitude(nbWorld);
1209  for(int i=0;i<nbWorld;i++)
1210  {
1211  switch(mAnnealingScheduleMutation.GetChoice())
1212  {
1213  case ANNEALING_BOLTZMANN:
1214  mutationAmplitude(i)=
1215  mMutationAmplitudeMin*log((REAL)(nbWorld-1))/log((REAL)(i+2));
1216  break;
1217  case ANNEALING_CAUCHY:
1218  mutationAmplitude(i)=mMutationAmplitudeMin*(REAL)(nbWorld-1)/(i+1);break;
1219  //case ANNEALING_QUENCHING:
1220  case ANNEALING_EXPONENTIAL:
1221  mutationAmplitude(i)=mMutationAmplitudeMax
1222  *pow(mMutationAmplitudeMin/mMutationAmplitudeMax,
1223  i/(REAL)(nbWorld-1));break;
1224  case ANNEALING_GAMMA:
1225  mutationAmplitude(i)=mMutationAmplitudeMax+(mMutationAmplitudeMin-mMutationAmplitudeMax)
1226  *pow(i/(REAL)(nbWorld-1),mMutationAmplitudeGamma);break;
1227  case ANNEALING_SMART:
1228  mutationAmplitude(i)=sqrt(mMutationAmplitudeMin*mMutationAmplitudeMax);break;
1229  default:
1230  mutationAmplitude(i)=sqrt(mMutationAmplitudeMin*mMutationAmplitudeMax);break;
1231  }
1232  }
1233  // Init the parameter sets for each World
1234  // All Worlds start from the same (current) configuration.
1235  CrystVector_long worldCurrentSetIndex(nbWorld);
1236  for(int i=nbWorld-1;i>=0;i--)
1237  {
1238  if((i!=(nbWorld-1))&&(i%2==0))
1239  for(int j=0;j<mRecursiveRefinedObjList.GetNb();j++)
1240  mRecursiveRefinedObjList.GetObj(j).RandomizeConfiguration();
1241  worldCurrentSetIndex(i)=mRefParList.CreateParamSet();
1242  mRefParList.RestoreParamSet(worldCurrentSetIndex(nbWorld-1));
1243  }
1244  TAU_PROFILE_STOP(timer0a);
1245  TAU_PROFILE_START(timer0b);
1246  //mNbTrial=nbSteps;;
1247  const long lastParSavedSetIndex=mRefParList.CreateParamSet("MonteCarloObj:Last parameters (PT)");
1248  const long runBestIndex=mRefParList.CreateParamSet("Best parameters for current run (PT)");
1249  CrystVector_REAL swapPar;
1250  //Keep track of how many trials are accepted for each World
1251  CrystVector_long worldNbAcceptedMoves(nbWorld);
1252  worldNbAcceptedMoves=0;
1253  //Do a report each... And check if mutation rate is OK (for annealing_smart)s
1254  const int nbTrialsReport=3000;
1255  // TEST : allow GENETIC mating of configurations
1256  //Get gene groups list :TODO: check for missing groups
1257  CrystVector_uint refParGeneGroupIndex(mRefParList.GetNbPar());
1258  unsigned int first=1;
1259  for(int i=0;i<mRecursiveRefinedObjList.GetNb();i++)
1260  mRecursiveRefinedObjList.GetObj(i).GetGeneGroup(mRefParList,refParGeneGroupIndex,first);
1261  #if 0
1262  if(!silent)
1263  for(int i=0;i<mRefParList.GetNbPar();i++)
1264  {
1265  cout << "Gene Group:"<<refParGeneGroupIndex(i)<<" :";
1266  mRefParList.GetPar(i).Print();
1267  }
1268  #endif
1269  // number of gene groups
1270  // to select which gene groups are exchanged in the mating
1271  //const unsigned int nbGeneGroup=refParGeneGroupIndex.max();
1272  //CrystVector_int crossoverGroupIndex(nbGeneGroup);
1273  //const long parSetOffspringA=mRefParList.CreateParamSet("Offspring A");
1274  //const long parSetOffspringB=mRefParList.CreateParamSet("Offspring B");
1275  // record the statistical distribution n=f(cost function) for each World
1276  //CrystMatrix_REAL trialsDensity(100,nbWorld+1);
1277  //trialsDensity=0;
1278  //for(int i=0;i<100;i++) trialsDensity(i,0)=i/(float)100;
1279  // Do we need to update the display ?
1280  bool needUpdateDisplay=false;
1281  //Do the refinement
1282  bool makeReport=false;
1283  Chronometer chrono;
1284  chrono.start();
1285  float lastUpdateDisplayTime=chrono.seconds();
1286  TAU_PROFILE_STOP(timer0b);
1287  for(;mNbTrial<nbSteps;)
1288  {
1289  for(int i=0;i<nbWorld;i++)
1290  {
1291  mContext=i;
1292  //mRefParList.RestoreParamSet(worldCurrentSetIndex(i));
1293  mMutationAmplitude=mutationAmplitude(i);
1294  mTemperature=simAnnealTemp(i);
1295  for(int j=0;j<nbTryPerWorld;j++)
1296  {
1297  //mRefParList.SaveParamSet(lastParSavedSetIndex);
1298  TAU_PROFILE_START(timer1);
1299  mRefParList.RestoreParamSet(worldCurrentSetIndex(i));
1300  this->NewConfiguration();
1301  accept=0;
1302  REAL cost=this->GetLogLikelihood();
1303  TAU_PROFILE_STOP(timer1);
1304  //trialsDensity((long)(cost*100.),i+1)+=1;
1305  if(cost<currentCost(i))
1306  {
1307  accept=1;
1308  currentCost(i)=cost;
1309  mRefParList.SaveParamSet(worldCurrentSetIndex(i));
1310  if(cost<runBestCost)
1311  {
1312  accept=2;
1313  runBestCost=currentCost(i);
1314  this->TagNewBestConfig();
1315  needUpdateDisplay=true;
1316 
1317  mRefParList.SaveParamSet(runBestIndex);
1318  if(runBestCost<mBestCost)
1319  {
1320  mBestCost=currentCost(i);
1321  mRefParList.SaveParamSet(mBestParSavedSetIndex);
1322  if(!silent) cout << "->Trial :" << mNbTrial
1323  << " World="<< worldSwapIndex(i)
1324  << " Temp="<< mTemperature
1325  << " Mutation Ampl.: "<<mMutationAmplitude
1326  << " NEW OVERALL Best Cost="<<mBestCost<< endl;
1327  }
1328  else if(!silent) cout << "->Trial :" << mNbTrial
1329  << " World="<< worldSwapIndex(i)
1330  << " Temp="<< mTemperature
1331  << " Mutation Ampl.: "<<mMutationAmplitude
1332  << " NEW RUN Best Cost="<<runBestCost<< endl;
1333  if(!silent) this->DisplayReport();
1334  }
1335  worldNbAcceptedMoves(i)++;
1336  }
1337  else
1338  {
1339  if(log((rand()+1)/(REAL)RAND_MAX)<(-(cost-currentCost(i))/mTemperature) )
1340  {
1341  accept=1;
1342  currentCost(i)=cost;
1343  mRefParList.SaveParamSet(worldCurrentSetIndex(i));
1344  worldNbAcceptedMoves(i)++;
1345  }
1346  }
1347  //if(accept==1 && i==(nbWorld-1)){this->UpdateDisplay();}
1348  if( ((mXMLAutoSave.GetChoice()==1)&&((chrono.seconds()-secondsWhenAutoSave)>86400))
1349  ||((mXMLAutoSave.GetChoice()==2)&&((chrono.seconds()-secondsWhenAutoSave)>3600))
1350  ||((mXMLAutoSave.GetChoice()==3)&&((chrono.seconds()-secondsWhenAutoSave)> 600))
1351  ||((mXMLAutoSave.GetChoice()==4)&&(accept==2)) )
1352  {
1353  secondsWhenAutoSave=(unsigned long)chrono.seconds();
1354  string saveFileName=this->GetName();
1355  time_t date=time(0);
1356  char strDate[40];
1357  strftime(strDate,sizeof(strDate),"%Y-%m-%d_%H-%M-%S",localtime(&date));//%Y-%m-%dT%H:%M:%S%Z
1358  char costAsChar[30];
1359  if(accept!=2) mRefParList.RestoreParamSet(mBestParSavedSetIndex);
1360  sprintf(costAsChar,"-Cost-%f",this->GetLogLikelihood());
1361  saveFileName=saveFileName+(string)strDate+(string)costAsChar+(string)".xml";
1362  XMLCrystFileSaveGlobal(saveFileName);
1363  //if(accept!=2) mRefParList.RestoreParamSet(lastParSavedSetIndex);
1364  }
1365  //if(accept==0) mRefParList.RestoreParamSet(lastParSavedSetIndex);
1366  mNbTrial++;nbStep--;
1367  if((mNbTrial%nbTrialsReport)==0) makeReport=true;
1368  }//nbTryPerWorld trials
1369  }//For each World
1370 
1371  if(mAutoLSQ.GetChoice()==2)
1372  if((mNbTrial%autoLSQPeriod)<(nbTryPerWorld*nbWorld))
1373  {// Try a quick LSQ ?
1374  for(int i=0;i<mRefinedObjList.GetNb();i++) mRefinedObjList.GetObj(i).SetApproximationFlag(false);
1375  for(int i=nbWorld-5;i<nbWorld;i++)
1376  {
1377  #ifdef __WX__CRYST__
1378  mMutexStopAfterCycle.Lock();
1379  if(mStopAfterCycle)
1380  {
1381  mMutexStopAfterCycle.Unlock();
1382  break;
1383  }
1384  mMutexStopAfterCycle.Unlock();
1385  #endif
1386 
1387  mRefParList.RestoreParamSet(worldCurrentSetIndex(i));
1388 
1389  #if 0
1390  // Report GoF values (Chi^2 / nbObs) values for all objects
1391  for(map<RefinableObj*,unsigned int>::iterator pos=mLSQ.GetRefinedObjMap().begin();pos!=mLSQ.GetRefinedObjMap().end();++pos)
1392  if(pos->first->GetNbLSQFunction()>0)
1393  {
1394  CrystVector_REAL tmp;
1395  tmp =pos->first->GetLSQCalc(pos->second);
1396  tmp-=pos->first->GetLSQObs (pos->second);
1397  tmp*=tmp;
1398  tmp*=pos->first->GetLSQWeight(pos->second);
1399  cout<<pos->first->GetClassName()<<":"<<pos->first->GetName()<<": GoF="<<tmp.sum()/tmp.numElements();
1400  }
1401  cout<<endl;
1402  #endif
1403 
1404  const REAL cost0=this->GetLogLikelihood();// cannot use currentCost(i), approximations changed...
1405  if(!silent) cout<<"LSQ: World="<<worldSwapIndex(i)<<": cost="<<cost0;
1406  try {mLSQ.Refine(-30,true,true,false,0.001);}
1407  catch(const ObjCrystException &except){};
1408  #if 0
1409  // Report GoF values (Chi^2 / nbObs) values for all objects
1410  for(map<RefinableObj*,unsigned int>::iterator pos=mLSQ.GetRefinedObjMap().begin();pos!=mLSQ.GetRefinedObjMap().end();++pos)
1411  if(pos->first->GetNbLSQFunction()>0)
1412  {
1413  CrystVector_REAL tmp;
1414  tmp =pos->first->GetLSQCalc(pos->second);
1415  tmp-=pos->first->GetLSQObs (pos->second);
1416  tmp*=tmp;
1417  tmp*=pos->first->GetLSQWeight(pos->second);
1418  cout<<pos->first->GetClassName()<<":"<<pos->first->GetName()<<": GoF="<<tmp.sum()/tmp.numElements();
1419  }
1420  cout<<endl;
1421  #endif
1422  const REAL cost=this->GetLogLikelihood();
1423  if(!silent) cout<<" -> "<<cost<<endl;
1424  if(cost<cost0) mRefParList.SaveParamSet(worldCurrentSetIndex(i));
1425  }
1426  // Need to go back to optimization with approximations allowed (they are not during LSQ)
1427  for(int i=0;i<mRefinedObjList.GetNb();i++) mRefinedObjList.GetObj(i).SetApproximationFlag(true);
1428  // And recompute LLK - since they will be lower
1429  for(int i=nbWorld-5;i<nbWorld;i++)
1430  {
1431  mRefParList.RestoreParamSet(worldCurrentSetIndex(i));
1432  const REAL cost=this->GetLogLikelihood();
1433  if(!silent) cout<<"LSQ2:"<<currentCost(i)<<"->"<<cost<<endl;
1434  if(cost<currentCost(i))
1435  {
1436  const REAL oldcost=currentCost(i);
1437  mRefParList.SaveParamSet(worldCurrentSetIndex(i));
1438  currentCost(i)=cost;
1439  if(cost<runBestCost)
1440  {
1441  runBestCost=currentCost(i);
1442  this->TagNewBestConfig();
1443  needUpdateDisplay=true;
1444 
1445  mRefParList.SaveParamSet(runBestIndex);
1446  if(runBestCost<mBestCost)
1447  {
1448  mBestCost=currentCost(i);
1449  mRefParList.SaveParamSet(mBestParSavedSetIndex);
1450  if(!silent) cout << "->Trial :" << mNbTrial
1451  << " World="<< worldSwapIndex(i)
1452  << " LSQ2: NEW OVERALL Best Cost="<<mBestCost<< endl;
1453  }
1454  else if(!silent) cout << "->Trial :" << mNbTrial
1455  << " World="<< worldSwapIndex(i)
1456  << " LSQ2: NEW RUN Best Cost="<<runBestCost<< endl;
1457  if(!silent) this->DisplayReport();
1458  }
1459  // KLUDGE - after a successful LSQ, we will be close to a minimum,
1460  // which will make most successive global optimization trials to
1461  // be rejected, until the temperature is increased a lot - this
1462  // is a problem as the temperature increases so much that the
1463  // benefit of the LSQ is essentially negated.
1464  // So we need to use a higher recorded cost, so that successive trials
1465  // may be accepted
1466  #if 0
1467  mMutationAmplitude=mutationAmplitude(i);
1468  for(unsigned int ii=0;ii<4;ii++) this->NewConfiguration(gpRefParTypeObjCryst,false);
1469  currentCost(i)=(this->GetLogLikelihood()+cost)/2;
1470  if(!silent) cout<<"LSQ3: #"<<worldSwapIndex(i)<<":"<<cost<<"->"<<currentCost(i)<<endl;
1471  #else
1472  currentCost(i)=oldcost;
1473  #endif
1474  }
1475  }
1476  }
1477 
1478  //Try swapping worlds
1479  for(int i=1;i<nbWorld;i++)
1480  {
1481  #if 0
1482  mRefParList.RestoreParamSet(worldCurrentSetIndex(i));
1483  mMutationAmplitude=mutationAmplitude(i);
1484  cout<<i<<":"<<currentCost(i)<<":"<<this->GetLogLikelihood()<<endl;
1485  #endif
1486  #if 1
1487  if( log((rand()+1)/(REAL)RAND_MAX)
1488  < (-(currentCost(i-1)-currentCost(i))/simAnnealTemp(i)))
1489  #else
1490  // Compare World (i-1) and World (i) with the same amplitude,
1491  // hence the same max likelihood error
1492  mRefParList.RestoreParamSet(worldCurrentSetIndex(i-1));
1493  mMutationAmplitude=mutationAmplitude(i);
1494  if( log((rand()+1)/(REAL)RAND_MAX)
1495  < (-(this->GetLogLikelihood()-currentCost(i))/simAnnealTemp(i)))
1496  #endif
1497  {
1498  /*
1499  if(i>2)
1500  {
1501  cout <<"->Swapping Worlds :" << i <<"(cost="<<currentCost(i)<<")"
1502  <<" with "<< (i-1) <<"(cost="<< currentCost(i-1)<<")"<<endl;
1503  }
1504  */
1505  swapPar=mRefParList.GetParamSet(worldCurrentSetIndex(i));
1506  mRefParList.GetParamSet(worldCurrentSetIndex(i))=
1507  mRefParList.GetParamSet(worldCurrentSetIndex(i-1));
1508  mRefParList.GetParamSet(worldCurrentSetIndex(i-1))=swapPar;
1509  const REAL tmp=currentCost(i);
1510  currentCost(i)=currentCost(i-1);
1511  currentCost(i-1)=tmp;
1512  const long tmpIndex=worldSwapIndex(i);
1513  worldSwapIndex(i)=worldSwapIndex(i-1);
1514  worldSwapIndex(i-1)=tmpIndex;
1515  #if 0
1516  // Compute correct costs in the case we use maximum likelihood
1517  mRefParList.RestoreParamSet(worldCurrentSetIndex(i));
1518  mMutationAmplitude=mutationAmplitude(i);
1519  currentCost(i)=this->GetLogLikelihood();
1520 
1521  mRefParList.RestoreParamSet(worldCurrentSetIndex(i-1));
1522  mMutationAmplitude=mutationAmplitude(i-1);
1523  currentCost(i-1)=this->GetLogLikelihood();
1524  #endif
1525  }
1526  }
1527  #if 0
1528  //Try mating worlds- NEW !
1529  TAU_PROFILE_TIMER(timer1,\
1530  "MonteCarloObj::Optimize (Try mating Worlds)"\
1531  ,"", TAU_FIELD);
1532  TAU_PROFILE_START(timer1);
1533  if( (rand()/(REAL)RAND_MAX)<.1)
1534  for(int k=nbWorld-1;k>nbWorld/2;k--)
1535  for(int i=k-nbWorld/3;i<k;i++)
1536  {
1537  #if 0
1538  // Random switching of gene groups
1539  for(unsigned int j=0;j<nbGeneGroup;j++)
1540  crossoverGroupIndex(j)= (int) floor(rand()/((REAL)RAND_MAX-1)*2);
1541  for(int j=0;j<mRefParList.GetNbPar();j++)
1542  {
1543  if(0==crossoverGroupIndex(refParGeneGroupIndex(j)-1))
1544  {
1545  mRefParList.GetParamSet(parSetOffspringA)(j)=
1546  mRefParList.GetParamSet(worldCurrentSetIndex(i))(j);
1547  mRefParList.GetParamSet(parSetOffspringB)(j)=
1548  mRefParList.GetParamSet(worldCurrentSetIndex(k))(j);
1549  }
1550  else
1551  {
1552  mRefParList.GetParamSet(parSetOffspringA)(j)=
1553  mRefParList.GetParamSet(worldCurrentSetIndex(k))(j);
1554  mRefParList.GetParamSet(parSetOffspringB)(j)=
1555  mRefParList.GetParamSet(worldCurrentSetIndex(i))(j);
1556  }
1557  }
1558  #endif
1559  #if 1
1560  // Switch gene groups in two parts
1561  unsigned int crossoverPoint1=
1562  (int)(1+floor(rand()/((REAL)RAND_MAX-1)*(nbGeneGroup)));
1563  unsigned int crossoverPoint2=
1564  (int)(1+floor(rand()/((REAL)RAND_MAX-1)*(nbGeneGroup)));
1565  if(crossoverPoint2<crossoverPoint1)
1566  {
1567  int tmp=crossoverPoint1;
1568  crossoverPoint1=crossoverPoint2;
1569  crossoverPoint2=tmp;
1570  }
1571  if(crossoverPoint1==crossoverPoint2) crossoverPoint2+=1;
1572  for(int j=0;j<mRefParList.GetNbPar();j++)
1573  {
1574  if((refParGeneGroupIndex(j)>crossoverPoint1)&&refParGeneGroupIndex(j)<crossoverPoint2)
1575  {
1576  mRefParList.GetParamSet(parSetOffspringA)(j)=
1577  mRefParList.GetParamSet(worldCurrentSetIndex(i))(j);
1578  mRefParList.GetParamSet(parSetOffspringB)(j)=
1579  mRefParList.GetParamSet(worldCurrentSetIndex(k))(j);
1580  }
1581  else
1582  {
1583  mRefParList.GetParamSet(parSetOffspringA)(j)=
1584  mRefParList.GetParamSet(worldCurrentSetIndex(k))(j);
1585  mRefParList.GetParamSet(parSetOffspringB)(j)=
1586  mRefParList.GetParamSet(worldCurrentSetIndex(i))(j);
1587  }
1588  }
1589  #endif
1590  // Try both offspring
1591  for(int junk=0;junk<2;junk++)
1592  {
1593  if(junk==0) mRefParList.RestoreParamSet(parSetOffspringA);
1594  else mRefParList.RestoreParamSet(parSetOffspringB);
1595  REAL cost=this->GetLogLikelihood();
1596  //if(log((rand()+1)/(REAL)RAND_MAX)
1597  // < (-(cost-currentCost(k))/simAnnealTemp(k)))
1598  if(cost<currentCost(k))
1599  {
1600  // Also exchange genes for higher-temperature World ?
1601  //if(junk==0)
1602  // mRefParList.GetParamSet(worldCurrentSetIndex(i))=
1603  // mRefParList.GetParamSet(parSetOffspringB);
1604  //else
1605  // mRefParList.GetParamSet(worldCurrentSetIndex(i))=
1606  // mRefParList.GetParamSet(parSetOffspringA);
1607  currentCost(k)=cost;
1608  mRefParList.SaveParamSet(worldCurrentSetIndex(k));
1609  //worldNbAcceptedMoves(k)++;
1610  if(!silent) cout << "Accepted mating :"<<k<<"(with"<<i<<")"
1611  <<" (crossoverGene1="<< crossoverPoint1<<","
1612  <<" crossoverGene2="<< crossoverPoint2<<")"
1613  <<endl;
1614  if(cost<runBestCost)
1615  {
1616  runBestCost=cost;
1617  this->TagNewBestConfig();
1618  needUpdateDisplay=true;
1619  mRefParList.SaveParamSet(runBestIndex);
1620  if(cost<mBestCost)
1621  {
1622  mBestCost=cost;
1623  mRefParList.SaveParamSet(mBestParSavedSetIndex);
1624  if(!silent) cout << "->Trial :" << mNbTrial
1625  << " World="<< worldSwapIndex(k)
1626  << " Temp="<< simAnnealTemp(k)
1627  << " Mutation Ampl.: "<<mMutationAmplitude
1628  << " NEW OVERALL Best Cost="<<mBestCost<< "(MATING !)"<<endl;
1629  }
1630  else if(!silent) cout << "->Trial :" << mNbTrial
1631  << " World="<< worldSwapIndex(k)
1632  << " Temp="<< simAnnealTemp(k)
1633  << " Mutation Ampl.: "<<mMutationAmplitude
1634  << " NEW RUN Best Cost="<<runBestCost<< "(MATING !)"<<endl;
1635  bestConfigNb=mNbTrial;
1636  if(!silent) this->DisplayReport();
1637  //for(int i=0;i<mRefinedObjList.GetNb();i++)
1638  // mRefinedObjList.GetObj(i).Print();
1639  }
1640  i=k;//Don't test other Worlds
1641  break;
1642  }
1643  //mNbTrial++;nbStep--;
1644  //if((mNbTrial%nbTrialsReport)==0) makeReport=true;
1645  }
1646  }
1647  TAU_PROFILE_STOP(timer1);
1648  #endif
1649  if(true==makeReport)
1650  {
1651  makeReport=false;
1652  worldNbAcceptedMoves*=nbWorld;
1653  if(!silent)
1654  {
1655  #if 0
1656  {// Experimental, dynamical weighting
1657  REAL max=0.;
1658  map<const RefinableObj*,REAL> ll,llvar;
1659  map<const RefinableObj*,LogLikelihoodStats>::iterator pos;
1660  for(pos=mvContextObjStats[0].begin();pos!=mvContextObjStats[0].end();++pos)
1661  {
1662  ll [pos->first]=0.;
1663  llvar[pos->first]=0.;
1664  }
1665  for(int i=0;i<nbWorld;i++)
1666  {
1667  for(pos=mvContextObjStats[0].begin();pos!=mvContextObjStats[0].end();++pos)
1668  {
1669  ll [pos->first] += pos->second.mTotalLogLikelihood;
1670  llvar[pos->first] += pos->second.mTotalLogLikelihoodDeltaSq;
1671  }
1672  }
1673  for(pos=mvContextObjStats[0].begin();pos!=mvContextObjStats[0].end();++pos)
1674  {
1675  cout << pos->first->GetName()
1676  << " " << llvar[pos->first]
1677  << " " << mvObjWeight[pos->first].mWeight
1678  << " " << max<<endl;
1679  llvar[pos->first] *= mvObjWeight[pos->first].mWeight;
1680  if(llvar[pos->first]>max) max=llvar[pos->first];
1681  }
1682  map<const RefinableObj*,REAL>::iterator pos2;
1683  for(pos2=llvar.begin();pos2!=llvar.end();++pos2)
1684  {
1685  const REAL d=pos2->second;
1686  if(d<(max/mvObjWeight.size()/10.))
1687  {
1688  if(d<1) continue;
1689  mvObjWeight[pos2->first].mWeight *=2;
1690  }
1691  }
1692  REAL ll1=0;
1693  REAL llt=0;
1694  for(pos2=ll.begin();pos2!=ll.end();++pos2)
1695  {
1696  llt += pos2->second;
1697  ll1 += pos2->second * mvObjWeight[pos2->first].mWeight;
1698  }
1699  map<const RefinableObj*,DynamicObjWeight>::iterator posw;
1700  for(posw=mvObjWeight.begin();posw!=mvObjWeight.end();++posw)
1701  {
1702  posw->second.mWeight *= llt/ll1;
1703  }
1704  }
1705  #endif //Experimental dynamical weighting
1706  #if 1 //def __DEBUG__
1707  for(int i=0;i<nbWorld;i++)
1708  {
1709  cout<<" World :"<<worldSwapIndex(i)<<":";
1710  map<const RefinableObj*,LogLikelihoodStats>::iterator pos;
1711  for(pos=mvContextObjStats[i].begin();pos!=mvContextObjStats[i].end();++pos)
1712  {
1713  cout << pos->first->GetName()
1714  << "(LLK="
1715  << pos->second.mLastLogLikelihood
1716  //<< "(<LLK>="
1717  //<< pos->second.mTotalLogLikelihood/nbTrialsReport
1718  //<< ", <delta(LLK)^2>="
1719  //<< pos->second.mTotalLogLikelihoodDeltaSq/nbTrialsReport
1720  << ", w="<<mvObjWeight[pos->first].mWeight
1721  <<") ";
1722  pos->second.mTotalLogLikelihood=0;
1723  pos->second.mTotalLogLikelihoodDeltaSq=0;
1724  }
1725  cout << endl;
1726  }
1727  #endif
1728  for(int i=0;i<nbWorld;i++)
1729  {
1730  //mRefParList.RestoreParamSet(worldCurrentSetIndex(i));
1731  cout <<" World :" << worldSwapIndex(i)
1732  <<" Temp.: " << simAnnealTemp(i)
1733  <<" Mutation Ampl.: " << mutationAmplitude(i)
1734  <<" Current Cost=" << currentCost(i)
1735  <<" Accepting "
1736  << (int)((REAL)worldNbAcceptedMoves(i)/nbTrialsReport*100)
1737  <<"% moves " <<endl;
1738  // <<"% moves " << mRefParList.GetPar("Pboccup").GetValue()<<endl;
1739  }
1740  }
1741  if(!silent) cout <<"Trial :" << mNbTrial << " Best Cost=" << runBestCost<< " ";
1742  if(!silent) chrono.print();
1743  //Change the mutation rate if necessary for each world
1744  if(ANNEALING_SMART==mAnnealingScheduleMutation.GetChoice())
1745  {
1746  for(int i=0;i<nbWorld;i++)
1747  {
1748  if((worldNbAcceptedMoves(i)/(REAL)nbTrialsReport)>0.30)
1749  mutationAmplitude(i)*=2.;
1750  if((worldNbAcceptedMoves(i)/(REAL)nbTrialsReport)<0.10)
1751  mutationAmplitude(i)/=2.;
1752  if(mutationAmplitude(i)>mMutationAmplitudeMax)
1753  mutationAmplitude(i)=mMutationAmplitudeMax;
1754  if(mutationAmplitude(i)<mMutationAmplitudeMin)
1755  mutationAmplitude(i)=mMutationAmplitudeMin;
1756  }
1757  }
1758  if(ANNEALING_SMART==mAnnealingScheduleTemp.GetChoice())
1759  {
1760  for(int i=0;i<nbWorld;i++)
1761  {
1762  if((worldNbAcceptedMoves(i)/(REAL)nbTrialsReport)>0.30)
1763  simAnnealTemp(i)/=1.5;
1764  if((worldNbAcceptedMoves(i)/(REAL)nbTrialsReport)>0.80)
1765  simAnnealTemp(i)/=1.5;
1766  if((worldNbAcceptedMoves(i)/(REAL)nbTrialsReport)>0.95)
1767  simAnnealTemp(i)/=1.5;
1768 
1769  if((worldNbAcceptedMoves(i)/(REAL)nbTrialsReport)<0.10)
1770  simAnnealTemp(i)*=1.5;
1771  if((worldNbAcceptedMoves(i)/(REAL)nbTrialsReport)<0.04)
1772  simAnnealTemp(i)*=1.5;
1773  //if((worldNbAcceptedMoves(i)/(REAL)nbTrialsReport)<0.01)
1774  // simAnnealTemp(i)*=1.5;
1775  //cout<<"World#"<<i<<":"<<worldNbAcceptedMoves(i)<<":"<<nbTrialsReport<<endl;
1776  //if(simAnnealTemp(i)>mTemperatureMax) simAnnealTemp(i)=mTemperatureMax;
1777  //if(simAnnealTemp(i)<mTemperatureMin) simAnnealTemp(i)=mTemperatureMin;
1778  }
1779  }
1780  worldNbAcceptedMoves=0;
1781  //this->DisplayReport();
1782 
1783  #ifdef __WX__CRYST__
1784  if(0!=mpWXCrystObj) mpWXCrystObj->UpdateDisplayNbTrial();
1785  #endif
1786  }
1787  if( (needUpdateDisplay&&(lastUpdateDisplayTime<(chrono.seconds()-1)))||(lastUpdateDisplayTime<(chrono.seconds()-10)))
1788  {
1789  mRefParList.RestoreParamSet(runBestIndex);
1790  this->UpdateDisplay();
1791  needUpdateDisplay=false;
1792  lastUpdateDisplayTime=chrono.seconds();
1793  }
1794  #ifdef __WX__CRYST__
1795  mMutexStopAfterCycle.Lock();
1796  #endif
1797  if((runBestCost<finalcost) || mStopAfterCycle ||( (maxTime>0)&&(chrono.seconds()>maxTime)))
1798  {
1799  #ifdef __WX__CRYST__
1800  mMutexStopAfterCycle.Unlock();
1801  #endif
1802  if(!silent) cout << endl <<endl << "Refinement Stopped:"<<mBestCost<<endl;
1803  break;
1804  }
1805  #ifdef __WX__CRYST__
1806  mMutexStopAfterCycle.Unlock();
1807  #endif
1808  }//Trials
1809 
1810  TAU_PROFILE_START(timerN);
1811  if(mAutoLSQ.GetChoice()>0)
1812  {// LSQ
1813  if(!silent) cout<<"Beginning final LSQ refinement"<<endl;
1814  for(int i=0;i<mRefinedObjList.GetNb();i++) mRefinedObjList.GetObj(i).SetApproximationFlag(false);
1815  mRefParList.RestoreParamSet(runBestIndex);
1816  mCurrentCost=this->GetLogLikelihood();
1817  try {mLSQ.Refine(-50,true,true,false,0.001);}
1818  catch(const ObjCrystException &except){};
1819  if(!silent) cout<<"LSQ cost: "<<mCurrentCost<<" -> "<<this->GetLogLikelihood()<<endl;
1820 
1821  // Need to go back to optimization with approximations allowed (they are not during LSQ)
1822  for(int i=0;i<mRefinedObjList.GetNb();i++) mRefinedObjList.GetObj(i).SetApproximationFlag(true);
1823 
1824  REAL cost=this->GetLogLikelihood();
1825  if(cost<mCurrentCost)
1826  {
1827  mCurrentCost=cost;
1828  mRefParList.SaveParamSet(lastParSavedSetIndex);
1829  if(mCurrentCost<runBestCost)
1830  {
1831  runBestCost=mCurrentCost;
1832  mRefParList.SaveParamSet(runBestIndex);
1833  if(runBestCost<mBestCost)
1834  {
1835  mBestCost=mCurrentCost;
1836  mRefParList.SaveParamSet(mBestParSavedSetIndex);
1837  if(!silent) cout << "LSQ : NEW OVERALL Best Cost="<<runBestCost<< endl;
1838  }
1839  else if(!silent) cout << " LSQ : NEW Run Best Cost="<<runBestCost<< endl;
1840  }
1841  }
1842  if(!silent) cout<<"Finished LSQ refinement"<<endl;
1843  }
1844 
1845  mLastOptimTime=chrono.seconds();
1846  //Restore Best values
1847  //mRefParList.Print();
1848  if(!silent) this->DisplayReport();
1849  mRefParList.RestoreParamSet(runBestIndex);
1850  //for(int i=0;i<mRefinedObjList.GetNb();i++) mRefinedObjList.GetObj(i).Print();
1851  mCurrentCost=this->GetLogLikelihood();
1852  if(!silent) cout<<"Run Best Cost:"<<mCurrentCost<<endl;
1853  if(!silent) chrono.print();
1854  //Save density of states
1855  //ofstream out("densityOfStates.txt");
1856  //out << trialsDensity<<endl;
1857  //out.close();
1858  // Clear temporary param set
1859  for(int i=0;i<nbWorld;i++)
1860  {
1861  mRefParList.ClearParamSet(worldCurrentSetIndex(i));
1862  //mvSavedParamSet.push_back(make_pair(worldCurrentSetIndex(i),currentCost(i)));
1863  }
1864  mRefParList.ClearParamSet(lastParSavedSetIndex);
1865  mRefParList.ClearParamSet(runBestIndex);
1866  TAU_PROFILE_STOP(timerN);
1867 }
1868 
1869 void MonteCarloObj::XMLOutput(ostream &os,int indent)const
1870 {
1871  VFN_DEBUG_ENTRY("MonteCarloObj::XMLOutput():"<<this->GetName(),5)
1872  for(int i=0;i<indent;i++) os << " " ;
1873  XMLCrystTag tag("GlobalOptimObj");
1874  tag.AddAttribute("Name",this->GetName());
1875  tag.AddAttribute("NbTrialPerRun",(boost::format("%d")%(this->NbTrialPerRun())).str());
1876 
1877  os <<tag<<endl;
1878  indent++;
1879 
1880  mGlobalOptimType.XMLOutput(os,indent);
1881  os<<endl;
1882 
1883  mAnnealingScheduleTemp.XMLOutput(os,indent);
1884  os<<endl;
1885 
1886  mXMLAutoSave.XMLOutput(os,indent);
1887  os<<endl;
1888 
1889  mAutoLSQ.XMLOutput(os,indent);
1890  os<<endl;
1891 
1892  {
1893  XMLCrystTag tag2("TempMaxMin");
1894  for(int i=0;i<indent;i++) os << " " ;
1895  os<<tag2<<mTemperatureMax << " "<< mTemperatureMin;
1896  tag2.SetIsEndTag(true);
1897  os<<tag2<<endl;
1898  }
1899 
1900  mAnnealingScheduleMutation.XMLOutput(os,indent);
1901  os<<endl;
1902 
1903  mSaveTrackedData.XMLOutput(os,indent);
1904  os<<endl;
1905 
1906  {
1907  XMLCrystTag tag2("MutationMaxMin");
1908  for(int i=0;i<indent;i++) os << " " ;
1909  os<<tag2<<mMutationAmplitudeMax << " "<< mMutationAmplitudeMin;
1910  tag2.SetIsEndTag(true);
1911  os<<tag2<<endl;
1912  }
1913 
1914  for(int j=0;j<mRefinedObjList.GetNb();j++)
1915  {
1916  XMLCrystTag tag2("RefinedObject",false,true);
1917  tag2.AddAttribute("ObjectType",mRefinedObjList.GetObj(j).GetClassName());
1918  tag2.AddAttribute("ObjectName",mRefinedObjList.GetObj(j).GetName());
1919  for(int i=0;i<indent;i++) os << " " ;
1920  os<<tag2<<endl;
1921  }
1922 
1923  indent--;
1924  tag.SetIsEndTag(true);
1925  for(int i=0;i<indent;i++) os << " " ;
1926  os <<tag<<endl;
1927  VFN_DEBUG_EXIT("MonteCarloObj::XMLOutput():"<<this->GetName(),5)
1928 }
1929 
1930 void MonteCarloObj::XMLInput(istream &is,const XMLCrystTag &tagg)
1931 {
1932  VFN_DEBUG_ENTRY("MonteCarloObj::XMLInput():"<<this->GetName(),5)
1933  for(unsigned int i=0;i<tagg.GetNbAttribute();i++)
1934  {
1935  if("Name"==tagg.GetAttributeName(i)) this->SetName(tagg.GetAttributeValue(i));
1936  if("NbTrialPerRun"==tagg.GetAttributeName(i))
1937  {
1938  stringstream ss(tagg.GetAttributeValue(i));
1939  long v;
1940  ss>>v;
1941  this->NbTrialPerRun()=v;
1942  }
1943  }
1944  while(true)
1945  {
1946  XMLCrystTag tag(is);
1947  if(("GlobalOptimObj"==tag.GetName())&&tag.IsEndTag())
1948  {
1949  VFN_DEBUG_EXIT("MonteCarloObj::Exit():"<<this->GetName(),5)
1950  this->UpdateDisplay();
1951  return;
1952  }
1953  if("Option"==tag.GetName())
1954  {
1955  for(unsigned int i=0;i<tag.GetNbAttribute();i++)
1956  if("Name"==tag.GetAttributeName(i))
1957  {
1958  if("Algorithm"==tag.GetAttributeValue(i))
1959  {
1960  mGlobalOptimType.XMLInput(is,tag);
1961  break;
1962  }
1963  if("Temperature Schedule"==tag.GetAttributeValue(i))
1964  {
1965  mAnnealingScheduleTemp.XMLInput(is,tag);
1966  break;
1967  }
1968  if("Displacement Amplitude Schedule"==tag.GetAttributeValue(i))
1969  {
1970  mAnnealingScheduleMutation.XMLInput(is,tag);
1971  break;
1972  }
1973  if("Save Best Config Regularly"==tag.GetAttributeValue(i))
1974  {
1975  mXMLAutoSave.XMLInput(is,tag);
1976  break;
1977  }
1978  if("Save Tracked Data"==tag.GetAttributeValue(i))
1979  {
1980  mSaveTrackedData.XMLInput(is,tag);
1981  break;
1982  }
1983  if("Automatic Least Squares Refinement"==tag.GetAttributeValue(i))
1984  {
1985  mAutoLSQ.XMLInput(is,tag);
1986  break;
1987  }
1988  }
1989  continue;
1990  }
1991  if("TempMaxMin"==tag.GetName())
1992  {
1993  is>>mTemperatureMax>> mTemperatureMin;
1994  if(false==tag.IsEmptyTag()) XMLCrystTag junk(is);//:KLUDGE: for first release
1995  continue;
1996  }
1997  if("MutationMaxMin"==tag.GetName())
1998  {
1999  is>>mMutationAmplitudeMax>>mMutationAmplitudeMin;
2000  if(false==tag.IsEmptyTag()) XMLCrystTag junk(is);//:KLUDGE: for first release
2001  continue;
2002  }
2003  if("RefinedObject"==tag.GetName())
2004  {
2005  string name,type;
2006  for(unsigned int i=0;i<tag.GetNbAttribute();i++)
2007  {
2008  if("ObjectName"==tag.GetAttributeName(i)) name=tag.GetAttributeValue(i);
2009  if("ObjectType"==tag.GetAttributeName(i)) type=tag.GetAttributeValue(i);
2010  }
2011  RefinableObj* obj=& (gRefinableObjRegistry.GetObj(name,type));
2012  this->AddRefinableObj(*obj);
2013  continue;
2014  }
2015  }
2016 }
2017 
2018 const string MonteCarloObj::GetClassName()const { return "MonteCarloObj";}
2019 
2020 LSQNumObj& MonteCarloObj::GetLSQObj() {return mLSQ;}
2021 
2022 const LSQNumObj& MonteCarloObj::GetLSQObj() const{return mLSQ;}
2023 
2024 void MonteCarloObj::NewConfiguration(const RefParType *type)
2025 {
2026  TAU_PROFILE("MonteCarloObj::NewConfiguration()","void ()",TAU_DEFAULT);
2027  VFN_DEBUG_ENTRY("MonteCarloObj::NewConfiguration()",4)
2028  for(int i=0;i<mRefinedObjList.GetNb();i++)
2029  mRefinedObjList.GetObj(i).BeginGlobalOptRandomMove();
2030  for(int i=0;i<mRefinedObjList.GetNb();i++)
2031  mRefinedObjList.GetObj(i).GlobalOptRandomMove(mMutationAmplitude,type);
2032  VFN_DEBUG_EXIT("MonteCarloObj::NewConfiguration()",4)
2033 }
2034 
2035 void MonteCarloObj::InitOptions()
2036 {
2037  VFN_DEBUG_MESSAGE("MonteCarloObj::InitOptions()",5)
2038  this->OptimizationObj::InitOptions();
2039  static string GlobalOptimTypeName;
2040  static string GlobalOptimTypeChoices[2];//:TODO: Add Genetic Algorithm
2041 
2042  static string AnnealingScheduleChoices[6];
2043 
2044  static string AnnealingScheduleTempName;
2045  static string AnnealingScheduleMutationName;
2046 
2047  static string runAutoLSQName;
2048  static string runAutoLSQChoices[3];
2049 
2050  static string saveTrackedDataName;
2051  static string saveTrackedDataChoices[2];
2052 
2053  static bool needInitNames=true;
2054  if(true==needInitNames)
2055  {
2056  GlobalOptimTypeName="Algorithm";
2057  GlobalOptimTypeChoices[0]="Simulated Annealing";
2058  GlobalOptimTypeChoices[1]="Parallel Tempering";
2059  //GlobalOptimTypeChoices[2]="Random-LSQ";
2060 
2061  AnnealingScheduleTempName="Temperature Schedule";
2062  AnnealingScheduleMutationName="Displacement Amplitude Schedule";
2063  AnnealingScheduleChoices[0]="Constant";
2064  AnnealingScheduleChoices[1]="Boltzmann";
2065  AnnealingScheduleChoices[2]="Cauchy";
2066  AnnealingScheduleChoices[3]="Exponential";
2067  AnnealingScheduleChoices[4]="Smart";
2068  AnnealingScheduleChoices[5]="Gamma";
2069 
2070  runAutoLSQName="Automatic Least Squares Refinement";
2071  runAutoLSQChoices[0]="Never";
2072  runAutoLSQChoices[1]="At the end of each run";
2073  runAutoLSQChoices[2]="Every 150000 trials, and at the end of each run";
2074 
2075  saveTrackedDataName="Save Tracked Data";
2076  saveTrackedDataChoices[0]="No (recommended!)";
2077  saveTrackedDataChoices[1]="Yes (for tests ONLY)";
2078 
2079  needInitNames=false;//Only once for the class
2080  }
2081  mGlobalOptimType.Init(2,&GlobalOptimTypeName,GlobalOptimTypeChoices);
2082  mAnnealingScheduleTemp.Init(6,&AnnealingScheduleTempName,AnnealingScheduleChoices);
2083  mAnnealingScheduleMutation.Init(6,&AnnealingScheduleMutationName,AnnealingScheduleChoices);
2084  mSaveTrackedData.Init(2,&saveTrackedDataName,saveTrackedDataChoices);
2085  mAutoLSQ.Init(3,&runAutoLSQName,runAutoLSQChoices);
2086  VFN_DEBUG_MESSAGE("MonteCarloObj::InitOptions():End",5)
2087 }
2088 
2089 void MonteCarloObj::InitLSQ(const bool useFullPowderPatternProfile)
2090 {
2091  mLSQ.SetRefinedObj(mRecursiveRefinedObjList.GetObj(0),0,true,true);
2092  for(unsigned int i=1;i<mRefinedObjList.GetNb();++i)
2093  mLSQ.SetRefinedObj(mRefinedObjList.GetObj(i),0,false,true);
2094 
2095  if(!useFullPowderPatternProfile)
2096  {// Use LSQ function #1 for powder patterns (integrated patterns - faster !)
2097  for(map<RefinableObj*,unsigned int>::iterator pos=mLSQ.GetRefinedObjMap().begin();pos!=mLSQ.GetRefinedObjMap().end();++pos)
2098  if(pos->first->GetClassName()=="PowderPattern") pos->second=1;
2099  }
2100  // Only refine structural parameters (excepting parameters already fixed) and scale factor
2101  mLSQ.PrepareRefParList(true);
2102  mLSQ.SetParIsFixed(gpRefParTypeScattData,true);
2103  mLSQ.SetParIsFixed(gpRefParTypeScattDataScale,false);
2104  mLSQ.SetParIsFixed(gpRefParTypeUnitCell,true);
2105  mLSQ.SetParIsFixed(gpRefParTypeScattPow,true);
2106  mLSQ.SetParIsFixed(gpRefParTypeRadiation,true);
2107 }
2108 
2109 #ifdef __WX__CRYST__
2110 WXCrystObjBasic* MonteCarloObj::WXCreate(wxWindow *parent)
2111 {
2112  mpWXCrystObj=new WXMonteCarloObj (parent,this);
2113  return mpWXCrystObj;
2114 }
2115 WXOptimizationObj* MonteCarloObj::WXGet()
2116 {
2117  return mpWXCrystObj;
2118 }
2119 void MonteCarloObj::WXDelete()
2120 {
2121  if(0!=mpWXCrystObj) delete mpWXCrystObj;
2122  mpWXCrystObj=0;
2123 }
2124 void MonteCarloObj::WXNotifyDelete()
2125 {
2126  mpWXCrystObj=0;
2127 }
2128 #endif
2129 
2130 }//namespace
ObjCryst::OptimizationObj::SetName
void SetName(const string &)
Set the name for this object.
Definition: GlobalOptimObj.cpp:248
ObjCryst::OptimizationObj::SetLimitsRelative
void SetLimitsRelative(const string &parName, const REAL min, const REAL max)
Change the relative limits for a parameter from its name.
Definition: GlobalOptimObj.cpp:162
ObjCryst::OptimizationObj::OptimizationObj
OptimizationObj(const string name="")
Constructor.
Definition: GlobalOptimObj.cpp:69
ObjCryst::OptimizationObj::XMLOutput
virtual void XMLOutput(ostream &os, int indent=0) const =0
Output a description of the object in XML format to a stream.
ObjCryst::MonteCarloObj::mAutoLSQ
RefObjOpt mAutoLSQ
Option to run automatic least-squares refinements.
Definition: GlobalOptimObj.h:510
ObjCryst::MonteCarloObj::RunSimulatedAnnealing
void RunSimulatedAnnealing(long &nbSteps, const bool silent=false, const REAL finalcost=0, const REAL maxTime=-1)
Definition: GlobalOptimObj.cpp:718
ObjCryst::RefinableObj::AddPar
void AddPar(const RefinablePar &newRefPar)
Add a refinable parameter.
Definition: RefinableObj.cpp:1416
ObjCryst::OptimizationObj::RandomizeStartingConfig
virtual void RandomizeStartingConfig()
Randomize starting configuration.
Definition: GlobalOptimObj.cpp:102
ObjCryst::RefinablePar::GetMin
REAL GetMin() const
Minimum value allowed (if limited or periodic)
Definition: RefinableObj.cpp:565
ObjCryst::MonteCarloObj::InitLSQ
virtual void InitLSQ(const bool useFullPowderPatternProfile=true)
Prepare mLSQ for least-squares refinement during the global optimization.
Definition: GlobalOptimObj.cpp:2089
ObjCryst::TrackerObject
Tracker for objects (RefinableObj, Crystal, PowderPattern, RefPar,...)
Definition: Tracker.h:109
ObjCryst::OptimizationObj::mName
string mName
Name of the GlobalOptimization object.
Definition: GlobalOptimObj.h:244
ObjCryst::MonteCarloObj::Optimize
virtual void Optimize(long &nbSteps, const bool silent=false, const REAL finalcost=0, const REAL maxTime=-1)
Launch optimization (a single run) for N steps.
Definition: GlobalOptimObj.cpp:499
ObjCryst::MonteCarloObj::mNbTrialRetry
long mNbTrialRetry
Number of trials before testing if we are below the given minimum cost.
Definition: GlobalOptimObj.h:504
ObjCryst::LSQNumObj::PrepareRefParList
void PrepareRefParList(const bool copy_param=false)
Prepare the full parameter list for the refinement.
Definition: LSQNumObj.cpp:816
ObjCryst::MonteCarloObj::mTemperature
REAL mTemperature
Current temperature for annealing.
Definition: GlobalOptimObj.h:478
ObjCryst::RefinableObjClock
We need to record exactly when refinable objects have been modified for the last time (to avoid re-co...
Definition: RefinableObj.h:138
ObjCryst::OptimizationObj::BuildRecursiveRefObjList
void BuildRecursiveRefObjList()
(Re)build OptimizationObj::mRecursiveRefinedObjList, if an object has been added or modified...
Definition: GlobalOptimObj.cpp:391
ObjCryst::MonteCarloObj::mCurrentCost
REAL mCurrentCost
Current value of the cost function.
Definition: GlobalOptimObj.h:469
ObjCryst::OptimizationObj::AddRefinableObj
void AddRefinableObj(RefinableObj &)
Add a refined object. All sub-objects are also added.
Definition: GlobalOptimObj.cpp:228
ObjCryst::AnnealingSchedule
AnnealingSchedule
Annealing schedule type.
Definition: GlobalOptimObj.h:66
ObjCryst::LSQNumObj::SetRefinedObj
void SetRefinedObj(RefinableObj &obj, const unsigned int LSQFuncIndex=0, const bool init=true, const bool recursive=false)
Choose the object to refine.
Definition: LSQNumObj.cpp:727
ObjCryst::RefinableObj::SaveParamSet
void SaveParamSet(const unsigned long id) const
Save the current set of refined values over a previously-created set of saved values.
Definition: RefinableObj.cpp:1527
ObjCryst::MonteCarloObj::mTemperatureMax
REAL mTemperatureMax
Beginning temperature for annealing.
Definition: GlobalOptimObj.h:480
ObjCryst::MonteCarloObj::mGlobalOptimType
RefObjOpt mGlobalOptimType
Method used for the global optimization.
Definition: GlobalOptimObj.h:465
Chronometer
Simple chronometer class, with microsecond precision.
Definition: Chronometer.h:33
ObjCryst::Crystal::GetBumpMergeCost
REAL GetBumpMergeCost() const
Get the Anti-bumping/pro-Merging cost function.
Definition: Crystal.cpp:820
ObjCryst::MonteCarloObj::GetLSQObj
LSQNumObj & GetLSQObj()
Access to the builtin LSQ optimization object.
Definition: GlobalOptimObj.cpp:2020
ObjCryst::MainTracker::UpdateDisplay
void UpdateDisplay() const
Update display, if any.
Definition: Tracker.cpp:133
ObjCryst::RefObjRegisterRecursive
void RefObjRegisterRecursive(T &obj, ObjRegistry< T > &reg)
Register a new object in a registry, and recursively include all included (sub)objects.
Definition: RefinableObj.cpp:1147
ObjCryst::gOptimizationObjRegistry
ObjRegistry< OptimizationObj > gOptimizationObjRegistry("List of all Optimization objects")
Global Registry for all OptimizationObj.
Definition: GlobalOptimObj.h:525
ObjCryst::RefinableObj::GetNbPar
long GetNbPar() const
Total number of refinable parameter in the object.
Definition: RefinableObj.cpp:1348
ObjCryst::OptimizationObj::SetParIsFixed
void SetParIsFixed(const string &parName, const bool fix)
Fix one parameter.
Definition: GlobalOptimObj.cpp:130
ObjCryst::RefinableObj::GetParamSet
const CrystVector_REAL & GetParamSet(const unsigned long setId) const
Access one save refpar set.
Definition: RefinableObj.cpp:1550
ObjCryst::RefinableObj::GetParNotFixed
RefinablePar & GetParNotFixed(const long i)
Access all parameters in the order they were inputted, skipping fixed parameters. ...
Definition: RefinableObj.cpp:1406
ObjCryst::RefinableObj::GetNbParNotFixed
long GetNbParNotFixed() const
Total number of non-fixed parameters. Is initialized by PrepareForRefinement()
Definition: RefinableObj.cpp:1350
ObjCryst::RefinableObj::GetPar
RefinablePar & GetPar(const long i)
Access all parameters in the order they were inputted.
Definition: RefinableObj.cpp:1352
ObjCryst::OptimizationObj::SetParIsUsed
void SetParIsUsed(const string &parName, const bool use)
Set a parameter to be used.
Definition: GlobalOptimObj.cpp:150
ObjCryst::OptimizationObj::DisplayReport
virtual void DisplayReport()
Show report to the user during refinement. Used for GUI update.
Definition: GlobalOptimObj.cpp:223
ObjCryst::RefObjOpt::XMLInput
void XMLInput(istream &is, const XMLCrystTag &tag)
XMLInput From stream.
Definition: RefinableObj/IO.cpp:313
ObjCryst::MainTracker::ClearValues
void ClearValues()
Removes all stored values.
Definition: Tracker.cpp:97
ObjCryst::OptimizationObj::GetLastOptimElapsedTime
REAL GetLastOptimElapsedTime() const
Get the elapsed time (in seconds) during the last optimization.
Definition: GlobalOptimObj.cpp:270
ObjCryst::OptimizationObj::FixAllPar
void FixAllPar()
Fix all parameters.
Definition: GlobalOptimObj.cpp:122
ObjCryst::RefinableObj
Generic Refinable Object.
Definition: RefinableObj.h:752
ObjCryst::LSQNumObj::GetRefinedObjMap
const std::map< RefinableObj *, unsigned int > & GetRefinedObjMap() const
Get the map of refined objects - this is a recursive list of all the objects that are taken into acco...
Definition: LSQNumObj.cpp:740
ObjCryst::OptimizationObj::LogLikelihoodStats::mTotalLogLikelihood
REAL mTotalLogLikelihood
Total Log(Likelihood), to compute the average.
Definition: GlobalOptimObj.h:266
ObjCryst::gpRefParTypeScattDataScale
const RefParType * gpRefParTypeScattDataScale
Type for scattering data scale factors.
Definition: ScatteringData.cpp:56
ObjCryst::MonteCarloObj::MonteCarloObj
MonteCarloObj(const string name="")
Constructor.
Definition: GlobalOptimObj.cpp:412
ObjCryst::RefinableObj::ResetParList
void ResetParList()
Re-init the list of refinable parameters, removing all parameters.
Definition: RefinableObj.cpp:1829
ObjCryst::OptimizationObj::mStopAfterCycle
bool mStopAfterCycle
If true, then stop at the end of the cycle. Used in multi-threaded environment.
Definition: GlobalOptimObj.h:294
ObjCryst::OptimizationObj::mMainTracker
MainTracker mMainTracker
MainTracker object to track the evolution of cost functions, likelihood, and individual parameters...
Definition: GlobalOptimObj.h:310
ObjCryst::gRefinableObjRegistry
ObjRegistry< RefinableObj > gRefinableObjRegistry("Global RefinableObj registry")
Global Registry for all RefinableObj.
Definition: RefinableObj.h:1312
ObjCryst::LSQNumObj::SetParIsFixed
void SetParIsFixed(const std::string &parName, const bool fix)
Fix one parameter.
ObjCryst::MonteCarloObj::mTemperatureMin
REAL mTemperatureMin
Lower temperature.
Definition: GlobalOptimObj.h:482
ObjCryst::WXCrystObjBasic
Abstract base class for all objects in wxCryst.
Definition: wxCryst.h:127
ObjCryst::MonteCarloObj::mMutationAmplitudeGamma
REAL mMutationAmplitudeGamma
Gamma for the 'gamma' Mutation amplitude schedule.
Definition: GlobalOptimObj.h:500
ObjCryst::MonteCarloObj::mMutationAmplitude
REAL mMutationAmplitude
Mutation amplitude.
Definition: GlobalOptimObj.h:492
ObjCryst::GetRefParListClockRecursive
void GetRefParListClockRecursive(ObjRegistry< RefinableObj > &reg, RefinableObjClock &clock)
Get the last time any RefinablePar was added in a recursive list of objects.
Definition: RefinableObj.cpp:2033
ObjCryst::OptimizationObj::GetName
const string & GetName() const
Get the name for this object.
Definition: GlobalOptimObj.cpp:247
ObjCryst::RefinableObj::PrepareForRefinement
void PrepareForRefinement() const
Find which parameters are used and not fixed, for a refinement /optimization.
Definition: RefinableObj.cpp:1273
ObjCryst::OptimizationObj::mBestCost
REAL mBestCost
Best value of the cost function so far.
Definition: GlobalOptimObj.h:254
ObjCryst::Crystal::GetBondValenceCost
REAL GetBondValenceCost() const
Get the Bond-Valence cost function, which compares the expected valence to the one computed from Bond...
Definition: Crystal.cpp:1205
ObjCryst::MonteCarloObj::mLSQ
LSQNumObj mLSQ
Least squares object.
Definition: GlobalOptimObj.h:508
ObjCryst::OptimizationObj::~OptimizationObj
virtual ~OptimizationObj()
Destructor.
Definition: GlobalOptimObj.cpp:95
ObjCryst::RefinableObj::EraseAllParamSet
const void EraseAllParamSet()
Erase all saved refpar sets.
Definition: RefinableObj.cpp:1572
ObjCryst::OptimizationObj::mBestParSavedSetIndex
long mBestParSavedSetIndex
Index of the 'best' saved parameter set.
Definition: GlobalOptimObj.h:256
ObjCryst::RefObjOpt
Base class for options.
Definition: RefinableObj.h:550
ObjCryst::OptimizationObj::Print
virtual void Print() const
Print some information about this object.
Definition: GlobalOptimObj.cpp:252
ObjCryst::XMLCrystFileSaveGlobal
void XMLCrystFileSaveGlobal(const string &filename)
Save all Objcryst++ objects.
Definition: ObjCryst/IO.cpp:106
ObjCryst::Molecule
Molecule : class for complex scatterer descriptions using cartesian coordinates with bond length/angl...
Definition: Molecule.h:731
ObjCryst::MonteCarloObj::NewConfiguration
virtual void NewConfiguration(const RefParType *type=gpRefParTypeObjCryst)
Make a random change in the configuration.
Definition: GlobalOptimObj.cpp:2024
ObjCryst::OptimizationObj::TagNewBestConfig
void TagNewBestConfig()
During a global optimization, tell all objects that the current config is the latest "best" config...
Definition: GlobalOptimObj.cpp:263
ObjCryst::MonteCarloObj::mMutationAmplitudeMin
REAL mMutationAmplitudeMin
Mutation amplitude at the end of the optimization.
Definition: GlobalOptimObj.h:496
ObjCryst::RefinableObj::RestoreParamSet
void RestoreParamSet(const unsigned long id)
Restore a saved set of values.
Definition: RefinableObj.cpp:1536
ObjCryst::MonteCarloObj::SetAlgorithmParallTempering
void SetAlgorithmParallTempering(const AnnealingSchedule scheduleTemp, const REAL tMax, const REAL tMin, const AnnealingSchedule scheduleMutation=ANNEALING_CONSTANT, const REAL mutMax=16., const REAL mutMin=.125)
Set the refinement method to Parallel Tempering.
Definition: GlobalOptimObj.cpp:481
ObjCryst::OptimizationObj::SetLimitsAbsolute
void SetLimitsAbsolute(const string &parName, const REAL min, const REAL max)
Change the absolute limits for a parameter from its name.
Definition: GlobalOptimObj.cpp:176
ObjCryst::MonteCarloObj::mAnnealingScheduleTemp
RefObjOpt mAnnealingScheduleTemp
Schedule for the annealing.
Definition: GlobalOptimObj.h:484
ObjCryst::RefinablePar::GetMax
REAL GetMax() const
Get the maximum value allowed (if limited)
Definition: RefinableObj.cpp:570
ObjCryst::RefObjOpt::XMLOutput
void XMLOutput(ostream &os, int indent=0) const
XMLOutput to stream in well-formed XML.
Definition: RefinableObj/IO.cpp:295
ObjCryst::RefinablePar::GetPeriod
REAL GetPeriod() const
Get the period (if periodic)
Definition: RefinableObj.cpp:575
ObjCryst::OptimizationObj::mContext
unsigned long mContext
The current 'context', in the case the optimization is run in different parallel contexts.
Definition: GlobalOptimObj.h:259
ObjCryst::OptimizationObj::GetFullRefinableObj
RefinableObj & GetFullRefinableObj(const bool rebuild=true)
Get the RefinableObj with all the parameters from all refined objects.
Definition: GlobalOptimObj.cpp:241
ObjCryst::OptimizationObj::IsOptimizing
bool IsOptimizing() const
Are we busy optimizing ?
Definition: GlobalOptimObj.cpp:261
ObjCryst::MonteCarloObj::MultiRunOptimize
virtual void MultiRunOptimize(long &nbCycle, long &nbSteps, const bool silent=false, const REAL finalcost=0, const REAL maxTime=-1)
Launch optimization for multiple runs of N steps.
Definition: GlobalOptimObj.cpp:576
ObjCryst::RefinableObj::GetRefParListClock
const RefinableObjClock & GetRefParListClock() const
What was the last time a RefinablePar was added/removed ?
Definition: RefinableObj.cpp:1878
ObjCryst::Crystal::GetScatt
Scatterer & GetScatt(const string &scattName)
Provides an access to the scatterers.
Definition: Crystal.cpp:198
ObjCryst::OptimizationObj::LogLikelihoodStats
Statistics about each object contributing to the overall Log(likelihood)
Definition: GlobalOptimObj.h:261
ObjCryst::OptimizationObj::mNbTrialPerRun
long mNbTrialPerRun
Number of trial per run, to be saved/restored in XML output.
Definition: GlobalOptimObj.h:250
ObjCryst::OptimizationObj::mRefParList
RefinableObj mRefParList
The refinable par list used during refinement.
Definition: GlobalOptimObj.h:242
ObjCryst::OptimizationObj::GetClassName
virtual const string GetClassName() const
Get the name for this class type.
Definition: GlobalOptimObj.cpp:250
ObjCryst::OptimizationObj::mXMLAutoSave
RefObjOpt mXMLAutoSave
Periodic save of complete environment as an xml file.
Definition: GlobalOptimObj.h:304
ObjCryst::MonteCarloObj::RunParallelTempering
void RunParallelTempering(long &nbSteps, const bool silent=false, const REAL finalcost=0, const REAL maxTime=-1)
Definition: GlobalOptimObj.cpp:1147
ObjCryst::Crystal::GetScattererRegistry
ObjRegistry< Scatterer > & GetScattererRegistry()
Get the registry of scatterers.
Definition: Crystal.cpp:218
ObjCryst::MainTracker
A class to hold all trackers.
Definition: Tracker.h:68
ObjCryst::OptimizationObj::RestoreBestConfiguration
void RestoreBestConfiguration()
Restore the Best configuration.
Definition: GlobalOptimObj.cpp:254
ObjCryst::MonteCarloObj::~MonteCarloObj
virtual ~MonteCarloObj()
Destructor.
Definition: GlobalOptimObj.cpp:454
ObjCryst::OptimizationObj::mRefinedObjList
ObjRegistry< RefinableObj > mRefinedObjList
The refined objects.
Definition: GlobalOptimObj.h:298
ObjCryst::OptimizationObj::mLastOptimTime
REAL mLastOptimTime
The time elapsed after the last optimization, in seconds.
Definition: GlobalOptimObj.h:307
ObjCryst::OptimizationObj::UnFixAllPar
void UnFixAllPar()
UnFix All parameters.
Definition: GlobalOptimObj.cpp:143
ObjCryst::OptimizationObj::GetMainTracker
MainTracker & GetMainTracker()
Get the MainTracker.
Definition: GlobalOptimObj.cpp:275
ObjCryst::RefinablePar::MutateTo
void MutateTo(const REAL newValue)
Change the current value to the given one.
Definition: RefinableObj.cpp:487
ObjCryst::RefinableObj::ClearParamSet
void ClearParamSet(const unsigned long id) const
Erase the param set with the given id, releasing memory.
Definition: RefinableObj.cpp:1520
ObjCryst::OptimizationObj::mvObjWeight
map< const RefinableObj *, DynamicObjWeight > mvObjWeight
Weights for each objects in each context (mutable for dynamic update during optimization) ...
Definition: GlobalOptimObj.h:281
ObjCryst::WXMonteCarloObj
Class for Graphical interface to Monte-Carlo objects (Simulated Annealing, Parallel Tempering) ...
Definition: wxGlobalOptimObj.h:96
ObjCryst::MonteCarloObj::mMutationAmplitudeMax
REAL mMutationAmplitudeMax
Mutation amplitude at the beginning of the optimization.
Definition: GlobalOptimObj.h:494
ObjCryst::OptimizationObj
Base object for Optimization methods.
Definition: GlobalOptimObj.h:102
ObjCryst::MainTracker::ClearTrackers
void ClearTrackers()
Removes all Trackers.
Definition: Tracker.cpp:88
ObjCryst::OptimizationObj::mRecursiveRefinedObjList
ObjRegistry< RefinableObj > mRecursiveRefinedObjList
The refined objects, recursively including all sub-objects.
Definition: GlobalOptimObj.h:301
ObjCryst::ObjCrystException
Exception class for ObjCryst++ library.
Definition: General.h:119
ObjCryst::MonteCarloObj::mMinCostRetry
REAL mMinCostRetry
Cost to reach unless an automatic randomization and retry is done.
Definition: GlobalOptimObj.h:506
ObjCryst::OptimizationObj::StopAfterCycle
void StopAfterCycle()
Stop after the current cycle. USed for interactive refinement.
Definition: GlobalOptimObj.cpp:211
ObjCryst
The namespace which includes all objects (crystallographic and algorithmic) in ObjCryst++.
Definition: Atom.cpp:47
ObjCryst::RefinableObj::GetLogLikelihood
virtual REAL GetLogLikelihood() const
Get -log(likelihood) of the current configuration for the object.
Definition: RefinableObj.cpp:1759
ObjCryst::MonteCarloObj::mAnnealingScheduleMutation
RefObjOpt mAnnealingScheduleMutation
Schedule for the annealing.
Definition: GlobalOptimObj.h:498
ObjCryst::OptimizationObj::EndOptimization
virtual void EndOptimization()
End optimization for all objects.
Definition: GlobalOptimObj.cpp:293
ObjCryst::MonteCarloObj::GetClassName
virtual const string GetClassName() const
Get the name for this class type.
Definition: GlobalOptimObj.cpp:2018
ObjCryst::MonteCarloObj::XMLInput
virtual void XMLInput(istream &is, const XMLCrystTag &tag)
Input in XML format from a stream, restoring the set of refined objects and the associated cost funct...
Definition: GlobalOptimObj.cpp:1930
ObjCryst::OptimizationObj::mvContextObjStats
map< unsigned long, map< const RefinableObj *, LogLikelihoodStats > > mvContextObjStats
Statistics for each context (mutable for dynamic update during optimization)
Definition: GlobalOptimObj.h:272
ObjCryst::RefinableObj::GetName
virtual const string & GetName() const
Name of the object.
Definition: RefinableObj.cpp:1249
ObjCryst::OptimizationObj::mIsOptimizing
bool mIsOptimizing
True if a refinement is being done. For multi-threaded environment.
Definition: GlobalOptimObj.h:292
ObjCryst::MainTracker::SaveAll
void SaveAll(std::ostream &out) const
Will save to a single file if all recorded trial numbers are the same Otherwise ? ...
Definition: Tracker.cpp:105
ObjCryst::GetSubRefObjListClockRecursive
void GetSubRefObjListClockRecursive(ObjRegistry< RefinableObj > &reg, RefinableObjClock &clock)
Get the last time any object was added in the recursive list of objects.
Definition: RefinableObj.cpp:1160
ObjCryst::OptimizationObj::InitOptions
virtual void InitOptions()
Initialization of options.
Definition: GlobalOptimObj.cpp:358
ObjCryst::Crystal
Crystal class: Unit cell, spacegroup, scatterers.
Definition: Crystal.h:97
ObjCryst::MonteCarloObj::XMLOutput
virtual void XMLOutput(ostream &os, int indent=0) const
Output a description of the object in XML format to a stream.
Definition: GlobalOptimObj.cpp:1869
ObjCryst::LSQNumObj
(Quick & dirty) Least-Squares Refinement Object with Numerical derivatives
Definition: LSQNumObj.h:37
ObjCryst::OptimizationObj::BeginOptimization
virtual void BeginOptimization(const bool allowApproximations=false, const bool enableRestraints=false)
Begin optimization for all objects.
Definition: GlobalOptimObj.cpp:285
ObjCryst::OptimizationObj::LogLikelihoodStats::mLastLogLikelihood
REAL mLastLogLikelihood
Previous log(likelihood)
Definition: GlobalOptimObj.h:264
ObjCryst::OptimizationObj::UpdateDisplay
virtual void UpdateDisplay()
Update Display (if any display is available), when a new 'relevant' configuration is reached...
Definition: GlobalOptimObj.cpp:381
ObjCryst::XMLCrystTag
class to input or output a well-formatted xml beginning or ending tag.
Definition: RefinableObj/IO.h:72
ObjCryst::MonteCarloObj::mTemperatureGamma
REAL mTemperatureGamma
Gamma for the 'gamma' temperature schedule.
Definition: GlobalOptimObj.h:486
ObjCryst::MonteCarloObj::mSaveTrackedData
RefObjOpt mSaveTrackedData
Option to save the evolution of tracked data (cost functions, likelihhod, individual parameters...
Definition: GlobalOptimObj.h:474
ObjCryst::RefinableObj::SetDeleteRefParInDestructor
void SetDeleteRefParInDestructor(const bool b)
Set this object not to delete its list of parameters when destroyed.
Definition: RefinableObj.cpp:1876
ObjCryst::RefParType
class of refinable parameter types.
Definition: RefinableObj.h:78
ObjCryst::MonteCarloObj::SetAlgorithmSimulAnnealing
void SetAlgorithmSimulAnnealing(const AnnealingSchedule scheduleTemp, const REAL tMax, const REAL tMin, const AnnealingSchedule scheduleMutation=ANNEALING_CONSTANT, const REAL mutMax=16., const REAL mutMin=.125, const long nbTrialRetry=0, const REAL minCostRetry=0.)
Set the refinement method to simulated Annealing.
Definition: GlobalOptimObj.cpp:460
ObjCryst::OptimizationObj::mNbTrial
long mNbTrial
Number of trials so far.
Definition: GlobalOptimObj.h:252
ObjCryst::OptimizationObj::NbTrialPerRun
virtual long & NbTrialPerRun()
Number of trial per run.
Definition: GlobalOptimObj.cpp:298
ObjCryst::gpRefParTypeScattData
const RefParType * gpRefParTypeScattData
Generic type for scattering data.
Definition: ScatteringData.cpp:55
ObjCryst::RefinableObj::CreateParamSet
unsigned long CreateParamSet(const string name="") const
Save the current set of refined values in a new set.
Definition: RefinableObj.cpp:1503
ObjCryst::OptimizationObj::mvSavedParamSet
std::vector< pair< long, REAL > > mvSavedParamSet
List of saved parameter sets.
Definition: GlobalOptimObj.h:289
ObjCryst::gpRefParTypeObjCryst
const RefParType * gpRefParTypeObjCryst
Top RefParType for the ObjCryst++ library.
Definition: RefinableObj.cpp:77
ObjCryst::OptimizationObj::GetLogLikelihood
virtual REAL GetLogLikelihood() const
The optimized (minimized, actually) function.
Definition: GlobalOptimObj.cpp:191
ObjCryst::OptimizationObj::GetBestCost
const REAL & GetBestCost() const
Access to current best cost.
Definition: GlobalOptimObj.cpp:282
ObjCryst::LSQNumObj::Refine
void Refine(int nbCycle=1, bool useLevenbergMarquardt=false, const bool silent=false, const bool callBeginEndOptimization=true, const float minChi2var=0.01)
Do the refinement.
Definition: LSQNumObj.cpp:104
ObjCryst::MonteCarloObj::InitOptions
virtual void InitOptions()
Initialization of options.
Definition: GlobalOptimObj.cpp:2035
ObjCryst::OptimizationObj::LogLikelihoodStats::mTotalLogLikelihoodDeltaSq
REAL mTotalLogLikelihoodDeltaSq
total of (Delta(Log(Likelihood)))^2 between successive trials
Definition: GlobalOptimObj.h:268