support for AeroMode instantiation with sampled mode type (incl. fixe…

…s in spec file logic) + sampled and mono mode types coverage in unit tests (in test_aero_[mode,state,dist].py) (#357)

Co-authored-by: Jeffrey Curtis <[email protected]>

src/aero_mode.F90

@@ -244,4 +244,40 @@ subroutine f_aero_mode_get_name(ptr_c, name_data, name_size) bind(C)
     name_data = c_loc(aero_mode%name)
     name_size = len_trim(aero_mode%name)
   end subroutine
+
+  subroutine f_aero_mode_get_sample_num_conc(ptr_c, arr_data, data_size) bind(C)
+    type(c_ptr), intent(in) :: ptr_c
+    type(aero_mode_t), pointer :: aero_mode => null()
+    integer(c_int), intent(in) :: data_size
+    real(c_double), dimension(data_size), intent(inout) :: arr_data
+
+    call c_f_pointer(ptr_c, aero_mode)
+
+    arr_data = aero_mode%sample_num_conc
+
+  end subroutine
+
+  subroutine f_aero_mode_get_sample_radius(ptr_c, arr_data, data_size) bind(C)
+    type(c_ptr), intent(in) :: ptr_c
+    type(aero_mode_t), pointer :: aero_mode => null()
+    integer(c_int), intent(in) :: data_size
+    real(c_double), dimension(data_size), intent(inout) :: arr_data
+
+    call c_f_pointer(ptr_c, aero_mode)
+
+    arr_data = aero_mode%sample_radius
+
+  end subroutine
+
+  subroutine f_aero_mode_get_sample_bins(ptr_c, n_bins) bind(c)
+    type(c_ptr), intent(in) :: ptr_c
+    type(aero_mode_t), pointer :: aero_mode => null()
+    integer(c_int), intent(out) :: n_bins
+
+    call c_f_pointer(ptr_c, aero_mode)
+
+    n_bins = size(aero_mode%sample_num_conc)
+
+  end subroutine
+
 end module

src/aero_mode.hpp

@@ -118,6 +118,20 @@ extern "C" void f_aero_mode_from_json(
     void *aero_data_ptr
 ) noexcept;
 
+extern "C" void f_aero_mode_get_sample_num_conc(
+    const void *ptr,
+    void *sample_num_conc_data,
+    const int *sample_num_conc_data_size
+) noexcept;
+
+extern "C" void f_aero_mode_get_sample_radius(
+    const void *ptr,
+    void *sample_radius_data,
+    const int *sample_radius_data_size
+) noexcept;
+
+extern "C" void f_aero_mode_get_sample_bins(const void *ptr, int *n_bins) noexcept;
+
 struct AeroMode {
     PMCResource ptr;
 
@@ -145,6 +159,24 @@ struct AeroMode {
             throw std::runtime_error("mass_frac value must be a list of single-element dicts");
         if (!InputJSONResource::unique_keys(mass_frac))
             throw std::runtime_error("mass_frac keys must be unique");
+        if (mode["mode_type"] == "sampled") {
+            if (mode.find("size_dist") == mode.end())
+                throw std::runtime_error("size_dist key must be set for mode_type=sampled");
+            auto sd = mode["size_dist"];
+            if (
+                sd.size() != 2 || 
+                !sd[0].is_object() || 
+                sd[0].size() != 1 || 
+                sd[1].size() != 1 ||
+                sd[0].find("diam") == sd[0].end() ||
+                sd[1].find("num_conc") == sd[1].end()
+            )
+                throw std::runtime_error("size_dist value must be an iterable of two single-element dicts (first with 'diam', second with 'num_conc' as keys)");
+            auto diam = *sd[0].find("diam");
+            auto num_conc = *sd[1].find("num_conc");
+            if (diam.size() != num_conc.size() + 1)
+                throw std::runtime_error("size_dist['num_conc'] must have len(size_dist['diam'])-1 elements");
+        }
     }
 
     static auto get_num_conc(const AeroMode &self){
@@ -291,7 +323,7 @@ struct AeroMode {
         int type;
         f_aero_mode_get_type(self.ptr.f_arg(), &type);
 
-        if (type < 0 || (unsigned int)type >= AeroMode::types().size())
+        if (type <= 0 || (unsigned int)type > AeroMode::types().size())
             throw std::logic_error("Unknown mode type.");
 
         return AeroMode::types()[type - 1];
@@ -312,4 +344,29 @@ struct AeroMode {
           name[i] = f_ptr[i];
         return name;
     }
+
+    static auto get_sample_radius(const AeroMode &self) {
+       int len;
+       f_aero_mode_get_sample_bins(self.ptr.f_arg(), &len);
+       len++;
+       std::valarray<double> sample_radius(len);
+       f_aero_mode_get_sample_radius(
+          self.ptr.f_arg(),
+          begin(sample_radius),
+          &len
+       );
+       return sample_radius;
+    }
+
+    static auto get_sample_num_conc(const AeroMode &self) {
+       int len;
+       f_aero_mode_get_sample_bins(self.ptr.f_arg(), &len);
+       std::valarray<double> sample_num_conc(len);
+       f_aero_mode_get_sample_num_conc(
+           self.ptr.f_arg(),
+           begin(sample_num_conc),
+           &len
+       );
+       return sample_num_conc;
+    }
 };

src/json_resource.hpp

@@ -26,7 +26,7 @@ struct JSONResource {
     }
 
   protected:
-    size_t index = 0;
+    size_t index = 0, named_array_read_count = 0;
 
     JSONResource() {}
 
@@ -71,6 +71,10 @@ struct JSONResource {
   public:
     virtual ~JSONResource() {}
 
+    auto n_named_array_read_count() noexcept {
+        return this->named_array_read_count;
+    }
+
     void zoom_in(const bpstd::string_view &sub) noexcept {
         auto it = this->json->is_array() 
             ? this->json->at(this->json->size()-1).find(sub)
@@ -84,6 +88,7 @@ struct JSONResource {
         else
             this->set_current_json_ptr(&(*it));
 
+        this->named_array_read_count = 0;
     }
 
     void zoom_out() noexcept {
@@ -101,13 +106,12 @@ struct JSONResource {
         return this->json->begin();
     }
 
-    // TODO #112: to be removed after initialising GasData with a list, and not JSON?
-    auto first_field_name() const noexcept {
+    auto first_field_name() noexcept {
         // TODO #112: handle errors
         std::string name = "";
         assert(this->json->size() > 0);
         assert(this->json->begin()->size() > 0);
-        for (auto &entry : this->json->at(0).items())
+        for (auto &entry : this->json->at(this->named_array_read_count++).items())
         {
             name = entry.key();
         }

src/pypartmc.cpp

@@ -453,6 +453,10 @@ PYBIND11_MODULE(_PyPartMC, m) {
         .def_property("gsd", &AeroMode::get_gsd,
              &AeroMode::set_gsd, "Geometric standard deviation")
         .def("set_sample", &AeroMode::set_sampled)
+        .def_property_readonly("sample_num_conc", &AeroMode::get_sample_num_conc,
+             "Sample bin number concentrations (m^{-3})")
+        .def_property_readonly("sample_radius", &AeroMode::get_sample_radius,
+             "Sample bin radii (m).")
         .def_property("type", &AeroMode::get_type, &AeroMode::set_type,
              "Mode type (given by module constants)")
         .def_property("name", &AeroMode::get_name, &AeroMode::set_name,

src/spec_file_pypartmc.F90

@@ -129,6 +129,8 @@ subroutine spec_file_read_real_named_array(file, max_lines, names, vals)
             n_rows = max_lines
         end if
 
+        if (allocated(names)) deallocate(names)
+        if (allocated(vals)) deallocate(vals)
         allocate(names(n_rows))
         allocate(vals(n_rows, n_cols))
         allocate(vals_row(max(1, n_cols)))

src/spec_file_pypartmc.cpp

@@ -145,10 +145,10 @@ void spec_file_read_real_named_array_size(
     int *n_rows,
     int *n_cols
 ) noexcept {
-    auto first_field = json_resource_ptr()->first_field_name();
     *n_rows = json_resource_ptr()->n_numeric_array_entries();
+
+    auto first_field = json_resource_ptr()->first_field_name();
     *n_cols = json_resource_ptr()->n_elements(first_field);
-    // TODO #112: check each line has the same number of elements as time
 }
 
 extern "C"
@@ -174,17 +174,17 @@ void spec_file_read_real_named_array_data(
         ++i, ++it
     ) {
         assert(it->is_object());
-        if (i == row-1) {
+        if (i == (row - 1) + (json_resource_ptr()->n_named_array_read_count() - 1)) {
             assert(it->size() == 1);
             for (auto &entry : it->items()) {
-                // TODO #112: use input name_size as limit param
+                assert(*name_size > (long)entry.key().size());
                 for (auto c=0u; c < entry.key().size(); ++c)
                     name_data[c] = entry.key()[c];
                 *name_size = entry.key().size();
                 for (auto idx=0u; idx < entry.value().size(); ++idx) {
                     vals[idx] = entry.value().at(idx).get<double>();
                 }
-                break; // TODO #112
+                break;
             }
         }
     }

tests/test_aero_dist.py

@@ -217,3 +217,22 @@ def test_ctor_error_on_repeated_massfrac_keys():
 
         # assert
         assert str(exc_info.value) == "mass_frac keys must be unique"
+
+    @staticmethod
+    def test_ctor_sampled_mode():
+        # arrange
+        aero_data = ppmc.AeroData(AERO_DATA_CTOR_ARG_MINIMAL)
+        ctor_arg = copy.deepcopy(AERO_DIST_CTOR_ARG_MINIMAL)
+        ctor_arg[0]["test_mode"]["mode_type"] = "sampled"
+        ctor_arg[0]["test_mode"]["size_dist"] = [
+            {"diam": [1, 2, 3, 4]},
+            {"num_conc": [1, 2, 3]},
+        ]
+
+        # act
+        sut = ppmc.AeroDist(aero_data, ctor_arg)
+
+        # assert
+        assert sut.mode(0).num_conc == sum(
+            ctor_arg[0]["test_mode"]["size_dist"][1]["num_conc"]
+        )

tests/test_aero_mode.py

@@ -48,6 +48,18 @@
     }
 }
 
+AERO_MODE_CTOR_SAMPLED = {
+    "test_mode": {
+        "mass_frac": [{"H2O": [1]}],
+        "diam_type": "geometric",
+        "mode_type": "sampled",
+        "size_dist": [
+            {"diam": [1, 2, 3, 4]},
+            {"num_conc": [100, 200, 300]},
+        ],
+    }
+}
+
 
 class TestAeroMode:
     @staticmethod
@@ -289,3 +301,127 @@ def test_segfault_case():  # TODO #319
         )
         print(fishy_ctor_arg)
         ppmc.AeroMode(aero_data, fishy_ctor_arg)
+
+    @staticmethod
+    @pytest.mark.skipif(platform.machine() == "arm64", reason="TODO #348")
+    def test_sampled_without_size_dist():
+        # arrange
+        aero_data = ppmc.AeroData(AERO_DATA_CTOR_ARG_MINIMAL)
+        fishy_ctor_arg = copy.deepcopy(AERO_MODE_CTOR_LOG_NORMAL)
+        fishy_ctor_arg["test_mode"]["mode_type"] = "sampled"
+
+        # act
+        with pytest.raises(Exception) as exc_info:
+            ppmc.AeroMode(aero_data, fishy_ctor_arg)
+
+        # assert
+        assert str(exc_info.value) == "size_dist key must be set for mode_type=sampled"
+
+    @staticmethod
+    @pytest.mark.parametrize(
+        "fishy",
+        (
+            None,
+            [],
+            [{}, {}, {}],
+            [{}, []],
+            [{"diam": None}, {}],
+            [{"num_conc": None}, {}],
+            [{"diam": None, "": None}, {}],
+            [{"num_conc": None, "": None}, {}],
+        ),
+    )
+    @pytest.mark.skipif(platform.machine() == "arm64", reason="TODO #348")
+    def test_sampled_with_fishy_size_dist(fishy):
+        # arrange
+        aero_data = ppmc.AeroData(AERO_DATA_CTOR_ARG_MINIMAL)
+        fishy_ctor_arg = copy.deepcopy(AERO_MODE_CTOR_LOG_NORMAL)
+        fishy_ctor_arg["test_mode"]["mode_type"] = "sampled"
+        fishy_ctor_arg["test_mode"]["size_dist"] = fishy
+
+        # act
+        with pytest.raises(Exception) as exc_info:
+            ppmc.AeroMode(aero_data, fishy_ctor_arg)
+
+        # assert
+        assert (
+            str(exc_info.value)
+            == "size_dist value must be an iterable of two single-element dicts"
+            + " (first with 'diam', second with 'num_conc' as keys)"
+        )
+
+    @staticmethod
+    @pytest.mark.skipif(platform.machine() == "arm64", reason="TODO #348")
+    def test_sampled_with_diam_of_different_len_than_num_conc():
+        # arrange
+        aero_data = ppmc.AeroData(AERO_DATA_CTOR_ARG_MINIMAL)
+        fishy_ctor_arg = copy.deepcopy(AERO_MODE_CTOR_LOG_NORMAL)
+        fishy_ctor_arg["test_mode"]["mode_type"] = "sampled"
+        fishy_ctor_arg["test_mode"]["size_dist"] = [
+            {"diam": [1, 2, 3]},
+            {"num_conc": [1, 2, 3]},
+        ]
+
+        # act
+        with pytest.raises(Exception) as exc_info:
+            ppmc.AeroMode(aero_data, fishy_ctor_arg)
+
+        # assert
+        assert (
+            str(exc_info.value)
+            == "size_dist['num_conc'] must have len(size_dist['diam'])-1 elements"
+        )
+
+    @staticmethod
+    def test_sampled():
+        # arrange
+        aero_data = ppmc.AeroData(AERO_DATA_CTOR_ARG_MINIMAL)
+
+        # act
+        sut = ppmc.AeroMode(aero_data, AERO_MODE_CTOR_SAMPLED)
+
+        # assert
+        assert sut.type == "sampled"
+        assert sut.num_conc == np.sum(
+            AERO_MODE_CTOR_SAMPLED["test_mode"]["size_dist"][1]["num_conc"]
+        )
+        assert (
+            sut.sample_num_conc
+            == AERO_MODE_CTOR_SAMPLED["test_mode"]["size_dist"][1]["num_conc"]
+        )
+        assert (
+            np.array(sut.sample_radius) * 2
+            == AERO_MODE_CTOR_SAMPLED["test_mode"]["size_dist"][0]["diam"]
+        ).all()
+
+    @staticmethod
+    def test_set_sample():
+        # arrange
+        aero_data = ppmc.AeroData(AERO_DATA_CTOR_ARG_MINIMAL)
+
+        diams = [1, 2, 3, 4]
+        num_concs = [100, 200, 300]
+        sut = ppmc.AeroMode(
+            aero_data,
+            {
+                "test_mode": {
+                    "mass_frac": [{"H2O": [1]}],
+                    "diam_type": "geometric",
+                    "mode_type": "sampled",
+                    "size_dist": [
+                        {"diam": diams},
+                        {"num_conc": num_concs},
+                    ],
+                }
+            },
+        )
+        num_conc_orig = sut.num_conc
+        # act
+        diams = [0.5 * x for x in diams]
+        num_concs = [2 * x for x in num_concs]
+        sut.set_sample(diams, num_concs)
+        # assert
+        assert sut.num_conc == np.sum(num_concs)
+        assert sut.sample_num_conc == num_concs
+        assert (np.array(sut.sample_radius) * 2 == diams).all()
+        assert sut.num_conc == num_conc_orig * 2