multifluid.hpp

    // Extract the set of possible identifiers to be used to match parameters
    auto identifierset = collect_identifiers(pureJSON);
    // Decide which identifier is to be used (Name, CAS, REFPROP name)
    auto identifiers = identifierset[select_identifier(BIPcollection, identifierset, flags)];

    // Things related to the mixture
    auto F = MultiFluidReducingFunction::get_F_matrix(BIPcollection, identifiers, flags);
    auto funcs = get_departure_function_matrix(depcollection, BIPcollection, identifiers, flags);
    auto [betaT, gammaT, betaV, gammaV] = MultiFluidReducingFunction::get_BIP_matrices(BIPcollection, identifiers, flags, Tc, vc);

    auto redfunc = MultiFluidReducingFunction(betaT, gammaT, betaV, gammaV, Tc, vc);

    return MultiFluid(
        std::move(redfunc),
        std::move(CorrespondingStatesContribution(std::move(EOSs))),
        std::move(DepartureContribution(std::move(F), std::move(funcs)))
    );
}

/// A builder function where the JSON-formatted strings are provided explicitly rather than file paths
inline auto build_multifluid_JSONstr(const std::vector<std::string>& componentJSON, const std::string& BIPJSON, const std::string& departureJSON, const nlohmann::json& flags = {}) {

    // Mixture things
    const auto BIPcollection = nlohmann::json::parse(BIPJSON);
    const auto depcollection = nlohmann::json::parse(departureJSON);

    // Pure fluids
    std::vector<nlohmann::json> pureJSON;
    for (auto& c : componentJSON) {
        pureJSON.emplace_back(nlohmann::json::parse(c));
    }
    return _build_multifluid_model(pureJSON, BIPcollection, depcollection, flags);
}

inline auto build_multifluid_model(const std::vector<std::string>& components, const std::string& coolprop_root, const std::string& BIPcollectionpath = {}, const nlohmann::json& flags = {}, const std::string& departurepath = {}) {

    std::string BIPpath = (BIPcollectionpath.empty()) ? coolprop_root + "/dev/mixtures/mixture_binary_pairs.json" : BIPcollectionpath;
    const auto BIPcollection = load_a_JSON_file(BIPpath);

    std::string deppath = (departurepath.empty()) ? coolprop_root + "/dev/mixtures/mixture_departure_functions.json" : departurepath;
    const auto depcollection = load_a_JSON_file(deppath);

    // Pure fluids
    std::vector<nlohmann::json> pureJSON;
    try {
        // Try the normal lookup, matching component name to a file in dev/fluids (case sensitive match on linux!)
        pureJSON = collect_component_json(components, coolprop_root);
    }
    catch(...){
        // Lookup the absolute paths for each component
        auto aliasmap = build_alias_map(coolprop_root);
        std::vector<std::string> abspaths;
        for (auto c : components) {
            // Allow matching of absolute paths first
            if (std::filesystem::is_regular_file(c)) {
                abspaths.push_back(c);
            }
            else {
                abspaths.push_back(aliasmap[c]);
            }
        }
        // Backup lookup with absolute paths resolved for each component
        pureJSON = collect_component_json(abspaths, coolprop_root);
    }
    return _build_multifluid_model(pureJSON, BIPcollection, depcollection, flags);
}

/**
* \brief Load a model from a JSON data structure
* 
* Required fields are: components, BIP, departure
* 
* BIP and departure can be either the data in JSON format, or a path to file with those contents
* components is an array, which either contains the paths to the JSON data, or the file path
*/
inline auto multifluidfactory(const nlohmann::json& spec) {
    
    auto JSON_from_path_or_contents = [](const nlohmann::json &path_or_contents) -> nlohmann::json {
        if (path_or_contents.is_string()) {
            return load_a_JSON_file(path_or_contents.get<std::string>());
        }
        else {
            return path_or_contents;
        }
    };

    auto components = spec.at("components"); 
    auto depcollection = JSON_from_path_or_contents(spec.at("departure"));
    auto BIPcollection = JSON_from_path_or_contents(spec.at("BIP"));
    nlohmann::json flags = (spec.contains("flags")) ? spec.at("flags") : nlohmann::json();

    // Pure components
    std::vector<nlohmann::json> pureJSON;
    for (auto c : components) {
        pureJSON.push_back(JSON_from_path_or_contents(c));
    }

    return _build_multifluid_model(pureJSON, BIPcollection, depcollection, flags);
}
/// An overload of multifluidfactory that takes in a string
inline auto multifluidfactory(const std::string& specstring) {
    return multifluidfactory(nlohmann::json::parse(specstring));
}

