RefinePowderInstrumentParameters3.cpp 45.6 KB
Newer Older
1
#include "MantidCurveFitting/Algorithms/RefinePowderInstrumentParameters3.h"
2
3
4
5

#include "MantidAPI/Axis.h"
#include "MantidAPI/TextAxis.h"
#include "MantidAPI/TableRow.h"
6
#include "MantidAPI/WorkspaceFactory.h"
7
8
#include "MantidKernel/ListValidator.h"

9
10
using namespace Mantid::API;
using namespace Mantid::CurveFitting;
11
using namespace Mantid::CurveFitting::Functions;
12
13
14
15
using namespace Mantid::DataObjects;
using namespace Mantid::Kernel;
using namespace std;

16
17
namespace Mantid {
namespace CurveFitting {
18
namespace Algorithms {
19

20
DECLARE_ALGORITHM(RefinePowderInstrumentParameters3)
21

22
23
24
//----------------------------------------------------------------------------------------------
/** Constructor
 */
25
RefinePowderInstrumentParameters3::RefinePowderInstrumentParameters3()
26
27
28
29
    : m_dataWS(), m_wsIndex(-1), m_paramTable(), m_fitMode(MONTECARLO),
      m_stdMode(CONSTANT), m_numWalkSteps(-1), m_randomSeed(-1),
      m_profileParameters(), m_positionFunc(), m_dampingFactor(0.),
      m_bestChiSq(0.), m_bestChiSqStep(-1), m_bestChiSqGroup(-1) {}
30

31
32
33
34
//----------------------------------------------------------------------------------------------
/** Destructor
 */
RefinePowderInstrumentParameters3::~RefinePowderInstrumentParameters3() {}
35

36
37
38
39
40
41
//----------------------------------------------------------------------------------------------
/** Declare properties
  */
void RefinePowderInstrumentParameters3::init() {
  // Peak position workspace
  declareProperty(
42
43
      Kernel::make_unique<WorkspaceProperty<Workspace2D>>(
          "InputPeakPositionWorkspace", "Anonymous", Direction::Input),
44
45
46
47
48
49
50
51
52
      "Data workspace containing workspace positions in TOF agains dSpacing.");

  // Workspace Index
  declareProperty(
      "WorkspaceIndex", 0,
      "Workspace Index of the peak positions in PeakPositionWorkspace.");

  // Output workspace
  declareProperty(
53
54
      Kernel::make_unique<WorkspaceProperty<Workspace2D>>(
          "OutputPeakPositionWorkspace", "Anonymous2", Direction::Output),
55
56
57
58
59
      "Output data workspace containing refined workspace positions in TOF "
      "agains dSpacing.");

  // Input Table workspace containing instrument profile parameters
  declareProperty(
60
61
      Kernel::make_unique<WorkspaceProperty<TableWorkspace>>(
          "InputInstrumentParameterWorkspace", "Anonymous3", Direction::Input),
62
63
64
65
      "INput tableWorkspace containg instrument's parameters.");

  // Output table workspace containing the refined parameters
  declareProperty(
66
      Kernel::make_unique<WorkspaceProperty<TableWorkspace>>(
67
68
69
70
71
          "OutputInstrumentParameterWorkspace", "Anonymous4",
          Direction::Output),
      "Output tableworkspace containing instrument's fitted parameters. ");

  // Refinement algorithm
72
  vector<string> algoptions{"OneStepFit", "MonteCarlo"};
73
74
75
76
77
78
79
80
81
82
83
84
85
  auto validator = boost::make_shared<Kernel::StringListValidator>(algoptions);
  declareProperty("RefinementAlgorithm", "MonteCarlo", validator,
                  "Algorithm to refine the instrument parameters.");

  // Random walk steps
  declareProperty("RandomWalkSteps", 10000,
                  "Number of Monte Carlo random walk steps. ");

  // Random seed
  declareProperty("MonteCarloRandomSeed", 0,
                  "Random seed for Monte Carlo simulation. ");

  // Method to calcualte the standard error of peaks
86
  vector<string> stdoptions{"ConstantValue", "UseInputValue"};
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
  auto listvalidator =
      boost::make_shared<Kernel::StringListValidator>(stdoptions);
  declareProperty(
      "StandardError", "ConstantValue", listvalidator,
      "Algorithm to calculate the standard error of peak positions.");

  // Damping factor
  declareProperty(
      "Damping", 1.0,
      "Damping factor for (1) minimizer 'damping'. (2) Monte Calro. ");

  // Anealing temperature
  declareProperty("AnnealingTemperature", 1.0,
                  "Starting aneealing temperature.");

  // Monte Carlo iterations
  declareProperty("MonteCarloIterations", 100,
                  "Number of iterations in Monte Carlo random walk.");

  // Output
  declareProperty("ChiSquare", DBL_MAX, Direction::Output);

  return;
}

//----------------------------------------------------------------------------------------------
/** Main execution body
  */
void RefinePowderInstrumentParameters3::exec() {
  // 1. Process input
  processInputProperties();

  // 2. Parse input table workspace
  parseTableWorkspaces();

  // 3. Set up main function for peak positions
  ThermalNeutronDtoTOFFunction rawfunc;
  m_positionFunc = boost::make_shared<ThermalNeutronDtoTOFFunction>(rawfunc);
  m_positionFunc->initialize();

  // 3. Fit
  // a) Set up parameter value
  setFunctionParameterValues(m_positionFunc, m_profileParameters);

  // b) Generate some global useful value and Calculate starting chi^2
  API::FunctionDomain1DVector domain(m_dataWS->readX(m_wsIndex));
  API::FunctionValues rawvalues(domain);
  m_positionFunc->function(domain, rawvalues);

  // d) Calcualte statistic
  double startchi2 =
      calculateFunctionError(m_positionFunc, m_dataWS, m_wsIndex);

  // b) Fit by type
  double finalchi2 = DBL_MAX;
  switch (m_fitMode) {
  case FIT:
    // Fit by non-Monte Carlo method
    g_log.notice("Fit by non Monte Carlo algorithm. ");
    finalchi2 = execFitParametersNonMC();
    break;

  case MONTECARLO:
    // Fit by Monte Carlo method
    g_log.notice("Fit by Monte Carlo algorithm.");
    finalchi2 = execFitParametersMC();
    break;

  default:
    // Unsupported
    throw runtime_error("Unsupported fit mode.");
    break;
159
160
  }

161
162
163
164
  // 4. Process the output
  TableWorkspace_sptr fitparamtable =
      genOutputProfileTable(m_profileParameters, startchi2, finalchi2);
  setProperty("OutputInstrumentParameterWorkspace", fitparamtable);
165

166
167
  Workspace2D_sptr outdataws = genOutputWorkspace(domain, rawvalues);
  setProperty("OutputPeakPositionWorkspace", outdataws);
168

169
  setProperty("ChiSquare", finalchi2);
170

171
172
  return;
}
173

174
175
176
177
178
179
180
181
182
183
184
//----------------------------------------------------------------------------------------------
/** Process input properties
  */
void RefinePowderInstrumentParameters3::processInputProperties() {
  // Data Workspace
  m_dataWS = getProperty("InputPeakPositionWorkspace");

  m_wsIndex = getProperty("WorkspaceIndex");
  if (m_wsIndex < 0 ||
      m_wsIndex >= static_cast<int>(m_dataWS->getNumberHistograms())) {
    throw runtime_error("Input workspace index is out of range.");
185
186
  }

187
188
189
190
191
192
193
194
195
196
197
198
199
  // Parameter TableWorkspace
  m_paramTable = getProperty("InputInstrumentParameterWorkspace");

  // Fit mode
  string fitmode = getProperty("RefinementAlgorithm");
  if (fitmode.compare("OneStepFit") == 0)
    m_fitMode = FIT;
  else if (fitmode.compare("MonteCarlo") == 0)
    m_fitMode = MONTECARLO;
  else {
    m_fitMode = FIT;
    throw runtime_error("Input RefinementAlgorithm is not supported.");
  }
200

201
202
203
204
205
206
207
208
209
  // Stanard error mode
  string stdmode = getProperty("StandardError");
  if (stdmode.compare("ConstantValue") == 0)
    m_stdMode = CONSTANT;
  else if (stdmode.compare("UseInputValue") == 0)
    m_stdMode = USEINPUT;
  else {
    m_stdMode = USEINPUT;
    throw runtime_error("Input StandardError (mode) is not supported.");
210
211
  }

212
213
214
215
216
  // Monte Carlo
  m_numWalkSteps = getProperty("RandomWalkSteps");
  if (m_numWalkSteps <= 0)
    throw runtime_error(
        "Monte Carlo walk steps cannot be less or equal to 0. ");
217

218
  m_randomSeed = getProperty("MonteCarloRandomSeed");
219

220
  m_dampingFactor = getProperty("Damping");
221

222
223
  return;
}
224

225
226
227
228
229
//----------------------------------------------------------------------------------------------
/** Parse TableWorkspaces
  */
void RefinePowderInstrumentParameters3::parseTableWorkspaces() {
  m_profileParameters.clear();
230

231
232
  parseTableWorkspace(m_paramTable, m_profileParameters);
}
233

234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
//----------------------------------------------------------------------------------------------
/** Parse table workspace to a map of Parameters
  */
void RefinePowderInstrumentParameters3::parseTableWorkspace(
    TableWorkspace_sptr tablews, map<string, Parameter> &parammap) {
  // 1. Process Table column names
  std::vector<std::string> colnames = tablews->getColumnNames();
  map<string, size_t> colnamedict;
  convertToDict(colnames, colnamedict);

  int iname = getStringIndex(colnamedict, "Name");
  int ivalue = getStringIndex(colnamedict, "Value");
  int ifit = getStringIndex(colnamedict, "FitOrTie");
  int imin = getStringIndex(colnamedict, "Min");
  int imax = getStringIndex(colnamedict, "Max");
  int istep = getStringIndex(colnamedict, "StepSize");

  if (iname < 0 || ivalue < 0 || ifit < 0)
    throw runtime_error(
        "TableWorkspace does not have column Name, Value and/or Fit.");

  // 3. Parse
  size_t numrows = tablews->rowCount();
  for (size_t irow = 0; irow < numrows; ++irow) {
    string parname = tablews->cell<string>(irow, iname);
    double parvalue = tablews->cell<double>(irow, ivalue);
    string fitq = tablews->cell<string>(irow, ifit);

    double minvalue;
    if (imin >= 0)
      minvalue = tablews->cell<double>(irow, imin);
    else
      minvalue = -DBL_MAX;
267

268
269
270
271
272
    double maxvalue;
    if (imax >= 0)
      maxvalue = tablews->cell<double>(irow, imax);
    else
      maxvalue = DBL_MAX;
273

274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
    double stepsize;
    if (istep >= 0)
      stepsize = tablews->cell<double>(irow, istep);
    else
      stepsize = 1.0;

    Parameter newpar;
    newpar.name = parname;
    newpar.curvalue = parvalue;
    newpar.minvalue = minvalue;
    newpar.maxvalue = maxvalue;
    newpar.stepsize = stepsize;

    // If empty string, fit is default to be false
    bool fit = false;
    if (fitq.size() > 0) {
      if (fitq[0] == 'F' || fitq[0] == 'f')
        fit = true;
    }
    newpar.fit = fit;
294

295
    parammap.emplace(parname, newpar);
296
297
  }

298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
  return;
}

//----------------------------------------------------------------------------------------------
/** Fit instrument parameters by non Monte Carlo algorithm
  * Requirement:  m_positionFunc should have the best fit result;
 */
double RefinePowderInstrumentParameters3::execFitParametersNonMC() {
  // 1. Set up constraints
  setFunctionParameterFitSetups(m_positionFunc, m_profileParameters);

  // 2. Fit function
  // FIXME powerfit should be a user option before freezing this algorithm
  // FIXME powdefit = True introduce segmentation fault
  bool powerfit = false;
  double chi2 = fitFunction(m_positionFunc, m_dataWS, m_wsIndex, powerfit);

  // 2. Summary
  stringstream sumss;
  sumss << "Non-Monte Carlo Results:  Best Chi^2 = " << chi2;
  g_log.notice(sumss.str());

  return chi2;
}

//----------------------------------------------------------------------------------------------
/** Refine instrument parameters by Monte Carlo/simulated annealing method
  */
double RefinePowderInstrumentParameters3::execFitParametersMC() {
  // 1. Monte Carlo simulation
  double chisq = doSimulatedAnnealing(m_profileParameters);

  // 2. Summary
  stringstream sumss;
  sumss << "Monte Carlo Results:  Best Chi^2 = " << chisq << " @ Step "
        << m_bestChiSqStep << ", Group " << m_bestChiSqGroup;
  g_log.notice(sumss.str());

  return chisq;
}

//----------------------------------------------------------------------------------------------
/** Do MC/simulated annealing to refine parameters
  *
  * Helpful:     double curchi2 = calculateD2TOFFunction(mFunction, domain,
  *values, rawY, rawE);
  */
double RefinePowderInstrumentParameters3::doSimulatedAnnealing(
    map<string, Parameter> inparammap) {
  // 1. Prepare/initialization
  //    Data structure
  const MantidVec &dataY = m_dataWS->readY(m_wsIndex);
350

351
  size_t numpts = dataY.size();
352

353
  vector<double> vecY(numpts, 0.0);
354

355
356
357
  //    Monte Carlo strategy and etc.
  vector<vector<string>> mcgroups;
  setupRandomWalkStrategy(inparammap, mcgroups);
358

359
360
  int randomseed = getProperty("MonteCarloRandomSeed");
  srand(randomseed);
361

362
363
364
  double temperature = getProperty("AnnealingTemperature");
  if (temperature < 1.0E-10)
    throw runtime_error("Annealing temperature is too low.");
365

366
367
368
  int maxiterations = getProperty("MonteCarloIterations");
  if (maxiterations <= 0)
    throw runtime_error("Max iteration cannot be 0 or less.");
369

370
371
372
373
374
375
  //    Book keeping
  map<string, Parameter> parammap;
  duplicateParameters(inparammap, parammap);
  // vector<pair<double, map<string, Parameter> > > bestresults;
  map<string, Parameter> bestresult;
  // size_t maxnumresults = 10;
376

377
378
379
380
  // 2. Set up parameters and get initial values
  m_bestChiSq = DBL_MAX;
  m_bestChiSqStep = -1;
  m_bestChiSqGroup = -1;
381

382
383
384
385
  double chisq0 = calculateFunction(parammap, vecY);
  double chisq0x = calculateFunctionError(m_positionFunc, m_dataWS, m_wsIndex);
  g_log.notice() << "[DBx510] Starting Chi^2 = " << chisq0 << " (homemade) "
                 << chisq0x << " (Levenber-marquadt)" << endl;
386

387
388
  bookKeepMCResult(parammap, chisq0, -1, -1,
                   bestresult); // bestresults, maxnumresults);
389

390
391
392
393
394
  // 3. Monte Carlo starts
  double chisqx = chisq0;
  int numtotalacceptance = 0;
  int numrecentacceptance = 0;
  int numrecentsteps = 0;
395

396
397
  map<string, Parameter> propparammap; // parameters with proposed value
  duplicateParameters(parammap, propparammap);
398

399
400
401
402
403
  for (int istep = 0; istep < maxiterations; ++istep) {
    for (int igroup = 0; igroup < static_cast<int>(mcgroups.size()); ++igroup) {
      // a) Propose value
      proposeNewValues(mcgroups[igroup], parammap, propparammap,
                       chisqx); // , prevbetterchi2);
404

405
406
      // b) Calcualte function and chi^2
      double propchisq = calculateFunction(propparammap, vecY);
407

408
409
410
411
412
      /*
      stringstream dbss;
      dbss << "[DBx541] New Chi^2 = " << propchisq << endl;
      vector<string> paramnames = m_positionFunc->getParameterNames();
      for (size_t i = 0; i < paramnames.size(); ++i)
413
      {
414
415
416
417
418
419
        string parname = paramnames[i];
        double curvalue = parammap[parname].value;
        double propvalue = propparammap[parname].value;
        dbss << parname << ":\t\t" << setw(20) << propvalue << "\t\t<-----\t\t"
      << curvalue << "\t Delta = "
             << curvalue-propvalue << endl;
420
      }
421
422
      g_log.notice(dbss.str());
      */
423

424
425
426
      // c) Determine to accept change
      bool acceptpropvalues =
          acceptOrDenyChange(propchisq, chisqx, temperature);
427

428
429
430
431
432
433
      // d) Change current chi^2, apply change, and book keep
      if (acceptpropvalues) {
        setFunctionParameterValues(m_positionFunc, propparammap);
        chisqx = propchisq;
        bookKeepMCResult(parammap, chisqx, istep, igroup,
                         bestresult); // s, maxnumresults);
434
435
      }

436
437
438
439
440
      // e) MC strategy control
      ++numtotalacceptance;
      ++numrecentacceptance;
      ++numrecentsteps;
    }
441

442
443
444
445
446
447
448
449
450
451
    // f) Annealing
    if (numrecentsteps >= 10) {
      double acceptratio = static_cast<double>(numrecentacceptance) /
                           static_cast<double>(numrecentsteps);
      if (acceptratio < 0.2) {
        // i) Low acceptance, need to raise temperature
        temperature *= 2.0;
      } else if (acceptratio >= 0.8) {
        // ii) Temperature too high to accept too much new change
        temperature /= 2.0;
452
453
      }

454
455
456
457
458
      // iii) Reset counters
      numrecentacceptance = 0;
      numrecentsteps = 0;
    }
  }
459

460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
  // 4. Apply the best result
  // sort(bestresults.begin(), bestresults.end());
  setFunctionParameterValues(m_positionFunc, bestresult);
  double chisqf = m_bestChiSq;

  g_log.warning() << "[DBx544] Best Chi^2 From MC = " << m_bestChiSq << endl;

  // 5. Use regular minimzer to try to get a better result
  string fitstatus;
  double fitchisq;
  bool goodfit =
      doFitFunction(m_positionFunc, m_dataWS, m_wsIndex,
                    "Levenberg-MarquardtMD", 1000, fitchisq, fitstatus);

  bool restoremcresult = false;
  if (goodfit) {
    map<string, Parameter> nullmap;
    fitchisq = calculateFunction(nullmap, vecY);
    if (fitchisq > chisqf) {
      // Fit is unable to achieve a better solution
      restoremcresult = true;
    } else {
      m_bestChiSq = fitchisq;
483
    }
484
485
486
487
  } else {
    // Fit is bad
    restoremcresult = true;
  }
488

489
  g_log.warning() << "[DBx545] Restore MC Result = " << restoremcresult << endl;
490

491
492
  if (restoremcresult) {
    setFunctionParameterValues(m_positionFunc, bestresult);
493
  }
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
  chisqf = m_bestChiSq;

  // 6. Final result
  double chisqfx = calculateFunctionError(m_positionFunc, m_dataWS, m_wsIndex);
  map<string, Parameter> emptymap;
  double chisqf0 = calculateFunction(emptymap, vecY);
  g_log.notice() << "Best Chi^2 (L-V) = " << chisqfx
                 << ", (homemade) = " << chisqf0 << endl;
  g_log.warning() << "Data Size = " << m_dataWS->readX(m_wsIndex).size()
                  << ", Number of parameters = "
                  << m_positionFunc->getParameterNames().size() << endl;

  return chisqf;
}

//----------------------------------------------------------------------------------------------
/** Propose new parameters
  *
  * @param mcgroup:     list of parameters to have new values proposed
  * @param currchisq:  present chi^2 (as a factor in step size)
  * @param curparammap: current parameter maps
  * @param newparammap: parameters map containing new/proposed value
  */
void RefinePowderInstrumentParameters3::proposeNewValues(
    vector<string> mcgroup, map<string, Parameter> &curparammap,
    map<string, Parameter> &newparammap, double currchisq) {
520
  for (auto paramname : mcgroup) {
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
    // random number between -1 and 1
    double randomnumber =
        2 * static_cast<double>(rand()) / static_cast<double>(RAND_MAX) - 1.0;

    // parameter information
    Parameter param = curparammap[paramname];
    double stepsize = m_dampingFactor * currchisq *
                      (param.curvalue * param.mcA1 + param.mcA0) *
                      randomnumber / m_bestChiSq;

    g_log.debug() << "Parameter " << paramname << " Step Size = " << stepsize
                  << " From " << param.mcA0 << ", " << param.mcA1 << ", "
                  << param.curvalue << ", " << m_dampingFactor << endl;

    // drunk walk or random walk
    double newvalue;
    // Random walk.  No preference on direction
    newvalue = param.curvalue + stepsize;
539

540
541
    /*
    if (m_walkStyle == RANDOMWALK)
542
543
    {

544
545
546
547
548
549
550
    }
    else if (m_walkStyle == DRUNKENWALK)
    {
      // Drunken walk.  Prefer to previous successful move direction
      int prevRightDirection;
      if (prevBetterRwp)
        prevRightDirection = 1;
551
      else
552
        prevRightDirection = -1;
553

554
555
      double randirint =
    static_cast<double>(rand())/static_cast<double>(RAND_MAX);
556

557
558
      // FIXME Here are some MAGIC numbers
      if (randirint < 0.1)
559
      {
560
561
562
        // Negative direction to previous direction
        stepsize =
    -1.0*fabs(stepsize)*static_cast<double>(param.movedirection*prevRightDirection);
563
      }
564
      else if (randirint < 0.4)
565
      {
566
567
        // No preferance
        stepsize = stepsize;
568
569
570
      }
      else
      {
571
572
573
        // Positive direction to previous direction
        stepsize =
    fabs(stepsize)*static_cast<double>(param.movedirection*prevRightDirection);
574
575
      }

576
      newvalue = param.value + stepsize;
577
578
579
    }
    else
    {
580
581
      newvalue = DBL_MAX;
      throw runtime_error("Unrecoganized walk style. ");
582
583
584
    }
    */

585
586
587
588
    // restriction
    if (param.nonnegative && newvalue < 0) {
      // If not allowed to be negative
      newvalue = fabs(newvalue);
589
590
    }

591
592
593
594
595
596
597
598
599
600
601
602
603
604
    // apply to new parameter map
    newparammap[paramname].curvalue = newvalue;

    // record some trace
    Parameter &p = curparammap[paramname];
    if (stepsize > 0) {
      p.movedirection = 1;
      ++p.numpositivemove;
    } else if (stepsize < 0) {
      p.movedirection = -1;
      ++p.numnegativemove;
    } else {
      p.movedirection = -1;
      ++p.numnomove;
605
    }
606
607
608
609
610
611
612
613
    p.sumstepsize += fabs(stepsize);
    if (fabs(stepsize) > p.maxabsstepsize)
      p.maxabsstepsize = fabs(stepsize);

    g_log.debug() << "[DBx257] " << paramname << "\t"
                  << "Proposed value = " << setw(15) << newvalue
                  << " (orig = " << param.curvalue << ",  step = " << stepsize
                  << "), totRwp = " << currchisq << endl;
614
615
  }

616
617
  return;
}
618

619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
//----------------------------------------------------------------------------------------------
/** Determine whether the proposed value should be accepted or denied
  *
  * @param curchisq:  present chi^2 (as a factor in step size)
  * @param newchisq:  new chi^2 (as a factor in step size)
  * @param temperature:  annealing temperature
  */
bool RefinePowderInstrumentParameters3::acceptOrDenyChange(double curchisq,
                                                           double newchisq,
                                                           double temperature) {
  bool accept;

  if (newchisq < curchisq) {
    // Lower Rwp.  Take the change
    accept = true;
  } else {
    // Higher Rwp. Take a chance to accept
    double dice = static_cast<double>(rand()) / static_cast<double>(RAND_MAX);
    double bar = exp(-(newchisq - curchisq) / (curchisq * temperature));
638
    accept = dice < bar;
639
640
  }

641
642
  return accept;
}
643

644
645
646
647
648
649
650
651
652
653
654
//----------------------------------------------------------------------------------------------
/** Book keep the best fitting result
  */
void RefinePowderInstrumentParameters3::bookKeepMCResult(
    map<string, Parameter> parammap, double chisq, int istep, int igroup,
    map<string, Parameter> &bestparammap)
// vector<pair<double, map<string, Parameter> > >& bestresults,
// size_t maxnumresults)
{
  // 1. Check whether input Chi^2 is the best Chi^2
  bool recordparameter = false;
655

656
657
658
659
  if (chisq < m_bestChiSq) {
    m_bestChiSq = chisq;
    m_bestChiSqStep = istep;
    m_bestChiSqGroup = igroup;
660

661
    recordparameter = true;
662
663
  }

664
  // 2. Record for the best parameters
665
  if (bestparammap.empty()) {
666
667
668
669
    // No record yet
    duplicateParameters(parammap, bestparammap);
  } else if (recordparameter) {
    // Replace the record
670
671
  }

672
673
674
675
676
  // 2. Determine whether to add this entry to records
  /*
  bool addentry = true;
  if (bestresults.size() >= maxnumresults && chisq > bestresults.back().first)
    addentry = false;
677

678
679
680
681
682
  // 3. Add entry
  if (addentry)
  {
    map<string, Parameter> storemap;
    duplicateParameters(parammap, storemap);
683
    bestresults.emplace_back(chisq, storemap);
684
685
686
    sort(bestresults.begin(), bestresults.end());
  }
  */
687

688
689
  return;
}
690

691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
//----------------------------------------------------------------------------------------------
/** Set up Monte Carlo random walk strategy
  */
void RefinePowderInstrumentParameters3::setupRandomWalkStrategy(
    map<string, Parameter> &parammap, vector<vector<string>> &mcgroups) {
  stringstream dboutss;
  dboutss << "Monte Carlo minimizer refines: ";

  // 1. Monte Carlo groups
  // a. Instrument gemetry
  vector<string> geomparams;
  addParameterToMCMinimize(geomparams, "Dtt1", parammap);
  addParameterToMCMinimize(geomparams, "Dtt1t", parammap);
  addParameterToMCMinimize(geomparams, "Dtt2t", parammap);
  addParameterToMCMinimize(geomparams, "Zero", parammap);
  addParameterToMCMinimize(geomparams, "Zerot", parammap);
  addParameterToMCMinimize(geomparams, "Width", parammap);
  addParameterToMCMinimize(geomparams, "Tcross", parammap);
  mcgroups.push_back(geomparams);

  dboutss << "Geometry parameters: ";
712
713
  for (auto &geomparam : geomparams)
    dboutss << geomparam << "\t\t";
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
  dboutss << endl;

  g_log.notice(dboutss.str());

  // 2. Dictionary for each parameter for non-negative, mcX0, mcX1
  parammap["Width"].mcA0 = 0.0;
  parammap["Width"].mcA1 = 1.0;
  parammap["Width"].nonnegative = true;

  parammap["Tcross"].mcA0 = 0.0;
  parammap["Tcross"].mcA1 = 1.0;
  parammap["Tcross"].nonnegative = true;

  parammap["Zero"].mcA0 = 5.0;
  parammap["Zero"].mcA1 = 0.0;
  parammap["Zero"].nonnegative = false;

  parammap["Zerot"].mcA0 = 5.0;
  parammap["Zerot"].mcA1 = 0.0;
  parammap["Zerot"].nonnegative = false;

  parammap["Dtt1"].mcA0 = 5.0;
  parammap["Dtt1"].mcA1 = 0.0;
  parammap["Dtt1"].nonnegative = true;

  parammap["Dtt1t"].mcA0 = 5.0;
  parammap["Dtt1t"].mcA1 = 0.0;
  parammap["Dtt1t"].nonnegative = true;

  parammap["Dtt2t"].mcA0 = 0.1;
  parammap["Dtt2t"].mcA1 = 1.0;
  parammap["Dtt2t"].nonnegative = false;

  // 3. Reset
  map<string, Parameter>::iterator mapiter;
  for (mapiter = parammap.begin(); mapiter != parammap.end(); ++mapiter) {
    mapiter->second.movedirection = 1;
    mapiter->second.sumstepsize = 0.0;
    mapiter->second.numpositivemove = 0;
    mapiter->second.numnegativemove = 0;
    mapiter->second.numnomove = 0;
    mapiter->second.maxabsstepsize = -0.0;
  }
757

758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
  return;
}

//----------------------------------------------------------------------------------------------
/** Add parameter (to a vector of string/name) for MC random walk
  * according to Fit in Parameter
  *
  * @param parnamesforMC: vector of parameter for MC minimizer
  * @param parname: name of parameter to check whether to put into refinement
  *list
  * @param parammap :: parammap
  */
void RefinePowderInstrumentParameters3::addParameterToMCMinimize(
    vector<string> &parnamesforMC, string parname,
    map<string, Parameter> parammap) {
  map<string, Parameter>::iterator pariter;
  pariter = parammap.find(parname);
  if (pariter == parammap.end()) {
    stringstream errss;
    errss << "Parameter " << parname
          << " does not exisit Le Bail function parameters. ";
    g_log.error(errss.str());
    throw runtime_error(errss.str());
781
782
  }

783
784
  if (pariter->second.fit)
    parnamesforMC.push_back(parname);
785

786
787
  return;
}
788

789
790
791
792
793
794
795
796
797
798
//----------------------------------------------------------------------------------------------
/** Implement parameter values, calculate function and its chi square.
  *
  * @param parammap:  if size = 0, there is no action to set function parameter.
  * @param vecY :: vecY
  * Return: chi^2
  */
double RefinePowderInstrumentParameters3::calculateFunction(
    map<string, Parameter> parammap, vector<double> &vecY) {
  // 1. Implement parameter values to m_positionFunc
799
  if (!parammap.empty())
800
    setFunctionParameterValues(m_positionFunc, parammap);
801

802
803
804
805
806
  // 2. Calculate
  const MantidVec &vecX = m_dataWS->readX(m_wsIndex);
  //    Check
  if (vecY.size() != vecX.size())
    throw runtime_error("vecY must be initialized with proper size!");
807

808
  m_positionFunc->function1D(vecY, vecX);
809

810
811
812
  // 3. Calcualte error
  double chisq = calculateFunctionChiSquare(vecY, m_dataWS->readY(m_wsIndex),
                                            m_dataWS->readE(m_wsIndex));
813

814
815
  return chisq;
}
816

817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
//----------------------------------------------------------------------------------------------
/** Calculate Chi^2
  */
double calculateFunctionChiSquare(const vector<double> &modelY,
                                  const vector<double> &dataY,
                                  const vector<double> &dataE) {
  // 1. Check
  if (modelY.size() != dataY.size() || dataY.size() != dataE.size())
    throw runtime_error("Input model, data and error have different size.");

  // 2. Calculation
  double chisq = 0.0;
  size_t numpts = modelY.size();
  for (size_t i = 0; i < numpts; ++i) {
    if (dataE[i] > 1.0E-5) {
      double temp = (modelY[i] - dataY[i]) / dataE[i];
      chisq += temp * temp;
834
835
836
    }
  }

837
838
  return chisq;
}
839

840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
//----------------------------------------------------------------------------------------------
/** Calculate Chi^2 of the a function with all parameters are fixed
  */
double RefinePowderInstrumentParameters3::calculateFunctionError(
    IFunction_sptr function, Workspace2D_sptr dataws, int wsindex) {
  // 1. Record the fitting information
  vector<string> parnames = function->getParameterNames();
  vector<bool> vecFix(parnames.size(), false);

  for (size_t i = 0; i < parnames.size(); ++i) {
    bool fixed = function->isFixed(i);
    vecFix[i] = fixed;
    if (!fixed)
      function->fix(i);
  }
855

856
857
858
  // 2. Fit with zero iteration
  double chi2;
  string fitstatus;
859
860
861
862
863
864
865
866
867
  const std::string minimizer = "Levenberg-MarquardtMD";
  bool fitOK =
      doFitFunction(function, dataws, wsindex, minimizer, 0, chi2, fitstatus);

  if (!fitOK) {
    g_log.warning() << "Fit by " << minimizer
                    << " with 0 iterations failed, with reason: " << fitstatus
                    << "\n";
  }
868

869
870
871
872
  // 3. Restore the fit/fix setup
  for (size_t i = 0; i < parnames.size(); ++i) {
    if (!vecFix[i])
      function->unfix(i);
873
874
  }

875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
  return chi2;
}

//----------------------------------------------------------------------------------------------
/** Fit a function by trying various minimizer or minimizer combination
  *
  * @param function :: an instance of a function to fit
  * @param dataws :: a workspace with the data
  * @param wsindex :: a histogram index
  * @param powerfit :: a flag to choose a robust algorithm to fit function
  *
  * Return: double chi2 of the final (best) solution.  If fitting fails, chi2
  *wil be maximum double
  */
double RefinePowderInstrumentParameters3::fitFunction(IFunction_sptr function,
                                                      Workspace2D_sptr dataws,
                                                      int wsindex,
                                                      bool powerfit) {
  // 1. Store original
  map<string, pair<double, double>> start_paramvaluemap, paramvaluemap1,
      paramvaluemap2, paramvaluemap3;
  storeFunctionParameterValue(function, start_paramvaluemap);

  // 2. Calculate starting chi^2
  double startchisq = calculateFunctionError(function, dataws, wsindex);
  g_log.notice() << "[DBx436] Starting Chi^2 = " << startchisq
                 << ", Power-Fit is " << powerfit << endl;

  // 3. Fitting
  int numiters;
  double final_chi2 = DBL_MAX;

  if (powerfit) {
    // a) Use Simplex to fit
    string minimizer = "Simplex";
    double chi2simplex;
    string fitstatussimplex;
    numiters = 10000;
    bool fitgood1 = doFitFunction(function, dataws, wsindex, minimizer,
                                  numiters, chi2simplex, fitstatussimplex);

    if (fitgood1)
      storeFunctionParameterValue(function, paramvaluemap1);
918
    else
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
      chi2simplex = DBL_MAX;

    // b) Continue Levenberg-Marquardt following Simplex
    minimizer = "Levenberg-MarquardtMD";
    double chi2lv2;
    string fitstatuslv2;
    numiters = 1000;
    bool fitgood2 = doFitFunction(function, dataws, wsindex, minimizer,
                                  numiters, chi2lv2, fitstatuslv2);
    if (fitgood2)
      storeFunctionParameterValue(function, paramvaluemap2);
    else
      chi2lv2 = DBL_MAX;

    // c) Fit by L.V. solely
    map<string, Parameter> tempparmap;
    restoreFunctionParameterValue(start_paramvaluemap, function, tempparmap);
    double chi2lv1;
    string fitstatuslv1;
    bool fitgood3 = doFitFunction(function, dataws, wsindex, minimizer,
                                  numiters, chi2lv1, fitstatuslv1);
    if (fitgood3)
      storeFunctionParameterValue(function, paramvaluemap3);
    else
      chi2lv1 = DBL_MAX;

    // 4. Compare best
    g_log.notice() << "Fit Result:  Chi2s: Simplex = " << chi2simplex << ", "
                   << "Levenberg 1 = " << chi2lv2
                   << ", Levenberg 2 = " << chi2lv1 << endl;

    if (fitgood1 || fitgood2 || fitgood3) {
      // At least one good fit
      if (fitgood1 && chi2simplex <= chi2lv2 && chi2simplex <= chi2lv1) {
        final_chi2 = chi2simplex;
        restoreFunctionParameterValue(paramvaluemap1, function,
                                      m_profileParameters);
      } else if (fitgood2 && chi2lv2 <= chi2lv1) {
        restoreFunctionParameterValue(paramvaluemap2, function,
                                      m_profileParameters);
        final_chi2 = chi2lv2;
      } else if (fitgood3) {
        final_chi2 = chi2lv1;
        restoreFunctionParameterValue(paramvaluemap3, function,
                                      m_profileParameters);
      } else {
        throw runtime_error("This situation is impossible to happen!");
      }
    } // END of Choosing Results
  } else {
    // 3B) Simple fit
    string minimizer = "Levenberg-MarquardtMD";
    string fitstatus;
    numiters = 1000;
    bool fitgood = doFitFunction(function, dataws, wsindex, minimizer, numiters,
                                 final_chi2, fitstatus);
    if (fitgood) {
      storeFunctionParameterValue(function, paramvaluemap1);
      restoreFunctionParameterValue(paramvaluemap1, function,
                                    m_profileParameters);
    } else {
      g_log.warning() << "Fit by " << minimizer
                      << " failed.  Reason: " << fitstatus << "\n";
982
983
    }
  }
984

985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
  return final_chi2;
}

//----------------------------------------------------------------------------------------------
/** Fit function
  * Minimizer: "Levenberg-MarquardtMD"/"Simplex"
 */
bool RefinePowderInstrumentParameters3::doFitFunction(
    IFunction_sptr function, Workspace2D_sptr dataws, int wsindex,
    string minimizer, int numiters, double &chi2, string &fitstatus) {
  // 0. Debug output
  stringstream outss;
  outss << "Fit function: " << m_positionFunc->asString() << endl
        << "Data To Fit: \n";
  for (size_t i = 0; i < dataws->readX(0).size(); ++i)
    outss << dataws->readX(wsindex)[i] << "\t\t" << dataws->readY(wsindex)[i]
For faster browsing, not all history is shown. View entire blame