Read/Write to plain 2D tables without swapping axes or rotating contents.

Device/Data/DataUtil.cpp

@@ -23,24 +23,6 @@
 #include <functional>
 #include <stdexcept>
 
-std::tuple<size_t, size_t, std::vector<double>>
-DataUtil::flatten2D(const std::vector<std::vector<double>>& vec)
-{
-    size_t nrows = vec.size();
-    size_t ncols(0);
-    if (nrows)
-        ncols = vec[0].size();
-    std::vector<double> outvec(nrows * ncols);
-    for (size_t row = nrows - 1; row != (size_t)-1; --row) {
-        ASSERT(vec[row].size() == ncols);
-        for (size_t col = 0; col < ncols; ++col)
-            // TODO: why do we need to shuffle cols and rows and reverse row order?
-            outvec[nrows - row - 1 + col * nrows] = vec[row][col];
-    }
-    ASSERT(outvec.size() == nrows * ncols);
-    return {nrows, ncols, outvec};
-}
-
 std::vector<std::vector<double>> DataUtil::createVector2D(const Datafield& data)
 {
     std::vector<std::vector<double>> result;
@@ -52,10 +34,8 @@ std::vector<std::vector<double>> DataUtil::createVector2D(const Datafield& data)
 
     for (size_t row = 0; row < nrows; ++row) {
         result[row].resize(ncols, 0.0);
-        for (size_t col = 0; col < ncols; ++col) {
-            size_t globalbin = nrows - row - 1 + col * nrows;
-            result[row][col] = data[globalbin];
-        }
+        for (size_t col = 0; col < ncols; ++col)
+            result[row][col] = data[row * ncols + col];
     }
 
     return result;

Device/Data/DataUtil.h

@@ -26,8 +26,6 @@ class Datafield;
 
 namespace DataUtil {
 
-std::tuple<size_t, size_t, std::vector<double>> flatten2D(const std::vector<std::vector<double>>&);
-
 //! 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

Device/Data/Datafield.cpp

@@ -43,31 +43,26 @@ PyObject* npExport(const Frame& frame, const std::vector<double>& flatData)
     std::vector<size_t> dimensions;
     for (size_t i = 0; i < frame.rank(); i++)
         dimensions.push_back(frame.axis(i).size());
-    // for rot90 of 2-dim arrays to conform with numpy
-    if (dimensions.size() == 2)
-        std::swap(dimensions[0], dimensions[1]);
+    std::reverse(dimensions.begin(), dimensions.end());
 
     // creating ndarray objects describing size of dimensions
     PyObjectPtr pyarray{PyInterpreter::Numpy::arrayND(dimensions)};
     ASSERT(pyarray.valid());
 
     // get the pointer to the data buffer of the array (assumed to be C-contiguous)
-    double* data{PyInterpreter::Numpy::getDataPtr(pyarray.get())};
-    ASSERT(data);
-
-    double* array_buffer = data;
+    double* array_buffer{PyInterpreter::Numpy::getDataPtr(pyarray.get())};
+    ASSERT(array_buffer);
 
     // filling numpy array with output_data
-    if (frame.rank() == 2) {
+    if (frame.rank() == 1) {
+        for (size_t i = 0; i < frame.size(); ++i)
+            *array_buffer++ = flatData[i];
+    } else if (frame.rank() == 2) {
         for (size_t i = 0; i < frame.size(); ++i) {
             std::vector<int> axes_indices = frame.allIndices(i);
-            size_t offset = axes_indices[0]
-                            + frame.axis(0).size() * (frame.axis(1).size() - 1 - axes_indices[1]);
+            size_t offset = axes_indices[0] * frame.axis(1).size() + axes_indices[1];
             array_buffer[offset] = flatData[i];
         }
-    } else if (frame.rank() == 1) {
-        for (size_t i = 0; i < frame.size(); ++i)
-            *array_buffer++ = flatData[i];
     } else
         ASSERT_NEVER;
 

Device/IO/ReadWrite2DTable.cpp

@@ -69,7 +69,7 @@ void write2DRepresentation(const Datafield& data, std::ostream& output_stream)
 Datafield* Util::RW::read2DTable(std::istream& input_stream)
 {
     std::string line;
-    std::vector<std::vector<double>> data;
+    std::vector<double> values;
 
     // Read numbers from input stream:
     size_t nrows = 0;
@@ -83,33 +83,24 @@ Datafield* Util::RW::read2DTable(std::istream& input_stream)
             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);
+        for (const double val : tmp)
+            values.push_back(val);
         ++nrows;
     }
     if (nrows == 0 || ncols == 0)
         throw std::runtime_error("No data found in table");
 
-    // Convert:
-    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) {
-        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 new Datafield(std::move(axes), vector1d);
-    }
-
-    std::vector<const Scale*> axes{newEquiDivision("axis0", ncols, 0.0, (double)ncols),
-                                   newEquiDivision("axis1", nrows, 0.0, (double)nrows)};
+    std::vector<const Scale*> axes;
 
-    const auto [nrows2, ncols2, outvec] = DataUtil::flatten2D(data);
-    ASSERT(nrows2 == nrows);
-    ASSERT(ncols2 == ncols);
+    if (nrows < 2)
+        axes = {newEquiDivision("axis0", ncols, 0.0, (double)ncols)};
+    else if (ncols < 2)
+        axes = {newEquiDivision("axis0", nrows, 0.0, (double)nrows)};
+    else
+        axes = {newEquiDivision("axis0", ncols, 0.0, (double)ncols),
+                newEquiDivision("axis1", nrows, 0.0, (double)nrows)};
 
-    return new Datafield(std::move(axes), outvec);
+    return new Datafield(std::move(axes), values);
 }
 
 void Util::RW::write2DTable(const Datafield& data, std::ostream& output_stream)

