IFunction.cpp 35.6 KB
Newer Older
1
2
3
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
4
#include "MantidAPI/Axis.h"
5
6
7
8
#include "MantidAPI/IFunction.h"
#include "MantidAPI/Jacobian.h"
#include "MantidAPI/IConstraint.h"
#include "MantidAPI/ParameterTie.h"
9
10
#include "MantidAPI/Expression.h"
#include "MantidAPI/ConstraintFactory.h"
11
#include "MantidAPI/FunctionFactory.h"
12
13
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/IFunctionWithLocation.h"
14
#include "MantidGeometry/Instrument.h"
15
16
17
18
#include "MantidGeometry/Instrument/ParameterMap.h"
#include "MantidGeometry/Instrument/Component.h"
#include "MantidGeometry/Instrument/DetectorGroup.h"
#include "MantidGeometry/Instrument/FitParameter.h"
19
20
21
#include "MantidGeometry/muParser_Silent.h"
#include "MantidKernel/Exception.h"
#include "MantidKernel/Logger.h"
22
#include "MantidKernel/UnitFactory.h"
23
#include "MantidKernel/MultiThreaded.h"
24
#include "MantidKernel/ProgressBase.h"
25
26
27

#include <boost/lexical_cast.hpp>

28
29
#include <Poco/StringTokenizer.h>

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

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

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

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

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

61
62
/**
 * Attach a progress reporter
63
64
 * @param reporter :: A pointer to a progress reporter that can be called during
 * function evaluation
65
 */
66
void IFunction::setProgressReporter(Kernel::ProgressBase *reporter) {
67
68
69
70
71
72
73
74
  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 = "")
 */
75
76
77
void IFunction::reportProgress(const std::string &msg) const {
  if (m_progReporter) {
    const_cast<Kernel::ProgressBase *>(m_progReporter)->report(msg);
78
79
80
81
82
  }
}

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

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

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

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

144
145
146
147
/** 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.
 */
148
void IFunction::removeTie(const std::string &parName) {
149
150
151
152
153
  size_t i = parameterIndex(parName);
  this->removeTie(i);
}

/**
154
155
 * Writes a string that can be used in Fit.IFunction to create a copy of this
 * IFunction
156
157
 * @return string representation of the function
 */
158
std::string IFunction::asString() const {
159
  std::ostringstream ostr;
160
  ostr << "name=" << this->name();
161
  // print the attributes
162
  std::vector<std::string> attr = this->getAttributeNames();
163
164
  for (const auto &attName : attr) {
    std::string attValue = this->getAttribute(attName).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
  for (const auto &expr : list) {
226
    IConstraint *c =
227
        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
const std::vector<std::string> IFunction::categories() const {
256
  Poco::StringTokenizer tokenizer(category(), categorySeparator(),
257
258
                                  Poco::StringTokenizer::TOK_TRIM |
                                      Poco::StringTokenizer::TOK_IGNORE_EMPTY);
259
  return std::vector<std::string>(tokenizer.begin(), tokenizer.end());
260
261
}

262
263
264
265
266
/**
 * Operator <<
 * @param ostr :: The output stream
 * @param f :: The IFunction
 */
267
std::ostream &operator<<(std::ostream &ostr, const IFunction &f) {
268
269
270
271
  ostr << f.asString();
  return ostr;
}

272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
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>";
  }
};
291
}
292

293
std::string IFunction::Attribute::type() const {
294
295
296
297
  AttType tmp;
  return apply(tmp);
}

298
299
300
301
302
303
namespace {
/**
 * Const attribute visitor returning the value of the attribute as a string
 */
class AttValue : public IFunction::ConstAttributeVisitor<std::string> {
public:
304
  explicit AttValue(bool quoteString = false)
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
      : 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());
331
    }
332
333
334
    res += ")";
    return res;
  }
335

336
337
338
339
private:
  /// Flag to quote a string value returned
  bool m_quoteString;
};
340
}
341

