Skip to content

Commit

Permalink
Fix a bug. Refs #6019.
Browse files Browse the repository at this point in the history
  • Loading branch information
@wdzhou
wdzhou committed Dec 31, 2012
1 parent 7f620af commit c77517c
Show file tree
Hide file tree
Showing 2 changed files with 46 additions and 45 deletions.
9 changes: 5 additions & 4 deletions Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/LeBailFit2.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -120,10 +120,10 @@ namespace CurveFitting
map<string, Parameter> parammap, bool recalpeakintesity);

/// LeBailFit
void execLeBailFit(size_t workspaceindex);
void execLeBailFit();

/// Do 1 iteration in Le Bail fit
bool do1StepLeBailFit(size_t workspaceindex, std::map<std::string, Parameter>& parammap);
bool do1StepLeBailFit(std::map<std::string, Parameter>& parammap);

/// Set up Lebail
void setLeBailFitParameters();
Expand Down Expand Up @@ -200,7 +200,8 @@ namespace CurveFitting
void writeFakedDataToOutputWS(size_t workspaceindex, int functionmode);

/// Write out (domain, values) to output workspace
void writeToOutputWorkspace(API::FunctionDomain1DVector domain, API::FunctionValues values);
void writeToOutputWorkspace(size_t wsindex, FunctionDomain1DVector domain, FunctionValues values);
// void writeToOutputWorkspace(API::FunctionDomain1DVector domain, API::FunctionValues values);

/// Write input data and difference to output workspace
void writeInputDataNDiff(size_t workspaceindex, API::FunctionDomain1DVector domain);
Expand Down Expand Up @@ -304,7 +305,7 @@ namespace CurveFitting
void doResultStatistics();

/// ============================= =========================== ///
size_t mWSIndexToWrite;
// size_t mWSIndexToWrite;

/// Map to store peak group information: key (int) = (hkl)^2; value = group ID
std::map<int, size_t> mPeakGroupMap;
Expand Down
82 changes: 41 additions & 41 deletions Code/Mantid/Framework/CurveFitting/src/LeBailFit2.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -127,7 +127,7 @@ namespace CurveFitting
this->declareProperty("BackgroundType", "Polynomial", bkgdvalidator, "Background type");

// Background function order
this->declareProperty("BackgroundFunctionOrder", 12, "Order of background function.");
// this->declareProperty("BackgroundFunctionOrder", 12, "Order of background function.");

