Skip to content

Commit

Permalink
Merge branch 'master' into 9198_add_diffsphere_support_to_convfit
Browse files Browse the repository at this point in the history 
Refs #9198
  • Loading branch information
@DanNixon
DanNixon committed Mar 25, 2015
2 parents ead902d + ca7daa3 commit 62af761
Show file tree
Hide file tree
Showing 26 changed files with 410 additions and 147 deletions.
5 changes: 4 additions & 1 deletion Code/Mantid/Framework/CurveFitting/src/DiffRotDiscreteCircle.cpp
View file
Expand Up @@ -54,6 +54,7 @@ InelasticDiffRotDiscreteCircle::InelasticDiffRotDiscreteCircle()
declareParameter("Decay", 1.0, "Inverse of transition rate, in nanoseconds "
"if energy in micro-ev, or picoseconds if "
"energy in mili-eV");
declareParameter("Shift", 0.0, "Shift in domain");

declareAttribute("Q", API::IFunction::Attribute(0.5));
declareAttribute("N", API::IFunction::Attribute(3));
Expand Down Expand Up @@ -82,6 +83,7 @@ void InelasticDiffRotDiscreteCircle::function1D(double *out,
const double rate = m_hbar / getParameter("Decay"); // micro-eV or mili-eV
const double Q = getAttribute("Q").asDouble();
const int N = getAttribute("N").asInt();
const double S = getParameter("Shift");

std::vector<double> sph(N);
for (int k = 1; k < N; k++) {
Expand All @@ -97,7 +99,7 @@ void InelasticDiffRotDiscreteCircle::function1D(double *out,
}

for (size_t i = 0; i < nData; i++) {
double w = xValues[i];
double w = xValues[i] - S;
double S = 0.0;
for (int l = 1; l < N; l++) // l goes up to N-1
{
Expand Down Expand Up @@ -166,6 +168,7 @@ void DiffRotDiscreteCircle::init() {
setAlias("f1.Intensity", "Intensity");
setAlias("f1.Radius", "Radius");
setAlias("f1.Decay", "Decay");
setAlias("f1.Shift", "Shift");

// Set the ties between Elastic and Inelastic parameters
addDefaultTies("f0.Height=f1.Intensity,f0.Radius=f1.Radius");
Expand Down
5 changes: 4 additions & 1 deletion Code/Mantid/Framework/CurveFitting/src/DiffSphere.cpp
View file
Expand Up @@ -145,6 +145,7 @@ InelasticDiffSphere::InelasticDiffSphere()
declareParameter("Diffusion", 0.05, "Diffusion coefficient, in units of "
"A^2*THz, if energy in meV, or A^2*PHz "
"if energy in ueV");
declareParameter("Shift", 0.0, "Shift in domain");

declareAttribute("Q", API::IFunction::Attribute(1.0));
}
Expand Down Expand Up @@ -194,6 +195,7 @@ void InelasticDiffSphere::function1D(double *out, const double *xValues,
const double R = getParameter("Radius");
const double D = getParameter("Diffusion");
const double Q = getAttribute("Q").asDouble();
const double S = getParameter("Shift");

// // Penalize negative parameters
if (I < std::numeric_limits<double>::epsilon() ||
Expand All @@ -216,7 +218,7 @@ void InelasticDiffSphere::function1D(double *out, const double *xValues,
std::vector<double> YJ;
YJ = LorentzianCoefficients(Q * R); // The (2l+1)*A_{n,l}
for (size_t i = 0; i < nData; i++) {
double energy = xValues[i]; // from meV to THz (or from micro-eV to PHz)
double energy = xValues[i] - S; // from meV to THz (or from micro-eV to PHz)
out[i] = 0.0;
for (size_t n = 0; n < ncoeff; n++) {
double L = (1.0 / M_PI) * HWHM[n] /
Expand Down Expand Up @@ -272,6 +274,7 @@ void DiffSphere::init() {
setAlias("f1.Intensity", "Intensity");
setAlias("f1.Radius", "Radius");
setAlias("f1.Diffusion", "Diffusion");
setAlias("f1.Shift", "Shift");

// Set the ties between Elastic and Inelastic parameters
addDefaultTies("f0.Height=f1.Intensity,f0.Radius=f1.Radius");
Expand Down
127 changes: 80 additions & 47 deletions Code/Mantid/Framework/CurveFitting/test/DiffRotDiscreteCircleTest.h
View file
Expand Up @@ -25,6 +25,7 @@ class DiffRotDiscreteCircleTest : public CxxTest::TestSuite
static DiffRotDiscreteCircleTest *createSuite() { return new DiffRotDiscreteCircleTest(); }
static void destroySuite( DiffRotDiscreteCircleTest *suite ) { delete suite; }


// convolve the elastic part with a resolution function, here a Gaussian
void testDiffRotDiscreteCircleElastic()
{
Expand Down Expand Up @@ -75,58 +76,17 @@ class DiffRotDiscreteCircleTest : public CxxTest::TestSuite
} // testDiffRotDiscreteCircleElastic


/// Fit the convolution of the inelastic part with a Gaussian resolution function
void testDiffRotDiscreteCircleInelastic()
{
/* Note: it turns out that parameters Intensity and Radius are highly covariant, so that more than one minimum exists.
* Thus, I tied parameter Radius. This is OK since one usually knows the radius of the circle of the jumping diffusion
*/

// initialize the fitting function in a Fit algorithm
// Parameter units are assumed in micro-eV, Angstroms, Angstroms**(-1), and nano-seconds. Intensities have arbitrary units
std::string funtion_string = "(composite=Convolution,FixResolution=true,NumDeriv=true;name=Gaussian,Height=1.0,PeakCentre=0.0,Sigma=20.0,ties=(Height=1.0,PeakCentre=0.0,Sigma=20.0);name=InelasticDiffRotDiscreteCircle,N=3,Q=0.5,Intensity=47.014,Radius=1.567,Decay=7.567)";

// Initialize the fit function in the Fit algorithm
Mantid::CurveFitting::Fit fitalg;
TS_ASSERT_THROWS_NOTHING( fitalg.initialize() );
TS_ASSERT( fitalg.isInitialized() );
fitalg.setProperty( "Function", funtion_string );

// create the data workspace by evaluating the fit function in the Fit algorithm
auto data_workspace = generateWorkspaceFromFitAlgorithm( fitalg );
//saveWorkspace( data_workspace, "/tmp/junk.nxs" ); // for debugging purposes only

//override the function with new parameters, then do the Fit
funtion_string = "(composite=Convolution,FixResolution=true,NumDeriv=true;name=Gaussian,Height=1.0,PeakCentre=0.0,Sigma=20.0,ties=(Height=1.0,PeakCentre=0.0,Sigma=20.0);name=InelasticDiffRotDiscreteCircle,N=3,Q=0.5,Intensity=10.0,Radius=1.567,Decay=20.0,ties=(Radius=1.567))";
fitalg.setProperty( "Function", funtion_string );
fitalg.setProperty( "InputWorkspace", data_workspace );
fitalg.setPropertyValue( "WorkspaceIndex", "0" );
TS_ASSERT_THROWS_NOTHING( TS_ASSERT( fitalg.execute() ) );
TS_ASSERT( fitalg.isExecuted() );
runDiffRotDiscreteCircleInelasticTest(0.0);
}

// check Chi-square is small
const double chi_squared = fitalg.getProperty("OutputChi2overDoF");
TS_ASSERT_LESS_THAN( chi_squared, 0.001 );
//std::cout << "\nchi_squared = " << chi_squared << "\n"; // only for debugging purposes

// check the parameters of the resolution did not change
Mantid::API::IFunction_sptr fitalg_function = fitalg.getProperty( "Function" );
auto fitalg_conv = boost::dynamic_pointer_cast<Mantid::CurveFitting::Convolution>( fitalg_function ) ;
Mantid::API::IFunction_sptr fitalg_resolution = fitalg_conv->getFunction( 0 );
TS_ASSERT_DELTA( fitalg_resolution -> getParameter( "PeakCentre" ), 0.0, 0.00001 ); // allow for a small percent variation
TS_ASSERT_DELTA( fitalg_resolution -> getParameter( "Height" ), 1.0, 1.0 * 0.001 ); // allow for a small percent variation
TS_ASSERT_DELTA( fitalg_resolution -> getParameter( "Sigma" ), 20.0, 20.0* 0.001 ); // allow for a small percent variation
//std::cout << "\nPeakCentre = " << fitalg_resolution->getParameter("PeakCentre") << " Height= " << fitalg_resolution->getParameter("Height") << " Sigma=" << fitalg_resolution->getParameter("Sigma") << "\n"; // only for debugging purposes

// check the parameters of the inelastic part
Mantid::API::IFunction_sptr fitalg_structure_factor = fitalg_conv->getFunction( 1 );
TS_ASSERT_DELTA( fitalg_structure_factor -> getParameter( "Intensity" ), 47.014, 47.014 * 0.05 ); // allow for a small percent variation
TS_ASSERT_DELTA( fitalg_structure_factor -> getParameter( "Radius" ), 1.567, 1.567 * 0.05 ); // allow for a small percent variation
TS_ASSERT_DELTA( fitalg_structure_factor -> getParameter( "Decay" ), 7.567, 7.567 * 0.05 ); // allow for a small percent variation
//std::cout << "\nGOAL: Intensity = 47.014, Radius = 1.567, Decay = 7.567\n"; // only for debugging purposes
//std::cout << "OPTIMIZED: Intensity = " << fitalg_structure_factor->getParameter("Intensity") << " Radius = " << fitalg_structure_factor->getParameter("Radius") << " Decay = " << fitalg_structure_factor->getParameter("Decay") << "\n"; // only for debugging purposes
void testDiffRotDiscreteCircleInelasticWithShift()
{
runDiffRotDiscreteCircleInelasticTest(0.5);
}

} // testDiffRotDiscreteCircleElastic

/* Check the particular case for N = 3
* In this case, the inelastic part should reduce to a single Lorentzian in 'w':
Expand Down Expand Up @@ -293,6 +253,79 @@ class DiffRotDiscreteCircleTest : public CxxTest::TestSuite


private:
/// Fit the convolution of the inelastic part with a Gaussian resolution function
void runDiffRotDiscreteCircleInelasticTest(const double S)
{
/* Note: it turns out that parameters Intensity and Radius are highly covariant, so that more than one minimum exists.
* Thus, I tied parameter Radius. This is OK since one usually knows the radius of the circle of the jumping diffusion
*/
const double I(47.014);
const double R(1.567);
const double tao(7.567);

// initialize the fitting function in a Fit algorithm
// Parameter units are assumed in micro-eV, Angstroms, Angstroms**(-1), and nano-seconds. Intensities have arbitrary units
std::ostringstream function_stream;
function_stream << "(composite=Convolution,FixResolution=true,NumDeriv=true;"
<< "name=Gaussian,Height=1.0,PeakCentre=0.0,Sigma=20.0,"
<< "ties=(Height=1.0,PeakCentre=0.0,Sigma=20.0);"
<< "name=InelasticDiffRotDiscreteCircle,N=3,Q=0.5,"
<< "Intensity=" << I
<< ",Radius=" << R
<< ",Decay=" << tao
<< ",Shift=" << S << ")";

// Initialize the fit function in the Fit algorithm
Mantid::CurveFitting::Fit fitalg;
TS_ASSERT_THROWS_NOTHING( fitalg.initialize() );
TS_ASSERT( fitalg.isInitialized() );
fitalg.setProperty( "Function", function_stream.str() );

function_stream.str( std::string() );
function_stream.clear();

// create the data workspace by evaluating the fit function in the Fit algorithm
auto data_workspace = generateWorkspaceFromFitAlgorithm( fitalg );
//saveWorkspace( data_workspace, "/tmp/junk.nxs" ); // for debugging purposes only

//override the function with new parameters, then do the Fit
function_stream << "(composite=Convolution,FixResolution=true,NumDeriv=true;"
<< "name=Gaussian,Height=1.0,PeakCentre=0.0,Sigma=20.0,"
<< "ties=(Height=1.0,PeakCentre=0.0,Sigma=20.0);"
<< "name=InelasticDiffRotDiscreteCircle,N=3,Q=0.5,"
<< "Intensity=10.0,Radius=1.567,Decay=20.0"
<< ",ties=(Radius=" << R << "))";
fitalg.setProperty( "Function", function_stream.str() );
fitalg.setProperty( "InputWorkspace", data_workspace );
fitalg.setPropertyValue( "WorkspaceIndex", "0" );
TS_ASSERT_THROWS_NOTHING( TS_ASSERT( fitalg.execute() ) );
TS_ASSERT( fitalg.isExecuted() );

// check Chi-square is small
const double chi_squared = fitalg.getProperty("OutputChi2overDoF");
TS_ASSERT_LESS_THAN( chi_squared, 0.001 );
//std::cout << "\nchi_squared = " << chi_squared << "\n"; // only for debugging purposes

// check the parameters of the resolution did not change
Mantid::API::IFunction_sptr fitalg_function = fitalg.getProperty( "Function" );
auto fitalg_conv = boost::dynamic_pointer_cast<Mantid::CurveFitting::Convolution>( fitalg_function ) ;
Mantid::API::IFunction_sptr fitalg_resolution = fitalg_conv->getFunction( 0 );
TS_ASSERT_DELTA( fitalg_resolution -> getParameter( "PeakCentre" ), 0.0, 0.00001 ); // allow for a small percent variation
TS_ASSERT_DELTA( fitalg_resolution -> getParameter( "Height" ), 1.0, 1.0 * 0.001 ); // allow for a small percent variation
TS_ASSERT_DELTA( fitalg_resolution -> getParameter( "Sigma" ), 20.0, 20.0* 0.001 ); // allow for a small percent variation
//std::cout << "\nPeakCentre = " << fitalg_resolution->getParameter("PeakCentre") << " Height= " << fitalg_resolution->getParameter("Height") << " Sigma=" << fitalg_resolution->getParameter("Sigma") << "\n"; // only for debugging purposes

// check the parameters of the inelastic part
Mantid::API::IFunction_sptr fitalg_structure_factor = fitalg_conv->getFunction( 1 );
TS_ASSERT_DELTA( fitalg_structure_factor -> getParameter( "Intensity" ), I, I * 0.05 ); // allow for a small percent variation
TS_ASSERT_DELTA( fitalg_structure_factor -> getParameter( "Radius" ), R, R * 0.05 ); // allow for a small percent variation
TS_ASSERT_DELTA( fitalg_structure_factor -> getParameter( "Decay" ), tao, tao * 0.05 ); // allow for a small percent variation
TS_ASSERT_DELTA( fitalg_structure_factor -> getParameter( "Shift" ), S, 0.00001 ); // allow for a small percent variation
//std::cout << "\nGOAL: Intensity = 47.014, Radius = 1.567, Decay = 7.567\n"; // only for debugging purposes
//std::cout << "OPTIMIZED: Intensity = " << fitalg_structure_factor->getParameter("Intensity") << " Radius = " << fitalg_structure_factor->getParameter("Radius") << " Decay = " << fitalg_structure_factor->getParameter("Decay") << "\n"; // only for debugging purposes

} // runDiffRotDiscreteCircleInelasticTest


/// returns a real value from a uniform distribution
double random_value(const double & a, const double & b)
Expand Down

0 comments on commit 62af761

Please sign in to comment.