342
std::string IFunction::Attribute::value() const {
343
344
345
346
  AttValue tmp(m_quoteValue);
  return apply(tmp);
}

347
348
349
/**
 * Return the attribute as a string if it is a string.
 */
350
351
352
353
354
std::string IFunction::Attribute::asString() const {
  if (m_quoteValue)
    return asQuotedString();

  try {
355
    return boost::get<std::string>(m_data);
356
357
358
  } catch (...) {
    throw std::runtime_error("Trying to access a " + type() + " attribute "
                                                              "as string");
359
360
361
  }
}

362
363
364
/**
 * Return the attribute as a quoted string if it is a string.
 */
365
std::string IFunction::Attribute::asQuotedString() const {
366
367
  std::string attr;

368
  try {
369
    attr = boost::get<std::string>(m_data);
370
371
372
  } catch (...) {
    throw std::runtime_error("Trying to access a " + type() + " attribute "
                                                              "as string");
373
  }
374

375
376
  if (attr.empty())
    return "\"\"";
377

378
  std::string quoted(attr);
379
380
381
382
  if (*(attr.begin()) != '\"')
    quoted = "\"" + attr;
  if (*(quoted.end() - 1) != '\"')
    quoted += "\"";
383
384
385
386

  return quoted;
}

387
388
389
/**
 * Return the attribute as an unquoted string if it is a string.
 */
390
std::string IFunction::Attribute::asUnquotedString() const {
391
392
  std::string attr;

393
  try {
394
    attr = boost::get<std::string>(m_data);
395
396
397
  } catch (...) {
    throw std::runtime_error("Trying to access a " + type() + " attribute "
                                                              "as string");
398
399
  }
  std::string unquoted(attr);
400
401
402
403
404
405
406
  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);

407
408
409
  return unquoted;
}

410
411
412
/**
 * Return the attribute as an int if it is a int.
 */
413
414
int IFunction::Attribute::asInt() const {
  try {
415
    return boost::get<int>(m_data);
416
417
418
  } catch (...) {
    throw std::runtime_error("Trying to access a " + type() + " attribute "
                                                              "as int");
419
420
421
  }
}

422
423
424
/**
 * Return the attribute as a double if it is a double.
 */
425
426
double IFunction::Attribute::asDouble() const {
  try {
427
    return boost::get<double>(m_data);
428
429
430
  } catch (...) {
    throw std::runtime_error("Trying to access a " + type() + " attribute "
                                                              "as double");
431
432
433
  }
}

434
435
436
/**
 * Return the attribute as a bool if it is a bool.
 */
437
438
bool IFunction::Attribute::asBool() const {
  try {
439
    return boost::get<bool>(m_data);
440
441
442
  } catch (...) {
    throw std::runtime_error("Trying to access a " + type() + " attribute "
                                                              "as bool");
443
  }
444
445
446
447
448
}

/**
 * Return the attribute as a bool if it is a vector.
 */
449
450
451
452
453
454
455
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");
  }
456
457
}

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

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

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

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

/**
 * Sets new value if attribute is a vector. If the type is wrong
 * throws an exception
 * @param v :: The new value
 */
515
516
517
518
519
520
521
522
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");
  }
523
524
}

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

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;
567
    }
568
569
570
571
572
573
    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);
      }
574
    }
575
576
577
578
579
    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]);
580
    }
581
  }
582

583
584
585
private:
  mutable std::string m_value; ///<the value as a string
};
586
}
587
588
589
590

/** Set value from a string. Throws exception if the string has wrong format
 * @param str :: String representation of the new value
 */
591
void IFunction::Attribute::fromString(const std::string &str) {
592
593
594
595
  SetValue tmp(str);
  apply(tmp);
}

596
597
598
599
/// 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)) {
600
601
602
603
604
    throw std::runtime_error("Attempt to use an inactive parameter");
  }
  return getParameter(i);
}

605
606
607
608
/// 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)) {
609
610
    throw std::runtime_error("Attempt to use an inactive parameter");
  }
