phigrape/config.cpp

#include "config.h"
#include <stdexcept>
#include <cstdio>
#include <fstream>
#include <algorithm>
#include <limits>

// Would be a bit more elegant to do the whole thing with std::variant.

using Dictionary = std::unordered_map<std::string,std::string>;

static constexpr double nix = -std::numeric_limits<double>::max(); // avoid nans

std::string Config::strip(const std::string str)
{
    std::string str_new = str;
    auto pos = str_new.find_first_not_of(" \t");
    if (pos != std::string::npos) str_new = str_new.substr(pos, str_new.size());
    pos = str_new.find_last_not_of(" \t");
    if (pos != std::string::npos) str_new = str_new.substr(0, pos+1);
    return str_new;
}

Dictionary Config::read_config_file(const std::string file_name)
{
    std::unordered_map<std::string,std::string> dictionary;
    std::ifstream file(file_name);
    if (!file.good()) throw std::runtime_error("File not found.");
    std::string str;
    int line_number = 0;
    while (std::getline(file, str)) {
        line_number++;
        auto pos = str.find('#');
        if (pos != std::string::npos) str = str.substr(0, pos);
        str = strip(str);
        if (str.size() == 0) continue;
        pos = str.find_first_of("=");
        if (pos == std::string::npos) throw std::runtime_error("Error: expected a key-value pair in line " + std::to_string(line_number) + " of file " + file_name);
        std::string key = strip(str.substr(0, pos));
        pos = str.find_first_not_of(" \t", pos+1);
        std::string val = strip(str.substr(pos, str.size()));
        dictionary[key] = val;
    }
    return dictionary;
}

template<> std::string Config::string_cast(const std::string str)
{
    return str;
}

template<> double Config::string_cast(const std::string str)
{
    size_t idx;
    auto value = std::stod(str, &idx);
    if (idx == str.size()) return value;
    else throw std::runtime_error("Cannot convert \"" + str + "\" into a double");
}

template<> int Config::string_cast(const std::string str)
{
    size_t idx;
    auto value = std::stoi(str, &idx); // WARNING stoi can throw
    if (idx == str.size()) return value;
    else throw std::runtime_error("Cannot convert \"" + str + "\" into an int");
}

template<> bool Config::string_cast(const std::string str)
{
    if      ((str=="true")  || (str=="True")  || (str=="yes") || (str=="Yes") || (str=="1")) return true;
    else if ((str=="false") || (str=="False") || (str=="no")  || (str=="No")  || (str=="0")) return false;
    throw std::runtime_error("Cannot convert \"" + str + "\" into a bool");
}

template<> std::vector<int> Config::string_cast(const std::string str)
{
    auto error = std::runtime_error("Cannot convert \"" + str + "\" into an integer array");
    if (!( (str.front()=='{') && (str.back()=='}'))) throw error;
    std::string new_str = strip(str.substr(1, str.length()-2));
    std::replace(new_str.begin(), new_str.end(), ',', ' ');

    std::vector<int> result;
    while (new_str.length() > 0) {
        size_t idx;
        auto value = std::stoi(new_str, &idx);
        result.push_back(value);
        new_str = new_str.substr(idx, new_str.length()-idx);
    }
    return result;
}

template<> std::vector<double> Config::string_cast(const std::string str)
{
    auto error = std::runtime_error("Cannot convert \"" + str + "\" into an integer array");
    if (!( (str.front()=='{') && (str.back()=='}'))) throw error;
    std::string new_str = strip(str.substr(1, str.length()-2));
    std::replace(new_str.begin(), new_str.end(), ',', ' ');

    std::vector<double> result;
    while (new_str.length() > 0) {
        size_t idx;
        auto value = std::stod(new_str, &idx);
        result.push_back(value);
        new_str = new_str.substr(idx, new_str.length()-idx);
    }
    return result;
}


// For mandatory parameters
template<typename T>
T Config::get_parameter(Dictionary dictionary, std::string name)
{
    auto item = dictionary.find(name);
    if (item==dictionary.end()) throw std::runtime_error("Mandatory parameter " + name + " must be defined");
    else return string_cast<T>((*item).second);
}

// For optional parameters
template<typename T>
T Config::get_parameter(Dictionary dictionary, std::string name, T default_value)
{
    auto item = dictionary.find(name);
    if (item==dictionary.end()) return default_value;
    else return string_cast<T>((*item).second);
}

#include <iostream>

