b2case_database.H

//------------------------------------------------------------------------
// b2case_database.H --
//
//
// written by Mathias Doreille
//            Harald Klimach <harald.klimach@dlr.de>
//            Thomas Blome <thomas.blome@dlr.de>
//
// (c) 2004-2012,2017 SMR Engineering & Development SA
//               2502 Bienne, Switzerland
//
// (c) 2023 Deutsches Zentrum für Luft- und Raumfahrt (DLR) e.V.
//          Linder Höhe, 51147 Köln
//
// All Rights Reserved.  Proprietary source code.  The contents of
// this file may not be disclosed to third parties, copied or
// duplicated in any form, in whole or in part, without the prior
// written permission of SMR.
//------------------------------------------------------------------------
 
#ifndef B2CASE_DATABASE_H_
#define B2CASE_DATABASE_H_
 
#include <list>
#include <string>
#include <vector>
 
#include "b2model_database.H"
#include "b2ppconfig.h"
#include "model/b2case.H"
#include "model_imp/b2case_component.H"
#include "muParser.h"
#include "utils/b2allocator.H"
#include "utils/b2dictionary.H"
#include "utils/b2logging.H"
#include "utils/b2splines.H"
#include "utils/b2threading.H"
 
namespace b2000::b2dbv3 {
 
class CaseAttributeList : public b2000::CaseComponent {
public:
    CaseAttributeList();
 
    std::string get_stage_id(int stage = -1) const override = 0;
 
    void push_rtable(const RTable& rtable) { rtable_list.push_back(&rtable); }
 
    void set_case_rtable(const RTable& rtable) { case_rtable = &rtable; }
 
    const RTable& get_case_rtable() { return *case_rtable; }
 
    void set_stage_rtable(const RTable& rtable) { stage_rtable = &rtable; }
 
    const RTable& top_rtable() { return *rtable_list.back(); }
 
    void pop_rtable() { rtable_list.pop_back(); }
 
    void set_cmdline_opt(Dictionary* cmdline_set_) { cmdline_set = cmdline_set_; }
 
    void copy(const CaseAttributeList& case_) {
        rtable_list = case_.rtable_list;
        cmdline_set = case_.cmdline_set;
    }
 
    bool has_key(const std::string& key) const override {
        if (cmdline_set && cmdline_set->has_key(key)) { return true; }
        std::string case_key("CASE");
        case_key += get_id();
        case_key += '.';
        std::string stage_key = case_key + "STAGE" + get_stage_id() + '.' + key;
        case_key += key;
        if (cmdline_set && (cmdline_set->has_key(case_key) || cmdline_set->has_key(stage_key))) {
            return true;
        }
        if (stage_rtable && stage_rtable->has_key(key)) { return true; }
        if (case_rtable && case_rtable->has_key(key)) { return true; }
        if (key_value.has_key(key)) { return true; }
        for (std::vector<const RTable*>::const_reverse_iterator i = rtable_list.rbegin();
             i != rtable_list.rend(); ++i) {
            if ((*i)->has_key(key)) { return true; }
        }
        keys_conversion_t::const_iterator res = keys_conversion.find(key);
        if (res != keys_conversion.end()) { return has_key(res->second); }
        return false;
    }
 
