Skip to content

Commit

Permalink
Refs #4394 Some refactoring
Browse files Browse the repository at this point in the history
  • Loading branch information
@jmborr
jmborr committed Dec 20, 2013
1 parent 4d5d49c commit 90cd767
Show file tree
Hide file tree
Showing 2 changed files with 103 additions and 96 deletions.
20 changes: 10 additions & 10 deletions Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/DiffSphere.h
View file
Expand Up @@ -67,8 +67,8 @@ class DLLExport ElasticDiffSphere : public DeltaFunction
/// structure to hold info on Volino's coefficients
struct xnlc {
double x;
unsigned int l;
unsigned int n;
size_t n;
size_t l;
};

/// simple structure to hold a linear interpolation of factor J around its numerical divergence point
Expand All @@ -85,8 +85,11 @@ class DLLExport InelasticDiffSphere : public API::ParamFunction, public API::IFu
public:

InelasticDiffSphere();

virtual ~InelasticDiffSphere() {}

virtual void init();

virtual std::string name()const{return "InelasticDiffSphere"; }

virtual const std::string category() const { return "QENS"; }
Expand All @@ -104,29 +107,26 @@ class DLLExport InelasticDiffSphere : public API::ParamFunction, public API::IFu
/// initialize the Xnl coefficients
void initXnlCoeff();

/// initialize the alpha coefficients
/// initialize the m_alpha coefficients
void initAlphaCoeff();

/// initialize the list of Linearized J values
void initLinJlist();

/// xnl coefficients
std::vector< xnlc > xnl;
std::vector< xnlc > m_xnl;

/// certain coefficients invariant during fitting
std::vector< double > alpha;
std::vector< double > m_alpha;

/// maximum value of l in xnlist
unsigned int lmax;

/// number of coefficients
unsigned int ncoeff;
size_t lmax;

/// linear interpolation zone around the numerical divergence of factor J
double m_divZone;

/// list of linearized J values
std::vector< linearJ > linearJlist;
std::vector< linearJ > m_linearJlist;

}; // end of class InelasticDiffSphere

Expand Down
179 changes: 93 additions & 86 deletions Code/Mantid/Framework/CurveFitting/src/DiffSphere.cpp
View file
Expand Up @@ -48,74 +48,85 @@ double ElasticDiffSphere::HeightPrefactor() const

DECLARE_FUNCTION(InelasticDiffSphere);