/**
This class holds a lightweight reference to the core parts of the model

The reducing and departure functions are moved into this class, while the donor class is used for the corresponding states portion
*/
template<typename ReducingFunction, typename DepartureFunction, typename BaseClass>
class MultiFluidAdapter {

private:
    std::string meta = "";

public:
    const BaseClass& base;
    const ReducingFunction redfunc;
    const DepartureFunction depfunc;

    template<class VecType>
    auto R(const VecType& molefrac) const { return base.R(molefrac); }

    MultiFluidAdapter(const BaseClass& base, ReducingFunction&& redfunc, DepartureFunction &&depfunc) : base(base), redfunc(redfunc), depfunc(depfunc) {};

    /// Store some sort of metadata in string form (perhaps a JSON representation of the model?)
    void set_meta(const std::string& m) { meta = m; }
    /// Get the metadata stored in string form
    auto get_meta() const { return meta; }

    template<typename TType, typename RhoType, typename MoleFracType>
    auto alphar(const TType& T,
        const RhoType& rho,
        const MoleFracType& molefrac) const
    {
        auto Tred = forceeval(redfunc.get_Tr(molefrac));
        auto rhored = forceeval(redfunc.get_rhor(molefrac));
        auto delta = forceeval(rho / rhored);
        auto tau = forceeval(Tred / T);
        auto val = base.corr.alphar(tau, delta, molefrac) + depfunc.alphar(tau, delta, molefrac);
        return forceeval(val);
    }
};

template<class Model>
auto build_multifluid_mutant(Model& model, const nlohmann::json& jj) {

    auto red = model.redfunc;
    auto N = red.Tc.size();

    auto betaT = red.betaT, gammaT = red.gammaT, betaV = red.betaV, gammaV = red.gammaV;
    auto Tc = red.Tc, vc = red.vc;

    // Allocate the matrices of default models and F factors
    Eigen::MatrixXd F(N, N); F.setZero();
    std::vector<std::vector<DepartureTerms>> funcs(N);
    for (auto i = 0; i < N; ++i) { funcs[i].resize(N); }

    for (auto i = 0; i < N; ++i) {
        for (auto j = i; j < N; ++j) {
            if (i == j) {
                funcs[i][i].add_term(NullEOSTerm());
            }
            else {
                // Extract the given entry
                auto entry = jj[std::to_string(i)][std::to_string(j)];

                auto BIP = entry["BIP"];
                // Set the reducing function parameters in the copy
                betaT(i, j) = BIP.at("betaT");
                betaT(j, i) = 1 / red.betaT(i, j);
                betaV(i, j) = BIP.at("betaV");
                betaV(j, i) = 1 / red.betaV(i, j);
                gammaT(i, j) = BIP.at("gammaT"); gammaT(j, i) = gammaT(i, j);
                gammaV(i, j) = BIP.at("gammaV"); gammaV(j, i) = gammaV(i, j);

                // Build the matrix of F and departure functions
                auto dep = entry["departure"];
                F(i, j) = BIP.at("Fij");
                F(j, i) = F(i, j);
                funcs[i][j] = build_departure_function(dep);
                funcs[j][i] = build_departure_function(dep);
            }
        }
    }

    auto newred = MultiFluidReducingFunction(betaT, gammaT, betaV, gammaV, Tc, vc);
    auto newdep = DepartureContribution(std::move(F), std::move(funcs));
    auto mfa = MultiFluidAdapter(model, std::move(newred), std::move(newdep));
    /// Store the departure term in the adapted multifluid class
    mfa.set_meta(jj.dump());
    return mfa;
}