    template <typename T>
    T get(const std::string& key) const {
        std::string case_key("CASE");
        case_key += get_id();
        case_key += '.';
        std::string stage_key = case_key + "STAGE" + get_stage_id() + '.' + key;
        case_key += key;
        if (cmdline_set && cmdline_set->has_key(stage_key)) {
            try {
                T r = cmdline_set->get<T>(stage_key);
                const_cast<RTable&>(solution_rtable).set(key, r);
                return r;
            } catch (Exception& e) {
                logging::get_logger("command_line")
                      << logging::warning << "Cannot parse the command line directive adir: " << e
                      << logging::LOGFLUSH;
            }
        }
        if (key_value.has_key(key)) {
            T r = key_value.get<T>(key);
            const_cast<RTable&>(solution_rtable).set(key, r);
            return r;
        }
        if (stage_rtable && stage_rtable->has_key(key)) {
            T r = stage_rtable->get<T>(key);
            const_cast<RTable&>(solution_rtable).set(key, r);
            return r;
        }
        if (cmdline_set && cmdline_set->has_key(case_key)) {
            try {
                T r = cmdline_set->get<T>(case_key);
                const_cast<RTable&>(solution_rtable).set(key, r);
                return r;
            } catch (Exception& e) {
                logging::get_logger("command_line")
                      << logging::warning << "Cannot parse the command line directive adir: " << e
                      << logging::LOGFLUSH;
            }
        }
        if (case_rtable && case_rtable->has_key(key)) {
            T r = case_rtable->get<T>(key);
            const_cast<RTable&>(solution_rtable).set(key, r);
            return r;
        }
        if (cmdline_set && cmdline_set->has_key(key)) {
            try {
                T r = cmdline_set->get<T>(key);
                const_cast<RTable&>(solution_rtable).set(key, r);
                return r;
            } catch (Exception& e) {
                logging::get_logger("command_line")
                      << logging::warning << "Cannot parse the command line directive adir: " << e
                      << logging::LOGFLUSH;
            }
        }
        for (std::vector<const RTable*>::const_reverse_iterator ii = rtable_list.rbegin();
             ii != rtable_list.rend(); ++ii) {
            if ((*ii)->has_key(key)) {
                T r = (*ii)->get<T>(key);
                const_cast<RTable&>(solution_rtable).set(key, r);
                return r;
            }
        }
        keys_conversion_t::const_iterator res = keys_conversion.find(key);
        if (res != keys_conversion.end()) {
            T r = get<T>(res->second);
            const_cast<RTable&>(solution_rtable).set(key, r);
            return r;
        }
        Exception() << "The key " << key << " of type " << typeid(T)
                    << " is needed but not defined for "
                       "the case "
                    << get_id() << "." << THROW;
        return T();
    }
 
    template <typename T>
    T get(const std::string& key, const T& default_) const {
        if (!has_key(key)) { return default_; }
        return get<T>(key);
    }
 
    std::string get_string(const std::string& key) const override { return get<std::string>(key); }
 
    bool get_bool(const std::string& key) const override { return get<bool>(key); }
 
    int get_int(const std::string& key) const override { return get<int>(key); }
 
    double get_double(const std::string& key) const override { return get<double>(key); }
 
    b2000::csda<double> get_csda_double(const std::string& key) const override {
        return get<b2000::csda<double>>(key);
    }
 
    std::complex<double> get_complex_double(const std::string& key) const override {
        return get<std::complex<double>>(key);
    }
 
    std::string get_string(const std::string& key, const std::string& default_) const override {
        return get<std::string>(key, default_);
    }
 
    bool get_bool(const std::string& key, const bool& default_) const {
        return get<bool>(key, default_);
    }
 
    int get_int(const std::string& key, const int& default_) const {
        return get<int>(key, default_);
    }
 
    double get_double(const std::string& key, const double& default_) const {
        return get<double>(key, default_);
    }
 
    b2000::csda<double> get_csda_double(
          const std::string& key, const b2000::csda<double>& default_) const {
        return get<b2000::csda<double>>(key, default_);
    }
 
    std::complex<double> get_complex_double(
          const std::string& key, const std::complex<double>& default_) const {
        return get<std::complex<double>>(key, default_);
    }
 
