From 403a0a59c79a0f5c2e32ac34eebe20138d5e78e2 Mon Sep 17 00:00:00 2001
From: Ian Bell <ian.bell@nist.gov>
Date: Wed, 9 Nov 2022 21:21:46 -0500
Subject: [PATCH] Move implementation of C++ interface into separate files
 (#28)

* Try to break up the C++ interface into smaller chunks

In order to make compilation less RAM intensive

* Inline the makers, also add the other sources

* Tidying up the splitting
---
 CMakeLists.txt                             |   7 +-
 include/teqp/models/mie/lennardjones.hpp   |   4 +-
 interface/CPP/teqp_impl_Arxy.cpp           |  59 +++++
 interface/CPP/teqp_impl_VLE.cpp            |  52 +++++
 interface/CPP/teqp_impl_VLLE.cpp           |  17 ++
 interface/CPP/teqp_impl_crit.cpp           |  61 +++++
 interface/CPP/teqp_impl_factory.cpp        |  15 ++
 interface/CPP/teqp_impl_isochoric.cpp      |  27 +++
 interface/CPP/teqp_impl_virials.cpp        |  29 +++
 interface/CPP/teqpcpp.cpp                  | 253 ++++-----------------
 interface/CPP/{ => test}/bench_teqpcpp.cpp |   0
 interface/CPP/{ => test}/test_teqpcpp.cpp  |   0
 12 files changed, 309 insertions(+), 215 deletions(-)
 create mode 100644 interface/CPP/teqp_impl_Arxy.cpp
 create mode 100644 interface/CPP/teqp_impl_VLE.cpp
 create mode 100644 interface/CPP/teqp_impl_VLLE.cpp
 create mode 100644 interface/CPP/teqp_impl_crit.cpp
 create mode 100644 interface/CPP/teqp_impl_factory.cpp
 create mode 100644 interface/CPP/teqp_impl_isochoric.cpp
 create mode 100644 interface/CPP/teqp_impl_virials.cpp
 rename interface/CPP/{ => test}/bench_teqpcpp.cpp (100%)
 rename interface/CPP/{ => test}/test_teqpcpp.cpp (100%)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 9447908..4926748 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -91,16 +91,17 @@ if (NOT TEQP_NO_TEQPCPP)
   # Add a static library with the C++ interface that uses only STL
   # types so that recompilation of a library that uses teqp 
   # doesn't require a full compile for a single LOC change
-  add_library(teqpcpp STATIC "${CMAKE_CURRENT_SOURCE_DIR}/interface/CPP/teqpcpp.cpp")
+  file(GLOB sources "${CMAKE_CURRENT_SOURCE_DIR}/interface/CPP/*.cpp")
+  add_library(teqpcpp STATIC ${sources})
   target_link_libraries(teqpcpp PUBLIC teqpinterface PUBLIC autodiff)
   target_include_directories(teqpcpp PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}/interface/CPP")
   set_property(TARGET teqpcpp PROPERTY POSITION_INDEPENDENT_CODE ON)
   target_compile_definitions(teqpcpp PUBLIC -DMULTICOMPLEX_NO_MULTIPRECISION)
 
   if (TEQP_TESTTEQPCPP)
-    add_executable(test_teqpcpp "${CMAKE_CURRENT_SOURCE_DIR}/interface/CPP/test_teqpcpp.cpp")
+    add_executable(test_teqpcpp "${CMAKE_CURRENT_SOURCE_DIR}/interface/CPP/test/test_teqpcpp.cpp")
     target_link_libraries(test_teqpcpp PUBLIC teqpcpp)
-    add_executable(bench_teqpcpp "${CMAKE_CURRENT_SOURCE_DIR}/interface/CPP/bench_teqpcpp.cpp")
+    add_executable(bench_teqpcpp "${CMAKE_CURRENT_SOURCE_DIR}/interface/CPP/test/bench_teqpcpp.cpp")
     target_link_libraries(bench_teqpcpp PUBLIC teqpcpp PRIVATE Catch2WithMain)
   endif()
 endif()
diff --git a/include/teqp/models/mie/lennardjones.hpp b/include/teqp/models/mie/lennardjones.hpp
index d97d9d7..4cdfe39 100644
--- a/include/teqp/models/mie/lennardjones.hpp
+++ b/include/teqp/models/mie/lennardjones.hpp
@@ -6,7 +6,7 @@ namespace teqp {
     /**
      * The EOS of Monika Thol and colleagues. DOI:10.1063/1.4945000
      */
-    auto build_LJ126_TholJPCRD2016() {
+    inline auto build_LJ126_TholJPCRD2016() {
         std::string contents = R"(
 
         {
@@ -65,4 +65,4 @@ namespace teqp {
 
         return teqp::build_multifluid_JSONstr({ contents }, "{}", "{}");
     }
-};
\ No newline at end of file
+};
diff --git a/interface/CPP/teqp_impl_Arxy.cpp b/interface/CPP/teqp_impl_Arxy.cpp
new file mode 100644
index 0000000..ff2fcd2
--- /dev/null
+++ b/interface/CPP/teqp_impl_Arxy.cpp
@@ -0,0 +1,59 @@
+#include "teqp/derivs.hpp"
+
+#include "teqpcpp.cpp"
+using MI = teqp::cppinterface::ModelImplementer;
+
+// Derivatives from isochoric thermodynamics (all have the same signature)
+#define X(f) \
+double MI::f(const double T, const REArrayd& rhovec) const  { \
+    return std::visit([&](const auto& model) { \
+        using id = IsochoricDerivatives<decltype(model), double, REArrayd>; \
+        return id::f(model, T, rhovec); \
+    }, m_model); \
+}
+ISOCHORIC_double_args
+#undef X
+
+#define X(f) \
+EArrayd MI::f(const double T, const REArrayd& rhovec) const  { \
+    return std::visit([&](const auto& model) { \
+        using id = IsochoricDerivatives<decltype(model), double, REArrayd>; \
+        return id::f(model, T, rhovec); \
+    }, m_model); \
+}
+ISOCHORIC_array_args
+#undef X
+
+#define X(f) \
+EMatrixd MI::f(const double T, const REArrayd& rhovec) const  { \
+    return std::visit([&](const auto& model) { \
+        using id = IsochoricDerivatives<decltype(model), double, REArrayd>; \
+        return id::f(model, T, rhovec); \
+    }, m_model); \
+}
+ISOCHORIC_matrix_args
+#undef X
+
+double MI::get_Arxy(const int NT, const int ND, const double T, const double rho, const EArrayd& molefracs) const {
+    return std::visit([&](const auto& model) {
+        using tdx = teqp::TDXDerivatives<decltype(model), double, EArrayd>;
+        return tdx::template get_Ar(NT, ND, model, T, rho, molefracs);
+    }, m_model);
+}
+
+double MI::get_neff(const double T, const double rho, const EArrayd& molefracs) const {
+    return std::visit([&](const auto& model) {
+        using tdx = teqp::TDXDerivatives<decltype(model), double, EArrayd>;
+        return tdx::template get_neff(model, T, rho, molefracs);
+    }, m_model);
+}
+
+EArray33d MI::get_deriv_mat2(const double T, double rho, const EArrayd& z) const {
+    return std::visit([&](const auto& model) {
+        // Although the template argument suggests that only residual terms
+        // are returned, also the ideal-gas ones are returned because the
+        // ideal-gas term is required to implement alphar which just redirects
+        // to alphaig
+        return DerivativeHolderSquare<2, AlphaWrapperOption::residual>(model, T, rho, z).derivs;
+    }, m_model);
+}
diff --git a/interface/CPP/teqp_impl_VLE.cpp b/interface/CPP/teqp_impl_VLE.cpp
new file mode 100644
index 0000000..017d761
--- /dev/null
+++ b/interface/CPP/teqp_impl_VLE.cpp
@@ -0,0 +1,52 @@
+#include "teqpcpp.cpp"
+#include "teqp/algorithms/VLE.hpp"
+
+using namespace teqp;
+using MI = teqp::cppinterface::ModelImplementer;
+
+std::tuple<EArrayd, EArrayd> MI::get_drhovecdp_Tsat(const double T, const REArrayd& rhovecL, const REArrayd& rhovecV) const {
+    return std::visit([&](const auto& model) {
+        return teqp::get_drhovecdp_Tsat(model, T, rhovecL, rhovecV);
+    }, m_model);
+}
+std::tuple<EArrayd, EArrayd> MI::get_drhovecdT_psat(const double T, const REArrayd& rhovecL, const REArrayd& rhovecV) const {
+    return std::visit([&](const auto& model) {
+        return teqp::get_drhovecdT_psat(model, T, rhovecL, rhovecV);
+    }, m_model);
+}
+double MI::get_dpsat_dTsat_isopleth(const double T, const REArrayd& rhovecL, const REArrayd& rhovecV) const {
+    return std::visit([&](const auto& model) {
+        return teqp::get_dpsat_dTsat_isopleth(model, T, rhovecL, rhovecV);
+    }, m_model);
+}
+
+nlohmann::json MI::trace_VLE_isotherm_binary(const double T0, const EArrayd& rhovecL0, const EArrayd& rhovecV0, const std::optional<TVLEOptions> &options) const{
+    return std::visit([&](const auto& model) {
+        return teqp::trace_VLE_isotherm_binary(model, T0, rhovecL0, rhovecV0, options);
+    }, m_model);
+}
+nlohmann::json MI::trace_VLE_isobar_binary(const double p, const double T0, const EArrayd& rhovecL0, const EArrayd& rhovecV0, const std::optional<PVLEOptions> &options) const{
+    return std::visit([&](const auto& model) {
+        return teqp::trace_VLE_isobar_binary(model, p, T0, rhovecL0, rhovecV0, options);
+    }, m_model);
+}
+std::tuple<VLE_return_code,EArrayd,EArrayd> MI::mix_VLE_Tx(const double T, const REArrayd& rhovecL0, const REArrayd& rhovecV0, const REArrayd& xspec, const double atol, const double reltol, const double axtol, const double relxtol, const int maxiter) const{
+    return std::visit([&](const auto& model) {
+        return teqp::mix_VLE_Tx(model, T, rhovecL0, rhovecV0, xspec, atol, reltol, axtol, relxtol, maxiter);
+    }, m_model);
+}
+MixVLEReturn MI::mix_VLE_Tp(const double T, const double pgiven, const REArrayd& rhovecL0, const REArrayd& rhovecV0, const std::optional<MixVLETpFlags> &flags) const{
+    return std::visit([&](const auto& model) {
+        return teqp::mix_VLE_Tp(model, T, pgiven, rhovecL0, rhovecV0, flags);
+    }, m_model);
+}
+std::tuple<VLE_return_code,double,EArrayd,EArrayd> MI::mixture_VLE_px(const double p_spec, const REArrayd& xmolar_spec, const double T0, const REArrayd& rhovecL0, const REArrayd& rhovecV0, const std::optional<MixVLEpxFlags>& flags ) const{
+    return std::visit([&](const auto& model) {
+        return teqp::mixture_VLE_px(model, p_spec, xmolar_spec, T0, rhovecL0, rhovecV0, flags);
+    }, m_model);
+}
+EArray2 MI::pure_VLE_T(const double T, const double rhoL, const double rhoV, int maxiter) const {
+    return std::visit([&](const auto& model) {
+        return teqp::pure_VLE_T(model, T, rhoL, rhoV, maxiter);
+    }, m_model);
+}
diff --git a/interface/CPP/teqp_impl_VLLE.cpp b/interface/CPP/teqp_impl_VLLE.cpp
new file mode 100644
index 0000000..eae5696
--- /dev/null
+++ b/interface/CPP/teqp_impl_VLLE.cpp
@@ -0,0 +1,17 @@
+#include "teqpcpp.cpp"
+#include "teqp/algorithms/VLLE.hpp"
+
+using namespace teqp;
+using MI = teqp::cppinterface::ModelImplementer;
+
+std::tuple<VLLE::VLLE_return_code,EArrayd,EArrayd,EArrayd> MI::mix_VLLE_T(const double T, const REArrayd& rhovecVinit, const REArrayd& rhovecL1init, const REArrayd& rhovecL2init, const double atol, const double reltol, const double axtol, const double relxtol, const int maxiter) const{
+    return std::visit([&](const auto& model) {
+        return VLLE::mix_VLLE_T(model, T, rhovecVinit, rhovecL1init, rhovecL2init, atol, reltol, axtol, relxtol, maxiter);
+    }, m_model);
+}
+
+std::vector<nlohmann::json> MI::find_VLLE_T_binary(const std::vector<nlohmann::json>& traces, const std::optional<VLLE::VLLEFinderOptions> options) const{
+    return std::visit([&](const auto& model) {
+        return VLLE::find_VLLE_T_binary(model, traces, options);
+    }, m_model);
+}
diff --git a/interface/CPP/teqp_impl_crit.cpp b/interface/CPP/teqp_impl_crit.cpp
new file mode 100644
index 0000000..8f43d96
--- /dev/null
+++ b/interface/CPP/teqp_impl_crit.cpp
@@ -0,0 +1,61 @@
+#include "teqp/algorithms/critical_tracing.hpp"
+
+#include "teqpcpp.cpp"
+
+using namespace teqp;
+using MI = teqp::cppinterface::ModelImplementer;
+
+nlohmann::json MI::trace_critical_arclength_binary(const double T0, const EArrayd& rhovec0, const std::optional<std::string>& filename, const std::optional<TCABOptions> &options) const {
+    return std::visit([&](const auto& model) {
+        using crit = teqp::CriticalTracing<decltype(model), double, std::decay_t<decltype(rhovec0)>>;
+        return crit::trace_critical_arclength_binary(model, T0, rhovec0, filename , options);
+    }, m_model);
+}
+EArrayd MI::get_drhovec_dT_crit(const double T, const REArrayd& rhovec) const {
+    return std::visit([&](const auto& model) {
+        using crit = teqp::CriticalTracing<decltype(model), double, std::decay_t<decltype(rhovec)>>;
+        return crit::get_drhovec_dT_crit(model, T, rhovec);
+    }, m_model);
+}
+double MI::get_dp_dT_crit(const double T, const REArrayd& rhovec) const {
+    return std::visit([&](const auto& model) {
+        using crit = teqp::CriticalTracing<decltype(model), double, std::decay_t<decltype(rhovec)>>;
+        return crit::get_dp_dT_crit(model, T, rhovec);
+    }, m_model);
+}
+EArray2 MI::get_criticality_conditions(const double T, const REArrayd& rhovec) const {
+    return std::visit([&](const auto& model) {
+        using crit = teqp::CriticalTracing<decltype(model), double, std::decay_t<decltype(rhovec)>>;
+        return crit::get_criticality_conditions(model, T, rhovec);
+    }, m_model);
+}
+EigenData MI::eigen_problem(const double T, const REArrayd& rhovec, const std::optional<REArrayd>& alignment_v0) const {
+    return std::visit([&](const auto& model) {
+        using crit = teqp::CriticalTracing<decltype(model), double, std::decay_t<decltype(rhovec)>>;
+        return crit::eigen_problem(model, T, rhovec, alignment_v0.value_or(Eigen::ArrayXd()));
+    }, m_model);
+}
+double MI::get_minimum_eigenvalue_Psi_Hessian(const double T, const REArrayd& rhovec) const {
+    return std::visit([&](const auto& model) {
+        using crit = teqp::CriticalTracing<decltype(model), double, std::decay_t<decltype(rhovec)>>;
+        return crit::get_minimum_eigenvalue_Psi_Hessian(model, T, rhovec);
+    }, m_model);
+}
+
+
+std::tuple<double, double> MI::solve_pure_critical(const double T, const double rho, const std::optional<nlohmann::json>& flags) const  {
+    return std::visit([&](const auto& model) {
+        return teqp::solve_pure_critical(model, T, rho, flags.value_or(nlohmann::json{}));
+    }, m_model);
+}
+std::tuple<EArrayd, EMatrixd> MI::get_pure_critical_conditions_Jacobian(const double T, const double rho, int alternative_pure_index, int alternative_length) const {
+    return std::visit([&](const auto& model) {
+        return teqp::get_pure_critical_conditions_Jacobian(model, T, rho, alternative_pure_index, alternative_length);
+    }, m_model);
+}
+std::tuple<double, double> MI::extrapolate_from_critical(const double Tc, const double rhoc, const double Tnew) const {
+    return std::visit([&](const auto& model) {
+        auto mat = teqp::extrapolate_from_critical(model, Tc, rhoc, Tnew);
+        return std::make_tuple(mat[0], mat[1]);
+    }, m_model);
+}
diff --git a/interface/CPP/teqp_impl_factory.cpp b/interface/CPP/teqp_impl_factory.cpp
new file mode 100644
index 0000000..ea2de8e
--- /dev/null
+++ b/interface/CPP/teqp_impl_factory.cpp
@@ -0,0 +1,15 @@
+#include "teqpcpp.cpp"
+
+using namespace teqp;
+using namespace teqp::cppinterface;
+
+std::shared_ptr<AbstractModel> teqp::cppinterface::make_model(const nlohmann::json& j) {
+    return std::make_shared<ModelImplementer>(build_model(j));
+}
+std::shared_ptr<AbstractModel> teqp::cppinterface::make_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) {
+    return std::make_shared<ModelImplementer>(build_multifluid_model(components, coolprop_root, BIPcollectionpath, flags, departurepath));
+}
+
+std::shared_ptr<AbstractModel> teqp::cppinterface::emplace_model(AllowedModels&& model){
+    return std::make_shared<ModelImplementer>(std::move(model));
+}
diff --git a/interface/CPP/teqp_impl_isochoric.cpp b/interface/CPP/teqp_impl_isochoric.cpp
new file mode 100644
index 0000000..0419937
--- /dev/null
+++ b/interface/CPP/teqp_impl_isochoric.cpp
@@ -0,0 +1,27 @@
+#include "teqp/derivs.hpp"
+
+#include "teqpcpp.cpp"
+using MI = teqp::cppinterface::ModelImplementer;
+
+// Here XMacros are used to create functions like get_Ar00, get_Ar01, ....
+#define X(i,j) \
+double MI::get_Ar ## i ## j(const double T, const double rho, const REArrayd& molefracs) const { \
+    return std::visit([&](const auto& model) { \
+        using tdx = teqp::TDXDerivatives<decltype(model), double, REArrayd>; \
+        return tdx::template get_Arxy<i,j>(model, T, rho, molefracs); \
+    }, m_model); \
+}
+ARXY_args
+#undef X
+
+// Here XMacros are used to create functions like get_Ar01n, get_Ar02n, ....
+#define X(i) \
+EArrayd MI::get_Ar0 ## i ## n(const double T, const double rho, const REArrayd& molefracs) const { \
+    return std::visit([&](const auto& model) { \
+        using tdx = teqp::TDXDerivatives<decltype(model), double, REArrayd>; \
+        auto vals = tdx::template get_Ar0n<i>(model, T, rho, molefracs); \
+        return Eigen::Map<Eigen::ArrayXd>(&(vals[0]), vals.size());\
+    }, m_model); \
+}
+AR0N_args
+#undef X
diff --git a/interface/CPP/teqp_impl_virials.cpp b/interface/CPP/teqp_impl_virials.cpp
new file mode 100644
index 0000000..35aaab0
--- /dev/null
+++ b/interface/CPP/teqp_impl_virials.cpp
@@ -0,0 +1,29 @@
+#include "teqp/derivs.hpp"
+#include "teqpcpp.cpp"
+
+using MI = teqp::cppinterface::ModelImplementer;
+
+double MI::get_B2vir(const double T, const EArrayd& molefrac) const  {
+    return std::visit([&](const auto& model) {
+        using vd = VirialDerivatives<decltype(model), double, RAX>;
+        return vd::get_B2vir(model, T, molefrac);
+    }, m_model);
+}
+double MI::get_B12vir(const double T, const EArrayd& molefrac) const {
+    return std::visit([&](const auto& model) {
+        using vd = VirialDerivatives<decltype(model), double, RAX>;
+        return vd::get_B12vir(model, T, molefrac);
+    }, m_model);
+}
+std::map<int, double> MI::get_Bnvir(const int Nderiv, const double T, const EArrayd& molefrac) const {
+    return std::visit([&](const auto& model) {
+        using vd = VirialDerivatives<decltype(model), double, RAX>;
+        return vd::get_Bnvir_runtime(Nderiv, model, T, molefrac);
+    }, m_model);
+}
+double MI::get_dmBnvirdTm(const int Nderiv, const int NTderiv, const double T, const EArrayd& molefrac) const {
+    return std::visit([&](const auto& model) {
+        using vd = VirialDerivatives<decltype(model), double, RAX>;
+        return vd::get_dmBnvirdTm_runtime(Nderiv, NTderiv, model, T, molefrac);
+    }, m_model);
+}
diff --git a/interface/CPP/teqpcpp.cpp b/interface/CPP/teqpcpp.cpp
index 30e5a8a..d23eaa9 100644
--- a/interface/CPP/teqpcpp.cpp
+++ b/interface/CPP/teqpcpp.cpp
@@ -1,12 +1,5 @@
 #include "teqp/cpp/teqpcpp.hpp"