611
  setParameter(i, value);
612
613
}

614
615
616
617
/**
 * Returns the name of an active parameter.
 * @param i :: Index of a parameter. The parameter must be active.
 */
618
619
std::string IFunction::nameOfActive(size_t i) const {
  if (!isActive(i)) {
620
621
622
623
624
    throw std::runtime_error("Attempt to use an inactive parameter");
  }
  return parameterName(i);
}

625
626
627
628
/**
 * Returns the description of an active parameter.
 * @param i :: Index of a parameter. The parameter must be active.
 */
629
630
std::string IFunction::descriptionOfActive(size_t i) const {
  if (!isActive(i)) {
631
632
633
634
635
    throw std::runtime_error("Attempt to use an inactive parameter");
  }
  return parameterDescription(i);
}

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

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

654
  // PARALLEL_CRITICAL(numeric_deriv)
655
656
  {
    applyTies(); // just in case
657
    function(domain, minusStep);
658
  }
659

660
661
  for (size_t iP = 0; iP < nParam; iP++) {
    if (isActive(iP)) {
662
      const double val = activeParameter(iP);
663
      if (fabs(val) < cutoff) {
664
        step = epsilon;
665
666
      } else {
        step = val * stepPercentage;
667
      }
668

669
      double paramPstep = val + step;
670

671
      // PARALLEL_CRITICAL(numeric_deriv)
672
673
      {
        setActiveParameter(iP, paramPstep);
674
675
        applyTies();
        function(domain, plusStep);
676
677
        setActiveParameter(iP, val);
      }
678

679
      step = paramPstep - val;
680
681
682
683
      for (size_t i = 0; i < nData; i++) {
        jacobian.set(i, iP,
                     (plusStep.getCalculated(i) - minusStep.getCalculated(i)) /
                         step);
684
685
      }
    }
686
  }
687
688
}

689
690
691
692
693
/** 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
694
 */
695
696
697
void IFunction::setMatrixWorkspace(
    boost::shared_ptr<const API::MatrixWorkspace> workspace, size_t wi,
    double startX, double endX) {
698
699
700
  UNUSED_ARG(startX);
  UNUSED_ARG(endX);

701
702
  if (!workspace)
    return; // unset the workspace
703

704
  try {
705
706
707

    // check if parameter are specified in instrument definition file

708
    const Geometry::ParameterMap &paramMap = workspace->instrumentParameters();
709
710

    Geometry::IDetector_const_sptr det;
711
712
713
714
715
716
    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.
717
      Instrument_const_sptr inst = workspace->getInstrument();
718
719
720
      det = inst->getDetector(
          *workspace->getSpectrum(wi)->getDetectorIDs().begin());
    } else
721
      // Get the detector (single) at this workspace index
722
723
724
725
726
727
728
729
      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()) {
730
          // get FitParameter
731
732
733
734
735
736
737
738
739
          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);
740
            if (testWithLocation == nullptr ||
741
742
                (fitParam.getLookUpTable().containData() == false &&
                 fitParam.getFormula().compare("") == 0)) {
743
              setParameter(i, fitParam.getValue());
744
            } else {
745
              double centreValue = testWithLocation->centre();
746
747
748
749
750
751
752
753
754
755
756
757
758
759
              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";
760
761
762
763
                  }
                }
              }

764
765
766
767
768
769
770
771
772
773
774
775
              // 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;
776
              }
777
778
779

              double paramValue = fitParam.getValue(centreValue);

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

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

