Merge branch 'develop' into 'main'

Develop

See merge request !14

config/config_keys.toml

@@ -236,11 +236,11 @@ is_internal = [false]
 
 [LAMBDA]
 default = [0]
-value_type = ['int | double']
-value_format = ['N']
-value_N = [1]
-description = """This KEY controls the regularisation parameter :math:`\\lambda` for both :cpp:class:`M_BLR` and :cpp:class:`M_KRR`. If :code:`N=0`, no regularisation is applied. If :code:`N>0`, :math:`\\lambda` is set to this value. If :code:`N<0`, an evidence-approximation algorithm estimates the value of :math:`\\lambda`."""
-examples = ['LAMBDA -1', 'LAMBDA 1e-4']
+value_type = ['int | double | int double']
+value_format = ['N | N M']
+value_N = [-2]
+description = """This KEY controls the regularisation parameter :math:`\\lambda` for both :cpp:class:`M_BLR` and :cpp:class:`M_KRR`. If :code:`N=0`, no regularisation is applied. If :code:`N>0`, :math:`\\lambda` is set to this value. If :code:`N<0`, an evidence-approximation algorithm estimates the value of :math:`\\lambda`. For LAMBDA 0, second number, M (double), can be specified, i.e. LAMBDA 0 1e-12, which helps to determine the effective rank of matrix ::math:`\\Phi` (default 1e-8)."""
+examples = ['LAMBDA -1', 'LAMBDA 1e-4', 'LAMBDA 0', 'LAMBDA 0 1e-12']
 is_internal = [false]
 
 [ALPHA]

include/tadah/core/aeds_type.h

