Skip to content

Commit

Permalink
Implement parallelization over events in box. Refs #5397
Browse files Browse the repository at this point in the history
  • Loading branch information
@martyngigg
martyngigg committed Aug 2, 2012
1 parent 5a175fb commit 7010f2f
Show file tree
Hide file tree
Showing 5 changed files with 78 additions and 43 deletions.
5 changes: 5 additions & 0 deletions ...id/Framework/MDAlgorithms/inc/MantidMDAlgorithms/Quantification/MDResolutionConvolution.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,6 +78,11 @@ namespace Mantid
* @param workspace :: The MD that will be used for the fit
*/
virtual void preprocess(const API::IMDEventWorkspace_const_sptr & workspace) { UNUSED_ARG(workspace); }
/**
* Called before any fit/simulation is started to tell the function how many threads will be used.
* Default does nothing.
*/
virtual void useNumberOfThreads(const int) {}

/// Declares the parameters. Overridden here to ensure that concrete models override it
void declareAttributes();
Expand Down
13 changes: 8 additions & 5 deletions ...rk/MDAlgorithms/inc/MantidMDAlgorithms/Quantification/Resolution/TobyFitResolutionModel.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,6 +83,9 @@ namespace Mantid

private:
DISABLE_COPY_AND_ASSIGN(TobyFitResolutionModel);

/// Informs the function how many threads will be processing it
void useNumberOfThreads(const int nthreads);
/// Cache detector observations
void preprocess(const API::IMDEventWorkspace_const_sptr & workspace);
/// Declare function attributes
Expand Down Expand Up @@ -128,15 +131,15 @@ namespace Mantid
int m_mosaicActive;

/// A pre-sized matrix for the resolution coefficients
mutable TobyFitBMatrix m_bmatrix;
mutable std::vector<TobyFitBMatrix> m_bmatrix;
/// A pre-sized vector for the randomly generated points
mutable TobyFitYVector m_yvector;
mutable std::vector<TobyFitYVector> m_yvector;
/// The generated value of the in-place mosaic (eta_2)
mutable double m_etaInPlane;
mutable std::vector<double> m_etaInPlane;
/// The generated value of the in-place mosaic (eta_3)
mutable double m_etaOutPlane;
mutable std::vector<double> m_etaOutPlane;
/// A pre-sized vector for the QE position to be evaluated
mutable std::vector<double> m_deltaQE;
mutable std::vector<std::vector<double>> m_deltaQE;

/// Cache of detector observations
std::map<std::pair<int, detid_t>, Observation*> m_observations;
Expand Down
2 changes: 0 additions & 2 deletions .../Framework/MDAlgorithms/inc/MantidMDAlgorithms/Quantification/Resolution/TobyFitYVector.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -91,8 +91,6 @@ namespace Mantid
const QOmegaPoint & qOmega);

private:
DISABLE_COPY_AND_ASSIGN(TobyFitYVector);

/// Sample from moderator time distribution
void calculateModeratorTime();
/// Aperature contribution
Expand Down
88 changes: 57 additions & 31 deletions Code/Mantid/Framework/MDAlgorithms/src/Quantification/Resolution/TobyFitResolutionModel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,8 @@ namespace Mantid
: MDResolutionConvolution(), m_randGen(new Kernel::SobolSequence(TobyFitYVector::variableCount() + 2)), // For eta
m_activeAttrValue(1),
m_mcLoopMin(100), m_mcLoopMax(1000), m_mcRelErrorTol(1e-5), m_mosaicActive(1),
m_bmatrix(), m_yvector(), m_etaInPlane(0.0), m_etaOutPlane(0.0), m_deltaQE(4, 0.0)
m_bmatrix(), m_yvector(), m_etaInPlane(), m_etaOutPlane(), m_deltaQE(),
m_observations()
{
}

Expand All @@ -51,7 +52,8 @@ namespace Mantid
: MDResolutionConvolution(fittedFunction, fgModel), m_randGen(new Kernel::SobolSequence(TobyFitYVector::variableCount())),
m_activeAttrValue(1),
m_mcLoopMin(100), m_mcLoopMax(1000), m_mcRelErrorTol(1e-5), m_mosaicActive(1),
m_bmatrix(), m_yvector(), m_observations()
m_bmatrix(), m_yvector(), m_etaInPlane(), m_etaOutPlane(), m_deltaQE(),
m_observations()
{
}

Expand Down Expand Up @@ -95,7 +97,8 @@ namespace Mantid
generateIntegrationVariables(detObservation, qOmega);
calculatePerturbedQE(detObservation, qOmega);
// Compute weight from the foreground at this point
const double weight = foregroundModel().scatteringIntensity(detObservation.experimentInfo(), m_deltaQE);
const double weight = foregroundModel().scatteringIntensity(detObservation.experimentInfo(),
m_deltaQE[PARALLEL_THREAD_NUMBER]);

