Fix WorkspaceCreation::create(const P &parent) for parent workspaces with varying bins #19247

Merged
merged 13 commits into from Mar 30, 2017

Conversation

Projects
None yet
3 participants
@soininen
Contributor

soininen commented Mar 24, 2017

WorkspaceCreation::create(const P &parent) always copied the X data of the first histogram in the parent workspace to the created one. This would lead to incorrect binning in the created workspace if the parent had varying bins.

To test:

The symptoms can be seen e.g. when using the MonteCarloAbsorption algorithm. In the following script, change the buildDir variable to point to Mantid's build dir and run.

buildDir = <path-to-Mantid-build-directory>
ws = Load(Filename=buildDir + 'ExternalData/Testing/Data/UnitTest/ILL/IN4/084446.nxs')
ws = ConvertUnits(ws, Target='Wavelength', EMode='Direct')
geometry = {
    'Shape': 'Cylinder',
    'Height': 10.0,
    'Radius': 0.5,
    'Center': [0.0, 0.0, 0.0]
}
material = {
    'ChemicalFormula': 'O2',
    'SampleMassDensity': 1
}
SetSample(ws, Geometry=geometry, Material=material)
corrections = MonteCarloAbsorption(ws, 
    NumberOfWavelengthPoints=2, 
    EventsPerPoint=1)

After the unit conversion ws has varying bins for workspace indices > 300. Binning should be the same for corrections after this PR has been applied.

Fixes #19242.

Release notes were updated accordingly.


Reviewer

Please comment on the following (full description):

Code Review
  • Is the code of an acceptable quality?
  • Does the code conform to the coding standards? Is it well structured with small focussed classes/methods/functions?
  • Are there unit/system tests in place? Are the unit tests small and test the a class in isolation?
  • If there are changes in the release notes then do they describe the changes appropriately?
Functional Tests
  • Do changes function as described? Add comments below that describe the tests performed?

  • How do the changes handle unexpected situations, e.g. bad input?

  • Has the relevant documentation been added/updated?

  • Is user-facing documentation written in a user-friendly manner?

  • Has developer documentation been updated if required?

  • Does everything look good? Comment with the ship it emoji but don't merge. A member of @mantidproject/gatekeepers will take care of it.

@SimonHeybrock SimonHeybrock self-assigned this Mar 27, 2017

@SimonHeybrock

This comment has been minimized.

Show comment
Hide comment
@SimonHeybrock

SimonHeybrock Mar 27, 2017

Contributor

I wrote these new factory methods a while ago, and considered exactly this problem. In the end, the reason I had decided against copying all the bins from the parent was the case of creating a smaller workspaces, e.g., with create<Workspace2D>(parent, parent->getNumberHistograms()/2). In that case it could be undefined which X (i.e., which spectra) should be used, so I decided to default.

As this is fixed now, I think we are running exactly into this problem again: Consider, with the changes from this PR:

ws1 = create<Workspace2D>(parent);
ws2 = create<Workspace2D>(parent, parent->getNumberHistograms());`

I think as it stands now, ws2 will have X equal to parent->x(0), whereas ws1 has X identical to parent, for all spectra. Is that a good behavior? Should there be as size check?

I am not trying to say that I do not like copying the bins (I would actually prefer it 👍), but we have to be careful how to define edge cases, and make sure that the differences in behavior are documented.

Contributor

SimonHeybrock commented Mar 27, 2017

I wrote these new factory methods a while ago, and considered exactly this problem. In the end, the reason I had decided against copying all the bins from the parent was the case of creating a smaller workspaces, e.g., with create<Workspace2D>(parent, parent->getNumberHistograms()/2). In that case it could be undefined which X (i.e., which spectra) should be used, so I decided to default.

As this is fixed now, I think we are running exactly into this problem again: Consider, with the changes from this PR:

ws1 = create<Workspace2D>(parent);
ws2 = create<Workspace2D>(parent, parent->getNumberHistograms());`

I think as it stands now, ws2 will have X equal to parent->x(0), whereas ws1 has X identical to parent, for all spectra. Is that a good behavior? Should there be as size check?

I am not trying to say that I do not like copying the bins (I would actually prefer it 👍), but we have to be careful how to define edge cases, and make sure that the differences in behavior are documented.