// initialize class attribute xnl with a list of coefficients in string format
// initialize class attribute m_xnl with a list of coefficients in string format
void InelasticDiffSphere::initXnlCoeff(){
/* List of 98 coefficients sorted by increasing value (F.Volino,Mol. Phys. 41,271-279,1980)
* For each coefficient, the triad (coeff, l, n) is defined
*/
ncoeff = 98;
double xnlist[]={
2.081576, 1, 0, 3.342094, 2, 0, 4.493409, 0, 1, 4.514100, 3, 0,
5.646704, 4, 0, 5.940370, 1, 1, 6.756456, 5, 0, 7.289932, 2, 1,
7.725252, 0, 2, 7.851078, 6, 0, 8.583755, 3, 1, 8.934839, 7, 0,
9.205840, 1, 2, 9.840446, 4, 1, 10.010371, 8, 0, 10.613855, 2, 2,
10.904122, 0, 3, 11.070207, 5, 1, 11.079418, 9, 0, 11.972730, 3, 2,
12.143204, 10, 0, 12.279334, 6, 1, 12.404445, 1, 3, 13.202620, 11, 0,
13.295564, 4, 2, 13.472030, 7, 1, 13.846112, 2, 3, 14.066194, 0, 4,
14.258341, 12, 0, 14.590552, 5, 2, 14.651263, 8, 1, 15.244514, 3, 3,
15.310887, 13, 0, 15.579236, 1, 4, 15.819216, 9, 1, 15.863222, 6, 2,
16.360674, 14, 0, 16.609346, 4, 3, 16.977550, 10, 1, 17.042902, 2, 4,
17.117506, 7, 2, 17.220755, 0, 5, 17.408034, 15, 0, 17.947180, 5, 3,
18.127564, 11, 1, 18.356318, 8, 2, 18.453241, 16, 0, 18.468148, 3, 4,
18.742646, 1, 5, 19.262710, 6, 3, 19.270294, 12, 1, 19.496524, 17, 0,
19.581889, 9, 2, 19.862424, 4, 4, 20.221857, 2, 5, 20.371303, 0, 6,
20.406581, 13, 1, 20.538074, 18, 0, 20.559428, 7, 3, 20.795967, 10, 2,
21.231068, 5, 4, 21.537120, 14, 1, 21.578053, 19, 0, 21.666607, 3, 5,
21.840012, 8, 3, 21.899697, 1, 6, 21.999955, 11, 2, 22.578058, 6, 4,
22.616601, 20, 0, 22.662493, 15, 1, 23.082796, 4, 5, 23.106568, 9, 3,
23.194996, 12, 2, 23.390490, 2, 6, 23.519453, 0, 7, 23.653839, 21, 0,
23.783192, 16, 1, 23.906450, 7, 4, 24.360789, 10, 3, 24.382038, 13, 2,
24.474825, 5, 5, 24.689873, 22, 0, 24.850085, 3, 6, 24.899636, 17, 1,
25.052825, 1, 7, 25.218652, 8, 4, 25.561873, 14, 2, 25.604057, 11, 3,
25.724794, 23, 0, 25.846084, 6, 5, 26.012188, 18, 1, 26.283265, 4, 6,
26.516603, 9, 4, 26.552589, 2, 7, 26.666054, 0, 8, 26.735177, 15, 2,
26.758685, 24, 0, 26.837518, 12, 3
size_t ncoeff = 98;

double xvalues[] =
{
2.081576, 3.342094, 4.493409, 4.514100, 5.646704, 5.940370, 6.756456, 7.289932,
7.725252, 7.851078, 8.583755, 8.934839, 9.205840, 9.840446, 10.010371, 10.613855,
10.904122, 11.070207, 11.079418, 11.972730, 12.143204, 12.279334, 12.404445,
13.202620, 13.295564, 13.472030, 13.846112, 14.066194, 14.258341, 14.590552,
14.651263, 15.244514, 15.310887, 15.579236, 15.819216, 15.863222, 16.360674,
16.609346, 16.977550, 17.042902, 17.117506, 17.220755, 17.408034, 17.947180,
18.127564, 18.356318, 18.453241, 18.468148, 18.742646, 19.262710, 19.270294,
19.496524, 19.581889, 19.862424, 20.221857, 20.371303, 20.406581, 20.538074,
20.559428, 20.795967, 21.231068, 21.537120, 21.578053, 21.666607, 21.840012,
21.899697, 21.999955, 22.578058, 22.616601, 22.662493, 23.082796, 23.106568,
23.194996, 23.390490, 23.519453, 23.653839, 23.783192, 23.906450, 24.360789,
24.382038, 24.474825, 24.689873, 24.850085, 24.899636, 25.052825, 25.218652,
25.561873, 25.604057, 25.724794, 25.846084, 26.012188, 26.283265, 26.516603,
26.552589, 26.666054, 26.735177, 26.758685, 26.837518
};

size_t lvalues[] =
{
1, 2, 0, 3, 4, 1, 5, 2, 0, 6, 3, 7, 1, 4, 8, 2, 0, 5, 9, 3, 10, 6, 1, 11, 4, 7,
2, 0, 12, 5, 8, 3, 13, 1, 9, 6, 14, 4, 10, 2, 7, 0, 15, 5, 11, 8, 16, 3, 1, 6, 12,
17, 9, 4, 2, 0, 13, 18, 7, 10, 5, 14, 19, 3, 8, 1, 11, 6, 20, 15, 4, 9, 12, 2, 0,
21, 16, 7, 10, 13, 5, 22, 3, 17, 1, 8, 14, 11, 23, 6, 18, 4, 9, 2, 0, 15, 24, 12
};

for( unsigned int i = 0; i < ncoeff; i += 3 )
size_t nvalues[] =
{
0, 0, 1, 0, 0, 1, 0, 1, 2, 0, 1, 0, 2, 1, 0, 2, 3, 1, 0, 2, 0, 1, 3, 0, 2, 1, 3,
4, 0, 2, 1, 3, 0, 4, 1, 2, 0, 3, 1, 4, 2, 5, 0, 3, 1, 2, 0, 4, 5, 3, 1, 0, 2, 4, 5,
6, 1, 0, 3, 2, 4, 1, 0, 5, 3, 6, 2, 4, 0, 1, 5, 3, 2, 6, 7, 0, 1, 4, 3, 2, 5, 0,
6, 1, 7, 4, 2, 3, 0, 5, 1, 6, 4, 7, 8, 2, 0, 3
};

for( size_t i = 0; i < ncoeff; i += 3 )
{
xnlc coeff;
coeff.x = xnlist[ i ]; //value of the coefficient
coeff.l = ( unsigned int )( xnlist[ i + 1 ] ); //corresponding n
coeff.n = ( unsigned int )( xnlist [ i + 2 ] ); //corresponding l
xnl.push_back( coeff );
coeff.x = xvalues[ i ];
coeff.l = lvalues[ i ];
coeff.n = nvalues[ i ];
m_xnl.push_back( coeff );
}
}

//initialize a set of coefficients that will remain constant during fitting
void InelasticDiffSphere::initAlphaCoeff(){
for( std::vector< xnlc >::const_iterator it = xnl.begin(); it != xnl.end(); ++it )
for( std::vector< xnlc >::const_iterator it = m_xnl.begin(); it != m_xnl.end(); ++it )
{
double x = it->x;
double x = it->x; // eigenvalue for a (n, l) pair
double l = ( double )( it->l );
alpha.push_back( ( 2.0 * l + 1 ) * ( 6.0 * x*x / ( x*x - l*(l + 1) ) ) / M_PI );
m_alpha.push_back( ( 2.0 * l + 1 ) * ( 6.0 * x*x / ( x*x - l*(l + 1) ) ) / M_PI );
}
}

//initialize linear interpolation of factor J around its numerical divergence point a = it->x
void InelasticDiffSphere::initLinJlist(){
for( std::vector< xnlc >::const_iterator it = xnl.begin(); it != xnl.end(); ++it )
/* Factor "J" is defined as [Q*a*j(l+1,Q*a) - l*j(l,Q*a)] / [(Q*a)^2 - x^2]
* Both numerator and denominator goes to zero when Q*a approaches x, giving rise to
* numerical indeterminacies. To avoid them, we will interpolate linearly.
*/
void InelasticDiffSphere::initLinJlist()
{
for( size_t i = 0; i < m_xnl.size(); i++ )
{
linearJ abJ;
double x = it->x;
unsigned int l = it->l;
double a = x - m_divZone; //left of the numerical divergence point
double J0 = ( a * boost::math::sph_bessel( l + 1, a ) - l * boost::math::sph_bessel( l, a ) ) / ( a*a - x*x );
a = x + m_divZone; //right of the numerical divergence point
double J1 = ( a * boost::math::sph_bessel( l + 1, a ) - l * boost::math::sph_bessel( l, a ) ) / ( a*a - x*x );
double x = m_xnl[ i ].x; // eigenvalue for a (n, l) pair
unsigned int l = ( unsigned int )( m_xnl[ i ].l );
double Qa = x - m_divZone; //left of the numerical divergence point
double J0 = ( Qa * boost::math::sph_bessel( l + 1, Qa ) - l * boost::math::sph_bessel( l, Qa ) ) / ( Qa*Qa - x*x );
Qa = x + m_divZone; //right of the numerical divergence point
double J1 = ( Qa * boost::math::sph_bessel( l + 1, Qa ) - l * boost::math::sph_bessel( l, Qa ) ) / ( Qa*Qa - x*x );
abJ.slope = ( J1 - J0 ) / ( 2 * m_divZone ); //slope of the linear interpolation
abJ.intercept = J0 - abJ.slope * ( x - m_divZone ); //intercept of the linear interpolation
linearJlist.push_back( abJ ); //store the parameters of the linear interpolation for this it->x
m_linearJlist.push_back( abJ ); //store the parameters of the linear interpolation for this it->x
}
}

Expand All @@ -126,7 +137,10 @@ InelasticDiffSphere::InelasticDiffSphere() : lmax( 24 ), m_divZone( 0.1 )
declareParameter( "Diffusion", 1.0, "Diffusion coefficient, in units of" );

declareAttribute( "Q", API::IFunction::Attribute( 1.0 ) );
}

void InelasticDiffSphere::init()
{
// Ensure positive values for Intensity, Radius, and Diffusion coefficient
BoundaryConstraint* IntensityConstraint = new BoundaryConstraint( this, "Intensity", 0, true );
addConstraint(IntensityConstraint);
Expand All @@ -137,52 +151,43 @@ InelasticDiffSphere::InelasticDiffSphere() : lmax( 24 ), m_divZone( 0.1 )
BoundaryConstraint* DiffusionConstraint = new BoundaryConstraint( this, "Diffusion", 0, true );
addConstraint( DiffusionConstraint );

initXnlCoeff(); // initialize this->xnl with the list of coefficients xnlist
initAlphaCoeff(); // initialize this->alpha, certain factors constant over the fit
initLinJlist(); // initialize this->linJlist, linear interpolation around numerical divergence
initXnlCoeff(); // initialize m_xnl with the list of coefficients xnlist
initAlphaCoeff(); // initialize m_alpha, certain factors constant over the fit
initLinJlist(); // initialize m_linearJlist, linear interpolation around numerical divergence
}

//calculate the coefficients for each Lorentzian
std::vector< double > InelasticDiffSphere::LorentzianCoefficients( double a ) const
{
//precompute the 2+lmax spherical bessel functions (26 in total)
double* jl = new double[ 2 + lmax ];
for( unsigned int l = 0; l <= 1 + lmax; l++ )
std::vector< double > jl( 2 + lmax );
for( size_t l = 0; l < 2 + lmax; l++ )
{
jl[ l ] = boost::math::sph_bessel( l, a);
jl[ l ] = boost::math::sph_bessel( ( unsigned int )( l ), a );
}

//store the coefficient of each Lorentzian in vector YJ(a,w)
size_t ncoeff = m_xnl.size();
std::vector< double > YJ( ncoeff );
std::vector< linearJ >::const_iterator itlinJ = linearJlist.begin();
//loop over all coefficients
std::vector< double >::const_iterator italpha = alpha.begin();
std::vector< double >::iterator itYJ = YJ.begin();
for( std::vector< xnlc >::const_iterator it = xnl.begin(); it != xnl.end(); ++it )

for( size_t i = 0; i < ncoeff; i++ )
{
//only to make expressions more readable
double x = it->x;
unsigned int l = it->l;
double x = m_xnl[ i ].x;
unsigned int l = ( unsigned int )( m_xnl[ i ].l );
//compute factors Y and J
double Y = *italpha; //Y is independent of parameters a and w, and independent of data x
/* J is dependent on parameter a, cannot be computed when active parameter a obeys a*a=c.
* Thus for each it->x we stored J(it->x-m_divZone) and J(it->x_m_divZone), and use linear
* interpolation
*/
double Y = m_alpha[ i ]; //Y is independent of parameters a and w, and independent of data x
double J;
if( fabs( a - x ) > m_divZone )
{
J = ( a * jl[l + 1] - l * jl[l] ) / ( a*a - x*x );
}else{
J = itlinJ->slope * a + itlinJ->intercept; //linear interpolation
}
*itYJ = Y * J * J;
++itYJ;
++italpha;
++itlinJ; //retrieve next linear interpolation
} // end of for(std::vector<xnlc>::const_iterator it=xnl.begin()
else
{
J = m_linearJlist[ i ].slope * a + m_linearJlist[ i ].intercept; // linear interpolation instead
}
YJ[ i ] = Y * J * J;
}

delete[] jl;
return YJ;
} // end of LorentzianCoefficients

