Merge pull request #110 from jcurtis2/add_histograms

add histograms

Author: slayoo

src/bin_grid.F90

@@ -70,4 +70,51 @@ subroutine f_bin_grid_centers(ptr_c, arr_data, arr_size) bind(C)
     arr_data = bin_grid%centers
   end subroutine
 
+  subroutine f_bin_grid_histogram_1d(x_bin_grid_ptr_c, x_data, weight_data, &
+       arr_size, output_data, bin_grid_size) bind(C)
+
+    type(c_ptr), intent(in) :: x_bin_grid_ptr_c
+    type(bin_grid_t), pointer :: bin_grid => null()
+    integer(c_int), intent(in) :: arr_size
+    integer(c_int), intent(in) :: bin_grid_size
+    real(c_double), dimension(bin_grid_size) :: output_data
+    real(c_double), dimension(arr_size) :: x_data
+    real(c_double), dimension(arr_size) :: weight_data
+
+    call c_f_pointer(x_bin_grid_ptr_c, bin_grid)
+    output_data = bin_grid_histogram_1d(bin_grid, x_data, weight_data)
+
+  end subroutine
+
+  subroutine f_bin_grid_histogram_2d(x_bin_grid_ptr_c, x_data, &
+       y_bin_grid_ptr_c, y_data, weight_data, &
+       arr_size, output_data, x_bin_grid_size, y_bin_grid_size) bind(C)
+
+    type(c_ptr), intent(in) :: x_bin_grid_ptr_c
+    type(c_ptr), intent(in) :: y_bin_grid_ptr_c
+    type(bin_grid_t), pointer :: x_bin_grid => null()
+    type(bin_grid_t), pointer :: y_bin_grid => null()
+    integer(c_int), intent(in) :: arr_size
+    integer(c_int), intent(in) :: x_bin_grid_size
+    integer(c_int), intent(in) :: y_bin_grid_size
+    real(c_double), dimension(x_bin_grid_size*y_bin_grid_size), intent(out) :: output_data
+    real(c_double), dimension(arr_size), intent(in) :: x_data, y_data
+    real(c_double), dimension(arr_size), intent(in) :: weight_data
+    real(c_double), allocatable :: output_data_local(:,:)
+    integer :: i, j
+
+    call c_f_pointer(x_bin_grid_ptr_c, x_bin_grid)
+    call c_f_pointer(y_bin_grid_ptr_c, y_bin_grid)
+    allocate(output_data_local(x_bin_grid_size,y_bin_grid_size))
+    output_data_local = bin_grid_histogram_2d(x_bin_grid, x_data, y_bin_grid, & 
+         y_data, weight_data)
+
+    do i = 1,x_bin_grid_size
+    do j = 1,y_bin_grid_size
+       output_data((i-1)*y_bin_grid_size + j) = output_data_local(i,j)
+    end do
+    end do
+
+  end subroutine
+
 end module

src/bin_grid.hpp