// Input background parameters (array)
this->declareProperty(new Kernel::ArrayProperty<double>("BackgroundParameters"),
Expand Down Expand Up @@ -244,7 +244,7 @@ namespace CurveFitting
g_log.notice() << "Function: Do LeBail Fit." << std::endl;
if (inputparamcorrect)
{
execLeBailFit(m_wsIndex);
execLeBailFit();
}
else
{
Expand Down Expand Up @@ -446,7 +446,7 @@ namespace CurveFitting
//----------------------------------------------------------------------------------------------
/** LeBail Fitting for one self-consistent iteration
*/
void LeBailFit2::execLeBailFit(size_t workspaceindex)
void LeBailFit2::execLeBailFit()
{
// FIXME: m_maxStep should be an instance variable and evaluate in exec as an input property
int m_maxSteps = 1;
Expand All @@ -458,7 +458,7 @@ namespace CurveFitting
// 2. Do 1 iteration of LeBail fit
for (int istep = 0; istep < m_maxSteps; ++istep)
{
this->do1StepLeBailFit(workspaceindex, parammap);
this->do1StepLeBailFit(parammap);
}

// 3. Output
Expand Down Expand Up @@ -536,50 +536,44 @@ namespace CurveFitting
* a) Calculate pattern for peak intensities
* b) Set peak intensities
*/
// TODO Release 2.4: Arguments: size_t workspaceindex, std::map<std::string, Parameter>& paramma, size_t istep
bool LeBailFit2::do1StepLeBailFit(size_t workspaceindex, std::map<std::string, Parameter>& parammap)
bool LeBailFit2::do1StepLeBailFit(map<string, Parameter>& parammap)
{
// 1. Generate domain and value
const std::vector<double> vecX = m_dataWS->readX(workspaceindex);
const std::vector<double> vecX = m_dataWS->readX(m_wsIndex);
API::FunctionDomain1DVector domain(vecX);
API::FunctionValues values(domain);

// 2. Calculate peak intensity and etc.
// FIXME Release 2.4 Single out the process to calculate intensity.
// FIXME Move calculation of peak out of UnitLeBailFit
bool calpeakintensity = true;
this->calculateDiffractionPattern(m_dataWS, workspaceindex, domain, values, parammap, calpeakintensity);
vector<double> allpeaksvalues(vecX.size(), 0.0);
calculatePeaksIntensities(m_dataWS, m_wsIndex, false, allpeaksvalues);
// calculateDiffractionPattern(m_dataWS, workspaceindex, domain, values, parammap, calpeakintensity);

// a) Apply initial calculated result to output workspace
mWSIndexToWrite = 5;
writeToOutputWorkspace(domain, values);
writeToOutputWorkspace(5, domain, values);

// b) Calculate input background
m_backgroundFunction->function(domain, values);
mWSIndexToWrite = 6;
writeToOutputWorkspace(domain, values);
writeToOutputWorkspace(6, domain, values);

// 3. Construct the tie. 2-level loop. (1) peak parameter (2) peak
// TODO Release 2.4: setLeBailFitParameters(istep)
this->setLeBailFitParameters();

// 4. Construct the Fit
this->fitLeBailFunction(workspaceindex, parammap);
this->fitLeBailFunction(m_wsIndex, parammap);

// TODO (1) Calculate Rwp, Chi^2, .... for the fitted pattern.

// 5. Do calculation again and set the output
// FIXME Move this part our of UnitLeBailFit
calpeakintensity = true;
bool calpeakintensity = true;
API::FunctionValues newvalues(domain);
this->calculateDiffractionPattern(m_dataWS, workspaceindex, domain, newvalues, parammap, calpeakintensity);
this->calculateDiffractionPattern(m_dataWS, m_wsIndex, domain, newvalues, parammap, calpeakintensity);

// Add final calculated value to output workspace
mWSIndexToWrite = 1;
writeToOutputWorkspace(domain, newvalues);
writeToOutputWorkspace(1, domain, newvalues);

// Add original data and
writeInputDataNDiff(workspaceindex, domain);
writeInputDataNDiff(m_wsIndex, domain);

return true;
}
Expand Down Expand Up @@ -2065,29 +2059,35 @@ void LeBailFit2::exportBraggPeakParameterToTable()
return;
}