-#include "teqp/derivs.hpp"
 #include "teqp/json_builder.hpp"
-#include "teqp/algorithms/critical_tracing.hpp"
-#include "teqp/algorithms/VLE.hpp"
-#include "teqp/algorithms/VLLE.hpp"
-
-using namespace teqp;
-using namespace teqp::cppinterface;
 
 namespace teqp {
     namespace cppinterface {
@@ -42,230 +35,70 @@ namespace teqp {
                 }, m_model);
             }
             
-            nlohmann::json trace_critical_arclength_binary(const double T0, const EArrayd& rhovec0, const std::optional<std::string>& filename, const std::optional<TCABOptions> &options) const override {
-                return std::visit([&](const auto& model) {
-                    using crit = teqp::CriticalTracing<decltype(model), double, std::decay_t<decltype(rhovec0)>>;
-                    return crit::trace_critical_arclength_binary(model, T0, rhovec0, filename , options);
-                }, m_model);
-            }
-            EArrayd get_drhovec_dT_crit(const double T, const REArrayd& rhovec) const override {
-                return std::visit([&](const auto& model) {
-                    using crit = teqp::CriticalTracing<decltype(model), double, std::decay_t<decltype(rhovec)>>;
-                    return crit::get_drhovec_dT_crit(model, T, rhovec);
-                }, m_model);
-            }
-            double get_dp_dT_crit(const double T, const REArrayd& rhovec) const override {
-                return std::visit([&](const auto& model) {
-                    using crit = teqp::CriticalTracing<decltype(model), double, std::decay_t<decltype(rhovec)>>;
-                    return crit::get_dp_dT_crit(model, T, rhovec);
-                }, m_model);
-            }
-            EArray2 get_criticality_conditions(const double T, const REArrayd& rhovec) const override {
-                return std::visit([&](const auto& model) {
-                    using crit = teqp::CriticalTracing<decltype(model), double, std::decay_t<decltype(rhovec)>>;
-                    return crit::get_criticality_conditions(model, T, rhovec);
-                }, m_model);
-            }
-            EigenData eigen_problem(const double T, const REArrayd& rhovec, const std::optional<REArrayd>& alignment_v0) const override {
-                return std::visit([&](const auto& model) {
-                    using crit = teqp::CriticalTracing<decltype(model), double, std::decay_t<decltype(rhovec)>>;
-                    return crit::eigen_problem(model, T, rhovec, alignment_v0.value_or(Eigen::ArrayXd()));
-                }, m_model);
-            }
-            double get_minimum_eigenvalue_Psi_Hessian(const double T, const REArrayd& rhovec) const override {
-                return std::visit([&](const auto& model) {
-                    using crit = teqp::CriticalTracing<decltype(model), double, std::decay_t<decltype(rhovec)>>;
-                    return crit::get_minimum_eigenvalue_Psi_Hessian(model, T, rhovec);
-                }, m_model);
-            }
-            
-            EArray2 pure_VLE_T(const double T, const double rhoL, const double rhoV, int maxiter) const override {
-                return std::visit([&](const auto& model) {
-                    return teqp::pure_VLE_T(model, T, rhoL, rhoV, maxiter);
-                }, m_model);
-            }
-            std::tuple<double, double> solve_pure_critical(const double T, const double rho, const std::optional<nlohmann::json>& flags) const override {
-                return std::visit([&](const auto& model) {
-                    return teqp::solve_pure_critical(model, T, rho, flags.value_or(nlohmann::json{}));
-                }, m_model);
-            }
-            std::tuple<EArrayd, EMatrixd> get_pure_critical_conditions_Jacobian(const double T, const double rho, int alternative_pure_index, int alternative_length) const override {
-                return std::visit([&](const auto& model) {
-                    return teqp::get_pure_critical_conditions_Jacobian(model, T, rho, alternative_pure_index, alternative_length);
-                }, m_model);
-            }
-            std::tuple<double, double> extrapolate_from_critical(const double Tc, const double rhoc, const double Tnew) const override {
-                return std::visit([&](const auto& model) {
-                    auto mat = teqp::extrapolate_from_critical(model, Tc, rhoc, Tnew);
-                    return std::make_tuple(mat[0], mat[1]);
-                }, m_model);
-            }
+            nlohmann::json trace_critical_arclength_binary(const double T0, const EArrayd& rhovec0, const std::optional<std::string>& filename, const std::optional<TCABOptions> &options) const override;
+            EArrayd get_drhovec_dT_crit(const double T, const REArrayd& rhovec) const override;
+            double get_dp_dT_crit(const double T, const REArrayd& rhovec) const override;
+            EArray2 get_criticality_conditions(const double T, const REArrayd& rhovec) const override;
+            EigenData eigen_problem(const double T, const REArrayd& rhovec, const std::optional<REArrayd>& alignment_v0) const override;
+            double get_minimum_eigenvalue_Psi_Hessian(const double T, const REArrayd& rhovec) const override ;
             