// Add on this contribution to the average
sumSigma += weight;
Expand All @@ -113,6 +116,20 @@ namespace Mantid
//---------------------------------------------------------------------------------
// Private members
//---------------------------------------------------------------------------------
/**
* Sets up the function to cope with the given number of threads processing it at once
* @param nthreads :: The maximum number of threads that will be used to evaluate the
* function
*/
void TobyFitResolutionModel::useNumberOfThreads(const int nthreads)
{
m_bmatrix = std::vector<TobyFitBMatrix>(nthreads, TobyFitBMatrix());
m_yvector = std::vector<TobyFitYVector>(nthreads, TobyFitYVector());
m_etaInPlane = std::vector<double>(nthreads, 0.0);
m_etaOutPlane = std::vector<double>(nthreads, 0.0);
m_deltaQE = std::vector<std::vector<double>>(nthreads, std::vector<double>(4, 0.0));
}

/**
* Called before any fit/simulation is started to allow caching of
* frequently used parameters
Expand Down Expand Up @@ -182,8 +199,13 @@ namespace Mantid
else if(name == MC_MAX_NAME) m_mcLoopMax = value.asInt();
else if(name == MC_LOOP_TOL) m_mcRelErrorTol = value.asDouble();
else if(name == CRYSTAL_MOSAIC) m_mosaicActive = value.asInt();
else m_yvector.setAttribute(name, value.asInt());

else
{
for(auto iter = m_yvector.begin(); iter != m_yvector.end(); ++iter)
{
iter->setAttribute(name, value.asInt());
}
}
}

/**
Expand All @@ -205,7 +227,7 @@ namespace Mantid
void TobyFitResolutionModel::calculateResolutionCoefficients(const Observation & observation,
const QOmegaPoint & eventPoint) const
{
m_bmatrix.recalculate(observation, eventPoint);
m_bmatrix[PARALLEL_THREAD_NUMBER].recalculate(observation, eventPoint);
}

/**
Expand All @@ -217,7 +239,7 @@ namespace Mantid
const QOmegaPoint & eventPoint) const
{
const std::vector<double> & randomNums = m_randGen->nextPoint();
const size_t nvars = m_yvector.recalculate(randomNums, observation, eventPoint);
const size_t nvars = m_yvector[PARALLEL_THREAD_NUMBER].recalculate(randomNums, observation, eventPoint);

// Calculate crystal mosaic contribution
if(m_mosaicActive)
Expand All @@ -227,12 +249,12 @@ namespace Mantid
const double r2 = randomNums[nvars+1];
const double etaSig = observation.experimentInfo().run().getLogAsSingleValue("eta_sigma");

m_etaInPlane = etaSig*prefactor*std::cos(2.0*M_PI*r2);
m_etaOutPlane = etaSig*prefactor*std::sin(2.0*M_PI*r2);
m_etaInPlane[PARALLEL_THREAD_NUMBER] = etaSig*prefactor*std::cos(2.0*M_PI*r2);
m_etaOutPlane[PARALLEL_THREAD_NUMBER] = etaSig*prefactor*std::sin(2.0*M_PI*r2);
}
else
{
m_etaInPlane = m_etaOutPlane = 0.0;
m_etaInPlane[PARALLEL_THREAD_NUMBER] = m_etaOutPlane[PARALLEL_THREAD_NUMBER] = 0.0;
}
}

Expand All @@ -255,15 +277,16 @@ namespace Mantid

double xVec0(0.0), xVec1(0.0), xVec2(0.0),
xVec3(0.0), xVec4(0.0), xVec5(0.0);
const std::vector<double> & yvalues = m_yvector.values();
const std::vector<double> & yvalues = m_yvector[PARALLEL_THREAD_NUMBER].values();
const TobyFitBMatrix & bmatrix = m_bmatrix[PARALLEL_THREAD_NUMBER];
for(unsigned int i =0; i < TobyFitYVector::variableCount(); ++i)
{
xVec0 += m_bmatrix[0][i] * yvalues[i];
xVec1 += m_bmatrix[1][i] * yvalues[i];
xVec2 += m_bmatrix[2][i] * yvalues[i];
xVec3 += m_bmatrix[3][i] * yvalues[i];
xVec4 += m_bmatrix[4][i] * yvalues[i];
xVec5 += m_bmatrix[5][i] * yvalues[i];
xVec0 += bmatrix[0][i] * yvalues[i];
xVec1 += bmatrix[1][i] * yvalues[i];
xVec2 += bmatrix[2][i] * yvalues[i];
xVec3 += bmatrix[3][i] * yvalues[i];
xVec4 += bmatrix[4][i] * yvalues[i];
xVec5 += bmatrix[5][i] * yvalues[i];
}

// Convert to dQ in lab frame
Expand All @@ -282,24 +305,27 @@ namespace Mantid
const double dqlab0 = (L00*xVec3 + L01*xVec4 + L02*xVec5)/determinant;
const double dqlab1 = (L10*xVec3 + L11*xVec4 + L12*xVec5)/determinant;
const double dqlab2 = (L20*xVec3 + L21*xVec4 + L22*xVec5)/determinant;
m_deltaQE[0] = (xVec0 - dqlab0);
m_deltaQE[1] = (xVec1 - dqlab1);
m_deltaQE[2] = (xVec2 - dqlab2);
std::vector<double> & deltaQE = m_deltaQE[PARALLEL_THREAD_NUMBER];
deltaQE[0] = (xVec0 - dqlab0);
deltaQE[1] = (xVec1 - dqlab1);
deltaQE[2] = (xVec2 - dqlab2);

const double efixed = observation.getEFixed();
const double wi = std::sqrt(efixed/PhysicalConstants::E_mev_toNeutronWavenumberSq);
const double wf = std::sqrt((efixed - qOmega.deltaE)/PhysicalConstants::E_mev_toNeutronWavenumberSq);
m_deltaQE[3] = 4.1442836 * (wi*xVec2 - wf*xVec5);
deltaQE[3] = 4.1442836 * (wi*xVec2 - wf*xVec5);

// Add on the nominal Q
m_deltaQE[0] += qOmega.qx;
m_deltaQE[1] += qOmega.qy;
m_deltaQE[2] += qOmega.qz;
m_deltaQE[3] += qOmega.deltaE;
deltaQE[0] += qOmega.qx;
deltaQE[1] += qOmega.qy;
deltaQE[2] += qOmega.qz;
deltaQE[3] += qOmega.deltaE;

if(m_mosaicActive)
{
const double qx(m_deltaQE[0]),qy(m_deltaQE[1]),qz(m_deltaQE[1]);
const double etaInPlane = m_etaInPlane[PARALLEL_THREAD_NUMBER];
const double etaOutPlane = m_etaOutPlane[PARALLEL_THREAD_NUMBER];
const double qx(deltaQE[0]),qy(deltaQE[1]),qz(deltaQE[1]);
const double qipmodSq = qy*qy + qz*qz;
const double qmod = std::sqrt(qx*qx + qipmodSq);
static const double small(1e-10);
Expand All @@ -308,14 +334,14 @@ namespace Mantid
const double qipmod = std::sqrt(qipmodSq);
if(qipmod > small)
{
m_deltaQE[0] -= qipmod*m_etaInPlane;
m_deltaQE[1] += ((qx*qy)*m_etaInPlane - (qz*qmod)*m_etaOutPlane)/qipmod;
m_deltaQE[2] += ((qx*qz)*m_etaInPlane + (qy*qmod)*m_etaOutPlane)/qipmod;
deltaQE[0] -= qipmod*etaInPlane;
deltaQE[1] += ((qx*qy)*etaInPlane - (qz*qmod)*etaOutPlane)/qipmod;
deltaQE[2] += ((qx*qz)*etaInPlane + (qy*qmod)*etaOutPlane)/qipmod;
}
else
{
m_deltaQE[1] += qmod*m_etaInPlane;
m_deltaQE[2] += qmod*m_etaOutPlane;
deltaQE[1] += qmod*etaInPlane;
deltaQE[2] += qmod*etaOutPlane;
}
}
}
Expand Down
13 changes: 8 additions & 5 deletions Code/Mantid/Framework/MDAlgorithms/src/Quantification/ResolutionConvolvedCrossSection.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -93,6 +93,7 @@ namespace Mantid
throw std::invalid_argument("ResolutionConvolvedCrossSection can only be used with MD event workspaces");
}
IFunctionMD::setWorkspace(workspace);
m_convolution->useNumberOfThreads(PARALLEL_GET_MAX_THREADS);
m_convolution->preprocess(m_workspace);
}

Expand All @@ -103,8 +104,9 @@ namespace Mantid
{
throw std::invalid_argument("Unexpected domain in IFunctionMD");
}

dmd->reset();
size_t i=0;
size_t i = 0;
for(const API::IMDIterator* r = dmd->getNextIterator(); r != NULL; r = dmd->getNextIterator())
{
values.setCalculated(i,functionMD(*r));
Expand All @@ -123,18 +125,19 @@ namespace Mantid
*/
double ResolutionConvolvedCrossSection::functionMD(const Mantid::API::IMDIterator& box) const
{
assert(m_convolution);
const size_t numEvents = box.getNumEvents();
if(numEvents == 0) return 0.0;

assert(m_convolution);

double signal(0.0);
// loop over each MDPoint in current MDBox
for(size_t j = 0; j < numEvents; j++)
PARALLEL_FOR_NO_WSP_CHECK()
for(size_t j = 0; j < numEvents; ++j)
{
double contribution = m_convolution->signal(box, box.getInnerRunIndex(j), j);
PARALLEL_ATOMIC
signal += contribution;
}

// Return the mean
return signal/static_cast<double>(numEvents);
}
Expand Down

0 comments on commit 7010f2f

Please sign in to comment.