@SimonHeybrock

This comment has been minimized.

Show comment
Hide comment
@SimonHeybrock

SimonHeybrock Mar 27, 2017

Contributor

Maybe one option would be to remove the variants that accept a size but not a Histogram?

  create<T>(ParentWS, NumSpectra)
  create<T>(ParentWS, IndexInfo)

That would force client code to be explicit, i.e., use create<T>(parent, parent->getNumberHistogram()/2, parent->histogram(0))?

Contributor

SimonHeybrock commented Mar 27, 2017

Maybe one option would be to remove the variants that accept a size but not a Histogram?

  create<T>(ParentWS, NumSpectra)
  create<T>(ParentWS, IndexInfo)

That would force client code to be explicit, i.e., use create<T>(parent, parent->getNumberHistogram()/2, parent->histogram(0))?

@soininen

This comment has been minimized.

Show comment
Hide comment
@soininen

soininen Mar 27, 2017

Contributor

These are very good points. I think forcing the clients to be explicit on how to fill the histogram data is the way to go, i.e. removing the size-only variants. I will try it out and see what will come.

Contributor

soininen commented Mar 27, 2017

These are very good points. I think forcing the clients to be explicit on how to fill the histogram data is the way to go, i.e. removing the size-only variants. I will try it out and see what will come.

soininen added some commits Mar 28, 2017

Remove WorkspaceCreation::create(parent, number) methods.
We now explicitly require clients to provide a template histogram.

Re #19242
@SimonHeybrock

Some comments below. I think most of the changes to the templated methods are not necessary, and you can get the same result with less code (and much fewer changes):

@@ -292,8 +292,7 @@ void DiffractionEventCalibrateDetectors::exec() {
// Get some stuff from the input workspace
// We make a copy of the instrument since we will be moving detectors in
// `inputW` but want to access original positions (etc.) via `detList` below.
- const auto &dummyW = create<EventWorkspace>(*inputW, 1);
- Instrument_const_sptr inst = dummyW->getInstrument();
+ Instrument_const_sptr inst(inputW->getInstrument()->baseInstrument()->clone());

This comment has been minimized.

@SimonHeybrock

SimonHeybrock Mar 28, 2017

Contributor

I don't think this is correct, this way you are losing all position information stored in parameter map and in DetectorInfo. DetectorInfo is currently owned by ExperimentInfo (since Instrument-2.0 is not complete yet), so I think you need to keep creating a workspace.

@SimonHeybrock

SimonHeybrock Mar 28, 2017

Contributor

I don't think this is correct, this way you are losing all position information stored in parameter map and in DetectorInfo. DetectorInfo is currently owned by ExperimentInfo (since Instrument-2.0 is not complete yet), so I think you need to keep creating a workspace.

+ }
+ if (histogram.ptrE()) {
+ histogram.setSharedE(nullptr);
+ }

This comment has been minimized.

@SimonHeybrock

SimonHeybrock Mar 28, 2017

Contributor

Why this change?

@SimonHeybrock

SimonHeybrock Mar 28, 2017

Contributor

Why this change?

This comment has been minimized.

@soininen

soininen Mar 28, 2017

Contributor

I was just being too careful here. I'll remove it.

@soininen

soininen Mar 28, 2017

Contributor

I was just being too careful here. I'll remove it.