+            std::tuple<double, double> solve_pure_critical(const double T, const double rho, const std::optional<nlohmann::json>& flags) const override ;
+            std::tuple<EArrayd, EMatrixd> get_pure_critical_conditions_Jacobian(const double T, const double rho, int alternative_pure_index, int alternative_length) const override ;
+            std::tuple<double, double> extrapolate_from_critical(const double Tc, const double rhoc, const double Tnew) const override ;
+
             // Derivatives from isochoric thermodynamics (all have the same signature)
             #define X(f) \
-            virtual double f(const double T, const REArrayd& rhovec) const override { \
-                return std::visit([&](const auto& model) { \
-                    using id = IsochoricDerivatives<decltype(model), double, REArrayd>; \
-                    return id::f(model, T, rhovec); \
-                }, m_model); \
-            }
+            double f(const double T, const REArrayd& rhovec) const override ;
             ISOCHORIC_double_args
             #undef X
-            
+
             #define X(f) \
-            virtual EArrayd f(const double T, const REArrayd& rhovec) const override { \
-                return std::visit([&](const auto& model) { \
-                    using id = IsochoricDerivatives<decltype(model), double, REArrayd>; \
-                    return id::f(model, T, rhovec); \
-                }, m_model); \
-            }
+            EArrayd f(const double T, const REArrayd& rhovec) const override ;
             ISOCHORIC_array_args
             #undef X
