IFunction.cpp 35.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidKernel/Exception.h"
#include "MantidKernel/Logger.h"
#include "MantidAPI/IFunction.h"
#include "MantidAPI/Jacobian.h"
#include "MantidAPI/IConstraint.h"
#include "MantidAPI/ParameterTie.h"
10
11
#include "MantidAPI/Expression.h"
#include "MantidAPI/ConstraintFactory.h"
12
#include "MantidAPI/FunctionFactory.h"
13
14
15
16
17
18
19
20

#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/IFunctionWithLocation.h"
#include "MantidGeometry/Instrument/ParameterMap.h"
#include "MantidGeometry/Instrument/Component.h"
#include "MantidGeometry/Instrument/DetectorGroup.h"
#include "MantidGeometry/Instrument/FitParameter.h"
#include "MantidKernel/UnitFactory.h"
21
#include "MantidKernel/MultiThreaded.h"
22
#include "MantidKernel/ProgressBase.h"
23
#include "MantidGeometry/muParser_Silent.h"
24
25
26

#include <boost/lexical_cast.hpp>

27
28
#include <Poco/StringTokenizer.h>

29
#include <limits>
30
#include <sstream>
31
#include <algorithm>
32

33
34
35
namespace Mantid {
namespace API {
using namespace Geometry;
Gigg, Martyn Anthony's avatar
Gigg, Martyn Anthony committed
36

37
38
39
40
namespace {
/// static logger
Kernel::Logger g_log("IFunction");
}
41
42
43
44

/**
 * Destructor
 */
45
46
47
48
49
IFunction::~IFunction() {
  m_attrs.clear();
  if (m_handler) {
    delete m_handler;
    m_handler = NULL;
50
  }
51
}
52

53
54
55
56
57
58
/**
 * Virtual copy constructor
 */
boost::shared_ptr<IFunction> IFunction::clone() const {
  return FunctionFactory::Instance().createInitialized(this->asString());
}
59

60
61
/**
 * Attach a progress reporter
62
63
 * @param reporter :: A pointer to a progress reporter that can be called during
 * function evaluation
64
 */
65
void IFunction::setProgressReporter(Kernel::ProgressBase *reporter) {
66
67
68
69
70
71
72
73
  m_progReporter = reporter;
  m_progReporter->setNotifyStep(0.01);
}

/**
 * If a reporter object is set, reports progress with an optional message
 * @param msg :: A message to display (default = "")
 */
74
75
76
void IFunction::reportProgress(const std::string &msg) const {
  if (m_progReporter) {
    const_cast<Kernel::ProgressBase *>(m_progReporter)->report(msg);
77
78
79
80
81
  }
}

/**
 *
82
83
 * @returns true if a progress reporter is set & evalaution has been requested
 *to stop
84
 */
85
86
87
88
89
bool IFunction::cancellationRequestReceived() const {
  if (m_progReporter)
    return m_progReporter->hasCancellationBeenRequested();
  else
    return false;
90
91
}

92
/** Base class implementation calculates the derivatives numerically.
93
 * @param domain :: The domain of the function
94
95
 * @param jacobian :: A Jacobian matrix. It is expected to have dimensions of
 * domain.size() by nParams().
96
 */
97
98
void IFunction::functionDeriv(const FunctionDomain &domain,
                              Jacobian &jacobian) {
99
  calNumericalDeriv(domain, jacobian);
100
101
102
103
104
}

/**
 * Ties a parameter to other parameters
 * @param parName :: The name of the parameter to tie.
105
 * @param expr :: A math expression
106
107
 * @param isDefault :: Flag to mark as default the value of an object associated
 * with this reference: a tie or a constraint.
108
109
 * @return newly ties parameters
 */
110
111
112
ParameterTie *IFunction::tie(const std::string &parName,
                             const std::string &expr, bool isDefault) {
  ParameterTie *ti = new ParameterTie(this, parName, expr, isDefault);
113
114
115
  addTie(ti);
  this->fix(getParameterIndex(*ti));
  return ti;
116
117
}

118
119
/**
 * Add ties to the function.
120
121
 * @param ties :: Comma-separated list of name=value pairs where name is a
 *parameter name and value
122
 *  is a math expression tying the parameter to other parameters or a constant.
123
124
 * @param isDefault :: Flag to mark as default the value of an object associated
 *with this reference: a tie or a constraint.
125
 *
126
 */
127
void IFunction::addTies(const std::string &ties, bool isDefault) {
128
129
130
  Expression list;
  list.parse(ties);
  list.toList();
131
132
  for (auto t = list.begin(); t != list.end(); ++t) {
    if (t->name() == "=" && t->size() >= 2) {
133
134
      size_t n = t->size() - 1;
      const std::string value = (*t)[n].str();
135
      for (size_t i = n; i != 0;) {
136
        --i;
137
        this->tie((*t)[i].name(), value, isDefault);
138
139
140
141
142
      }
    }
  }
}

143
144
145
146
/** Removes the tie off a parameter. The parameter becomes active
 * This method can be used when constructing and editing the IFunction in a GUI
 * @param parName :: The name of the parameter which ties will be removed.
 */
147
void IFunction::removeTie(const std::string &parName) {
148
149
150
151
152
  size_t i = parameterIndex(parName);
  this->removeTie(i);
}

/**
153
154
 * Writes a string that can be used in Fit.IFunction to create a copy of this
 * IFunction
155
156
 * @return string representation of the function
 */
157
std::string IFunction::asString() const {
158
  std::ostringstream ostr;
159
  ostr << "name=" << this->name();
160
  // print the attributes
161
  std::vector<std::string> attr = this->getAttributeNames();
162
  for (size_t i = 0; i < attr.size(); i++) {
163
164
    std::string attName = attr[i];
    std::string attValue = this->getAttribute(attr[i]).value();
165
166
    if (!attValue.empty() && attValue != "\"\"") {
      ostr << ',' << attName << '=' << attValue;
167
168
    }
  }
169
  // print the parameters
170
171
172
173
  for (size_t i = 0; i < nParams(); i++) {
    const ParameterTie *tie = getTie(i);
    if (!tie || !tie->isDefault()) {
      ostr << ',' << parameterName(i) << '=' << getParameter(i);
174
    }
175
  }
176
  // collect non-default constraints
177
  std::string constraints;
178
179
180
  for (size_t i = 0; i < nParams(); i++) {
    const IConstraint *c = getConstraint(i);
    if (c && !c->isDefault()) {
181
      std::string tmp = c->asString();
182
183
      if (!tmp.empty()) {
        if (!constraints.empty()) {
184
185
186
187
188
189
          constraints += ",";
        }
        constraints += tmp;
      }
    }
  }
190
  // print constraints
191
  if (!constraints.empty()) {
192
193
    ostr << ",constraints=(" << constraints << ")";
  }
194
  // collect the non-default ties
195
  std::string ties;
196
197
198
  for (size_t i = 0; i < nParams(); i++) {
    const ParameterTie *tie = getTie(i);
    if (tie && !tie->isDefault()) {
199
      std::string tmp = tie->asString(this);
200
201
      if (!tmp.empty()) {
        if (!ties.empty()) {
202
203
204
205
206
207
          ties += ",";
        }
        ties += tmp;
      }
    }
  }
208
  // print the ties
209
  if (!ties.empty()) {
210
211
212
213
214
    ostr << ",ties=(" << ties << ")";
  }
  return ostr.str();
}

215
/** Add a list of constraints from a string
216
 * @param str :: A comma-separated list of constraint expressions.
217
218
 * @param isDefault :: Flag to mark as default the value of an object associated
 *with this reference.
219
 *
220
 */
221
void IFunction::addConstraints(const std::string &str, bool isDefault) {
222
223
224
  Expression list;
  list.parse(str);
  list.toList();
225
226
227
  for (auto expr = list.begin(); expr != list.end(); ++expr) {
    IConstraint *c =
        ConstraintFactory::Instance().createInitialized(this, *expr, isDefault);
228
229
230
231
232
233
234
    this->addConstraint(c);
  }
}

/**
 * Return a vector with all parameter names.
 */
235
std::vector<std::string> IFunction::getParameterNames() const {
236
  std::vector<std::string> out;
237
  for (size_t i = 0; i < nParams(); ++i) {
238
239
240
241
242
    out.push_back(parameterName(i));
  }
  return out;
}

243
244
245
/** Set a function handler
 * @param handler :: A new handler
 */
246
void IFunction::setHandler(FunctionHandler *handler) {
247
  m_handler = handler;
248
  if (handler && handler->function().get() != this) {
249
250
251
252
253
    throw std::runtime_error("Function handler points to a different function");
  }
  m_handler->init();
}

254
/// Function to return all of the categories that contain this function
255
256
const std::vector<std::string> IFunction::categories() const {
  std::vector<std::string> res;
257
  Poco::StringTokenizer tokenizer(category(), categorySeparator(),
258
259
                                  Poco::StringTokenizer::TOK_TRIM |
                                      Poco::StringTokenizer::TOK_IGNORE_EMPTY);
260
261
  Poco::StringTokenizer::Iterator h = tokenizer.begin();

262
  for (; h != tokenizer.end(); ++h) {
263
264
265
266
267
268
    res.push_back(*h);
  }

  return res;
}

269
270
271
272
273
/**
 * Operator <<
 * @param ostr :: The output stream
 * @param f :: The IFunction
 */
274
std::ostream &operator<<(std::ostream &ostr, const IFunction &f) {
275
276
277
278
  ostr << f.asString();
  return ostr;
}

279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
namespace {
/**
 * Const attribute visitor returning the type of the attribute
 */
class AttType : public IFunction::ConstAttributeVisitor<std::string> {
protected:
  /// Apply if string
  std::string apply(const std::string &) const { return "std::string"; }
  /// Apply if int
  std::string apply(const int &) const { return "int"; }
  /// Apply if double
  std::string apply(const double &) const { return "double"; }
  /// Apply if bool
  std::string apply(const bool &) const { return "bool"; }
  /// Apply if vector
  std::string apply(const std::vector<double> &) const {
    return "std::vector<double>";
  }
};
298
}
299

300
std::string IFunction::Attribute::type() const {
301
302
303
304
  AttType tmp;
  return apply(tmp);
}

305
306
307
308
309
310
namespace {
/**
 * Const attribute visitor returning the value of the attribute as a string
 */
class AttValue : public IFunction::ConstAttributeVisitor<std::string> {
public:
311
  explicit AttValue(bool quoteString = false)
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
      : IFunction::ConstAttributeVisitor<std::string>(),
        m_quoteString(quoteString) {}

protected:
  /// Apply if string
  std::string apply(const std::string &str) const {
    return (m_quoteString) ? std::string("\"" + str + "\"") : str;
  }
  /// Apply if int
  std::string apply(const int &i) const {
    return boost::lexical_cast<std::string>(i);
  }
  /// Apply if double
  std::string apply(const double &d) const {
    return boost::lexical_cast<std::string>(d);
  }
  /// Apply if bool
  std::string apply(const bool &b) const { return b ? "true" : "false"; }
  /// Apply if vector
  std::string apply(const std::vector<double> &v) const {
    std::string res = "(";
    if (v.size() > 0) {
      for (size_t i = 0; i < v.size() - 1; ++i) {
        res += boost::lexical_cast<std::string>(v[i]) + ",";
      }
      res += boost::lexical_cast<std::string>(v.back());
338
    }
339
340
341
    res += ")";
    return res;
  }
342

343
344
345
346
private:
  /// Flag to quote a string value returned
  bool m_quoteString;
};
347
}
348

349
std::string IFunction::Attribute::value() const {
350
351
352
353
  AttValue tmp(m_quoteValue);
  return apply(tmp);
}

354
355
356
/**
 * Return the attribute as a string if it is a string.
 */
357
358
359
360
361
std::string IFunction::Attribute::asString() const {
  if (m_quoteValue)
    return asQuotedString();

  try {
362
    return boost::get<std::string>(m_data);
363
364
365
  } catch (...) {
    throw std::runtime_error("Trying to access a " + type() + " attribute "
                                                              "as string");
366
367
368
  }
}

369
370
371
/**
 * Return the attribute as a quoted string if it is a string.
 */
372
std::string IFunction::Attribute::asQuotedString() const {
373
374
  std::string attr;

375
  try {
376
    attr = boost::get<std::string>(m_data);
377
378
379
  } catch (...) {
    throw std::runtime_error("Trying to access a " + type() + " attribute "
                                                              "as string");
380
  }
381

382
383
  if (attr.empty())
    return "\"\"";
384

385
  std::string quoted(attr);
386
387
388
389
  if (*(attr.begin()) != '\"')
    quoted = "\"" + attr;
  if (*(quoted.end() - 1) != '\"')
    quoted += "\"";
390
391
392
393

  return quoted;
}

394
395
396
/**
 * Return the attribute as an unquoted string if it is a string.
 */
397
std::string IFunction::Attribute::asUnquotedString() const {
398
399
  std::string attr;

400
  try {
401
    attr = boost::get<std::string>(m_data);
402
403
404
  } catch (...) {
    throw std::runtime_error("Trying to access a " + type() + " attribute "
                                                              "as string");
405
406
  }
  std::string unquoted(attr);
407
408
409
410
411
412
413
  if (attr.empty())
    return "";
  if (*(attr.begin()) == '\"')
    unquoted = std::string(attr.begin() + 1, attr.end() - 1);
  if (*(unquoted.end() - 1) == '\"')
    unquoted = std::string(unquoted.begin(), unquoted.end() - 1);

414
415
416
  return unquoted;
}

417
418
419
/**
 * Return the attribute as an int if it is a int.
 */
420
421
int IFunction::Attribute::asInt() const {
  try {
422
    return boost::get<int>(m_data);
423
424
425
  } catch (...) {
    throw std::runtime_error("Trying to access a " + type() + " attribute "
                                                              "as int");
426
427
428
  }
}

429
430
431
/**
 * Return the attribute as a double if it is a double.
 */
432
433
double IFunction::Attribute::asDouble() const {
  try {
434
    return boost::get<double>(m_data);
435
436
437
  } catch (...) {
    throw std::runtime_error("Trying to access a " + type() + " attribute "
                                                              "as double");
438
439
440
  }
}

441
442
443
/**
 * Return the attribute as a bool if it is a bool.
 */
444
445
bool IFunction::Attribute::asBool() const {
  try {
446
    return boost::get<bool>(m_data);
447
448
449
  } catch (...) {
    throw std::runtime_error("Trying to access a " + type() + " attribute "
                                                              "as bool");
450
  }
451
452
453
454
455
}

/**
 * Return the attribute as a bool if it is a vector.
 */
456
457
458
459
460
461
462
std::vector<double> IFunction::Attribute::asVector() const {
  try {
    return boost::get<std::vector<double>>(m_data);
  } catch (...) {
    throw std::runtime_error("Trying to access a " + type() + " attribute "
                                                              "as vector");
  }
463
464
}

465
/** Sets new value if attribute is a string. If the type is wrong
466
467
468
 * throws an exception
 * @param str :: The new value
 */
469
470
void IFunction::Attribute::setString(const std::string &str) {
  try {
471
    boost::get<std::string>(m_data) = str;
472
473
474
  } catch (...) {
    throw std::runtime_error("Trying to access a " + type() + " attribute "
                                                              "as string");
475
476
477
  }
}

478
/** Sets new value if attribute is a double. If the type is wrong
479
480
481
 * throws an exception
 * @param d :: The new value
 */
482
483
void IFunction::Attribute::setDouble(const double &d) {
  try {
484
    boost::get<double>(m_data) = d;
485
486
487
  } catch (...) {
    throw std::runtime_error("Trying to access a " + type() + " attribute "
                                                              "as double");
488
489
490
  }
}

491
/** Sets new value if attribute is an int. If the type is wrong
492
493
494
 * throws an exception
 * @param i :: The new value
 */
495
496
void IFunction::Attribute::setInt(const int &i) {
  try {
497
    boost::get<int>(m_data) = i;
498
499
500
  } catch (...) {
    throw std::runtime_error("Trying to access a " + type() + " attribute "
                                                              "as int");
501
502
503
  }
}

504
/** Sets new value if attribute is an bool. If the type is wrong
505
506
 * throws an exception
 * @param b :: The new value
507
 */
508
509
void IFunction::Attribute::setBool(const bool &b) {
  try {
510
    boost::get<bool>(m_data) = b;
511
512
513
  } catch (...) {
    throw std::runtime_error("Trying to access a " + type() + " attribute "
                                                              "as bool");
514
  }
515
516
517
518
519
520
521
}

/**
 * Sets new value if attribute is a vector. If the type is wrong
 * throws an exception
 * @param v :: The new value
 */
522
523
524
525
526
527
528
529
void IFunction::Attribute::setVector(const std::vector<double> &v) {
  try {
    auto &value = boost::get<std::vector<double>>(m_data);
    value.assign(v.begin(), v.end());
  } catch (...) {
    throw std::runtime_error("Trying to access a " + type() + " attribute "
                                                              "as vector");
  }
530
531
}

532
533
534
535
536
537
namespace {
/**
 * Attribute visitor setting new value to an attribute
 */
class SetValue : public IFunction::AttributeVisitor<> {
public:
538
  /**
539
540
   * Constructor
   * @param value :: The value to set
541
   */
542
  explicit SetValue(const std::string &value) : m_value(value) {}
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573

protected:
  /// Apply if string
  void apply(std::string &str) const { str = m_value; }
  /// Apply if int
  void apply(int &i) const {
    std::istringstream istr(m_value + " ");
    istr >> i;
    if (!istr.good())
      throw std::invalid_argument("Failed to set int attribute "
                                  "from string " +
                                  m_value);
  }
  /// Apply if double
  void apply(double &d) const {
    std::istringstream istr(m_value + " ");
    istr >> d;
    if (!istr.good())
      throw std::invalid_argument("Failed to set double attribute "
                                  "from string " +
                                  m_value);
  }
  /// Apply if bool
  void apply(bool &b) const {
    b = (m_value == "true" || m_value == "TRUE" || m_value == "1");
  }
  /// Apply if vector
  void apply(std::vector<double> &v) const {
    if (m_value.empty()) {
      v.clear();
      return;
574
    }
575
576
577
578
579
580
    if (m_value.size() > 2) {
      // check if the value is in barckets (...)
      if (m_value[0] == '(' && m_value[m_value.size() - 1] == ')') {
        m_value.erase(0, 1);
        m_value.erase(m_value.size() - 1);
      }
581
    }
582
583
584
585
586
    Poco::StringTokenizer tokenizer(m_value, ",",
                                    Poco::StringTokenizer::TOK_TRIM);
    v.resize(tokenizer.count());
    for (size_t i = 0; i < v.size(); ++i) {
      v[i] = boost::lexical_cast<double>(tokenizer[i]);
587
    }
588
  }
589

590
591
592
private:
  mutable std::string m_value; ///<the value as a string
};
593
}
594
595
596
597

/** Set value from a string. Throws exception if the string has wrong format
 * @param str :: String representation of the new value
 */
598
void IFunction::Attribute::fromString(const std::string &str) {
599
600
601
602
  SetValue tmp(str);
  apply(tmp);
}

603
604
605
606
/// Value of i-th active parameter. Override this method to make fitted
/// parameters different from the declared
double IFunction::activeParameter(size_t i) const {
  if (!isActive(i)) {
607
608
609
610
611
    throw std::runtime_error("Attempt to use an inactive parameter");
  }
  return getParameter(i);
}

612
613
614
615
/// Set new value of i-th active parameter. Override this method to make fitted
/// parameters different from the declared
void IFunction::setActiveParameter(size_t i, double value) {
  if (!isActive(i)) {
616
617
    throw std::runtime_error("Attempt to use an inactive parameter");
  }
618
  setParameter(i, value);
619
620
}

621
622
623
624
/**
 * Returns the name of an active parameter.
 * @param i :: Index of a parameter. The parameter must be active.
 */
625
626
std::string IFunction::nameOfActive(size_t i) const {
  if (!isActive(i)) {
627
628
629
630
631
    throw std::runtime_error("Attempt to use an inactive parameter");
  }
  return parameterName(i);
}

632
633
634
635
/**
 * Returns the description of an active parameter.
 * @param i :: Index of a parameter. The parameter must be active.
 */
636
637
std::string IFunction::descriptionOfActive(size_t i) const {
  if (!isActive(i)) {
638
639
640
641
642
    throw std::runtime_error("Attempt to use an inactive parameter");
  }
  return parameterDescription(i);
}

643
/** Calculate numerical derivatives.
644
 * @param domain :: The domain of the function
645
646
 * @param jacobian :: A Jacobian matrix. It is expected to have dimensions of
 * domain.size() by nParams().
647
 */
648
649
void IFunction::calNumericalDeriv(const FunctionDomain &domain,
                                  Jacobian &jacobian) {
650
651
652
  const double minDouble = std::numeric_limits<double>::min();
  const double epsilon = std::numeric_limits<double>::epsilon() * 100;
  double stepPercentage = 0.001; // step percentage
653
654
  double step;                   // real step
  double cutoff = 100.0 * minDouble / stepPercentage;
655
656
657
658
659
660
  size_t nParam = nParams();
  size_t nData = domain.size();

  FunctionValues minusStep(domain);
  FunctionValues plusStep(domain);

661
  // PARALLEL_CRITICAL(numeric_deriv)
662
663
  {
    applyTies(); // just in case
664
    function(domain, minusStep);
665
  }
666

667
668
  for (size_t iP = 0; iP < nParam; iP++) {
    if (isActive(iP)) {
669
      const double val = activeParameter(iP);
670
      if (fabs(val) < cutoff) {
671
        step = epsilon;
672
673
      } else {
        step = val * stepPercentage;
674
      }
675

676
      double paramPstep = val + step;
677

678
      // PARALLEL_CRITICAL(numeric_deriv)
679
680
      {
        setActiveParameter(iP, paramPstep);
681
682
        applyTies();
        function(domain, plusStep);
683
684
        setActiveParameter(iP, val);
      }
685

686
      step = paramPstep - val;
687
688
689
690
      for (size_t i = 0; i < nData; i++) {
        jacobian.set(i, iP,
                     (plusStep.getCalculated(i) - minusStep.getCalculated(i)) /
                         step);
691
692
      }
    }
693
  }
694
695
}

696
697
698
699
700
/** Initialize the function providing it the workspace
 * @param workspace :: The workspace to set
 * @param wi :: The workspace index
 * @param startX :: The lower bin index
 * @param endX :: The upper bin index
701
 */
702
703
704
void IFunction::setMatrixWorkspace(
    boost::shared_ptr<const API::MatrixWorkspace> workspace, size_t wi,
    double startX, double endX) {
705
706
707
  UNUSED_ARG(startX);
  UNUSED_ARG(endX);

708
709
  if (!workspace)
    return; // unset the workspace
710

711
  try {
712
713
714

    // check if parameter are specified in instrument definition file

715
    const Geometry::ParameterMap &paramMap = workspace->instrumentParameters();
716
717

    Geometry::IDetector_const_sptr det;
718
719
720
721
722
723
    size_t numDetectors = workspace->getSpectrum(wi)->getDetectorIDs().size();
    if (numDetectors > 1) {
      // If several detectors are on this workspace index, just use the ID of
      // the first detector
      // Note JZ oct 2011 - I'm not sure why the code uses the first detector
      // and not the group. Ask Roman.
724
      Instrument_const_sptr inst = workspace->getInstrument();
725
726
727
      det = inst->getDetector(
          *workspace->getSpectrum(wi)->getDetectorIDs().begin());
    } else
728
      // Get the detector (single) at this workspace index
729
730
731
732
733
734
735
736
      det = workspace->getDetector(wi);
    ;

    for (size_t i = 0; i < nParams(); i++) {
      if (!isExplicitlySet(i)) {
        Geometry::Parameter_sptr param =
            paramMap.getRecursive(&(*det), parameterName(i), "fitting");
        if (param != Geometry::Parameter_sptr()) {
737
          // get FitParameter
738
739
740
741
742
743
744
745
746
747
748
749
          const Geometry::FitParameter &fitParam =
              param->value<Geometry::FitParameter>();

          // check first if this parameter is actually specified for this
          // function
          if (name().compare(fitParam.getFunction()) == 0) {
            // update value
            IFunctionWithLocation *testWithLocation =
                dynamic_cast<IFunctionWithLocation *>(this);
            if (testWithLocation == NULL ||
                (fitParam.getLookUpTable().containData() == false &&
                 fitParam.getFormula().compare("") == 0)) {
750
              setParameter(i, fitParam.getValue());
751
            } else {
752
              double centreValue = testWithLocation->centre();
753
754
755
756
757
758
759
760
761
762
763
764
765
766
              Kernel::Unit_sptr centreUnit; // unit of value used in formula or
                                            // to look up value in lookup table
              if (fitParam.getFormula().compare("") == 0)
                centreUnit = fitParam.getLookUpTable().getXUnit(); // from table
              else {
                if (!fitParam.getFormulaUnit().empty()) {
                  try {
                    centreUnit = Kernel::UnitFactory::Instance().create(
                        fitParam.getFormulaUnit()); // from formula
                  } catch (...) {
                    g_log.warning()
                        << fitParam.getFormulaUnit()
                        << " Is not an recognised formula unit for parameter "
                        << fitParam.getName() << "\n";
767
768
769
770
                  }
                }
              }

771
772
773
774
775
776
777
778
779
780
781
782
              // if unit specified convert centre value to unit required by
              // formula or look-up-table
              if (centreUnit) {
                g_log.debug()
                    << "For FitParameter " << parameterName(i)
                    << " centre of peak before any unit convertion is "
                    << centreValue << std::endl;
                centreValue =
                    convertValue(centreValue, centreUnit, workspace, wi);
                g_log.debug() << "For FitParameter " << parameterName(i)
                              << " centre of peak after any unit convertion is "
                              << centreValue << std::endl;
783
              }
784
785
786

              double paramValue = fitParam.getValue(centreValue);

787
788
              // this returned param value by a formula or a look-up-table may
              // have
789
              // a unit of its own. If set convert param value
790
791
792
              // See section 'Using fitting parameters in
              // www.mantidproject.org/IDF
              if (fitParam.getFormula().compare("") == 0) {
793
                // so from look up table
794
795
796
797
798
                Kernel::Unit_sptr resultUnit =
                    fitParam.getLookUpTable().getYUnit(); // from table
                g_log.debug() << "The FitParameter " << parameterName(i)
                              << " = " << paramValue
                              << " before y-unit convertion" << std::endl;
799
                paramValue /= convertValue(1.0, resultUnit, workspace, wi);
800
801
802
803
804
                g_log.debug() << "The FitParameter " << parameterName(i)
                              << " = " << paramValue
                              << " after y-unit convertion" << std::endl;
              } else {
                // so from formula
805
806
807

                std::string resultUnitStr = fitParam.getResultUnit();

808
809
810
811
                if (!resultUnitStr.empty()) {
                  std::vector<std::string> allUnitStr =
                      Kernel::UnitFactory::Instance().getKeys();
                  for (unsigned iUnit = 0; iUnit < allUnitStr.size(); iUnit++) {
812
                    size_t found = resultUnitStr.find(allUnitStr[iUnit]);
813
                    if (found != std::string::npos) {
814
815
                      size_t len = allUnitStr[iUnit].size();
                      std::stringstream readDouble;
816
817
818
819
                      Kernel::Unit_sptr unt =
                          Kernel::UnitFactory::Instance().create(
                              allUnitStr[iUnit]);
                      readDouble << 1.0 / convertValue(1.0, unt, workspace, wi);
820
821
                      resultUnitStr.replace(found, len, readDouble.str());
                    }
822
                  } // end for
823

824
                  try {
825
826
                    mu::Parser p;
                    p.SetExpr(resultUnitStr);
827
828
829
830
                    g_log.debug() << "The FitParameter " << parameterName(i)
                                  << " = " << paramValue
                                  << " before result-unit convertion (using "
                                  << resultUnitStr << ")" << std::endl;
831
                    paramValue *= p.Eval();
832
833
834
835
836
837
838
839
840
841
842
                    g_log.debug() << "The FitParameter " << parameterName(i)
                                  << " = " << paramValue
                                  << " after result-unit convertion"
                                  << std::endl;
                  } catch (mu::Parser::exception_type &e) {
                    g_log.error()
                        << "Cannot convert formula unit to workspace unit"
                        << " Formula unit which cannot be passed is "
                        << resultUnitStr
                        << ". Muparser error message is: " << e.GetMsg()
                        << std::endl;
843
844
                  }
                } // end if
845
              } // end trying to convert result-unit from formula or y-unit for
846
              // lookuptable
847
848

              setParameter(i, paramValue);
849
            } // end of update parameter value
850

851
852
853
            // add tie if specified for this parameter in instrument definition
            // file
            if (fitParam.getTie().compare("")) {
854
855
856
857
858
              std::ostringstream str;
              str << getParameter(i);
              tie(parameterName(i), str.str());
            }

859
860
861
862
863
864
865
866
867
868
            // add constraint if specified for this parameter in instrument
            // definition file
            if (fitParam.getConstraint().compare("")) {
              IConstraint *constraint =
                  ConstraintFactory::Instance().createInitialized(
                      this, fitParam.getConstraint());
              if (fitParam.getConstraintPenaltyFactor().compare("")) {
                try {
                  double penalty =
                      atof(fitParam.getConstraintPenaltyFactor().c_str());
869
                  constraint->setPenaltyFactor(penalty);
870
871
872
                } catch (...) {
                  g_log.warning()
                      << "Can't set penalty factor for constraint\n";
873
874
875
876
877
878
879
880
                }
              }
              addConstraint(constraint);
            }
          }
        }
      }
    }
881
  } catch (...) {
882
  }
883
884
}

885
886
887
888
889
890
891
/** Convert a value from unit defined in workspace (ws) to outUnit
 *
 *  @param value ::   assumed to be in unit of workspace
 *  @param outUnit ::  unit to convert to
 *  @param ws ::      workspace
 *  @param wsIndex :: workspace index
 *  @return converted value
892
 */
893
double IFunction::convertValue(double value, Kernel::Unit_sptr &outUnit,
894
                               boost::shared_ptr<const MatrixWorkspace> ws,
895
                               size_t wsIndex) const {
896
  // only required if formula or look-up-table different from ws unit
897
898
899
  const auto &wsUnit = ws->getAxis(0)->unit();
  if (outUnit->unitID().compare(wsUnit->unitID()) == 0)
    return value;
900
901
902
903

  // first check if it is possible to do a quick conversion and convert
  // slight duplication to below to avoid instantiating vector unless necessary
  double factor(0.0), power(0.0);
904
905
906
  if (wsUnit->quickConversion(*outUnit, factor, power)) {
    return factor * std::pow(value, power);
  } else {
907
908
909
910
    std::vector<double> singleValue(1, value);
    convertValue(singleValue, outUnit, ws, wsIndex);
    return singleValue.front();
  }
911
912
913
914
}

/** Convert values from unit defined in workspace (ws) to outUnit
 *
915
 *  @param values ::   As input: assumed to be in unit of workspace.
916
917
918
919
920
 *                  As output: in unit of outUnit
 *  @param outUnit ::  unit to convert to
 *  @param ws ::      workspace
 *  @param wsIndex :: workspace index
 */
921
922
923
924
void IFunction::convertValue(std::vector<double> &values,
                             Kernel::Unit_sptr &outUnit,
                             boost::shared_ptr<const MatrixWorkspace> ws,
                             size_t wsIndex) const {
925
  // only required if  formula or look-up-table different from ws unit
926
927
928
  const auto &wsUnit = ws->getAxis(0)->unit();
  if (outUnit->unitID().compare(wsUnit->unitID()) == 0)
    return;
929

930
931
  // first check if it is possible to do a quick conversion convert
  double factor, power;
932
  if (wsUnit->quickConversion(*outUnit, factor, power)) {
933
    auto iend = values.end();
934
    for (auto itr = values.begin(); itr != iend; ++itr)
935
      (*itr) = factor * std::pow(*itr, power);
936
  } else {
937
938
    // Get l1, l2 and theta  (see also RemoveBins.calculateDetectorPosition())
    Instrument_const_sptr instrument = ws->getInstrument();
939
    Geometry::IComponent_const_sptr sample = instrument->getSample();
940
941
942
    if (sample == NULL) {
      g_log.error()
          << "No sample defined instrument. Cannot convert units for function\n"
943
944
          << "Ignore convertion.";
      return;
945
    }
946
947
948
    double l1 = instrument->getSource()->getDistance(*sample);
    Geometry::IDetector_const_sptr det = ws->getDetector(wsIndex);
    double l2(-1.0), twoTheta(0.0);
949
    if (!det->isMonitor()) {
950
951
      l2 = det->getDistance(*sample);
      twoTheta = ws->detectorTwoTheta(det);
952
953
    } else // If this is a monitor then make l1+l2 = source-detector distance
           // and twoTheta=0
954
955
956
957
    {
      l2 = det->getDistance(*(instrument->getSource()));
      l2 = l2 - l1;
      twoTheta = 0.0;
958
    }
959
960
    int emode = static_cast<int>(ws->getEMode());
    double efixed(0.0);
961
    try {
962
      efixed = ws->getEFixed(det);
963
    } catch (std::exception &) {
964
965
966
      // assume elastic
      efixed = 0.0;
      emode = 0;
967
    }
968
969
970

    std::vector<double> emptyVec;
    wsUnit->toTOF(values, emptyVec, l1, l2, twoTheta, emode, efixed, 0.0);
971
    outUnit->fromTOF(values, emptyVec, l1, l2, twoTheta, emode, efixed, 0.0);
972
  }
973
974
}

975
976
977
/**
* Returns the number of attributes associated with the function
*/
978
size_t IFunction::nAttributes() const { return m_attrs.size(); }
979
980

/// Check if attribute named exists
981
982
bool IFunction::hasAttribute(const std::string &name) const {
  return m_attrs.find(name) != m_attrs.end();
983
984
}

985
986
987
988
989
/**
  * Overload for const char* values.
  * @param attName :: Attribute name
  * @param value :: New attribute value to set
  */
990
991
992
993
void IFunction::setAttributeValue(const std::string &attName,
                                  const char *value) {
  std::string str(value);
  setAttributeValue(attName, str);
994
995
996
997
998
999
1000
}

/**
  * Set string attribute by value. Make sure that quoted style doesn't change.
  * @param attName :: Attribute name
  * @param value :: New attribute value to set
  */
For faster browsing, not all history is shown. View entire blame