@@ -21,6 +21,12 @@ extern "C" void f_bin_grid_init(
 extern "C" void f_bin_grid_size(const void *ptr, int *val) noexcept;
 extern "C" void f_bin_grid_edges(const void *ptr, void *arr_data, const int *arr_size) noexcept;
 extern "C" void f_bin_grid_centers(const void *ptr, void *arr_data, const int *arr_size) noexcept;
+extern "C" void f_bin_grid_histogram_1d(const void *x_bin_grid_ptr_c, const void *x_data,
+    const void *weight_data, const int *arr_size, void *output_data, const int *grid_size) noexcept;
+extern "C" void f_bin_grid_histogram_2d(const void *x_bin_grid_ptr_c, const void *x_data,
+    const void *y_bin_grid_ptr_c, const void *y_data,
+    const void *weight_data, const int *arr_size, void *output_data, const int *x_grid_size,
+    const int *y_grid_size) noexcept;
 
 struct BinGrid {
     PMCResource ptr;
@@ -61,3 +67,48 @@ struct BinGrid {
         return data;
     }
 };
+
+static std::valarray<double> histogram_1d(
+    const BinGrid &bin_grid,
+    std::valarray<double> values,
+    std::valarray<double> weights){
+
+        int len;
+        f_bin_grid_size(&bin_grid.ptr, &len);
+        int data_size = values.size();
+        std::valarray<double> data(len);
+        f_bin_grid_histogram_1d(&bin_grid.ptr, begin(values), begin(weights),
+            &data_size, begin(data), &len);
+
+        return data;
+}
+
+static std::vector<std::vector<double>> histogram_2d(
+    const BinGrid &x_bin_grid,
+    std::valarray<double> x_values,
+    const BinGrid &y_bin_grid,
+    std::valarray<double> y_values,
+    std::valarray<double> weights){
+
+        int x_len;
+        f_bin_grid_size(&x_bin_grid.ptr, &x_len);
+        int y_len;
+        f_bin_grid_size(&y_bin_grid.ptr, &y_len);
+        int data_size = x_values.size();
+
+        std::vector<std::vector<double>> data( x_len , std::vector<double>(y_len,0));
+        std::valarray<double> data_flat(x_len*y_len);
+        f_bin_grid_histogram_2d(&x_bin_grid.ptr, begin(x_values),
+            &y_bin_grid.ptr, begin(y_values), begin(weights),
+            &data_size, begin(data_flat), &x_len, &y_len);
+
+        for(int i = 0; i < x_len; i++)
+        {
+            for(int j = 0; j < y_len; j++)
+            {
+            data[i][j] = data_flat[i*y_len + j];
+            }
+        }
+
+        return data;
+}

src/pypartmc.cpp

@@ -176,9 +176,20 @@ PYBIND11_MODULE(_PyPartMC, m) {
     py::class_<BinGrid>(m,"BinGrid")
         .def(py::init<const double, const py::str, const double, const double>())
         .def("__len__", BinGrid::__len__)
-        .def_property_readonly("edges", BinGrid::edges)
-        .def_property_readonly("centers", BinGrid::centers)
+        .def_property_readonly("edges", BinGrid::edges, "Bin edges")
+        .def_property_readonly("centers", BinGrid::centers, "Bin centers")
     ;
+
+    m.def(
+        "histogram_1d", &histogram_1d, py::return_value_policy::copy,
+        "Return a 1D histogram with of the given weighted data, scaled by the bin sizes."
+    );
+
+    m.def(
+        "histogram_2d", &histogram_2d, py::return_value_policy::copy,
+        "Return a 2D histogram with of the given weighted data, scaled by the bin sizes."
+    );
+
     //  TODO #120: auto util = m.def_submodule("util", "...");
     m.def(
         "pow2_above", &pow2_above, py::return_value_policy::copy,
@@ -199,6 +210,8 @@ PYBIND11_MODULE(_PyPartMC, m) {
         "Scenario",
         "condense_equilib_particles",
         "run_part",
-        "pow2_above"
+        "pow2_above",
+        "histogram_1d",
+        "histogram_2d"
     );
 }

tests/test_bin_grid.py

@@ -97,3 +97,94 @@ def test_invalid_grid():
             _ = ppmc.BinGrid(grid_size, "X", 1, 100)
         except ValueError as error:
             assert str(error) == "Invalid grid spacing."
+
+    @staticmethod
+    def test_histogram_1d():
+        # arrange
+        n_data = 1000
+        grid = ppmc.BinGrid(100,"linear",0,1000)
+        vals = np.random.random(n_data)*1000
+        weights = np.ones(n_data)
+        hist, bin_edges = np.histogram(vals, bins=grid.edges)
+
+        # act
+        data = ppmc.histogram_1d(grid, vals, weights)
+
+        # assert
+        np.testing.assert_array_almost_equal(data, hist/(bin_edges[1]-bin_edges[0]))
+
+    @staticmethod
+    def test_histogram_2d_linear_linear():
+        # arrange
+        n_data = 1000
+        x_grid = ppmc.BinGrid(15, "linear", 0, 1000)
+        y_grid = ppmc.BinGrid(12, "linear", 0, 500)
+        x_vals = np.random.random(n_data)*1000
+        y_vals = np.random.random(n_data)*500
+        weights = np.random.random(n_data)
+        data_numpy, bin_edges_x, bin_edges_y = np.histogram2d(
+            x_vals, y_vals,
+            bins=[x_grid.edges, y_grid.edges],
+            weights=weights
+        )
+        cell_size = (bin_edges_x[1]-bin_edges_x[0]) * (bin_edges_y[1]-bin_edges_y[0])
+
+        # act
+        data = ppmc.histogram_2d(x_grid, x_vals, y_grid, y_vals, weights)
+
+        # assert
+        np.testing.assert_array_almost_equal(np.array(data), data_numpy/cell_size, decimal=15)
+
+    @staticmethod
+    def test_histogram_2d_linear_log():
+        # arrange
+        n_data = 100
+        y_data_min = .1
+        y_data_max = 10.
+        x_grid = ppmc.BinGrid(15,"linear",0,1000)
+        y_grid = ppmc.BinGrid(12,"log",y_data_min,y_data_max)
+        x_vals = np.random.random(n_data)*1000
+        y_vals = y_data_min*10**(np.log10(y_data_max/y_data_min) * np.random.random(n_data))
+        weights = np.random.random(n_data)
+        data_numpy, bin_edges_x, bin_edges_y = np.histogram2d(
+            x_vals, y_vals,
+            bins=[x_grid.edges, y_grid.edges],
+            weights=weights
+        )
+        cell_size = (bin_edges_x[1]-bin_edges_x[0]) * np.log(bin_edges_y[1]/bin_edges_y[0])
+
+        # act
+        data = ppmc.histogram_2d(x_grid, x_vals, y_grid, y_vals, weights)
+
+        # assert
+        np.testing.assert_array_almost_equal(np.array(data), data_numpy/cell_size, decimal=15)
+
+    @staticmethod
+    def test_histogram_2d_log_log():
+        # arrange
+        n_data = 100
+        x_data_min = 1
+        x_data_max = 1000
+        y_data_min = .1
+        y_data_max = 10.
+        x_grid = ppmc.BinGrid(15,"log", x_data_min, x_data_max)
+        y_grid = ppmc.BinGrid(12,"log", y_data_min, y_data_max)
+        x_vals = x_data_min * 10**(np.log10(x_data_max / x_data_min) * np.random.random(n_data))
+        y_vals = y_data_min * 10**(np.log10(y_data_max / y_data_min) * np.random.random(n_data))
+        weights = np.random.random(n_data)
+        data_numpy, bin_edges_x, bin_edges_y = np.histogram2d(
+            x_vals, y_vals,
+            bins=[x_grid.edges, y_grid.edges],
+            weights=weights
+        )
+        cell_size = np.log(
+            bin_edges_x[1] / bin_edges_x[0]
+        ) * np.log(
+            bin_edges_y[1] / bin_edges_y[0]
+        )
+
+        # act
+        data = ppmc.histogram_2d(x_grid, x_vals, y_grid, y_vals, weights)
+
+        # assert
+        np.testing.assert_array_almost_equal(np.array(data), data_numpy/cell_size, decimal=13)