#16 Added unit test for io

src/pybind_kernels/Makefile

@@ -46,5 +46,8 @@ $(foreach PLATFORM, $(PLATFORMS), \
 	) \
 )
 
+test: all
+	python3 -m pytest test
+
 clean:
 	rm -f $(TARGETS)

src/pybind_kernels/cpu/morphology.cc

@@ -33,7 +33,7 @@ void morphology_3d_sphere(
                             // TODO exact match with ndimage
                             bool within = px*px + py*py + pz*pz <= sqradius; // sphere kernel
                             int64_t offset = pz*strides[0] + py*strides[1] + px*strides[2];
-                            value = within? op(value, voxels[flat_index+offset]) : value;
+                            value = within ? op(value, voxels[flat_index+offset]) : value;
                         }
                     }
                 }

src/pybind_kernels/cpu_seq/morphology.cc

@@ -25,14 +25,13 @@ void morphology_3d_sphere(
 
                 // Apply the spherical kernel
                 bool value = neutral;
-                //#pragma omp simd collapse(3) reduction(op:value)
                 for (int64_t pz = limits[0]; pz <= limits[1]; pz++) {
                     for (int64_t py = limits[2]; py <= limits[3]; py++) {
                         for (int64_t px = limits[4]; px <= limits[5]; px++) {
                             // TODO exact match with ndimage
                             bool within = px*px + py*py + pz*pz <= sqradius; // sphere kernel
                             int64_t offset = pz*strides[0] + py*strides[1] + px*strides[2];
-                            value = within? op(value, voxels[flat_index+offset]) : value;
+                            value = within ? op(value, voxels[flat_index+offset]) : value;
                         }
                     }
                 }

src/pybind_kernels/gpu/morphology.cc

@@ -12,7 +12,7 @@ void morphology_3d_sphere(
     Op op;
     int64_t sqradius = radius * radius;
 
-    #pragma acc data copyin(voxels[:Nz*Ny*Nx], N[:3], strides[:3], sqradius) copyout(result[:Nz*Ny*Nx])
+    #pragma acc data copyin(voxels[:N[0]*N[1]*N[2]], N[:3], strides[:3], sqradius) copyout(result[:N[0]*N[1]*N[2]])
     {
         #pragma acc parallel loop collapse(3)
         for (int64_t z = 0; z < N[0]; z++) {
@@ -35,7 +35,7 @@ void morphology_3d_sphere(
                             for (int64_t px = limits[4]; px <= limits[5]; px++) {
                                 bool within = px*px + py*py + pz*pz <= sqradius; // sphere kernel
                                 int64_t offset = pz*strides[0] + py*strides[1] + px*strides[2];
-                                value = within? op(value, voxels[flat_index+offset]) : value;
+                                value = within ? op(value, voxels[flat_index+offset]) : value;
                             }
                         }
                     }

src/pybind_kernels/pybind/io-pybind.cc

@@ -34,15 +34,17 @@ void write_slice(const py::array_t<T> &np_data,
 
 PYBIND11_MODULE(io, m) {
     m.doc() = "I/O functions for handling flat binary format files."; // optional module docstring
-    m.def("load_slice", &load_slice<mask_type>);
-    m.def("load_slice", &load_slice<voxel_type>);
-    m.def("load_slice", &load_slice<field_type>);
-    m.def("load_slice", &load_slice<gauss_type>);
-    m.def("load_slice", &load_slice<real_t>);
+    m.def("load_slice", &load_slice<uint8_t>);
+    m.def("load_slice", &load_slice<uint16_t>);
+    m.def("load_slice", &load_slice<uint32_t>);
+    m.def("load_slice", &load_slice<uint64_t>);
+    m.def("load_slice", &load_slice<float>);
+    m.def("load_slice", &load_slice<double>);
 
-    m.def("write_slice", &write_slice<mask_type>);
-    m.def("write_slice", &write_slice<voxel_type>);
-    m.def("write_slice", &write_slice<field_type>);
-    m.def("write_slice", &write_slice<gauss_type>);
-    m.def("write_slice", &write_slice<real_t>);
+    m.def("write_slice", &write_slice<uint8_t>);
+    m.def("write_slice", &write_slice<uint16_t>);
+    m.def("write_slice", &write_slice<uint32_t>);
+    m.def("write_slice", &write_slice<uint64_t>);
+    m.def("write_slice", &write_slice<float>);
+    m.def("write_slice", &write_slice<double>);
 }

src/pybind_kernels/test/test_io.py

@@ -0,0 +1,66 @@
+'''
+Unittests for the I/O pybind kernels.
+'''
+import sys
+sys.path.append(sys.path[0]+"/../")
+import cpu_seq.io as io
+import numpy as np
+import tempfile
+import os
+import pytest
+
+# TODO np.bool doesn't work. It works when writing, but numpy doesn't recognize that the memory has been updated. It works fine if data_read is a np.uint8 array, even though an np.bool array has been written.
+dtypes_to_test = [np.uint8, np.uint16, np.uint32, np.uint64, np.float32, np.float64]
+tmp_folder = tempfile._get_default_tempdir()
+tmp_filename = next(tempfile._get_candidate_names())
+tmp_file = f'{tmp_folder}/{tmp_filename}'
+dim_size = 16
+dim_shape = (dim_size, dim_size, dim_size)
+partial_factor = 4
+
+def random(shape, dtype):
+    rnds = np.random.random(shape) * 100
+    return rnds > .5 if dtype == np.bool else rnds.astype(dtype)
+
+@pytest.mark.parametrize("dtype", dtypes_to_test)
+def test_dtype(dtype):
+    individual_tmp_file = f'{tmp_file}.{dtype.__name__}'
+    data = random(dim_shape, dtype)
+    data[0,0,1] = False
+    partial = dim_size // partial_factor
+
+    # Write out a new file
+    io.write_slice(data, individual_tmp_file, (0,0,0), dim_shape)
+    assert os.path.getsize(individual_tmp_file) == data.nbytes
+
+    # Read back and verify in chunks
+    read_data = np.zeros((partial, dim_size, dim_size), dtype=dtype)
+    for i in range(partial_factor):
+        io.load_slice(read_data, individual_tmp_file, (i*partial,0,0), read_data.shape)
+        assert np.allclose(data[i*partial:(i+1)*partial], read_data)
+
+    # Append another layer
+    data = np.append(data, random((partial, dim_size, dim_size), dtype), axis=0)
+    io.write_slice(data[dim_size:], individual_tmp_file, (dim_size,0,0), data.shape)
+    assert os.path.getsize(individual_tmp_file) == data.nbytes
+
+    # Read back and verify in chunks
+    for i in range(partial_factor+1):
+        io.load_slice(read_data, individual_tmp_file, (i*partial,0,0), read_data.shape)
+        assert np.allclose(data[i*partial:(i+1)*partial], read_data)
+
+    # Overwrite one of the "middle" chunks
+    data[partial:2*partial] = random((partial, dim_size, dim_size), dtype)
+    io.write_slice(data[partial:partial*2], individual_tmp_file, (partial,0,0), data.shape)
+
+    # Read back and verify in chunks
+    for i in range(partial_factor+1):
+        io.load_slice(read_data, individual_tmp_file, (i*partial,0,0), read_data.shape)
+        assert np.allclose(data[i*partial:(i+1)*partial], read_data)
+
+    os.remove(individual_tmp_file)
+
+if __name__ == '__main__':
+    for dtype in dtypes_to_test:
+        print (f'Testing {dtype.__name__}')
+        test_dtype(dtype)