void Config::error_checking()
{
#ifndef HAS_HDF5
    if (output_hdf5)
        throw std::runtime_error("HDF5 output format (output_hdf5=true) requires the code to be compiled with HAS_HDF5");
#endif
    if (live_smbh_count < 0)
        throw std::runtime_error("The number of live black holes (live_smbh_count) has to be greater than or equals to zero");
    if (binary_smbh_pn && (live_smbh_count!=2))
        throw std::runtime_error("Post-Newtonian gravity (binary_smbh_pn=true) requires live_smbh_count=2");
    if (binary_smbh_pn && (pn_c <= 0))
        throw std::runtime_error("Post-Newtonian gravity (binary_smbh_pn=true) requires pn_c > 0");
    if (live_smbh_custom_eps == eps) live_smbh_custom_eps = -1;
    if (live_smbh_output && (live_smbh_count == 0))
        throw std::runtime_error("Black hole output (live_smbh_output=true) requires at least one live black hole (live_smbh_count)");
    if (live_smbh_neighbor_output && (live_smbh_count == 0))
        throw std::runtime_error("Black hole neighbour output (live_smbh_neighbor_output=true) requires at least one live black hole (live_smbh_count)");
    if (binary_smbh_influence_sphere_output && (live_smbh_count != 2))
        throw std::runtime_error("Binary black hole influence sphere output (binary_smbh_influence_sphere_output=true) requires exactly two live black holes (live_smbh_count=2)");
    if (pn_usage.size() != 7)
        throw std::runtime_error("PN usage array (pn_usage) must have exactly seven components");
    if (binary_smbh_pn)
        for (int i=0; i<7; i++)
            if (!((pn_usage[i] == 0) || (pn_usage[i] == 1)))
                throw std::runtime_error("PN usage array (pn_usage) must be a 7-component vector filled with ones and zeros only");
    if ((smbh1_spin.size()!=3) || (smbh2_spin.size()!=3))
        throw std::runtime_error("Spins must be three-component vectors");
    if ((pn_usage[6]==1) && ((smbh1_spin[0]==nix) || (smbh2_spin[0]==nix)))
        throw std::runtime_error("Please define smbh1_spin and smbh2_spin or disable the spin by setting the last component of pn_usage to zero");
    std::cout << smbh1_spin[0] << std::endl;
        std::cout << smbh1_spin[1] << std::endl;
            std::cout << smbh1_spin[2] << std::endl;
    if ((pn_usage[6]==0) && ((smbh1_spin[0]!=nix) || (smbh2_spin[0]!=nix)))
        throw std::runtime_error("Spins (smbh1_spin and smbh2_spin) may not be defined if the last element of pn_usage is set to zero");


    if (ext_units_physical && ((unit_mass == 0) || (unit_length == 0)))
        throw std::runtime_error("Physical units for external gravity (ext_units_physical) requires ext_unit_mass and ext_unit_length to be positive numbers");
    if ((ext_m_bulge > 0) && (ext_b_bulge < 0))
        throw std::runtime_error("To use external bulge gravity, please specify positive ext_m_bulge and ext_b_bulge");
    if ((ext_m_halo_plummer > 0) && (ext_b_halo_plummer < 0))
        throw std::runtime_error("To use external Plummer halo gravity, please specify positive ext_m_halo_plummer and ext_b_halo_plummer");
    if ((ext_m_disk > 0) && ((ext_a_disk < 0) || (ext_b_disk < 0)))
        throw std::runtime_error("To use external disk gravity, please specify positive ext_m_disk, ext_a_disk and ext_b_disk");
    if (((ext_log_halo_r > 0) && (ext_log_halo_v <= 0)) || ((ext_log_halo_r <= 0) && (ext_log_halo_v > 0)))
        throw std::runtime_error("To use external logarithmic halo gravity, please specify positive ext_log_halo_r and ext_log_halo_v");
    if ((ext_dehnen_m > 0) && ((ext_dehnen_r <= 0) || (ext_dehnen_gamma <= 0)))
        throw std::runtime_error("To use external Dehnen model, please specify positive ext_dehnen_r and ext_dehnen_gamma");
}

