LCOV - code coverage report
Current view: top level - src/database - Database.cpp (source / functions) Hit Total Coverage
Test: coverage.info Lines: 444 560 79.3 %
Date: 2026-10-14 13:16:16 Functions: 31 37 83.8 %

          Line data    Source code
       1             : // SPDX-FileCopyrightText: 2024 PairInteraction Developers
       2             : // SPDX-License-Identifier: LGPL-3.0-or-later
       3             : 
       4             : #include "pairinteraction/database/Database.hpp"
       5             : 
       6             : #include "pairinteraction/basis/BasisAtom.hpp"
       7             : #include "pairinteraction/database/AtomDescriptionByParameters.hpp"
       8             : #include "pairinteraction/database/AtomDescriptionByRanges.hpp"
       9             : #include "pairinteraction/database/GitHubDownloader.hpp"
      10             : #include "pairinteraction/database/ParquetManager.hpp"
      11             : #include "pairinteraction/enums/OperatorType.hpp"
      12             : #include "pairinteraction/ket/KetAtom.hpp"
      13             : #include "pairinteraction/ket/KetNotUniqueError.hpp"
      14             : #include "pairinteraction/utils/TaskControl.hpp"
      15             : #include "pairinteraction/utils/hash.hpp"
      16             : #include "pairinteraction/utils/ket_id.hpp"
      17             : #include "pairinteraction/utils/paths.hpp"
      18             : #include "pairinteraction/utils/streamed.hpp"
      19             : 
      20             : #include <algorithm>
      21             : #include <cmath>
      22             : #include <cpptrace/cpptrace.hpp>
      23             : #include <duckdb.hpp>
      24             : #include <fmt/core.h>
      25             : #include <fmt/format.h>
      26             : #include <fmt/ranges.h>
      27             : #include <fstream>
      28             : #include <iterator>
      29             : #include <nlohmann/json.hpp>
      30             : #include <oneapi/tbb.h>
      31             : #include <spdlog/spdlog.h>
      32             : #include <string>
      33             : #include <system_error>
      34             : #include <unordered_map>
      35             : #include <unordered_set>
      36             : #include <utility>
      37             : 
      38             : namespace pairinteraction {
      39             : 
      40             : namespace {
      41             : 
      42             : // Absolute tolerance by which quantum number ranges are widened so that states are not dropped
      43             : // because of the numerical noise of the (expectation) values stored in the database. This noise
      44             : // stems from the expectation value calculation as well as from the nullspace solver tolerances
      45             : // (MQDT in rydstate).
      46             : constexpr double database_numerical_precision = 1e-9;
      47             : 
      48         542 : std::string format_quantum_number_range(const std::string &column, const Range<double> &range) {
      49           0 :     return fmt::format("{} BETWEEN {} AND {}", column, range.min() - database_numerical_precision,
      50        1084 :                        range.max() + database_numerical_precision);
      51             : }
      52             : 
      53         406 : std::string format_expectation_value_range(const std::string &value_column,
      54             :                                            const std::string &std_column,
      55             :                                            const Range<double> &range,
      56             :                                            double standard_deviation_factor) {
      57           0 :     return fmt::format("{} BETWEEN {}-{}*{}-{} AND {}+{}*{}+{}", value_column, range.min(),
      58             :                        standard_deviation_factor, std_column, database_numerical_precision,
      59         406 :                        range.max(), standard_deviation_factor, std_column,
      60         812 :                        database_numerical_precision);
      61             : }
      62             : 
      63             : // Find the index of the result column with the given name.
      64        4367 : size_t get_column_index(const std::vector<std::string> &names, const std::string &name) {
      65        4367 :     auto it = std::find(names.begin(), names.end(), name);
      66        4367 :     if (it == names.end()) {
      67           0 :         throw std::runtime_error("Missing database column '" + name + "'.");
      68             :     }
      69        4367 :     return static_cast<size_t>(std::distance(names.begin(), it));
      70             : }
      71             : 
      72        2335 : bool is_integer_like(duckdb::LogicalTypeId id) {
      73        2335 :     switch (id) {
      74         903 :     case duckdb::LogicalTypeId::BOOLEAN:
      75             :     case duckdb::LogicalTypeId::TINYINT:
      76             :     case duckdb::LogicalTypeId::SMALLINT:
      77             :     case duckdb::LogicalTypeId::INTEGER:
      78             :     case duckdb::LogicalTypeId::BIGINT:
      79             :     case duckdb::LogicalTypeId::UTINYINT:
      80             :     case duckdb::LogicalTypeId::USMALLINT:
      81             :     case duckdb::LogicalTypeId::UINTEGER:
      82             :     case duckdb::LogicalTypeId::UBIGINT:
      83         903 :         return true;
      84        1432 :     default:
      85        1432 :         return false;
      86             :     }
      87             : }
      88             : 
      89             : template <typename T>
      90     1532439 : T get_entry_as(duckdb::Vector &vector, const duckdb::LogicalType &type, size_t row) {
      91     1532439 :     switch (type.id()) {
      92     1388120 :     case duckdb::LogicalTypeId::DOUBLE:
      93     1388120 :         return static_cast<T>(duckdb::FlatVector::GetData<double>(vector)[row]);
      94           0 :     case duckdb::LogicalTypeId::FLOAT:
      95           0 :         return static_cast<T>(duckdb::FlatVector::GetData<float>(vector)[row]);
      96           0 :     case duckdb::LogicalTypeId::BOOLEAN:
      97           0 :         return static_cast<T>(duckdb::FlatVector::GetData<bool>(vector)[row] ? 1 : 0);
      98           0 :     case duckdb::LogicalTypeId::TINYINT:
      99           0 :         return static_cast<T>(duckdb::FlatVector::GetData<int8_t>(vector)[row]);
     100           0 :     case duckdb::LogicalTypeId::SMALLINT:
     101           0 :         return static_cast<T>(duckdb::FlatVector::GetData<int16_t>(vector)[row]);
     102           0 :     case duckdb::LogicalTypeId::INTEGER:
     103           0 :         return static_cast<T>(duckdb::FlatVector::GetData<int32_t>(vector)[row]);
     104      144318 :     case duckdb::LogicalTypeId::BIGINT:
     105      144318 :         return static_cast<T>(duckdb::FlatVector::GetData<int64_t>(vector)[row]);
     106           0 :     case duckdb::LogicalTypeId::UTINYINT:
     107           0 :         return static_cast<T>(duckdb::FlatVector::GetData<uint8_t>(vector)[row]);
     108           0 :     case duckdb::LogicalTypeId::USMALLINT:
     109           0 :         return static_cast<T>(duckdb::FlatVector::GetData<uint16_t>(vector)[row]);
     110           0 :     case duckdb::LogicalTypeId::UINTEGER:
     111           0 :         return static_cast<T>(duckdb::FlatVector::GetData<uint32_t>(vector)[row]);
     112           0 :     case duckdb::LogicalTypeId::UBIGINT:
     113           0 :         return static_cast<T>(duckdb::FlatVector::GetData<uint64_t>(vector)[row]);
     114           2 :     default:
     115           2 :         throw std::runtime_error("Cannot read database column of type " + type.ToString() +
     116             :                                  " as a quantum number.");
     117             :     }
     118             : }
     119             : 
     120     1484596 : double get_entry_as_double(duckdb::Vector &vector, const duckdb::LogicalType &type, size_t row) {
     121     1484596 :     return get_entry_as<double>(vector, type, row);
     122             : }
     123             : 
     124       47842 : int64_t get_entry_as_int64(duckdb::Vector &vector, const duckdb::LogicalType &type, size_t row) {
     125       47842 :     return get_entry_as<int64_t>(vector, type, row);
     126             : }
     127             : 
     128             : struct QuantumNumbers {
     129             :     std::unordered_map<std::string, double> values;
     130             :     std::unordered_map<std::string, double> stds;
     131             : };
     132             : 
     133       47842 : QuantumNumbers get_quantum_numbers_from_row(duckdb::DataChunk &chunk,
     134             :                                             const std::vector<duckdb::LogicalType> &types,
     135             :                                             const std::vector<std::string> &names,
     136             :                                             const std::unordered_set<std::string> &excluded_columns,
     137             :                                             size_t row) {
     138       47842 :     QuantumNumbers quantum_numbers;
     139     1626109 :     for (size_t col = 0; col < names.size(); ++col) {
     140     1578265 :         const std::string &name = names[col];
     141     1578265 :         if (excluded_columns.contains(name)) {
     142      143526 :             continue;
     143             :         }
     144     1434735 :         double value = get_entry_as_double(chunk.data[col], types[col], row);
     145     1434731 :         if (name.starts_with("std_")) {
     146      574101 :             quantum_numbers.stds[name.substr(4)] = value;
     147      860629 :         } else if (name.starts_with("exp_")) {
     148      574101 :             quantum_numbers.values[name.substr(4)] = value;
     149             :         } else {
     150      286536 :             quantum_numbers.values[name] = value;
     151             :         }
     152             :     }
     153       47842 :     return quantum_numbers;
     154           0 : }
     155             : 
     156        2487 : bool must_match_exactly(const std::unordered_map<std::string, duckdb::LogicalTypeId> &columns,
     157             :                         const std::string &column) {
     158        2487 :     if (column == "f" || column == "m" || column == "parity") {
     159         152 :         return true;
     160             :     }
     161        2335 :     auto it = columns.find(column);
     162        2335 :     return it != columns.end() && is_integer_like(it->second);
     163             : }
     164             : } // namespace
     165             : 
     166       47840 : void ensure_consistent_quantum_numbers(double quantum_number_f, double quantum_number_m) {
     167       47840 :     if (2 * quantum_number_m != std::rint(2 * quantum_number_m)) {
     168           0 :         throw std::runtime_error("The quantum number m must be an integer or half-integer.");
     169             :     }
     170       47840 :     if (2 * quantum_number_f != std::rint(2 * quantum_number_f)) {
     171           0 :         throw std::runtime_error("The quantum number f must be an integer or half-integer.");
     172             :     }
     173       47840 :     if (quantum_number_f + quantum_number_m != std::rint(quantum_number_f + quantum_number_m)) {
     174           0 :         throw std::invalid_argument(
     175           0 :             "The quantum numbers f and m must be both either integers or half-integers.");
     176             :     }
     177       47840 :     if (std::abs(quantum_number_m) > quantum_number_f) {
     178           1 :         throw std::invalid_argument(
     179           2 :             "The absolute value of the quantum number m must be less than or equal to f.");
     180             :     }
     181       47839 : }
     182             : 
     183           0 : Database::Database() : Database(default_download_missing) {}
     184             : 
     185           0 : Database::Database(bool download_missing)
     186           0 :     : Database(download_missing, default_use_cache, default_database_dir) {}
     187             : 
     188           0 : Database::Database(std::filesystem::path database_dir)
     189           0 :     : Database(default_download_missing, default_use_cache, std::move(database_dir)) {}
     190             : 
     191           3 : Database::Database(bool download_missing, bool use_cache, std::filesystem::path database_dir)
     192           3 :     : download_missing_(download_missing), use_cache_(use_cache),
     193           3 :       database_dir_(std::move(database_dir)), db(std::make_unique<duckdb::DuckDB>(nullptr)),
     194          15 :       con(std::make_unique<duckdb::Connection>(*db)) {
     195             : 
     196           3 :     if (database_dir_.empty()) {
     197           0 :         database_dir_ = default_database_dir;
     198             :     }
     199             : 
     200             :     // Ensure the database directory exists
     201           3 :     if (!std::filesystem::exists(database_dir_)) {
     202           0 :         std::filesystem::create_directories(database_dir_);
     203             :     }
     204           3 :     database_dir_ = std::filesystem::canonical(database_dir_);
     205           3 :     if (!std::filesystem::is_directory(database_dir_)) {
     206           0 :         throw std::filesystem::filesystem_error("Cannot access database", database_dir_.string(),
     207           0 :                                                 std::make_error_code(std::errc::not_a_directory));
     208             :     }
     209           3 :     SPDLOG_INFO("Using database directory: {}", database_dir_.string());
     210             : 
     211             :     // Ensure that the config directory exists
     212           3 :     std::filesystem::path configdir = paths::get_config_directory();
     213           3 :     if (!std::filesystem::exists(configdir)) {
     214           1 :         std::filesystem::create_directories(configdir);
     215           2 :     } else if (!std::filesystem::is_directory(configdir)) {
     216           0 :         throw std::filesystem::filesystem_error("Cannot access config directory ",
     217           0 :                                                 configdir.string(),
     218           0 :                                                 std::make_error_code(std::errc::not_a_directory));
     219             :     }
     220             : 
     221             :     // Read in the database_repo_paths if a config file exists, otherwise use the default and
     222             :     // write it to the config file
     223           3 :     std::filesystem::path configfile = configdir / "database.json";
     224           3 :     std::string database_repo_host;
     225           3 :     std::vector<std::string> database_repo_paths;
     226           3 :     if (std::filesystem::exists(configfile)) {
     227           2 :         std::ifstream file(configfile);
     228           2 :         nlohmann::json doc = nlohmann::json::parse(file, nullptr, false);
     229             : 
     230           4 :         if (!doc.is_discarded() && doc.contains("hash") && doc.contains("database_repo_host") &&
     231           2 :             doc.contains("database_repo_paths")) {
     232           2 :             database_repo_host = doc["database_repo_host"].get<std::string>();
     233           2 :             database_repo_paths = doc["database_repo_paths"].get<std::vector<std::string>>();
     234             : 
     235             :             // If the values are not equal to the default values but the hash is consistent (i.e.,
     236             :             // the user has not changed anything manually), clear the values so that they can be
     237             :             // updated
     238           4 :             if (database_repo_host != default_database_repo_host ||
     239           2 :                 database_repo_paths != default_database_repo_paths) {
     240           0 :                 std::size_t seed = 0;
     241           0 :                 utils::hash_combine(seed, database_repo_paths);
     242           0 :                 utils::hash_combine(seed, database_repo_host);
     243           0 :                 if (seed == doc["hash"].get<std::size_t>()) {
     244           0 :                     database_repo_host.clear();
     245           0 :                     database_repo_paths.clear();
     246             :                 } else {
     247           0 :                     SPDLOG_INFO("The database repository host and paths have been changed "
     248             :                                 "manually. Thus, they will not be updated automatically. To reset "
     249             :                                 "them, delete the file '{}'.",
     250             :                                 configfile.string());
     251             :                 }
     252             :             }
     253             :         }
     254           2 :     }
     255             : 
     256             :     // Read in and store the default values if necessary
     257           3 :     if (database_repo_host.empty() || database_repo_paths.empty()) {
     258           2 :         SPDLOG_INFO("Updating the database repository host and paths:");
     259             : 
     260           1 :         database_repo_host = default_database_repo_host;
     261           1 :         database_repo_paths = default_database_repo_paths;
     262           1 :         std::ofstream file(configfile);
     263           1 :         nlohmann::json doc;
     264             : 
     265           1 :         SPDLOG_INFO("* New host: {}", default_database_repo_host);
     266           2 :         SPDLOG_INFO("* New paths: {}", fmt::join(default_database_repo_paths, ", "));
     267             : 
     268           1 :         doc["database_repo_host"] = default_database_repo_host;
     269           1 :         doc["database_repo_paths"] = database_repo_paths;
     270             : 
     271           1 :         std::size_t seed = 0;
     272           1 :         utils::hash_combine(seed, default_database_repo_paths);
     273           1 :         utils::hash_combine(seed, default_database_repo_host);
     274           1 :         doc["hash"] = seed;
     275             : 
     276           1 :         file << doc.dump(4);
     277           1 :     }
     278             : 
     279             :     // Limit the memory usage of duckdb's buffer manager
     280             :     {
     281           3 :         auto result = con->Query("PRAGMA max_memory = '8GB';");
     282           3 :         if (result->HasError()) {
     283           0 :             throw cpptrace::runtime_error("Error setting the memory limit: " + result->GetError());
     284             :         }
     285           3 :     }
     286             : 
     287             :     // Instantiate a database manager that provides access to database tables. If a table
     288             :     // is outdated/not available locally, it will be downloaded if download_missing_ is true.
     289           3 :     if (!download_missing_) {
     290           3 :         database_repo_paths.clear();
     291             :     }
     292           3 :     downloader = std::make_unique<GitHubDownloader>();
     293           6 :     manager = std::make_unique<ParquetManager>(database_dir_, *downloader, database_repo_paths,
     294           6 :                                                *con, use_cache_);
     295           3 :     manager->scan_local();
     296           3 :     manager->scan_remote();
     297             : 
     298             :     // Print versions of tables
     299           3 :     std::istringstream iss(manager->get_versions_info());
     300          36 :     for (std::string line; std::getline(iss, line);) {
     301          33 :         SPDLOG_INFO(line);
     302           3 :     }
     303           9 : }
     304             : 
     305           2 : Database::~Database() = default;
     306             : 
     307        2187 : const Database::column_info_t &Database::get_states_column_info(const std::string &species) {
     308        2187 :     std::string table_path = manager->get_path(species, "states");
     309        2187 :     if (auto it = column_info_cache.find(table_path); it != column_info_cache.end()) {
     310        2177 :         return it->second;
     311             :     }
     312          20 :     auto result = con->Query(fmt::format(R"(SELECT * FROM '{}' LIMIT 0)", table_path));
     313          10 :     if (result->HasError()) {
     314           0 :         throw cpptrace::runtime_error("Error querying the database columns: " + result->GetError());
     315             :     }
     316          10 :     column_info_t columns;
     317          10 :     columns.reserve(result->names.size());
     318         320 :     for (size_t col = 0; col < result->names.size(); ++col) {
     319         310 :         columns.emplace(result->names[col], result->types[col].id());
     320             :     }
     321             : 
     322             :     // Every states table must provide the columns required for constructing kets. All other
     323             :     // columns are treated generically as quantum numbers, whatever they are named.
     324          50 :     for (const auto &required : {"id", "energy", "f", "parity"}) {
     325          40 :         if (!columns.contains(required)) {
     326           0 :             throw std::runtime_error(
     327           0 :                 fmt::format("The table of states of species '{}' is missing the required column "
     328             :                             "'{}'.",
     329           0 :                             species, required));
     330             :         }
     331             :     }
     332          40 :     for (const auto &forbidden : {"m", "exp_m", "std_m"}) {
     333          30 :         if (columns.contains(forbidden)) {
     334           0 :             throw std::runtime_error(
     335           0 :                 fmt::format("The table of states of species '{}' must not contain a column '{}', "
     336             :                             "since the quantum number m is generated by the software itself.",
     337           0 :                             species, forbidden));
     338             :         }
     339             :     }
     340             : 
     341             :     // If another thread inserted the same entry concurrently, insert keeps the existing value and
     342             :     // returns an iterator to it, so the reference stays valid (elements are never erased).
     343          10 :     return column_info_cache.insert({table_path, std::move(columns)}).first->second;
     344        2187 : }
     345             : 
     346         520 : std::shared_ptr<const KetAtom> Database::get_ket(const std::string &species,
     347             :                                                  const AtomDescriptionByParameters &description) {
     348             :     // Check that the specifications are valid
     349         520 :     if (!description.quantum_numbers.contains("m")) {
     350           0 :         throw std::invalid_argument("The quantum number m must be specified.");
     351             :     }
     352        2581 :     for (const auto &[name, value] : description.quantum_numbers) {
     353        2061 :         if ((name == "f" || name == "m") && 2 * value != std::rint(2 * value)) {
     354           0 :             throw std::invalid_argument("The quantum number " + name +
     355           0 :                                         " must be an integer or half-integer.");
     356             :         }
     357        2061 :         if (name == "f" && value < 0) {
     358           0 :             throw std::invalid_argument("The quantum number " + name + " must be positive.");
     359             :         }
     360             :     }
     361             : 
     362         520 :     const auto &columns = get_states_column_info(species);
     363             : 
     364             :     // Describe the state. Integer-valued quantum numbers and (f, m, parity) are matched exactly,
     365             :     // while all other quantum numbers are matched within a +-0.5 window (they can deviate from the
     366             :     // requested value, e.g. expectation values in MQDT). The result is ordered by the distance to
     367             :     // the requested values, so the nearest state is returned.
     368         520 :     std::string where;
     369         520 :     std::string where_separator;
     370         520 :     std::string orderby;
     371         520 :     std::string orderby_separator;
     372         520 :     if (description.energy.has_value()) {
     373             :         // The following condition derives from demanding that quantum number n that corresponds to
     374             :         // the energy "E_n = -1/(2*n^2)" is not off by more than 1 from the actual quantum number n,
     375             :         // i.e., "sqrt(-1/(2*E_n)) - sqrt(-1/(2*E_{n-1})) = 1"
     376           0 :         double n_from_energy = std::sqrt(-1 / (2 * description.energy.value()));
     377           0 :         where += where_separator +
     378           0 :             fmt::format("SQRT(-1/(2*energy)) BETWEEN {} AND {}", n_from_energy - 0.5,
     379           0 :                         n_from_energy + 0.5);
     380           0 :         where_separator = " AND ";
     381           0 :         orderby += orderby_separator + fmt::format("(SQRT(-1/(2*energy)) - {})^2", n_from_energy);
     382           0 :         orderby_separator = " + ";
     383             :     }
     384        2579 :     for (const auto &[name, value] : description.quantum_numbers) {
     385        2060 :         if (name == "m") {
     386         520 :             continue; // m is encoded into the id, not stored as a queryable column
     387             :         }
     388        1540 :         std::string column = columns.contains("exp_" + name) ? "exp_" + name : name;
     389        1540 :         if (!columns.contains(column)) {
     390           2 :             throw std::invalid_argument(
     391           2 :                 fmt::format("The quantum number '{}' is not stored in the database table for "
     392             :                             "species '{}'.",
     393           2 :                             name, species));
     394             :         }
     395        1539 :         double tolerance = must_match_exactly(columns, column) ? 0.0 : 0.5;
     396        1539 :         where += where_separator +
     397        4617 :             fmt::format("{} BETWEEN {} AND {}", column, value - tolerance, value + tolerance);
     398        1539 :         where_separator = " AND ";
     399        3078 :         orderby += orderby_separator + fmt::format("({} - {})^2", column, value);
     400        1539 :         orderby_separator = " + ";
     401        1540 :     }
     402         519 :     if (where_separator.empty()) {
     403           0 :         where += "FALSE";
     404             :     }
     405         519 :     if (orderby_separator.empty()) {
     406           0 :         orderby += "id";
     407             :     }
     408             : 
     409             :     // Ask the database for the described state
     410         519 :     set_task_status("Loading atomic ket from database...");
     411        1038 :     auto result = con->Query(fmt::format(
     412             :         R"(SELECT *, {} AS order_val FROM '{}' WHERE {} ORDER BY order_val ASC LIMIT 2)", orderby,
     413        1557 :         manager->get_path(species, "states"), where));
     414             : 
     415         519 :     if (result->HasError()) {
     416           0 :         throw cpptrace::runtime_error("Error querying the database: " + result->GetError());
     417             :     }
     418             : 
     419         519 :     if (result->RowCount() == 0) {
     420           5 :         throw std::invalid_argument("No state found.");
     421             :     }
     422             : 
     423             :     // Get the first chunk of the results (the first chunk is sufficient as we need two rows at
     424             :     // most). Every column except energy, id and the synthetic order_val is treated as a quantum
     425             :     // number; m is not a database column and is injected from the description.
     426         514 :     const auto &types = result->types;
     427         514 :     const auto &names = result->names;
     428        2056 :     const std::unordered_set<std::string> excluded_columns = {"energy", "id", "order_val"};
     429         514 :     auto chunk = result->Fetch();
     430             : 
     431         514 :     size_t energy_column = get_column_index(names, "energy");
     432         514 :     size_t id_column = get_column_index(names, "id");
     433         514 :     double quantum_number_m = description.quantum_numbers.at("m");
     434             : 
     435         515 :     auto make_ket = [&](size_t row) {
     436             :         auto quantum_numbers =
     437         515 :             get_quantum_numbers_from_row(*chunk, types, names, excluded_columns, row);
     438         515 :         quantum_numbers.values["m"] = quantum_number_m;
     439         515 :         double energy = get_entry_as_double(chunk->data[energy_column], types[energy_column], row);
     440        1030 :         auto id = utils::encode_as_ket_id({.id = static_cast<size_t>(get_entry_as_int64(
     441         515 :                                                chunk->data[id_column], types[id_column], row)),
     442         515 :                                            .m = quantum_number_m});
     443         515 :         return KetAtom(typename KetAtom::Private(), energy, species,
     444        1030 :                        std::move(quantum_numbers.values), std::move(quantum_numbers.stds), *this,
     445        1545 :                        id);
     446         515 :     };
     447             : 
     448             :     // Check that the ket is uniquely specified. If not, throw a KetNotUniqueError.
     449         514 :     if (chunk->size() > 1) {
     450           7 :         size_t order_val_column = get_column_index(names, "order_val");
     451             :         auto order_val_0 =
     452           7 :             get_entry_as_double(chunk->data[order_val_column], types[order_val_column], 0);
     453             :         auto order_val_1 =
     454           7 :             get_entry_as_double(chunk->data[order_val_column], types[order_val_column], 1);
     455             : 
     456           7 :         if (order_val_1 - order_val_0 <= order_val_0) {
     457           4 :             throw KetNotUniqueError({
     458           0 :                 std::make_shared<const KetAtom>(make_ket(0)),
     459           1 :                 std::make_shared<const KetAtom>(make_ket(1)),
     460           5 :             });
     461             :         }
     462             :     }
     463             : 
     464             :     // Construct the state
     465         513 :     auto ket = make_ket(0);
     466             : 
     467             :     // Check database consistency
     468         515 :     ensure_consistent_quantum_numbers(ket.get_quantum_number("f"), quantum_number_m);
     469             : 
     470        1024 :     return std::make_shared<const KetAtom>(std::move(ket));
     471        1585 : }
     472             : 
     473             : template <typename Scalar>
     474             : std::shared_ptr<const BasisAtom<Scalar>>
     475        1667 : Database::get_basis(const std::string &species, const AtomDescriptionByRanges &description,
     476             :                     const std::vector<size_t> &additional_ket_ids) {
     477             :     // The quantum number m is restricted separately because it is generated by UNNEST below.
     478        3334 :     auto range_quantum_number_m = [&description]() {
     479        1667 :         auto it = description.quantum_number_ranges.find("m");
     480        3334 :         return it != description.quantum_number_ranges.end() ? it->second : Range<double>{};
     481        1667 :     }();
     482             : 
     483        1667 :     const auto &columns = get_states_column_info(species);
     484             : 
     485             :     // Describe the states by all restrictions that do not involve the quantum number m
     486        1667 :     std::string where = "(";
     487        1667 :     std::string separator;
     488        1667 :     if (description.range_energy.is_finite()) {
     489          57 :         where += separator +
     490          57 :             fmt::format("energy BETWEEN {} AND {}", description.range_energy.min(),
     491          57 :                         description.range_energy.max());
     492          57 :         separator = " AND ";
     493             :     }
     494        3446 :     for (const auto &[name, range] : description.quantum_number_ranges) {
     495         949 :         if (name == "m" || !range.is_finite()) {
     496         119 :             continue;
     497             :         }
     498         830 :         std::string exp_column = "exp_" + name;
     499         830 :         std::string std_column = "std_" + name;
     500         830 :         if (columns.contains(exp_column) && columns.contains(std_column)) {
     501         406 :             where += separator +
     502             :                 format_expectation_value_range(
     503             :                          exp_column, std_column, range,
     504         406 :                          description.quantum_number_standard_deviation_factor);
     505         424 :         } else if (columns.contains(name)) {
     506         423 :             where += separator + format_quantum_number_range(name, range);
     507             :         } else {
     508           2 :             throw std::invalid_argument(
     509             :                 fmt::format("The quantum number '{}' is not stored in the database table for "
     510             :                             "species '{}'.",
     511             :                             name, species));
     512             :         }
     513         829 :         separator = " AND ";
     514             :     }
     515        1666 :     if (separator.empty()) {
     516             :         // If the description contains no restrictions at all, it describes no states
     517        1247 :         where += range_quantum_number_m.is_finite() ? "TRUE" : "FALSE";
     518             :     }
     519        1666 :     where += ")";
     520             : 
     521             :     // Describe the restriction of the quantum number m
     522        1666 :     std::string where_m = "(";
     523        1666 :     if (range_quantum_number_m.is_finite()) {
     524         119 :         where_m += format_quantum_number_range("m", range_quantum_number_m);
     525             :     } else {
     526        1547 :         where_m += "TRUE";
     527             :     }
     528        1666 :     where_m += ")";
     529             : 
     530             :     // Create a table containing the described states
     531        1666 :     std::string canonical_basis_id;
     532             :     {
     533        1666 :         auto result = con->Query(R"(SELECT UUID()::varchar)");
     534        1666 :         if (result->HasError()) {
     535           0 :             throw cpptrace::runtime_error("Error selecting canonical_basis_id: " +
     536           0 :                                           result->GetError());
     537             :         }
     538        1666 :         canonical_basis_id =
     539        3332 :             duckdb::FlatVector::GetData<duckdb::string_t>(result->Fetch()->data[0])[0].GetString();
     540        1666 :     }
     541             :     {
     542        1666 :         set_task_status("Selecting atomic basis states...");
     543        4998 :         auto result = con->Query(fmt::format(
     544             :             R"(CREATE TEMP TABLE '{}' AS SELECT *, id*{}+(2*m+{})::bigint AS ketid FROM (
     545             :                 SELECT *,
     546             :                 UNNEST(list_transform(generate_series(0,(2*f)::bigint),
     547             :                 x -> x::double-f)) AS m FROM (
     548             :                     SELECT * FROM '{}' WHERE {}
     549             :                 )
     550             :             ) WHERE {})",
     551             :             canonical_basis_id, utils::KET_ID_STRIDE, utils::M_OFFSET,
     552        3332 :             manager->get_path(species, "states"), where, where_m));
     553             : 
     554        1666 :         if (result->HasError()) {
     555           0 :             throw cpptrace::runtime_error("Error creating table: " + result->GetError());
     556             :         }
     557        1666 :     }
     558             :     {
     559        3332 :         auto result = con->Query(
     560             :             fmt::format(R"(ALTER TABLE '{}' ADD PRIMARY KEY (ketid))", canonical_basis_id));
     561             : 
     562        1666 :         if (result->HasError()) {
     563           0 :             throw cpptrace::runtime_error("Error adding primary key: " + result->GetError());
     564             :         }
     565        1666 :     }
     566             : 
     567             :     // Add the additional kets to the table if they are not already contained in it
     568        1666 :     if (!additional_ket_ids.empty()) {
     569        1321 :         set_task_status("Selecting additional kets...");
     570        1321 :         std::vector<std::string> values;
     571        1321 :         values.reserve(additional_ket_ids.size());
     572        3527 :         for (size_t ket_id : additional_ket_ids) {
     573        2206 :             auto [id, m] = utils::decode_from_ket_id(ket_id);
     574        4412 :             values.push_back(fmt::format("({}, {}::double, {})", id, m, ket_id));
     575             :         }
     576        5284 :         auto result = con->Query(fmt::format(
     577             :             R"(INSERT OR IGNORE INTO '{0}'
     578             :             SELECT s.*, v.m, v.ketid FROM '{1}' AS s
     579             :             JOIN (VALUES {2}) AS v(id, m, ketid) ON s.id = v.id)",
     580        3963 :             canonical_basis_id, manager->get_path(species, "states"), fmt::join(values, ",")));
     581             : 
     582        1321 :         if (result->HasError()) {
     583           0 :             throw cpptrace::runtime_error("Error adding additional kets: " + result->GetError());
     584             :         }
     585        1321 :     }
     586             : 
     587             :     // Ask the table for the extreme values of the quantum numbers
     588             :     {
     589        1666 :         set_task_status("Validating atomic basis coverage...");
     590             : 
     591             :         // Collect the finite quantum-number ranges to validate against the loaded basis. Energy is
     592             :         // handled separately because it is compared via its corresponding effective quantum number.
     593             :         struct CoverageCheck {
     594             :             std::string name;
     595             :             std::string column;
     596             :             double tolerance{}; // allowed slack between the requested and the available range
     597             :             Range<double> range;
     598             :         };
     599        1666 :         std::vector<CoverageCheck> checks;
     600        2614 :         for (const auto &[name, range] : description.quantum_number_ranges) {
     601         948 :             if (!range.is_finite()) {
     602           0 :                 continue;
     603             :             }
     604         948 :             bool is_expectation_value =
     605         948 :                 columns.contains("exp_" + name) && columns.contains("std_" + name);
     606         948 :             std::string column = is_expectation_value ? "exp_" + name : name;
     607         948 :             double tolerance = must_match_exactly(columns, column) ? 0.0 : 1.0;
     608         948 :             checks.push_back({name, column, tolerance, range});
     609             :         }
     610             : 
     611        1666 :         std::string select;
     612        1666 :         std::string separator;
     613        1666 :         if (description.range_energy.is_finite()) {
     614          57 :             select += separator + "MIN(energy) AS min_energy, MAX(energy) AS max_energy";
     615          57 :             separator = ", ";
     616             :         }
     617        2614 :         for (const auto &check : checks) {
     618         948 :             select += separator +
     619         948 :                 fmt::format("MIN({0}) AS min_{1}, MAX({0}) AS max_{1}", check.column, check.name);
     620         948 :             separator = ", ";
     621             :         }
     622             : 
     623        1666 :         if (!separator.empty()) {
     624         838 :             auto result =
     625             :                 con->Query(fmt::format(R"(SELECT {} FROM '{}')", select, canonical_basis_id));
     626             : 
     627         419 :             if (result->HasError()) {
     628           0 :                 throw cpptrace::runtime_error("Error querying the database: " + result->GetError());
     629             :             }
     630             : 
     631         419 :             auto chunk = result->Fetch();
     632         419 :             const auto &types = result->types;
     633             : 
     634        2429 :             for (size_t i = 0; i < chunk->ColumnCount(); i++) {
     635        2010 :                 if (duckdb::FlatVector::IsNull(chunk->data[i], 0)) {
     636           0 :                     throw std::invalid_argument("No state found.");
     637             :                 }
     638             :             }
     639             : 
     640         419 :             size_t idx = 0;
     641         419 :             if (description.range_energy.is_finite()) {
     642          57 :                 auto min_energy = get_entry_as_double(chunk->data[idx], types[idx], 0);
     643          57 :                 idx++;
     644         114 :                 if (std::sqrt(-1 / (2 * min_energy)) - 1 >
     645          57 :                     std::sqrt(-1 / (2 * description.range_energy.min()))) {
     646           0 :                     SPDLOG_DEBUG("No state found with the requested minimum energy. Requested: {}, "
     647             :                                  "found: {}.",
     648             :                                  description.range_energy.min(), min_energy);
     649             :                 }
     650          57 :                 auto max_energy = get_entry_as_double(chunk->data[idx], types[idx], 0);
     651          57 :                 idx++;
     652         114 :                 if (std::sqrt(-1 / (2 * max_energy)) + 1 <
     653          57 :                     std::sqrt(-1 / (2 * description.range_energy.max()))) {
     654           0 :                     SPDLOG_DEBUG("No state found with the requested maximum energy. Requested: {}, "
     655             :                                  "found: {}.",
     656             :                                  description.range_energy.max(), max_energy);
     657             :                 }
     658             :             }
     659        1367 :             for (const auto &check : checks) {
     660         948 :                 auto min_value = get_entry_as_double(chunk->data[idx], types[idx], 0);
     661         948 :                 idx++;
     662         948 :                 if (min_value - check.tolerance > check.range.min()) {
     663         144 :                     SPDLOG_DEBUG("No state found with the requested minimum quantum number {}. "
     664             :                                  "Requested: {}, found: {}.",
     665             :                                  check.name, check.range.min(), min_value);
     666             :                 }
     667         948 :                 auto max_value = get_entry_as_double(chunk->data[idx], types[idx], 0);
     668         948 :                 idx++;
     669         948 :                 if (max_value + check.tolerance < check.range.max()) {
     670          52 :                     SPDLOG_DEBUG("No state found with the requested maximum quantum number {}. "
     671             :                                  "Requested: {}, found: {}.",
     672             :                                  check.name, check.range.max(), max_value);
     673             :                 }
     674             :             }
     675         419 :         }
     676        1666 :     }
     677             : 
     678             :     // Ask the table for the described states
     679        1666 :     set_task_status("Loading atomic basis states...");
     680        3332 :     auto result =
     681             :         con->Query(fmt::format(R"(SELECT * FROM '{}' ORDER BY ketid ASC)", canonical_basis_id));
     682             : 
     683        1666 :     if (result->HasError()) {
     684           0 :         throw cpptrace::runtime_error("Error querying the database: " + result->GetError());
     685             :     }
     686             : 
     687        1666 :     if (result->RowCount() == 0) {
     688           0 :         throw std::invalid_argument("No state found.");
     689             :     }
     690             : 
     691             :     // Construct the states. Every column except energy and the raw/encoded id is treated as a
     692             :     // quantum number ("id" is the raw states-table id, "ketid" the m-encoded id of the basis
     693             :     // state).
     694        1666 :     const auto &types = result->types;
     695        1666 :     const auto &names = result->names;
     696        6664 :     const std::unordered_set<std::string> excluded_columns = {"energy", "id", "ketid"};
     697        1666 :     size_t energy_column = get_column_index(names, "energy");
     698        1666 :     size_t ketid_column = get_column_index(names, "ketid");
     699             : 
     700        1666 :     std::vector<std::shared_ptr<const KetAtom>> kets;
     701        1666 :     kets.reserve(result->RowCount());
     702        1666 :     double last_energy = std::numeric_limits<double>::lowest();
     703             : 
     704        3332 :     for (auto chunk = result->Fetch(); chunk; chunk = result->Fetch()) {
     705        1666 :         set_task_status("Constructing atomic basis...");
     706             : 
     707       48993 :         for (size_t i = 0; i < chunk->size(); i++) {
     708       94654 :             auto quantum_numbers =
     709       47327 :                 get_quantum_numbers_from_row(*chunk, types, names, excluded_columns, i);
     710       47327 :             double energy =
     711       47327 :                 get_entry_as_double(chunk->data[energy_column], types[energy_column], i);
     712       47327 :             auto id = static_cast<size_t>(
     713       47327 :                 get_entry_as_int64(chunk->data[ketid_column], types[ketid_column], i));
     714             : 
     715             :             // Check database consistency
     716       47327 :             ensure_consistent_quantum_numbers(quantum_numbers.values.at("f"),
     717       94654 :                                               quantum_numbers.values.at("m"));
     718       47327 :             if (energy < last_energy) {
     719           0 :                 throw std::runtime_error("The states are not sorted by energy.");
     720             :             }
     721       47327 :             last_energy = energy;
     722             : 
     723             :             // Append a new state
     724       47327 :             kets.push_back(std::make_shared<const KetAtom>(
     725       94654 :                 typename KetAtom::Private(), energy, species, std::move(quantum_numbers.values),
     726       47327 :                 std::move(quantum_numbers.stds), *this, id));
     727             :         }
     728             :     }
     729             : 
     730           0 :     return std::make_shared<const BasisAtom<Scalar>>(typename BasisAtom<Scalar>::Private(),
     731        1666 :                                                      std::move(kets), std::move(canonical_basis_id),
     732        3332 :                                                      *this);
     733        5948 : }
     734             : 
     735             : template <typename Scalar>
     736        5888 : Eigen::SparseMatrix<Scalar, Eigen::RowMajor> Database::get_matrix_elements_in_canonical_basis(
     737             :     std::shared_ptr<const BasisAtom<Scalar>> initial_basis,
     738             :     std::shared_ptr<const BasisAtom<Scalar>> final_basis, OperatorType type, int q) {
     739             :     using real_t = typename traits::NumTraits<Scalar>::real_t;
     740             :     using cached_matrix_ptr_t = std::shared_ptr<const cached_matrix_t>;
     741             : 
     742        5888 :     if (&initial_basis->get_database() != this || &final_basis->get_database() != this) {
     743           1 :         throw std::invalid_argument(
     744             :             "The initial and final bases must belong to the Database instance used for the "
     745             :             "matrix element calculation.");
     746             :     }
     747        5916 :     if (initial_basis->get_species() != final_basis->get_species()) {
     748           0 :         throw std::invalid_argument(fmt::format(
     749             :             "The initial and final bases must have the same species, but got '{}' and '{}'.",
     750           0 :             initial_basis->get_species(), final_basis->get_species()));
     751             :     }
     752             : 
     753        5909 :     std::string specifier;
     754        5971 :     int kappa{};
     755        5971 :     switch (type) {
     756         129 :     case OperatorType::ELECTRIC_MONOPOLE:
     757         129 :         specifier = "matrix_elements_m";
     758         129 :         kappa = 0;
     759         129 :         break;
     760        4023 :     case OperatorType::ELECTRIC_DIPOLE:
     761        4023 :         specifier = "matrix_elements_d";
     762        4062 :         kappa = 1;
     763        4062 :         break;
     764         499 :     case OperatorType::ELECTRIC_QUADRUPOLE:
     765         499 :         specifier = "matrix_elements_q";
     766         499 :         kappa = 2;
     767         499 :         break;
     768          95 :     case OperatorType::ELECTRIC_QUADRUPOLE_ZERO:
     769          95 :         specifier = "matrix_elements_q0";
     770          94 :         kappa = 0;
     771          94 :         break;
     772           0 :     case OperatorType::ELECTRIC_OCTUPOLE:
     773           0 :         specifier = "matrix_elements_o";
     774           0 :         kappa = 3;
     775           0 :         break;
     776          56 :     case OperatorType::MAGNETIC_DIPOLE:
     777          56 :         specifier = "matrix_elements_mu";
     778          57 :         kappa = 1;
     779          57 :         break;
     780           4 :     case OperatorType::ENERGY:
     781           4 :         specifier = "energy";
     782           4 :         kappa = 0;
     783           4 :         break;
     784        1137 :     case OperatorType::IDENTITY:
     785        1137 :         specifier = "identity";
     786        1137 :         kappa = 0;
     787        1137 :         break;
     788          28 :     default:
     789          28 :         throw std::invalid_argument("Unknown operator type.");
     790             :     }
     791             : 
     792        5982 :     std::string canonical_basis_id_initial = initial_basis->get_canonical_basis_id();
     793        5952 :     std::string canonical_basis_id_final = final_basis->get_canonical_basis_id();
     794        5968 :     std::string cache_key = fmt::format("{}_{}_{}_{}", specifier, q, canonical_basis_id_initial,
     795             :                                         canonical_basis_id_final);
     796        5973 :     auto &matrix_elements_cache = get_matrix_elements_cache();
     797        5885 :     std::promise<cached_matrix_ptr_t> matrix_promise;
     798        5907 :     auto [cache_it, inserted] =
     799             :         matrix_elements_cache.insert({cache_key, matrix_promise.get_future().share()});
     800             : 
     801        5884 :     if (inserted) {
     802             :         try {
     803        1775 :             Eigen::Index num_rows = final_basis->get_number_of_kets();
     804        1774 :             Eigen::Index num_cols = initial_basis->get_number_of_kets();
     805             : 
     806        1774 :             std::vector<int> outerIndexPtr;
     807        1774 :             std::vector<int> innerIndices;
     808        1774 :             std::vector<real_t> values;
     809             : 
     810             :             // Check that the specifications are valid
     811        1774 :             if (std::abs(q) > kappa) {
     812           0 :                 throw std::invalid_argument("Invalid q.");
     813             :             }
     814             : 
     815             :             // Ask the database for the operator
     816        1774 :             set_task_status("Loading matrix elements from database...");
     817        1777 :             std::string species = initial_basis->get_species();
     818        1777 :             duckdb::unique_ptr<duckdb::MaterializedQueryResult> result;
     819        1777 :             if (specifier == "identity") {
     820        1972 :                 result = con->Query(fmt::format(
     821             :                     R"(SELECT s2.ketid AS row, s1.ketid AS col, 1.0::DOUBLE AS val
     822             :                     FROM '{}' AS s1
     823             :                     INNER JOIN '{}' AS s2 ON s1.ketid = s2.ketid
     824             :                     ORDER BY row ASC)",
     825             :                     canonical_basis_id_initial, canonical_basis_id_final));
     826         790 :             } else if (specifier == "energy") {
     827           6 :                 result = con->Query(fmt::format(
     828             :                     R"(SELECT s2.ketid AS row, s1.ketid AS col, s1.energy AS val
     829             :                     FROM '{}' AS s1
     830             :                     INNER JOIN '{}' AS s2 ON s1.ketid = s2.ketid
     831             :                     ORDER BY row ASC)",
     832             :                     canonical_basis_id_initial, canonical_basis_id_final));
     833             :             } else {
     834        3151 :                 result = con->Query(fmt::format(
     835             :                     R"(WITH s1 AS (
     836             :                         SELECT id, f, m, ketid FROM '{}'
     837             :                     ),
     838             :                     s2 AS (
     839             :                         SELECT id, f, m, ketid FROM '{}'
     840             :                     ),
     841             :                     b AS (
     842             :                         SELECT MIN(f) AS min_f, MAX(f) AS max_f,
     843             :                         MIN(id) AS min_id, MAX(id) AS max_id
     844             :                         FROM (SELECT f, id FROM s1 UNION ALL SELECT f, id FROM s2)
     845             :                     ),
     846             :                     w_filtered AS (
     847             :                         SELECT *
     848             :                         FROM '{}'
     849             :                         WHERE kappa = {} AND q = {} AND
     850             :                         f_initial BETWEEN (SELECT min_f FROM b) AND (SELECT max_f FROM b) AND
     851             :                         f_final BETWEEN (SELECT min_f FROM b) AND (SELECT max_f FROM b)
     852             :                     ),
     853             :                     e_filtered AS (
     854             :                         SELECT *
     855             :                         FROM '{}'
     856             :                         WHERE
     857             :                         id_initial BETWEEN (SELECT min_id FROM b) AND (SELECT max_id FROM b) AND
     858             :                         id_final BETWEEN (SELECT min_id FROM b) AND (SELECT max_id FROM b)
     859             :                     )
     860             :                     SELECT
     861             :                     s2.ketid AS row,
     862             :                     s1.ketid AS col,
     863             :                     e.val*w.val AS val
     864             :                     FROM e_filtered AS e
     865             :                     JOIN s1 ON e.id_initial = s1.id
     866             :                     JOIN s2 ON e.id_final = s2.id
     867             :                     JOIN w_filtered AS w ON
     868             :                     w.f_initial = s1.f AND w.m_initial = s1.m AND
     869             :                     w.f_final = s2.f AND w.m_final = s2.m
     870             :                     ORDER BY row ASC, col ASC)",
     871             :                     canonical_basis_id_initial, canonical_basis_id_final,
     872        1575 :                     manager->get_path("misc", "wigner"), kappa, q,
     873             :                     manager->get_path(species, specifier)));
     874             :             }
     875             : 
     876        1777 :             if (result->HasError()) {
     877           0 :                 throw cpptrace::runtime_error("Error querying the database: " + result->GetError());
     878             :             }
     879             : 
     880             :             // Check the types of the columns
     881        1777 :             const auto &types = result->types;
     882        1777 :             const auto &labels = result->names;
     883        7108 :             const std::vector<duckdb::LogicalType> ref_types = {duckdb::LogicalType::BIGINT,
     884             :                                                                 duckdb::LogicalType::BIGINT,
     885             :                                                                 duckdb::LogicalType::DOUBLE};
     886        7108 :             for (size_t i = 0; i < types.size(); i++) {
     887        5331 :                 if (types[i] != ref_types[i]) {
     888           0 :                     throw std::runtime_error("Wrong type for '" + labels[i] + "'.");
     889             :                 }
     890             :             }
     891             : 
     892        1777 :             set_task_status("Constructing matrix elements...");
     893             : 
     894             :             // Construct the matrix
     895        1777 :             int num_entries = static_cast<int>(result->RowCount());
     896        1777 :             outerIndexPtr.reserve(num_rows + 1);
     897        1777 :             innerIndices.reserve(num_entries);
     898        1777 :             values.reserve(num_entries);
     899             : 
     900        1777 :             int last_row = -1;
     901             : 
     902        3723 :             for (auto chunk = result->Fetch(); chunk; chunk = result->Fetch()) {
     903        1946 :                 auto *chunk_row = duckdb::FlatVector::GetData<int64_t>(chunk->data[0]);
     904        1946 :                 auto *chunk_col = duckdb::FlatVector::GetData<int64_t>(chunk->data[1]);
     905        1946 :                 auto *chunk_val = duckdb::FlatVector::GetData<double>(chunk->data[2]);
     906             : 
     907      878163 :                 for (size_t i = 0; i < chunk->size(); i++) {
     908      876217 :                     int row = final_basis->get_ket_index_from_id(chunk_row[i]);
     909      876217 :                     if (row != last_row) {
     910       58970 :                         if (row < last_row) {
     911           0 :                             throw std::runtime_error("The rows are not sorted.");
     912             :                         }
     913      124143 :                         for (; last_row < row; last_row++) {
     914       65173 :                             outerIndexPtr.push_back(static_cast<int>(innerIndices.size()));
     915             :                         }
     916             :                     }
     917      876217 :                     innerIndices.push_back(initial_basis->get_ket_index_from_id(chunk_col[i]));
     918      876217 :                     values.push_back(chunk_val[i]);
     919             :                 }
     920             :             }
     921             : 
     922        6685 :             for (; last_row < num_rows + 1; last_row++) {
     923        4908 :                 outerIndexPtr.push_back(static_cast<int>(innerIndices.size()));
     924             :             }
     925             : 
     926        1777 :             Eigen::Map<const cached_matrix_t> matrix_map(num_rows, num_cols, values.size(),
     927        1777 :                                                          outerIndexPtr.data(), innerIndices.data(),
     928        1777 :                                                          values.data());
     929             : 
     930        1777 :             auto cached_matrix = std::make_shared<const cached_matrix_t>(matrix_map);
     931        1777 :             matrix_promise.set_value(std::move(cached_matrix));
     932        1777 :         } catch (...) {
     933           0 :             matrix_promise.set_exception(std::current_exception());
     934           0 :             throw;
     935             :         }
     936             :     }
     937             : 
     938        5886 :     set_task_status("Returning matrix elements in canonical basis...");
     939             : 
     940       11995 :     return cache_it->second.get()->template cast<Scalar>();
     941        7747 : }
     942             : 
     943        1266 : int Database::get_minimal_kappa(const std::string &species) {
     944             :     // The monopole operator is only stored for charged species. Thus, if the table exists, the
     945             :     // species is charged and the monopole operator has to be taken into account.
     946        1266 :     return manager->has_table(species, "matrix_elements_m") ? 0 : 1;
     947             : }
     948             : 
     949           4 : bool Database::get_download_missing() const { return download_missing_; }
     950             : 
     951           2 : bool Database::get_use_cache() const { return use_cache_; }
     952             : 
     953          10 : std::filesystem::path Database::get_database_dir() const { return database_dir_; }
     954             : 
     955           0 : std::string Database::get_versions_info() const { return manager->get_versions_info(); }
     956             : 
     957        5985 : Database::matrix_elements_cache_t &Database::get_matrix_elements_cache() {
     958        5985 :     static matrix_elements_cache_t matrix_elements_cache;
     959        5946 :     return matrix_elements_cache;
     960             : }
     961             : 
     962          63 : Database &Database::get_global_instance() {
     963         126 :     return get_global_instance_without_checks(default_download_missing, default_use_cache,
     964         126 :                                               default_database_dir);
     965             : }
     966             : 
     967           0 : Database &Database::get_global_instance(bool download_missing) {
     968           0 :     Database &database = get_global_instance_without_checks(download_missing, default_use_cache,
     969             :                                                             default_database_dir);
     970           0 :     if (download_missing != database.download_missing_) {
     971           0 :         throw std::invalid_argument(
     972           0 :             "The 'download_missing' argument must not change between calls to the method.");
     973             :     }
     974           0 :     return database;
     975             : }
     976             : 
     977           0 : Database &Database::get_global_instance(std::filesystem::path database_dir) {
     978           0 :     if (database_dir.empty()) {
     979           0 :         database_dir = default_database_dir;
     980             :     }
     981           0 :     Database &database = get_global_instance_without_checks(default_download_missing,
     982             :                                                             default_use_cache, database_dir);
     983           0 :     if (!std::filesystem::exists(database_dir) ||
     984           0 :         std::filesystem::canonical(database_dir) != database.database_dir_) {
     985           0 :         throw std::invalid_argument(
     986           0 :             "The 'database_dir' argument must not change between calls to the method.");
     987             :     }
     988           0 :     return database;
     989             : }
     990             : 
     991           1 : Database &Database::get_global_instance(bool download_missing, bool use_cache,
     992             :                                         std::filesystem::path database_dir) {
     993           1 :     if (database_dir.empty()) {
     994           0 :         database_dir = default_database_dir;
     995             :     }
     996             :     Database &database =
     997           1 :         get_global_instance_without_checks(download_missing, use_cache, database_dir);
     998           1 :     if (download_missing != database.download_missing_ || use_cache != database.use_cache_ ||
     999           3 :         !std::filesystem::exists(database_dir) ||
    1000           2 :         std::filesystem::canonical(database_dir) != database.database_dir_) {
    1001           0 :         throw std::invalid_argument(
    1002             :             "The 'download_missing', 'use_cache' and 'database_dir' arguments must not "
    1003           0 :             "change between calls to the method.");
    1004             :     }
    1005           1 :     return database;
    1006             : }
    1007             : 
    1008          64 : Database &Database::get_global_instance_without_checks(bool download_missing, bool use_cache,
    1009             :                                                        std::filesystem::path database_dir) {
    1010          64 :     static Database database(download_missing, use_cache, std::move(database_dir));
    1011          64 :     return database;
    1012             : }
    1013             : 
    1014             : struct database_dir_noexcept : std::filesystem::path {
    1015           2 :     explicit database_dir_noexcept() noexcept try : std
    1016           2 :         ::filesystem::path(paths::get_cache_directory() / "database") {}
    1017           0 :     catch (...) {
    1018           0 :         SPDLOG_ERROR("Error getting the PairInteraction cache directory.");
    1019           0 :         std::terminate();
    1020           2 :     }
    1021             : };
    1022             : 
    1023             : const std::filesystem::path Database::default_database_dir = database_dir_noexcept();
    1024             : 
    1025             : // Explicit instantiations
    1026             : // NOLINTBEGIN(bugprone-macro-parentheses, cppcoreguidelines-macro-usage)
    1027             : #define INSTANTIATE_GETTERS(SCALAR)                                                                \
    1028             :     template std::shared_ptr<const BasisAtom<SCALAR>> Database::get_basis<SCALAR>(                 \
    1029             :         const std::string &species, const AtomDescriptionByRanges &description,                    \
    1030             :         const std::vector<size_t> &additional_ket_ids);                                            \
    1031             :     template Eigen::SparseMatrix<SCALAR, Eigen::RowMajor>                                          \
    1032             :     Database::get_matrix_elements_in_canonical_basis<SCALAR>(                                      \
    1033             :         std::shared_ptr<const BasisAtom<SCALAR>> initial_basis,                                    \
    1034             :         std::shared_ptr<const BasisAtom<SCALAR>> final_basis, OperatorType type, int q);
    1035             : // NOLINTEND(bugprone-macro-parentheses, cppcoreguidelines-macro-usage)
    1036             : 
    1037             : INSTANTIATE_GETTERS(double)
    1038             : INSTANTIATE_GETTERS(std::complex<double>)
    1039             : 
    1040             : #undef INSTANTIATE_GETTERS
    1041             : } // namespace pairinteraction

Generated by: LCOV version 1.16