Skip to content

Commit

Permalink
Refactored to use centralized version of much of the code. Refs #4208.
Browse files Browse the repository at this point in the history 
The constructs introduced in the event branch that were copied into the histogram
version are now centralized. The new central versions are now back into the event
branch.
  • Loading branch information
@peterfpeterson
peterfpeterson committed Dec 2, 2011
1 parent 7ad2d9d commit c9d4c0d
Show file tree
Hide file tree
Showing 2 changed files with 114 additions and 128 deletions.
7 changes: 6 additions & 1 deletion Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/DiffractionFocussing2.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -100,6 +100,8 @@ class DLLExport DiffractionFocussing2: public API::Algorithm
void exec();
void cleanup();

std::size_t setupGroupToWSIndices();

// For events
void execEvent();

Expand All @@ -117,7 +119,6 @@ class DLLExport DiffractionFocussing2: public API::Algorithm
/// Shared pointer to the event workspace
DataObjects::EventWorkspace_const_sptr m_eventW;


// This map does not need to be ordered, just a lookup for udet
/// typedef for the storage of the UDET-group mapping
typedef std::map<detid_t, int> udet2groupmap;
Expand All @@ -139,6 +140,10 @@ class DLLExport DiffractionFocussing2: public API::Algorithm
int nHist;
/// Number of points in the 2D workspace
int nPoints;
/// Mapping of group number to vector of inputworkspace indices.
std::vector< std::vector<std::size_t> > m_wsIndices;
/// List of valid group numbers
std::vector<int> m_validGroups;

};

Expand Down
235 changes: 108 additions & 127 deletions Code/Mantid/Framework/Algorithms/src/DiffractionFocussing2.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -90,6 +90,8 @@ void DiffractionFocussing2::cleanup()
std::vector<int>().swap(groupAtWorkspaceIndex);
group2xvector.clear();
group2wgtvector.clear();
this->m_validGroups.clear();
this->m_wsIndices.clear();
}


Expand Down Expand Up @@ -145,6 +147,8 @@ void DiffractionFocussing2::exec()
// It also initializes the groupAtWorkspaceIndex[] array.
determineRebinParameters();

size_t totalHistProcess = this->setupGroupToWSIndices();

// determine event workspace min/max tof
double eventXMin = 0.;
double eventXMax = 0.;
Expand Down Expand Up @@ -173,50 +177,14 @@ void DiffractionFocussing2::exec()
// Now the real work
MantidVec limits(2), weights_default(1,1.0), emptyVec(1,0.0), EOutDummy(nPoints); // Vectors for use with the masking stuff

// The output spectrum, will be set at the first group
Progress * prog;
prog = new API::Progress(this,0.2,0.25,nHist);

// determine the workspace indices for each group
vector< vector<size_t> > ws_indices(nGroups+1);
for (size_t wi=0;wi<static_cast<size_t>(nHist);wi++)
{
//i is the workspace index (of the input)
const int group = groupAtWorkspaceIndex[wi];
if (group < 1) // Not in a group, or invalid group #
continue;

// resize the ws_indices if it is not big enough
if (ws_indices.size() < static_cast<size_t>(group+1))
{
ws_indices.resize(group + 1);
}

// Also record a list of workspace indices
ws_indices[group].push_back(wi);
prog->reportIncrement(1, "Pre-counting");
}

// initialize a vector of the valid group numbers
vector<int> valid_groups;
size_t totalHistProcess = 0;
for (size_t i = 0; i < ws_indices.size(); i++)
{
if (!(ws_indices[i].empty()))
{
valid_groups.push_back(static_cast<int>(i));
totalHistProcess += ws_indices[i].size();
}
}