-            
+
             #define X(f) \
-            virtual EMatrixd f(const double T, const REArrayd& rhovec) const override { \
-                return std::visit([&](const auto& model) { \
-                    using id = IsochoricDerivatives<decltype(model), double, REArrayd>; \
-                    return id::f(model, T, rhovec); \
-                }, m_model); \
-            }
+            EMatrixd f(const double T, const REArrayd& rhovec) const override ;
             ISOCHORIC_matrix_args
             #undef X
+
+            EArray33d get_deriv_mat2(const double T, double rho, const EArrayd& z) const override ;
             
-            EArray33d get_deriv_mat2(const double T, double rho, const EArrayd& z) const override {
-                return std::visit([&](const auto& model) {
-                    // Although the template argument suggests that only residual terms
-                    // are returned, also the ideal-gas ones are returned because the
-                    // ideal-gas term is required to implement alphar which just redirects
-                    // to alphaig
-                    return DerivativeHolderSquare<2, AlphaWrapperOption::residual>(model, T, rho, z).derivs;
-                }, m_model);
-            }
-            double get_B2vir(const double T, const EArrayd& molefrac) const override {
-                return std::visit([&](const auto& model) {
-                    using vd = VirialDerivatives<decltype(model), double, RAX>;
-                    return vd::get_B2vir(model, T, molefrac);
-                }, m_model);
-            }
-            double get_B12vir(const double T, const EArrayd& molefrac) const override {
-                return std::visit([&](const auto& model) {
-                    using vd = VirialDerivatives<decltype(model), double, RAX>;
-                    return vd::get_B12vir(model, T, molefrac);
-                }, m_model);
-            }
-            std::map<int, double> get_Bnvir(const int Nderiv, const double T, const EArrayd& molefrac) const override {
-                return std::visit([&](const auto& model) {
-                    using vd = VirialDerivatives<decltype(model), double, RAX>;
-                    return vd::get_Bnvir_runtime(Nderiv, model, T, molefrac);
-                }, m_model);
-            }
-            double get_dmBnvirdTm(const int Nderiv, const int NTderiv, const double T, const EArrayd& molefrac) const override {
-                return std::visit([&](const auto& model) {
-                    using vd = VirialDerivatives<decltype(model), double, RAX>;
-                    return vd::get_dmBnvirdTm_runtime(Nderiv, NTderiv, model, T, molefrac);
-                }, m_model);
-            }
-            
-            double get_Arxy(const int NT, const int ND, const double T, const double rho, const EArrayd& molefracs) const override {
-                return std::visit([&](const auto& model) {
-                    using tdx = teqp::TDXDerivatives<decltype(model), double, EArrayd>;
-                    return tdx::template get_Ar(NT, ND, model, T, rho, molefracs);
-                }, m_model);
-            }
+            double get_B2vir(const double T, const EArrayd& molefrac) const override;
+            double get_B12vir(const double T, const EArrayd& molefrac) const override;
+            std::map<int, double> get_Bnvir(const int Nderiv, const double T, const EArrayd& molefrac) const override;
+            double get_dmBnvirdTm(const int Nderiv, const int NTderiv, const double T, const EArrayd& molefrac) const override;
+
+            double get_Arxy(const int NT, const int ND, const double T, const double rho, const EArrayd& molefracs) const override;
             // Here XMacros are used to create functions like get_Ar00, get_Ar01, ....
             #define X(i,j) \
-            double get_Ar ## i ## j(const double T, const double rho, const REArrayd& molefracs) const override { \
-                return std::visit([&](const auto& model) { \
-                    using tdx = teqp::TDXDerivatives<decltype(model), double, REArrayd>; \
-                    return tdx::template get_Arxy<i,j>(model, T, rho, molefracs); \
-                }, m_model); \
-            }
+            double get_Ar ## i ## j(const double T, const double rho, const REArrayd& molefracs) const override;
             ARXY_args
             #undef X
-            
             // Here XMacros are used to create functions like get_Ar01n, get_Ar02n, ....
             #define X(i) \
-            EArrayd get_Ar0 ## i ## n(const double T, const double rho, const REArrayd& molefracs) const override { \
-                return std::visit([&](const auto& model) { \
-                    using tdx = teqp::TDXDerivatives<decltype(model), double, REArrayd>; \
-                    auto vals = tdx::template get_Ar0n<i>(model, T, rho, molefracs); \
-                    return Eigen::Map<Eigen::ArrayXd>(&(vals[0]), vals.size());\
-                }, m_model); \
-            }
+            EArrayd get_Ar0 ## i ## n(const double T, const double rho, const REArrayd& molefracs) const override;
             AR0N_args
             #undef X
-            
-            double get_neff(const double T, const double rho, const EArrayd& molefracs) const override {
-                return std::visit([&](const auto& model) {
-                    using tdx = teqp::TDXDerivatives<decltype(model), double, EArrayd>;
-                    return tdx::template get_neff(model, T, rho, molefracs);
-                }, m_model);
-            }
-            
-            std::tuple<EArrayd, EArrayd> get_drhovecdp_Tsat(const double T, const REArrayd& rhovecL, const REArrayd& rhovecV) const override {
-                return std::visit([&](const auto& model) {
-                    return teqp::get_drhovecdp_Tsat(model, T, rhovecL, rhovecV);
-                }, m_model);
-            }
-            std::tuple<EArrayd, EArrayd> get_drhovecdT_psat(const double T, const REArrayd& rhovecL, const REArrayd& rhovecV) const override {
-                return std::visit([&](const auto& model) {
-                    return teqp::get_drhovecdT_psat(model, T, rhovecL, rhovecV);
-                }, m_model);
-            }
-            double get_dpsat_dTsat_isopleth(const double T, const REArrayd& rhovecL, const REArrayd& rhovecV) const override {
-                return std::visit([&](const auto& model) {
-                    return teqp::get_dpsat_dTsat_isopleth(model, T, rhovecL, rhovecV);
-                }, m_model);
-            }
-            
-            nlohmann::json trace_VLE_isotherm_binary(const double T0, const EArrayd& rhovecL0, const EArrayd& rhovecV0, const std::optional<TVLEOptions> &options) const override{
-                return std::visit([&](const auto& model) {
-                    return teqp::trace_VLE_isotherm_binary(model, T0, rhovecL0, rhovecV0, options);
-                }, m_model);
-            }
-            nlohmann::json trace_VLE_isobar_binary(const double p, const double T0, const EArrayd& rhovecL0, const EArrayd& rhovecV0, const std::optional<PVLEOptions> &options) const override{
-                return std::visit([&](const auto& model) {
-                    return teqp::trace_VLE_isobar_binary(model, p, T0, rhovecL0, rhovecV0, options);
-                }, m_model);
-            }
-            std::tuple<VLE_return_code,EArrayd,EArrayd> mix_VLE_Tx(const double T, const REArrayd& rhovecL0, const REArrayd& rhovecV0, const REArrayd& xspec, const double atol, const double reltol, const double axtol, const double relxtol, const int maxiter) const override{
-                return std::visit([&](const auto& model) {
-                    return teqp::mix_VLE_Tx(model, T, rhovecL0, rhovecV0, xspec, atol, reltol, axtol, relxtol, maxiter);
-                }, m_model);
-            }
-            MixVLEReturn mix_VLE_Tp(const double T, const double pgiven, const REArrayd& rhovecL0, const REArrayd& rhovecV0, const std::optional<MixVLETpFlags> &flags) const override{
-                return std::visit([&](const auto& model) {
-                    return teqp::mix_VLE_Tp(model, T, pgiven, rhovecL0, rhovecV0, flags);
-                }, m_model);
-            }
-            std::tuple<VLE_return_code,double,EArrayd,EArrayd> mixture_VLE_px(const double p_spec, const REArrayd& xmolar_spec, const double T0, const REArrayd& rhovecL0, const REArrayd& rhovecV0, const std::optional<MixVLEpxFlags>& flags ) const override{
-                return std::visit([&](const auto& model) {
-                    return teqp::mixture_VLE_px(model, p_spec, xmolar_spec, T0, rhovecL0, rhovecV0, flags);
-                }, m_model);
-            }
-            
-            std::tuple<VLLE::VLLE_return_code,EArrayd,EArrayd,EArrayd> mix_VLLE_T(const double T, const REArrayd& rhovecVinit, const REArrayd& rhovecL1init, const REArrayd& rhovecL2init, const double atol, const double reltol, const double axtol, const double relxtol, const int maxiter) const override{
-                return std::visit([&](const auto& model) {
-                    return teqp::VLLE::mix_VLLE_T(model, T, rhovecVinit, rhovecL1init, rhovecL2init, atol, reltol, axtol, relxtol, maxiter);
-                }, m_model);
-            }
-            
-            std::vector<nlohmann::json> find_VLLE_T_binary(const std::vector<nlohmann::json>& traces, const std::optional<VLLE::VLLEFinderOptions> options) const override{
-                return std::visit([&](const auto& model) {
-                    return teqp::VLLE::find_VLLE_T_binary(model, traces, options);
-                }, m_model);
-            }
+
+            double get_neff(const double T, const double rho, const EArrayd& molefracs) const override;
+
+            EArray2 pure_VLE_T(const double T, const double rhoL, const double rhoV, int maxiter) const override;
+            std::tuple<EArrayd, EArrayd> get_drhovecdp_Tsat(const double T, const REArrayd& rhovecL, const REArrayd& rhovecV) const override;
+            std::tuple<EArrayd, EArrayd> get_drhovecdT_psat(const double T, const REArrayd& rhovecL, const REArrayd& rhovecV) const override;
+            double get_dpsat_dTsat_isopleth(const double T, const REArrayd& rhovecL, const REArrayd& rhovecV) const override;
+
+            nlohmann::json trace_VLE_isotherm_binary(const double T0, const EArrayd& rhovecL0, const EArrayd& rhovecV0, const std::optional<TVLEOptions> &options) const override;
+            nlohmann::json trace_VLE_isobar_binary(const double p, const double T0, const EArrayd& rhovecL0, const EArrayd& rhovecV0, const std::optional<PVLEOptions> &options) const override;
+            std::tuple<VLE_return_code,EArrayd,EArrayd> mix_VLE_Tx(const double T, const REArrayd& rhovecL0, const REArrayd& rhovecV0, const REArrayd& xspec, const double atol, const double reltol, const double axtol, const double relxtol, const int maxiter) const override;
+            MixVLEReturn mix_VLE_Tp(const double T, const double pgiven, const REArrayd& rhovecL0, const REArrayd& rhovecV0, const std::optional<MixVLETpFlags> &flags) const override;
+            std::tuple<VLE_return_code,double,EArrayd,EArrayd> mixture_VLE_px(const double p_spec, const REArrayd& xmolar_spec, const double T0, const REArrayd& rhovecL0, const REArrayd& rhovecV0, const std::optional<MixVLEpxFlags>& flags ) const override;
+
+            std::tuple<VLLE::VLLE_return_code,EArrayd,EArrayd,EArrayd> mix_VLLE_T(const double T, const REArrayd& rhovecVinit, const REArrayd& rhovecL1init, const REArrayd& rhovecL2init, const double atol, const double reltol, const double axtol, const double relxtol, const int maxiter) const override;
+
+            std::vector<nlohmann::json> find_VLLE_T_binary(const std::vector<nlohmann::json>& traces, const std::optional<VLLE::VLLEFinderOptions> options) const override;
         };
 
-        std::shared_ptr<AbstractModel> make_model(const nlohmann::json& j) {
-            return std::make_shared<ModelImplementer>(build_model(j));
-        }
-        std::shared_ptr<AbstractModel> make_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) {
-            return std::make_shared<ModelImplementer>(build_multifluid_model(components, coolprop_root, BIPcollectionpath, flags, departurepath));
-        }
-    
-        std::shared_ptr<AbstractModel> emplace_model(AllowedModels&& model){
-            return std::make_shared<ModelImplementer>(std::move(model));
-        }
+        
     }
 }
diff --git a/interface/CPP/bench_teqpcpp.cpp b/interface/CPP/test/bench_teqpcpp.cpp
similarity index 100%
rename from interface/CPP/bench_teqpcpp.cpp
rename to interface/CPP/test/bench_teqpcpp.cpp
diff --git a/interface/CPP/test_teqpcpp.cpp b/interface/CPP/test/test_teqpcpp.cpp
similarity index 100%
rename from interface/CPP/test_teqpcpp.cpp
rename to interface/CPP/test/test_teqpcpp.cpp
-- 
GitLab