    template <typename T>
    void get(const std::string& key, std::vector<T>& v) const {
        std::string case_key("CASE");
        case_key += get_id();
        case_key += '.';
        std::string stage_key = case_key + "STAGE" + get_stage_id() + '.' + key;
        case_key += key;
        if (cmdline_set && cmdline_set->has_key(stage_key)) {
            try {
                cmdline_set->get(stage_key, v);
                const_cast<RTable&>(solution_rtable).set(key, v);
                return;
            } catch (Exception& e) {
                logging::get_logger("command_line")
                      << logging::warning << "Cannot parse the command line directive adir: " << e
                      << logging::LOGFLUSH;
            }
        }
        if (stage_rtable && stage_rtable->has_key(key)) {
            stage_rtable->get(key, v);
            const_cast<RTable&>(solution_rtable).set(key, v);
            return;
        }
        if (cmdline_set && cmdline_set->has_key(case_key)) {
            try {
                cmdline_set->get(case_key, v);
                const_cast<RTable&>(solution_rtable).set(key, v);
                return;
            } catch (Exception& e) {
                logging::get_logger("command_line")
                      << logging::warning << "Cannot parse the command line directive adir: " << e
                      << logging::LOGFLUSH;
            }
        }
        if (case_rtable && case_rtable->has_key(key)) {
            case_rtable->get(key, v);
            const_cast<RTable&>(solution_rtable).set(key, v);
            return;
        }
        if (cmdline_set && cmdline_set->has_key(key)) {
            try {
                cmdline_set->get(key, v);
                const_cast<RTable&>(solution_rtable).set(key, v);
                return;
            } catch (Exception& e) {
                logging::get_logger("command_line")
                      << logging::warning << "Cannot parse the command line directive adir: " << e
                      << logging::LOGFLUSH;
            }
        }
        for (std::vector<const RTable*>::const_reverse_iterator i = rtable_list.rbegin();
             i != rtable_list.rend(); ++i) {
            if ((*i)->has_key(key)) {
                (*i)->get(key, v);
                const_cast<RTable&>(solution_rtable).set(key, v);
                return;
            }
        }
        keys_conversion_t::const_iterator res = keys_conversion.find(key);
        if (res != keys_conversion.end()) {
            get(res->second, v);
            const_cast<RTable&>(solution_rtable).set(key, v);
            return;
        }
        Exception() << "The key " << key << " of type array of " << typeid(T)
                    << " is needed but not "
                       "defined for the case "
                    << get_id() << "." << THROW;
    }
 
    void set(const std::string& key, int value) override { key_value.set(key, value); }
 
    void set(const std::string& key, double value) override { key_value.set(key, value); }
 
    void set(const std::string& key, b2000::csda<double> value) override {
        key_value.set(key, value);
    }
 
    void set(const std::string& key, std::complex<double> value) override {
        key_value.set(key, value);
    }
 
    void set(const std::string& key, const std::string& value) override {
        key_value.set(key, value);
    }
 
    void get(const std::string& key, std::vector<int>& v) const override { get<int>(key, v); }
 
    void get(const std::string& key, std::vector<double>& v) const override { get<double>(key, v); }
 
    void get(const std::string& key, std::vector<b2000::csda<double>>& v) const override {
        get<b2000::csda<double>>(key, v);
    }
 
    void get(const std::string& key, std::vector<std::complex<double>>& v) const override {
        get<std::complex<double>>(key, v);
    }
 
    void warn_on_non_used_key() const override {
        if (cmdline_set) { cmdline_set->warn_on_non_used_key(); }
    }
 
    RTable& get_solution_rtable() { return solution_rtable; }
 
private:
    std::vector<const RTable*> rtable_list;
    const RTable* case_rtable;
    const RTable* stage_rtable;
    Dictionary* cmdline_set;
    using keys_conversion_t = std::map<std::string, std::string>;
    keys_conversion_t keys_conversion;
    RTable key_value;
    RTable solution_rtable;
};
 
class Case : public CaseAttributeList {
public:
    ~Case() override {
        for (auto& comp : list_comp) { comp->decref_or_free(); }
        for (auto& time_func : list_time_function) { delete time_func; }
    }
    void init(Model& model, const std::string& db_id, Allocator& allocator);
 
    // Add component to list_comp with information from RTable and subtype_name
    void set_component(
          Model& model, const std::string& case_id, const RTable& rt,
          const std::string& subtype_name, double scalef, Allocator& allocator);
 
    // Get component from list_comp using component id
    Component& get_component(const std::string& comp_id);
 
    std::string get_id() const override { return id; }
 
    void get_subcase_ids(std::set<int>& subcase_ids_) const override { subcase_ids_ = subcase_ids; }
 