+ class = typename std::enable_if<
+ std::is_base_of<API::MatrixWorkspace, P>::value>::type>
+std::unique_ptr<T> create(const P &parent, const IndexArg &indexArg,
+ const HistogramData::Histogram &histogram) {

This comment has been minimized.

@SimonHeybrock

SimonHeybrock Mar 28, 2017

Contributor

Why this new variant? I think it will actually break things since there is code that creates workspaces based on a histogram including data (via the templated version above).

@SimonHeybrock

SimonHeybrock Mar 28, 2017

Contributor

Why this new variant? I think it will actually break things since there is code that creates workspaces based on a histogram including data (via the templated version above).

This comment has been minimized.

@soininen

soininen Mar 28, 2017

Contributor

EventWorkspace::initialize() won't accept histograms with non-NULL y or e data. Check the new unit test test_create_event_from_event(): it will throw at create() when using the original implementation above.
I added this variant to fix the problem. However, as you note, it erroneously strips the data from histogram workspaces.
Maybe this variant (and the other one below) could be merged to the corresponding histArg accepting create():

template <class T, class P, class IndexArg, class HistArg,
          class = typename std::enable_if<
              std::is_base_of<API::MatrixWorkspace, P>::value>::type>
std::unique_ptr<T> create(const P &parent, const IndexArg &indexArg,
                          const HistArg &histArg) {
  std::unique_ptr<T> ws = detail::createUninitialized<T, P>(parent);
  HistogramData::Histogram templateHist(histArg);
  if (std::is_base_of<DataObjects::EventWorkspace, T>::value) {
    templateHisto = detail::stripData(templateHist);
  }  
  ws->initialize(indexArg, templateHist);
  detail::initializeFromParent(parent, *ws);
  return ws;
}

This way we would get rid of the variants and handle EventWorkspaces gracefully. Or am I missing something?

@soininen

soininen Mar 28, 2017

Contributor

EventWorkspace::initialize() won't accept histograms with non-NULL y or e data. Check the new unit test test_create_event_from_event(): it will throw at create() when using the original implementation above.
I added this variant to fix the problem. However, as you note, it erroneously strips the data from histogram workspaces.
Maybe this variant (and the other one below) could be merged to the corresponding histArg accepting create():

template <class T, class P, class IndexArg, class HistArg,
          class = typename std::enable_if<
              std::is_base_of<API::MatrixWorkspace, P>::value>::type>
std::unique_ptr<T> create(const P &parent, const IndexArg &indexArg,
                          const HistArg &histArg) {
  std::unique_ptr<T> ws = detail::createUninitialized<T, P>(parent);
  HistogramData::Histogram templateHist(histArg);
  if (std::is_base_of<DataObjects::EventWorkspace, T>::value) {
    templateHisto = detail::stripData(templateHist);
  }  
  ws->initialize(indexArg, templateHist);
  detail::initializeFromParent(parent, *ws);
  return ws;
}

This way we would get rid of the variants and handle EventWorkspaces gracefully. Or am I missing something?

This comment has been minimized.

@SimonHeybrock

SimonHeybrock Mar 28, 2017

Contributor

I think client code should simply use a Histogram with null y and e data in this case. Removing it silently might actually hide issues.

@SimonHeybrock

SimonHeybrock Mar 28, 2017

Contributor

I think client code should simply use a Histogram with null y and e data in this case. Removing it silently might actually hide issues.

This comment has been minimized.

@SimonHeybrock

SimonHeybrock Mar 28, 2017

Contributor

@soininen Note that it is trivial to do that, since client code can simply use

create<EventWorkspace>(parent, BinEdges{1,2,3,4});

The histogram with null Y and E is created automatically from BinEdges, so there is no burden on the client code.

@SimonHeybrock

SimonHeybrock Mar 28, 2017

Contributor

@soininen Note that it is trivial to do that, since client code can simply use

create<EventWorkspace>(parent, BinEdges{1,2,3,4});

The histogram with null Y and E is created automatically from BinEdges, so there is no burden on the client code.

This comment has been minimized.

@soininen

soininen Mar 28, 2017

Contributor

In the code above, the e and y data are removed only for EventWorkspaces, which don't use this data anyhow. However, if it is OK that

create<WorkspaceType>(parent, 2, parent.histogram(0));

works only if WorkspaceType != EventWorkspace, then I think I can completely revert the variants. Personally, I don't care either way, it is just that somebody might be surprised when EventWorkspaces suddenly cause runtime errors with something that works fine with Workspace2D. Let me know what you think.

@soininen

soininen Mar 28, 2017

Contributor

In the code above, the e and y data are removed only for EventWorkspaces, which don't use this data anyhow. However, if it is OK that

create<WorkspaceType>(parent, 2, parent.histogram(0));

works only if WorkspaceType != EventWorkspace, then I think I can completely revert the variants. Personally, I don't care either way, it is just that somebody might be surprised when EventWorkspaces suddenly cause runtime errors with something that works fine with Workspace2D. Let me know what you think.

This comment has been minimized.

@SimonHeybrock

SimonHeybrock Mar 28, 2017

Contributor

I prefer runtime errors over behavior that silently removes data -- correctness is the most important aspect I think!

@SimonHeybrock

SimonHeybrock Mar 28, 2017

Contributor

I prefer runtime errors over behavior that silently removes data -- correctness is the most important aspect I think!

+ !std::is_base_of<API::MatrixWorkspace, IndexArg>::value>::type * =
+ nullptr>
+std::unique_ptr<T> create(const IndexArg &indexArg,
+ const HistogramData::Histogram &histogram) {

This comment has been minimized.

@SimonHeybrock

SimonHeybrock Mar 28, 2017

Contributor

Why this new variant? (see above)

@SimonHeybrock

SimonHeybrock Mar 28, 2017

Contributor

Why this new variant? (see above)

- detail::stripData(parent.histogram(0)));
-}
+ const auto numHistograms = parent.getNumberHistograms();
+ if (parent.isCommonBins()) {

This comment has been minimized.

@SimonHeybrock

SimonHeybrock Mar 28, 2017

Contributor

I think the isCommonBins method is a bit broken, I would not use it. Simply remove this if and set up sharing X below, that should be both safer and faster.

@SimonHeybrock

SimonHeybrock Mar 28, 2017

Contributor

I think the isCommonBins method is a bit broken, I would not use it. Simply remove this if and set up sharing X below, that should be both safer and faster.

This comment has been minimized.

@soininen

soininen Mar 28, 2017

Contributor

Good idea! isCommonBins is indeed not the most reliable of tests.

@soininen

soininen Mar 28, 2017

Contributor

Good idea! isCommonBins is indeed not the most reliable of tests.

+ const auto XLength = parent.isHistogramData() ? YLength + 1 : YLength;
+ ws->initialize(numHistograms, XLength, YLength);
+ for (size_t i = 0; i < numHistograms; ++i) {
+ ws->mutableX(i) = parent.x(i);

This comment has been minimized.

@SimonHeybrock

SimonHeybrock Mar 28, 2017

Contributor

ws->setSharedX(i, parent.sharedX(i)); sets up sharing and is much faster.

@SimonHeybrock

SimonHeybrock Mar 28, 2017

Contributor

ws->setSharedX(i, parent.sharedX(i)); sets up sharing and is much faster.

+ std::unique_ptr<T> ws = detail::createUninitialized<T, P>(parent);
+ const auto YLength = parent.blocksize();
+ const auto XLength = parent.isHistogramData() ? YLength + 1 : YLength;
+ ws->initialize(numHistograms, XLength, YLength);

This comment has been minimized.

@SimonHeybrock

SimonHeybrock Mar 28, 2017

Contributor

I think you can simply call the create(parent, numHistograms, detail::stripData(parent.histogram(0))); method here, and then set the correct X afterwards, the result is identical. As a result, you do not need to split out the createUninitialized method from the original create method. There should be no performance penalty to this.

@SimonHeybrock

SimonHeybrock Mar 28, 2017

Contributor

I think you can simply call the create(parent, numHistograms, detail::stripData(parent.histogram(0))); method here, and then set the correct X afterwards, the result is identical. As a result, you do not need to split out the createUninitialized method from the original create method. There should be no performance penalty to this.

soininen added some commits Mar 28, 2017

Remove two create() variants introduced earlier.
We now strip y and e data for EventWorkspaces in the original methods.

Re #19242
Remove special treatment of EventWorkspaces in workspaceCreation.
Clients should explicitly use BinEdges when providing histArg.

Re #19242
@@ -130,6 +128,13 @@ template <>
MANTID_DATAOBJECTS_DLL std::unique_ptr<API::HistoWorkspace>
createConcreteHelper();
+template <class T, class P, class = typename std::enable_if<std::is_base_of<
+ API::MatrixWorkspace, P>::value>::type>
+std::unique_ptr<T> createUninitialized(const P &parent) {

This comment has been minimized.

@SimonHeybrock

SimonHeybrock Mar 28, 2017

Contributor

Now unused?

@SimonHeybrock

SimonHeybrock Mar 28, 2017

Contributor

Now unused?

@@ -57,7 +57,8 @@ As a consequence of these changes, :ref:`CopyInstrumentParmeters <algm-CopyInstr
Bugs
----
-- We have fixed a bug where Mantid could crash when deleteing a large number of workspaces.
+- We have fixed a bug where Mantid could crash when deleting a large number of workspaces.
+- A bug was fixed where sometimes a workspace created using another one as a template would have incorrect X data. This happened only when the template workspace had varying bins between the histograms. In these cases the binning was not copied correctly from the template. The bug was discovered in :ref:`algm-MonteCarloAbsorption` but other algorithms may have been affected as well.

This comment has been minimized.

@SimonHeybrock

SimonHeybrock Mar 28, 2017

Contributor

These workspace factory methods are not exposed to Python, so I am not sure it makes sense to be so specific in the release notes. Maybe only report the fix for the affected algorithm?

@SimonHeybrock

SimonHeybrock Mar 28, 2017

Contributor

These workspace factory methods are not exposed to Python, so I am not sure it makes sense to be so specific in the release notes. Maybe only report the fix for the affected algorithm?

This comment has been minimized.

@soininen

soininen Mar 28, 2017

Contributor

I wasn't sure if this affected other algorithms/user scripts, so I wrote the entire story here. But perhaps you are right, this is not so big deal. I'll revise the release notes.

@soininen

soininen Mar 28, 2017

Contributor

I wasn't sure if this affected other algorithms/user scripts, so I wrote the entire story here. But perhaps you are right, this is not so big deal. I'll revise the release notes.

- const auto &dummyW = create<EventWorkspace>(*inputW, 1);
- Instrument_const_sptr inst = dummyW->getInstrument();
+ Instrument_const_sptr inst(
+ inputW->getInstrument()->baseInstrument()->clone());

This comment has been minimized.

@SimonHeybrock

SimonHeybrock Mar 30, 2017

Contributor

@soininen I had previously commented on this, but I think it was lost in clang-format:

I don't think this is correct, this way you are losing all position information stored in parameter map and in DetectorInfo. DetectorInfo is currently owned by ExperimentInfo (since Instrument-2.0 is not complete yet), so I think you need to keep creating a workspace.

@SimonHeybrock

SimonHeybrock Mar 30, 2017

Contributor

@soininen I had previously commented on this, but I think it was lost in clang-format:

I don't think this is correct, this way you are losing all position information stored in parameter map and in DetectorInfo. DetectorInfo is currently owned by ExperimentInfo (since Instrument-2.0 is not complete yet), so I think you need to keep creating a workspace.

This comment has been minimized.

@soininen

soininen Mar 30, 2017

Contributor

I see - I'll fix it. Sorry for missing your comment!

@soininen

soininen Mar 30, 2017

Contributor

I see - I'll fix it. Sorry for missing your comment!

@SimonHeybrock

This comment has been minimized.

Show comment
Hide comment
@SimonHeybrock

SimonHeybrock Mar 30, 2017

Contributor

Latest changes look good. Overall I think this is a nice improvement of these factory methods 👍
:shipit: once builds pass.

Contributor

SimonHeybrock commented Mar 30, 2017

Latest changes look good. Overall I think this is a nice improvement of these factory methods 👍
:shipit: once builds pass.

@peterfpeterson peterfpeterson merged commit 9a16638 into master Mar 30, 2017

9 checks passed

ClangFormat Jenkins build pull_requests-clang-format 12453 has succeeded
Details
Doxygen Jenkins build pull_requests-doxygen 11778 has succeeded
Details
Flake8 Jenkins build pull_requests-flake8 3102 has succeeded
Details
OSX Jenkins build pull_requests-osx 12967 has succeeded
Details
RHEL7 + System Tests Jenkins build pull_requests-rhel7 12849 has succeeded
Details
Ubuntu + Doc Tests Jenkins build pull_requests-ubuntu 13444 has succeeded
Details
Ubuntu Python 3 Jenkins build pull_requests-ubuntu-python3 985 has succeeded
Details
Windows Jenkins build pull_requests-win7 13710 has succeeded
Details
cppcheck Jenkins build pull_requests-cppcheck 13437 has succeeded
Details

@peterfpeterson peterfpeterson deleted the 19242_WorkspaceCreation_fix branch Mar 30, 2017

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment