diff --git a/Device/Analyze/Fourier.cpp b/Device/Analyze/Fourier.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..bb8cd39a6cbf3628dfbb53979e48353dd35e4c51
--- /dev/null
+++ b/Device/Analyze/Fourier.cpp
@@ -0,0 +1,38 @@
+// ************************************************************************************************
+//
+// BornAgain: simulate and fit reflection and scattering
+//
+//! @file Device/Analyze/Fourier.cpp
+//! @brief Implements Fourier transform function in namespace Analyze.
+//!
+//! @homepage http://www.bornagainproject.org
+//! @license GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2018
+//! @authors Scientific Computing Group at MLZ (see CITATION, AUTHORS)
+//
+// ************************************************************************************************
+
+#include "Device/Analyze/Fourier.h"
+#include "Base/Math/FourierTransform.h"
+#include "Device/Data/DataUtil.h"
+#include "Device/Data/Datafield.h"
+
+namespace {
+
+std::vector<std::vector<double>> FT2DArray(const std::vector<std::vector<double>>& signal)
+{
+ FourierTransform ft;
+ std::vector<std::vector<double>> result;
+ ft.fft(signal, result);
+ ft.fftshift(result); // low frequency to center of array
+ return result;
+}
+
+} // namespace
+
+Datafield Analyze::createFFT(const Datafield& data)
+{
+ auto array_2d = DataUtil::create2DArrayfromDatafield(data);
+ auto fft_array_2d = FT2DArray(array_2d);
+ return DataUtil::vecvecToDatafield(fft_array_2d);
+}
diff --git a/Device/Analyze/Fourier.h b/Device/Analyze/Fourier.h
new file mode 100644
index 0000000000000000000000000000000000000000..d4101cd3349a08529879dd4375309cb6723e6636
--- /dev/null
+++ b/Device/Analyze/Fourier.h
@@ -0,0 +1,30 @@
+// ************************************************************************************************
+//
+// BornAgain: simulate and fit reflection and scattering
+//
+//! @file Device/Analyze/Fourier.h
+//! @brief Defines Fourier transform function in namespace Analyze.
+//!
+//! @homepage http://www.bornagainproject.org
+//! @license GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2018
+//! @authors Scientific Computing Group at MLZ (see CITATION, AUTHORS)
+//
+// ************************************************************************************************
+
+#ifdef SWIG
+#error no need to expose this header to Swig
+#endif // SWIG
+#ifndef BORNAGAIN_DEVICE_ANALYZE_FOURIER_H
+#define BORNAGAIN_DEVICE_ANALYZE_FOURIER_H
+
+class Datafield;
+
+namespace Analyze {
+
+//! Creates Fourier Transform (Datafield format) of intensity map (Datafield format).
+Datafield createFFT(const Datafield& data);
+
+} // namespace Analyze
+
+#endif // BORNAGAIN_DEVICE_ANALYZE_FOURIER_H
diff --git a/Device/Analyze/Peaks.cpp b/Device/Analyze/Peaks.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a7b216d34ccfa4bc28a6e65775d3f2b9736feb16
--- /dev/null
+++ b/Device/Analyze/Peaks.cpp
@@ -0,0 +1,50 @@
+// ************************************************************************************************
+//
+// BornAgain: simulate and fit reflection and scattering
+//
+//! @file Device/Analyze/Peaks.cpp
+//! @brief Implements function FindPeaks namespace Analyze.
+//!
+//! @homepage http://www.bornagainproject.org
+//! @license GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2018
+//! @authors Scientific Computing Group at MLZ (see CITATION, AUTHORS)
+//
+// ************************************************************************************************
+
+#include "Device/Analyze/Peaks.h"
+#include "Base/Axis/Scale.h"
+#include "Device/Data/DataUtil.h"
+#include "Device/Data/Datafield.h"
+#include <tspectrum.h> // third-party code, extracted from CERN ROOT (class TSpectrum2)
+
+std::vector<std::pair<double, double>>
+Analyze::FindPeaks(const Datafield& data, double sigma, const std::string& option, double threshold)
+{
+ const std::vector<std::vector<double>> arr = DataUtil::createVector2D(data);
+ tspectrum::Spectrum2D spec;
+ const auto peaks = spec.find_peaks(arr, sigma, option, threshold);
+
+ // coordinates of peaks in histogram axes units
+ std::vector<std::pair<double, double>> result;
+
+ for (const auto& p : peaks) {
+ const double row_value = p.first;
+ const double col_value = p.second;
+
+ const auto xaxis_index = static_cast<size_t>(col_value);
+ const size_t yaxis_index = data.yAxis().size() - 1 - static_cast<size_t>(row_value);
+
+ const Bin1D xbin = data.xAxis().bin(xaxis_index);
+ const Bin1D ybin = data.yAxis().bin(yaxis_index);
+
+ const double dx = col_value - static_cast<size_t>(col_value);
+ const double dy = -1.0 * (row_value - static_cast<size_t>(row_value));
+
+ const double x = xbin.center() + xbin.binSize() * dx;
+ const double y = ybin.center() + ybin.binSize() * dy;
+
+ result.emplace_back(x, y);
+ }
+ return result;
+}
diff --git a/Device/Analyze/Peaks.h b/Device/Analyze/Peaks.h
new file mode 100644
index 0000000000000000000000000000000000000000..c379577693774cfc868c0f5068353b07245987d3
--- /dev/null
+++ b/Device/Analyze/Peaks.h
@@ -0,0 +1,32 @@
+// ************************************************************************************************
+//
+// BornAgain: simulate and fit reflection and scattering
+//
+//! @file Device/Analyze/Peaks.h
+//! @brief Defines function FindPeaks namespace Analyze.
+//!
+//! @homepage http://www.bornagainproject.org
+//! @license GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2018
+//! @authors Scientific Computing Group at MLZ (see CITATION, AUTHORS)
+//
+// ************************************************************************************************
+
+#ifndef BORNAGAIN_DEVICE_ANALYZE_PEAKS_H
+#define BORNAGAIN_DEVICE_ANALYZE_PEAKS_H
+
+#include <string>
+#include <vector>
+
+class Datafield;
+
+namespace Analyze {
+
+//! Returns vector of peak center coordinates, for peaks in given histogram.
+std::vector<std::pair<double, double>> FindPeaks(const Datafield& data, double sigma = 2,
+ const std::string& option = {},
+ double threshold = 0.05);
+
+} // namespace Analyze
+
+#endif // BORNAGAIN_DEVICE_ANALYZE_PEAKS_H
diff --git a/Device/Data/DataUtil.cpp b/Device/Data/DataUtil.cpp
index fc35f5ad56c98ad82671ebd68295d5c80f9e0623..43c20ea62ddf5ccaba77c5ec697e2d3e79754f6b 100644
--- a/Device/Data/DataUtil.cpp
+++ b/Device/Data/DataUtil.cpp
@@ -16,27 +16,12 @@
#include "Base/Axis/Frame.h"
#include "Base/Axis/MakeScale.h"
#include "Base/Axis/Scale.h"
-#include "Base/Math/FourierTransform.h"
#include "Base/Util/Assert.h"
#include "Device/Data/Datafield.h"
#include <algorithm>
#include <cmath>
#include <functional>
#include <stdexcept>
-#include <tspectrum.h> // third-party code, extracted from CERN ROOT (class TSpectrum2)
-
-namespace {
-
-std::vector<std::vector<double>> FT2DArray(const std::vector<std::vector<double>>& signal)
-{
- FourierTransform ft;
- std::vector<std::vector<double>> result;
- ft.fft(signal, result);
- ft.fftshift(result); // low frequency to center of array
- return result;
-}
-
-} // namespace
std::tuple<size_t, size_t, std::vector<double>>
DataUtil::flatten2D(const std::vector<std::vector<double>>& vec)
@@ -56,33 +41,6 @@ DataUtil::flatten2D(const std::vector<std::vector<double>>& vec)
return {nrows, ncols, outvec};
}
-std::unique_ptr<Datafield> DataUtil::createPField1D(const std::vector<double>& vec)
-{
- const size_t N = vec.size();
- std::vector<const Scale*> axes{newEquiDivision("axis0", N, 0.0, (double)N)};
- return std::make_unique<Datafield>(std::move(axes), vec);
-}
-
-std::unique_ptr<Datafield> DataUtil::createPField2D(const std::vector<std::vector<double>>& vec,
- const std::vector<std::vector<double>>& stdv)
-{
- const auto [nrows, ncols, outvec] = flatten2D(vec);
- ASSERT(nrows > 0);
- ASSERT(ncols > 0);
-
- std::vector<const Scale*> axes{newEquiDivision("axis0", ncols, 0.0, (double)ncols),
- newEquiDivision("axis1", nrows, 0.0, (double)nrows)};
-
- if (stdv.empty())
- return std::make_unique<Datafield>(std::move(axes), outvec);
-
- const auto [nrows2, ncols2, outvec2] = flatten2D(stdv);
- ASSERT(nrows2 == nrows);
- ASSERT(ncols2 == ncols);
-
- return std::make_unique<Datafield>(std::move(axes), outvec, outvec2);
-}
-
std::vector<std::vector<double>> DataUtil::createVector2D(const Datafield& data)
{
std::vector<std::vector<double>> result;
@@ -219,42 +177,3 @@ Datafield DataUtil::vecvecToDatafield(const std::vector<std::vector<double>>& ar
}
return {std::move(axes), out};
}
-
-Datafield DataUtil::createFFT(const Datafield& data)
-{
- auto array_2d = DataUtil::create2DArrayfromDatafield(data);
- auto fft_array_2d = FT2DArray(array_2d);
- return DataUtil::vecvecToDatafield(fft_array_2d);
-}
-
-std::vector<std::pair<double, double>> DataUtil::FindPeaks(const Datafield& data, double sigma,
- const std::string& option,
- double threshold)
-{
- const std::vector<std::vector<double>> arr = DataUtil::createVector2D(data);
- tspectrum::Spectrum2D spec;
- const auto peaks = spec.find_peaks(arr, sigma, option, threshold);
-
- // coordinates of peaks in histogram axes units
- std::vector<std::pair<double, double>> result;
-
- for (const auto& p : peaks) {
- const double row_value = p.first;
- const double col_value = p.second;
-
- const auto xaxis_index = static_cast<size_t>(col_value);
- const size_t yaxis_index = data.yAxis().size() - 1 - static_cast<size_t>(row_value);
-
- const Bin1D xbin = data.xAxis().bin(xaxis_index);
- const Bin1D ybin = data.yAxis().bin(yaxis_index);
-
- const double dx = col_value - static_cast<size_t>(col_value);
- const double dy = -1.0 * (row_value - static_cast<size_t>(row_value));
-
- const double x = xbin.center() + xbin.binSize() * dx;
- const double y = ybin.center() + ybin.binSize() * dy;
-
- result.emplace_back(x, y);
- }
- return result;
-}
diff --git a/Device/Data/DataUtil.h b/Device/Data/DataUtil.h
index 0c8524fd06c756e65d883340d82d7ca52fabb271..107b8d62f9233fbf3ecc218fef7c7a7ea60bc178 100644
--- a/Device/Data/DataUtil.h
+++ b/Device/Data/DataUtil.h
@@ -12,6 +12,9 @@
//
// ************************************************************************************************
+#ifdef SWIG
+#error no need to expose this header to Swig
+#endif // SWIG
#ifndef BORNAGAIN_DEVICE_DATA_DATAUTIL_H
#define BORNAGAIN_DEVICE_DATA_DATAUTIL_H
@@ -23,13 +26,8 @@ class Datafield;
namespace DataUtil {
-#ifndef SWIG
std::tuple<size_t, size_t, std::vector<double>> flatten2D(const std::vector<std::vector<double>>&);
-std::unique_ptr<Datafield> createPField1D(const std::vector<double>& vec);
-std::unique_ptr<Datafield> createPField2D(const std::vector<std::vector<double>>& vec,
- const std::vector<std::vector<double>>& stdv = {});
-
//! Returns new object with input data rotated by
//! n*90 deg counterclockwise (n > 0) or clockwise (n < 0)
//! Axes are swapped if the data is effectively rotated by 90 or 270 degrees
@@ -39,9 +37,6 @@ Datafield createRearrangedDataSet(const Datafield& data, int n);
//! Creates Datafield from a 2D Array.
Datafield vecvecToDatafield(const std::vector<std::vector<double>>& array_2d);
-//! Creates Fourier Transform (Datafield format) of intensity map (Datafield format).
-Datafield createFFT(const Datafield& data);
-
//! Creates 2D vector from Datafield.
std::vector<std::vector<double>> createVector2D(const Datafield& data);
@@ -54,14 +49,6 @@ std::vector<std::vector<double>> transpose(const std::vector<std::vector<double>
//! Creates a vector of vectors of double (2D Array) from Datafield.
std::vector<std::vector<double>> create2DArrayfromDatafield(const Datafield& data);
-
-#endif // SWIG
-
-//! Returns vector of peak center coordinates, for peaks in given histogram.
-std::vector<std::pair<double, double>> FindPeaks(const Datafield& data, double sigma = 2,
- const std::string& option = {},
- double threshold = 0.05);
-
} // namespace DataUtil
#endif // BORNAGAIN_DEVICE_DATA_DATAUTIL_H
diff --git a/Device/IO/ReadWrite2DTable.cpp b/Device/IO/ReadWrite2DTable.cpp
index 91abd3223675f8b645b003a2172fe60979e24788..476e8e24c0205b335943a9425db70e5c13b893e3 100644
--- a/Device/IO/ReadWrite2DTable.cpp
+++ b/Device/IO/ReadWrite2DTable.cpp
@@ -13,6 +13,7 @@
// ************************************************************************************************
#include "Device/IO/ReadWrite2DTable.h"
+#include "Base/Axis/MakeScale.h"
#include "Base/Axis/Scale.h"
#include "Base/Util/Assert.h"
#include "Base/Util/StringUtil.h"
@@ -71,37 +72,44 @@ Datafield* Util::RW::read2DTable(std::istream& input_stream)
std::vector<std::vector<double>> data;
// Read numbers from input stream:
+ size_t nrows = 0;
+ size_t ncols = 0;
while (std::getline(input_stream, line)) {
line = Base::String::trim(line);
if (line.empty() || !isDoubleStartChar(line[0]))
continue;
- std::vector<double> dataInRow = Util::Parse::parse_doubles(line);
- data.push_back(dataInRow);
- }
-
- // Validate:
- size_t nrows = data.size();
- size_t ncols(0);
- if (nrows)
- ncols = data[0].size();
- if (ncols == 0)
- throw std::runtime_error("Text file is empty or has invalid content");
- for (size_t row = 0; row < nrows; row++)
- if (data[row].size() != ncols)
+ std::vector<double> tmp = Util::Parse::parse_doubles(line);
+ if (nrows == 0)
+ ncols = tmp.size();
+ else if (tmp.size() != ncols)
throw std::runtime_error("Number of elements is not the same for all rows");
+ data.push_back(tmp);
+ ++nrows;
+ }
+ if (nrows == 0 || ncols == 0)
+ throw std::runtime_error("No data found in table");
// Convert:
- if (nrows < 2)
- return DataUtil::createPField1D(data[0]).release();
+ if (nrows < 2) {
+ std::vector<const Scale*> axes{newEquiDivision("axis0", ncols, 0.0, (double)ncols)};
+ return new Datafield(std::move(axes), data[0]);
+ }
if (ncols < 2) {
- const size_t size = data.size();
- std::vector<double> vector1d(size);
- for (size_t i = 0; i < size; ++i)
+ std::vector<const Scale*> axes{newEquiDivision("axis0", nrows, 0.0, (double)nrows)};
+ std::vector<double> vector1d(nrows);
+ for (size_t i = 0; i < nrows; ++i)
vector1d[i] = data[i][0];
- return DataUtil::createPField1D(vector1d).release();
+ return new Datafield(std::move(axes), vector1d);
}
- return DataUtil::createPField2D(data).release();
+ std::vector<const Scale*> axes{newEquiDivision("axis0", ncols, 0.0, (double)ncols),
+ newEquiDivision("axis1", nrows, 0.0, (double)nrows)};
+
+ const auto [nrows2, ncols2, outvec] = DataUtil::flatten2D(data);
+ ASSERT(nrows2 == nrows);
+ ASSERT(ncols2 == ncols);
+
+ return new Datafield(std::move(axes), outvec);
}
void Util::RW::write2DTable(const Datafield& data, std::ostream& output_stream)
diff --git a/GUI/View/Plot2D/IntensityDataFFTPresenter.cpp b/GUI/View/Plot2D/IntensityDataFFTPresenter.cpp
index 9f280cb036da1736425056f0123bde971f2fdf85..d20da4f26a062e9b825a78839acc512e6145e246 100644
--- a/GUI/View/Plot2D/IntensityDataFFTPresenter.cpp
+++ b/GUI/View/Plot2D/IntensityDataFFTPresenter.cpp
@@ -14,7 +14,7 @@
#include "GUI/View/Plot2D/IntensityDataFFTPresenter.h"
#include "Base/Util/Assert.h"
-#include "Device/Data/DataUtil.h"
+#include "Device/Analyze/Fourier.h"
#include "Device/Data/Datafield.h"
#include "GUI/Model/Data/IntensityDataItem.h"
#include <QApplication>
@@ -43,7 +43,7 @@ IntensityDataItem* IntensityDataFFTPresenter::fftItem(IntensityDataItem* origIte
ASSERT(origItem);
QApplication::setOverrideCursor(Qt::WaitCursor);
- m_fftItem->setDatafield(DataUtil::createFFT(*origItem->c_field()));
+ m_fftItem->setDatafield(Analyze::createFFT(*origItem->c_field()));
QApplication::restoreOverrideCursor();
diff --git a/Wrap/Swig/libBornAgainDevice.i b/Wrap/Swig/libBornAgainDevice.i
index 00f17c0afc020595146ff11f153ed46b3a7c5d71..2c85392fca1772db6040149af5ec50e72067796d 100644
--- a/Wrap/Swig/libBornAgainDevice.i
+++ b/Wrap/Swig/libBornAgainDevice.i
@@ -22,10 +22,10 @@
#include "Base/Axis/Frame.h"
#include "Base/Axis/Scale.h"
#include "Base/Spin/SpinMatrix.h"
+#include "Device/Analyze/Peaks.h"
#include "Device/Beam/Beam.h"
#include "Device/Beam/FootprintGauss.h"
#include "Device/Beam/FootprintSquare.h"
-#include "Device/Data/DataUtil.h"
#include "Device/Detector/OffspecDetector.h"
#include "Device/Detector/FlatDetector.h"
#include "Device/Detector/SphericalDetector.h"
@@ -52,9 +52,7 @@
%include "Base/Axis/Coordinate.h" // required by ImportSettings
-%include "Device/Data/DataUtil.h"
%include "Device/Data/Datafield.h"
-%include "Device/Data/DataUtil.h"
%include "Device/Beam/Beam.h"
%include "Device/Beam/IFootprint.h"
@@ -80,3 +78,5 @@
%include "Device/IO/DiffUtil.h"
%include "Device/IO/ImportSettings.h"
%include "Device/IO/IOFactory.h"
+
+%include "Device/Analyze/Peaks.h"
diff --git a/auto/Examples/bayesian/likelihood_sampling.py b/auto/Examples/bayesian/likelihood_sampling.py
index 49ccddbce6b207bafbd3529d7dbd86b16532e1da..1fa6401c3e32d614afc415eaac60cdbde15eb209 100755
--- a/auto/Examples/bayesian/likelihood_sampling.py
+++ b/auto/Examples/bayesian/likelihood_sampling.py
@@ -57,9 +57,12 @@ def run_simulation(points, ni_thickness, ti_thickness):
if __name__ == '__main__':
filepath = os.path.join(datadir, "specular/genx_alternating_layers.dat.gz")
- data = ba.readData2D(filepath).npArray()
- data[:, 0] *= np.pi/360 # convert incident angles from deg to rad
- data[:, 2] = data[:, 1]*0.1 # arbitrary uncertainties of 10%
+ flags = ba.ImportSettings1D("2alpha (deg)", "#", "", 1, 2)
+ data = ba.readData1D(filepath, ba.csv1D, flags)
+
+ q = data.npXcenters()
+ y = data.npArray()
+ dy = y * 0.1 # arbitrary uncertainties
def log_likelihood(P):
"""
@@ -70,9 +73,6 @@ if __name__ == '__main__':
:array yerr: the ordinate uncertainty (dR-values)
:return: log-likelihood
"""
- q = data[:, 0]
- y = data[:, 1]
- dy = data[:, 2]
y_sim = run_simulation(q, *P)
sigma2 = dy**2 + y_sim**2
return -0.5*np.sum((y - y_sim)**2/sigma2 + np.log(sigma2))
@@ -104,14 +104,10 @@ if __name__ == '__main__':
plt.show()
# Plot and show MLE and data of reflectivity
- plt.errorbar(data[:, 0],
- data[:, 1],
- data[:, 2],
- marker='.',
- ls='')
+ plt.errorbar(q, y, dy, marker='.', ls='')
plt.plot(
- data[:, 0],
- run_simulation(data[:, 0], *flat_samples.mean(axis=0)),
+ q,
+ run_simulation(q, *flat_samples.mean(axis=0)),
'-')
plt.xlabel('$\\alpha$/rad')
plt.ylabel('$R$')
diff --git a/auto/MiniExamples/bayesian/likelihood_sampling.py b/auto/MiniExamples/bayesian/likelihood_sampling.py
index 9144431229378469ed89ad12dba6bc3a0a0165be..c5094bd65147f84e277d0ac112cc5e56ae385f3d 100755
--- a/auto/MiniExamples/bayesian/likelihood_sampling.py
+++ b/auto/MiniExamples/bayesian/likelihood_sampling.py
@@ -57,9 +57,12 @@ def run_simulation(points, ni_thickness, ti_thickness):
if __name__ == '__main__':
filepath = os.path.join(datadir, "specular/genx_alternating_layers.dat.gz")
- data = ba.readData2D(filepath).npArray()
- data[:, 0] *= np.pi/360 # convert incident angles from deg to rad
- data[:, 2] = data[:, 1]*0.1 # arbitrary uncertainties of 10%
+ flags = ba.ImportSettings1D("2alpha (deg)", "#", "", 1, 2)
+ data = ba.readData1D(filepath, ba.csv1D, flags)
+
+ q = data.npXcenters()
+ y = data.npArray()
+ dy = y * 0.1 # arbitrary uncertainties
def log_likelihood(P):
"""
@@ -70,9 +73,6 @@ if __name__ == '__main__':
:array yerr: the ordinate uncertainty (dR-values)
:return: log-likelihood
"""
- q = data[:, 0]
- y = data[:, 1]
- dy = data[:, 2]
y_sim = run_simulation(q, *P)
sigma2 = dy**2 + y_sim**2
return -0.5*np.sum((y - y_sim)**2/sigma2 + np.log(sigma2))
@@ -104,14 +104,10 @@ if __name__ == '__main__':
# plt.show()
# Plot and show MLE and data of reflectivity
- plt.errorbar(data[:, 0],
- data[:, 1],
- data[:, 2],
- marker='.',
- ls='')
+ plt.errorbar(q, y, dy, marker='.', ls='')
plt.plot(
- data[:, 0],
- run_simulation(data[:, 0], *flat_samples.mean(axis=0)),
+ q,
+ run_simulation(q, *flat_samples.mean(axis=0)),
'-')
plt.xlabel('$\\alpha$/rad')
plt.ylabel('$R$')
diff --git a/auto/Wrap/libBornAgainDevice.py b/auto/Wrap/libBornAgainDevice.py
index b0322647e3bb1453892ec813e28a26bfd0665a18..c483cce9e5fb350dec836562e647ce83ea746ee2 100644
--- a/auto/Wrap/libBornAgainDevice.py
+++ b/auto/Wrap/libBornAgainDevice.py
@@ -2066,10 +2066,6 @@ class Coordinate(object):
# Register Coordinate in _libBornAgainDevice:
_libBornAgainDevice.Coordinate_swigregister(Coordinate)
-
-def FindPeaks(*args):
- r"""FindPeaks(Datafield data, double sigma=2, std::string const & option={}, double threshold=0.05) -> vector_pvacuum_double_t"""
- return _libBornAgainDevice.FindPeaks(*args)
class Datafield(object):
r"""Proxy of C++ Datafield class."""
@@ -3022,3 +3018,7 @@ def writeDatafield(data, fname):
def dataMatchesFile(data, refFileName, tol):
r"""dataMatchesFile(Datafield data, std::string const & refFileName, double tol) -> bool"""
return _libBornAgainDevice.dataMatchesFile(data, refFileName, tol)
+
+def FindPeaks(*args):
+ r"""FindPeaks(Datafield data, double sigma=2, std::string const & option={}, double threshold=0.05) -> vector_pvacuum_double_t"""
+ return _libBornAgainDevice.FindPeaks(*args)
diff --git a/auto/Wrap/libBornAgainDevice_wrap.cpp b/auto/Wrap/libBornAgainDevice_wrap.cpp
index db81f4eec3a761758921b55909bf7795042418dc..c1e4290de48e56955523474cfea6b456ce5e3b57 100644
--- a/auto/Wrap/libBornAgainDevice_wrap.cpp
+++ b/auto/Wrap/libBornAgainDevice_wrap.cpp
@@ -7017,10 +7017,10 @@ SWIGINTERN void std_vector_Sl_std_pair_Sl_double_Sc_double_Sg__Sg__insert__SWIG_
#include "Base/Axis/Frame.h"
#include "Base/Axis/Scale.h"
#include "Base/Spin/SpinMatrix.h"
+#include "Device/Analyze/Peaks.h"
#include "Device/Beam/Beam.h"
#include "Device/Beam/FootprintGauss.h"
#include "Device/Beam/FootprintSquare.h"
-#include "Device/Data/DataUtil.h"
#include "Device/Detector/OffspecDetector.h"
#include "Device/Detector/FlatDetector.h"
#include "Device/Detector/SphericalDetector.h"
@@ -28729,288 +28729,6 @@ SWIGINTERN PyObject *Coordinate_swiginit(PyObject *SWIGUNUSEDPARM(self), PyObjec
return SWIG_Python_InitShadowInstance(args);
}
-SWIGINTERN PyObject *_wrap_FindPeaks__SWIG_0(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
- PyObject *resultobj = 0;
- Datafield *arg1 = 0 ;
- double arg2 ;
- std::string *arg3 = 0 ;
- double arg4 ;
- void *argp1 = 0 ;
- int res1 = 0 ;
- double val2 ;
- int ecode2 = 0 ;
- int res3 = SWIG_OLDOBJ ;
- double val4 ;
- int ecode4 = 0 ;
- std::vector< std::pair< double,double >,std::allocator< std::pair< double,double > > > result;
-
- if ((nobjs < 4) || (nobjs > 4)) SWIG_fail;
- res1 = SWIG_ConvertPtr(swig_obj[0], &argp1, SWIGTYPE_p_Datafield, 0 | 0);
- if (!SWIG_IsOK(res1)) {
- SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "FindPeaks" "', argument " "1"" of type '" "Datafield const &""'");
- }
- if (!argp1) {
- SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "FindPeaks" "', argument " "1"" of type '" "Datafield const &""'");
- }
- arg1 = reinterpret_cast< Datafield * >(argp1);
- ecode2 = SWIG_AsVal_double(swig_obj[1], &val2);
- if (!SWIG_IsOK(ecode2)) {
- SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "FindPeaks" "', argument " "2"" of type '" "double""'");
- }
- arg2 = static_cast< double >(val2);
- {
- std::string *ptr = (std::string *)0;
- res3 = SWIG_AsPtr_std_string(swig_obj[2], &ptr);
- if (!SWIG_IsOK(res3)) {
- SWIG_exception_fail(SWIG_ArgError(res3), "in method '" "FindPeaks" "', argument " "3"" of type '" "std::string const &""'");
- }
- if (!ptr) {
- SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "FindPeaks" "', argument " "3"" of type '" "std::string const &""'");
- }
- arg3 = ptr;
- }
- ecode4 = SWIG_AsVal_double(swig_obj[3], &val4);
- if (!SWIG_IsOK(ecode4)) {
- SWIG_exception_fail(SWIG_ArgError(ecode4), "in method '" "FindPeaks" "', argument " "4"" of type '" "double""'");
- }
- arg4 = static_cast< double >(val4);
- {
- try {
- result = DataUtil::FindPeaks((Datafield const &)*arg1,arg2,(std::string const &)*arg3,arg4);
- } catch (const std::exception& ex) {
- // message shown in the Python interpreter
- const std::string msg {
- "BornAgain C++ Exception: " + std::string(ex.what())
- };
- SWIG_exception(SWIG_RuntimeError, msg.c_str());
- }
- }
- resultobj = swig::from(static_cast< std::vector< std::pair< double,double >,std::allocator< std::pair< double,double > > > >(result));
- if (SWIG_IsNewObj(res3)) delete arg3;
- return resultobj;
-fail:
- if (SWIG_IsNewObj(res3)) delete arg3;
- return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_FindPeaks__SWIG_1(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
- PyObject *resultobj = 0;
- Datafield *arg1 = 0 ;
- double arg2 ;
- std::string *arg3 = 0 ;
- void *argp1 = 0 ;
- int res1 = 0 ;
- double val2 ;
- int ecode2 = 0 ;
- int res3 = SWIG_OLDOBJ ;
- std::vector< std::pair< double,double >,std::allocator< std::pair< double,double > > > result;
-
- if ((nobjs < 3) || (nobjs > 3)) SWIG_fail;
- res1 = SWIG_ConvertPtr(swig_obj[0], &argp1, SWIGTYPE_p_Datafield, 0 | 0);
- if (!SWIG_IsOK(res1)) {
- SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "FindPeaks" "', argument " "1"" of type '" "Datafield const &""'");
- }
- if (!argp1) {
- SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "FindPeaks" "', argument " "1"" of type '" "Datafield const &""'");
- }
- arg1 = reinterpret_cast< Datafield * >(argp1);
- ecode2 = SWIG_AsVal_double(swig_obj[1], &val2);
- if (!SWIG_IsOK(ecode2)) {
- SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "FindPeaks" "', argument " "2"" of type '" "double""'");
- }
- arg2 = static_cast< double >(val2);
- {
- std::string *ptr = (std::string *)0;
- res3 = SWIG_AsPtr_std_string(swig_obj[2], &ptr);
- if (!SWIG_IsOK(res3)) {
- SWIG_exception_fail(SWIG_ArgError(res3), "in method '" "FindPeaks" "', argument " "3"" of type '" "std::string const &""'");
- }
- if (!ptr) {
- SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "FindPeaks" "', argument " "3"" of type '" "std::string const &""'");
- }
- arg3 = ptr;
- }
- {
- try {
- result = DataUtil::FindPeaks((Datafield const &)*arg1,arg2,(std::string const &)*arg3);
- } catch (const std::exception& ex) {
- // message shown in the Python interpreter
- const std::string msg {
- "BornAgain C++ Exception: " + std::string(ex.what())
- };
- SWIG_exception(SWIG_RuntimeError, msg.c_str());
- }
- }
- resultobj = swig::from(static_cast< std::vector< std::pair< double,double >,std::allocator< std::pair< double,double > > > >(result));
- if (SWIG_IsNewObj(res3)) delete arg3;
- return resultobj;
-fail:
- if (SWIG_IsNewObj(res3)) delete arg3;
- return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_FindPeaks__SWIG_2(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
- PyObject *resultobj = 0;
- Datafield *arg1 = 0 ;
- double arg2 ;
- void *argp1 = 0 ;
- int res1 = 0 ;
- double val2 ;
- int ecode2 = 0 ;
- std::vector< std::pair< double,double >,std::allocator< std::pair< double,double > > > result;
-
- if ((nobjs < 2) || (nobjs > 2)) SWIG_fail;
- res1 = SWIG_ConvertPtr(swig_obj[0], &argp1, SWIGTYPE_p_Datafield, 0 | 0);
- if (!SWIG_IsOK(res1)) {
- SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "FindPeaks" "', argument " "1"" of type '" "Datafield const &""'");
- }
- if (!argp1) {
- SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "FindPeaks" "', argument " "1"" of type '" "Datafield const &""'");
- }
- arg1 = reinterpret_cast< Datafield * >(argp1);
- ecode2 = SWIG_AsVal_double(swig_obj[1], &val2);
- if (!SWIG_IsOK(ecode2)) {
- SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "FindPeaks" "', argument " "2"" of type '" "double""'");
- }
- arg2 = static_cast< double >(val2);
- {
- try {
- result = DataUtil::FindPeaks((Datafield const &)*arg1,arg2);
- } catch (const std::exception& ex) {
- // message shown in the Python interpreter
- const std::string msg {
- "BornAgain C++ Exception: " + std::string(ex.what())
- };
- SWIG_exception(SWIG_RuntimeError, msg.c_str());
- }
- }
- resultobj = swig::from(static_cast< std::vector< std::pair< double,double >,std::allocator< std::pair< double,double > > > >(result));
- return resultobj;
-fail:
- return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_FindPeaks__SWIG_3(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
- PyObject *resultobj = 0;
- Datafield *arg1 = 0 ;
- void *argp1 = 0 ;
- int res1 = 0 ;
- std::vector< std::pair< double,double >,std::allocator< std::pair< double,double > > > result;
-
- if ((nobjs < 1) || (nobjs > 1)) SWIG_fail;
- res1 = SWIG_ConvertPtr(swig_obj[0], &argp1, SWIGTYPE_p_Datafield, 0 | 0);
- if (!SWIG_IsOK(res1)) {
- SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "FindPeaks" "', argument " "1"" of type '" "Datafield const &""'");
- }
- if (!argp1) {
- SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "FindPeaks" "', argument " "1"" of type '" "Datafield const &""'");
- }
- arg1 = reinterpret_cast< Datafield * >(argp1);
- {
- try {
- result = DataUtil::FindPeaks((Datafield const &)*arg1);
- } catch (const std::exception& ex) {
- // message shown in the Python interpreter
- const std::string msg {
- "BornAgain C++ Exception: " + std::string(ex.what())
- };
- SWIG_exception(SWIG_RuntimeError, msg.c_str());
- }
- }
- resultobj = swig::from(static_cast< std::vector< std::pair< double,double >,std::allocator< std::pair< double,double > > > >(result));
- return resultobj;
-fail:
- return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_FindPeaks(PyObject *self, PyObject *args) {
- Py_ssize_t argc;
- PyObject *argv[5] = {
- 0
- };
-
- if (!(argc = SWIG_Python_UnpackTuple(args, "FindPeaks", 0, 4, argv))) SWIG_fail;
- --argc;
- if (argc == 1) {
- int _v = 0;
- int res = SWIG_ConvertPtr(argv[0], 0, SWIGTYPE_p_Datafield, SWIG_POINTER_NO_NULL | 0);
- _v = SWIG_CheckState(res);
- if (_v) {
- return _wrap_FindPeaks__SWIG_3(self, argc, argv);
- }
- }
- if (argc == 2) {
- int _v = 0;
- int res = SWIG_ConvertPtr(argv[0], 0, SWIGTYPE_p_Datafield, SWIG_POINTER_NO_NULL | 0);
- _v = SWIG_CheckState(res);
- if (_v) {
- {
- int res = SWIG_AsVal_double(argv[1], NULL);
- _v = SWIG_CheckState(res);
- }
- if (_v) {
- return _wrap_FindPeaks__SWIG_2(self, argc, argv);
- }
- }
- }
- if (argc == 3) {
- int _v = 0;
- int res = SWIG_ConvertPtr(argv[0], 0, SWIGTYPE_p_Datafield, SWIG_POINTER_NO_NULL | 0);
- _v = SWIG_CheckState(res);
- if (_v) {
- {
- int res = SWIG_AsVal_double(argv[1], NULL);
- _v = SWIG_CheckState(res);
- }
- if (_v) {
- int res = SWIG_AsPtr_std_string(argv[2], (std::string**)(0));
- _v = SWIG_CheckState(res);
- if (_v) {
- return _wrap_FindPeaks__SWIG_1(self, argc, argv);
- }
- }
- }
- }
- if (argc == 4) {
- int _v = 0;
- int res = SWIG_ConvertPtr(argv[0], 0, SWIGTYPE_p_Datafield, SWIG_POINTER_NO_NULL | 0);
- _v = SWIG_CheckState(res);
- if (_v) {
- {
- int res = SWIG_AsVal_double(argv[1], NULL);
- _v = SWIG_CheckState(res);
- }
- if (_v) {
- int res = SWIG_AsPtr_std_string(argv[2], (std::string**)(0));
- _v = SWIG_CheckState(res);
- if (_v) {
- {
- int res = SWIG_AsVal_double(argv[3], NULL);
- _v = SWIG_CheckState(res);
- }
- if (_v) {
- return _wrap_FindPeaks__SWIG_0(self, argc, argv);
- }
- }
- }
- }
- }
-
-fail:
- SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'FindPeaks'.\n"
- " Possible C/C++ prototypes are:\n"
- " DataUtil::FindPeaks(Datafield const &,double,std::string const &,double)\n"
- " DataUtil::FindPeaks(Datafield const &,double,std::string const &)\n"
- " DataUtil::FindPeaks(Datafield const &,double)\n"
- " DataUtil::FindPeaks(Datafield const &)\n");
- return 0;
-}
-
-
SWIGINTERN PyObject *_wrap_new_Datafield__SWIG_0(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
PyObject *resultobj = 0;
std::string arg1 ;
@@ -43315,6 +43033,288 @@ fail:
}
+SWIGINTERN PyObject *_wrap_FindPeaks__SWIG_0(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
+ PyObject *resultobj = 0;
+ Datafield *arg1 = 0 ;
+ double arg2 ;
+ std::string *arg3 = 0 ;
+ double arg4 ;
+ void *argp1 = 0 ;
+ int res1 = 0 ;
+ double val2 ;
+ int ecode2 = 0 ;
+ int res3 = SWIG_OLDOBJ ;
+ double val4 ;
+ int ecode4 = 0 ;
+ std::vector< std::pair< double,double >,std::allocator< std::pair< double,double > > > result;
+
+ if ((nobjs < 4) || (nobjs > 4)) SWIG_fail;
+ res1 = SWIG_ConvertPtr(swig_obj[0], &argp1, SWIGTYPE_p_Datafield, 0 | 0);
+ if (!SWIG_IsOK(res1)) {
+ SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "FindPeaks" "', argument " "1"" of type '" "Datafield const &""'");
+ }
+ if (!argp1) {
+ SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "FindPeaks" "', argument " "1"" of type '" "Datafield const &""'");
+ }
+ arg1 = reinterpret_cast< Datafield * >(argp1);
+ ecode2 = SWIG_AsVal_double(swig_obj[1], &val2);
+ if (!SWIG_IsOK(ecode2)) {
+ SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "FindPeaks" "', argument " "2"" of type '" "double""'");
+ }
+ arg2 = static_cast< double >(val2);
+ {
+ std::string *ptr = (std::string *)0;
+ res3 = SWIG_AsPtr_std_string(swig_obj[2], &ptr);
+ if (!SWIG_IsOK(res3)) {
+ SWIG_exception_fail(SWIG_ArgError(res3), "in method '" "FindPeaks" "', argument " "3"" of type '" "std::string const &""'");
+ }
+ if (!ptr) {
+ SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "FindPeaks" "', argument " "3"" of type '" "std::string const &""'");
+ }
+ arg3 = ptr;
+ }
+ ecode4 = SWIG_AsVal_double(swig_obj[3], &val4);
+ if (!SWIG_IsOK(ecode4)) {
+ SWIG_exception_fail(SWIG_ArgError(ecode4), "in method '" "FindPeaks" "', argument " "4"" of type '" "double""'");
+ }
+ arg4 = static_cast< double >(val4);
+ {
+ try {
+ result = Analyze::FindPeaks((Datafield const &)*arg1,arg2,(std::string const &)*arg3,arg4);
+ } catch (const std::exception& ex) {
+ // message shown in the Python interpreter
+ const std::string msg {
+ "BornAgain C++ Exception: " + std::string(ex.what())
+ };
+ SWIG_exception(SWIG_RuntimeError, msg.c_str());
+ }
+ }
+ resultobj = swig::from(static_cast< std::vector< std::pair< double,double >,std::allocator< std::pair< double,double > > > >(result));
+ if (SWIG_IsNewObj(res3)) delete arg3;
+ return resultobj;
+fail:
+ if (SWIG_IsNewObj(res3)) delete arg3;
+ return NULL;
+}
+
+
+SWIGINTERN PyObject *_wrap_FindPeaks__SWIG_1(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
+ PyObject *resultobj = 0;
+ Datafield *arg1 = 0 ;
+ double arg2 ;
+ std::string *arg3 = 0 ;
+ void *argp1 = 0 ;
+ int res1 = 0 ;
+ double val2 ;
+ int ecode2 = 0 ;
+ int res3 = SWIG_OLDOBJ ;
+ std::vector< std::pair< double,double >,std::allocator< std::pair< double,double > > > result;
+
+ if ((nobjs < 3) || (nobjs > 3)) SWIG_fail;
+ res1 = SWIG_ConvertPtr(swig_obj[0], &argp1, SWIGTYPE_p_Datafield, 0 | 0);
+ if (!SWIG_IsOK(res1)) {
+ SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "FindPeaks" "', argument " "1"" of type '" "Datafield const &""'");
+ }
+ if (!argp1) {
+ SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "FindPeaks" "', argument " "1"" of type '" "Datafield const &""'");
+ }
+ arg1 = reinterpret_cast< Datafield * >(argp1);
+ ecode2 = SWIG_AsVal_double(swig_obj[1], &val2);
+ if (!SWIG_IsOK(ecode2)) {
+ SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "FindPeaks" "', argument " "2"" of type '" "double""'");
+ }
+ arg2 = static_cast< double >(val2);
+ {
+ std::string *ptr = (std::string *)0;
+ res3 = SWIG_AsPtr_std_string(swig_obj[2], &ptr);
+ if (!SWIG_IsOK(res3)) {
+ SWIG_exception_fail(SWIG_ArgError(res3), "in method '" "FindPeaks" "', argument " "3"" of type '" "std::string const &""'");
+ }
+ if (!ptr) {
+ SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "FindPeaks" "', argument " "3"" of type '" "std::string const &""'");
+ }
+ arg3 = ptr;
+ }
+ {
+ try {
+ result = Analyze::FindPeaks((Datafield const &)*arg1,arg2,(std::string const &)*arg3);
+ } catch (const std::exception& ex) {
+ // message shown in the Python interpreter
+ const std::string msg {
+ "BornAgain C++ Exception: " + std::string(ex.what())
+ };
+ SWIG_exception(SWIG_RuntimeError, msg.c_str());
+ }
+ }
+ resultobj = swig::from(static_cast< std::vector< std::pair< double,double >,std::allocator< std::pair< double,double > > > >(result));
+ if (SWIG_IsNewObj(res3)) delete arg3;
+ return resultobj;
+fail:
+ if (SWIG_IsNewObj(res3)) delete arg3;
+ return NULL;
+}
+
+
+SWIGINTERN PyObject *_wrap_FindPeaks__SWIG_2(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
+ PyObject *resultobj = 0;
+ Datafield *arg1 = 0 ;
+ double arg2 ;
+ void *argp1 = 0 ;
+ int res1 = 0 ;
+ double val2 ;
+ int ecode2 = 0 ;
+ std::vector< std::pair< double,double >,std::allocator< std::pair< double,double > > > result;
+
+ if ((nobjs < 2) || (nobjs > 2)) SWIG_fail;
+ res1 = SWIG_ConvertPtr(swig_obj[0], &argp1, SWIGTYPE_p_Datafield, 0 | 0);
+ if (!SWIG_IsOK(res1)) {
+ SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "FindPeaks" "', argument " "1"" of type '" "Datafield const &""'");
+ }
+ if (!argp1) {
+ SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "FindPeaks" "', argument " "1"" of type '" "Datafield const &""'");
+ }
+ arg1 = reinterpret_cast< Datafield * >(argp1);
+ ecode2 = SWIG_AsVal_double(swig_obj[1], &val2);
+ if (!SWIG_IsOK(ecode2)) {
+ SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "FindPeaks" "', argument " "2"" of type '" "double""'");
+ }
+ arg2 = static_cast< double >(val2);
+ {
+ try {
+ result = Analyze::FindPeaks((Datafield const &)*arg1,arg2);
+ } catch (const std::exception& ex) {
+ // message shown in the Python interpreter
+ const std::string msg {
+ "BornAgain C++ Exception: " + std::string(ex.what())
+ };
+ SWIG_exception(SWIG_RuntimeError, msg.c_str());
+ }
+ }
+ resultobj = swig::from(static_cast< std::vector< std::pair< double,double >,std::allocator< std::pair< double,double > > > >(result));
+ return resultobj;
+fail:
+ return NULL;
+}
+
+
+SWIGINTERN PyObject *_wrap_FindPeaks__SWIG_3(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
+ PyObject *resultobj = 0;
+ Datafield *arg1 = 0 ;
+ void *argp1 = 0 ;
+ int res1 = 0 ;
+ std::vector< std::pair< double,double >,std::allocator< std::pair< double,double > > > result;
+
+ if ((nobjs < 1) || (nobjs > 1)) SWIG_fail;
+ res1 = SWIG_ConvertPtr(swig_obj[0], &argp1, SWIGTYPE_p_Datafield, 0 | 0);
+ if (!SWIG_IsOK(res1)) {
+ SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "FindPeaks" "', argument " "1"" of type '" "Datafield const &""'");
+ }
+ if (!argp1) {
+ SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "FindPeaks" "', argument " "1"" of type '" "Datafield const &""'");
+ }
+ arg1 = reinterpret_cast< Datafield * >(argp1);
+ {
+ try {
+ result = Analyze::FindPeaks((Datafield const &)*arg1);
+ } catch (const std::exception& ex) {
+ // message shown in the Python interpreter
+ const std::string msg {
+ "BornAgain C++ Exception: " + std::string(ex.what())
+ };
+ SWIG_exception(SWIG_RuntimeError, msg.c_str());
+ }
+ }
+ resultobj = swig::from(static_cast< std::vector< std::pair< double,double >,std::allocator< std::pair< double,double > > > >(result));
+ return resultobj;
+fail:
+ return NULL;
+}
+
+
+SWIGINTERN PyObject *_wrap_FindPeaks(PyObject *self, PyObject *args) {
+ Py_ssize_t argc;
+ PyObject *argv[5] = {
+ 0
+ };
+
+ if (!(argc = SWIG_Python_UnpackTuple(args, "FindPeaks", 0, 4, argv))) SWIG_fail;
+ --argc;
+ if (argc == 1) {
+ int _v = 0;
+ int res = SWIG_ConvertPtr(argv[0], 0, SWIGTYPE_p_Datafield, SWIG_POINTER_NO_NULL | 0);
+ _v = SWIG_CheckState(res);
+ if (_v) {
+ return _wrap_FindPeaks__SWIG_3(self, argc, argv);
+ }
+ }
+ if (argc == 2) {
+ int _v = 0;
+ int res = SWIG_ConvertPtr(argv[0], 0, SWIGTYPE_p_Datafield, SWIG_POINTER_NO_NULL | 0);
+ _v = SWIG_CheckState(res);
+ if (_v) {
+ {
+ int res = SWIG_AsVal_double(argv[1], NULL);
+ _v = SWIG_CheckState(res);
+ }
+ if (_v) {
+ return _wrap_FindPeaks__SWIG_2(self, argc, argv);
+ }
+ }
+ }
+ if (argc == 3) {
+ int _v = 0;
+ int res = SWIG_ConvertPtr(argv[0], 0, SWIGTYPE_p_Datafield, SWIG_POINTER_NO_NULL | 0);
+ _v = SWIG_CheckState(res);
+ if (_v) {
+ {
+ int res = SWIG_AsVal_double(argv[1], NULL);
+ _v = SWIG_CheckState(res);
+ }
+ if (_v) {
+ int res = SWIG_AsPtr_std_string(argv[2], (std::string**)(0));
+ _v = SWIG_CheckState(res);
+ if (_v) {
+ return _wrap_FindPeaks__SWIG_1(self, argc, argv);
+ }
+ }
+ }
+ }
+ if (argc == 4) {
+ int _v = 0;
+ int res = SWIG_ConvertPtr(argv[0], 0, SWIGTYPE_p_Datafield, SWIG_POINTER_NO_NULL | 0);
+ _v = SWIG_CheckState(res);
+ if (_v) {
+ {
+ int res = SWIG_AsVal_double(argv[1], NULL);
+ _v = SWIG_CheckState(res);
+ }
+ if (_v) {
+ int res = SWIG_AsPtr_std_string(argv[2], (std::string**)(0));
+ _v = SWIG_CheckState(res);
+ if (_v) {
+ {
+ int res = SWIG_AsVal_double(argv[3], NULL);
+ _v = SWIG_CheckState(res);
+ }
+ if (_v) {
+ return _wrap_FindPeaks__SWIG_0(self, argc, argv);
+ }
+ }
+ }
+ }
+ }
+
+fail:
+ SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'FindPeaks'.\n"
+ " Possible C/C++ prototypes are:\n"
+ " Analyze::FindPeaks(Datafield const &,double,std::string const &,double)\n"
+ " Analyze::FindPeaks(Datafield const &,double,std::string const &)\n"
+ " Analyze::FindPeaks(Datafield const &,double)\n"
+ " Analyze::FindPeaks(Datafield const &)\n");
+ return 0;
+}
+
+
static PyMethodDef SwigMethods[] = {
{ "delete_SwigPyIterator", _wrap_delete_SwigPyIterator, METH_O, "delete_SwigPyIterator(SwigPyIterator self)"},
{ "SwigPyIterator_value", _wrap_SwigPyIterator_value, METH_O, "SwigPyIterator_value(SwigPyIterator self) -> PyObject *"},
@@ -44019,7 +44019,6 @@ static PyMethodDef SwigMethods[] = {
{ "delete_Coordinate", _wrap_delete_Coordinate, METH_O, "delete_Coordinate(Coordinate self)"},
{ "Coordinate_swigregister", Coordinate_swigregister, METH_O, NULL},
{ "Coordinate_swiginit", Coordinate_swiginit, METH_VARARGS, NULL},
- { "FindPeaks", _wrap_FindPeaks, METH_VARARGS, "FindPeaks(Datafield data, double sigma=2, std::string const & option={}, double threshold=0.05) -> vector_pvacuum_double_t"},
{ "new_Datafield", _wrap_new_Datafield, METH_VARARGS, "\n"
"Datafield(std::string title, Frame frame, vdouble1d_t values, vdouble1d_t errSigmas={})\n"
"Datafield(std::string title, Frame frame)\n"
@@ -44342,6 +44341,7 @@ static PyMethodDef SwigMethods[] = {
""},
{ "writeDatafield", _wrap_writeDatafield, METH_VARARGS, "writeDatafield(Datafield data, std::string const & fname)"},
{ "dataMatchesFile", _wrap_dataMatchesFile, METH_VARARGS, "dataMatchesFile(Datafield data, std::string const & refFileName, double tol) -> bool"},
+ { "FindPeaks", _wrap_FindPeaks, METH_VARARGS, "FindPeaks(Datafield data, double sigma=2, std::string const & option={}, double threshold=0.05) -> vector_pvacuum_double_t"},
{ NULL, NULL, 0, NULL }
};
diff --git a/rawEx/bayesian/likelihood_sampling.py b/rawEx/bayesian/likelihood_sampling.py
index ce7f70578fd2690e35a7b13fc586ef6ab8af96c9..be2f4332cfe7ecf3f6dd552706f20c3a571be465 100755
--- a/rawEx/bayesian/likelihood_sampling.py
+++ b/rawEx/bayesian/likelihood_sampling.py
@@ -57,9 +57,12 @@ def run_simulation(points, ni_thickness, ti_thickness):
if __name__ == '__main__':
filepath = os.path.join(datadir, "specular/genx_alternating_layers.dat.gz")
- data = ba.readData2D(filepath).npArray()
- data[:, 0] *= np.pi/360 # convert incident angles from deg to rad
- data[:, 2] = data[:, 1]*0.1 # arbitrary uncertainties of 10%
+ flags = ba.ImportSettings1D("2alpha (deg)", "#", "", 1, 2)
+ data = ba.readData1D(filepath, ba.csv1D, flags)
+
+ q = data.npXcenters()
+ y = data.npArray()
+ dy = y * 0.1 # arbitrary uncertainties
def log_likelihood(P):
"""
@@ -70,9 +73,6 @@ if __name__ == '__main__':
:array yerr: the ordinate uncertainty (dR-values)
:return: log-likelihood
"""
- q = data[:, 0]
- y = data[:, 1]
- dy = data[:, 2]
y_sim = run_simulation(q, *P)
sigma2 = dy**2 + y_sim**2
return -0.5*np.sum((y - y_sim)**2/sigma2 + np.log(sigma2))
@@ -104,14 +104,10 @@ if __name__ == '__main__':
<%= sm ? "# " : "" %>plt.show()
# Plot and show MLE and data of reflectivity
- plt.errorbar(data[:, 0],
- data[:, 1],
- data[:, 2],
- marker='.',
- ls='')
+ plt.errorbar(q, y, dy, marker='.', ls='')
plt.plot(
- data[:, 0],
- run_simulation(data[:, 0], *flat_samples.mean(axis=0)),
+ q,
+ run_simulation(q, *flat_samples.mean(axis=0)),
'-')
plt.xlabel('$\\alpha$/rad')
plt.ylabel('$R$')