    size_t get_number_of_stage() const override { return 1; }
 
    std::string get_stage_id(int stage = -1) const override { return id; }
 
    int get_stage_number() const override { return 0; }
 
    double get_stage_size(int stage = -1) const override {
        double pamax = get_double("PAMAX", 1.0);
        if (pamax == 0.0) { pamax = 1.0; }
        return pamax;
    }
 
    bool next_stage() override { return false; }
 
    void set_stage(int stage) override { return; }
 
    class CaseComponentIterator : public ComponentIterator {
    public:
        Component* next() override {
            ensure_matching_subcase_id();
            if (begin == end) { return nullptr; }
            return *begin++;
        }
 
    private:
        friend class b2dbv3::Case;
        CaseComponentIterator(Component** begin_, Component** end_, const int subcase_id_)
            : begin(begin_), end(end_), subcase_id(subcase_id_) {
            ensure_matching_subcase_id();
        }
 
        void ensure_matching_subcase_id() {
            if (subcase_id < 0) { return; }
            for (; begin != end; ++begin) {
                if (*begin != nullptr && (*begin)->get_subcase_id() == subcase_id) { break; }
            }
        }
 
        using list_comp_t = std::vector<Component*>;
        Component** begin;
        Component** end;
        int subcase_id;
    };
 
    component_iterator get_component_iterator() override {
        return component_iterator(
              new CaseComponentIterator(&list_comp.front(), &list_comp.back() + 1, -1));
    }
 
    component_iterator get_component_iterator_for_stage(size_t stage) override {
        return get_component_iterator();
    }
 
    component_iterator get_component_iterator_for_stage_and_subcase(
          size_t stage_, int subcase_id_) override {
        if (stage_ != 0 && subcase_id_ != 0) {
            UnimplementedError() << "Stages and subcases cannot be "
                                    "combined."
                                 << THROW;
        }
        return component_iterator(
              new CaseComponentIterator(&list_comp.front(), &list_comp.back() + 1, subcase_id_));
    }
 
    using type_t = ObjectTypeComplete<Case>;
    static type_t type;
 
protected:
    class TimeFunction : public b2000::TimeFunction {
    public:
        TimeFunction(double scale_) : scale(scale_) {}
 
        double d_value(double t, int order) override {
            // Numerical computation of the derivative of the
            // expression.
            static const int binomial_coefficient[7][7] = {
                  {1},
                  {1, -1},
                  {1, -2, 1},
                  {1, -3, 3, -1},
                  {1, -4, 6, -4, 1},
                  {1, -5, 10, -10, 5, -1},
                  {1, -6, 15, -20, 15, -6, 1}};
            if (order > 6) { UnimplementedError() << THROW; }
            double h = 0.0001;
            double vt = value(t);
            if (vt != vt) {
                // The evaluation of the python expression is
                // not a number.
                return vt;
            }
            double r = 0;
            while (1) {
                double new_r = 0;
                for (int i = 0; i <= order; ++i) {
                    new_r += binomial_coefficient[order][i] * value(t + (order / 2.0 - i) * h);
                }
                for (int i = 0; i != order; ++i) { new_r /= h; }
                if (std::abs(new_r - r) < 0.0000001) { return new_r; }
                r = new_r;
                h /= 100;
            }
            return -1;
        }
 
    protected:
        double scale;
    };
 
    class TimeFunctionIdentity : public TimeFunction {
    public:
        TimeFunctionIdentity(double scale_) : TimeFunction(scale_) {}
 
        double value(double t) override { return scale * t; }
 
        double d_value(double t, int order) override {
            if (order == 1) { return scale; }
            if (order > 1) { return 0; }
            Exception() << THROW;
            return 0;
        }
    };
 
    class TimeFunctionExpression : public TimeFunction {
    public:
        TimeFunctionExpression(double scale_, const std::string& expression_string)
            : TimeFunction(scale_) {
            expression.SetExpr(expression_string);
            expression.DefineVar("t", &var_t);
            expression.SetVarFactory(AddVariable, this);
        }
 