//-----------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------
/** Write data (domain, values) to one specified spectrum of output workspace
*/
void LeBailFit2::writeToOutputWorkspace(API::FunctionDomain1DVector domain, API::FunctionValues values)
void LeBailFit2::writeToOutputWorkspace(size_t wsindex, FunctionDomain1DVector domain, FunctionValues values)
{
if (m_outputWS->getNumberHistograms() <= mWSIndexToWrite)
// Check workspace index
if (m_outputWS->getNumberHistograms() <= wsindex)
{
g_log.error() << "LeBailFit.writeToOutputWorkspace. Try to write to spectrum " << mWSIndexToWrite << " out of range = "
<< m_outputWS->getNumberHistograms() << std::endl;
throw std::invalid_argument("Try to write to a spectrum out of range.");
stringstream errss;
errss << "LeBailFit.writeToOutputWorkspace. Try to write to spectrum " << wsindex << " out of range = "
<< m_outputWS->getNumberHistograms();
g_log.error(errss.str());
throw runtime_error(errss.str());
}

// Write X
for (size_t i = 0; i < domain.size(); ++i)
{
m_outputWS->dataX(mWSIndexToWrite)[i] = domain[i];
m_outputWS->dataX(wsindex)[i] = domain[i];
}

// Write Y & E (square root of Y)
for (size_t i = 0; i < values.size(); ++i)
{
m_outputWS->dataY(mWSIndexToWrite)[i] = values[i];
m_outputWS->dataY(wsindex)[i] = values[i];
if (fabs(values[i]) > 1.0)
m_outputWS->dataE(mWSIndexToWrite)[i] = std::sqrt(fabs(values[i]));
m_outputWS->dataE(wsindex)[i] = std::sqrt(fabs(values[i]));
else
m_outputWS->dataE(mWSIndexToWrite)[i] = 1.0;
m_outputWS->dataE(wsindex)[i] = 1.0;
}

return;
Expand Down Expand Up @@ -2380,14 +2380,13 @@ void LeBailFit2::createOutputDataWorkspace(size_t workspaceindex, FunctionMode f
}


// ============================ Random Walk Suite ==============================
//-----------------------------------------------------------------------------
// ====================================== Random Walk Suite ======================================
//------------------------------------------------------------------------------------------------
/** Refine instrument parameters by random walk algorithm (MC)
*
* @param wsindex : workspace index of the diffraction data (observed)
*/
void LeBailFit2::execRandomWalkMinimizer(size_t maxcycles, size_t wsindex,
map<string, Parameter>& parammap)
void LeBailFit2::execRandomWalkMinimizer(size_t maxcycles, size_t wsindex, map<string, Parameter>& parammap)
{
// 1. Initialization
const MantidVec& vecInY = m_dataWS->readY(wsindex);
Expand Down Expand Up @@ -2468,8 +2467,8 @@ void LeBailFit2::execRandomWalkMinimizer(size_t maxcycles, size_t wsindex,

// 4. Random walk loops
// generate some MC trace structure
vector<double> vecIndex(maxcycles);
vector<double> vecRwp(maxcycles);
vector<double> vecIndex(maxcycles+1);
vector<double> vecRwp(maxcycles+1);
size_t numinvalidmoves = 0;
size_t numacceptance = 0;
bool prevcyclebetterrwp = true;
Expand Down Expand Up @@ -2590,7 +2589,7 @@ void LeBailFit2::execRandomWalkMinimizer(size_t maxcycles, size_t wsindex,
// 5. Sum up
// a) Summary output
g_log.notice() << "[SUMMARY] Random-walk Rwp: Starting = " << startrwp << ", Best = " << m_bestRwp
<< " @ Step = " << m_bestMCStep << ", Final Rwp = " << currwp
<< " @ Step = " << m_bestMCStep << "; Last Step Rwp = " << currwp
<< ", Rp = " << currp
<< ", Acceptance ratio = " << double(numacceptance)/double(maxcycles*m_numMCGroups) << endl;

Expand Down Expand Up @@ -3175,6 +3174,7 @@ void LeBailFit2::bookKeepBestMCResult(map<string, Parameter> parammap, vector<do
{
// If not be initialized, initialize it!
m_bestParameters = parammap;
// copyParameterMap(parammap, m_bestParameters);
}
else
{
Expand All @@ -3191,7 +3191,7 @@ void LeBailFit2::bookKeepBestMCResult(map<string, Parameter> parammap, vector<do
return;
}

//-----------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------
/** Apply the value of parameters in the source to target
*/
void LeBailFit2::applyParameterValues(map<string, Parameter>& srcparammap, map<string, Parameter>& tgtparammap)
Expand Down Expand Up @@ -3221,7 +3221,7 @@ void LeBailFit2::applyParameterValues(map<string, Parameter>& srcparammap, map<s

//=============== Background Functions ========================================

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------
/** Re-fit background according to the new values
*
* @param wsindex raw data's workspace index
Expand Down

0 comments on commit c77517c

Please sign in to comment.