@@ -8,9 +8,13 @@
 #include <iostream>
 
 class aeds_type2 {
+  /* Column major order
+   *
+   * Each column represents single aed_type of dimension dim_
+   */
     double* data_ = nullptr;
-    size_t dim_ = 0;
-    size_t natoms_ = 0;
+    size_t dim_ = 0;  // nrows
+    size_t natoms_ = 0; // ncols
     aed_type** aptr = nullptr;
 
 public:
@@ -24,27 +28,27 @@ public:
     aeds_type2& operator=(aeds_type2&& other) noexcept;
 
     inline double& operator()(size_t r, size_t c) {
-        return data_[r * natoms_ + c];
+        return data_[c * dim_ + r];
     }
 
     inline const double& operator()(size_t r, size_t c) const {
-        return data_[r * natoms_ + c];
+        return data_[c * dim_ + r];
     }
 
-    inline aed_type& operator[](size_t i) {
-        return *aptr[i];
+    inline aed_type& operator[](size_t c) {
+        return *aptr[c];
     }
 
-    inline const aed_type& operator[](size_t i) const {
-        return *aptr[i];
+    inline const aed_type& operator[](size_t c) const {
+        return *aptr[c];
     }
 
-    inline aed_type& col(size_t i) {
-        return *aptr[i];
+    inline aed_type& col(size_t c) {
+        return *aptr[c];
     }
 
-    inline const aed_type& col(size_t i) const {
-        return *aptr[i];
+    inline const aed_type& col(size_t c) const {
+        return *aptr[c];
     }
 
     inline size_t cols() const {

include/tadah/core/utils.h

@@ -9,6 +9,8 @@
 #include <sstream>
 #include <iomanip>
 #include <random>
+#include <tuple>
+#include <cstddef>
 
 
 /** Vector printing to streams */
@@ -283,4 +285,33 @@ void reverse_columns(Matrix_Base<T> &matrix, const std::vector<size_t>& column_m
  * @param original_map A mapping of the rearranged indices to original positions.
  */
 void reverse_vector(t_type& vec, const std::vector<size_t>& original_map);
+
+/**
+ * @typedef histogram_bin
+ * @brief Type alias for a histogram bin containing the bin center coordinates and the count.
+ *
+ * Each histogram bin is represented as a tuple containing:
+ * - `double x_center`: The center coordinate of the bin along the x-axis.
+ * - `double y_center`: The center coordinate of the bin along the y-axis.
+ * - `size_t count`: The number of data points that fall into the bin.
+ */
+typedef std::tuple<double, double, size_t> histogram_bin;
+
+/**
+ * @brief Generates a 2D histogram from input data.
+ *
+ * This function takes a 2D array of doubles representing (x, y) data points and generates
+ * a 2D histogram. The histogram is returned as a vector of tuples, where each tuple contains
+ * the x and y coordinates of the bin center and the count of data points within that bin.
+ * The minimum and maximum values along each axis are calculated from the input data.
+ *
+ * @param data_2d A vector containing the (x, y) data points.
+ * @param x_bins The number of bins along the x-axis.
+ * @param y_bins The number of bins along the y-axis.
+ * @return A vector of histogram bins, each containing the bin center coordinates and the count.
+ */
+std::vector<histogram_bin> generate_2d_histogram(
+    const std::vector<std::pair<double,double>>& data_2d,
+    int x_bins,
+    int y_bins);
 #endif

src/aeds_type.cpp

@@ -109,7 +109,7 @@ void aeds_type2::set_zero() {
 std::ostream& operator<<(std::ostream& os, const aeds_type2& v) {
     for (size_t i = 0; i < v.dim_; ++i) {
         for (size_t j = 0; j < v.natoms_; ++j) {
-            os << v.data_[i * v.natoms_ + j] << " ";
+            os << v.data_[j * v.dim_ + i] << " ";
         }
         os << std::endl;
     }

src/utils.cpp

@@ -5,6 +5,8 @@
 #include <iostream>
 #include <algorithm>
 #include <cctype>
+#include <cmath>
+#include <limits>
 
 v_type logspace(double start, double stop, int num, double base) {
     if (start==0) {
@@ -235,3 +237,72 @@ void reverse_vector(t_type& vec, const std::vector<size_t>& original_map) {
     }
 }
 
+std::vector<histogram_bin> generate_2d_histogram(
+    const std::vector<std::pair<double,double>>& data_2d,
+    int x_bins,
+    int y_bins)
+{
+    if (data_2d.empty()) {
+        return {};
+    }
+
+    // Initialize minimum and maximum values
+    double x_min = std::numeric_limits<double>::max();
+    double x_max = std::numeric_limits<double>::lowest();
+    double y_min = std::numeric_limits<double>::max();
+    double y_max = std::numeric_limits<double>::lowest();
+
+    // Find min and max values along x and y axes from the data
+    for (const auto& point : data_2d) {
+        double x_value = point.first;
+        double y_value = point.second;
+
+        if (x_value < x_min) x_min = x_value;
+        if (x_value > x_max) x_max = x_value;
+        if (y_value < y_min) y_min = y_value;
+        if (y_value > y_max) y_max = y_value;
+    }
+
+    // Check if min and max values are valid
+    if (x_min >= x_max || y_min >= y_max) {
+        // Invalid data range
+        return {};
+    }
+
+    // Initialize the histogram bins with zero counts
+    std::vector<std::vector<size_t>> histogram(x_bins, std::vector<size_t>(y_bins, 0));
+
+    // Calculate the width of each bin along x and y axes
+    double x_bin_width = (x_max - x_min) / x_bins;
+    double y_bin_width = (y_max - y_min) / y_bins;
+
+    // Iterate over the data points and increment the corresponding histogram bins
+    for (const auto& point : data_2d) {
+        double x_value = point.first;
+        double y_value = point.second;
+
+        // Calculate the bin indices for the current data point
+        int x_bin = static_cast<int>((x_value - x_min) / x_bin_width);
+        int y_bin = static_cast<int>((y_value - y_min) / y_bin_width);
+
+        // Handle the edge case where the data point is exactly at the maximum value
+        if (x_bin == x_bins) x_bin = x_bins - 1;
+        if (y_bin == y_bins) y_bin = y_bins - 1;
+
+        // Increment the count in the appropriate bin
+        ++histogram[x_bin][y_bin];
+    }
+
+    // Prepare the output vector containing histogram bins with centers and counts
+    std::vector<histogram_bin> histogram_bins;
+    for (int ix = 0; ix < x_bins; ++ix) {
+        double x_center = x_min + (ix + 0.5) * x_bin_width;
+        for (int iy = 0; iy < y_bins; ++iy) {
+            double y_center = y_min + (iy + 0.5) * y_bin_width;
+            size_t count = histogram[ix][iy];
+            histogram_bins.push_back(std::make_tuple(x_center, y_center, count));
+        }
+    }
+
+    return histogram_bins;
+}

tests/test_periodic_table.cpp

+#include "catch2/catch.hpp"
+#include <tadah/core/periodic_table.h>
+
+TEST_CASE("PeriodicTable Initialization", "[initialize]") {
+    // Ensure that the periodic table is properly initialized
+    REQUIRE(PeriodicTable::symbol.empty());
+    REQUIRE(PeriodicTable::name.empty());
+    REQUIRE(PeriodicTable::Z.empty());
+    REQUIRE(PeriodicTable::masses.empty());
+
+    PeriodicTable::initialize();
+
+    REQUIRE_FALSE(PeriodicTable::symbol.empty());
+    REQUIRE_FALSE(PeriodicTable::name.empty());
+    REQUIRE_FALSE(PeriodicTable::Z.empty());
+    REQUIRE_FALSE(PeriodicTable::masses.empty());
+}
+
+TEST_CASE("Find Element by Symbol", "[find_by_symbol]") {
+    PeriodicTable::initialize();
+
+    SECTION("Known symbol") {
+        const Element& oxygen = PeriodicTable::find_by_symbol("O");
+        REQUIRE(oxygen.symbol == "O");
+        REQUIRE(oxygen.Z == 8);
+    }
+
+    SECTION("Unknown symbol") {
+        REQUIRE_THROWS(
+            PeriodicTable::find_by_symbol("Xx")
+        );
+    }
+}
+
+TEST_CASE("Find Element by Name", "[find_by_name]") {
+    PeriodicTable::initialize();
+
+    SECTION("Known name") {
+        const Element& carbon = PeriodicTable::find_by_name("Carbon");
+        REQUIRE(carbon.symbol == "C");
+        REQUIRE(carbon.Z == 6);
+    }
+
+    SECTION("Unknown name") {
+        REQUIRE_THROWS(
+            PeriodicTable::find_by_name("Unobtainium")
+        );
+    }
+}
+
+TEST_CASE("Find Element by Atomic Number (Z)", "[find_by_Z]") {
+    PeriodicTable::initialize();
+
+    SECTION("Known atomic number") {
+        const Element& iron = PeriodicTable::find_by_Z(26);
+        REQUIRE(iron.symbol == "Fe");
+        REQUIRE(iron.Z == 26);
+    }
+
+    SECTION("Unknown atomic number") {
+        REQUIRE_THROWS(
+            PeriodicTable::find_by_Z(999)
+        );
+    }
+}
+
+TEST_CASE("Get Atomic Mass", "[get_mass]") {
+    PeriodicTable::initialize();
+
+    SECTION("By atomic number") {
+        double mass_O = PeriodicTable::get_mass(8);
+        REQUIRE(mass_O == Approx(15.9990).epsilon(0.0001));
+    }
+
+    SECTION("By Element object") {
+        const Element& gold = PeriodicTable::find_by_symbol("Au");
+        double mass_Au = PeriodicTable::get_mass(gold);
+        REQUIRE(mass_Au == Approx(196.967).epsilon(0.0001));
+    }
+
+    SECTION("Unknown atomic number") {
+        REQUIRE_THROWS(
+            PeriodicTable::get_mass(999)
+        );
+    }
+}
+