Tests/Py/Unit/intensitydata_io.py

 # Functional test: tests of IO operations with the IntensityData object
 
-import math, unittest, numpy
+import unittest, numpy
 import bornagain as ba
-from bornagain import deg
-
-
-def is_the_same_data(data1, data2):
-    """
-    Checks if two data are identical
-    """
-    if data1.size() != data2.size():
-        return False
-    if data1.rank() != data2.rank():
-        return False
-    for i in range(0, data1.rank()):
-        if data1.axis(i) != data2.axis(i):
-            return False
-    for i in range(0, data1.size()):
-        if data1[i] != data2[i]:
-            return False
-    return True
-
-
-def get_boundaries_flat_in_sin(nbins, start, end):
-    """
-    Returns flat_in_sin binning of angle axis
-    """
-    result = []
-    start_sin = math.sin(start)
-    end_sin = math.sin(end)
-    step = (end_sin - start_sin) / nbins
-    for i in range(0, nbins + 1):
-        result.append(math.degrees(math.asin(start_sin + step * i)))
-    return result
-
 
 class IOTest(unittest.TestCase):
     """
@@ -43,18 +11,12 @@ class IOTest(unittest.TestCase):
     def test_SaveNumpyArray_ReadDatafield(self):
         arr = numpy.array([[0, 1, 2, 3.0], [4, 5, 6, 7.0], [8, 9, 10, 11.0]])
         numpy.savetxt('tmp.txt', arr)
-        newdata = ba.readData2D("tmp.txt")
-        self.assertTrue(numpy.array_equal(newdata.npArray(), arr))
-
-    def test_SaveNumpyArray_ReadRawDataVector(self):
-        arr = numpy.array([[0, 1, 2, 3.0], [4, 5, 6, 7.0], [8, 9, 10, 11.0]])
-        numpy.savetxt('tmp.txt', arr)
-        newdata = numpy.array(
-            ba.readData2D("tmp.txt").flatVector())
-        expected = numpy.array(
-            [8., 4., 0., 9., 5., 1., 10., 6., 2., 11., 7., 3.])
-        self.assertTrue(numpy.array_equal(newdata, expected))
-
+        data2 = ba.readData2D("tmp.txt")
+        arr2 = data2.npArray()
+        print("arr = ", arr)
+        print("arr2 = ", arr2)
+        print("data2.flat = ", data2.flatVector())
+        self.assertTrue(numpy.array_equal(arr2, arr))
 
 if __name__ == '__main__':
     unittest.main()

Tests/Unit/Device/ArrayUtilsTest.cpp

@@ -8,14 +8,12 @@
 TEST(ArrayUtilsTest, DatafieldToVector2D)
 {
     Datafield data(
-        {newEquiDivision("axis0", 4, 10.0, 20.0), newEquiDivision("axis1", 3, 30.0, 40.0)});
-    const std::vector<double> values = {8.0,  4.0, 0.0, 9.0,  5.0, 1.0,
-                                        10.0, 6.0, 2.0, 11.0, 7.0, 3.0};
+        {newEquiDivision("axis0", 3, 10.0, 20.0), newEquiDivision("axis1", 2, 30.0, 40.0)});
+    const std::vector<double> values = {0., 1., 2., 10., 11., 12.};
 
     data.setVector(values);
 
     auto vec = DataUtil::createVector2D(data);
-    const std::vector<std::vector<double>> expected = {
-        {0.0, 1.0, 2.0, 3.0}, {4.0, 5.0, 6.0, 7.0}, {8.0, 9.0, 10.0, 11.0}};
+    const std::vector<std::vector<double>> expected = {{0., 1., 2.}, {10., 11., 12.}};
     EXPECT_EQ(vec, expected);
 }