801
802
803
                if (!resultUnitStr.empty()) {
                  std::vector<std::string> allUnitStr =
                      Kernel::UnitFactory::Instance().getKeys();
804
805
                  for (auto &iUnit : allUnitStr) {
                    size_t found = resultUnitStr.find(iUnit);
806
                    if (found != std::string::npos) {
807
                      size_t len = iUnit.size();
808
                      std::stringstream readDouble;
809
                      Kernel::Unit_sptr unt =
810
                          Kernel::UnitFactory::Instance().create(iUnit);
811
                      readDouble << 1.0 / convertValue(1.0, unt, workspace, wi);
812
813
                      resultUnitStr.replace(found, len, readDouble.str());
                    }
814
                  } // end for
815

816
                  try {
817
818
                    mu::Parser p;
                    p.SetExpr(resultUnitStr);
819
820
821
822
                    g_log.debug() << "The FitParameter " << parameterName(i)
                                  << " = " << paramValue
                                  << " before result-unit convertion (using "
                                  << resultUnitStr << ")" << std::endl;
823
                    paramValue *= p.Eval();
824
825
826
827
828
829
830
831
832
833
834
                    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;
835
836
                  }
                } // end if
837
              } // end trying to convert result-unit from formula or y-unit for
838
              // lookuptable
839
840

              setParameter(i, paramValue);
841
            } // end of update parameter value
842

843
844
845
            // add tie if specified for this parameter in instrument definition
            // file
            if (fitParam.getTie().compare("")) {
846
847
848
849
850
              std::ostringstream str;
              str << getParameter(i);
              tie(parameterName(i), str.str());
            }

851
852
853
854
855
856
857
858
859
860
            // 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());
861
                  constraint->setPenaltyFactor(penalty);
862
863
864
                } catch (...) {
                  g_log.warning()
                      << "Can't set penalty factor for constraint\n";
865
866
867
868
869
870
871
872
                }
              }
              addConstraint(constraint);
            }
          }
        }
      }
    }
873
  } catch (...) {
874
  }
875
876
}

877
878
879
880
881
882
883
/** 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
884
 */
885
double IFunction::convertValue(double value, Kernel::Unit_sptr &outUnit,
886
                               boost::shared_ptr<const MatrixWorkspace> ws,
887
                               size_t wsIndex) const {
888
  // only required if formula or look-up-table different from ws unit
889
890
891
  const auto &wsUnit = ws->getAxis(0)->unit();
  if (outUnit->unitID().compare(wsUnit->unitID()) == 0)
    return value;
892
893
894
895

  // 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);
896
897
898
  if (wsUnit->quickConversion(*outUnit, factor, power)) {
    return factor * std::pow(value, power);
  } else {
899
900
901
902
    std::vector<double> singleValue(1, value);
    convertValue(singleValue, outUnit, ws, wsIndex);
    return singleValue.front();
  }
903
904
905
906
}

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

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

    std::vector<double> emptyVec;
    wsUnit->toTOF(values, emptyVec, l1, l2, twoTheta, emode, efixed, 0.0);
963
    outUnit->fromTOF(values, emptyVec, l1, l2, twoTheta, emode, efixed, 0.0);
964
  }
965
966
}

967
968
969
/**
* Returns the number of attributes associated with the function
*/
970
size_t IFunction::nAttributes() const { return m_attrs.size(); }
971
972

/// Check if attribute named exists
973
974
bool IFunction::hasAttribute(const std::string &name) const {
  return m_attrs.find(name) != m_attrs.end();
975
976
}

977
978
979
980
981
/**
  * Overload for const char* values.
  * @param attName :: Attribute name
  * @param value :: New attribute value to set
  */
982
983
984
985
void IFunction::setAttributeValue(const std::string &attName,
                                  const char *value) {
  std::string str(value);
  setAttributeValue(attName, str);
986
987
988
989
990
991
992
}

/**
  * Set string attribute by value. Make sure that quoted style doesn't change.
  * @param attName :: Attribute name
  * @param value :: New attribute value to set
  */
993
994
995
996
997
void IFunction::setAttributeValue(const std::string &attName,
                                  const std::string &value) {
  Attribute att = getAttribute(attName);
  att.setString(value);
  setAttribute(attName, att);
998
999
}

1000
/// Returns a list of attribute names
For faster browsing, not all history is shown. View entire blame