Mantid
Loading...
Searching...
No Matches
BinaryOperation.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
8#include "MantidAPI/Axis.h"
21#include "MantidHistogramData/Histogram.h"
22#include "MantidKernel/Timer.h"
23#include "MantidKernel/Unit.h"
24#include <memory>
25
26using namespace Mantid::Geometry;
27using namespace Mantid::API;
28using namespace Mantid::Kernel;
29using namespace Mantid::DataObjects;
30using namespace Mantid::HistogramData;
31using std::size_t;
32
33namespace Mantid::Algorithms {
40 "The name of the input workspace on the left hand side of the operation");
42 "The name of the input workspace on the right hand side of "
43 "the operation");
45 "The name to call the output workspace");
46 declareProperty(std::make_unique<PropertyWithValue<bool>>("AllowDifferentNumberSpectra", false, Direction::Input),
47 "Are workspaces with different number of spectra allowed? "
48 "For example, the LHSWorkspace might have one spectrum per detector, "
49 "but the RHSWorkspace could have its spectra averaged per bank. If true, "
50 "then matching between the LHS and RHS spectra is performed (all "
51 "detectors "
52 "in a LHS spectrum have to be in the corresponding RHS) in order to "
53 "apply the RHS spectrum to the LHS.");
54
55 declareProperty(std::make_unique<PropertyWithValue<bool>>("ClearRHSWorkspace", false, Direction::Input),
56 "For EventWorkspaces only. This will clear out event lists "
57 "from the RHS workspace as the binary operation is applied. "
58 "This can prevent excessive memory use, e.g. when subtracting "
59 "an EventWorkspace from another: memory use will be approximately "
60 "constant instead of increasing by 50%. At completion, the RHS workspace "
61 "will be empty.");
62}
63
69 // Is the LHS operand a single number?
70 WorkspaceSingleValue_const_sptr lhs_singleVal = std::dynamic_pointer_cast<const WorkspaceSingleValue>(m_lhs);
71 WorkspaceSingleValue_const_sptr rhs_singleVal = std::dynamic_pointer_cast<const WorkspaceSingleValue>(m_rhs);
72
73 if (lhs_singleVal) {
74 MatrixWorkspace_sptr out = getProperty("OutputWorkspace");
75 if (this->name() == "Divide" && !bool(rhs_singleVal)) {
76 // x / workspace = Power(workspace, -1) * x
77 // workspace ^ -1
78 auto pow = createChildAlgorithm("Power", 0.0, 0.5, true);
79 pow->setProperty("InputWorkspace", std::const_pointer_cast<MatrixWorkspace>(m_rhs));
80 pow->setProperty("Exponent", -1.0);
81 pow->setProperty("OutputWorkspace", out);
82 pow->executeAsChildAlg();
83 out = pow->getProperty("OutputWorkspace");
84
85 // Multiply by x
86 auto mult = createChildAlgorithm("Multiply", 0.5, 1.0, true);
87 mult->setProperty(inputPropName1(), out); //(workspace^-1)
88 mult->setProperty(inputPropName2(),
89 std::const_pointer_cast<MatrixWorkspace>(m_lhs)); // (1.0) or other number
90 mult->setProperty(outputPropName(), out);
91 mult->executeAsChildAlg();
92 out = mult->getProperty("OutputWorkspace");
93 setProperty("OutputWorkspace", out);
94 return true;
95 } else if (this->name() == "Minus") {
96 // x - workspace = x + (workspace * -1)
97 MatrixWorkspace_sptr minusOne = create<WorkspaceSingleValue>(1, Points(1));
98 minusOne->dataY(0)[0] = -1.0;
99 minusOne->dataE(0)[0] = 0.0;
100
101 // workspace * -1
102 auto mult = createChildAlgorithm("Multiply", 0.0, 0.5, true);
103 mult->setProperty(inputPropName1(), std::const_pointer_cast<MatrixWorkspace>(m_rhs));
104 mult->setProperty(inputPropName2(), minusOne);
105 mult->setProperty("OutputWorkspace", out);
106 mult->executeAsChildAlg();
107 out = mult->getProperty("OutputWorkspace");
108
109 // Multiply by x
110 auto plus = createChildAlgorithm("Plus", 0.5, 1.0, true);
111 plus->setProperty(inputPropName1(), out); //(workspace^-1)
112 plus->setProperty(inputPropName2(),
113 std::const_pointer_cast<MatrixWorkspace>(m_lhs)); // (1.0) or other number
114 plus->setProperty(outputPropName(), out);
115 plus->executeAsChildAlg();
116 out = plus->getProperty("OutputWorkspace");
117 setProperty("OutputWorkspace", out);
118 return true;
119 }
120
121 } // lhs_singleVal
122
123 // Process normally
124 return false;
125}
126
132 // get input workspace, dynamic cast not needed
135 m_AllowDifferentNumberSpectra = getProperty("AllowDifferentNumberSpectra");
136
137 m_lhsBlocksize = m_lhs->blocksize();
138 m_rhsBlocksize = m_rhs->blocksize();
139
140 // Special handling for 1-WS and 1/WS.
141 if (this->handleSpecialDivideMinus())
142 return;
143
144 // Cast to EventWorkspace pointers
145 m_elhs = std::dynamic_pointer_cast<const EventWorkspace>(m_lhs);
146 m_erhs = std::dynamic_pointer_cast<const EventWorkspace>(m_rhs);
147
148 // We can clear the RHS workspace if it is an event,
149 // and we are not doing mismatched spectra (in which case you might clear it
150 // too soon!)
151 // and lhs is not rhs.
152 // and out is not rhs.
153 m_ClearRHSWorkspace = getProperty("ClearRHSWorkspace");
155 if (m_AllowDifferentNumberSpectra || (!m_erhs) || (m_rhs == m_lhs) || (m_out == m_rhs)) {
156 // std::cout << "m_ClearRHSWorkspace = false\n";
157 m_ClearRHSWorkspace = false;
158 }
159 }
160
161 // Make a check of what will be needed to setup the workspaces, based on the
162 // input types.
163 this->checkRequirements();
164
165 if (m_flipSides) {
166 // Flip the workspaces left and right
167 std::swap(m_lhs, m_rhs);
168 std::swap(m_elhs, m_erhs);
169 std::swap(m_lhsBlocksize, m_rhsBlocksize);
170 }
171
172 // Check that the input workspaces are compatible
174 std::ostringstream ostr;
175 ostr << "The two workspaces are not compatible for algorithm " << this->name();
176 g_log.error() << ostr.str() << '\n';
177 throw std::invalid_argument(ostr.str());
178 }
179
180 // Get the output workspace
182 if (m_elhs) {
183 m_eout = std::dynamic_pointer_cast<EventWorkspace>(m_out);
184 }
185
186 // Is the output going to be an EventWorkspace?
188 // The output WILL be EventWorkspace (this implies lhs is EW or rhs is EW +
189 // it gets flipped)
190 if (!m_elhs)
191 throw std::runtime_error("BinaryOperation:: the output was set to be an "
192 "EventWorkspace (m_keepEventWorkspace == true), "
193 "but the lhs is not an EventWorkspace. There "
194 "must be a mistake in the algorithm. Contact "
195 "the developers.");
196
197 if (m_out == m_lhs) {
198 // Will be modifying the EventWorkspace in-place on the lhs. Good.
199 if (!m_eout)
200 throw std::runtime_error("BinaryOperation:: the output was set to be lhs, and to be an "
201 "EventWorkspace (m_keepEventWorkspace == true), but the output is "
202 "not an EventWorkspace. There must be a mistake in the algorithm. "
203 "Contact the developers.");
204 } else {
205 // You HAVE to copy the data from lhs to to the output!
206 m_out = m_lhs->clone();
207 // Make sure m_eout still points to the same as m_out;
208 m_eout = std::dynamic_pointer_cast<EventWorkspace>(m_out);
209 }
210
211 // Always clear the MRUs.
212 m_eout->clearMRU();
213 if (m_elhs)
214 m_elhs->clearMRU();
215 if (m_erhs)
216 m_erhs->clearMRU();
217
218 } else {
219 // ---- Output will be WS2D -------
220
221 // We need to create a new workspace for the output if:
222 // (a) the output workspace hasn't been set to one of the input ones, or
223 // (b) it has been, but it's not the correct dimensions
224 if ((m_out != m_lhs && m_out != m_rhs) || (m_out == m_rhs && (m_lhs->size() > m_rhs->size()))) {
225 // if the input workspace are specialworkspace2d, then we need to ensure
226 // the map is set
227 auto specialLHS = dynamic_cast<const SpecialWorkspace2D *>(m_lhs.get());
228 auto specialRHS = dynamic_cast<const SpecialWorkspace2D *>(m_rhs.get());
229 if (specialLHS && specialRHS) {
230 m_out = create<SpecialWorkspace2D>(*specialLHS);
231 } else {
232 m_out = create<HistoWorkspace>(*m_lhs);
233 }
234 }
235 }
236
237 // only overridden for some operations (plus and minus at the time of writing)
238 operateOnRun(m_lhs->run(), m_rhs->run(), m_out->mutableRun());
239
240 // Initialise the progress reporting object
241 m_progress = std::make_unique<Progress>(this, 0.0, 1.0, m_lhs->getNumberHistograms());
242
243 // There are now 4 possible scenarios, shown schematically here:
244 // xxx x xxx xxx xxx xxx xxx x
245 // xxx , xxx xxx , xxx , xxx x
246 // xxx , xxx xxx xxx xxx x
247 // So work out which one we have and call the appropriate function
248
249 // Single value workspace on the right : if it is an EventWorkspace with 1
250 // spectrum, 1 bin, it is treated as a scalar
251 if ((m_rhs->size() == 1) && !m_do2D_even_for_SingleColumn_on_rhs) {
252 doSingleValue(); // m_lhs,m_rhs,m_out
253 } else if (m_rhs->getNumberHistograms() == 1) // Single spectrum on rhs
254 {
256 }
257 // Single column on rhs; if the RHS is an event workspace with one bin, it is
258 // treated as a scalar.
261 } else // The two are both 2D and should be the same size (except if LHS is an
262 // event workspace)
263 {
264 bool mismatchedSpectra =
265 (m_AllowDifferentNumberSpectra && (m_rhs->getNumberHistograms() != m_lhs->getNumberHistograms()));
266 do2D(mismatchedSpectra);
267 }
268
270
271 // Assign the result to the output workspace property
273}
274
282 UNUSED_ARG(lhs);
284 // This should never happen
285 throw Exception::NotImplementedError("BinaryOperation::execEvent() is not implemented for this operation.");
286}
287
297 Unit_const_sptr lhs_unit;
298 Unit_const_sptr rhs_unit;
299 if (lhs->axes() && rhs->axes()) // If one of these is a WorkspaceSingleValue
300 // then we don't want to check units match
301 {
302 lhs_unit = lhs->getAxis(0)->unit();
303 rhs_unit = rhs->getAxis(0)->unit();
304 }
305
306 const std::string lhs_unitID = (lhs_unit ? lhs_unit->unitID() : "");
307 const std::string rhs_unitID = (rhs_unit ? rhs_unit->unitID() : "");
308
309 // Check the workspaces have the same units and distribution flag
310 if (lhs_unitID != rhs_unitID && m_lhsBlocksize > 1 && m_rhsBlocksize > 1) {
311 g_log.error("The two workspace are not compatible because they have "
312 "different units on the X axis.");
313 return false;
314 }
315
316 // Check the size compatibility
317 const std::string checkSizeCompatibilityResult = checkSizeCompatibility(lhs, rhs);
318 if (!checkSizeCompatibilityResult.empty()) {
319 throw std::invalid_argument(checkSizeCompatibilityResult);
320 }
321
322 return true;
323}
324
334 UNUSED_ARG(lhs);
336 return false;
337}
338
351 const size_t lhsSize = lhs->size();
352 const size_t rhsSize = rhs->size();
353 // A SingleValueWorkspace on the right matches anything
354 if (rhsSize == 1)
355 return "";
356 // The lhs must not be smaller than the rhs
357 if (lhsSize < rhsSize)
358 return "Left hand side smaller than right hand side.";
359
360 // Did checkRequirements() tell us that the X histogram size did not matter?
361 if (!m_matchXSize) {
362 // If so, only the vertical # needs to match
363
364 if (lhs->getNumberHistograms() == rhs->getNumberHistograms()) {
365 return "";
366 } else {
367 return "Number of histograms not identical.";
368 }
369 }
370 // Otherwise they must match both ways, or horizontally or vertically with the
371 // other rhs dimension=1
372 if (m_rhsBlocksize == 1 && lhs->getNumberHistograms() == rhs->getNumberHistograms())
373 return "";
374 // Past this point, we require the X arrays to match. Note this only checks
375 // the first spectrum
376 if (!WorkspaceHelpers::matchingBins(*lhs, *rhs, true)) {
377 return "X arrays must match when performing this operation on a 2D "
378 "workspaces.";
379 }
380
381 const size_t rhsSpec = rhs->getNumberHistograms();
382
384 if (rhsSpec == 1 || lhs->getNumberHistograms() == rhsSpec) {
385 return "";
386 } else {
387 // can't be more specific as if this is reached both failed and only one
388 // or both are needed
389 return "Left and right sides should contain the same amount of spectra "
390 "or the right side should contain only one spectra.";
391 }
392 } else {
393 // blocksize check failed, but still check the number of spectra to see if
394 // that was wrong too
395 if (rhsSpec == 1 || lhs->getNumberHistograms() == rhsSpec) {
396 return "Number of y values not equal on left and right sides.";
397 } else {
398 // can't be more specific as if this is reached both failed and only one
399 // or both are needed
400 return "Number of y values not equal on left and right sides and the "
401 "right side contained neither only one spectra or the same amount "
402 "of spectra as the left.";
403 }
404 }
405}
406
419bool BinaryOperation::propagateSpectraMask(const SpectrumInfo &lhsSpectrumInfo, const SpectrumInfo &rhsSpectrumInfo,
420 const int64_t index, MatrixWorkspace &out, SpectrumInfo &outSpectrumInfo) {
421 bool continueOp(true);
422
423 if ((lhsSpectrumInfo.hasDetectors(index) && lhsSpectrumInfo.isMasked(index)) ||
424 (rhsSpectrumInfo.hasDetectors(index) && rhsSpectrumInfo.isMasked(index))) {
425 continueOp = false;
427 PARALLEL_CRITICAL(setMasked) { outSpectrumInfo.setMasked(index, true); }
428 }
429 return continueOp;
430}
431
438 // Don't propate masking from the rhs here - it would be decidedly odd if the
439 // single value was masked
440
441 // Pull out the single value and its error
442 const double rhsY = m_rhs->y(0)[0];
443 const double rhsE = m_rhs->e(0)[0];
444
445 // Now loop over the spectra of the left hand side calling the virtual
446 // function
447 const int64_t numHists = m_lhs->getNumberHistograms();
448
449 if (m_eout) {
450 // ---- The output is an EventWorkspace ------
452 for (int64_t i = 0; i < numHists; ++i) {
454 m_out->setSharedX(i, m_lhs->sharedX(i));
455 performEventBinaryOperation(m_eout->getSpectrum(i), rhsY, rhsE);
456 m_progress->report(this->name());
458 }
460 } else {
461 // ---- Histogram Output -----
463 for (int64_t i = 0; i < numHists; ++i) {
465 m_out->setSharedX(i, m_lhs->sharedX(i));
466 // Get reference to output vectors here to break any sharing outside the
467 // function call below
468 // where the order of argument evaluation is not guaranteed (if it's L->R
469 // there would be a data race)
470 HistogramData::HistogramY &outY = m_out->mutableY(i);
471 HistogramData::HistogramE &outE = m_out->mutableE(i);
472 performBinaryOperation(m_lhs->histogram(i), rhsY, rhsE, outY, outE);
473 m_progress->report(this->name());
475 }
477 }
478}
479
485 // Don't propate masking from the m_rhs here - it would be decidedly odd if
486 // the single bin was masked
487
488 // Now loop over the spectra of the left hand side pulling m_out the single
489 // value from each m_rhs 'spectrum'
490 // and then calling the virtual function
491 const int64_t numHists = m_lhs->getNumberHistograms();
492 auto &outSpectrumInfo = m_out->mutableSpectrumInfo();
493 auto &lhsSpectrumInfo = m_lhs->spectrumInfo();
494 auto &rhsSpectrumInfo = m_rhs->spectrumInfo();
495 if (m_eout) {
496 // ---- The output is an EventWorkspace ------
498 for (int64_t i = 0; i < numHists; ++i) {
500 const double rhsY = m_rhs->y(i)[0];
501 const double rhsE = m_rhs->e(i)[0];
502 if (propagateSpectraMask(lhsSpectrumInfo, rhsSpectrumInfo, i, *m_out, outSpectrumInfo)) {
503 performEventBinaryOperation(m_eout->getSpectrum(i), rhsY, rhsE);
504 }
505 m_progress->report(this->name());
507 }
509 } else {
510 // ---- Histogram Output -----
512 for (int64_t i = 0; i < numHists; ++i) {
514 const double rhsY = m_rhs->y(i)[0];
515 const double rhsE = m_rhs->e(i)[0];
516
517 m_out->setSharedX(i, m_lhs->sharedX(i));
518 if (propagateSpectraMask(lhsSpectrumInfo, rhsSpectrumInfo, i, *m_out, outSpectrumInfo)) {
519 // Get reference to output vectors here to break any sharing outside the
520 // function call below
521 // where the order of argument evaluation is not guaranteed (if it's
522 // L->R there would be a data race)
523 HistogramData::HistogramY &outY = m_out->mutableY(i);
524 HistogramData::HistogramE &outE = m_out->mutableE(i);
525 performBinaryOperation(m_lhs->histogram(i), rhsY, rhsE, outY, outE);
526 }
527 m_progress->report(this->name());
529 }
531 }
532}
533
538
539 // Propagate any masking first or it could mess up the numbers
540 // TODO: Check if this works for event workspaces...
542
543 if (m_eout) {
544 // ----------- The output is an EventWorkspace -------------
545
546 if (m_erhs) {
547 // -------- The rhs is ALSO an EventWorkspace --------
548
549 // Pull out the single eventList on the right
550 const EventList &rhs_spectrum = m_erhs->getSpectrum(0);
551
552 // Now loop over the spectra of the left hand side calling the virtual
553 // function
554 const int64_t numHists = m_lhs->getNumberHistograms();
556 for (int64_t i = 0; i < numHists; ++i) {
558 // Perform the operation on the event list on the output (== lhs)
559 performEventBinaryOperation(m_eout->getSpectrum(i), rhs_spectrum);
560 m_progress->report(this->name());
562 }
564 } else {
565 // -------- The rhs is a histogram ---------
566 // Pull m_out the m_rhs spectrum
567 const MantidVec &rhsX = m_rhs->readX(0);
568 const MantidVec &rhsY = m_rhs->readY(0);
569 const MantidVec &rhsE = m_rhs->readE(0);
570
571 // Now loop over the spectra of the left hand side calling the virtual
572 // function
573 const int64_t numHists = m_lhs->getNumberHistograms();
574
576 for (int64_t i = 0; i < numHists; ++i) {
578 // Perform the operation on the event list on the output (== lhs)
579 performEventBinaryOperation(m_eout->getSpectrum(i), rhsX, rhsY, rhsE);
580 m_progress->report(this->name());
582 }
584 }
585
586 } else {
587 // -------- The output is a histogram ----------
588 // (inputs can be EventWorkspaces, but their histogram representation
589 // will be used instead)
590
591 // Pull m_out the m_rhs spectrum
592 const auto rhs = m_rhs->histogram(0);
593
594 // Now loop over the spectra of the left hand side calling the virtual
595 // function
596 const int64_t numHists = m_lhs->getNumberHistograms();
597
599 for (int64_t i = 0; i < numHists; ++i) {
601 m_out->setSharedX(i, m_lhs->sharedX(i));
602 // Get reference to output vectors here to break any sharing outside the
603 // function call below
604 // where the order of argument evaluation is not guaranteed (if it's L->R
605 // there would be a data race)
606 HistogramData::HistogramY &outY = m_out->mutableY(i);
607 HistogramData::HistogramE &outE = m_out->mutableE(i);
608 performBinaryOperation(m_lhs->histogram(i), rhs, outY, outE);
609 m_progress->report(this->name());
611 }
613 }
614}
615
624void BinaryOperation::do2D(bool mismatchedSpectra) {
626 if (mismatchedSpectra) {
628 }
629
630 // Propagate any masking first or it could mess up the numbers
631 // TODO: Check if this works for event workspaces...
633
634 auto &outSpectrumInfo = m_out->mutableSpectrumInfo();
635 auto &lhsSpectrumInfo = m_lhs->spectrumInfo();
636 auto &rhsSpectrumInfo = m_rhs->spectrumInfo();
637
638 if (m_eout) {
639 // ----------- The output is an EventWorkspace -------------
640
642 // ------------ The rhs is ALSO an EventWorkspace ---------------
643 // Now loop over the spectra of each one calling the virtual function
644 const int64_t numHists = m_lhs->getNumberHistograms();
646 for (int64_t i = 0; i < numHists; ++i) {
648 m_progress->report(this->name());
649
650 int64_t rhs_wi = i;
651 if (mismatchedSpectra && table) {
652 rhs_wi = (*table)[i];
653 if (rhs_wi < 0)
654 continue;
655 } else {
656 // Check for masking except when mismatched sizes
657 if (!propagateSpectraMask(lhsSpectrumInfo, rhsSpectrumInfo, i, *m_out, outSpectrumInfo))
658 continue;
659 }
660 // Reach here? Do the division
661 // Perform the operation on the event list on the output (== lhs)
662 performEventBinaryOperation(m_eout->getSpectrum(i), m_erhs->getSpectrum(rhs_wi));
663
664 // Free up memory on the RHS if that is possible
666 const_cast<EventList &>(m_erhs->getSpectrum(rhs_wi)).clear();
668 }
670 } else {
671 // -------- The rhs is a histogram, or we want to use the histogram
672 // representation of it ---------
673
674 // Now loop over the spectra of each one calling the virtual function
675 const int64_t numHists = m_lhs->getNumberHistograms();
676
678 for (int64_t i = 0; i < numHists; ++i) {
680 m_progress->report(this->name());
681 int64_t rhs_wi = i;
682 if (mismatchedSpectra && table) {
683 rhs_wi = (*table)[i];
684 if (rhs_wi < 0)
685 continue;
686 } else {
687 // Check for masking except when mismatched sizes
688 if (!propagateSpectraMask(lhsSpectrumInfo, rhsSpectrumInfo, i, *m_out, outSpectrumInfo))
689 continue;
690 }
691
692 // Reach here? Do the division
693 performEventBinaryOperation(m_eout->getSpectrum(i), m_rhs->readX(rhs_wi), m_rhs->readY(rhs_wi),
694 m_rhs->readE(rhs_wi));
695
696 // Free up memory on the RHS if that is possible
698 const_cast<EventList &>(m_erhs->getSpectrum(rhs_wi)).clear();
699
701 }
703 }
704
705 } else {
706 // -------- The output is a histogram ----------
707 // (inputs can be EventWorkspaces, but their histogram representation
708 // will be used instead)
709
710 // Now loop over the spectra of each one calling the virtual function
711 const int64_t numHists = m_lhs->getNumberHistograms();
712
714 for (int64_t i = 0; i < numHists; ++i) {
716 m_progress->report(this->name());
717 m_out->setSharedX(i, m_lhs->sharedX(i));
718 int64_t rhs_wi = i;
719 if (mismatchedSpectra && table) {
720 rhs_wi = (*table)[i];
721 if (rhs_wi < 0)
722 continue;
723 } else {
724 // Check for masking except when mismatched sizes
725 if (!propagateSpectraMask(lhsSpectrumInfo, rhsSpectrumInfo, i, *m_out, outSpectrumInfo))
726 continue;
727 }
728 // Reach here? Do the division
729 // Get reference to output vectors here to break any sharing outside the
730 // function call below
731 // where the order of argument evaluation is not guaranteed (if it's L->R
732 // there would be a data race)
733 HistogramData::HistogramY &outY = m_out->mutableY(i);
734 HistogramData::HistogramE &outE = m_out->mutableE(i);
735 performBinaryOperation(m_lhs->histogram(i), m_rhs->histogram(rhs_wi), outY, outE);
736
737 // Free up memory on the RHS if that is possible
739 const_cast<EventList &>(m_erhs->getSpectrum(rhs_wi)).clear();
740
742 }
744 }
745 // Make sure we don't use outdated MRU
747 m_erhs->clearMRU();
748}
749
757 const API::MatrixWorkspace_sptr &out) {
758 const int64_t outHists = out->getNumberHistograms();
759 const int64_t rhsHists = rhs->getNumberHistograms();
760 for (int64_t i = 0; i < outHists; ++i) {
761 // Copy over masks from the rhs, if any exist.
762 // If rhs is single spectrum, copy masks from that to all spectra in the
763 // output.
764 if (rhs->hasMaskedBins((rhsHists == 1) ? 0 : i)) {
765 const MatrixWorkspace::MaskList &masks = rhs->maskedBins((rhsHists == 1) ? 0 : i);
766 MatrixWorkspace::MaskList::const_iterator it;
767 for (it = masks.begin(); it != masks.end(); ++it) {
768 out->flagMasked(i, it->first, it->second);
769 }
770 }
771 }
772}
773
774// ------- Default implementations of Event binary operations --------
775
785 UNUSED_ARG(lhs);
787 throw Exception::NotImplementedError("BinaryOperation::performEventBinaryOperation() not implemented.");
788}
789
800 const MantidVec &rhsY, const MantidVec &rhsE) {
801 UNUSED_ARG(lhs);
802 UNUSED_ARG(rhsX);
803 UNUSED_ARG(rhsY);
804 UNUSED_ARG(rhsE);
805 throw Exception::NotImplementedError("BinaryOperation::performEventBinaryOperation() not implemented.");
806}
807
816void BinaryOperation::performEventBinaryOperation(DataObjects::EventList &lhs, const double &rhsY, const double &rhsE) {
817 UNUSED_ARG(lhs);
818 UNUSED_ARG(rhsY);
819 UNUSED_ARG(rhsE);
820 throw Exception::NotImplementedError("BinaryOperation::performEventBinaryOperation() not implemented.");
821}
822
829 // An event workspace?
830 EventWorkspace_const_sptr ews = std::dynamic_pointer_cast<const EventWorkspace>(ws);
831 if (ews)
832 return eEventList;
833
834 // If the workspace has no axes, then it is a WorkspaceSingleValue
835 if (!ws->axes())
836 return eNumber;
837
838 // TODO: Check if it is a single-colum one, then
839 // return Number;
840
841 // Otherwise, we take it as a histogram (workspace 2D)
842 return eHistogram;
843}
844
852
853 // In general, the X sizes have to match.
854 // (only for some EventWorkspace cases does it NOT matter...)
855 m_matchXSize = true;
856
857 // Operations are not always commutative! Don't flip sides.
858 m_flipSides = false;
859
860 // And in general, EventWorkspaces get turned to Workspace2D
861 m_keepEventWorkspace = false;
862
863 // This will be set to true for Divide/Multiply
865}
866
881 // An addition table is a list of pairs:
882 // First int = workspace index in the EW being added
883 // Second int = workspace index to which it will be added in the OUTPUT EW.
884 // -1 if it should add a new entry at the end.
885 auto table = std::make_shared<BinaryOperationTable>();
886
887 auto rhs_nhist = static_cast<int>(rhs->getNumberHistograms());
888 auto lhs_nhist = static_cast<int>(lhs->getNumberHistograms());
889
890 // Initialize the table; filled with -1 meaning no match
891 table->resize(lhs_nhist, -1);
892
893 const detid2index_map rhs_det_to_wi = rhs->getDetectorIDToWorkspaceIndexMap();
894
896 for (int lhsWI = 0; lhsWI < lhs_nhist; lhsWI++) {
897 bool done = false;
898
899 // List of detectors on lhs side
900 const auto &lhsDets = lhs->getSpectrum(lhsWI).getDetectorIDs();
901
902 // ----------------- Matching Workspace Indices and Detector IDs
903 // --------------------------------------
904 // First off, try to match the workspace indices. Most times, this will be
905 // ok right away.
906 int64_t rhsWI = lhsWI;
907 if (rhsWI < rhs_nhist) // don't go out of bounds
908 {
909 // Get the detector IDs at that workspace index.
910 const auto &rhsDets = rhs->getSpectrum(rhsWI).getDetectorIDs();
911
912 // Checks that lhsDets is a subset of rhsDets
913 if (std::includes(rhsDets.begin(), rhsDets.end(), lhsDets.begin(), lhsDets.end())) {
914 // We found the workspace index right away. No need to keep looking
915 (*table)[lhsWI] = rhsWI;
916 done = true;
917 }
918 }
919
920 // ----------------- Scrambled Detector IDs with one Detector per Spectrum
921 // --------------------------------------
922 if (!done && (lhsDets.size() == 1)) {
923 // Didn't find it. Try to use the RHS map.
924
925 // First, we have to get the (single) detector ID of the LHS
926 auto lhsDets_it = lhsDets.cbegin();
927 detid_t lhs_detector_ID = *lhsDets_it;
928
929 // Now we use the RHS map to find it. This only works if both the lhs and
930 // rhs have 1 detector per pixel
931 auto map_it = rhs_det_to_wi.find(lhs_detector_ID);
932 if (map_it != rhs_det_to_wi.end()) {
933 rhsWI = map_it->second; // This is the workspace index in the RHS that
934 // matched lhs_detector_ID
935 } else {
936 // Did not find it!
937 rhsWI = -1; // Marker to mean its not in the LHS.
938
939 // std::ostringstream mess;
940 // mess << "BinaryOperation: cannot find a RHS spectrum that
941 // contains the detectors in LHS workspace index " << lhsWI
942 // << "\n";
943 // throw std::runtime_error(mess.str());
944 }
945 (*table)[lhsWI] = rhsWI;
946 done = true; // Great, we did it.
947 }
948
949 // ----------------- LHS detectors are subset of RHS, which are Grouped
950 // --------------------------------------
951 if (!done) {
952
953 // Didn't find it? Now we need to iterate through the output workspace to
954 // match the detector ID.
955 // NOTE: This can be SUPER SLOW!
956 for (rhsWI = 0; rhsWI < static_cast<int64_t>(rhs_nhist); rhsWI++) {
957 const auto &rhsDets = rhs->getSpectrum(rhsWI).getDetectorIDs();
958
959 // Checks that lhsDets is a subset of rhsDets
960 if (std::includes(rhsDets.begin(), rhsDets.end(), lhsDets.begin(), lhsDets.end())) {
961 // This one is right. Now we can stop looking.
962 (*table)[lhsWI] = rhsWI;
963 done = true;
964 continue;
965 }
966 }
967 }
968
969 // ------- Still nothing ! -----------
970 if (!done) {
971 (*table)[lhsWI] = -1;
972
973 // std::ostringstream mess;
974 // mess << "BinaryOperation: cannot find a RHS spectrum that
975 // contains the detectors in LHS workspace index " << lhsWI <<
976 // "\n";
977 // throw std::runtime_error(mess.str());
978 }
979 }
980
981 return table;
982}
983
984Parallel::ExecutionMode
985BinaryOperation::getParallelExecutionMode(const std::map<std::string, Parallel::StorageMode> &storageModes) const {
986 if (static_cast<bool>(getProperty("AllowDifferentNumberSpectra")))
987 return Parallel::ExecutionMode::Invalid;
988 auto lhs = storageModes.find(inputPropName1())->second;
989 auto rhs = storageModes.find(inputPropName2())->second;
990 // Two identical modes is ok
991 if (lhs == rhs)
992 return getCorrespondingExecutionMode(storageModes.begin()->second);
993 // Mode <X> times Cloned is ok if the cloned workspace is WorkspaceSingleValue
994 if (rhs == Parallel::StorageMode::Cloned) {
996 if (std::dynamic_pointer_cast<const WorkspaceSingleValue>(ws))
997 return getCorrespondingExecutionMode(lhs);
998 }
999 // Other options are not ok (e.g., MasterOnly times Distributed)
1000 return Parallel::ExecutionMode::Invalid;
1001}
1002
1003} // namespace Mantid::Algorithms
const std::vector< double > & rhs
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
#define PARALLEL_START_INTERRUPT_REGION
Begins a block to skip processing is the algorithm has been interupted Note the end of the block if n...
Definition: MultiThreaded.h:94
#define PARALLEL_FOR_NO_WSP_CHECK()
#define PARALLEL_CRITICAL(name)
#define PARALLEL_END_INTERRUPT_REGION
Ends a block to skip processing is the algorithm has been interupted Note the start of the block if n...
#define PARALLEL_FOR_IF(condition)
Empty definitions - to enable set your complier to enable openMP.
#define PARALLEL_CHECK_INTERRUPT_REGION
Adds a check after a Parallel region to see if it was interupted.
#define UNUSED_ARG(x)
Function arguments are sometimes unused in certain implmentations but are required for documentation ...
Definition: System.h:64
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
Definition: Algorithm.cpp:1913
void clear() override
Clears all properties under management.
Definition: Algorithm.cpp:2125
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
virtual std::shared_ptr< Algorithm > createChildAlgorithm(const std::string &name, const double startProgress=-1., const double endProgress=-1., const bool enableLogging=true, const int &version=-1)
Create a Child Algorithm.
Definition: Algorithm.cpp:842
Kernel::Logger & g_log
Definition: Algorithm.h:451
const std::string name() const override=0
function to return a name of the algorithm, must be overridden in all algorithms
virtual void clearData()=0
Base MatrixWorkspace Abstract Class.
virtual ISpectrum & getSpectrum(const size_t index)=0
Return the underlying ISpectrum ptr at the given workspace index.
std::map< size_t, double > MaskList
Masked bins for each spectrum are stored as a set of pairs containing <bin index, weight>
API::SpectrumInfo is an intermediate step towards a SpectrumInfo that is part of Instrument-2....
Definition: SpectrumInfo.h:53
void setMasked(const size_t index, bool masked)
Set the mask flag of the spectrum with given index.
bool hasDetectors(const size_t index) const
Returns true if the spectrum is associated with detectors in the instrument.
bool isMasked(const size_t index) const
Returns true if the detector(s) associated with the spectrum are masked.
A property class for workspaces.
bool m_do2D_even_for_SingleColumn_on_rhs
Special case for plus/minus: if there is only one bin on the RHS, use the 2D method (appending event ...
void doSingleSpectrum()
Called when the m_rhs operand is a single spectrum.
virtual void operateOnRun(const API::Run &lhs, const API::Run &rhs, API::Run &ans) const
Only overridden by operations that affect the properties of the run (e.g.
size_t m_rhsBlocksize
Cache for RHS workspace's blocksize.
virtual bool propagateSpectraMask(const API::SpectrumInfo &lhsSpectrumInfo, const API::SpectrumInfo &rhsSpectrumInfo, const int64_t index, API::MatrixWorkspace &out, API::SpectrumInfo &outSpectrumInfo)
Checks if the spectra at the given index of either input workspace is masked.
void doSingleValue()
Called when the rhs operand is a single value.
void propagateBinMasks(const API::MatrixWorkspace_const_sptr &rhs, const API::MatrixWorkspace_sptr &out)
Copies any bin masking from the smaller/rhs input workspace to the output.
virtual void setOutputUnits(const API::MatrixWorkspace_const_sptr lhs, const API::MatrixWorkspace_const_sptr rhs, API::MatrixWorkspace_sptr out)
Should be overridden by operations that need to manipulate the units of the output workspace.
static BinaryOperationTable_sptr buildBinaryOperationTable(const API::MatrixWorkspace_const_sptr &lhs, const API::MatrixWorkspace_const_sptr &rhs)
Build up an BinaryOperationTable for performing a binary operation e.g.
bool m_keepEventWorkspace
Variable set to true if the operation allows the output to stay as an EventWorkspace.
std::shared_ptr< BinaryOperationTable > BinaryOperationTable_sptr
virtual bool checkCompatibility(const API::MatrixWorkspace_const_sptr lhs, const API::MatrixWorkspace_const_sptr rhs) const
Checks the compatibility of the two workspaces.
DataObjects::EventWorkspace_sptr m_eout
Output EventWorkspace.
virtual std::string checkSizeCompatibility(const API::MatrixWorkspace_const_sptr lhs, const API::MatrixWorkspace_const_sptr rhs) const
Checks the overall size compatibility of two workspaces.
void init() override
Initialisation method.
void doSingleColumn()
Called when the m_rhs operand is a 2D workspace of single values.
API::MatrixWorkspace_sptr m_out
Output workspace.
virtual void execEvent(DataObjects::EventWorkspace_const_sptr lhs, DataObjects::EventWorkspace_const_sptr rhs)
Execution method for event workspaces, to be overridden as needed.
bool handleSpecialDivideMinus()
Special handling for 1-WS and 1/WS.
OperandType getOperandType(const API::MatrixWorkspace_const_sptr &ws)
Get the type of operand from a workspace.
bool m_AllowDifferentNumberSpectra
The property value.
bool m_useHistogramForRhsEventWorkspace
Are we going to use the histogram representation of the RHS event list when performing the operation?...
bool m_matchXSize
matchXSize set to true if the X sizes of histograms must match.
bool m_flipSides
flipSides set to true if the rhs and lhs operands should be flipped - for commutative binary operatio...
virtual std::string outputPropName() const
The name of the output workspace property.
DataObjects::EventWorkspace_const_sptr m_erhs
Right-hand side EventWorkspace.
virtual bool checkEventCompatibility(const API::MatrixWorkspace_const_sptr lhs, const API::MatrixWorkspace_const_sptr rhs)
Checks the compatibility of event-based processing of the two workspaces.
std::unique_ptr< API::Progress > m_progress
Progress reporting.
virtual void performEventBinaryOperation(DataObjects::EventList &lhs, const DataObjects::EventList &rhs)
Carries out the binary operation IN-PLACE on a single EventList, with another EventList as the right-...
virtual void performBinaryOperation(const HistogramData::Histogram &lhs, const HistogramData::Histogram &rhs, HistogramData::HistogramY &YOut, HistogramData::HistogramE &EOut)=0
Carries out the binary operation on a single spectrum, with another spectrum as the right-hand operan...
void do2D(bool mismatchedSpectra)
Called when the two workspaces are the same size.
size_t m_lhsBlocksize
Cache for LHS workspace's blocksize.
void exec() override
Executes the algorithm.
virtual std::string inputPropName1() const
The name of the first input workspace property.
DataObjects::EventWorkspace_const_sptr m_elhs
Left-hand side EventWorkspace.
API::MatrixWorkspace_const_sptr m_rhs
Right-hand side workspace.
bool m_ClearRHSWorkspace
Flag to clear RHS workspace in binary operation.
Parallel::ExecutionMode getParallelExecutionMode(const std::map< std::string, Parallel::StorageMode > &storageModes) const override
Get correct execution mode based on input storage modes for an MPI run.
virtual void checkRequirements()
Check what operation will be needed in order to apply the operation to these two types of workspaces.
virtual std::string inputPropName2() const
The name of the second input workspace property.
API::MatrixWorkspace_const_sptr m_lhs
Left-hand side workspace.
A class for holding :
Definition: EventList.h:56
Marks code as not implemented yet.
Definition: Exception.h:138
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void error(const std::string &msg)
Logs at error level.
Definition: Logger.cpp:77
The concrete, templated class for properties.
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::shared_ptr< const EventWorkspace > EventWorkspace_const_sptr
shared pointer to a const Workspace2D
std::shared_ptr< const WorkspaceSingleValue > WorkspaceSingleValue_const_sptr
std::shared_ptr< const Unit > Unit_const_sptr
Shared pointer to the Unit base class (const version)
Definition: Unit.h:231
std::enable_if< std::is_pointer< Arg >::value, bool >::type threadSafe(Arg workspace)
Thread-safety check Checks the workspace to ensure it is suitable for multithreaded access.
Definition: MultiThreaded.h:22
int32_t detid_t
Typedef for a detector ID.
Definition: SpectrumInfo.h:21
std::unordered_map< detid_t, size_t > detid2index_map
Map with key = detector ID, value = workspace index.
std::vector< double > MantidVec
typedef for the data storage used in Mantid matrix workspaces
Definition: cow_ptr.h:172
static bool matchingBins(const MatrixWorkspace &ws1, const MatrixWorkspace &ws2, const bool firstOnly=false)
Checks whether the bins (X values) of two workspace are the same.
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54