diff --git a/Core/Instrument/IntensityDataFunctions.cpp b/Core/Instrument/IntensityDataFunctions.cpp
index 338f267f81b70ce2047f491e4172dc160be052a5..4f354ca75ef08d0236c4285a7a1bfa536d3f8591 100644
--- a/Core/Instrument/IntensityDataFunctions.cpp
+++ b/Core/Instrument/IntensityDataFunctions.cpp
@@ -18,8 +18,26 @@
#include "FourierTransform.h"
#include "IHistogram.h"
#include "Numeric.h"
+#include "SimulationResult.h"
#include <math.h>
+//! Returns sum of relative differences between each pair of elements:
+//! (a, b) -> 2*abs(a - b)/(|a| + |b|) ( and zero if a=b=0 within epsilon )
+double IntensityDataFunctions::RelativeDifference(const SimulationResult& dat,
+ const SimulationResult& ref)
+{
+ if (dat.size() != ref.size())
+ throw std::runtime_error("Error in IntensityDataFunctions::RelativeDifference: "
+ "different number of elements");
+ if (dat.size()==0) return 0.0;
+ double sum_of_diff = 0.0;
+ for (size_t i=0; i<dat.size(); ++i) {
+ sum_of_diff += Numeric::get_relative_difference(dat[i], ref[i]);
+ }
+ return sum_of_diff / dat.size();
+}
+
+
//! Returns relative difference between two data sets sum(dat[i] - ref[i])/ref[i]).
double IntensityDataFunctions::getRelativeDifference(const OutputData<double>& dat,
const OutputData<double>& ref)
diff --git a/Core/Instrument/SimulationResult.cpp b/Core/Instrument/SimulationResult.cpp
index c66490aecfec75f7217681cd47162eac415ccf11..faa83f6a8778a6601a5df5f4f06cb9f4018e0a28 100644
--- a/Core/Instrument/SimulationResult.cpp
+++ b/Core/Instrument/SimulationResult.cpp
@@ -69,6 +69,27 @@ Histogram2D* SimulationResult::histogram2d(AxesUnits units_type) const
return new Histogram2D(*P_data);
}
+double& SimulationResult::operator[](size_t i)
+{
+ if (mP_data) return (*mP_data)[i];
+ throw std::runtime_error("Error in SimulationResult::operator[]: "
+ "no data initialized");
+}
+
+const double& SimulationResult::operator[](size_t i) const
+{
+ if (mP_data) return (*mP_data)[i];
+ throw std::runtime_error("Error in SimulationResult::operator[]: "
+ "no data initialized");
+}
+
+size_t SimulationResult::size() const
+{
+ if (mP_data) return mP_data->getAllocatedSize();
+ throw std::runtime_error("Error in SimulationResult::size(): "
+ "no data initialized");
+}
+
PyObject* SimulationResult::array() const
{
return mP_data->getArray();
diff --git a/Core/Instrument/SimulationResult.h b/Core/Instrument/SimulationResult.h
index 897c28a253f982edcf8e89555639e57a0677316c..4b6fbc2d04b6dcffeb019ce20ee29d53e16789a8 100644
--- a/Core/Instrument/SimulationResult.h
+++ b/Core/Instrument/SimulationResult.h
@@ -48,6 +48,11 @@ public:
OutputData<double>* data(AxesUnits units_type = AxesUnits::DEFAULT) const;
Histogram2D* histogram2d(AxesUnits units_type = AxesUnits::DEFAULT) const;
+ //! Data element access
+ double& operator[](size_t i);
+ const double& operator[](size_t i) const;
+ size_t size() const;
+
//! returns data as Python numpy array
#ifdef BORNAGAIN_PYTHON
PyObject* array() const;
diff --git a/Core/Tools/Numeric.cpp b/Core/Tools/Numeric.cpp
index dd9ffbb70e847724bf117443f6bab0e78f8705b8..2f8ff7c886ad807607fcf7e0e02bbc4f8e9773b1 100644
--- a/Core/Tools/Numeric.cpp
+++ b/Core/Tools/Numeric.cpp
@@ -28,17 +28,15 @@ bool areAlmostEqual(double a, double b, double tolerance)
return std::abs(a-b) <= eps * std::max( tolerance*eps, std::max(1., tolerance)*std::abs(b) );
}
-//! Returns the safe relative difference, which is |(a-b)/b| except in special cases.
+//! Returns the safe relative difference, which is 2(|a-b|)/(|a|+|b|) except in special cases.
double get_relative_difference(double a, double b)
{
constexpr double eps = std::numeric_limits<double>::epsilon();
+ const double avg_abs = (std::abs(a) + std::abs(b))/2.0;
// return 0.0 if relative error smaller than epsilon
- if (std::abs(a-b) <= eps*std::abs(b))
+ if (std::abs(a-b) <= eps*avg_abs)
return 0.0;
- // for small numbers, divide by epsilon (to avoid catastrophic cancellation)
- if (std::abs(b) <= eps)
- return std::abs((a-b)/eps);
- return std::abs((a-b)/b);
+ return std::abs(a-b)/avg_abs;
}
} // Numeric namespace
diff --git a/Tests/Functional/Python/PyCore/sliced_composition.py b/Tests/Functional/Python/PyCore/sliced_composition.py
index e5779f6369467e442f2fb908f36ab375e759cbf4..265274227a5743517488ae57e982822549cae19e 100644
--- a/Tests/Functional/Python/PyCore/sliced_composition.py
+++ b/Tests/Functional/Python/PyCore/sliced_composition.py
@@ -42,11 +42,11 @@ class SlicedSpheresTest(unittest.TestCase):
multi_layer.addLayer(substrate)
return multi_layer
- def get_intensity_data(self, particle_to_air=None, particle_to_substrate=None):
+ def get_result(self, particle_to_air=None, particle_to_substrate=None):
sample = self.get_sample(particle_to_air, particle_to_substrate)
simulation = utils.get_simulation_MiniGISAS(sample)
simulation.runSimulation()
- return simulation.getIntensityData()
+ return simulation.result()
def get_composition(self, top_material, bottom_material):
"""
@@ -91,13 +91,13 @@ class SlicedSpheresTest(unittest.TestCase):
# spherical particle
sphere = ba.Particle(mParticle, ba.FormFactorFullSphere(sphere_radius))
- reference = self.get_intensity_data(sphere)
+ reference = self.get_result(sphere)
# spherical composition
composition = self.get_composition(mParticle, mParticle)
- data = self.get_intensity_data(composition)
+ data = self.get_result(composition)
- diff = ba.getRelativeDifference(data, reference)
+ diff = ba.RelativeDifference(data, reference)
print(diff)
self.assertLess(diff, 1e-10)
@@ -109,13 +109,13 @@ class SlicedSpheresTest(unittest.TestCase):
"""
composition = self.get_composition(mParticle, mSubstrate)
- reference = self.get_intensity_data(composition)
+ reference = self.get_result(composition)
# spherical composition
composition2 = self.get_composition_v2(mParticle, mSubstrate)
- data = self.get_intensity_data(composition2)
+ data = self.get_result(composition2)
- diff = ba.getRelativeDifference(data, reference)
+ diff = ba.RelativeDifference(data, reference)
print(diff)
self.assertLess(diff, 1e-10)
@@ -131,14 +131,14 @@ class SlicedSpheresTest(unittest.TestCase):
shift = bottom_cup_height
# Scattering from empty multilayer
- reference = self.get_intensity_data()
+ reference = self.get_result()
# spherical composition
composition = self.get_composition(mAmbience, mSubstrate)
composition.setPosition(0, 0, -shift)
- data = self.get_intensity_data(composition)
+ data = self.get_result(composition)
- diff = ba.getRelativeDifference(data, reference)
+ diff = ba.RelativeDifference(data, reference)
print(diff)
self.assertLess(diff, 1e-10)
@@ -155,14 +155,14 @@ class SlicedSpheresTest(unittest.TestCase):
# spherical particle
sphere = ba.Particle(mParticle, ba.FormFactorFullSphere(sphere_radius))
sphere.setPosition(0, 0, -shift)
- reference = self.get_intensity_data(sphere)
+ reference = self.get_result(sphere)
# spherical composition
composition = self.get_composition(mParticle, mParticle)
composition.setPosition(0, 0, -shift)
- data = self.get_intensity_data(composition)
+ data = self.get_result(composition)
- diff = ba.getRelativeDifference(data, reference)
+ diff = ba.RelativeDifference(data, reference)
print(diff)
self.assertLess(diff, 1e-10)
@@ -177,33 +177,19 @@ class SlicedSpheresTest(unittest.TestCase):
# truncated sphere on top of substrate with height 16nm
truncatedSphere = ba.Particle(mParticle, ba.FormFactorTruncatedSphere(sphere_radius, sphere_radius*2 - bottom_cup_height))
- reference = self.get_intensity_data(truncatedSphere)
+ reference = self.get_result(truncatedSphere)
# Particle composition, top part made of same material, as particle. Bottom part made of same material as substrate.
composition = self.get_composition(mParticle, mSubstrate)
composition_shift = bottom_cup_height
composition.setPosition(0, 0, -composition_shift)
- data = self.get_intensity_data(composition)
+ data = self.get_result(composition)
- diff = ba.getRelativeDifference(data, reference)
+ diff = ba.RelativeDifference(data, reference)
print(diff)
self.assertLess(diff, 1e-10)
- # def testExample(self):
- # """
- # Temporary method to generate reference data
- # """
- # m_top_cup = ba.HomogeneousMaterial("Ag", 1.245e-5, 5.419e-7)
- # m_bottom_cup = ba.HomogeneousMaterial("Teflon", 2.900e-6, 6.019e-9)
- #
- # composition = self.get_composition(m_top_cup, m_bottom_cup)
- # shift = 4*nm
- # composition.setPosition(0, 0, -shift)
- # data = self.get_intensity_data(composition)
- # ba.IntensityDataIOFactory.writeIntensityData(data, "SlicedComposition.int")
-
-
if __name__ == '__main__':
unittest.main()