template<typename Model>
auto build_multifluid_mutant_invariant(Model& model, const nlohmann::json& jj) {

    auto red = model.redfunc;
    auto N = red.Tc.size();
    if (N != 2) {
        throw std::invalid_argument("Only binary mixtures are currently supported with invariant departure functions");
    }

    auto phiT = red.betaT, lambdaT = red.gammaT, phiV = red.betaV, lambdaV = red.gammaV;
    phiT.setOnes(); phiV.setOnes();
    lambdaV.setZero(); lambdaV.setZero();
    auto Tc = red.Tc, vc = red.vc;

    // Allocate the matrices of default models and F factors
    Eigen::MatrixXd F(N, N); F.setZero();
    std::vector<std::vector<DepartureTerms>> funcs(N);
    for (auto i = 0; i < N; ++i) { funcs[i].resize(N); }

    for (auto i = 0; i < N; ++i) {
        for (auto j = i; j < N; ++j) {
            if (i == j) {
                funcs[i][i].add_term(NullEOSTerm());
            }
            else {
                // Extract the given entry
                auto entry = jj[std::to_string(i)][std::to_string(j)];

                auto BIP = entry["BIP"];
                // Set the reducing function parameters in the copy
                phiT(i, j) = BIP.at("phiT");
                phiT(j, i) = phiT(i, j);
                lambdaT(i, j) = BIP.at("lambdaT"); 
                lambdaT(j, i) = -lambdaT(i, j);

                phiV(i, j) = BIP.at("phiV");
                phiV(j, i) = phiV(i, j);
                lambdaV(i, j) = BIP.at("lambdaV"); 
                lambdaV(j, i) = -lambdaV(i, j);

                // Build the matrix of F and departure functions
                auto dep = entry["departure"];
                F(i, j) = BIP.at("Fij");
                F(j, i) = F(i, j);
                funcs[i][j] = build_departure_function(dep);
                funcs[j][i] = build_departure_function(dep);
            }
        }
    }

    auto newred = MultiFluidInvariantReducingFunction(phiT, lambdaT, phiV, lambdaV, Tc, vc);
    auto newdep = DepartureContribution(std::move(F), std::move(funcs));
    auto mfa = MultiFluidAdapter(model, std::move(newred), std::move(newdep));
    /// Store the departure term in the adapted multifluid class
    mfa.set_meta(jj.dump());
    return mfa;
}


class DummyEOS {
public:
    template<typename TType, typename RhoType> auto alphar(TType tau, const RhoType& delta) const { return tau * delta; }
};
class DummyReducingFunction {
public:
    template<typename MoleFractions> auto get_Tr(const MoleFractions& molefracs) const { return molefracs[0]; }
    template<typename MoleFractions> auto get_rhor(const MoleFractions& molefracs) const { return molefracs[0]; }
};
inline auto build_dummy_multifluid_model(const std::vector<std::string>& components) {
    std::vector<DummyEOS> EOSs(2);
    std::vector<std::vector<DummyEOS>> funcs(2); for (auto i = 0; i < funcs.size(); ++i) { funcs[i].resize(funcs.size()); }
    std::vector<std::vector<double>> F(2); for (auto i = 0; i < F.size(); ++i) { F[i].resize(F.size()); }

    struct Fwrapper {
    private: 
        const std::vector<std::vector<double>> F_;
    public:
        Fwrapper(const std::vector<std::vector<double>> &F) : F_(F){};
        auto operator ()(std::size_t i, std::size_t j) const{ return F_[i][j]; }
    };
    auto ff = Fwrapper(F);
    auto redfunc = DummyReducingFunction();
    return MultiFluid(std::move(redfunc), std::move(CorrespondingStatesContribution(std::move(EOSs))), std::move(DepartureContribution(std::move(ff), std::move(funcs))));
}
inline void test_dummy() {
    auto model = build_dummy_multifluid_model({ "A", "B" });
    std::valarray<double> rhovec = { 1.0, 2.0 };
    auto alphar = model.alphar(300.0, rhovec);
}

}; // namespace teqp