Skip to content

Commit

Permalink
Update inelastic tests to run with shiift param
Browse files Browse the repository at this point in the history 
Refs #10189
  • Loading branch information
@DanNixon
DanNixon committed Mar 24, 2015
1 parent c0ec79c commit 30e5e3a
Show file tree
Hide file tree
Showing 2 changed files with 148 additions and 104 deletions.
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 30e5e3a

Please sign in to comment.