Expand All @@ -196,27 +201,25 @@ void InelasticDiffSphere::function1D( double* out, const double* xValues, const

std::vector<double> YJ;
YJ = LorentzianCoefficients( Q * R );
for (unsigned int i = 0; i < nData; i++)
for (size_t i = 0; i < nData; i++)
{
double x = xValues[ i ];
//loop over all coefficients
out[ i ] = 0.0;
std::vector< double >::const_iterator itYJ = YJ.begin();
for( std::vector< xnlc >::const_iterator it = xnl.begin(); it != xnl.end(); ++it )
size_t ncoeff = m_xnl.size();
for ( size_t n = 0; n < ncoeff; n++ )
{
double z = it->x;
double zw = z * z * D / ( R * R );
double z = m_xnl[ n ].x; // eigenvalue
double zw = z * z * D / ( R * R ); // HWHM
double L = zw / ( zw * zw + x * x ); //Lorentzian
out[i] += I * ( *itYJ ) * L;
++itYJ; //retrieve next coefficient
} // end of for(std::vector<xnlc>::const_iterator it....
} // end of for (unsigned int i...
} // end of void DiffSphere::functionMW
out[i] += I * YJ[ n ] * L;
}
}
}

/** Calculate numerical derivatives.
* @param domain :: The domain of the function
* @param jacobian :: A Jacobian matrix. It is expected to have dimensions of domain.size() by nParams().
*/
void InelasticDiffSphere::calNumericalDeriv2( const API::FunctionDomain& domain, API::Jacobian& jacobian )
{
const double minDouble = std::numeric_limits< double >::min();
Expand Down Expand Up @@ -267,21 +270,25 @@ void InelasticDiffSphere::calNumericalDeriv2( const API::FunctionDomain& domain,
}
}
} // calNumericalDeriv()
*/

/*
/// Using numerical derivative
void InelasticDiffSphere::functionDeriv( const API::FunctionDomain& domain, API::Jacobian& jacobian )
{
this->calNumericalDeriv( domain, jacobian );
return;
}
*/

/*
/// Using numerical derivative
void InelasticDiffSphere::functionDeriv1D( Jacobian* jacobian, const double* xValues, const size_t nData )
{
FunctionDomain1DView domain( xValues, nData );
this->calNumericalDeriv( domain, *jacobian );
}

*/

/* Propagate the attribute to its member functions.
* NOTE: we pass this->getAttribute(name) by reference, thus the same
Expand Down

0 comments on commit 90cd767

Please sign in to comment.