IFitFunction.cpp 12.4 KB
Newer Older
1
2
3
4
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidKernel/Exception.h"
5
#include "MantidKernel/Logger.h"
6
#include "MantidAPI/IFitFunction.h"
7
#include "MantidAPI/Jacobian.h"
8
9
10
#include "MantidAPI/IConstraint.h"
#include "MantidAPI/ParameterTie.h"

Campbell, Stuart's avatar
Campbell, Stuart committed
11
#include <boost/lexical_cast.hpp>
12
13
14
15
16
17
18
19
20
21

#include <sstream>
#include <iostream> 

namespace Mantid
{
namespace API
{
  
/** Base class implementation of derivative IFitFunction throws error. This is to check if such a function is provided
22
    by derivative class. In the derived classes this method must return the derivatives of the function
23
    with respect to the fit parameters. If this method is not reimplemented the derivative free simplex minimization
24
    algorithm is used or the derivatives are computed numerically.
25
 * @param out :: Pointer to a Jacobian matrix. If it is NULL the method is called in order to check whether it's implemented or not.
26
 *      If the derivatives are implemented the method must simply return, otherwise it must throw Kernel::Exception::NotImplementedError.
27
28
29
 */
void IFitFunction::functionDeriv(Jacobian* out)
{
Gigg, Martyn Anthony's avatar
Gigg, Martyn Anthony committed
30
  UNUSED_ARG(out);
31
  throw ("No derivative IFitFunction provided");
32
33
34
35
}



36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
///** Update active parameters. Ties are applied.
// *  @param in :: Pointer to an array with active parameters values. Must be at least nActive() doubles long.
// */
//void IFitFunction::updateActive(const double* in)
//{
//  if (in)
//    for(size_t i=0;i<nActive();i++)
//    {
//      setActiveParameter(i,in[i]);
//    }
//  applyTies();
//}
//
///**
// * Sets active parameter i to value. Ties are not applied.
// * @param i :: The index of active parameter to set
// * @param value :: The new value for the parameter
// */
//void IFitFunction::setActiveParameter(size_t i,double value)
//{
//  size_t j = indexOfActive(i);
//  setParameter(j,value,false);
//}
//
//double IFitFunction::activeParameter(size_t i)const
//{
//  size_t j = indexOfActive(i);
//  return getParameter(j);
//}

///** Create a new tie. IFitFunctions can have their own types of ties.
// * @param parName :: The parameter name for this tie
// * @return a new parameter tie
// */
//ParameterTie* IFitFunction::createTie(const std::string& parName)
//{
//  return new ParameterTie(this,parName);
//}
//
///**
// * Ties a parameter to other parameters
// * @param parName :: The name of the parameter to tie.
// * @param expr ::    A math expression 
// * @return newly ties parameters
// */
//ParameterTie* IFitFunction::tie(const std::string& parName,const std::string& expr)
//{
//  ParameterTie* tie = this->createTie(parName);
//  size_t i = getParameterIndex(*tie);
//  //if (i < 0)
//  //{
//  //  delete tie;
//  //  throw std::logic_error("Parameter "+parName+" was not found.");
//  //}
//
//  tie->set(expr);
//  addTie(tie);
//  this->removeActive(i);
//  return tie;
//}
//
///** Removes the tie off a parameter. The parameter becomes active
// * This method can be used when constructing and editing the IFitFunction in a GUI
// * @param parName :: The name of the parameter which ties will be removed.
// */
//void IFitFunction::removeTie(const std::string& parName)
//{
//  size_t i = parameterIndex(parName);
//  this->removeTie(i);
//}
106

107
108
109
/**
  * If any of the parameters do not satisfy a constraint penalty values will be added to the function output.
  * This method is called by Fit algorithm after calling function(double*out)
110
  * @param out :: The output form function(double* out) to which the penalty will be added.
111
112
113
114
  */
void IFitFunction::addPenalty(double *out)const
{
    double penalty = 0.;
115
    for(size_t i=0;i<nParams();++i)
116
117
118
119
120
121
122
123
    {
      API::IConstraint* c = getConstraint(i);
      if (c)
      {
        penalty += c->check();
      }
    }

124
    size_t n = dataSize() - 1;
125
126
127
128
129
130
    // add penalty to first and last point and every 10th point in between
    if ( penalty != 0.0 )
    {
      out[0] += penalty;
      out[n] += penalty;

131
      for (size_t i = 9; i < n; i+=10)
132
133
134
135
136
137
138
139
140
      {
        out[i] += penalty;
      }
    }
}

/**
  * If a penalty was added to the function output the derivatives are modified accordingly.
  * This method is called by Fit algorithm after calling functionDeriv(Jacobian *out)
141
  * @param out :: The Jacobian to be corrected
142
143
144
  */
void IFitFunction::addPenaltyDeriv(Jacobian *out)const
{
145
146
  size_t n = dataSize() - 1;
  for(size_t i=0;i<nParams();++i)
147
148
149
150
151
152
  {  
    API::IConstraint* c = getConstraint(i);
    if (c)
    {
      double penalty = c->checkDeriv();
      if ( penalty != 0.0 )
153
      {
154
        size_t ia = activeIndex(i);
155
156
157
158
159
        double deriv = out->get(0,ia);
        out->set(0,ia,deriv + penalty);
        deriv = out->get(n,ia);
        out->set(n,ia,deriv+penalty);

160
        for (size_t j = 9; j < n; j+=10)
161
162
163
164
        {
          deriv = out->get(j,ia);
          out->set(j,ia,deriv+penalty);
        }
165
      }
166
167
    } // if (c)
  }
168
169
170

}

171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
///**
// * Writes a string that can be used in Fit.IFitFunction to create a copy of this IFitFunction
// * @return string representation of the function
// */
//std::string IFitFunction::asString()const
//{
//  std::ostringstream ostr;
//  ostr << "name="<<this->name();
//  std::vector<std::string> attr = this->getAttributeNames();
//  for(size_t i=0;i<attr.size();i++)
//  {
//    std::string attName = attr[i];
//    std::string attValue = this->getAttribute(attr[i]).value();
//    if (!attValue.empty())
//    {
//      ostr<<','<<attName<<'='<<attValue;
//    }
//  }
//  for(size_t i=0;i<nParams();i++)
//  {
//    ostr<<','<<parameterName(i)<<'='<<getParameter(i);
//  }
//  std::string constraints;
//  for(size_t i=0;i<nParams();i++)
//  {
//    const IConstraint* c = getConstraint(i);
//    if (c)
//    {
//      std::string tmp = c->asString();
//      if (!tmp.empty())
//      {
//        if (!constraints.empty())
//        {
//          constraints += ",";
//        }
//        constraints += tmp;
//      }
//    }
//  }
//  if (!constraints.empty())
//  {
//    ostr << ",constraints=(" << constraints << ")";
//  }
//
//  std::string ties;
//  for(size_t i=0;i<nParams();i++)
//  {
//    const ParameterTie* tie = getTie(i);
//    if (tie)
//    {
//      std::string tmp = tie->asString(this);
//      if (!tmp.empty())
//      {
//        if (!ties.empty())
//        {
//          ties += ",";
//        }
//        ties += tmp;
//      }
//    }
//  }
//  if (!ties.empty())
//  {
//    ostr << ",ties=(" << ties << ")";
//  }
//  return ostr.str();
//}
//
///** Set a function handler
// * @param handler :: A new handler
// */
//void IFitFunction::setHandler(FitFunctionHandler* handler)
//{
//  m_handler = handler;
//  if (handler && handler->function() != this)
//  {
//    throw std::runtime_error("Function handler points to a different function");
//  }
//  m_handler->init();
//}
//
252
/**
253
 * Operator \<\<
254
255
 * @param ostr :: The output stream
 * @param f :: The IFitFunction
256
257
258
259
260
261
262
 */
std::ostream& operator<<(std::ostream& ostr,const IFitFunction& f)
{
  ostr << f.asString();
  return ostr;
}

263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
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
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
///**
// * Const attribute visitor returning the type of the attribute
// */
//class AttType: public IFitFunction::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";}
//};
//
//std::string IFitFunction::Attribute::type()const
//{
//  AttType tmp;
//  return apply(tmp);
//}
//
///**
// * Const attribute visitor returning the type of the attribute
// */
//class AttValue: public IFitFunction::ConstAttributeVisitor<std::string>
//{
//public:
//  AttValue(bool quoteString=false) : 
//    IFitFunction::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);}
//
//private:
//  /// Flag to quote a string value returned
//  bool m_quoteString;
//};
//
//std::string IFitFunction::Attribute::value()const
//{
//  AttValue tmp(m_quoteValue);
//  return apply(tmp);
//}
//
//std::string IFitFunction::Attribute::asString()const
//{
//  if( m_quoteValue ) return asQuotedString();
//  
//  try
//  {
//    return boost::get<std::string>(m_data);
//  }
//  catch(...)
//  {
//    throw std::runtime_error("Trying to access a "+type()+" attribute "
//      "as string");
//  }
//}
//
//std::string IFitFunction::Attribute::asQuotedString()const
//{
//  std::string attr;
//
//  try
//  {
//    attr = boost::get<std::string>(m_data);
//  }
//  catch(...)
//  {
//    throw std::runtime_error("Trying to access a "+type()+" attribute "
//      "as string");
//  }
//  std::string quoted(attr);
//  if( *(attr.begin()) != '\"' ) quoted = "\"" + attr;
//  if( *(quoted.end() - 1) != '\"' ) quoted += "\"";
//
//  return quoted;
//}
//
//std::string IFitFunction::Attribute::asUnquotedString()const
//{
//  std::string attr;
//
//  try
//  {
//    attr = boost::get<std::string>(m_data);
//  }
//  catch(...)
//  {
//    throw std::runtime_error("Trying to access a "+type()+" attribute "
//      "as string");
//  }
//  std::string unquoted(attr);
//  if( *(attr.begin()) == '\"' ) unquoted = std::string(attr.begin() + 1, attr.end() - 1);
//  if( *(unquoted.end() - 1) == '\"' ) unquoted = std::string(unquoted.begin(), unquoted.end() - 1);
//  
//  return unquoted;
//}
//
//int IFitFunction::Attribute::asInt()const
//{
//  try
//  {
//    return boost::get<int>(m_data);
//  }
//  catch(...)
//  {
//    throw std::runtime_error("Trying to access a "+type()+" attribute "
//      "as int");
//  }
//}
//
//double IFitFunction::Attribute::asDouble()const
//{
//  try
//  {
//    return boost::get<double>(m_data);
//  }
//  catch(...)
//  {
//    throw std::runtime_error("Trying to access a "+type()+" attribute "
//      "as double");
//  }
//}
//
///** Sets new value if attribute is a string. If the type is wrong 
// * throws an exception
// * @param str :: The new value
// */
//void IFitFunction::Attribute::setString(const std::string& str)
//{
//  try
//  {
//    boost::get<std::string>(m_data) = str;
//  }
//  catch(...)
//  {
//    throw std::runtime_error("Trying to access a "+type()+" attribute "
//      "as string");
//  }
//}
//
///** Sets new value if attribute is a double. If the type is wrong 
// * throws an exception
// * @param d :: The new value
// */
//void IFitFunction::Attribute::setDouble(const double& d)
//{
//  try
//  {
//    boost::get<double>(m_data) = d;
//  }
//  catch(...)
//  {
//    throw std::runtime_error("Trying to access a "+type()+" attribute "
//      "as double");
//  }
//}
//
///** Sets new value if attribute is an int. If the type is wrong 
// * throws an exception
// * @param i :: The new value
// */
//void IFitFunction::Attribute::setInt(const int& i)
//{
//  try
//  {
//    boost::get<int>(m_data) = i;
//  }
//  catch(...)
//  {
//    throw std::runtime_error("Trying to access a "+type()+" attribute "
//      "as int");
//  }
//}
//
///**
// * Attribute visitor setting new value to an attribute
// */
//class SetValue: public IFitFunction::AttributeVisitor<>
//{
//public:
//  /**
//   * Constructor
//   * @param value :: The value to set
//   */
//  SetValue(const std::string& value):m_value(value){}
//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);
//  }
//private:
//  std::string m_value; ///<the value as a string
//};
//
///** Set value from a string. Throws exception if the string has wrong format
// * @param str :: String representation of the new value
// */
//void IFitFunction::Attribute::fromString(const std::string& str)
//{
//  SetValue tmp(str);
//  apply(tmp);
//}
491
492
493

} // namespace API
} // namespace Mantid