// loop over groups
delete prog;
prog = new API::Progress(this, 0.25, 0.95, static_cast<int>(totalHistProcess) + nGroups);
Progress * prog;
prog = new API::Progress(this, 0.2, 0.95, static_cast<int>(totalHistProcess) + nGroups);
PARALLEL_FOR1(m_matrixInputW)
for (int outWorkspaceIndex = 0; outWorkspaceIndex < static_cast<int>(valid_groups.size()); outWorkspaceIndex++)
for (int outWorkspaceIndex = 0; outWorkspaceIndex < static_cast<int>(m_validGroups.size()); outWorkspaceIndex++)
{
PARALLEL_START_INTERUPT_REGION
int group = valid_groups[outWorkspaceIndex];
int group = m_validGroups[outWorkspaceIndex];

// Get the group
group2vectormap::iterator it=group2xvector.find(group);
Expand All @@ -242,7 +210,7 @@ void DiffractionFocussing2::exec()
MantidVec groupWgt(nPoints,0.0);

// loop through the contributing histograms
std::vector<size_t> indices = ws_indices[group];
std::vector<size_t> indices = m_wsIndices[group];
const size_t groupSize = indices.size();
for (size_t i=0; i<groupSize; i++)
{
Expand Down Expand Up @@ -381,55 +349,69 @@ void DiffractionFocussing2::execEvent()
g_log.debug("Focussing EventWorkspace in-place.");
g_log.debug() << nGroups << " groups found in .cal file (counting group 0).\n";

Geometry::Instrument_const_sptr instrument = m_eventW->getInstrument();
Geometry::IObjComponent_const_sptr source;
Geometry::IObjComponent_const_sptr sample;
if (instrument != NULL)
EventType eventWtype = m_eventW->getEventType();

Progress * prog;
prog = new Progress(this,0.2,0.25,nGroups);

// determine precount size
vector<size_t> size_required(this->m_validGroups.size(),0);
int totalHistProcess = 0;
for (size_t iGroup = 0; iGroup < this->m_validGroups.size(); iGroup++)
{
source = instrument->getSource();
sample = instrument->getSample();
const int group = this->m_validGroups[iGroup];
const vector<size_t> &indices = this->m_wsIndices[group];

totalHistProcess += static_cast<int>(indices.size());
for (vector<size_t>::const_iterator index = indices.begin();
index != indices.end(); ++index)
{
size_required[iGroup] += m_eventW->getEventList(*index).getNumberEvents();
}
prog->report(1, "Pre-counting");
}

Progress * prog;
prog = new Progress(this,0.2,0.35,nHist);
// ------------- Pre-allocate Event Lists ----------------------------
delete prog; prog = new Progress(this,0.25,0.3,totalHistProcess);

bool checkForMask = (instrument != NULL) && (source != NULL) && (sample != NULL);
EventType eventWtype = m_eventW->getEventType();
// This creates and reserves the space required
for (size_t iGroup = 0; iGroup < this->m_validGroups.size(); iGroup++)
{
const int group = this->m_validGroups[iGroup];
EventList & groupEL = out->getOrAddEventList(iGroup);
groupEL.switchTo(eventWtype);
groupEL.reserve(size_required[iGroup]);
groupEL.clearDetectorIDs();
groupEL.setSpectrumNo(group);
prog->reportIncrement(1, "Allocating");
}

// ----------- Focus ---------------
delete prog; prog = new Progress(this,0.3,0.9,totalHistProcess);

if (nGroups == 1)
if (this->m_validGroups.size() == 1)
{
g_log.information() << "Performing focussing on a single group\n";
// Special case of a single group - parallelize differently
EventList & groupEL = out->getOrAddEventList(0);
groupEL.switchTo(eventWtype);
groupEL.reserve(m_eventW->getNumberEvents()); // does this slow it all down?

// Only one group, spec # is = 1
groupEL.setSpectrumNo(1);
// Make sure you start with no detector IDs.
groupEL.clearDetectorIDs();
const int group = m_validGroups[0];
std::vector<size_t> indices = this->m_wsIndices[group];

int chunkSize = 200;

PRAGMA_OMP(parallel for schedule(dynamic, 1) )
for (int wiChunk=0;wiChunk<(nHist/chunkSize)+1;wiChunk++)
for (int wiChunk=0;wiChunk<(totalHistProcess/chunkSize)+1;wiChunk++)
{
PARALLEL_START_INTERUPT_REGION

// Perform in chunks for more efficiency
int max = (wiChunk+1)*chunkSize;
if (max > nHist) max = nHist;
if (max > totalHistProcess) max = totalHistProcess;

// precalculate output size
/* size_t numEventsInChunk(0);
for (int wi=wiChunk*chunkSize; wi < max; wi++)
{
// Check for masking. TODO: Most of the pointer checks are redundant
if (checkForMask)
{
Geometry::IDetector_const_sptr det = eventW->getDetector(static_cast<size_t>(wi));
if ( det->isMasked() ) continue; // should be cached
}
const int group = groupAtWorkspaceIndex[wi];
if (group == 1)
{
Expand All @@ -446,18 +428,8 @@ void DiffractionFocussing2::execEvent()
// process the chunk
for (int wi=wiChunk*chunkSize; wi < max; wi++)
{
// Check for masking. TODO: Most of the pointer checks are redundant
if (checkForMask)
{
Geometry::IDetector_const_sptr det = m_eventW->getDetector(static_cast<size_t>(wi));
if ( det->isMasked() ) continue;
}
const int group = groupAtWorkspaceIndex[wi];
if (group == 1)
{
// Accumulate the chunk
chunkEL += m_eventW->getEventList(wi);
}
// Accumulate the chunk
chunkEL += m_eventW->getEventList(wi);
}

// Rejoin the chunk with the rest.
Expand All @@ -469,56 +441,24 @@ void DiffractionFocussing2::execEvent()
PARALLEL_END_INTERUPT_REGION
}
PARALLEL_CHECK_INTERUPT_REGION

out->doneAddingEventLists();
}
else
{
// ------ PARALLELIZE BY GROUPS -------------------------
// ------------- Pre-allocate Event Lists ----------------------------
std::vector< std::vector<int> > ws_indices(nGroups+1);
std::vector<size_t> size_required(nGroups+1,0);

for (int wi=0;wi<nHist;wi++)
{
if (checkForMask)
{
Geometry::IDetector_const_sptr det = m_eventW->getDetector(static_cast<size_t>(wi));
if ( det->isMasked() ) continue;
}
//i is the workspace index (of the input)
const int group = groupAtWorkspaceIndex[wi];
if (group < 1 || group > nGroups) // Not in a group, or invalid group #
continue;
size_required[group] += m_eventW->getEventList(wi).getNumberEvents();
// Also record a list of workspace indices
ws_indices[group].push_back(wi);
prog->reportIncrement(1, "Pre-counting");
}

delete prog; prog = new Progress(this,0.15,0.3,nGroups);
// This creates and reserves the space required
for (int group=1; group<nGroups+1; group++)
{
out->getOrAddEventList(group-1).reserve(size_required[group]);
out->getOrAddEventList(group-1).clearDetectorIDs();
out->getOrAddEventList(group-1).setSpectrumNo(group);
prog->reportIncrement(1, "Allocating");
}

// ----------- Focus ---------------
delete prog; prog = new Progress(this,0.40,0.9,nHist);
int nValidGroups = static_cast<int>(this->m_validGroups.size());
PARALLEL_FOR1(m_eventW)
for (int group=1; group<nGroups+1; group++)
for (int iGroup = 0; iGroup < nValidGroups; iGroup++)
{
PARALLEL_START_INTERUPT_REGION
std::vector<int> indices = ws_indices[group];
const int group = this->m_validGroups[iGroup];
std::vector<size_t> indices = this->m_wsIndices[group];
for (size_t i=0; i<indices.size(); i++)
{
int wi = indices[i];
size_t wi = indices[i];

//In workspace index group-1, put what was in the OLD workspace index wi
out->getOrAddEventList(group-1) += m_eventW->getEventList(wi);
//In workspace index iGroup, put what was in the OLD workspace index wi
out->getOrAddEventList(iGroup) += m_eventW->getEventList(wi);

prog->reportIncrement(1, "Appending Lists");

Expand All @@ -532,30 +472,29 @@ void DiffractionFocussing2::execEvent()
PARALLEL_END_INTERUPT_REGION
}
PARALLEL_CHECK_INTERUPT_REGION

//Finalize the maps
out->doneAddingEventLists();

} // (done with parallel by groups)

//Finalize the maps
out->doneAddingEventLists();

//Now that the data is cleaned up, go through it and set the X vectors to the input workspace we first talked about.
delete prog; prog = new Progress(this,0.9,1.0,nGroups);
for (int g=1; g < nGroups+1; g++)
for (size_t workspaceIndex = 0; workspaceIndex < this->m_validGroups.size(); workspaceIndex++)
{
const int group = this->m_validGroups[workspaceIndex];
//Now this is the workspace index of that group; simply 1 offset
int workspaceIndex = g-1;
prog->reportIncrement(1, "Setting X");

if (workspaceIndex >= static_cast<int>(out->getNumberHistograms()))
if (workspaceIndex >= out->getNumberHistograms())
{
g_log.warning() << "Warning! Invalid workspace index found for group # " << g << ". Histogram will be empty.\n";
g_log.warning() << "Warning! Invalid workspace index found for group # " << group << ". Histogram will be empty.\n";
continue;
}

//Now you set the X axis to the X you saved before.
if (group2xvector.size() > 0)
{
group2vectormap::iterator git = group2xvector.find(g);
group2vectormap::iterator git = group2xvector.find(group);
if (git != group2xvector.end())
out->setX(workspaceIndex, (git->second) );
else
Expand Down Expand Up @@ -722,6 +661,48 @@ void DiffractionFocussing2::determineRebinParameters()
return;
}

/***
* Configure the mapping of output group to list of input workspace
* indices, and the list of valid group numbers.
*
* @return the total number of input histograms that will be read.
*/
size_t DiffractionFocussing2::setupGroupToWSIndices()
{
// set up the mapping of group to input workspace index
this->m_wsIndices.reserve(this->nGroups+1);
size_t nHist_st = static_cast<size_t>(nHist);
for (size_t wi=0;wi<nHist_st;wi++)
{
//wi is the workspace index (of the input)
const int group = groupAtWorkspaceIndex[wi];
if (group < 1) // Not in a group, or invalid group #
continue;

// resize the ws_indices if it is not big enough
if (this->m_wsIndices.size() < static_cast<size_t>(group+1))
{
this->m_wsIndices.resize(group + 1);
}

// Also record a list of workspace indices
this->m_wsIndices[group].push_back(wi);
}

// initialize a vector of the valid group numbers
this->m_validGroups.reserve(nGroups);
size_t totalHistProcess = 0;
for (size_t i = 0; i <this->m_wsIndices.size(); i++)
{
if (!(this->m_wsIndices[i].empty()))
{
this->m_validGroups.push_back(static_cast<int>(i));
totalHistProcess += this->m_wsIndices[i].size();
}
}

return totalHistProcess;
}

} // namespace Algorithm
} // namespace Mantid

0 comments on commit c9d4c0d

Please sign in to comment.