Provide data accessor and RelativeDifference for SimulationResult (2)

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)

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();

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;

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

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()