Config::Config(std::string file_name)
{
    auto dictionary = read_config_file(file_name);

    // TODO check if dt_disk and dt_contr are powers of two
    eps = get_parameter<double>(dictionary, "eps");
    t_end = get_parameter<double>(dictionary, "t_end");
    dt_disk = get_parameter<double>(dictionary, "dt_disk");
    dt_contr = get_parameter<double>(dictionary, "dt_contr");
    dt_bh = get_parameter<double>(dictionary, "dt_bh", dt_contr);
    eta = get_parameter<double>(dictionary, "eta");
    input_file_name = get_parameter<std::string>(dictionary, "input_file_name", "data.con");
    devices_per_node = get_parameter<int>(dictionary, "devices_per_node", 0);
    dt_max_power = get_parameter<int>(dictionary, "dt_max_power", -3);
    dt_min_power = get_parameter<int>(dictionary, "dt_min_power", -36);
    eta_s_corr = get_parameter<double>(dictionary, "eta_s_corr", 4);
    eta_bh_corr = get_parameter<double>(dictionary, "eta_bh_corr", 4);

    output_hdf5 = get_parameter<bool>(dictionary, "output_hdf5", false);
    output_hdf5_double_precision = get_parameter<bool>(dictionary, "output_hdf5_double_precision", true);
    output_ascii_precision = get_parameter<int>(dictionary, "output_ascii_precision", 10);
    output_extra_mode = get_parameter<int>(dictionary, "output_extra_mode", 10);
    dt_min_warning = get_parameter<bool>(dictionary, "dt_min_warning", false);

    live_smbh_count = get_parameter<int>(dictionary, "live_smbh_count", 0);
    live_smbh_custom_eps = get_parameter<double>(dictionary, "live_smbh_custom_eps", -1);
    live_smbh_output = get_parameter<bool>(dictionary, "live_smbh_output", false);
    live_smbh_neighbor_output = get_parameter<bool>(dictionary, "live_smbh_neighbor_output", false);
    live_smbh_neighbor_number = get_parameter<int>(dictionary, "live_smbh_neighbor_number", 10);

    binary_smbh_influence_sphere_output = get_parameter<bool>(dictionary, "binary_smbh_influence_sphere_output", false);
    binary_smbh_influence_radius_factor = get_parameter<double>(dictionary, "binary_smbh_influence_radius_factor", 10.);

    binary_smbh_pn = get_parameter<bool>(dictionary, "binary_smbh_pn", false);
    pn_usage = get_parameter<std::vector<int>>(dictionary, "pn_usage", std::vector<int>({-1,-1,-1,-1,-1,-1,-1}));
    pn_c = get_parameter<double>(dictionary, "pn_c", 0);
    smbh1_spin = get_parameter<std::vector<double>>(dictionary, "smbh1_spin", std::vector<double>({nix,nix,nix}));
    smbh2_spin = get_parameter<std::vector<double>>(dictionary, "smbh2_spin", std::vector<double>({nix,nix,nix}));

    ext_units_physical = get_parameter<bool>(dictionary, "ext_units_physical", false);
    unit_mass = get_parameter<double>(dictionary, "unit_mass", !ext_units_physical);
    unit_length = get_parameter<double>(dictionary, "unit_length", !ext_units_physical);
    ext_m_bulge = get_parameter<double>(dictionary, "ext_m_bulge", 0);
    ext_b_bulge = get_parameter<double>(dictionary, "ext_b_bulge", -1);
    ext_m_disk = get_parameter<double>(dictionary, "ext_m_disk", 0);
    ext_a_disk = get_parameter<double>(dictionary, "ext_a_disk", -1);
    ext_b_disk = get_parameter<double>(dictionary, "ext_b_disk", -1);
    ext_m_halo_plummer = get_parameter<double>(dictionary, "ext_m_halo_plummer", 0);
    ext_b_halo_plummer = get_parameter<double>(dictionary, "ext_b_halo_plummer", -1);
    ext_log_halo_v = get_parameter<double>(dictionary, "ext_log_halo_v", 0);
    ext_log_halo_r = get_parameter<double>(dictionary, "ext_log_halo_r", 0);
    ext_dehnen_m = get_parameter<double>(dictionary, "ext_dehnen_m", 0);
    ext_dehnen_r = get_parameter<double>(dictionary, "ext_dehnen_r", -1);
    ext_dehnen_gamma = get_parameter<double>(dictionary, "ext_dehnen_gamma", -1);

#ifdef ETICS
    dt_scf = get_parameter<double>(dictionary, "dt_scf");
    grapite_mask_file_name = get_parameter<std::string>(dictionary, "grapite_mask_file_name", "grapite.mask");
    etics_dump_coeffs = get_parameter<bool>(dictionary, "etics_dump_coeffs", false);
    grapite_active_search = get_parameter<bool>(dictionary, "grapite_active_search", false);
    grapite_smbh_star_eps = get_parameter<double>(dictionary, "grapite_smbh_star_eps", -1);
    grapite_dev_exec_threshold = get_parameter<int>(dictionary, "grapite_dev_exec_threshold", 32);
#endif

    error_checking();
}