From e6ca2c0d81c8bcc6ee6bcb83c0ba20bef950abdd Mon Sep 17 00:00:00 2001
From: Zamaan Raza <z.raza@fz-juelich.de>
Date: Thu, 12 Oct 2023 10:40:50 +0200
Subject: [PATCH 01/22] Replace shell with d max
---
gui/subframe_merge/SubframeMerge.cpp | 5 +++--
1 file changed, 3 insertions(+), 2 deletions(-)
diff --git a/gui/subframe_merge/SubframeMerge.cpp b/gui/subframe_merge/SubframeMerge.cpp
index d6e98525f..af54a362b 100644
--- a/gui/subframe_merge/SubframeMerge.cpp
+++ b/gui/subframe_merge/SubframeMerge.cpp
@@ -739,7 +739,8 @@ void SubframeMerge::refreshGraph(int column)
QVector<double> sum_yvals;
QVector<double> profile_yvals;
for (int i = 0; i < nshells; ++i) {
- xvals.push_back(double(i));
+ double d_max = _sum_shell_model->item(i, 1)->data(Qt::DisplayRole).value<double>();
+ xvals.push_back(d_max);
double sum_val = _sum_shell_model->item(i, column)->data(Qt::DisplayRole).value<double>();
double profile_val =
_profile_shell_model->item(i, column)->data(Qt::DisplayRole).value<double>();
@@ -765,7 +766,7 @@ void SubframeMerge::refreshGraph(int column)
_statistics_plot->graph(1)->setName("Profile");
_statistics_plot->legend->setVisible(true);
- _statistics_plot->xAxis->setLabel("shell");
+ _statistics_plot->xAxis->setLabel("d max");
_statistics_plot->yAxis->setLabel(_plottable_statistics->itemText(column));
_statistics_plot->setNotAntialiasedElements(QCP::aeAll);
--
GitLab
From 02e8d5c5eb2a8caacf01deae13e8f8c3fecf9a07 Mon Sep 17 00:00:00 2001
From: Zamaan Raza <z.raza@fz-juelich.de>
Date: Tue, 14 May 2024 08:46:37 +0200
Subject: [PATCH 02/22] Buffer all frames by default in GUI
---
core/data/DataSet.cpp | 7 ++++++-
core/data/DataSet.h | 6 +++++-
gui/models/Session.cpp | 2 ++
3 files changed, 13 insertions(+), 2 deletions(-)
diff --git a/core/data/DataSet.cpp b/core/data/DataSet.cpp
index ee52ffebc..58d00a753 100644
--- a/core/data/DataSet.cpp
+++ b/core/data/DataSet.cpp
@@ -53,7 +53,11 @@
namespace ohkl {
DataSet::DataSet(const std::string& dataset_name, Diffractometer* diffractometer)
- : _diffractometer{diffractometer}, _states(nullptr), _total_histogram(nullptr), _buffered(false)
+ : _diffractometer{diffractometer}
+ , _states(nullptr)
+ , _total_histogram(nullptr)
+ , _buffered(false)
+ , _all_frames_buffered(true)
{
setName(dataset_name);
if (!_diffractometer)
@@ -513,6 +517,7 @@ void DataSet::initBuffer(bool bufferAll)
_frame_buffer.at(idx) = std::make_unique<Eigen::MatrixXi>(_reader->data(idx));
}
_buffered = true;
+ _all_frames_buffered = bufferAll;
ohklLog(Level::Debug, "DataSet::initBuffer: ", _name, " buffered");
}
diff --git a/core/data/DataSet.h b/core/data/DataSet.h
index 54a401aa9..d6c2201d8 100644
--- a/core/data/DataSet.h
+++ b/core/data/DataSet.h
@@ -174,6 +174,8 @@ class DataSet {
void initBuffer(bool bufferAll = true);
//! Clear the frame buffer
void clearBuffer();
+ //! Check whether all frames are buffered
+ bool allFramesBuffered() const { return _all_frames_buffered; };
virtual void setNFrames(std::size_t nframes) { std::ignore = nframes; };
@@ -204,8 +206,10 @@ class DataSet {
//! Buffer for image data
std::vector<std::unique_ptr<Eigen::MatrixXi>> _frame_buffer;
- //! Whether or not the buffer is active
+ //! Buffer is active
bool _buffered;
+ //! All frames are buffered
+ bool _all_frames_buffered;
};
/*! @}*/
diff --git a/gui/models/Session.cpp b/gui/models/Session.cpp
index beb8828be..fbf3db7a1 100644
--- a/gui/models/Session.cpp
+++ b/gui/models/Session.cpp
@@ -572,6 +572,8 @@ void Session::onDataChanged()
if (!gSession->currentProject()->hasDataSet())
return;
+ _data_combo->currentData()->initBuffer(true);
+
double x_offset =
_data_combo->currentData()->diffractometer()->source().selectedMonochromator().xOffset();
double y_offset =
--
GitLab
From e70260d1c306bc463b929e9158bfb84c417d7b1d Mon Sep 17 00:00:00 2001
From: Zamaan Raza <z.raza@fz-juelich.de>
Date: Tue, 14 May 2024 08:49:50 +0200
Subject: [PATCH 03/22] Remove OpenMP pragmas
---
core/algo/RefinementBatch.cpp | 2 --
core/algo/Refiner.cpp | 2 --
core/experiment/PeakFinder.cpp | 16 ++++++----------
core/integration/IIntegrator.cpp | 2 --
core/peak/Qs2Events.cpp | 16 ++++------------
core/shape/ShapeModel.cpp | 1 -
6 files changed, 10 insertions(+), 29 deletions(-)
diff --git a/core/algo/RefinementBatch.cpp b/core/algo/RefinementBatch.cpp
index 37e721354..72f3ece4f 100644
--- a/core/algo/RefinementBatch.cpp
+++ b/core/algo/RefinementBatch.cpp
@@ -224,7 +224,6 @@ int RefinementBatch::qSpaceResiduals(Eigen::VectorXd& fvec) const
UnitCell uc = _cell->fromParameters(_u0, _uOffsets, _cellParameters);
const Eigen::Matrix3d& UB = uc.reciprocalBasis();
- // #pragma omp parallel for
for (unsigned int i = 0; i < _peaks.size(); ++i) {
const Eigen::RowVector3d q0 = _peaks[i]->q().rowVector();
const Eigen::RowVector3d q1 = _hkls[i] * UB;
@@ -249,7 +248,6 @@ int RefinementBatch::realSpaceResiduals(Eigen::VectorXd& fvec) const
UnitCell uc = _cell->fromParameters(_u0, _uOffsets, _cellParameters);
const Eigen::Matrix3d& UB = uc.reciprocalBasis();
- // #pragma omp parallel for
for (unsigned int i = 0; i < _peaks.size(); ++i) {
const Eigen::RowVector3d x0 = _peaks[i]->shape().center();
auto data = _peaks[i]->dataSet();
diff --git a/core/algo/Refiner.cpp b/core/algo/Refiner.cpp
index bb8eded89..1afb5a376 100644
--- a/core/algo/Refiner.cpp
+++ b/core/algo/Refiner.cpp
@@ -273,10 +273,8 @@ bool Refiner::refine(sptrDataSet data, const std::vector<Peak3D*> peaks, sptrUni
return false;
unsigned int failed_batches = 0;
- // #pragma omp parallel for
for (auto&& batch : _batches) {
if (!batch.refine(_params->max_iter)) {
- // #pragma omp atomic
++failed_batches;
}
}
diff --git a/core/experiment/PeakFinder.cpp b/core/experiment/PeakFinder.cpp
index 7080b460f..e3034d18b 100644
--- a/core/experiment/PeakFinder.cpp
+++ b/core/experiment/PeakFinder.cpp
@@ -238,7 +238,6 @@ void PeakFinder::findPrimaryBlobs(
// Iterate on all pixels in the image
int nframes = 0;
-#pragma omp for schedule(dynamic, DYNAMIC_CHUNK)
for (size_t idx = begin; idx < end; ++idx) {
++nframes;
@@ -318,15 +317,12 @@ void PeakFinder::findPrimaryBlobs(
labels[index2D] = labels2[index2D] = label;
index2D++;
auto value = frame_data(row, col);
-#pragma omp critical(dataupdate)
- {
- // Create a new blob if necessary
- if (newlabel)
- blobs.insert(std::make_pair(label, Blob3D(col, row, idx, value)));
- else {
- auto it = blobs.find(label);
- it->second.addPoint(col, row, idx, value);
- }
+ // Create a new blob if necessary
+ if (newlabel)
+ blobs.insert(std::make_pair(label, Blob3D(col, row, idx, value)));
+ else {
+ auto it = blobs.find(label);
+ it->second.addPoint(col, row, idx, value);
}
}
}
diff --git a/core/integration/IIntegrator.cpp b/core/integration/IIntegrator.cpp
index 2b18162c3..8a74f1a0a 100644
--- a/core/integration/IIntegrator.cpp
+++ b/core/integration/IIntegrator.cpp
@@ -172,7 +172,6 @@ void IIntegrator::integrate(
current_peak->updateMask(mask, idx);
}
-#pragma omp parallel for
for (auto peak : peaks) {
auto* current_peak = regions.at(peak).get();
// Check whether the peak intersects a mask
@@ -221,7 +220,6 @@ void IIntegrator::integrate(
peak->setIntegrationFlag(
RejectionFlag::SaturatedPixel, _params.integrator_type);
} else {
-#pragma omp atomic
++nfailures;
// This is a fallback. The RejectionFlag should have been set by this point.
peak->setIntegrationFlag(
diff --git a/core/peak/Qs2Events.cpp b/core/peak/Qs2Events.cpp
index ee2ce955f..eb1064611 100644
--- a/core/peak/Qs2Events.cpp
+++ b/core/peak/Qs2Events.cpp
@@ -54,15 +54,11 @@ std::vector<std::pair<MillerIndex, DetectorEvent>> algo::qMap2Events(
}
// for each sample q, determine the rotation that makes it intersect the Ewald sphere
-#pragma omp parallel for
for (const auto& [hkl, sample_q] : sample_qs) {
std::vector<DetectorEvent> new_events =
qVector2Events(sample_q, states, detector, n_intervals);
-#pragma omp critical(dataupdate)
- {
- for (auto event : new_events)
- events.push_back({hkl, event});
- }
+ for (auto event : new_events)
+ events.push_back({hkl, event});
if (handler)
handler->setProgress(++count * 100.0 / sample_qs.size());
}
@@ -93,15 +89,11 @@ std::vector<DetectorEvent> algo::qVectorList2Events(
}
// for each sample q, determine the rotation that makes it intersect the Ewald sphere
-#pragma omp parallel for
for (const ReciprocalVector& sample_q : sample_qs) {
std::vector<DetectorEvent> new_events =
qVector2Events(sample_q, states, detector, n_intervals);
-#pragma omp critical(dataupdate)
- {
- for (auto event : new_events)
- events.emplace_back(event);
- }
+ for (auto event : new_events)
+ events.emplace_back(event);
if (handler)
handler->setProgress(++count * 100.0 / sample_qs.size());
}
diff --git a/core/shape/ShapeModel.cpp b/core/shape/ShapeModel.cpp
index f591508be..942deec56 100644
--- a/core/shape/ShapeModel.cpp
+++ b/core/shape/ShapeModel.cpp
@@ -446,7 +446,6 @@ void ShapeModel::setPredictedShapes(PeakCollection* peaks)
_handler->setProgress(0);
}
-#pragma omp parallel for
for (auto peak : peaks->getPeakList()) {
// Skip the peak if any error occur when computing its mean covariance (e.g.
// too few or no neighbouring peaks found)
--
GitLab
From 04bee7344b8894dd804366133f72d3d931caefc4 Mon Sep 17 00:00:00 2001
From: Zamaan Raza <z.raza@fz-juelich.de>
Date: Tue, 14 May 2024 10:58:01 +0200
Subject: [PATCH 04/22] Correct reference values for Profile3D integration test
---
.../functional/TestProfile3DIntegration.cpp | 34 +++++++++++++------
1 file changed, 23 insertions(+), 11 deletions(-)
diff --git a/test/cpp/functional/TestProfile3DIntegration.cpp b/test/cpp/functional/TestProfile3DIntegration.cpp
index 6569dcd5d..c0d9c3321 100644
--- a/test/cpp/functional/TestProfile3DIntegration.cpp
+++ b/test/cpp/functional/TestProfile3DIntegration.cpp
@@ -38,16 +38,19 @@ TEST_CASE("test/data/TestProfile3DIntegration.cpp", "")
{
// reference values
const int ref_n_peaks = 59209;
- const int ref_n_integrated_peaks = 56160;
- const double ref_rpim_overall = 0.0798;
- const double ref_ccstar_overall = 0.9972;
+ const int ref_n_integrated_peaks = 55986;
+ const double ref_rpim_overall = 0.0763;
+ const double ref_ccstar_overall = 0.9973;
+ const double ref_cchalf_overall = 0.9892;
const std::vector<double> ref_rpim =
- {0.0267, 0.0605, 0.1028, 0.1362, 0.1766, 0.2219, 0.2649, 0.3084, 0.3632, 0.4278};
+ {0.0267, 0.0605, 0.1019, 0.1325, 0.1729, 0.2113, 0.2474, 0.2790, 0.3272, 0.3974};
const std::vector<double> ref_ccstar =
- {0.9988, 0.9957, 0.9871, 0.9778, 0.9670, 0.9406, 0.9223, 0.9002, 0.8395, 0.7645};
+ {0.9988, 0.9956, 0.9873, 0.9786, 0.9685, 0.9464, 0.9330, 0.9283, 0.8837, 0.8071};
+ const std::vector<double> ref_cchalf =
+ {0.9953, 0.9826, 0.9507, 0.9187, 0.8833, 0.8114, 0.7707, 0.7570, 0.6406, 0.4831};
const int eps_peaks = 100;
- const double eps_stats = 0.02;
+ const double eps_stats = 0.01;
const std::string filename = "Trypsin-pxsum.ohkl";
ohkl::Experiment experiment("Trypsin", "BioDiff");
@@ -66,7 +69,7 @@ TEST_CASE("test/data/TestProfile3DIntegration.cpp", "")
shape_params.sigma_d = found->sigmaD();
shape_params.sigma_m = found->sigmaM();
shape_params.strength_min = 1.0;
- shape_params.neighbour_range_pixels = 200;
+ shape_params.neighbour_range_pixels = 50;
shape_params.neighbour_range_frames = 10;
found->buildShapeModel(data, shape_params);
ohkl::ShapeModel* shapes = found->shapeModel();
@@ -80,9 +83,9 @@ TEST_CASE("test/data/TestProfile3DIntegration.cpp", "")
params->discard_saturated = true;
params->skip_masked = true;
params->use_max_strength = true;
- params->max_strength = 5.0;
+ params->max_strength = 1.0;
params->use_max_d = true;
- params->max_d = 2.56;
+ params->max_d = 2.5;
integrator->integratePeaks(data, peaks, params, shapes);
@@ -111,19 +114,28 @@ TEST_CASE("test/data/TestProfile3DIntegration.cpp", "")
auto* shell_quality = merger->profileShellQuality();
const double rpim_overall = overall_quality->shells[0].Rpim;
const double ccstar_overall = overall_quality->shells[0].CCstar;
+ const double cchalf_overall = overall_quality->shells[0].CChalf;
std::vector<double> rpim;
std::vector<double> ccstar;
+ std::vector<double> cchalf;
for (int i = 0; i < merge_params->n_shells; ++i) {
rpim.push_back(shell_quality->shells[i].Rpim);
ccstar.push_back(shell_quality->shells[i].CCstar);
+ cchalf.push_back(shell_quality->shells[i].CChalf);
}
CHECK(rpim_overall < ref_rpim_overall + eps_stats);
- CHECK(ccstar_overall < ref_ccstar_overall + eps_stats);
+ CHECK(ccstar_overall > ref_ccstar_overall - eps_stats);
+ CHECK(cchalf_overall > ref_cchalf_overall - eps_stats);
+ std::cout << "rpim_overall = " << rpim_overall << std::endl;
+ std::cout << "ccstar_overall = " << ccstar_overall << std::endl;
+ std::cout << "cchalf_overall = " << cchalf_overall << std::endl;
for (int i = 0; i < rpim.size(); ++i) {
CHECK(rpim[i] < ref_rpim[i] + eps_stats);
- CHECK(ccstar[i] < ref_ccstar[i] + eps_stats);
+ CHECK(ccstar[i] > ref_ccstar[i] - eps_stats);
+ CHECK(cchalf[i] > ref_cchalf[i] - eps_stats);
std::cout << rpim[i] << " " << ref_rpim[i] << std::endl;
std::cout << ccstar[i] << " " << ref_ccstar[i] << std::endl;
+ std::cout << cchalf[i] << " " << ref_cchalf[i] << std::endl << std::endl;
}
}
--
GitLab
From 418647e2b12e076a4fbfa356db0bb64aae2974d2 Mon Sep 17 00:00:00 2001
From: Zamaan Raza <z.raza@fz-juelich.de>
Date: Tue, 14 May 2024 12:55:18 +0200
Subject: [PATCH 05/22] Generate Integrators as-needed using IntegratorFactory
---
core/experiment/Integrator.cpp | 32 +++--------
core/experiment/Integrator.h | 9 +--
core/integration/GaussianIntegrator.cpp | 6 +-
core/integration/GaussianIntegrator.h | 2 +-
core/integration/ISigmaIntegrator.cpp | 2 +-
core/integration/IntegratorFactory.cpp | 57 +++++++++++++++++++
core/integration/IntegratorFactory.h | 37 ++++++++++++
core/integration/PixelSumIntegrator.cpp | 6 +-
core/integration/PixelSumIntegrator.h | 2 +-
core/integration/ShapeIntegrator.cpp | 2 +-
.../unit/integrate/TestPeakIntegration.cpp | 4 +-
11 files changed, 121 insertions(+), 38 deletions(-)
create mode 100644 core/integration/IntegratorFactory.cpp
create mode 100644 core/integration/IntegratorFactory.h
diff --git a/core/experiment/Integrator.cpp b/core/experiment/Integrator.cpp
index 4b2f8ef76..74bc47bb0 100644
--- a/core/experiment/Integrator.cpp
+++ b/core/experiment/Integrator.cpp
@@ -17,11 +17,7 @@
#include "base/utils/Logger.h"
#include "core/experiment/DataHandler.h"
#include "core/experiment/PeakFinder.h"
-#include "core/integration/GaussianIntegrator.h"
-#include "core/integration/ISigmaIntegrator.h"
-#include "core/integration/PixelSumIntegrator.h"
-#include "core/integration/Profile1DIntegrator.h"
-#include "core/integration/Profile3DIntegrator.h"
+#include "core/integration/IIntegrator.h"
#include "core/integration/ShapeIntegrator.h"
#include "core/peak/Peak3D.h"
#include "core/shape/PeakCollection.h"
@@ -30,19 +26,8 @@
namespace ohkl {
Integrator::Integrator(std::shared_ptr<DataHandler> data_handler)
- : _handler(nullptr), _data_handler(data_handler)
+ : _handler(nullptr), _data_handler(data_handler), _factory()
{
- _integrator_map.clear();
- _integrator_map.insert(
- std::make_pair(IntegratorType::PixelSum, std::make_unique<PixelSumIntegrator>(true, true)));
- _integrator_map.insert(
- std::make_pair(IntegratorType::Gaussian, std::make_unique<GaussianIntegrator>(true, true)));
- _integrator_map.insert(
- std::make_pair(IntegratorType::ISigma, std::make_unique<ISigmaIntegrator>()));
- _integrator_map.insert(
- std::make_pair(IntegratorType::Profile1D, std::make_unique<Profile1DIntegrator>()));
- _integrator_map.insert(
- std::make_pair(IntegratorType::Profile3D, std::make_unique<Profile3DIntegrator>()));
_params = std::make_unique<IntegrationParameters>();
}
@@ -53,12 +38,7 @@ DataHandler* Integrator::getDataHandler()
IIntegrator* Integrator::getIntegrator(const IntegratorType integrator_type) const
{
- std::map<IntegratorType, std::unique_ptr<IIntegrator>>::const_iterator it;
- for (it = _integrator_map.begin(); it != _integrator_map.end(); ++it) {
- if (it->first == integrator_type)
- return it->second.get();
- }
- return nullptr;
+ return _factory.create(integrator_type);
}
void Integrator::integratePeaks(
@@ -82,6 +62,8 @@ void Integrator::integratePeaks(
if (peak->enabled())
++_n_valid;
}
+
+ delete integrator;
}
void Integrator::integratePeaks(
@@ -105,6 +87,8 @@ void Integrator::integratePeaks(
if (peak->enabled())
++_n_valid;
}
+
+ delete integrator;
}
void Integrator::integrateFoundPeaks(PeakFinder* peak_finder)
@@ -125,6 +109,8 @@ void Integrator::integrateFoundPeaks(PeakFinder* peak_finder)
peak_finder->setIntegrated(true);
if (_params->use_gradient)
peak_finder->setBkgGradient(true);
+
+ delete integrator;
}
void Integrator::integrateShapeModel(
diff --git a/core/experiment/Integrator.h b/core/experiment/Integrator.h
index c09b2dcb0..270450b85 100644
--- a/core/experiment/Integrator.h
+++ b/core/experiment/Integrator.h
@@ -15,19 +15,19 @@
#define OPENHKL_CORE_EXPERIMENT_INTEGRATOR_H
#include "core/data/DataTypes.h"
-#include "core/integration/IIntegrator.h"
+#include "core/integration/IntegratorFactory.h"
#include <map>
#include <string>
namespace ohkl {
-using IntegratorMap = std::map<IntegratorType, std::unique_ptr<ohkl::IIntegrator>>;
-
class PeakCollection;
class PeakFinder;
class DataHandler;
+class IIntegrator;
struct ShapeModelParameters;
+struct IntegrationParameters;
/*! \addtogroup python_api
* @{*/
@@ -75,12 +75,13 @@ class Integrator {
private:
sptrProgressHandler _handler;
- IntegratorMap _integrator_map;
std::shared_ptr<DataHandler> _data_handler;
std::shared_ptr<IntegrationParameters> _params;
unsigned int _n_peaks;
unsigned int _n_valid;
+
+ IntegratorFactory _factory;
};
/*! @}*/
diff --git a/core/integration/GaussianIntegrator.cpp b/core/integration/GaussianIntegrator.cpp
index d19bf5420..4571fbcf4 100644
--- a/core/integration/GaussianIntegrator.cpp
+++ b/core/integration/GaussianIntegrator.cpp
@@ -25,10 +25,10 @@
namespace ohkl {
-GaussianIntegrator::GaussianIntegrator(bool fit_center, bool fit_cov)
+GaussianIntegrator::GaussianIntegrator() : IIntegrator()
{
- _params.fit_center = fit_center;
- _params.fit_cov = fit_cov;
+ _params.fit_center = true;
+ _params.fit_cov = true;
}
static Eigen::Matrix3d from_cholesky(const Eigen::VectorXd a)
diff --git a/core/integration/GaussianIntegrator.h b/core/integration/GaussianIntegrator.h
index ed531abec..2ef5e3d66 100644
--- a/core/integration/GaussianIntegrator.h
+++ b/core/integration/GaussianIntegrator.h
@@ -25,7 +25,7 @@ namespace ohkl {
/*! \brief Compute integrated intensity by fitting to an analytic 3D Gaussian.*/
class GaussianIntegrator : public IIntegrator {
public:
- GaussianIntegrator(bool fit_center, bool fit_cov);
+ GaussianIntegrator();
protected:
//! Integrate a peak
diff --git a/core/integration/ISigmaIntegrator.cpp b/core/integration/ISigmaIntegrator.cpp
index 36a7bdbb7..11322aa82 100644
--- a/core/integration/ISigmaIntegrator.cpp
+++ b/core/integration/ISigmaIntegrator.cpp
@@ -23,7 +23,7 @@
namespace ohkl {
-ISigmaIntegrator::ISigmaIntegrator() : PixelSumIntegrator(false, false) { }
+ISigmaIntegrator::ISigmaIntegrator() : PixelSumIntegrator() { }
bool ISigmaIntegrator::compute(
Peak3D* peak, ShapeModel* shape_model, const IntegrationRegion& region)
diff --git a/core/integration/IntegratorFactory.cpp b/core/integration/IntegratorFactory.cpp
new file mode 100644
index 000000000..cf8e4440e
--- /dev/null
+++ b/core/integration/IntegratorFactory.cpp
@@ -0,0 +1,57 @@
+// ***********************************************************************************************
+//
+// OpenHKL: data reduction for single crystal diffraction
+//
+//! @file core/convolve/IntegratorFactory.cpp
+//! @brief Implements class IntegratorFactory
+//!
+//! @homepage https://openhkl.org
+//! @license GNU General Public License v3 or higher (see COPYING)
+//! @copyright Institut Laue-Langevin and Forschungszentrum Jülich GmbH 2016-
+//! @authors see CITATION, MAINTAINER
+//
+// ***********************************************************************************************
+
+#include "core/integration/IntegratorFactory.h"
+
+#include "core/integration/GaussianIntegrator.h"
+#include "core/integration/IIntegrator.h"
+#include "core/integration/ISigmaIntegrator.h"
+#include "core/integration/PixelSumIntegrator.h"
+#include "core/integration/Profile1DIntegrator.h"
+#include "core/integration/Profile3DIntegrator.h"
+#include "core/peak/Peak3D.h"
+
+
+namespace ohkl {
+
+template <typename T> T* create_integrator()
+{
+ return new T();
+}
+
+IntegratorFactory::IntegratorFactory() : _callbacks()
+{
+ _callbacks[IntegratorType::PixelSum] = &create_integrator<PixelSumIntegrator>;
+ _callbacks[IntegratorType::Profile3D] = &create_integrator<Profile3DIntegrator>;
+ _callbacks[IntegratorType::Profile1D] = &create_integrator<Profile1DIntegrator>;
+ _callbacks[IntegratorType::ISigma] = &create_integrator<ISigmaIntegrator>;
+ _callbacks[IntegratorType::Gaussian] = &create_integrator<GaussianIntegrator>;
+}
+
+IIntegrator* IntegratorFactory::create(const IntegratorType& integrator_type) const
+{
+ const auto it = _callbacks.find(integrator_type);
+
+ if (it == _callbacks.end())
+ throw std::runtime_error("IntegratorFactory::create: Invalid integrator type");
+
+ return it->second();
+}
+
+const std::map<IntegratorType, IntegratorFactory::callback>& IntegratorFactory::callbacks() const
+{
+ return _callbacks;
+}
+
+} // namespace ohkl
diff --git a/core/integration/IntegratorFactory.h b/core/integration/IntegratorFactory.h
new file mode 100644
index 000000000..80695a056
--- /dev/null
+++ b/core/integration/IntegratorFactory.h
@@ -0,0 +1,37 @@
+// ***********************************************************************************************
+//
+// OpenHKL: data reduction for single crystal diffraction
+//
+//! @file core/integration/IntegratorFactory.h
+//! @brief Defines class IntegratorFactory
+//!
+//! @homepage https://openhkl.org
+//! @license GNU General Public License v3 or higher (see COPYING)
+//! @copyright Institut Laue-Langevin and Forschungszentrum Jülich GmbH 2016-
+//! @authors see CITATION, MAINTAINER
+//
+// ***********************************************************************************************
+
+#ifndef OHKL_CORE_INTEGRATE_INTEGRATORFACTORY_H
+#define OHKL_CORE_INTEGRATE_INTEGRATORFACTORY_H
+
+#include "core/integration/IIntegrator.h"
+
+namespace ohkl {
+
+class IntegratorFactory {
+ public:
+ using callback = std::function<IIntegrator*()>;
+
+ IntegratorFactory();
+ IIntegrator* create(const IntegratorType& integrator_type) const;
+
+ const std::map<IntegratorType, callback>& callbacks() const;
+
+ private:
+ std::map<IntegratorType, callback> _callbacks;
+};
+
+} // namespace ohkl
+
+#endif // OHKL_CORE_INTEGRATE_INTEGRATORFACTORY_H
diff --git a/core/integration/PixelSumIntegrator.cpp b/core/integration/PixelSumIntegrator.cpp
index 592f7aa8b..ed5aae273 100644
--- a/core/integration/PixelSumIntegrator.cpp
+++ b/core/integration/PixelSumIntegrator.cpp
@@ -85,10 +85,10 @@ std::pair<Intensity, Intensity> compute_background(
} // namespace
-PixelSumIntegrator::PixelSumIntegrator(bool fit_center, bool fit_covariance) : IIntegrator()
+PixelSumIntegrator::PixelSumIntegrator() : IIntegrator()
{
- _params.fit_center = fit_center;
- _params.fit_cov = fit_covariance;
+ _params.fit_center = true;
+ _params.fit_cov = true;
}
bool PixelSumIntegrator::compute(Peak3D* peak, ShapeModel*, const IntegrationRegion& region)
diff --git a/core/integration/PixelSumIntegrator.h b/core/integration/PixelSumIntegrator.h
index e19d05704..60139d9a9 100644
--- a/core/integration/PixelSumIntegrator.h
+++ b/core/integration/PixelSumIntegrator.h
@@ -29,7 +29,7 @@ class PixelSumIntegrator : public IIntegrator {
//! Construct the pixel sum integrator
//! @param fit_center update the peak center as part of integration
//! @param fit_covariance update the peak shape covariance matrix as part of integration
- PixelSumIntegrator(bool fit_center, bool fit_covariance);
+ PixelSumIntegrator();
protected:
//! Integrate a peak
diff --git a/core/integration/ShapeIntegrator.cpp b/core/integration/ShapeIntegrator.cpp
index 07b547f2e..e3bc6a1da 100644
--- a/core/integration/ShapeIntegrator.cpp
+++ b/core/integration/ShapeIntegrator.cpp
@@ -26,7 +26,7 @@ namespace ohkl {
ShapeIntegrator::ShapeIntegrator(
ShapeModel* lib, const AABB& aabb, int nx, int ny, int nz, int subdiv)
- : PixelSumIntegrator(false, false)
+ : PixelSumIntegrator()
, _shape_model(lib)
, _aabb(aabb)
, _nx(nx)
diff --git a/test/cpp/unit/integrate/TestPeakIntegration.cpp b/test/cpp/unit/integrate/TestPeakIntegration.cpp
index 73f7e71d5..4a1564819 100644
--- a/test/cpp/unit/integrate/TestPeakIntegration.cpp
+++ b/test/cpp/unit/integrate/TestPeakIntegration.cpp
@@ -78,8 +78,10 @@ TEST_CASE("test/integrate/Test_6_12_38.cpp", "")
std::vector<ohkl::Peak3D*> peaks;
peaks.push_back(&peak);
- ohkl::PixelSumIntegrator integrator(false, false);
+ ohkl::PixelSumIntegrator integrator;
ohkl::IntegrationParameters params{};
+ params.fit_center = false;
+ params.fit_cov = false;
params.peak_end = 2.7;
params.bkg_begin = 3.0;
params.bkg_end = 4.0;
--
GitLab
From 59f73632bf946bc99c145ffbe308a9d321e8a791 Mon Sep 17 00:00:00 2001
From: Zamaan Raza <z.raza@fz-juelich.de>
Date: Tue, 14 May 2024 13:04:16 +0200
Subject: [PATCH 06/22] Add ParallelFor
---
base/utils/ParallelFor.h | 79 ++++++++++++++++++++++++++++++++++++++++
1 file changed, 79 insertions(+)
create mode 100644 base/utils/ParallelFor.h
diff --git a/base/utils/ParallelFor.h b/base/utils/ParallelFor.h
new file mode 100644
index 000000000..3245fd170
--- /dev/null
+++ b/base/utils/ParallelFor.h
@@ -0,0 +1,79 @@
+// ***********************************************************************************************
+//
+// OpenHKL: data reduction for single crystal diffraction
+//
+//! @file base/utils/ParallelFor.h
+//! @brief Defines class ParallelFor
+//!
+//! @homepage https://openhkl.org
+//! @license GNU General Public License v3 or higher (see COPYING)
+//! @copyright Institut Laue-Langevin and Forschungszentrum Jülich GmbH 2016-
+//! @authors see CITATION, MAINTAINER
+//
+// ***********************************************************************************************
+
+#ifndef OHKL_BASE_UTILS_PARALLELFOR_H
+#define OHKL_BASE_UTILS_PARALLELFOR_H
+
+#include <algorithm>
+#include <functional>
+#include <thread>
+#include <vector>
+
+#include <iostream>
+
+namespace ohkl {
+
+/// @param[in] nb_elements : size of your for loop
+/// @param[in] functor(start, end) :
+/// your function processing a sub chunk of the for loop.
+/// "start" is the first index to process (included) until the index "end"
+/// (excluded)
+/// @code
+/// for(int i = start; i < end; ++i)
+/// computation(i);
+/// @endcode
+/// @param use_threads : enable / disable threads.
+///
+///
+static void parallel_for(
+ unsigned int nb_elements, std::function<void (int start, int end)> functor,
+ bool use_threads = true, unsigned int max_threads = 0)
+{
+ unsigned int nb_threads_hint = std::thread::hardware_concurrency();
+ unsigned int nb_threads = nb_threads_hint == 0 ? 8 : nb_threads_hint;
+ if (max_threads > 0) // Do not exceed the maximum number of threads (prevents OOM)
+ nb_threads = nb_threads > max_threads ? max_threads : nb_threads;
+
+ unsigned int batch_size = nb_elements / nb_threads;
+ unsigned int batch_remainder = nb_elements % nb_threads;
+
+ std::vector<std::thread> my_threads(nb_threads);
+
+ if (use_threads) {
+ // Multithread execution
+ for(unsigned int i = 0; i < nb_threads; ++i) {
+ int start = i * batch_size;
+ my_threads[i] = std::thread(functor, start, start+batch_size);
+ }
+ }
+ else {
+ // Single thread execution (for easy debugging)
+ for(unsigned int i = 0; i < nb_threads; ++i){
+ int start = i * batch_size;
+ functor( start, start+batch_size );
+ }
+ }
+
+ // Deform the elements left
+ int start = nb_threads * batch_size;
+ functor( start, start+batch_remainder);
+
+ // Wait for the other thread to finish their task
+ if (use_threads)
+ std::for_each(my_threads.begin(), my_threads.end(), std::mem_fn(&std::thread::join));
+}
+
+} // namespace ohkl
+
+#endif // OHKL_BASE_UTILS_PARALLELFOR_H
--
GitLab
From 8808afca53471c4389f7d88e70aa7e928b728906 Mon Sep 17 00:00:00 2001
From: Zamaan Raza <z.raza@fz-juelich.de>
Date: Tue, 14 May 2024 15:47:35 +0200
Subject: [PATCH 07/22] Make serial integration thread-safe
add missing definition
---
core/integration/GaussianIntegrator.cpp | 78 ++++++-----
core/integration/GaussianIntegrator.h | 3 +-
core/integration/IIntegrator.cpp | 161 ++++++++++-------------
core/integration/IIntegrator.h | 77 +++++++----
core/integration/ISigmaIntegrator.cpp | 46 +++----
core/integration/ISigmaIntegrator.h | 3 +-
core/integration/PixelSumIntegrator.cpp | 75 ++++++-----
core/integration/PixelSumIntegrator.h | 2 +-
core/integration/Profile1DIntegrator.cpp | 60 +++++----
core/integration/Profile1DIntegrator.h | 3 +-
core/integration/Profile3DIntegrator.cpp | 68 +++++-----
core/integration/Profile3DIntegrator.h | 3 +-
core/integration/ShapeIntegrator.cpp | 27 ++--
core/integration/ShapeIntegrator.h | 3 +-
core/peak/IntegrationRegion.cpp | 94 +++++++------
core/peak/IntegrationRegion.h | 10 +-
core/peak/Peak3D.cpp | 116 ++++++++--------
core/peak/Peak3D.h | 8 +-
core/peak/PeakData.cpp | 17 ++-
core/peak/PeakData.h | 6 +
20 files changed, 464 insertions(+), 396 deletions(-)
diff --git a/core/integration/GaussianIntegrator.cpp b/core/integration/GaussianIntegrator.cpp
index 4571fbcf4..5821f463f 100644
--- a/core/integration/GaussianIntegrator.cpp
+++ b/core/integration/GaussianIntegrator.cpp
@@ -13,6 +13,7 @@
// ***********************************************************************************************
#include "core/integration/GaussianIntegrator.h"
+
#include "base/fit/FitParameters.h"
#include "base/fit/Minimizer.h"
#include "base/geometry/Ellipsoid.h"
@@ -23,15 +24,9 @@
#include <Eigen/Cholesky>
-namespace ohkl {
-
-GaussianIntegrator::GaussianIntegrator() : IIntegrator()
-{
- _params.fit_center = true;
- _params.fit_cov = true;
-}
+namespace {
-static Eigen::Matrix3d from_cholesky(const Eigen::VectorXd a)
+Eigen::Matrix3d from_cholesky(const Eigen::VectorXd a)
{
// Reconstruct Cholesky L factor
Eigen::Matrix3d L;
@@ -46,7 +41,7 @@ static Eigen::Matrix3d from_cholesky(const Eigen::VectorXd a)
return L * L.transpose();
}
-static void residuals(
+void residuals(
Eigen::VectorXd& res, double B, double I, const Eigen::Vector3d x0, const Eigen::VectorXd& a,
const std::vector<Eigen::Vector3d>& x, const std::vector<double>& M, double* pearson)
{
@@ -59,7 +54,7 @@ static void residuals(
double u = 0, v = 0, uu = 0, vv = 0, uv = 0;
for (size_t i = 0; i < n; ++i) {
- Eigen::Vector3d dx = x[i] - x0;
+ const Eigen::Vector3d dx = x[i] - x0;
const double xAx = dx.dot(A * dx);
const double M_pred = B + I * std::exp(-0.5 * xAx);
const double M_obs = M[i];
@@ -85,11 +80,21 @@ static void residuals(
}
}
-bool GaussianIntegrator::compute(
+} // namespace
+
+namespace ohkl {
+
+GaussianIntegrator::GaussianIntegrator() : IIntegrator()
+{
+ _params.fit_center = true;
+ _params.fit_cov = true;
+}
+
+ComputeResult GaussianIntegrator::compute(
Peak3D* peak, ShapeModel* /*shape_model*/, const IntegrationRegion& region)
{
- if (!peak)
- return false;
+ ComputeResult result;
+ result.integrator_type = IntegratorType::Gaussian;
const size_t N = region.peakData().events().size();
@@ -149,39 +154,41 @@ bool GaussianIntegrator::compute(
min.setWeights(wts);
try {
- bool success = min.fit(100);
- if (!success)
- return false;
+ const bool success = min.fit(100);
+ if (!success) {
+ result.integration_flag = RejectionFlag::BadGaussianFit;
+ return result;
+ }
} catch (std::exception& e) {
- ohklLog(Level::Warning, __FUNCTION__, ": Gaussian fit failed: ", e.what());
- return false;
+ result.integration_flag = RejectionFlag::BadGaussianFit;
+ return result;
}
// consistency check: center should still be in dataset!
if (x0(0) < 0 || x0(0) >= peak->dataSet()->nCols()) {
- peak->setIntegrationFlag(RejectionFlag::CentreOutOfBounds, IntegratorType::Gaussian);
- return false;
+ result.integration_flag = RejectionFlag::CentreOutOfBounds;
+ return result;
}
if (x0(1) < 0 || x0(1) >= peak->dataSet()->nRows()) {
- peak->setIntegrationFlag(RejectionFlag::CentreOutOfBounds, IntegratorType::Gaussian);
- return false;
+ result.integration_flag = RejectionFlag::CentreOutOfBounds;
+ return result;
}
if (x0(2) < 0 || x0(2) >= peak->dataSet()->nFrames()) {
- peak->setIntegrationFlag(RejectionFlag::CentreOutOfBounds, IntegratorType::Gaussian);
- return false;
+ result.integration_flag = RejectionFlag::CentreOutOfBounds;
+ return result;
}
// consistency check: covariance matrix should be positive definite
- Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> solver(from_cholesky(a));
+ const Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> solver(from_cholesky(a));
if (solver.eigenvalues().minCoeff() <= 0) {
- peak->setIntegrationFlag(RejectionFlag::InvalidCovariance, IntegratorType::Gaussian);
- return false;
+ result.integration_flag = RejectionFlag::InvalidCovariance;
+ return result;
}
const auto& covar = min.covariance();
- _profileBackground = {B, covar(0, 0)};
- _profileIntensity = {I, covar(1, 1)};
+ result.profile_background = {B, covar(0, 0)};
+ result.profile_intensity = {I, covar(1, 1)};
// Gets pearson coefficient of fit
double pearson;
@@ -189,14 +196,15 @@ bool GaussianIntegrator::compute(
residuals(r, B, I, x0, a, x, counts, &pearson);
if (pearson <= 0.75) {
- peak->setIntegrationFlag(RejectionFlag::BadIntegrationFit, IntegratorType::Gaussian);
- return false;
+ result.integration_flag = RejectionFlag::BadIntegrationFit;
+ return result;
}
- _sumIntensity = {};
- _sumBackground = {};
+ result.sum_background = {};
+ result.sum_intensity = {};
- peak->setShape({x0, from_cholesky(a)});
- return true;
+ result.shape = {x0, from_cholesky(a)};
+
+ return result;
}
std::vector<double> GaussianIntegrator::profile(Peak3D* peak, const IntegrationRegion& region)
diff --git a/core/integration/GaussianIntegrator.h b/core/integration/GaussianIntegrator.h
index 2ef5e3d66..32764a266 100644
--- a/core/integration/GaussianIntegrator.h
+++ b/core/integration/GaussianIntegrator.h
@@ -29,7 +29,8 @@ class GaussianIntegrator : public IIntegrator {
protected:
//! Integrate a peak
- bool compute(Peak3D* peak, ShapeModel* shape_model, const IntegrationRegion& region) override;
+ ComputeResult compute(
+ Peak3D* peak, ShapeModel* shape_model, const IntegrationRegion& region) override;
//! Returns the analytic profile computed over the given integration region
std::vector<double> profile(Peak3D* peak, const IntegrationRegion& region);
};
diff --git a/core/integration/IIntegrator.cpp b/core/integration/IIntegrator.cpp
index 8a74f1a0a..9ba640452 100644
--- a/core/integration/IIntegrator.cpp
+++ b/core/integration/IIntegrator.cpp
@@ -54,23 +54,25 @@ void IntegrationParameters::log(const Level& level) const
ohklLog(level, "max_width = ", max_width);
}
-IIntegrator::IIntegrator()
- : _sumBackground()
- , _profileBackground()
- , _meanBkgGradient()
- , _sumIntensity()
- , _profileIntensity()
- , _rockingCurve()
- , _handler(nullptr)
- , _params{}
+ComputeResult::ComputeResult()
+ : integration_flag(RejectionFlag::NotRejected)
+ , sum_intensity()
+ , profile_intensity()
+ , sum_background()
+ , profile_background()
+ , bkg_gradient()
+ , rocking_curve()
+ , integrator_type(IntegratorType::PixelSum)
+ , shape()
{
}
-IIntegrator::~IIntegrator() = default;
-
-const std::vector<Intensity>& IIntegrator::rockingCurve() const
+IIntegrator::IIntegrator()
+ : _handler(nullptr)
+ , _params{}
+ , _thread_parallel(true)
+ , _max_threads(8)
{
- return _rockingCurve;
}
void IIntegrator::integrate(
@@ -93,38 +95,40 @@ void IIntegrator::integrate(
_handler->setProgress(0);
}
- bool profile_integration = false;
+ _profile_integration = false;
if (_params.integrator_type == IntegratorType::Profile1D
|| _params.integrator_type == IntegratorType::Profile3D
- || _params.integrator_type == IntegratorType::Gaussian)
- profile_integration = true;
+ || _params.integrator_type == IntegratorType::Gaussian
+ || _params.integrator_type == IntegratorType::ISigma)
+ _profile_integration = true;
size_t idx = 0;
- int num_frames_done = 0;
+ _n_frames_done = 0;
std::map<Peak3D*, std::unique_ptr<IntegrationRegion>> regions;
std::map<Peak3D*, bool> integrated;
- double peak_end, bkg_begin, bkg_end;
if (_params.region_type == ohkl::RegionType::VariableEllipsoid) {
- peak_end = _params.peak_end;
- bkg_begin = _params.bkg_begin;
- bkg_end = _params.bkg_end;
+ _peak_end = _params.peak_end;
+ _bkg_begin = _params.bkg_begin;
+ _bkg_end = _params.bkg_end;
} else {
- peak_end = _params.fixed_peak_end;
- bkg_begin = _params.fixed_bkg_begin;
- bkg_end = _params.fixed_bkg_end;
+ _peak_end = _params.fixed_peak_end;
+ _bkg_begin = _params.fixed_bkg_begin;
+ _bkg_end = _params.fixed_bkg_end;
}
for (auto peak : peaks) {
- if (!peak->enabled() && peak->isRejectedFor(RejectionFlag::Extinct))
- continue;
+ if (peak->isRejectedFor(RejectionFlag::Extinct)) {
+ if (!peak->enabled())
+ continue;
+ }
regions.emplace(std::make_pair(
peak,
std::make_unique<IntegrationRegion>(
- peak, peak_end, bkg_begin, bkg_end, _params.region_type)));
+ peak, _peak_end, _bkg_begin, _bkg_end, _params.region_type)));
integrated.emplace(std::make_pair(peak, false));
// ignore partials
@@ -154,7 +158,7 @@ void IIntegrator::integrate(
removeOverlaps(regions);
ohklLog(Level::Debug, "IIntegrator::integrate: frames loop");
- int nfailures = 0;
+ _n_failures = 0;
for (idx = 0; idx < data->nFrames(); ++idx) {
Eigen::MatrixXd current_frame, gradient;
Eigen::MatrixXi mask;
@@ -167,42 +171,43 @@ void IIntegrator::integrate(
for (auto peak : peaks) {
assert(peak != nullptr);
- auto* current_peak = regions.at(peak).get();
- // assert(current_peak != regions.end());
- current_peak->updateMask(mask, idx);
+ auto* current_region = regions.at(peak).get();
+ current_region->updateMask(mask, idx);
}
for (auto peak : peaks) {
- auto* current_peak = regions.at(peak).get();
+ auto* current_region = regions.at(peak).get();
// Check whether the peak intersects a mask
if (_params.skip_masked) {
Ellipsoid shape = peak->shape();
- shape.scale(bkg_end);
+ shape.scale(_bkg_end);
for (const auto* mask : data->masks()) {
if (mask->collide(shape)) {
- peak->setRejectionFlag(RejectionFlag::Masked);
- ++nfailures;
+ peak->setMasked();
continue;
}
}
}
// Check for saturated pixels
- const auto& counts = current_peak->peakData().counts();
+ const auto& counts = current_region->peakData().counts();
double max = 0;
if (!counts.empty()) // std::max on empty vector segfaults under MacOS
max = *std::max_element(counts.begin(), counts.end());
- bool saturated = _params.discard_saturated && (max > _params.max_counts);
+ const bool saturated = _params.discard_saturated && (max > _params.max_counts);
- bool result;
+ bool result = false;
if (_params.use_gradient)
- result = current_peak->advanceFrame(current_frame, mask, idx, gradient);
+ result = current_region->advanceFrame(current_frame, mask, idx, &gradient);
else
- result = current_peak->advanceFrame(current_frame, mask, idx);
-
- // Skip strong and low resolution peaks during profile integration
- if (profile_integration && !reintegrate(peak)) {
- result = false;
- integrated[peak] = true;
+ result = current_region->advanceFrame(current_frame, mask, idx);
+
+ bool reintegrate = true;
+ if (_profile_integration) {
+ if (_params.use_max_strength && peak->sumIntensity().strength() >
+ _params.max_strength)
+ reintegrate = false;
+ if (_params.use_max_d && peak->d() > _params.max_d)
+ reintegrate = false;
}
// this allows for partials at end of data
@@ -210,34 +215,38 @@ void IIntegrator::integrate(
// done reading peak data
if (result && !integrated[peak]) {
- current_peak->peakData().standardizeCoords();
- if (compute(peak, shape_model, *current_peak)) {
- peak->updateIntegration(
- _rockingCurve, _sumBackground, _profileBackground, _meanBkgGradient,
- _sumIntensity, _profileIntensity, _params.peak_end, _params.bkg_begin,
- _params.bkg_end, _params.region_type);
- if (saturated)
- peak->setIntegrationFlag(
- RejectionFlag::SaturatedPixel, _params.integrator_type);
- } else {
- ++nfailures;
- // This is a fallback. The RejectionFlag should have been set by this point.
- peak->setIntegrationFlag(
- RejectionFlag::IntegrationFailure, _params.integrator_type);
+ ComputeResult compute_result;
+ if (!reintegrate) { // do not reintegrate using profile fitting
+ compute_result.integrator_type = _params.integrator_type;
+ compute_result.rocking_curve = peak->rockingCurve();
+ compute_result.profile_intensity = peak->sumIntensity();
+ compute_result.profile_background = peak->sumBackground();
+ compute_result.integration_flag = peak->sumRejectionFlag();
+ integrated[peak] = true;
+ } else { // do the integration
+ current_region->peakData().standardizeCoords();
+ compute_result = compute(peak, shape_model, *current_region);
}
- // free memory (important!!)
- current_peak->reset();
+
+ if (saturated)
+ compute_result.integration_flag = RejectionFlag::SaturatedPixel;
+ if (compute_result.integration_flag != RejectionFlag::NotRejected)
+ ++_n_failures;
integrated[peak] = true;
+ peak->updateIntegration(
+ compute_result, _params.peak_end, _params.bkg_begin, _params.bkg_end,
+ _params.region_type);
+ current_region->reset();
}
}
if (_handler) {
- ++num_frames_done;
- double progress = num_frames_done * 100.0 / data->nFrames();
+ ++_n_frames_done;
+ double progress = _n_frames_done * 100.0 / data->nFrames();
_handler->setProgress(progress);
}
}
- ohklLog(Level::Info, "IIntegrator::integrate: end; ", nfailures, " failures");
+ ohklLog(Level::Info, "IIntegrator::integrate: end; ", _n_failures, " failures");
}
void IIntegrator::setHandler(sptrProgressHandler handler)
@@ -295,28 +304,4 @@ void IIntegrator::removeOverlaps(
ohklLog(Level::Info, "IIntegrator::removeOverlaps: ", nrejected, " overlapping peaks rejected");
}
-bool IIntegrator::reintegrate(Peak3D* peak)
-{
- bool reintegrate = true;
- if (_params.use_max_strength) // Skip strong peaks
- if (peak->sumIntensity().strength() > _params.max_strength)
- reintegrate = false;
- if (_params.use_max_d) // Skip low resolution peaks
- if (peak->d() > _params.max_d)
- reintegrate = false;
-
- if (!reintegrate) {
- if (peak->sumRejectionFlag() == RejectionFlag::NotRejected) {
- peak->updateIntegration(
- _rockingCurve, peak->sumBackground(), peak->sumBackground(),
- peak->meanBkgGradient(), peak->sumIntensity(), peak->sumIntensity(),
- _params.peak_end, _params.bkg_begin, _params.bkg_end, _params.region_type);
- } else {
- peak->setIntegrationFlag(peak->sumRejectionFlag(), IntegratorType::Profile3D);
- }
- }
-
- return reintegrate;
-}
-
} // namespace ohkl
diff --git a/core/integration/IIntegrator.h b/core/integration/IIntegrator.h
index 242ad0dd4..c5ad87964 100644
--- a/core/integration/IIntegrator.h
+++ b/core/integration/IIntegrator.h
@@ -20,15 +20,47 @@
#include "core/peak/IntegrationRegion.h"
#include "core/peak/Intensity.h"
+#include <optional>
+
namespace ohkl {
enum class Level;
+enum class RejectionFlag;
-enum class IntegratorType { PixelSum = 0, Gaussian, ISigma, Profile1D, Profile3D, Count };
+enum class IntegratorType { PixelSum = 0, Gaussian, ISigma, Profile1D, Profile3D, Shape, Count };
/*! \addtogroup python_api
* @{*/
+ //! Return value of all integrator compute methods. Contains any mutable data from integration.
+struct ComputeResult {
+ //! Construct an invalid result by default
+ ComputeResult();
+ ComputeResult(const ComputeResult& other) = default;
+ ComputeResult& operator=(const ComputeResult& other) = default;
+
+ //! Rejection flag for integration
+ RejectionFlag integration_flag;
+ //! Pixel sum integrated intensity after background subtraction
+ Intensity sum_intensity;
+ //! Profile integrated intensity after background correction
+ Intensity profile_intensity;
+ //! Mean local pixel sum background of peak. The uncertainty is the uncertainty of the
+ //! _estimate_ of the background.
+ Intensity sum_background;
+ //! Profile integrated background of peak. The uncertainty is the uncertainty of the
+ //! _estimate_ of the background.
+ Intensity profile_background;
+ //! Mean gradient of background (Gaussian statistics)
+ Intensity bkg_gradient;
+ //! The rocking curve of the peak.
+ std::vector<Intensity> rocking_curve;
+ //! Type of integrator used
+ IntegratorType integrator_type;
+ //! Shape of peak when centre/covariance are changed during integration
+ std::optional<Ellipsoid> shape;
+};
+
//! Structure containing parameters for all integrators
struct IntegrationParameters {
//! End of peak region for RegionType::VariableEllipsoid (sigmas)
@@ -90,46 +122,37 @@ class ShapeModel;
class IIntegrator {
public:
IIntegrator();
- virtual ~IIntegrator();
- //! Integrate all peaks in the list which are contained in the specified data set.
+ virtual ~IIntegrator() = default;
+ //! Integrate a vector of peaks
void integrate(std::vector<Peak3D*> peaks, ShapeModel* shape_model, sptrDataSet data);
- const std::vector<Intensity>& rockingCurve() const;
//! Set the progress handler.
void setHandler(sptrProgressHandler handler);
+ //! Assign a parameter set to the integrator
+ void setParameters(const IntegrationParameters& params);
+ //! Remove overlapping peak intensity regions
+ void removeOverlaps(const std::map<Peak3D*, std::unique_ptr<IntegrationRegion>>& regions);
protected:
- //! If IntegratorType != PixelSum, check whether we should profile integrate this peak
- bool reintegrate(Peak3D* peak);
-
- //! Mean local pixel sum background of peak. The uncertainty is the uncertainty of the
- //! _estimate_ of the background.
- Intensity _sumBackground;
- //! Profile integrated background of peak. The uncertainty is the uncertainty of the
- //! _estimate_ of the background.
- Intensity _profileBackground;
- //! Mean gradient of background (Gaussian statistics)
- Intensity _meanBkgGradient;
- //! Pixel sum integrated intensity, after background correction.
- Intensity _sumIntensity;
- //! Profile integrated intensity, after background correction.
- Intensity _profileIntensity;
- //! The rocking curve of the peak.
- std::vector<Intensity> _rockingCurve;
//! Optional pointer to progress handler.
sptrProgressHandler _handler;
//! Container for user-defined integration parameters
IntegrationParameters _params;
- public:
- //! Assign a parameter set to the integrator
- void setParameters(const IntegrationParameters& params);
+ sptrDataSet _data;
+ std::atomic_int _n_failures;
+ std::atomic_int _n_frames_done;
+ std::atomic_int _n_peaks_done;
+ bool _profile_integration;
+ double _peak_end;
+ double _bkg_begin;
+ double _bkg_end;
- //! Remove overlapping peak intensity regions
- void removeOverlaps(const std::map<Peak3D*, std::unique_ptr<IntegrationRegion>>& regions);
+ bool _thread_parallel;
+ unsigned int _max_threads;
private:
//! Compute the integrated intensity of the peak given the integration region.
- virtual bool compute(
+ virtual ComputeResult compute(
Peak3D* peak, ShapeModel* shape_model, const IntegrationRegion& region) = 0;
};
diff --git a/core/integration/ISigmaIntegrator.cpp b/core/integration/ISigmaIntegrator.cpp
index 11322aa82..e75ac4950 100644
--- a/core/integration/ISigmaIntegrator.cpp
+++ b/core/integration/ISigmaIntegrator.cpp
@@ -25,33 +25,33 @@ namespace ohkl {
ISigmaIntegrator::ISigmaIntegrator() : PixelSumIntegrator() { }
-bool ISigmaIntegrator::compute(
+ComputeResult ISigmaIntegrator::compute(
Peak3D* peak, ShapeModel* shape_model, const IntegrationRegion& region)
{
+ ComputeResult result;
+ result.integrator_type = IntegratorType::ISigma;
+
if (!shape_model) {
- peak->setIntegrationFlag(RejectionFlag::NoShapeModel, IntegratorType::ISigma);
- return false;
+ result.integration_flag = RejectionFlag::NoShapeModel;
+ return result;
}
- if (!peak)
- return false;
-
// first get mean background
- PixelSumIntegrator::compute(peak, shape_model, region);
- const double mean_bkg = _sumBackground.value();
- const double var_bkg = _sumBackground.variance();
+ const ComputeResult pxsum_result = PixelSumIntegrator::compute(peak, shape_model, region);
+ const double mean_bkg = pxsum_result.sum_background.value();
+ const double var_bkg = pxsum_result.sum_background.variance();
const auto& events = region.peakData().events();
const auto& counts = region.peakData().counts();
// TODO: should this be hard-coded??
if (events.size() < 29) {
- peak->setIntegrationFlag(RejectionFlag::TooFewPoints, IntegratorType::ISigma);
- return false;
+ result.integration_flag = RejectionFlag::TooFewPoints;
+ return result;
}
- Eigen::Vector3d c = peak->shape().center();
- Eigen::Matrix3d A = peak->shape().metric();
+ const Eigen::Vector3d c = peak->shape().center();
+ const Eigen::Matrix3d A = peak->shape().metric();
Profile1D profile;
std::vector<Intensity> mean_profile = shape_model->meanProfile1D(DetectorEvent(c));
@@ -87,26 +87,26 @@ bool ISigmaIntegrator::compute(
// something went wrong (nans?)
if (best_idx < 0) {
- peak->setIntegrationFlag(RejectionFlag::NoISigmaMinimum, IntegratorType::ISigma);
- return false;
+ result.integration_flag = RejectionFlag::NoISigmaMinimum;
+ return result;
}
const double M = profile.counts()[best_idx];
const int n = profile.npoints()[best_idx];
- _profileIntensity = Intensity(M - n * mean_bkg, M + n * n * var_bkg);
- _profileIntensity = _profileIntensity / mean_profile[best_idx];
+ result.profile_intensity = Intensity(M - n * mean_bkg, M + n * n * var_bkg);
+ result.profile_intensity = result.profile_intensity / mean_profile[best_idx];
- double sigma = _profileIntensity.sigma();
+ const double sigma = result.profile_intensity.sigma();
if (std::isnan(sigma) && sigma > 0) {
- peak->setIntegrationFlag(RejectionFlag::InvalidSigma, IntegratorType::ISigma);
- return false;
+ result.integration_flag = RejectionFlag::InvalidSigma;
+ return result;
}
- _sumIntensity = {};
- _sumBackground = {};
+ result.sum_intensity = {};
+ result.sum_background = {};
- return true;
+ return result;
}
} // namespace ohkl
diff --git a/core/integration/ISigmaIntegrator.h b/core/integration/ISigmaIntegrator.h
index 799133986..397fdc319 100644
--- a/core/integration/ISigmaIntegrator.h
+++ b/core/integration/ISigmaIntegrator.h
@@ -42,7 +42,8 @@ class ISigmaIntegrator : public PixelSumIntegrator {
ISigmaIntegrator();
protected:
- bool compute(Peak3D* peak, ShapeModel* shape_model, const IntegrationRegion& region) override;
+ ComputeResult compute(
+ Peak3D* peak, ShapeModel* shape_model, const IntegrationRegion& region) override;
};
/*! @}*/
diff --git a/core/integration/PixelSumIntegrator.cpp b/core/integration/PixelSumIntegrator.cpp
index ed5aae273..3cffcfaf1 100644
--- a/core/integration/PixelSumIntegrator.cpp
+++ b/core/integration/PixelSumIntegrator.cpp
@@ -13,6 +13,7 @@
// ***********************************************************************************************
#include "core/integration/PixelSumIntegrator.h"
+
#include "base/geometry/Ellipsoid.h"
#include "base/utils/Logger.h"
#include "core/data/DataSet.h"
@@ -91,31 +92,34 @@ PixelSumIntegrator::PixelSumIntegrator() : IIntegrator()
_params.fit_cov = true;
}
-bool PixelSumIntegrator::compute(Peak3D* peak, ShapeModel*, const IntegrationRegion& region)
+ComputeResult PixelSumIntegrator::compute(Peak3D* peak, ShapeModel*, const IntegrationRegion& region)
{
+ ComputeResult result;
+ result.integrator_type = IntegratorType::PixelSum;
+
const auto& events = region.peakData().events();
const auto& counts = region.peakData().counts();
if (events.empty()) {
- peak->setIntegrationFlag(RejectionFlag::TooFewPoints, IntegratorType::PixelSum);
- return false;
+ result.integration_flag = RejectionFlag::TooFewPoints;
+ return result;
}
auto [meanBackground, bkgGradient] = compute_background(region, _params.use_gradient);
if (!meanBackground.isValid()) {
- peak->setIntegrationFlag(RejectionFlag::TooFewPoints, IntegratorType::PixelSum);
- return false;
+ result.integration_flag = RejectionFlag::TooFewPoints;
+ return result;
}
- _sumBackground = meanBackground;
- _meanBkgGradient = bkgGradient;
+ result.sum_background = meanBackground;
+ result.bkg_gradient = bkgGradient;
- PeakCoordinateSystem frame(peak);
+ const PeakCoordinateSystem frame(peak);
double sum_peak = 0.0;
- const double mean_bkg = _sumBackground.value();
+ const double mean_bkg = result.sum_background.value();
// note that this is the std of the _estimate_ of the background
// should be approximately mean_bkg / num_bkg for Poisson statistics
- const double std_bkg = _sumBackground.sigma();
+ const double std_bkg = result.sum_background.sigma();
size_t npeak = 0.0;
Blob3D blob;
@@ -126,9 +130,8 @@ bool PixelSumIntegrator::compute(Peak3D* peak, ShapeModel*, const IntegrationReg
const auto& ev = events[i];
auto ev_type = region.classify(ev);
- if (ev_type == IntegrationRegion::EventType::BACKGROUND) {
+ if (ev_type == IntegrationRegion::EventType::BACKGROUND)
continue;
- }
if (ev_type == IntegrationRegion::EventType::PEAK) {
sum_peak += counts[i];
@@ -140,11 +143,11 @@ bool PixelSumIntegrator::compute(Peak3D* peak, ShapeModel*, const IntegrationReg
}
}
- sum_peak -= npeak * _sumBackground.value();
+ sum_peak -= npeak * result.sum_background.value();
// TODO: ERROR ESTIMATE!!
// This INCORRECTLY assumes Poisson statistics (no gain or baseline)
- _sumIntensity =
+ result.sum_intensity =
Intensity(sum_peak, sum_peak + npeak * mean_bkg + npeak * npeak * std_bkg * std_bkg);
// TODO: compute rocking curve
@@ -164,8 +167,8 @@ bool PixelSumIntegrator::compute(Peak3D* peak, ShapeModel*, const IntegrationReg
if (blob.isValid()) {
center = blob.center();
} else {
- peak->setIntegrationFlag(RejectionFlag::InvalidCentroid, IntegratorType::PixelSum);
- return false;
+ result.integration_flag = RejectionFlag::InvalidCentroid;
+ return result;
}
} else {
center = peak->shape().center();
@@ -175,8 +178,8 @@ bool PixelSumIntegrator::compute(Peak3D* peak, ShapeModel*, const IntegrationReg
if (blob.isValid()) {
cov = blob.covariance();
} else {
- peak->setIntegrationFlag(RejectionFlag::InvalidCentroid, IntegratorType::PixelSum);
- return false;
+ result.integration_flag = RejectionFlag::InvalidCentroid;
+ return result;
}
} else {
cov = peak->shape().inverseMetric();
@@ -184,37 +187,37 @@ bool PixelSumIntegrator::compute(Peak3D* peak, ShapeModel*, const IntegrationReg
// center of mass is consistent
if (std::isnan(center.norm())) {
- peak->setIntegrationFlag(RejectionFlag::InvalidCentroid, IntegratorType::PixelSum);
- return false;
+ result.integration_flag = RejectionFlag::InvalidCentroid;
+ return result;
}
if (!peak->shape().isInside(center)) {
- peak->setIntegrationFlag(RejectionFlag::InvalidCentroid, IntegratorType::PixelSum);
- return false;
+ result.integration_flag = RejectionFlag::InvalidCentroid;
+ return result;
}
- Eigen::Matrix3d A0 = peak->shape().metric();
- Eigen::Matrix3d A1 = cov.inverse();
+ const Eigen::Matrix3d A0 = peak->shape().metric();
+ const Eigen::Matrix3d A1 = cov.inverse();
const double dA = (A1 - A0).norm() / A0.norm();
// check that the covariance is consistent
if (!(dA < 2.0)) {
- peak->setIntegrationFlag(RejectionFlag::InvalidCovariance, IntegratorType::PixelSum);
- return false;
+ result.integration_flag = RejectionFlag::InvalidCovariance;
+ return result;
}
// shape is not too small or too large
- Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> solver(cov);
+ const Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> solver(cov);
const auto& w = solver.eigenvalues();
if (w.minCoeff() < 0.1 || w.maxCoeff() > 100) {
- peak->setIntegrationFlag(RejectionFlag::InvalidShape, IntegratorType::PixelSum);
- return false;
+ result.integration_flag = RejectionFlag::InvalidShape;
+ return result;
}
- peak->setShape(Ellipsoid(center, A1));
+ result.shape = Ellipsoid(center, A1);
- size_t nframes = size_t(f_max - f_min) + 1;
- _rockingCurve.resize(nframes);
+ const size_t nframes = size_t(f_max - f_min) + 1;
+ result.rocking_curve.resize(nframes);
std::vector<double> intensity_per_frame(nframes, 0.0);
std::vector<double> n_peak_points_per_frame(nframes, 0.0);
@@ -236,12 +239,12 @@ bool PixelSumIntegrator::compute(Peak3D* peak, ShapeModel*, const IntegrationReg
for (int i = 0; i < nframes; ++i) {
const double corrected_intensity =
intensity_per_frame[i] - n_peak_points_per_frame[i] * mean_bkg;
- _rockingCurve[i] = Intensity(corrected_intensity, sqrt(corrected_intensity));
+ result.rocking_curve[i] = Intensity(corrected_intensity, sqrt(corrected_intensity));
}
- _profileIntensity = {};
- _profileBackground = {};
+ result.profile_intensity = {};
+ result.profile_background = {};
- return true;
+ return result;
}
} // namespace ohkl
diff --git a/core/integration/PixelSumIntegrator.h b/core/integration/PixelSumIntegrator.h
index 60139d9a9..d42e659d3 100644
--- a/core/integration/PixelSumIntegrator.h
+++ b/core/integration/PixelSumIntegrator.h
@@ -33,7 +33,7 @@ class PixelSumIntegrator : public IIntegrator {
protected:
//! Integrate a peak
- bool compute(Peak3D* peak, ShapeModel*, const IntegrationRegion& region) override;
+ ComputeResult compute(Peak3D* peak, ShapeModel*, const IntegrationRegion& region) override;
};
/*! @}*/
diff --git a/core/integration/Profile1DIntegrator.cpp b/core/integration/Profile1DIntegrator.cpp
index feb17e6eb..26fee4cc5 100644
--- a/core/integration/Profile1DIntegrator.cpp
+++ b/core/integration/Profile1DIntegrator.cpp
@@ -21,13 +21,11 @@
#include "core/peak/PeakCoordinateSystem.h"
#include "core/shape/ShapeModel.h"
-namespace ohkl {
-
-Profile1DIntegrator::Profile1DIntegrator() : IIntegrator() { }
-
-static void updateFit(
- Intensity& I, Intensity& B, const std::vector<double>& dp, const std::vector<double>& dM,
- const std::vector<int>& dn)
+namespace {
+
+void updateFit(
+ ohkl::Intensity& I, ohkl::Intensity& B, const std::vector<double>& dp,
+ const std::vector<double>& dM, const std::vector<int>& dn)
{
Eigen::Matrix2d A;
A.setZero();
@@ -72,36 +70,42 @@ static void updateFit(
// Note: this error estimate assumes the variances are correct (i.e., gain and
// baseline accounted for)
- B = Intensity(new_B, cov(0, 0));
- I = Intensity(new_I, cov(1, 1));
+ B = ohkl::Intensity(new_B, cov(0, 0));
+ I = ohkl::Intensity(new_I, cov(1, 1));
}
-bool Profile1DIntegrator::compute(
+} // namespace
+
+namespace ohkl {
+
+Profile1DIntegrator::Profile1DIntegrator() : IIntegrator() { }
+
+ComputeResult Profile1DIntegrator::compute(
Peak3D* peak, ShapeModel* shape_model, const IntegrationRegion& region)
{
+ ComputeResult result;
+ result.integrator_type = IntegratorType::Profile1D;
+
if (!shape_model) {
- peak->setIntegrationFlag(RejectionFlag::NoShapeModel, IntegratorType::Profile1D);
- return false;
+ result.integration_flag = RejectionFlag::NoShapeModel;
+ return result;
}
- if (!peak)
- return false;
-
const auto& events = region.peakData().events();
const auto& counts = region.peakData().counts();
// TODO: should this be hard-coded??
if (events.size() < 29) {
- peak->setIntegrationFlag(RejectionFlag::TooFewPoints, IntegratorType::Profile1D);
- return false;
+ result.integration_flag = RejectionFlag::TooFewPoints;
+ return result;
}
- Eigen::Vector3d c = peak->shape().center();
- Eigen::Matrix3d A = peak->shape().metric();
+ const Eigen::Vector3d c = peak->shape().center();
+ const Eigen::Matrix3d A = peak->shape().metric();
Profile1D profile(0.0, region.peakEnd());
- std::vector<Intensity> mean_profile = shape_model->meanProfile1D(DetectorEvent(c));
+ const std::vector<Intensity> mean_profile = shape_model->meanProfile1D(DetectorEvent(c));
// construct the observed profile
for (size_t i = 0; i < events.size(); ++i) {
@@ -134,22 +138,22 @@ bool Profile1DIntegrator::compute(
for (auto i = 0; i < 10 && I.value() > 0; ++i)
updateFit(I, B, dp, dm, dn);
- double sigma = I.sigma();
+ const double sigma = I.sigma();
if (std::isnan(sigma) || sigma <= 0.0) {
- peak->setIntegrationFlag(RejectionFlag::InvalidSigma, IntegratorType::Profile1D);
- return false;
+ result.integration_flag = RejectionFlag::InvalidSigma;
+ return result;
}
- _profileIntensity = I;
- _profileBackground = B;
+ result.profile_intensity = I;
+ result.profile_background = B;
- _sumIntensity = {};
- _sumBackground = {};
+ result.sum_intensity = {};
+ result.sum_background = {};
// TODO: rocking curve!
- return true;
+ return result;
}
} // namespace ohkl
diff --git a/core/integration/Profile1DIntegrator.h b/core/integration/Profile1DIntegrator.h
index 306e4b8bc..2a503a2cb 100644
--- a/core/integration/Profile1DIntegrator.h
+++ b/core/integration/Profile1DIntegrator.h
@@ -41,7 +41,8 @@ class Profile1DIntegrator : public IIntegrator {
protected:
//! Integrate a peak
- bool compute(Peak3D* peak, ShapeModel* shape_model, const IntegrationRegion& region) override;
+ ComputeResult compute(
+ Peak3D* peak, ShapeModel* shape_model, const IntegrationRegion& region) override;
};
/*! @}*/
diff --git a/core/integration/Profile3DIntegrator.cpp b/core/integration/Profile3DIntegrator.cpp
index 4a5c57ef2..ff2c616f1 100644
--- a/core/integration/Profile3DIntegrator.cpp
+++ b/core/integration/Profile3DIntegrator.cpp
@@ -18,10 +18,10 @@
#include "core/peak/Peak3D.h"
#include "core/shape/ShapeModel.h"
-namespace ohkl {
+namespace {
-static void updateFit(
- Intensity& I, Intensity& B, const std::vector<double>& profile,
+void updateFit(
+ ohkl::Intensity& I, ohkl::Intensity& B, const std::vector<double>& profile,
const std::vector<double>& counts)
{
Eigen::Matrix2d A;
@@ -43,7 +43,7 @@ static void updateFit(
b(1) += M * p / var;
}
- Eigen::Matrix2d AI = A.inverse();
+ const Eigen::Matrix2d AI = A.inverse();
const Eigen::Vector2d& x = AI * b;
const double new_B = x(0);
@@ -54,33 +54,37 @@ static void updateFit(
// Note: this error estimate assumes the variances are correct (i.e., gain and
// baseline accounted for)
- B = Intensity(new_B, cov(0, 0));
- I = Intensity(new_I, cov(1, 1));
+ B = ohkl::Intensity(new_B, cov(0, 0));
+ I = ohkl::Intensity(new_I, cov(1, 1));
+}
+
}
-bool Profile3DIntegrator::compute(
+namespace ohkl {
+
+ComputeResult Profile3DIntegrator::compute(
Peak3D* peak, ShapeModel* shape_model, const IntegrationRegion& region)
{
+ ComputeResult result;
+ result.integrator_type = IntegratorType::Profile3D;
+
if (!shape_model) {
- peak->setIntegrationFlag(RejectionFlag::NoShapeModel, IntegratorType::Profile3D);
- return false;
+ result.integration_flag = RejectionFlag::NoShapeModel;
+ return result;
}
- if (!peak)
- return false;
-
const auto& events = region.peakData().events();
const auto& counts = region.peakData().counts();
// TODO: should this be hard-coded??
if (events.size() < 29) {
- peak->setIntegrationFlag(RejectionFlag::TooFewPoints, IntegratorType::Profile3D);
- return false;
+ result.integration_flag = RejectionFlag::TooFewPoints;
+ return result;
}
// dummy value for initial guess
- _profileBackground = Intensity(1.0, 1.0);
- _profileIntensity = Intensity(0.0, 0.0);
+ result.profile_background = Intensity(1.0, 1.0);
+ result.profile_intensity = Intensity(0.0, 0.0);
std::vector<double> profile;
std::vector<double> obs_counts;
@@ -90,11 +94,9 @@ bool Profile3DIntegrator::compute(
const double tolerance = 1e-5;
- DetectorEvent event(peak->shape().center());
-
- Profile3D model_profile = shape_model->meanProfile(event);
-
- PeakCoordinateSystem coord(peak);
+ const DetectorEvent event(peak->shape().center());
+ const Profile3D model_profile = shape_model->meanProfile(event);
+ const PeakCoordinateSystem coord(peak);
// evaluate the model profile at the given events
for (int i = 0; i < events.size(); ++i) {
@@ -116,13 +118,13 @@ bool Profile3DIntegrator::compute(
// todo: stopping criterion
for (auto i = 0; i < 20; ++i) {
- Intensity old_intensity = _profileIntensity;
- const double I0 = _profileIntensity.value();
- updateFit(_profileIntensity, _profileBackground, profile, obs_counts);
- const double I1 = _profileIntensity.value();
+ const Intensity old_intensity = result.profile_intensity;
+ const double I0 = result.profile_intensity.value();
+ updateFit(result.profile_intensity, result.profile_background, profile, obs_counts);
+ const double I1 = result.profile_intensity.value();
- if (std::isnan(I1) || std::isnan(_profileBackground.value())) {
- _profileIntensity = old_intensity;
+ if (std::isnan(I1) || std::isnan(result.profile_background.value())) {
+ result.profile_intensity = old_intensity;
break;
}
@@ -130,16 +132,16 @@ bool Profile3DIntegrator::compute(
break;
}
- double sigma = _profileIntensity.sigma();
+ const double sigma = result.profile_intensity.sigma();
if (std::isnan(sigma) && sigma > 0) {
- peak->setIntegrationFlag(RejectionFlag::InvalidSigma, IntegratorType::Profile3D);
- return false;
+ result.integration_flag = RejectionFlag::InvalidSigma;
+ return result;
}
- _sumIntensity = {};
- _sumBackground = {};
+ result.sum_intensity = {};
+ result.sum_background = {};
- return true;
+ return result;
}
} // namespace ohkl
diff --git a/core/integration/Profile3DIntegrator.h b/core/integration/Profile3DIntegrator.h
index 04a9eede9..bae5085f6 100644
--- a/core/integration/Profile3DIntegrator.h
+++ b/core/integration/Profile3DIntegrator.h
@@ -54,7 +54,8 @@ class Profile3DIntegrator : public IIntegrator {
protected:
//! Do the integration
- bool compute(Peak3D* peak, ShapeModel* shape_model, const IntegrationRegion& region) override;
+ ComputeResult compute(
+ Peak3D* peak, ShapeModel* shape_model, const IntegrationRegion& region) override;
};
/*! @}*/
diff --git a/core/integration/ShapeIntegrator.cpp b/core/integration/ShapeIntegrator.cpp
index e3bc6a1da..daeb331ba 100644
--- a/core/integration/ShapeIntegrator.cpp
+++ b/core/integration/ShapeIntegrator.cpp
@@ -15,6 +15,7 @@
#include "core/integration/ShapeIntegrator.h"
#include "base/geometry/Ellipsoid.h"
#include "core/data/DataSet.h"
+#include "core/integration/IIntegrator.h"
#include "core/peak/Intensity.h"
#include "core/peak/Peak3D.h"
#include "core/peak/PeakCoordinateSystem.h"
@@ -36,34 +37,32 @@ ShapeIntegrator::ShapeIntegrator(
{
}
-bool ShapeIntegrator::compute(
+ComputeResult ShapeIntegrator::compute(
Peak3D* peak, ShapeModel* shape_model, const IntegrationRegion& region)
{
+ ComputeResult result;
+ result.integrator_type = IntegratorType::Shape;
+
const UnitCell* uc = peak->unitCell();
auto data = peak->dataSet();
if (!uc) {
- peak->setIntegrationFlag(RejectionFlag::NoUnitCell, IntegratorType::PixelSum);
- return false;
- }
-
- if (!data) {
- peak->setIntegrationFlag(RejectionFlag::NoDataSet, IntegratorType::PixelSum);
- return false;
+ result.integration_flag = RejectionFlag::NoUnitCell;
+ return result;
}
- PixelSumIntegrator::compute(peak, shape_model, region);
+ const ComputeResult pxsum_result = PixelSumIntegrator::compute(peak, shape_model, region);
- const double mean_bkg = _sumBackground.value();
+ const double mean_bkg = pxsum_result.sum_background.value();
const auto& events = region.peakData().events();
const auto& counts = region.peakData().counts();
Profile3D profile(_aabb, _nx, _ny, _nz);
// todo: don't use default constructor!
- Profile1D integrated_profile(_sumBackground, region.peakEnd());
- PeakCoordinateSystem frame(peak);
+ Profile1D integrated_profile(pxsum_result.sum_background, region.peakEnd());
+ const PeakCoordinateSystem frame(peak);
- Ellipsoid e = peak->shape();
+ const Ellipsoid e = peak->shape();
for (size_t i = 0; i < events.size(); ++i) {
const auto& ev = events[i];
@@ -87,7 +86,7 @@ bool ShapeIntegrator::compute(
}
if (profile.normalize())
_shape_model->addPeak(peak, std::move(profile), std::move(integrated_profile));
- return true;
+ return result;
}
const ShapeModel* ShapeIntegrator::shapeModel() const
diff --git a/core/integration/ShapeIntegrator.h b/core/integration/ShapeIntegrator.h
index f1e79d846..ea7b07d73 100644
--- a/core/integration/ShapeIntegrator.h
+++ b/core/integration/ShapeIntegrator.h
@@ -32,7 +32,8 @@ class ShapeIntegrator : public PixelSumIntegrator {
protected:
//! Integrate a peak
- bool compute(Peak3D* peak, ShapeModel* shape_model, const IntegrationRegion& region) override;
+ ComputeResult compute(
+ Peak3D* peak, ShapeModel* shape_model, const IntegrationRegion& region) override;
//! Returns the collection of cached peak shapes
const ShapeModel* shapeModel() const;
diff --git a/core/peak/IntegrationRegion.cpp b/core/peak/IntegrationRegion.cpp
index efd986e47..83d0ca5a8 100644
--- a/core/peak/IntegrationRegion.cpp
+++ b/core/peak/IntegrationRegion.cpp
@@ -19,6 +19,7 @@
#include "core/instrument/Diffractometer.h"
#include "core/integration/IIntegrator.h"
#include "core/peak/Peak3D.h"
+#include "core/peak/PeakData.h"
#include "tables/crystal/UnitCell.h"
namespace ohkl {
@@ -26,14 +27,14 @@ namespace ohkl {
IntegrationRegion::IntegrationRegion(
Peak3D* peak, double peak_end, double bkg_begin, double bkg_end,
RegionType region_type /* = RegionType::VariableEllipsoid */)
- : _bkgBegin(bkg_begin)
+ : _shape(peak->shape())
+ , _bkgBegin(bkg_begin)
, _bkgEnd(bkg_end)
, _data(peak)
, _regionType(region_type)
, _peak(peak)
, _valid(true)
{
- _shape = peak->shape();
if (_shape.aabb().lower().hasNaN() || _shape.aabb().upper().hasNaN()) {
peak->setIntegrationFlag(RejectionFlag::InvalidShape, IntegratorType::PixelSum);
peak->setIntegrationFlag(RejectionFlag::InvalidShape, IntegratorType::Profile3D);
@@ -132,7 +133,7 @@ void IntegrationRegion::updateMask(Eigen::MatrixXi& mask, double z) const
if (val == EventType::FORBIDDEN)
continue;
- DetectorEvent ev(x, y, z);
+ const DetectorEvent ev(x, y, z);
auto ev_type = classify(ev);
switch (ev_type) {
@@ -168,8 +169,8 @@ RegionData* IntegrationRegion::getRegion()
xmax = std::min(xmax, long(data->nCols()));
ymax = std::min(ymax, long(data->nRows()));
- int zmin = std::ceil(lower[2]);
- int zmax = std::floor(upper[2]);
+ const int zmin = std::ceil(lower[2]);
+ const int zmax = std::floor(upper[2]);
_region_data = RegionData(this, xmin, xmax, ymin, ymax, zmin, zmax);
@@ -185,7 +186,7 @@ RegionData* IntegrationRegion::getRegion()
for (auto y = ymin; y < ymax; ++y) {
EventType val = EventType::EXCLUDED;
- DetectorEvent ev(x, y, z);
+ const DetectorEvent ev(x, y, z);
auto ev_type = classify(ev);
switch (ev_type) {
@@ -237,7 +238,8 @@ IntegrationRegion::EventType IntegrationRegion::classify(const DetectorEvent& ev
}
bool IntegrationRegion::advanceFrame(
- const Eigen::MatrixXd& image, const Eigen::MatrixXi& mask, double frame)
+ const Eigen::MatrixXd& image, const Eigen::MatrixXi& mask, double frame,
+ const Eigen::MatrixXd* gradient /* = nullptr */)
{
const auto aabb = _hull.aabb();
auto lower = aabb.lower();
@@ -249,28 +251,30 @@ bool IntegrationRegion::advanceFrame(
if (frame > upper[2])
return true;
- long xmin = std::max(0L, std::lround(lower[0]));
- long ymin = std::max(0L, std::lround(lower[1]));
+ const long xmin = std::max(0L, std::lround(lower[0]));
+ const long ymin = std::max(0L, std::lround(lower[1]));
- long xmax = std::min(long(image.cols()), std::lround(upper[0]));
- long ymax = std::min(long(image.rows()), std::lround(upper[1]));
+ const long xmax = std::min(long(image.cols()), std::lround(upper[0]));
+ const long ymax = std::min(long(image.rows()), std::lround(upper[1]));
for (auto x = xmin; x < xmax; ++x) {
for (auto y = ymin; y < ymax; ++y) {
- EventType mask_type = EventType(mask(y, x));
+ const EventType mask_type = EventType(mask(y, x));
if (mask_type == EventType::FORBIDDEN)
continue;
- DetectorEvent ev(x, y, frame);
- Eigen::Vector3d p(x, y, frame);
+ const DetectorEvent ev(x, y, frame);
+ const Eigen::Vector3d p(x, y, frame);
auto event_type = classify(ev);
- if (event_type == EventType::PEAK) {
- _data.addEvent(ev, image(y, x));
- }
+ double grad = 0;
+ if (gradient)
+ grad = (*gradient)(y, x);
- if (event_type == EventType::BACKGROUND && mask_type == EventType::BACKGROUND) {
- _data.addEvent(ev, image(y, x));
- }
+ if (event_type == EventType::PEAK)
+ _data.addEvent(ev, image(y, x), grad);
+
+ if (event_type == EventType::BACKGROUND && mask_type == EventType::BACKGROUND)
+ _data.addEvent(ev, image(y, x), grad);
// check if point is in Brillouin zone (or AABB if no UC available)
// if (_hull.contains(p)) {
@@ -281,53 +285,51 @@ bool IntegrationRegion::advanceFrame(
return false;
}
-bool IntegrationRegion::advanceFrame(
+PeakData IntegrationRegion::threadSafeAdvanceFrame(
const Eigen::MatrixXd& image, const Eigen::MatrixXi& mask, double frame,
- const Eigen::MatrixXd& gradient)
+ const Eigen::MatrixXd* gradient /* = nullptr */)
{
+ PeakData peak_data(_peak);
+
const auto aabb = _hull.aabb();
auto lower = aabb.lower();
auto upper = aabb.upper();
if (frame < lower[2])
- return false;
+ return peak_data;
if (frame > upper[2])
- return true;
+ return peak_data;
- long xmin = std::max(0L, std::lround(lower[0]));
- long ymin = std::max(0L, std::lround(lower[1]));
+ const long xmin = std::max(0L, std::lround(lower[0]));
+ const long ymin = std::max(0L, std::lround(lower[1]));
- long xmax = std::min(long(image.cols()), std::lround(upper[0]));
- long ymax = std::min(long(image.rows()), std::lround(upper[1]));
+ const long xmax = std::min(long(image.cols()), std::lround(upper[0]));
+ const long ymax = std::min(long(image.rows()), std::lround(upper[1]));
for (auto x = xmin; x < xmax; ++x) {
for (auto y = ymin; y < ymax; ++y) {
- EventType mask_type = EventType(mask(y, x));
+ const EventType mask_type = EventType(mask(y, x));
if (mask_type == EventType::FORBIDDEN)
continue;
- DetectorEvent ev(x, y, frame);
- Eigen::Vector3d p(x, y, frame);
+ const DetectorEvent ev(x, y, frame);
+ const Eigen::Vector3d p(x, y, frame);
auto event_type = classify(ev);
- if (event_type == EventType::PEAK) {
- _data.addEvent(ev, image(y, x), gradient(y, x));
- }
+ double grad = 0;
+ if (gradient)
+ grad = (*gradient)(y, x);
- if (event_type == EventType::BACKGROUND && mask_type == EventType::BACKGROUND) {
- _data.addEvent(ev, image(y, x), gradient(y, x));
- }
+ if (event_type == EventType::PEAK)
+ peak_data.addEvent(ev, image(y, x), grad);
- // check if point is in Brillouin zone (or AABB if no UC available)
- // if (_hull.contains(p)) {
- // _data.addEvent(ev, image(y,x));
- //}
+ if (event_type == EventType::BACKGROUND && mask_type == EventType::BACKGROUND)
+ peak_data.addEvent(ev, image(y, x), grad);
}
}
- return false;
+ return peak_data;
}
-
void IntegrationRegion::reset()
{
_data.reset();
@@ -343,6 +345,12 @@ PeakData& IntegrationRegion::peakData()
return _data;
}
+void IntegrationRegion::appendPeakData(const PeakData& peak_data)
+{
+ if (!peak_data.empty())
+ _data.append(peak_data);
+}
+
const Ellipsoid& IntegrationRegion::shape() const
{
return _shape;
diff --git a/core/peak/IntegrationRegion.h b/core/peak/IntegrationRegion.h
index 0c4e5858b..09ba78152 100644
--- a/core/peak/IntegrationRegion.h
+++ b/core/peak/IntegrationRegion.h
@@ -57,17 +57,21 @@ class IntegrationRegion {
//! Classify a detector event (peak, background, forbidden, etc.)
EventType classify(const DetectorEvent& ev) const;
//! Update the region with the next frame
- bool advanceFrame(const Eigen::MatrixXd& image, const Eigen::MatrixXi& mask, double frame);
- //! Update the region with the next frame
bool advanceFrame(
const Eigen::MatrixXd& image, const Eigen::MatrixXi& mask, double frame,
- const Eigen::MatrixXd& gradient);
+ const Eigen::MatrixXd* gradient = nullptr);
+ //! Update the region with the next frame (thread safe)
+ PeakData threadSafeAdvanceFrame(
+ const Eigen::MatrixXd& image, const Eigen::MatrixXi& mask, double frame,
+ const Eigen::MatrixXd* gradient = nullptr);
//! Reset the integration region (i.e. free memory)
void reset();
//! Returns the underlying data stored by the region
const PeakData& peakData() const;
//! Returns the data stored by the region
PeakData& peakData();
+ //! Append PeakData to this instance
+ void appendPeakData(const PeakData& peak_data);
//! Returns the peak shape used by the region
const Ellipsoid& shape() const;
//! Returns the convex hull of the region (e.g. BrillouinZone)
diff --git a/core/peak/Peak3D.cpp b/core/peak/Peak3D.cpp
index 9aeb062b8..3da3d2519 100644
--- a/core/peak/Peak3D.cpp
+++ b/core/peak/Peak3D.cpp
@@ -61,6 +61,7 @@ const std::map<RejectionFlag, std::string> Peak3D::_rejection_map{
{RejectionFlag::NoShapeModel, "No shape model found"},
{RejectionFlag::NoISigmaMinimum, "Failed to find minimum of I/Sigma"},
{RejectionFlag::TooWide, "Bounding box spans too many images"},
+ {RejectionFlag::BadGaussianFit, "Unable to fit Gaussian to peak"},
{RejectionFlag::PredictionUpdateFailure, "Failure updating prediction post-refinement"},
{RejectionFlag::ManuallyRejected, "Manually unselected by user"},
{RejectionFlag::OutsideIndexingTol, "Outside indexing tolerance"},
@@ -110,26 +111,26 @@ Peak3D::Peak3D(sptrDataSet data, const MillerIndex& hkl)
}
Peak3D::Peak3D(std::shared_ptr<ohkl::Peak3D> peak)
+ : _sumIntensity(peak->sumIntensity())
+ , _profileIntensity(peak->profileIntensity())
+ , _sumBackground(peak->sumBackground())
+ , _profileBackground(peak->profileBackground())
+ , _meanBkgGradient(peak->meanBkgGradient())
+ , _peakEnd(peak->peakEnd())
+ , _bkgBegin(peak->bkgBegin())
+ , _bkgEnd(peak->bkgEnd())
+ , _regionType(peak->regionType())
+ , _hkl(peak->hkl())
+ , _unitCell(peak->_unitCell)
+ , _scale(peak->scale())
+ , _caught_by_filter(false)
+ , _rejected_by_filter(false)
+ , _transmission(peak->transmission())
+ , _data(peak->dataSet())
+ , _rockingCurve(peak->rockingCurve())
+
{
setShape(peak->shape());
- _sumIntensity = peak->sumIntensity();
- _profileIntensity = peak->profileIntensity();
- _sumBackground = peak->sumBackground();
- _profileBackground = peak->profileBackground();
- _meanBkgGradient = peak->meanBkgGradient();
- _peakEnd = peak->peakEnd();
- _bkgBegin = peak->bkgBegin();
- _bkgEnd = peak->bkgEnd();
- _regionType = peak->regionType();
- _unitCell = peak->_unitCell;
- _scale = peak->scale();
- _transmission = peak->transmission();
- _data = peak->dataSet();
- _rockingCurve = peak->rockingCurve();
- _hkl = peak->hkl();
-
- _caught_by_filter = false;
- _rejected_by_filter = false;
}
void Peak3D::setShape(const Ellipsoid& shape)
@@ -164,8 +165,7 @@ const UnitCell* Peak3D::unitCell() const
{
if (auto uc = _unitCell.lock())
return uc.get();
- else
- return nullptr;
+ return nullptr;
}
Intensity Peak3D::correctedIntensity(const Intensity& intensity) const
@@ -173,16 +173,16 @@ Intensity Peak3D::correctedIntensity(const Intensity& intensity) const
auto c = _shape.center();
auto state = InterpolatedState::interpolate(_data->instrumentStates(), c[2]);
if (!state.isValid()) // Interpolation error
- return Intensity();
+ return {};
auto diff = state.diffractometer();
if (diff == nullptr) {
- return Intensity();
+ return {};
}
auto detector = diff->detector();
if (detector == nullptr) {
- return Intensity();
+ return {};
}
const double lorentz = state.lorentzFactor(c[0], c[1]);
@@ -235,35 +235,40 @@ void Peak3D::reject(RejectionFlag flag)
}
void Peak3D::updateIntegration(
- const std::vector<Intensity>& rockingCurve, const Intensity& sumBkg, const Intensity& profBkg,
- const Intensity& meanBkgGradient, const Intensity& sumInt, const Intensity& profInt,
- double peakEnd, double bkgBegin, double bkgEnd, const RegionType& regionType)
-{
- _rockingCurve = rockingCurve;
- _meanBkgGradient = meanBkgGradient;
- if (sumBkg.isValid())
- _sumBackground = sumBkg;
- if (sumInt.isValid()) { // Default intensity constructor is zero, _valid = false
- if (sumInt.sigma() < _sigma2_eps) // NaN sigma handled by Intensity constructor
- setIntegrationFlag(RejectionFlag::InvalidSigma, IntegratorType::PixelSum);
- else
- _sumIntensity = sumInt;
- }
- if (profBkg.isValid())
- _profileBackground = profBkg;
- if (profInt.isValid()) {
- if (profInt.sigma() < _sigma2_eps) // NaN sigma handled by Intensity constructor
- setIntegrationFlag(RejectionFlag::InvalidSigma, IntegratorType::Profile3D);
- else
- _profileIntensity = profInt;
- }
+ const ComputeResult& result, double peakEnd, double bkgBegin, double bkgEnd,
+ const RegionType& regionType)
+{
+ if (result.integrator_type == IntegratorType::Shape)
+ return;
- //_sumIntensity = integrator.peakIntensity(); // TODO: test, reactivate ???
- //_shape = integrator.fitShape(); // TODO: test, reactivate ???
_peakEnd = peakEnd;
_bkgBegin = bkgBegin;
_bkgEnd = bkgEnd;
_regionType = regionType;
+ _rockingCurve = result.rocking_curve;
+
+ _rockingCurve = result.rocking_curve;
+ if (result.integrator_type == IntegratorType::PixelSum)
+ _meanBkgGradient = result.bkg_gradient;
+
+ if (result.shape)
+ setShape(result.shape.value());
+
+ setIntegrationFlag(result.integration_flag, result.integrator_type, true);
+
+ if (result.sum_background.isValid())
+ _sumBackground = result.sum_background;
+
+ if (result.sum_intensity.isValid())
+ _sumIntensity = result.sum_intensity;
+
+ if (result.profile_background.isValid())
+ _profileBackground = result.profile_background;
+
+ if (result.profile_intensity.isValid())
+ _profileIntensity = result.profile_intensity;
+
+ //_shape = integrator.fitShape(); // TODO: test, reactivate ???
}
ReciprocalVector Peak3D::q() const
@@ -276,7 +281,7 @@ ReciprocalVector Peak3D::q() const
double Peak3D::d() const
{
- double modq = q().rowVector().norm();
+ const double modq = q().rowVector().norm();
if (modq < _sigma2_eps)
return 0.0;
return 1.0 / modq;
@@ -311,15 +316,15 @@ Ellipsoid Peak3D::qShape() const
throw std::range_error("Interpolation error");
}
- Eigen::Vector3d q0 = q().rowVector();
+ const Eigen::Vector3d q0 = q().rowVector();
// Jacobian of map from detector coords to sample q space
- Eigen::Matrix3d J = state.jacobianQ(p[0], p[1]);
+ const Eigen::Matrix3d J = state.jacobianQ(p[0], p[1]);
const Eigen::Matrix3d JI = J.inverse();
// inverse covariance matrix in sample q space
const Eigen::Matrix3d q_inv_cov = JI.transpose() * _shape.metric() * JI;
- return Ellipsoid(q0, q_inv_cov);
+ return {q0, q_inv_cov};
}
bool Peak3D::caughtByFilter() const
@@ -389,7 +394,7 @@ const MillerIndex& Peak3D::hkl() const
void Peak3D::setMillerIndices()
{
if (unitCell()) {
- ReciprocalVector rv = q();
+ const ReciprocalVector rv = q();
if (rv.isValid()) {
_hkl = MillerIndex(q(), *unitCell());
} else {
@@ -429,6 +434,13 @@ void Peak3D::setIntegrationFlag(
}
}
+void Peak3D::setMasked()
+{
+ _rejection_flag = RejectionFlag::Masked;
+ _sum_integration_flag = RejectionFlag::Masked;
+ _profile_integration_flag = RejectionFlag::Masked;
+}
+
RejectionFlag Peak3D::rejectionFlag() const
{
if (_profile_integration_flag != RejectionFlag::NotRejected)
diff --git a/core/peak/Peak3D.h b/core/peak/Peak3D.h
index bbf136785..3eec28dc7 100644
--- a/core/peak/Peak3D.h
+++ b/core/peak/Peak3D.h
@@ -29,6 +29,7 @@ enum class RegionType;
enum class IntegratorType;
class ReciprocalVector;
class UnitCell;
+struct ComputeResult;
/*! \addtogroup python_api
* @{*/
@@ -61,6 +62,7 @@ enum class RejectionFlag {
NoShapeModel,
NoISigmaMinimum,
TooWide,
+ BadGaussianFit,
PredictionUpdateFailure, // from refiner
ManuallyRejected,
OutsideIndexingTol,
@@ -175,9 +177,7 @@ class Peak3D {
//! Update the integration parameters for this peak
void updateIntegration(
- const std::vector<Intensity>& rockingCurve, const Intensity& sumBkg,
- const Intensity& profBkg, const Intensity& meanBkgGradient, const Intensity& sumInt,
- const Intensity& profInt, double peakEnd, double bkgBegin, double bkgEnd,
+ const ComputeResult& result, double peakEnd, double bkgBegin, double bkgEnd,
const RegionType& regionType);
//! Return the q vector of the peak, transformed into sample coordinates.
ReciprocalVector q() const;
@@ -203,6 +203,8 @@ class Peak3D {
//! Set the reason for rejection during integration
void setIntegrationFlag(
const RejectionFlag& flag, const IntegratorType& integrator, bool overwrite = false);
+ //! Label the peak as masked
+ void setMasked();
//! Return the rejection flag only
RejectionFlag getRejectionFlag() const { return _rejection_flag; };
//! Return the sum integration flag only
diff --git a/core/peak/PeakData.cpp b/core/peak/PeakData.cpp
index d52a7bc80..be9d11d9f 100644
--- a/core/peak/PeakData.cpp
+++ b/core/peak/PeakData.cpp
@@ -46,11 +46,6 @@ const std::deque<double>& PeakData::gradients() const
void PeakData::standardizeCoords()
{
- if (_peak == nullptr) {
- throw std::runtime_error(
- "PeakData::standardizeCoords() cannot be called if _peak is nullptr");
- }
-
_coords.resize(_events.size());
for (size_t i = 0; i < _events.size(); ++i)
@@ -64,6 +59,13 @@ void PeakData::addEvent(const DetectorEvent& ev, double count, double gradient /
_gradients.push_back(gradient);
}
+void PeakData::append(const PeakData& other)
+{
+ _events.insert(_events.end(), other._events.begin(), other._events.end());
+ _counts.insert(_counts.end(), other._counts.begin(), other._counts.end());
+ _gradients.insert(_gradients.end(), other._gradients.begin(), other._gradients.end());
+}
+
void PeakData::reset()
{
std::deque<DetectorEvent> e;
@@ -77,4 +79,9 @@ void PeakData::reset()
std::swap(_coords, crds);
}
+bool PeakData::empty() const
+{
+ return _events.empty();
+}
+
} // namespace ohkl
diff --git a/core/peak/PeakData.h b/core/peak/PeakData.h
index 1817a3379..335f18985 100644
--- a/core/peak/PeakData.h
+++ b/core/peak/PeakData.h
@@ -37,8 +37,14 @@ class PeakData {
void standardizeCoords();
//! Add an event to the list of events.
void addEvent(const DetectorEvent& ev, double count, double gradient = 0);
+ //! Append a PeakData object to this one
+ void append(const PeakData& other);
//! Clear the events
void reset();
+ //! Is this PeakData empty?:
+ bool empty() const;
+ //! Get the peak pointer
+ Peak3D* peak() const { return _peak; };
private:
Peak3D* _peak;
--
GitLab
From 04b3ebae1e130755771bd931d1b94520ada23a88 Mon Sep 17 00:00:00 2001
From: Zamaan Raza <z.raza@fz-juelich.de>
Date: Wed, 22 May 2024 14:56:16 +0200
Subject: [PATCH 08/22] Parallelize integration
---
core/experiment/Integrator.cpp | 18 +-
core/experiment/Integrator.h | 7 +-
core/integration/IIntegrator.cpp | 219 ++++++++++++++++++
core/integration/IIntegrator.h | 4 +
core/peak/Peak3D.h | 2 +-
.../functional/TestProfile3DIntegration.cpp | 6 +-
test/cpp/functional/TestPxsumIntegration.cpp | 2 +
test/cpp/functional/TestShapeModel.cpp | 1 +
8 files changed, 247 insertions(+), 12 deletions(-)
diff --git a/core/experiment/Integrator.cpp b/core/experiment/Integrator.cpp
index 74bc47bb0..41307e94b 100644
--- a/core/experiment/Integrator.cpp
+++ b/core/experiment/Integrator.cpp
@@ -50,7 +50,7 @@ void Integrator::integratePeaks(
peaks->resetIntegrationFlags(integrator_type);
IIntegrator* integrator = getIntegrator(integrator_type);
integrator->setParameters(*_params);
- integrator->integrate(peaks->getPeakList(), peaks->shapeModel(), data);
+ integrator->parallelIntegrate(peaks->getPeakList(), peaks->shapeModel(), data);
peaks->setIntegrated(true);
if (_params->use_gradient)
peaks->setBkgGradient(true);
@@ -67,15 +67,17 @@ void Integrator::integratePeaks(
}
void Integrator::integratePeaks(
- sptrDataSet data, PeakCollection* peaks, IntegrationParameters* params, ShapeModel* shapes)
+ sptrDataSet data, PeakCollection* peaks, IntegrationParameters* params, ShapeModel* shapes,
+ bool parallel)
{
ohklLog(
Level::Info,
"Integrator::integratePeaks: integrating PeakCollection '" + peaks->name() + "'");
peaks->resetIntegrationFlags(_params->integrator_type);
IIntegrator* integrator = getIntegrator(_params->integrator_type);
+ integrator->setParallel(parallel);
integrator->setParameters(*params);
- integrator->integrate(peaks->getPeakList(), shapes, data);
+ integrator->parallelIntegrate(peaks->getPeakList(), shapes, data);
peaks->setIntegrated(true);
if (params->use_gradient)
peaks->setBkgGradient(true);
@@ -91,13 +93,14 @@ void Integrator::integratePeaks(
delete integrator;
}
-void Integrator::integrateFoundPeaks(PeakFinder* peak_finder)
+ void Integrator::integrateFoundPeaks(PeakFinder* peak_finder, bool parallel)
{
ohklLog(Level::Info, "Integrator::integrateFoundPeaks");
IIntegrator* integrator = getIntegrator(IntegratorType::PixelSum);
+ integrator->setParallel(parallel);
integrator->setParameters(*_params);
- integrator->integrate(peak_finder->currentPeaks(), nullptr, peak_finder->currentData());
+ integrator->parallelIntegrate(peak_finder->currentPeaks(), nullptr, peak_finder->currentData());
_n_peaks = 0;
_n_valid = 0;
@@ -115,14 +118,15 @@ void Integrator::integrateFoundPeaks(PeakFinder* peak_finder)
void Integrator::integrateShapeModel(
std::vector<Peak3D*> fit_peaks, sptrDataSet data, ShapeModel* shape_model, const AABB& aabb,
- const ShapeModelParameters& params)
+ const ShapeModelParameters& params, bool parallel)
{
ohklLog(Level::Info, "Integrator::integrateShapeModel");
ShapeIntegrator integrator{shape_model, aabb, params.nbins_x, params.nbins_y, params.nbins_z};
+ integrator.setParallel(parallel);
if (_handler)
integrator.setHandler(_handler);
integrator.setParameters(params);
- integrator.integrate(fit_peaks, shape_model, data);
+ integrator.parallelIntegrate(fit_peaks, shape_model, data);
}
void Integrator::setParameters(std::shared_ptr<IntegrationParameters> params)
diff --git a/core/experiment/Integrator.h b/core/experiment/Integrator.h
index 270450b85..897d9c560 100644
--- a/core/experiment/Integrator.h
+++ b/core/experiment/Integrator.h
@@ -54,13 +54,14 @@ class Integrator {
void integratePeaks(IntegratorType integrator_type, sptrDataSet data, PeakCollection* peaks);
//! Integrate a peak collection
void integratePeaks(
- sptrDataSet data, PeakCollection* peaks, IntegrationParameters* params, ShapeModel* shapes);
+ sptrDataSet data, PeakCollection* peaks, IntegrationParameters* params, ShapeModel* shapes,
+ bool parallel = true);
//! Integrate peaks found by _peak_finder
- void integrateFoundPeaks(PeakFinder* peak_finder);
+ void integrateFoundPeaks(PeakFinder* peak_finder, bool parallel = true);
//! Integrate the shape collection
void integrateShapeModel(
std::vector<Peak3D*> peaks, sptrDataSet data, ShapeModel* shape_model, const AABB& aabb,
- const ShapeModelParameters& params);
+ const ShapeModelParameters& params, bool parallel = true);
//! Set the parameters
void setParameters(std::shared_ptr<IntegrationParameters> params);
//! Get the parameters
diff --git a/core/integration/IIntegrator.cpp b/core/integration/IIntegrator.cpp
index 9ba640452..5061ae568 100644
--- a/core/integration/IIntegrator.cpp
+++ b/core/integration/IIntegrator.cpp
@@ -17,6 +17,7 @@
#include "base/geometry/Ellipsoid.h"
#include "base/mask/IMask.h"
#include "base/utils/Logger.h"
+#include "base/utils/ParallelFor.h"
#include "base/utils/ProgressHandler.h"
#include "core/data/DataSet.h"
#include "core/peak/Peak3D.h"
@@ -25,6 +26,7 @@
#include "tables/crystal/UnitCell.h"
#include <algorithm>
+#include <mutex>
namespace ohkl {
@@ -249,6 +251,223 @@ void IIntegrator::integrate(
ohklLog(Level::Info, "IIntegrator::integrate: end; ", _n_failures, " failures");
}
+void IIntegrator::parallelIntegrate(
+ std::vector<Peak3D*> peaks, ShapeModel* shape_model, sptrDataSet data)
+{
+ ohklLog(Level::Info, "IIntegrator::parallelIntegrate: integrating ", peaks.size(), " peaks");
+ _params.log(Level::Info);
+ if (shape_model)
+ shape_model->parameters()->log(Level::Info);
+
+ // integrate only those peaks that belong to the specified dataset
+ auto it = std::remove_if(peaks.begin(), peaks.end(), [&](const Peak3D* peak) {
+ return peak->dataSet() != data;
+ });
+ peaks.erase(it, peaks.end());
+ std::ostringstream oss;
+ oss << "Integrating " << peaks.size() << " peaks";
+ ohklLog(Level::Info, "IIntegrator::parallelIntegrate: integrating ", peaks.size(), " peaks");
+ if (_handler) {
+ _handler->setStatus(oss.str().c_str());
+ _handler->setProgress(0);
+ }
+
+ _profile_integration = false;
+ if (_params.integrator_type == IntegratorType::Profile1D
+ || _params.integrator_type == IntegratorType::Profile3D
+ || _params.integrator_type == IntegratorType::Gaussian
+ || _params.integrator_type == IntegratorType::ISigma)
+ _profile_integration = true;
+
+
+ size_t idx = 0;
+ _n_frames_done = 0;
+
+ std::map<Peak3D*, std::unique_ptr<IntegrationRegion>> regions;
+ std::map<Peak3D*, bool> integrated;
+
+ if (_params.region_type == ohkl::RegionType::VariableEllipsoid) {
+ _peak_end = _params.peak_end;
+ _bkg_begin = _params.bkg_begin;
+ _bkg_end = _params.bkg_end;
+ } else {
+ _peak_end = _params.fixed_peak_end;
+ _bkg_begin = _params.fixed_bkg_begin;
+ _bkg_end = _params.fixed_bkg_end;
+ }
+
+ for (auto peak : peaks) {
+ if (peak->isRejectedFor(RejectionFlag::Extinct)) {
+ if (!peak->enabled())
+ continue;
+ }
+
+ regions.emplace(std::make_pair(
+ peak,
+ std::make_unique<IntegrationRegion>(
+ peak, _peak_end, _bkg_begin, _bkg_end, _params.region_type)));
+ integrated.emplace(std::make_pair(peak, false));
+
+ // ignore partials
+ auto bb = regions.at(peak)->peakBB();
+ auto data = peak->dataSet();
+ auto lo = bb.lower();
+ auto hi = bb.upper();
+
+ double width = hi[2] - lo[2];
+ if (_params.use_max_width && width > _params.max_width)
+ peak->setIntegrationFlag(RejectionFlag::TooWide, _params.integrator_type);
+
+ if (lo[0] < 0 || lo[1] < 0 || lo[2] < 0 || hi[0] >= data->nCols() || hi[1] >= data->nRows()
+ || hi[2] >= data->nFrames())
+ peak->setIntegrationFlag(RejectionFlag::InvalidRegion, _params.integrator_type);
+ }
+
+ // only integrate the peaks with valid integration regions
+ ohklLog(Level::Debug, "IIntegrator::integrate: remove invalid regions");
+ it = std::remove_if(peaks.begin(), peaks.end(), [&](Peak3D*& p) {
+ return regions.find(p) == regions.end();
+ });
+ peaks.erase(it, peaks.end());
+
+ // check for overlaps if requested
+ if (_params.remove_overlaps)
+ removeOverlaps(regions);
+
+ oss << "Building " << peaks.size() << " integration regions";
+ if (_handler) {
+ _handler->setStatus(oss.str().c_str());
+ _handler->setProgress(0);
+ }
+ ohklLog(Level::Debug, "IIntegrator::parallelIntegrate: integration region loop");
+ std::mutex mut1;
+
+ parallel_for(data->nFrames(), [&](int start, int end) {
+ for (int idx = start; idx < end; ++idx) {
+
+ Eigen::MatrixXd current_frame, gradient;
+ Eigen::MatrixXi mask;
+ current_frame = data->transformedFrame(idx);
+ if (_params.use_gradient)
+ gradient =
+ data->gradientFrame(idx, _params.gradient_type, !_params.fft_gradient);
+
+ mask.resize(data->nRows(), data->nCols());
+ mask.setConstant(int(IntegrationRegion::EventType::EXCLUDED));
+
+ for (auto peak : peaks)
+ regions.at(peak)->updateMask(mask, idx);
+
+ std::vector<std::unique_ptr<PeakData>> peak_data_vector;
+ for (auto peak : peaks) {
+ auto* current_region = regions.at(peak).get();
+ PeakData peak_data(peak);
+ if (_params.use_gradient)
+ peak_data =
+ current_region->threadSafeAdvanceFrame(current_frame, mask, idx, &gradient);
+ else
+ peak_data =
+ current_region->threadSafeAdvanceFrame(current_frame, mask, idx);
+ if (!peak_data.empty())
+ peak_data_vector.push_back(std::make_unique<PeakData>(peak_data));
+ }
+
+ {
+ const std::lock_guard<std::mutex> lock(mut1);
+ for (auto& peak_data : peak_data_vector)
+ regions[peak_data->peak()]->appendPeakData(*peak_data);
+ }
+
+ if (_handler) {
+ const double progress = _n_frames_done++ * 100.0 / data->nFrames();
+ _handler->setProgress(progress);
+ }
+ }
+ }, _thread_parallel, _max_threads);
+
+ if (_handler)
+ _handler->setProgress(100);
+
+ _n_peaks_done = 0;
+ if (_handler)
+ _handler->setProgress(0);
+
+ ohklLog(Level::Debug, "IIntegrator::parallelIntegrate: compute loop");
+ _n_failures = 0;
+ std::mutex mut2;
+ parallel_for(peaks.size(), [&](int start, int end) {
+ for (int idx = start; idx < end; ++idx) {
+ auto* peak = peaks.at(idx);
+ auto* current_region = regions.at(peak).get();
+ const auto& counts = current_region->peakData().counts();
+
+ RejectionFlag reject(RejectionFlag::NotRejected);
+
+ if (_params.discard_saturated) {
+ double max = 0;
+ if (!counts.empty()) // std::max on empty vector segfaults under MacOS
+ max = *std::max_element(counts.begin(), counts.end());
+ if (max > _params.max_counts)
+ reject = RejectionFlag::SaturatedPixel;
+ }
+
+ if (_params.use_max_width) {
+ auto aabb = current_region->peakBB();
+ const auto lower = aabb.lower();
+ const auto upper = aabb.upper();
+ const double width = upper[2] - lower[2];
+ if (width > static_cast<double>(_params.max_width))
+ reject = RejectionFlag::TooWide;
+ }
+
+ bool reintegrate = true;
+ if (_profile_integration) {
+ if (_params.use_max_strength
+ && peak->sumIntensity().strength() > _params.max_strength)
+ reintegrate = false;
+ if (_params.use_max_d && peak->d() > _params.max_d)
+ reintegrate = false;
+ }
+
+ ComputeResult compute_result;
+ if (!reintegrate) {
+ compute_result.integrator_type = _params.integrator_type;
+ compute_result.rocking_curve = peak->rockingCurve();
+ compute_result.profile_intensity = peak->sumIntensity();
+ compute_result.profile_background = peak->sumBackground();
+ compute_result.integration_flag = peak->sumRejectionFlag();
+ } else {
+ if (reject == RejectionFlag::NotRejected) {
+ current_region->peakData().standardizeCoords();
+ compute_result = compute(peak, shape_model, *current_region);
+ } else
+ compute_result.integration_flag = reject;
+
+ if (compute_result.integration_flag != RejectionFlag::NotRejected)
+ ++_n_failures;
+ }
+
+ {
+ const std::lock_guard<std::mutex> lock(mut2);
+ peak->updateIntegration(
+ compute_result, _params.peak_end, _params.bkg_begin, _params.bkg_end,
+ _params.region_type);
+ current_region->reset();
+ }
+
+ if (_handler) {
+ const double progress = _n_peaks_done++ * 100.0 / peaks.size();
+ _handler->setProgress(progress);
+ }
+ }
+ }, _thread_parallel, _max_threads);
+
+ if (_handler)
+ _handler->setProgress(100);
+
+ ohklLog(Level::Info, "IIntegrator::parallelIntegrate: end; ", _n_failures, " failures");
+}
+
void IIntegrator::setHandler(sptrProgressHandler handler)
{
_handler = handler;
diff --git a/core/integration/IIntegrator.h b/core/integration/IIntegrator.h
index c5ad87964..fdc988718 100644
--- a/core/integration/IIntegrator.h
+++ b/core/integration/IIntegrator.h
@@ -125,10 +125,14 @@ class IIntegrator {
virtual ~IIntegrator() = default;
//! Integrate a vector of peaks
void integrate(std::vector<Peak3D*> peaks, ShapeModel* shape_model, sptrDataSet data);
+ //! Thread-parallel integrate
+ void parallelIntegrate(std::vector<Peak3D*> peaks, ShapeModel* shape_model, sptrDataSet data);
//! Set the progress handler.
void setHandler(sptrProgressHandler handler);
//! Assign a parameter set to the integrator
void setParameters(const IntegrationParameters& params);
+ //! Toggle parallel integration
+ void setParallel(bool parallel) { _thread_parallel = parallel; };
//! Remove overlapping peak intensity regions
void removeOverlaps(const std::map<Peak3D*, std::unique_ptr<IntegrationRegion>>& regions);
diff --git a/core/peak/Peak3D.h b/core/peak/Peak3D.h
index 3eec28dc7..c3675ec74 100644
--- a/core/peak/Peak3D.h
+++ b/core/peak/Peak3D.h
@@ -36,7 +36,7 @@ struct ComputeResult;
//! Specifies reason why a peak was rejected
enum class RejectionFlag {
- NotRejected,
+ NotRejected = 0,
Masked,
OutsideThreshold,
OutsideFrames,
diff --git a/test/cpp/functional/TestProfile3DIntegration.cpp b/test/cpp/functional/TestProfile3DIntegration.cpp
index c0d9c3321..ca8bdfdca 100644
--- a/test/cpp/functional/TestProfile3DIntegration.cpp
+++ b/test/cpp/functional/TestProfile3DIntegration.cpp
@@ -57,7 +57,11 @@ TEST_CASE("test/data/TestProfile3DIntegration.cpp", "")
experiment.loadFromFile(filename);
const ohkl::sptrUnitCell cell = experiment.getSptrUnitCell("indexed");
- const ohkl::sptrDataSet data = experiment.getData("Scan I"); ohkl::PeakCollection* peaks = experiment.getPeakCollection("predicted");
+ const ohkl::sptrDataSet data = experiment.getData("Scan I");
+ ohkl::PeakCollection* peaks = experiment.getPeakCollection("predicted");
+
+ data->initBuffer(true);
+
auto* integrator = experiment.integrator();
auto* params = integrator->parameters();
diff --git a/test/cpp/functional/TestPxsumIntegration.cpp b/test/cpp/functional/TestPxsumIntegration.cpp
index e3dbe993f..72e3d2e60 100644
--- a/test/cpp/functional/TestPxsumIntegration.cpp
+++ b/test/cpp/functional/TestPxsumIntegration.cpp
@@ -57,6 +57,8 @@ TEST_CASE("test/data/TestPxsumIntegration.cpp", "")
const ohkl::sptrDataSet data = experiment.getData("Scan I");
ohkl::PeakCollection* peaks = experiment.getPeakCollection("predicted");
+ data->initBuffer(true);
+
auto* integrator = experiment.integrator();
auto* params = integrator->parameters();
diff --git a/test/cpp/functional/TestShapeModel.cpp b/test/cpp/functional/TestShapeModel.cpp
index 3c79be9a8..909b017d6 100644
--- a/test/cpp/functional/TestShapeModel.cpp
+++ b/test/cpp/functional/TestShapeModel.cpp
@@ -40,6 +40,7 @@ TEST_CASE("test/data/TestShapeModel.cpp", "")
ohkl::sptrUnitCell cell = experiment.getSptrUnitCell("accepted");
ohkl::sptrDataSet data = experiment.getData("testdata");
+ data->initBuffer(true);
auto* found_peaks = experiment.getPeakCollection("found");
found_peaks->computeSigmas();
--
GitLab
From 05b5ea81cd2f42f3e11068ab0e0931f86e1e3694 Mon Sep 17 00:00:00 2001
From: Zamaan Raza <z.raza@fz-juelich.de>
Date: Wed, 22 May 2024 14:58:27 +0200
Subject: [PATCH 09/22] Add missing integration YAML nodes
---
core/experiment/ExperimentYAML.cpp | 10 ++++++++++
1 file changed, 10 insertions(+)
diff --git a/core/experiment/ExperimentYAML.cpp b/core/experiment/ExperimentYAML.cpp
index 7843935bc..4ede768af 100644
--- a/core/experiment/ExperimentYAML.cpp
+++ b/core/experiment/ExperimentYAML.cpp
@@ -135,6 +135,11 @@ void ExperimentYAML::grabIntegrationParameters(IntegrationParameters* params)
params->skip_masked = getNode<bool>(branch, "skip_masked");
params->remove_overlaps = getNode<bool>(branch, "remove_overlaps");
params->use_max_strength = getNode<bool>(branch, "use_max_strength");
+ params->max_strength = getNode<double>(branch, "max_strength");
+ params->use_max_d = getNode<bool>(branch, "use_max_d");
+ params->max_d = getNode<double>(branch, "max_d");
+ params->discard_saturated = getNode<bool>(branch, "discard_saturated");
+ params->max_counts = getNode<double>(branch, "max_counts");
}
void ExperimentYAML::setIntegrationParameters(IntegrationParameters* params)
@@ -162,6 +167,11 @@ void ExperimentYAML::setIntegrationParameters(IntegrationParameters* params)
int_node["skip_masked"] = params->skip_masked;
int_node["remove_overlaps"] = params->remove_overlaps;
int_node["use_max_strength"] = params->use_max_strength;
+ int_node["max_strength"] = params->max_strength;
+ int_node["use_max_d"] = params->use_max_d;
+ int_node["max_d"] = params->max_d;
+ int_node["discard_saturated"] = params->discard_saturated;
+ int_node["max_counts"] = params->max_counts;
}
void ExperimentYAML::grabPeakFinderParameters(PeakFinderParameters* params)
--
GitLab
From 1abd94d7c32b9b74ec61f1394faf464e79529329 Mon Sep 17 00:00:00 2001
From: Zamaan Raza <z.raza@fz-juelich.de>
Date: Wed, 22 May 2024 15:26:27 +0200
Subject: [PATCH 10/22] Fix masking loop in parallel
---
core/integration/IIntegrator.cpp | 15 ++++++++++++++-
1 file changed, 14 insertions(+), 1 deletion(-)
diff --git a/core/integration/IIntegrator.cpp b/core/integration/IIntegrator.cpp
index 5061ae568..0b57a6ac2 100644
--- a/core/integration/IIntegrator.cpp
+++ b/core/integration/IIntegrator.cpp
@@ -280,7 +280,6 @@ void IIntegrator::parallelIntegrate(
_profile_integration = true;
- size_t idx = 0;
_n_frames_done = 0;
std::map<Peak3D*, std::unique_ptr<IntegrationRegion>> regions;
@@ -297,6 +296,7 @@ void IIntegrator::parallelIntegrate(
}
for (auto peak : peaks) {
+ peak->setIntegrationFlag(RejectionFlag::NotRejected, _params.integrator_type, true);
if (peak->isRejectedFor(RejectionFlag::Extinct)) {
if (!peak->enabled())
continue;
@@ -392,6 +392,19 @@ void IIntegrator::parallelIntegrate(
if (_handler)
_handler->setProgress(0);
+ for (auto peak : peaks) {
+ if (_params.skip_masked) {
+ Ellipsoid shape = peak->shape();
+ shape.scale(_bkg_end);
+ for (const auto* mask : data->masks()) {
+ if (mask->collide(shape)) {
+ peak->setMasked();
+ continue;
+ }
+ }
+ }
+ }
+
ohklLog(Level::Debug, "IIntegrator::parallelIntegrate: compute loop");
_n_failures = 0;
std::mutex mut2;
--
GitLab
From db53f8720695b0710b3157a32c110bd910d51427 Mon Sep 17 00:00:00 2001
From: Zamaan Raza <z.raza@fz-juelich.de>
Date: Wed, 22 May 2024 15:47:11 +0200
Subject: [PATCH 11/22] Parallelise shape model integration
---
core/shape/ShapeModel.cpp | 20 +++++++++++++------
core/shape/ShapeModel.h | 7 +++++--
.../functional/TestProfile3DIntegration.cpp | 2 --
3 files changed, 19 insertions(+), 10 deletions(-)
diff --git a/core/shape/ShapeModel.cpp b/core/shape/ShapeModel.cpp
index 942deec56..9986cec3a 100644
--- a/core/shape/ShapeModel.cpp
+++ b/core/shape/ShapeModel.cpp
@@ -121,8 +121,14 @@ struct FitData {
}
};
-ShapeModel::ShapeModel()
- : _id(0), _profiles(), _choleskyD(), _choleskyM(), _choleskyS(), _handler(nullptr)
+ShapeModel::ShapeModel(bool thread_parallel)
+ : _id(0)
+ , _profiles()
+ , _choleskyD()
+ , _choleskyM()
+ , _choleskyS()
+ , _handler(nullptr)
+ , _thread_parallel(thread_parallel)
{
_choleskyD.fill(1e-6);
_choleskyM.fill(1e-6);
@@ -135,7 +141,7 @@ ShapeModel::ShapeModel(const std::string& name) : ShapeModel()
_name = name;
}
-ShapeModel::ShapeModel(std::shared_ptr<ShapeModelParameters> params)
+ShapeModel::ShapeModel(std::shared_ptr<ShapeModelParameters> params, bool thread_parallel)
: _id(0)
, _profiles()
, _choleskyD()
@@ -143,6 +149,7 @@ ShapeModel::ShapeModel(std::shared_ptr<ShapeModelParameters> params)
, _choleskyS()
, _params(params)
, _handler(nullptr)
+ , _thread_parallel(thread_parallel)
{
_choleskyD.fill(1e-6);
@@ -522,9 +529,9 @@ void ShapeModel::integrate(
int_params.fixed_bkg_begin = _params->fixed_bkg_begin;
int_params.fixed_bkg_end = _params->fixed_bkg_end;
integrator.setHandler(handler);
+ integrator.setParallel(_thread_parallel);
integrator.setParameters(int_params);
-
- integrator.integrate(peaks, this, data);
+ integrator.parallelIntegrate(peaks, this, data);
ohklLog(Level::Info, "ShapeModel::integrate: finished integrating shapes");
}
@@ -564,7 +571,8 @@ void ShapeModel::build(PeakCollection* peaks, sptrDataSet data)
int_params.fixed_bkg_end = _params->fixed_bkg_end;
int_params.region_type = _params->region_type;
integrator.setParameters(int_params);
- integrator.integrate(fit_peaks, this, data);
+ integrator.setParallel(_thread_parallel);
+ integrator.parallelIntegrate(fit_peaks, this, data);
ohklLog(Level::Info, "ShapeModel::build: finished integrating shapes");
// ohklLog(Level::Info, "ShapeModel::build: updating fit");
// updateFit(1000);
diff --git a/core/shape/ShapeModel.h b/core/shape/ShapeModel.h
index ca1b7d259..12c35fb47 100644
--- a/core/shape/ShapeModel.h
+++ b/core/shape/ShapeModel.h
@@ -77,9 +77,9 @@ class ShapeModel {
//! Construct an empty collection.
//! @param detector_coords if true, store profiles in detector coordinates;
//! otherwise store in Kabsch coordinates
- ShapeModel();
+ ShapeModel(bool thread_parallel = true);
ShapeModel(const std::string& name);
- ShapeModel(std::shared_ptr<ShapeModelParameters> params);
+ ShapeModel(std::shared_ptr<ShapeModelParameters> params, bool thread_parallel = true);
//! Get the integer id
unsigned int id() const { return _id; };
@@ -185,6 +185,9 @@ class ShapeModel {
//! Progress handler
sptrProgressHandler _handler;
+
+ //! Parallelise integration
+ bool _thread_parallel;
};
/*! @}*/
diff --git a/test/cpp/functional/TestProfile3DIntegration.cpp b/test/cpp/functional/TestProfile3DIntegration.cpp
index ca8bdfdca..fba52cbdf 100644
--- a/test/cpp/functional/TestProfile3DIntegration.cpp
+++ b/test/cpp/functional/TestProfile3DIntegration.cpp
@@ -60,8 +60,6 @@ TEST_CASE("test/data/TestProfile3DIntegration.cpp", "")
const ohkl::sptrDataSet data = experiment.getData("Scan I");
ohkl::PeakCollection* peaks = experiment.getPeakCollection("predicted");
- data->initBuffer(true);
-
auto* integrator = experiment.integrator();
auto* params = integrator->parameters();
--
GitLab
From 0064f3aea4b0aaf792d65b809bf25d6bad71af47 Mon Sep 17 00:00:00 2001
From: Zamaan Raza <z.raza@fz-juelich.de>
Date: Thu, 23 May 2024 08:39:12 +0200
Subject: [PATCH 12/22] Parallelise prediction
- Fix GUI enum class items
---
core/loader/XFileHandler.cpp | 4 +--
core/peak/Qs2Events.cpp | 32 +++++++++++++-------
core/peak/Qs2Events.h | 2 +-
gui/subframe_integrate/SubframeIntegrate.cpp | 3 ++
gui/subframe_integrate/SubframeIntegrate.h | 3 +-
5 files changed, 29 insertions(+), 15 deletions(-)
diff --git a/core/loader/XFileHandler.cpp b/core/loader/XFileHandler.cpp
index 369f136c1..671522665 100644
--- a/core/loader/XFileHandler.cpp
+++ b/core/loader/XFileHandler.cpp
@@ -14,6 +14,7 @@
#include "core/loader/XFileHandler.h"
+#include <string>
#include <fstream>
#include <iterator>
#include <sstream>
@@ -88,9 +89,8 @@ void XFileHandler::readXFile(int frame)
if (stop) {
_mask = Eigen::MatrixXi(refl.maskx, refl.masky);
break;
- } else {
- _reflections.emplace_back(refl);
}
+ _reflections.emplace_back(refl);
}
// The integration mask
diff --git a/core/peak/Qs2Events.cpp b/core/peak/Qs2Events.cpp
index eb1064611..bda9b2666 100644
--- a/core/peak/Qs2Events.cpp
+++ b/core/peak/Qs2Events.cpp
@@ -16,12 +16,15 @@
#include "base/geometry/DirectVector.h"
#include "base/utils/Logger.h"
+#include "base/utils/ParallelFor.h"
#include "base/utils/ProgressHandler.h"
#include "core/data/DataSet.h" // TODO mv interpolatedState
#include "core/detector/Detector.h"
#include "core/instrument/InterpolatedState.h"
#include "tables/crystal/MillerIndex.h"
+#include <mutex>
+
namespace ohkl {
// Tolerance for bisection search
@@ -37,7 +40,7 @@ auto compute_sign = [](const Eigen::RowVector3d& q, const InterpolatedState& sta
std::vector<std::pair<MillerIndex, DetectorEvent>> algo::qMap2Events(
const std::vector<std::pair<MillerIndex, ReciprocalVector>>& sample_qs,
const InstrumentStateList& states, const Detector& detector, const int n_intervals,
- sptrProgressHandler handler /* = nullptr */)
+ sptrProgressHandler handler /* = nullptr */, bool thread_parallel /* = true */)
{
ohklLog(
Level::Debug, "algo::Qs2Events::qVectorList2Events: processing ", sample_qs.size(),
@@ -45,7 +48,6 @@ std::vector<std::pair<MillerIndex, DetectorEvent>> algo::qMap2Events(
std::vector<std::pair<MillerIndex, DetectorEvent>> events;
- int count = 1;
if (handler) {
std::ostringstream oss;
oss << "Transforming " << sample_qs.size() << " q-vectors to detector events";
@@ -53,15 +55,23 @@ std::vector<std::pair<MillerIndex, DetectorEvent>> algo::qMap2Events(
handler->setProgress(0);
}
-// for each sample q, determine the rotation that makes it intersect the Ewald sphere
- for (const auto& [hkl, sample_q] : sample_qs) {
- std::vector<DetectorEvent> new_events =
- qVector2Events(sample_q, states, detector, n_intervals);
- for (auto event : new_events)
- events.push_back({hkl, event});
- if (handler)
- handler->setProgress(++count * 100.0 / sample_qs.size());
- }
+ // for each sample q, determine the rotation that makes it intersect the Ewald sphere
+ std::mutex mutex;
+ std::atomic_int count = 1;
+ parallel_for(sample_qs.size(), [&](int start, int end) {
+ for (int idx = start; idx < end; ++idx) {
+ std::vector<DetectorEvent> new_events =
+ qVector2Events(sample_qs.at(idx).second, states, detector, n_intervals);
+ {
+ std::lock_guard<std::mutex> lock(mutex);
+ for (auto event : new_events)
+ events.push_back({sample_qs.at(idx).first, event});
+ }
+ if (handler)
+ handler->setProgress(++count * 100.0 / sample_qs.size());
+ }
+ }, thread_parallel);
+
if (handler)
handler->setProgress(100);
ohklLog(
diff --git a/core/peak/Qs2Events.h b/core/peak/Qs2Events.h
index aa41d7921..164e8179a 100644
--- a/core/peak/Qs2Events.h
+++ b/core/peak/Qs2Events.h
@@ -48,7 +48,7 @@ namespace algo {
std::vector<std::pair<MillerIndex, DetectorEvent>> qMap2Events(
const std::vector<std::pair<MillerIndex, ReciprocalVector>>& sample_qs,
const InstrumentStateList& states, const Detector& detector, const int n_intervals,
- sptrProgressHandler handler = nullptr);
+ sptrProgressHandler handler = nullptr, bool thread_parallel = true);
//! Returns detector events corresponding to the list of q values.
std::vector<DetectorEvent> qVectorList2Events(
const std::vector<ReciprocalVector>& sample_qs, const InstrumentStateList& states,
diff --git a/gui/subframe_integrate/SubframeIntegrate.cpp b/gui/subframe_integrate/SubframeIntegrate.cpp
index 981383edd..2220098dc 100644
--- a/gui/subframe_integrate/SubframeIntegrate.cpp
+++ b/gui/subframe_integrate/SubframeIntegrate.cpp
@@ -17,6 +17,7 @@
#include "core/experiment/Experiment.h"
#include "core/experiment/Integrator.h"
+#include "core/integration/IIntegrator.h"
#include "core/shape/PeakCollection.h"
#include "core/shape/ShapeModel.h"
#include "gui/MainWin.h" // gGui
@@ -441,6 +442,8 @@ void SubframeIntegrate::setIntegrateUp()
for (std::size_t idx = 0; idx < static_cast<int>(ohkl::IntegratorType::Count); ++idx) {
const std::string integrator_type =
_integrator_strings.at(static_cast<ohkl::IntegratorType>(idx));
+ if (static_cast<ohkl::IntegratorType>(idx) == ohkl::IntegratorType::Shape)
+ continue;
_integrator_combo->addItem(QString::fromStdString(integrator_type));
}
diff --git a/gui/subframe_integrate/SubframeIntegrate.h b/gui/subframe_integrate/SubframeIntegrate.h
index a019fdc18..954b9d893 100644
--- a/gui/subframe_integrate/SubframeIntegrate.h
+++ b/gui/subframe_integrate/SubframeIntegrate.h
@@ -165,7 +165,8 @@ class SubframeIntegrate : public QWidget {
{ohkl::IntegratorType::Gaussian, "Gaussian integrator"},
{ohkl::IntegratorType::ISigma, "I/Sigma integrator"},
{ohkl::IntegratorType::Profile1D, "1D Profile integrator"},
- {ohkl::IntegratorType::Profile3D, "3D Profile integrator"}};
+ {ohkl::IntegratorType::Profile3D, "3D Profile integrator"},
+ {ohkl::IntegratorType::Shape, "Shape model integrator"}};
const std::map<ohkl::GradientKernel, QString> _kernel_description{
{ohkl::GradientKernel::CentralDifference, "Central difference"},
--
GitLab
From a2d04481a5ff5026fb7a2ceb5dae3544c0c6122c Mon Sep 17 00:00:00 2001
From: Zamaan Raza <z.raza@fz-juelich.de>
Date: Mon, 27 May 2024 11:00:45 +0200
Subject: [PATCH 13/22] Fix incorrect masked peak count
---
core/integration/IIntegrator.cpp | 27 +++++++++++++--------------
1 file changed, 13 insertions(+), 14 deletions(-)
diff --git a/core/integration/IIntegrator.cpp b/core/integration/IIntegrator.cpp
index 0b57a6ac2..8c3fb7bd7 100644
--- a/core/integration/IIntegrator.cpp
+++ b/core/integration/IIntegrator.cpp
@@ -296,7 +296,6 @@ void IIntegrator::parallelIntegrate(
}
for (auto peak : peaks) {
- peak->setIntegrationFlag(RejectionFlag::NotRejected, _params.integrator_type, true);
if (peak->isRejectedFor(RejectionFlag::Extinct)) {
if (!peak->enabled())
continue;
@@ -392,19 +391,6 @@ void IIntegrator::parallelIntegrate(
if (_handler)
_handler->setProgress(0);
- for (auto peak : peaks) {
- if (_params.skip_masked) {
- Ellipsoid shape = peak->shape();
- shape.scale(_bkg_end);
- for (const auto* mask : data->masks()) {
- if (mask->collide(shape)) {
- peak->setMasked();
- continue;
- }
- }
- }
- }
-
ohklLog(Level::Debug, "IIntegrator::parallelIntegrate: compute loop");
_n_failures = 0;
std::mutex mut2;
@@ -433,6 +419,17 @@ void IIntegrator::parallelIntegrate(
reject = RejectionFlag::TooWide;
}
+ if (_params.skip_masked) {
+ Ellipsoid shape = peak->shape();
+ shape.scale(_bkg_end);
+ for (const auto* mask : data->masks()) {
+ if (mask->collide(shape)) {
+ reject = RejectionFlag::Masked;
+ break;
+ }
+ }
+ }
+
bool reintegrate = true;
if (_profile_integration) {
if (_params.use_max_strength
@@ -465,6 +462,8 @@ void IIntegrator::parallelIntegrate(
peak->updateIntegration(
compute_result, _params.peak_end, _params.bkg_begin, _params.bkg_end,
_params.region_type);
+ if (reject == RejectionFlag::Masked)
+ peak->setMasked();
current_region->reset();
}
--
GitLab
From b8403d411a126c6a73bc1d26d6fb53f550896509 Mon Sep 17 00:00:00 2001
From: Zamaan Raza <z.raza@fz-juelich.de>
Date: Mon, 27 May 2024 15:40:37 +0200
Subject: [PATCH 14/22] Parallelize shape assignment
---
core/integration/IIntegrator.cpp | 14 +++++++-------
core/shape/ShapeModel.cpp | 26 +++++++++++++++-----------
2 files changed, 22 insertions(+), 18 deletions(-)
diff --git a/core/integration/IIntegrator.cpp b/core/integration/IIntegrator.cpp
index 8c3fb7bd7..3396fd737 100644
--- a/core/integration/IIntegrator.cpp
+++ b/core/integration/IIntegrator.cpp
@@ -265,12 +265,7 @@ void IIntegrator::parallelIntegrate(
});
peaks.erase(it, peaks.end());
std::ostringstream oss;
- oss << "Integrating " << peaks.size() << " peaks";
ohklLog(Level::Info, "IIntegrator::parallelIntegrate: integrating ", peaks.size(), " peaks");
- if (_handler) {
- _handler->setStatus(oss.str().c_str());
- _handler->setProgress(0);
- }
_profile_integration = false;
if (_params.integrator_type == IntegratorType::Profile1D
@@ -387,9 +382,14 @@ void IIntegrator::parallelIntegrate(
if (_handler)
_handler->setProgress(100);
- _n_peaks_done = 0;
- if (_handler)
+ oss.str(std::string()); // clear the stream
+ oss << "Integrating " << peaks.size() << " integration regions";
+ if (_handler) {
+ _handler->setStatus(oss.str().c_str());
_handler->setProgress(0);
+ }
+
+ _n_peaks_done = 0;
ohklLog(Level::Debug, "IIntegrator::parallelIntegrate: compute loop");
_n_failures = 0;
diff --git a/core/shape/ShapeModel.cpp b/core/shape/ShapeModel.cpp
index 9986cec3a..86726429f 100644
--- a/core/shape/ShapeModel.cpp
+++ b/core/shape/ShapeModel.cpp
@@ -19,6 +19,7 @@
#include "base/geometry/Ellipsoid.h"
#include "base/geometry/ReciprocalVector.h"
#include "base/utils/Logger.h"
+#include "base/utils/ParallelFor.h"
#include "base/utils/ProgressHandler.h"
#include "core/data/DataSet.h"
#include "core/data/DataTypes.h"
@@ -453,18 +454,21 @@ void ShapeModel::setPredictedShapes(PeakCollection* peaks)
_handler->setProgress(0);
}
- for (auto peak : peaks->getPeakList()) {
- // Skip the peak if any error occur when computing its mean covariance (e.g.
- // too few or no neighbouring peaks found)
- auto cov = meanCovariance(peak);
- Eigen::Vector3d center = peak->shape().center();
- peak->setShape(Ellipsoid(center, cov.inverse()));
-
- if (_handler) {
- double progress = ++count * 100.0 / npeaks;
- _handler->setProgress(progress);
+ auto peaklist = peaks->getPeakList();
+ parallel_for(peaklist.size(), [&](int start, int end) {
+ for (int idx = start; idx < end; ++idx) {
+ Peak3D* peak = peaklist.at(idx);
+ const auto cov = meanCovariance(peak);
+ const Eigen::Vector3d center = peak->shape().center();
+ peak->setShape(Ellipsoid(center, cov.inverse()));
+
+ if (_handler) {
+ const double progress = ++count * 100.0 / npeaks;
+ _handler->setProgress(progress);
+ }
}
- }
+ }, _thread_parallel);
+
ohklLog(Level::Info, "ShapeModel: finished computing shapes");
}
--
GitLab
From 5f9421ff51961d9006da09370743248ef2fc0c98 Mon Sep 17 00:00:00 2001
From: Zamaan Raza <z.raza@fz-juelich.de>
Date: Mon, 3 Jun 2024 10:22:01 +0200
Subject: [PATCH 15/22] Parallelise peak processing for GUI
---
gui/graphics/PeakCollectionGraphics.cpp | 49 +++++++++++++++----------
gui/graphics/PeakCollectionGraphics.h | 3 +-
gui/items/PeakCollectionItem.cpp | 27 +++++++++++---
gui/items/PeakCollectionItem.h | 2 +-
4 files changed, 54 insertions(+), 27 deletions(-)
diff --git a/gui/graphics/PeakCollectionGraphics.cpp b/gui/graphics/PeakCollectionGraphics.cpp
index 85430bac6..0f0d5cfd7 100644
--- a/gui/graphics/PeakCollectionGraphics.cpp
+++ b/gui/graphics/PeakCollectionGraphics.cpp
@@ -14,6 +14,7 @@
#include "gui/graphics/PeakCollectionGraphics.h"
+#include "base/utils/ParallelFor.h"
#include "core/data/DataSet.h"
#include "core/loader/XFileHandler.h"
#include "core/peak/IntegrationRegion.h"
@@ -150,32 +151,42 @@ QVector<PeakCenterGraphic*> PeakCollectionGraphics::detectorSpots(std::size_t fr
return graphics;
}
-void PeakCollectionGraphics::getIntegrationMask(Eigen::MatrixXi& mask, std::size_t frame_idx)
+void PeakCollectionGraphics::getIntegrationMask(
+ Eigen::MatrixXi& mask, std::size_t frame_idx, bool thread_parallel /* = true */)
{
std::vector<PeakItem*> peak_items = _peak_model->root()->peakItems();
std::mutex mutex;
- for (auto* peak_item : peak_items) {
- ohkl::Peak3D* peak = peak_item->peak();
- double peak_end, bkg_begin, bkg_end;
- if (_preview_int_regions) {
- if (_int_params.region_type == ohkl::RegionType::VariableEllipsoid) {
- peak_end = _int_params.peak_end;
- bkg_begin = _int_params.bkg_begin;
- bkg_end = _int_params.bkg_end;
+
+ ohkl::parallel_for(peak_items.size(), [&](int start, int end) {
+ for (int idx = start; idx < end; ++idx) {
+ PeakItem* peak_item = peak_items.at(idx);
+ ohkl::Peak3D* peak = peak_item->peak();
+
+ double peak_end, bkg_begin, bkg_end;
+ if (_preview_int_regions) {
+ if (_int_params.region_type == ohkl::RegionType::VariableEllipsoid) {
+ peak_end = _int_params.peak_end;
+ bkg_begin = _int_params.bkg_begin;
+ bkg_end = _int_params.bkg_end;
+ } else {
+ peak_end = _int_params.fixed_peak_end;
+ bkg_begin = _int_params.fixed_bkg_begin;
+ bkg_end = _int_params.fixed_bkg_end;
+ }
} else {
- peak_end = _int_params.fixed_peak_end;
- bkg_begin = _int_params.fixed_bkg_begin;
- bkg_end = _int_params.fixed_bkg_end;
+ peak_end = peak_item->peak()->peakEnd();
+ bkg_begin = peak_item->peak()->bkgBegin();
+ bkg_end = peak_item->peak()->bkgEnd();
+ }
+ ohkl::IntegrationRegion region(
+ peak, peak_end, bkg_begin, bkg_end, _int_params.region_type);
+ {
+ std::lock_guard<std::mutex> lock(mutex);
+ region.updateMask(mask, frame_idx);
}
- } else {
- peak_end = peak_item->peak()->peakEnd();
- bkg_begin = peak_item->peak()->bkgBegin();
- bkg_end = peak_item->peak()->bkgEnd();
}
- ohkl::IntegrationRegion region(peak, peak_end, bkg_begin, bkg_end, _int_params.region_type);
- region.updateMask(mask, frame_idx);
- }
+ }, thread_parallel);
}
void PeakCollectionGraphics::getSinglePeakIntegrationMask(
diff --git a/gui/graphics/PeakCollectionGraphics.h b/gui/graphics/PeakCollectionGraphics.h
index eeee416ec..615577007 100644
--- a/gui/graphics/PeakCollectionGraphics.h
+++ b/gui/graphics/PeakCollectionGraphics.h
@@ -87,7 +87,8 @@ class PeakCollectionGraphics {
private:
//! Generate a mask of integration regions (a matrix of integers classifying pixels)
- void getIntegrationMask(Eigen::MatrixXi& mask, std::size_t frame_idx);
+ void getIntegrationMask(
+ Eigen::MatrixXi& mask, std::size_t frame_idx, bool thread_parallel = true);
//! Generate a mask for a single peak only
void getSinglePeakIntegrationMask(
ohkl::Peak3D* peak, Eigen::MatrixXi& mask, std::size_t frame_idx);
diff --git a/gui/items/PeakCollectionItem.cpp b/gui/items/PeakCollectionItem.cpp
index 0368caae5..b1bd8c0f2 100644
--- a/gui/items/PeakCollectionItem.cpp
+++ b/gui/items/PeakCollectionItem.cpp
@@ -15,10 +15,12 @@
#include "gui/items/PeakCollectionItem.h"
#include "base/geometry/ReciprocalVector.h"
+#include "base/utils/ParallelFor.h"
#include "core/data/DataSet.h"
#include "core/detector/Detector.h"
#include "core/instrument/Diffractometer.h"
#include "core/instrument/InstrumentState.h"
+#include "core/peak/Peak3D.h"
#include "core/raw/DataKeys.h"
#include "core/raw/MetaData.h"
#include "core/shape/PeakCollection.h"
@@ -28,6 +30,10 @@
#include "tables/crystal/MillerIndex.h"
#include "tables/crystal/UnitCell.h"
+#include <Eigen/Dense>
+
+#include <mutex>
+
PeakCollectionItem::PeakCollectionItem()
{
_peak_collection = nullptr;
@@ -46,7 +52,8 @@ PeakCollectionItem::PeakCollectionItem(const ohkl::PeakCollection* peak_collecti
}
}
-void PeakCollectionItem::setPeakCollection(const ohkl::PeakCollection* peak_collection)
+void PeakCollectionItem::setPeakCollection(
+ const ohkl::PeakCollection* peak_collection, bool thread_parallel /* = true */)
{
if (!peak_collection) {
throw std::runtime_error(
@@ -55,11 +62,19 @@ void PeakCollectionItem::setPeakCollection(const ohkl::PeakCollection* peak_coll
_peak_collection = peak_collection;
std::vector<ohkl::Peak3D*> peak_list = _peak_collection->getPeakList();
- _peak_items.clear();
- for (auto* peak : peak_list) {
- auto item = std::make_unique<PeakItem>(peak);
- _peak_items.push_back(std::move(item));
- }
+ std::mutex mutex;
+
+ ohkl::parallel_for(peak_list.size(), [&](int start, int end) {
+ for (int idx = start; idx < end; ++idx) {
+ ohkl::Peak3D* peak = peak_list.at(idx);
+ auto item = std::make_unique<PeakItem>(peak);
+
+ {
+ const std::lock_guard<std::mutex> lock(mutex);
+ _peak_items.push_back(std::move(item));
+ }
+ }
+ }, thread_parallel);
}
std::string PeakCollectionItem::name() const
diff --git a/gui/items/PeakCollectionItem.h b/gui/items/PeakCollectionItem.h
index 90f27e13f..8cd9911dc 100644
--- a/gui/items/PeakCollectionItem.h
+++ b/gui/items/PeakCollectionItem.h
@@ -34,7 +34,7 @@ class PeakCollectionItem {
public:
//! Set the peak collection
- void setPeakCollection(const ohkl::PeakCollection* peak_collection);
+ void setPeakCollection(const ohkl::PeakCollection* peak_collection, bool thread_parallel = true);
//! Retrieve the name of the collection if present
std::string name() const;
//! Retrieve the number of children of this item
--
GitLab
From 4eb73f82b3cfb2e3972e582c96267156b1cc7ea9 Mon Sep 17 00:00:00 2001
From: Zamaan Raza <z.raza@fz-juelich.de>
Date: Mon, 3 Jun 2024 10:23:42 +0200
Subject: [PATCH 16/22] Dump list of peaks sorted by symmetry relations
---
core/integration/IIntegrator.h | 2 +-
core/peak/Peak3D.cpp | 17 ++++----
core/statistics/MergedPeak.cpp | 11 +++++
core/statistics/MergedPeak.h | 3 ++
core/statistics/MergedPeakCollection.cpp | 11 +++++
core/statistics/MergedPeakCollection.h | 3 ++
core/statistics/PeakMerger.cpp | 17 ++++++++
core/statistics/PeakMerger.h | 3 ++
scripts/merge.py | 54 ++++++++++++++++++++++++
9 files changed, 112 insertions(+), 9 deletions(-)
create mode 100755 scripts/merge.py
diff --git a/core/integration/IIntegrator.h b/core/integration/IIntegrator.h
index fdc988718..ed30c9b71 100644
--- a/core/integration/IIntegrator.h
+++ b/core/integration/IIntegrator.h
@@ -32,7 +32,7 @@ enum class IntegratorType { PixelSum = 0, Gaussian, ISigma, Profile1D, Profile3D
/*! \addtogroup python_api
* @{*/
- //! Return value of all integrator compute methods. Contains any mutable data from integration.
+//! Return value of all integrator compute methods. Contains any mutable data from integration.
struct ComputeResult {
//! Construct an invalid result by default
ComputeResult();
diff --git a/core/peak/Peak3D.cpp b/core/peak/Peak3D.cpp
index 3da3d2519..000019610 100644
--- a/core/peak/Peak3D.cpp
+++ b/core/peak/Peak3D.cpp
@@ -483,14 +483,15 @@ std::string Peak3D::rejectionString() const
std::string Peak3D::toString() const
{
std::ostringstream oss;
- // h, k, l, x, y, frame, intensity, sigma
- oss << std::fixed << std::setw(5) << _hkl.h() << std::fixed << std::setw(5) << _hkl.k()
- << std::fixed << std::setw(5) << _hkl.l() << std::fixed << std::setw(10)
- << std::setprecision(2) << shape().center()[0] << std::fixed << std::setw(10)
- << std::setprecision(2) << shape().center()[1] << std::fixed << std::setw(10)
- << std::setprecision(2) << shape().center()[2] << std::fixed << std::setw(10)
- << std::setprecision(2) << correctedSumIntensity().value() << std::fixed << std::setw(10)
- << std::setprecision(2) << correctedSumIntensity().sigma();
+ // d, h, k, l, x, y, frame, intensity, sigma
+ oss << std::fixed << std::setw(8) << std::setprecision(3) << d()
+ << std::fixed << std::setw(5) << _hkl.h()
+ << std::fixed << std::setw(5) << _hkl.k() << std::fixed << std::setw(5) << _hkl.l()
+ << std::fixed << std::setw(10) << std::setprecision(2) << shape().center()[0]
+ << std::fixed << std::setw(10) << std::setprecision(2) << shape().center()[1]
+ << std::fixed << std::setw(10) << std::setprecision(2) << shape().center()[2]
+ << std::fixed << std::setw(14) << std::setprecision(2) << correctedSumIntensity().value()
+ << std::fixed << std::setw(10) << std::setprecision(2) << correctedSumIntensity().sigma();
return oss.str();
}
diff --git a/core/statistics/MergedPeak.cpp b/core/statistics/MergedPeak.cpp
index 716e55ee6..c1ee044e3 100644
--- a/core/statistics/MergedPeak.cpp
+++ b/core/statistics/MergedPeak.cpp
@@ -208,4 +208,15 @@ double MergedPeak::pValue() const
return gsl_cdf_chisq_P(x, k);
}
+std::string MergedPeak::toString() const
+{
+ std::string peaklist = "";
+ for (auto* peak : _peaks) {
+ peaklist += peak->toString();
+ peaklist += "\n";
+ }
+
+ return peaklist;
+}
+
} // namespace ohkl
diff --git a/core/statistics/MergedPeak.h b/core/statistics/MergedPeak.h
index a29ff0127..b74b5cf87 100644
--- a/core/statistics/MergedPeak.h
+++ b/core/statistics/MergedPeak.h
@@ -61,6 +61,9 @@ class MergedPeak {
//! split the merged peak randomly into two, for calculation of CC
std::pair<MergedPeak, MergedPeak> split(bool sum_intensity) const;
+ //! Return list of symmetry-related peaks as a string
+ std::string toString() const;
+
private:
//! Update the hkl that represents the set of equivalences.
void determineRepresentativeHKL();
diff --git a/core/statistics/MergedPeakCollection.cpp b/core/statistics/MergedPeakCollection.cpp
index 4f34bc132..261135552 100644
--- a/core/statistics/MergedPeakCollection.cpp
+++ b/core/statistics/MergedPeakCollection.cpp
@@ -208,4 +208,15 @@ double MergedPeakCollection::dMax() const
return _d_max;
}
+std::string MergedPeakCollection::toStringUnmerged() const
+{
+ std::string peaklist = "";
+ for (auto merged_peak : _merged_peak_set) {
+ peaklist += merged_peak.toString();
+ peaklist += "\n";
+ }
+
+ return peaklist;
+}
+
} // namespace ohkl
diff --git a/core/statistics/MergedPeakCollection.h b/core/statistics/MergedPeakCollection.h
index 6a0136eda..717099d9c 100644
--- a/core/statistics/MergedPeakCollection.h
+++ b/core/statistics/MergedPeakCollection.h
@@ -71,6 +71,9 @@ class MergedPeakCollection {
//! Get maximum d
double dMax() const;
+ //! Output a list of unmerged peaks as a string, sorted by symmetry relations
+ std::string toStringUnmerged() const;
+
SpaceGroup spaceGroup() const { return _group; };
private:
diff --git a/core/statistics/PeakMerger.cpp b/core/statistics/PeakMerger.cpp
index 6cfc291c0..2cde0d407 100644
--- a/core/statistics/PeakMerger.cpp
+++ b/core/statistics/PeakMerger.cpp
@@ -328,6 +328,23 @@ bool PeakMerger::saveStatistics(std::string filename)
return true;
}
+bool PeakMerger::savePeaks(std::string filename, bool sum_intensities) const
+{
+ std::string data;
+ if (sum_intensities)
+ data = _sum_merged_data->toStringUnmerged();
+ else
+ data = _profile_merged_data->toStringUnmerged();
+
+ std::fstream file(filename, std::ios::out);
+ if (!file.is_open())
+ return false;
+
+ file << data;
+ file.close();
+ return true;
+}
+
std::vector<double> PeakMerger::getFigureOfMerit(FigureOfMerit fom, IntegratorType integrator)
{
DataResolution* resolution;
diff --git a/core/statistics/PeakMerger.h b/core/statistics/PeakMerger.h
index 378f05d70..054444838 100644
--- a/core/statistics/PeakMerger.h
+++ b/core/statistics/PeakMerger.h
@@ -97,6 +97,9 @@ class PeakMerger {
//! Saves the shell information to file.
bool saveStatistics(std::string filename);
+ //! Save peaks grouped by symmetry relations to file
+ bool savePeaks(std::string filename, bool sum_intensities) const;
+
//! Return a vector containing the resolution-dependent figure of merit
std::vector<double> getFigureOfMerit(FigureOfMerit fom, IntegratorType integrator);
diff --git a/scripts/merge.py b/scripts/merge.py
new file mode 100755
index 000000000..dcef23667
--- /dev/null
+++ b/scripts/merge.py
@@ -0,0 +1,54 @@
+#!/usr/bin/python
+## ***********************************************************************************************
+##
+## OpenHKL: data reduction for single-crystal diffraction
+##
+##! @file test/python/TestWorkFlow.py
+##! @brief Test ...
+##!
+##! @homepage https://openhkl.org
+##! @license GNU General Public License v3 or higher (see COPYING)
+##! @copyright Institut Laue-Langevin and Forschungszentrum Jülich GmbH 2016-
+##! @authors see CITATION, MAINTAINER
+##
+## ***********************************************************************************************
+
+import sys
+import numpy as np
+from pathlib import Path
+
+lib_dir = "@SWIG_INSTALL_PATH@" # Path to pyohkl.py
+
+sys.path.append(lib_dir)
+import pyohkl as ohkl
+
+file = Path('@CMAKE_BINARY_DIR@/test/data/Trypsin-pxsum.ohkl')
+experiment = 'Trypsin'
+diffractometer = 'BioDiff'
+data_name = 'Scan I'
+space_group = "P 21 21 21"
+found_peaks_name = 'found'
+predicted_peaks_name = 'predicted'
+
+expt = ohkl.Experiment(experiment, diffractometer)
+expt.loadFromFile(str(file))
+data = expt.getData(data_name)
+found_peaks = expt.getPeakCollection(found_peaks_name)
+predicted_peaks = expt.getPeakCollection(predicted_peaks_name)
+
+
+# Merge parameters
+merger = expt.peakMerger()
+params = merger.parameters()
+merger.reset()
+params.d_min = 1.5
+params.d_max = 50.0
+params.n_shells = 10
+params.friedel = True
+merger.addPeakCollection(predicted_peaks)
+merger.setSpaceGroup(ohkl.SpaceGroup(space_group))
+merger.mergePeaks()
+merger.computeQuality()
+print(merger.summary())
+
+merger.savePeaks("unmerged_sum.txt", True)
--
GitLab
From 602f918184f12c76a0fdf87b14b8f273fcba4ef9 Mon Sep 17 00:00:00 2001
From: Zamaan Raza <z.raza@fz-juelich.de>
Date: Tue, 4 Jun 2024 14:51:16 +0200
Subject: [PATCH 17/22] Fix GUI hang on integration
---
base/utils/ParallelFor.h | 2 --
gui/items/PeakCollectionItem.cpp | 16 +++++-----------
2 files changed, 5 insertions(+), 13 deletions(-)
diff --git a/base/utils/ParallelFor.h b/base/utils/ParallelFor.h
index 3245fd170..8a3a5523c 100644
--- a/base/utils/ParallelFor.h
+++ b/base/utils/ParallelFor.h
@@ -20,8 +20,6 @@
#include <thread>
#include <vector>
-#include <iostream>
-
namespace ohkl {
/// @param[in] nb_elements : size of your for loop
diff --git a/gui/items/PeakCollectionItem.cpp b/gui/items/PeakCollectionItem.cpp
index b1bd8c0f2..d423c19fb 100644
--- a/gui/items/PeakCollectionItem.cpp
+++ b/gui/items/PeakCollectionItem.cpp
@@ -62,19 +62,13 @@ void PeakCollectionItem::setPeakCollection(
_peak_collection = peak_collection;
std::vector<ohkl::Peak3D*> peak_list = _peak_collection->getPeakList();
- std::mutex mutex;
- ohkl::parallel_for(peak_list.size(), [&](int start, int end) {
- for (int idx = start; idx < end; ++idx) {
- ohkl::Peak3D* peak = peak_list.at(idx);
- auto item = std::make_unique<PeakItem>(peak);
+ _peak_items.clear();
- {
- const std::lock_guard<std::mutex> lock(mutex);
- _peak_items.push_back(std::move(item));
- }
- }
- }, thread_parallel);
+ for (auto* peak : peak_list) {
+ auto item = std::make_unique<PeakItem>(peak);
+ _peak_items.push_back(std::move(item));
+ }
}
std::string PeakCollectionItem::name() const
--
GitLab
From 065f3d4c3869dd0efc3519ec4c9033002b625f54 Mon Sep 17 00:00:00 2001
From: Zamaan Raza <z.raza@fz-juelich.de>
Date: Tue, 4 Jun 2024 14:52:22 +0200
Subject: [PATCH 18/22] Always save DataReaderParameters section of YAML file
---
core/experiment/Experiment.cpp | 4 ++
core/experiment/Experiment.h | 7 +++
gui/models/Session.cpp | 78 +++++++++++++++++-----------------
3 files changed, 50 insertions(+), 39 deletions(-)
diff --git a/core/experiment/Experiment.cpp b/core/experiment/Experiment.cpp
index d664dc7e7..80ad738c2 100644
--- a/core/experiment/Experiment.cpp
+++ b/core/experiment/Experiment.cpp
@@ -34,6 +34,7 @@
#include "core/instrument/InstrumentStateSet.h"
#include "core/instrument/Monochromator.h"
#include "core/raw/DataKeys.h"
+#include "core/loader/IDataReader.h"
#include "core/raw/MetaData.h"
#include "core/shape/PeakCollection.h"
#include "core/shape/PeakFilter.h"
@@ -72,6 +73,7 @@ Experiment::Experiment() : _diffractometer(nullptr), _strategy(false)
_peak_merger = std::make_unique<PeakMerger>();
_shape_params = std::make_shared<ShapeModelParameters>();
+ _data_reader_params = std::make_unique<DataReaderParameters>();
}
Experiment::~Experiment() = default;
@@ -106,6 +108,7 @@ void Experiment::readFromYaml(const std::string& filename)
{
ohklLog(Level::Info, "Experiment::readFromYaml: ", filename);
ExperimentYAML yaml(filename);
+ yaml.grabDataReaderParameters(_data_reader_params.get());
yaml.grabPeakFinderParameters(_peak_finder->parameters());
yaml.grabAutoindexerParameters(_auto_indexer->parameters());
yaml.grabShapeParameters(_shape_params.get());
@@ -118,6 +121,7 @@ void Experiment::saveToYaml(const std::string& filename)
{
ohklLog(Level::Info, "Experiment::saveToYaml: ", filename);
ExperimentYAML yaml(filename);
+ yaml.setDataReaderParameters(_data_reader_params.get());
yaml.setPeakFinderParameters(_peak_finder->parameters());
yaml.setAutoindexerParameters(_auto_indexer->parameters());
yaml.setShapeParameters(_shape_params.get());
diff --git a/core/experiment/Experiment.h b/core/experiment/Experiment.h
index bb1981b40..bc3d44d51 100644
--- a/core/experiment/Experiment.h
+++ b/core/experiment/Experiment.h
@@ -43,6 +43,7 @@ class Refiner;
class ShapeHandler;
class ShapeModel;
class UnitCellHandler;
+struct DataReaderParameters;
struct ShapeModelParameters;
using DataMap = std::map<std::string, sptrDataSet>;
@@ -91,6 +92,9 @@ class Experiment {
//! get pointer to ShapeModel parameters
std::shared_ptr<ShapeModelParameters> shapeParameters() const { return _shape_params; };
+ //! Get a pointer to the data reader parameters
+ DataReaderParameters* dataReaderParameters() const {return _data_reader_params.get(); };
+
// Data handler
//! Get the map of data sets from the handler
const DataMap* getDataMap() const;
@@ -296,6 +300,9 @@ class Experiment {
// ShapeModel parameters, since there is no experiment-level container
std::shared_ptr<ShapeModelParameters> _shape_params;
+ // Data reader parameters, since there is nowhere else to store this
+ std::unique_ptr<DataReaderParameters> _data_reader_params;
+
bool _strategy;
};
diff --git a/gui/models/Session.cpp b/gui/models/Session.cpp
index fbf3db7a1..013d83cb2 100644
--- a/gui/models/Session.cpp
+++ b/gui/models/Session.cpp
@@ -366,54 +366,54 @@ bool Session::loadRawData(bool single_file /* = false */)
// Get metadata from readme file, then edit them in dialog.
const QStringList& extant_dataset_names = currentProject()->getDataNames();
- ohkl::DataReaderParameters parameters;
- QString yml_file = QString::fromStdString(currentProject()->experiment()->name() + ".yml");
+ currentProject()->setDirectory(path);
+ ohkl::Experiment* exp = currentProject()->experiment();
+
+ ohkl::DataReaderParameters* parameters = exp->dataReaderParameters();
+ QString yml_file = QString::fromStdString(exp->name() + ".yml");
QString yml_path = QDir(path).filePath(yml_file);
QFileInfo info(yml_path);
if (info.exists() && info.isFile())
- parameters.loadFromYAML(yml_path.toStdString());
- ImageReaderDialog dialog(filenames, ¶meters, ohkl::DataFormat::RAW);
+ parameters->loadFromYAML(yml_path.toStdString());
+ ImageReaderDialog dialog(filenames, parameters, ohkl::DataFormat::RAW);
dialog.setWindowTitle("Raw data parameters");
if (single_file)
dialog.setSingleImageMode();
if (!dialog.exec())
return false;
-
- currentProject()->setDirectory(path);
- ohkl::Experiment* exp = currentProject()->experiment();
- parameters = dialog.dataReaderParameters();
+ *parameters = dialog.dataReaderParameters();
// Transfer metadata to diffractometer.
ohkl::Detector* detector = exp->getDiffractometer()->detector();
- detector->setBaseline(parameters.baseline);
- detector->setGain(parameters.gain);
+ detector->setBaseline(parameters->baseline);
+ detector->setGain(parameters->gain);
// Transfer metadata to dataset, and load the raw data.
if (single_file) {
const std::shared_ptr<ohkl::DataSet> dataset{std::make_shared<ohkl::SingleFrame>(
- parameters.dataset_name, exp->getDiffractometer())};
- dataset->setImageReaderParameters(parameters);
+ parameters->dataset_name, exp->getDiffractometer())};
+ dataset->setImageReaderParameters(*parameters);
dataset->addRawFrame(filenames[0].toStdString());
dataset->finishRead();
- parameters.cols = dataset->nCols();
- parameters.rows = dataset->nRows();
+ parameters->cols = dataset->nCols();
+ parameters->rows = dataset->nRows();
exp->addData(dataset);
} else {
const std::shared_ptr<ohkl::DataSet> dataset{
- std::make_shared<ohkl::DataSet>(parameters.dataset_name, exp->getDiffractometer())};
- dataset->setImageReaderParameters(parameters);
+ std::make_shared<ohkl::DataSet>(parameters->dataset_name, exp->getDiffractometer())};
+ dataset->setImageReaderParameters(*parameters);
for (const auto& filename : filenames)
dataset->addRawFrame(filename.toStdString());
dataset->finishRead();
- parameters.cols = dataset->nCols();
- parameters.rows = dataset->nRows();
+ parameters->cols = dataset->nCols();
+ parameters->rows = dataset->nRows();
exp->addData(dataset);
}
ohkl::ExperimentYAML yaml(yml_path.toStdString());
- yaml.setDataReaderParameters(¶meters);
+ yaml.setDataReaderParameters(parameters);
yaml.writeFile(yml_path.toStdString());
onDataChanged();
} catch (std::exception& e) {
@@ -431,7 +431,7 @@ bool Session::loadTiffData(bool single_file /* = false */)
ohkl::Experiment* exp = currentProject()->experiment();
ohkl::Detector* detector = exp->getDiffractometer()->detector();
- ohkl::DataReaderParameters params;
+ ohkl::DataReaderParameters* params = exp->dataReaderParameters();
try {
// Get input filenames from dialog.
@@ -450,9 +450,9 @@ bool Session::loadTiffData(bool single_file /* = false */)
QString yml_path = QDir(path).filePath(yml_file);
QFileInfo info(yml_path);
if (info.exists() && info.isFile())
- params.loadFromYAML(yml_path.toStdString());
+ params->loadFromYAML(yml_path.toStdString());
ImageReaderDialog dialog(
- filenames, static_cast<ohkl::DataReaderParameters*>(¶ms), ohkl::DataFormat::TIFF);
+ filenames, static_cast<ohkl::DataReaderParameters*>(params), ohkl::DataFormat::TIFF);
dialog.setWindowTitle("Tiff data parameters");
if (single_file)
@@ -462,21 +462,21 @@ bool Session::loadTiffData(bool single_file /* = false */)
return false;
currentProject()->setDirectory(path);
- params = dialog.dataReaderParameters();
- detector->setBaseline(params.baseline);
- detector->setGain(params.gain);
+ *params = dialog.dataReaderParameters();
+ detector->setBaseline(params->baseline);
+ detector->setGain(params->gain);
if (single_file) {
const std::shared_ptr<ohkl::DataSet> dataset{
- std::make_shared<ohkl::SingleFrame>(params.dataset_name, exp->getDiffractometer())};
- dataset->setImageReaderParameters(params);
+ std::make_shared<ohkl::SingleFrame>(params->dataset_name, exp->getDiffractometer())};
+ dataset->setImageReaderParameters(*params);
dataset->addTiffFrame(filenames[0].toStdString());
dataset->finishRead();
exp->addData(dataset);
} else {
const std::shared_ptr<ohkl::DataSet> dataset{
- std::make_shared<ohkl::DataSet>(params.dataset_name, exp->getDiffractometer())};
- dataset->setImageReaderParameters(params);
+ std::make_shared<ohkl::DataSet>(params->dataset_name, exp->getDiffractometer())};
+ dataset->setImageReaderParameters(*params);
for (const auto& filename : filenames)
dataset->addTiffFrame(filename.toStdString());
dataset->finishRead();
@@ -499,7 +499,7 @@ bool Session::loadPlainTextData(bool single_file /* = false */)
ohkl::Experiment* exp = currentProject()->experiment();
ohkl::Detector* detector = exp->getDiffractometer()->detector();
- ohkl::DataReaderParameters params;
+ ohkl::DataReaderParameters* params = exp->dataReaderParameters();
try {
// Get input filenames from dialog.
@@ -518,9 +518,9 @@ bool Session::loadPlainTextData(bool single_file /* = false */)
QString yml_path = QDir(path).filePath(yml_file);
QFileInfo info(yml_path);
if (info.exists() && info.isFile())
- params.loadFromYAML(yml_path.toStdString());
+ params->loadFromYAML(yml_path.toStdString());
ImageReaderDialog dialog(
- filenames, static_cast<ohkl::DataReaderParameters*>(¶ms),
+ filenames, static_cast<ohkl::DataReaderParameters*>(params),
ohkl::DataFormat::PLAINTEXT);
dialog.setWindowTitle("Plain text data parameters");
@@ -531,21 +531,21 @@ bool Session::loadPlainTextData(bool single_file /* = false */)
return false;
currentProject()->setDirectory(path);
- params = dialog.dataReaderParameters();
- detector->setBaseline(params.baseline);
- detector->setGain(params.gain);
+ *params = dialog.dataReaderParameters();
+ detector->setBaseline(params->baseline);
+ detector->setGain(params->gain);
if (single_file) {
const std::shared_ptr<ohkl::DataSet> dataset{
- std::make_shared<ohkl::SingleFrame>(params.dataset_name, exp->getDiffractometer())};
- dataset->setImageReaderParameters(params);
+ std::make_shared<ohkl::SingleFrame>(params->dataset_name, exp->getDiffractometer())};
+ dataset->setImageReaderParameters(*params);
dataset->addPlainTextFrame(filenames[0].toStdString());
dataset->finishRead();
exp->addData(dataset);
} else {
const std::shared_ptr<ohkl::DataSet> dataset{
- std::make_shared<ohkl::DataSet>(params.dataset_name, exp->getDiffractometer())};
- dataset->setImageReaderParameters(params);
+ std::make_shared<ohkl::DataSet>(params->dataset_name, exp->getDiffractometer())};
+ dataset->setImageReaderParameters(*params);
for (const auto& filename : filenames)
dataset->addPlainTextFrame(filename.toStdString());
dataset->finishRead();
--
GitLab
From fc9fd6772bd1142e1df17c0ed333f1dce32f1ad9 Mon Sep 17 00:00:00 2001
From: Zamaan Raza <z.raza@fz-juelich.de>
Date: Thu, 6 Jun 2024 08:43:17 +0200
Subject: [PATCH 19/22] Disable full integrator tests
---
cmake/CTestCustom.cmake | 1 +
1 file changed, 1 insertion(+)
diff --git a/cmake/CTestCustom.cmake b/cmake/CTestCustom.cmake
index 626014259..45f0e4844 100644
--- a/cmake/CTestCustom.cmake
+++ b/cmake/CTestCustom.cmake
@@ -1,4 +1,5 @@
# Tests specified in the following variable will not be executed
set(CTEST_CUSTOM_TESTS_IGNORE
TestProfile3DIntegration
+ TestPxsumIntegration
)
--
GitLab
From a8846386c1b998623760f1896a94d83dcbe0f71a Mon Sep 17 00:00:00 2001
From: Zamaan Raza <z.raza@fz-juelich.de>
Date: Thu, 6 Jun 2024 09:00:19 +0200
Subject: [PATCH 20/22] CMake variable to toggle integration tests (off by
default)
---
CMakeLists.txt | 5 ++++-
1 file changed, 4 insertions(+), 1 deletion(-)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 8c2b569f8..58c5e8c22 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -33,6 +33,7 @@ option(BUILD_PDF_DOCUMENTATION "Build Sphinx documentation using LaTeX" ON)
option(OHKL_TIDY "Upon 'make', clang-tidy will print warnings as specified in .clang-tidy" OFF)
option(OHKL_FULL_WORKFLOW_TEST "Add full workflow test to ctest" OFF)
option(BUILD_WITH_QT6 "Build using QT6" OFF)
+option(TEST_INTEGRATORS "Build intensive integrator tests" OFF)
# Path to test data directory must provided via -DOHKL_TESTDATA_DIR or from cache.
set(OHKL_TESTDATA_URL https://computing.mlz-garching.de/download/testdata/nsx)
@@ -226,7 +227,9 @@ if(OHKL_PYTHON AND NOT MSVC AND BUILD_TESTING)
endif()
# Ignore tests specified in CTestCustom.cmake
-configure_file(${CMAKE_SOURCE_DIR}/cmake/CTestCustom.cmake ${CMAKE_BINARY_DIR})
+if (NOT TEST_INTEGRATORS)
+ configure_file(${CMAKE_SOURCE_DIR}/cmake/CTestCustom.cmake ${CMAKE_BINARY_DIR})
+endif()
###################################################################################################
# Build documentation
--
GitLab
From 77938ddd21402f93240ec9fdb16ae8edde79180a Mon Sep 17 00:00:00 2001
From: Zamaan Raza <z.raza@fz-juelich.de>
Date: Wed, 24 Jul 2024 13:01:11 +0200
Subject: [PATCH 21/22] Code review changes
---
core/data/DataSet.h | 2 +-
core/experiment/ExperimentYAML.cpp | 10 ++++++----
core/loader/XFileHandler.cpp | 12 ++++++------
core/peak/Peak3D.cpp | 1 -
core/peak/Qs2Events.cpp | 4 ++--
5 files changed, 15 insertions(+), 14 deletions(-)
diff --git a/core/data/DataSet.h b/core/data/DataSet.h
index d6c2201d8..e4d9d20e2 100644
--- a/core/data/DataSet.h
+++ b/core/data/DataSet.h
@@ -155,7 +155,7 @@ class DataSet {
void clearHistograms();
//! Maximum per pixel count for whole DataSet
double maxCount();
- //! getting number of available
+ //! Get number of available histograms
size_t getNumberHistograms() { return _histograms.size(); }
//! accessing created histograms
diff --git a/core/experiment/ExperimentYAML.cpp b/core/experiment/ExperimentYAML.cpp
index 4ede768af..6cf1be920 100644
--- a/core/experiment/ExperimentYAML.cpp
+++ b/core/experiment/ExperimentYAML.cpp
@@ -2,8 +2,8 @@
//
// OpenHKL: data reduction for single crystal diffraction
//
-//! @file base/parser/ExperimentYAML.cpp
-//! @brief Defines function eigenToVector
+//! @file core/experiment/ExperimentYAML.cpp
+//! @brief Implements class ExperimentYAML
//!
//! @homepage https://openhkl.org
//! @license GNU General Public License v3 or higher (see COPYING)
@@ -124,7 +124,6 @@ void ExperimentYAML::grabIntegrationParameters(IntegrationParameters* params)
params->fixed_peak_end = getNode<double>(branch, "fixed_peak_end");
params->fixed_bkg_begin = getNode<double>(branch, "fixed_bkg_begin");
params->fixed_bkg_end = getNode<double>(branch, "fixed_bkg_end");
- params->max_strength = getNode<double>(branch, "max_strength");
params->fit_center = getNode<bool>(branch, "fit_center");
params->fit_cov = getNode<bool>(branch, "fit_cov");
params->integrator_type = static_cast<IntegratorType>(getNode<int>(branch, "integrator_type"));
@@ -138,6 +137,8 @@ void ExperimentYAML::grabIntegrationParameters(IntegrationParameters* params)
params->max_strength = getNode<double>(branch, "max_strength");
params->use_max_d = getNode<bool>(branch, "use_max_d");
params->max_d = getNode<double>(branch, "max_d");
+ params->use_max_width = getNode<bool>(branch, "use_max_width");
+ params->max_width = getNode<double>(branch, "max_width");
params->discard_saturated = getNode<bool>(branch, "discard_saturated");
params->max_counts = getNode<double>(branch, "max_counts");
}
@@ -156,7 +157,6 @@ void ExperimentYAML::setIntegrationParameters(IntegrationParameters* params)
int_node["fixed_peak_end"] = params->fixed_peak_end;
int_node["fixed_bkg_begin"] = params->fixed_bkg_begin;
int_node["fixed_bkg_end"] = params->fixed_bkg_end;
- int_node["max_strength"] = params->max_strength;
int_node["fit_center"] = params->fit_center;
int_node["fit_cov"] = params->fit_cov;
int_node["integrator_type"] = static_cast<int>(params->integrator_type);
@@ -170,6 +170,8 @@ void ExperimentYAML::setIntegrationParameters(IntegrationParameters* params)
int_node["max_strength"] = params->max_strength;
int_node["use_max_d"] = params->use_max_d;
int_node["max_d"] = params->max_d;
+ int_node["use_max_width"] = params->use_max_width;
+ int_node["max_width"] = params->max_width;
int_node["discard_saturated"] = params->discard_saturated;
int_node["max_counts"] = params->max_counts;
}
diff --git a/core/loader/XFileHandler.cpp b/core/loader/XFileHandler.cpp
index 671522665..e72cc05d3 100644
--- a/core/loader/XFileHandler.cpp
+++ b/core/loader/XFileHandler.cpp
@@ -51,24 +51,24 @@ void XFileHandler::readXFile(int frame)
_a(1, 0) = std::stod(tokens[1]);
_a(2, 0) = std::stod(tokens[2]);
_u(0, 0) = std::stod(tokens[3]);
- _u(1, 0) = std::stod(tokens[3]);
- _u(2, 0) = std::stod(tokens[4]);
+ _u(1, 0) = std::stod(tokens[4]);
+ _u(2, 0) = std::stod(tokens[5]);
std::getline(ifs, line);
tokens = tokenize(line);
_a(0, 1) = std::stod(tokens[0]);
_a(1, 1) = std::stod(tokens[1]);
_a(2, 1) = std::stod(tokens[2]);
_u(0, 1) = std::stod(tokens[3]);
- _u(2, 1) = std::stod(tokens[3]);
- _u(1, 1) = std::stod(tokens[4]);
+ _u(2, 1) = std::stod(tokens[4]);
+ _u(1, 1) = std::stod(tokens[5]);
std::getline(ifs, line);
tokens = tokenize(line);
_a(0, 2) = std::stod(tokens[0]);
_a(1, 2) = std::stod(tokens[1]);
_a(2, 2) = std::stod(tokens[2]);
_u(0, 2) = std::stod(tokens[3]);
- _u(1, 2) = std::stod(tokens[3]);
- _u(2, 2) = std::stod(tokens[4]);
+ _u(1, 2) = std::stod(tokens[4]);
+ _u(2, 2) = std::stod(tokens[5]);
// Other metadata (line 5)
std::getline(ifs, line);
diff --git a/core/peak/Peak3D.cpp b/core/peak/Peak3D.cpp
index 000019610..27368e1eb 100644
--- a/core/peak/Peak3D.cpp
+++ b/core/peak/Peak3D.cpp
@@ -247,7 +247,6 @@ void Peak3D::updateIntegration(
_regionType = regionType;
_rockingCurve = result.rocking_curve;
- _rockingCurve = result.rocking_curve;
if (result.integrator_type == IntegratorType::PixelSum)
_meanBkgGradient = result.bkg_gradient;
diff --git a/core/peak/Qs2Events.cpp b/core/peak/Qs2Events.cpp
index bda9b2666..3d18620da 100644
--- a/core/peak/Qs2Events.cpp
+++ b/core/peak/Qs2Events.cpp
@@ -43,7 +43,7 @@ std::vector<std::pair<MillerIndex, DetectorEvent>> algo::qMap2Events(
sptrProgressHandler handler /* = nullptr */, bool thread_parallel /* = true */)
{
ohklLog(
- Level::Debug, "algo::Qs2Events::qVectorList2Events: processing ", sample_qs.size(),
+ Level::Debug, "algo::Qs2Events::qMap2Events: processing ", sample_qs.size(),
" q-vectors");
std::vector<std::pair<MillerIndex, DetectorEvent>> events;
@@ -75,7 +75,7 @@ std::vector<std::pair<MillerIndex, DetectorEvent>> algo::qMap2Events(
if (handler)
handler->setProgress(100);
ohklLog(
- Level::Debug, "algo::Qs2Events::qVectorList2Events: finished; generated ", events.size(),
+ Level::Debug, "algo::Qs2Events::qMap2Events: finished; generated ", events.size(),
" events");
return events;
}
--
GitLab
From a0c9bc6d3349285494147b3d8a8d2579d3c6cbc9 Mon Sep 17 00:00:00 2001
From: Zamaan Raza <z.raza@fz-juelich.de>
Date: Wed, 24 Jul 2024 14:52:07 +0200
Subject: [PATCH 22/22] Add integration maximum width to SubframeFindPeaks
---
gui/subframe_find/SubframeFindPeaks.cpp | 21 ++++++++++++++++++++
gui/subframe_find/SubframeFindPeaks.h | 3 +++
gui/subframe_integrate/SubframeIntegrate.cpp | 2 +-
gui/subframe_shapes/SubframeShapes.cpp | 2 +-
4 files changed, 26 insertions(+), 2 deletions(-)
diff --git a/gui/subframe_find/SubframeFindPeaks.cpp b/gui/subframe_find/SubframeFindPeaks.cpp
index 0df88de8b..5591cb4df 100644
--- a/gui/subframe_find/SubframeFindPeaks.cpp
+++ b/gui/subframe_find/SubframeFindPeaks.cpp
@@ -216,6 +216,23 @@ void SubframeFindPeaks::setIntegrateUp()
"Background end",
"(" + QString(QChar(0x03C3)) + ") - scaling factor for upper limit of background");
+ _use_max_width = new QGroupBox("Maximum width for integration");
+ _use_max_width->setAlignment(Qt::AlignLeft);
+ _use_max_width->setCheckable(true);
+ _use_max_width->setChecked(false);
+ _use_max_width->setToolTip("Skip integration of peaks spanning too many images");
+
+ _max_width = new SafeSpinBox();
+ _max_width->setMaximum(100);
+
+ QLabel* label = new QLabel("Maximum width");
+ label->setToolTip("Maximum width for peak to be integrated");
+ QGridLayout* grid = new QGridLayout();
+ _use_max_width->setLayout(grid);
+ grid->addWidget(label, 0, 0, 1, 1);
+ grid->addWidget(_max_width, 0, 1, 1, 1);
+ f.addWidget(_use_max_width);
+
_gradient_check = f.addCheckBox(
"Compute gradient", "Compute mean gradient and sigma of background region", 1);
@@ -450,6 +467,8 @@ void SubframeFindPeaks::grabIntegrationParameters()
_bkg_end->setValue(params->fixed_bkg_end);
}
+ _use_max_width->setChecked(params->use_max_width);
+ _max_width->setValue(params->max_width);
_gradient_check->setChecked(params->use_gradient);
_fft_gradient_check->setChecked(params->fft_gradient);
_gradient_kernel->setCurrentIndex(static_cast<int>(params->gradient_type));
@@ -472,6 +491,8 @@ void SubframeFindPeaks::setIntegrationParameters()
params->fixed_bkg_begin = _bkg_begin->value();
params->fixed_bkg_end = _bkg_end->value();
}
+ params->use_max_width = _use_max_width->isChecked();
+ params->max_width = _max_width->value();
params->use_gradient = _gradient_check->isChecked();
params->fft_gradient = _fft_gradient_check->isChecked();
params->gradient_type = static_cast<ohkl::GradientKernel>(_gradient_kernel->currentIndex());
diff --git a/gui/subframe_find/SubframeFindPeaks.h b/gui/subframe_find/SubframeFindPeaks.h
index d9a1d3945..2d57304bc 100644
--- a/gui/subframe_find/SubframeFindPeaks.h
+++ b/gui/subframe_find/SubframeFindPeaks.h
@@ -28,6 +28,7 @@ class PeakTableView;
class PeakViewWidget;
class QCheckBox;
class QComboBox;
+class QGroupBox;
class QPushButton;
class QSplitter;
class QTableWidget;
@@ -128,6 +129,8 @@ class SubframeFindPeaks : public QWidget {
SafeDoubleSpinBox* _peak_end;
SafeDoubleSpinBox* _bkg_begin;
SafeDoubleSpinBox* _bkg_end;
+ SafeSpinBox* _max_width;
+ QGroupBox* _use_max_width;
QCheckBox* _gradient_check;
QCheckBox* _fft_gradient_check;
QComboBox* _gradient_kernel;
diff --git a/gui/subframe_integrate/SubframeIntegrate.cpp b/gui/subframe_integrate/SubframeIntegrate.cpp
index 2220098dc..8921edabf 100644
--- a/gui/subframe_integrate/SubframeIntegrate.cpp
+++ b/gui/subframe_integrate/SubframeIntegrate.cpp
@@ -411,7 +411,7 @@ void SubframeIntegrate::setIntegrateUp()
_use_max_width->setToolTip("Skip integration of peaks spanning too many images");
_max_width = new SafeSpinBox();
- _max_width->setMaximum(20);
+ _max_width->setMaximum(100);
label = new QLabel("Maximum width");
label->setToolTip("Maximum width for peak to be integrated");
diff --git a/gui/subframe_shapes/SubframeShapes.cpp b/gui/subframe_shapes/SubframeShapes.cpp
index 15d52f9eb..c80069718 100644
--- a/gui/subframe_shapes/SubframeShapes.cpp
+++ b/gui/subframe_shapes/SubframeShapes.cpp
@@ -184,7 +184,7 @@ void SubframeShapes::setInputUp()
_min_strength->setMaximum(10000);
_min_d->setMaximum(1000);
_max_d->setMaximum(10000);
- _max_width->setMaximum(20);
+ _max_width->setMaximum(100);
_peak_end->setMinimum(0.1);
_peak_end->setMaximum(50);
_peak_end->setSingleStep(0.1);
--
GitLab