        void set_case(b2000::Case& case_) override {
            expression.ClearConst();
            expression.DefineConst("s", case_.get_stage_number() + 1);
            expression.DefineStrConst("sid", case_.get_stage_id());
            expression.DefineConst("ssize", case_.get_stage_size());
        }
 
        double value(double t) override {
            try {
                b2spin_rw_mutex::scoped_lock lock(mutex);
 
                var_t = t;
                return scale * expression.Eval();
            } catch (mu::Parser::exception_type& e) {
                Exception() << "In the expression \"" << expression.GetExpr()
                            << "\": " << e.GetMsg() << THROW;
            }
            return 0;
        }
 
    private:
        b2spin_rw_mutex mutex;
 
        static double* AddVariable(const char* a_szName, void* pUserData) {
            TimeFunctionExpression* this_ = static_cast<TimeFunctionExpression*>(pUserData);
            std::deque<double>& dyn_var = this_->dyn_var;
 
            b2spin_rw_mutex::scoped_lock lock(this_->mutex);
 
            dyn_var.push_back(0);
            return &(dyn_var.back());
        }
        double var_t;
        mu::Parser expression;
        std::deque<double> dyn_var;
    };
 
    class TimeFunctionTabulated : public TimeFunction {
    public:
        TimeFunctionTabulated(double scale_, const std::vector<double>& tabulated_vector)
            : TimeFunction(scale_),
              spline(tabulated_vector.begin(), tabulated_vector.end(), true) {}
 
        double value(double t) override { return scale * spline(t); }
 
    private:
        CubicSpline spline;
    };
 
    RTable rtable;
    std::string id;
    using list_comp_t = std::vector<Component*>;
    list_comp_t list_comp;
 
    using list_time_function_t = std::vector<TimeFunction*>;
    list_time_function_t list_time_function;
    std::set<int> subcase_ids;
 
    void add_case(Model& model, const std::string& set_id, double scalef, Allocator& allocator);
};
 
class CaseStage : public Case {
public:
    void init(Model& model, const std::string& db_id, Allocator& allocator);
 
    size_t get_number_of_stage() const override { return list_stage.size(); }
 
    std::string get_stage_id(int stage_ = -1) const override {
        if (stage_ == -1) {
            return stage->first.get_id();
        } else {
            return get_stage(stage_)->first.get_id();
        }
    }
 
    int get_stage_number() const override { return stage_number; }
 
    double get_stage_size(int stage_ = -1) const override {
        if (stage_ == -1) {
            return stage->second;
        } else {
            return get_stage(stage_)->second;
        }
    }
 
    bool next_stage() override;
 
    void set_stage(int stage_) override;
 
    component_iterator get_component_iterator() override {
        return stage->first.get_component_iterator();
    }
 
    component_iterator get_component_iterator_for_stage(size_t stage_) override {
        return get_stage(stage_)->first.get_component_iterator();
    }
 
    component_iterator get_component_iterator_for_stage_and_subcase(
          size_t stage_, int subcase_id_) override {
        if (stage_ != 0 && subcase_id_ != 0) {
            UnimplementedError() << "Stages and subcases cannot be "
                                    "combined."
                                 << THROW;
        }
        return get_stage(stage_)->first.get_component_iterator();
    }
 
    using type_t = ObjectTypeComplete<CaseStage>;
    static type_t type;
 
private:
    void add_stage(Model& model, const std::string& set_id, double scalef, Allocator& allocator);
 
    using list_stage_t = std::list<std::pair<Case, double>>;
    list_stage_t list_stage;
    list_stage_t::iterator stage;
    int stage_number;
 
    list_stage_t::const_iterator get_stage(int stage_) const {
        if (stage_ < 0 || size_t(stage_) >= list_stage.size()) { Exception() << THROW; }
        list_stage_t::const_iterator i = list_stage.begin();
        while (stage_--) { ++i; }
        return i;
    }
 
