Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
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
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
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidAPI/IPeakFunction.h"
#include "MantidAPI/Jacobian.h"
#include "MantidKernel/Exception.h"
#include "MantidKernel/ConfigService.h"
#include <boost/lexical_cast.hpp>
#include <cmath>
namespace Mantid
{
namespace API
{
/** A Jacobian for individual functions
*/
class PartialJacobian1: public Jacobian
{
Jacobian* m_J; ///< pointer to the overall Jacobian
int m_iY0; ///< offset in the overall Jacobian for a particular function
public:
/** Constructor
* @param J :: A pointer to the overall Jacobian
* @param iY0 :: The data offset for a particular function
*/
PartialJacobian1(Jacobian* J,int iY0):m_J(J),m_iY0(iY0)
{}
/**
* Overridden Jacobian::set(...).
* @param iY :: The index of the data point
* @param iP :: The parameter index of an individual function.
* @param value :: The derivative value
*/
void set(size_t iY, size_t iP, double value)
{
m_J->set(m_iY0 + iY,iP,value);
}
/**
* Overridden Jacobian::get(...).
* @param iY :: The index of the data point
* @param iP :: The parameter index of an individual function.
*/
double get(size_t iY, size_t iP)
{
return m_J->get(m_iY0 + iY,iP);
}
};
/// Default value for the peak radius
int IPeakFunction::s_peakRadius = 5;
/**
* Constructor. Sets peak radius to the value of curvefitting.peakRadius property
*/
IPeakFunction::IPeakFunction()
{
int peakRadius;
if ( Kernel::ConfigService::Instance().getValue("curvefitting.peakRadius",peakRadius) )
{
if ( peakRadius != s_peakRadius )
{
setPeakRadius(peakRadius);
}
}
}
/**
* General implementation of the method for all peaks. Limits the peak evaluation to
* a certain number of FWHMs around the peak centre. The outside points are set to 0.
* Calls functionLocal() to compute the actual values
* @param out :: Output function values
* @param xValues :: X values for data points
* @param nData :: Number of data points
*/
void IPeakFunction::function1D(double* out, const double* xValues, const size_t nData)const
{
double c = this->centre();
double dx = fabs(s_peakRadius*this->fwhm());
int i0 = -1;
int n = 0;
for(size_t i = 0; i < nData; ++i)
{
if (fabs(xValues[i] - c) < dx)
{
if (i0 < 0) i0 = static_cast<int>(i);
++n;
}
else
{
out[i] = 0.0;
}
}
if (i0 < 0 || n == 0) return;
this->functionLocal(out+i0, xValues+i0, n);
}
/**
* General implementation of the method for all peaks. Calculates derivatives only
* for a range of x values limited to a certain number of FWHM around the peak centre.
* For the points outside the range all derivatives are set to 0.
* Calls functionDerivLocal() to compute the actual values
* @param out :: Derivatives
* @param xValues :: X values for data points
* @param nData :: Number of data points
*/
void IPeakFunction::functionDeriv1D(Jacobian* out, const double* xValues, const size_t nData)
{
double c = this->centre();
double dx = fabs(s_peakRadius*this->fwhm());
int i0 = -1;
int n = 0;
for(size_t i = 0; i < nData; ++i)
{
if (fabs(xValues[i] - c) < dx)
{
if (i0 < 0) i0 = static_cast<int>(i);
++n;
}
else
{
for(size_t ip = 0; ip < this->nParams(); ++ip)
{
out->set(i,ip, 0.0);
}
}
}
if (i0 < 0 || n == 0) return;
PartialJacobian1 J(out,i0);
this->functionDerivLocal(&J,xValues+i0,n);
}
void IPeakFunction::setPeakRadius(const int& r)
{
if (r > 0)
{
s_peakRadius = r;
std::string setting = boost::lexical_cast<std::string>(r);
Kernel::ConfigService::Instance().setString("curvefitting.peakRadius",setting);
}
}
} // namespace API
} // namespace Mantid