    list_stage_t::iterator get_stage(int stage_) {
        if (stage_ < 0 || size_t(stage_) >= list_stage.size()) { Exception() << THROW; }
        list_stage_t::iterator i = list_stage.begin();
        while (stage_--) { ++i; }
        return i;
    }
};
 
class CaseList : public b2000::CaseList, protected std::vector<Case*> {
public:
    void set_model(b2000::Model& model) override {}
 
    size_t get_number_of_case() const override { return size(); }
 
    class CaseIterator : public b2000::CaseList::CaseIterator {
    public:
        Case* next() override {
            if (begin != end) { return *begin++; }
            return nullptr;
        }
        size_t get_number_of_case() const override { return end - begin; }
 
    private:
        friend class b2dbv3::CaseList;
        CaseIterator(std::vector<Case*>& cl) : begin(cl.begin()), end(cl.end()) {}
        std::vector<Case*>::iterator begin;
        std::vector<Case*>::iterator end;
    };
 
    b2000::CaseList::case_iterator get_case_iterator() override {
        return case_iterator(new CaseIterator(*this));
    }
 
    b2000::CaseList::case_list_iterator get_case_sublist_iterator(
          const b2000::Domain* domain, b2000::CaseList::SortType st) override {
        return case_list_iterator(new CaseListIterator(*this, std::less<Case*>()));
    }
 
    b2000::CaseList::case_list_iterator get_case_sublist_iterator(
          const NaturalBoundaryCondition* nbc, b2000::CaseList::SortType st) override {
        return case_list_iterator(new CaseListIterator(*this, std::less<Case*>()));
    }
 
    b2000::CaseList::case_list_iterator get_case_sublist_iterator(
          const EssentialBoundaryCondition* ebc, b2000::CaseList::SortType st) override {
        return case_list_iterator(new CaseListIterator(*this, std::less<Case*>()));
    }
 
    b2000::CaseList::case_list_iterator get_case_sublist_iterator(
          const ModelReductionBoundaryCondition* mrbc, b2000::CaseList::SortType st) override {
        return case_list_iterator(new CaseListIterator(*this, std::less<Case*>()));
    }
 
    case_list_iterator get_case_sublist_iterator(
          const InitialCondition* ic, b2000::CaseList::SortType st) override {
        return case_list_iterator(new CaseListIterator(*this, std::less<Case*>()));
    }
 
    class CaseListIterator : public b2000::CaseList::CaseListIterator {
    public:
        CaseList* next() override {
            if (iter != list_of_case_list.end()) { return &(*(iter++)); }
            return nullptr;
        }
        size_t get_number_of_case_list() const override { return list_of_case_list.size(); }
 
    private:
        friend class b2dbv3::CaseList;
 
        template <typename CASE_COMP>
        CaseListIterator(std::vector<Case*>& clist, const CASE_COMP& comp) {
            std::vector<Case*> tmp(clist);
            std::vector<Case*>::iterator i = tmp.begin();
            std::sort(i, tmp.end(), comp);
            list_of_case_list.push_back(CaseList());
            list_of_case_list.back().push_back(*i++);
            for (; i != tmp.end(); ++i) {
                if (comp(list_of_case_list.back().front(), *i)) {
                    list_of_case_list.push_back(CaseList());
                }
                list_of_case_list.back().push_back(*i);
            }
            iter = list_of_case_list.begin();
        }
        std::list<CaseList>::iterator iter;
        std::list<CaseList> list_of_case_list;
    };
 
    friend class CaseListIterator;
 
    using type_t = ObjectTypeComplete<CaseList, CaseList::type_t>;
    static type_t type;
};
 
class CaseListRoot : public CaseList {
public:
    ~CaseListRoot() override {
        for (std::vector<Case*>::iterator i = begin(); i != end(); ++i) { (*i)->decref_or_free(); }
    }
 
    void set_model(b2000::Model& model) override {}
 
    void init(Model& model, Dictionary* cmdline_set_, const RTable& adir, const RTable& strat);
 
    using type_t = ObjectTypeComplete<CaseListRoot, CaseList::type_t>;
    static type_t type;
 
private:
    Allocator allocator;
};
 
}  // namespace b2000::b2dbv3
 
#endif  // B2CASE_DATABASE_H_