More generic derivative functions

include/teqp/derivs.hpp

@@ -6,37 +6,29 @@ caller(const FuncType& f, TType T, const ContainerType& rho) {
     return f(T, rho);
 }
 
-/// Given a 
+/// Given a function, use complex step derivatives to calculate the derivative with respect to the first variable
+/// which here is temperature
 template <typename TType, typename ContainerType, typename FuncType>
 typename std::enable_if<is_container<ContainerType>::value, typename ContainerType::value_type>::type
-deriv1(const FuncType& f, TType T, const ContainerType& rho) {
+derivT(const FuncType& f, TType T, const ContainerType& rho) {
     double h = 1e-100;
     return f(std::complex<TType>(T, h), rho).imag() / h;
 }
 
-template <typename TType, typename ContainerType, typename FuncType>
+/// Given a function, use complex step derivatives to calculate the derivative with respect to the first composition variable
+template <typename TType, typename ContainerType, typename FuncType, typename Integer>
 typename std::enable_if<is_container<ContainerType>::value, typename ContainerType::value_type>::type
-deriv2(const FuncType& f, TType T, const ContainerType& rho) {
+derivrhoi(const FuncType& f, TType T, const ContainerType& rho, Integer j) {
     double h = 1e-100;
     using comtype = std::complex<ContainerType::value_type>;
     std::valarray<comtype> rhocom(rho.size());
-    rhocom[0] = comtype(rho[0], h);
-    for (auto i = 1; i < rho.size(); ++i) {
+    for (auto i = 0; i < rho.size(); ++i) {
         rhocom[i] = comtype(rho[i], 0.0);
     }
+    rhocom[j] = comtype(rho[j], h);
     return f(T, rhocom).imag() / h;
 }
 
-template <typename TType, typename ContainerType, typename FuncType>
-typename std::enable_if<is_container<ContainerType>::value, typename ContainerType::value_type>::type
-deriv3(const FuncType& f, TType T, const ContainerType& rho) {
-    double h = 1e-100;
-    using comtype = std::complex<ContainerType::value_type>;
-    std::valarray<comtype> rhocom(rho.size());
-    rhocom[0] = comtype(rho[0], 0.0);
-    rhocom[1] = comtype(rho[1], h);
-    return f(T, rhocom).imag() / h;
-}
 
 template<typename Model, typename TType, typename RhoType>
 auto build_Psir_Hessian(const Model& model, const TType T, const RhoType& rho) {

src/main.cpp

@@ -33,8 +33,8 @@ void test_vdW() {
     auto t21 = std::chrono::steady_clock::now();
     
     auto Psir = vdW.Psir(T, rhovec);
-    auto dPsirdrho0 = rhovec[0] * deriv2([&vdW](const auto& T, const auto& rhovec) { return vdW.Psir(T, rhovec); }, T, rhovec);
-    auto dPsirdrho1 = rhovec[1] * deriv3([&vdW](const auto& T, const auto& rhovec) { return vdW.Psir(T, rhovec); }, T, rhovec);
+    auto dPsirdrho0 = rhovec[0] * derivrhoi([&vdW](const auto& T, const auto& rhovec) { return vdW.Psir(T, rhovec); }, T, rhovec, 0);
+    auto dPsirdrho1 = rhovec[1] * derivrhoi([&vdW](const auto& T, const auto& rhovec) { return vdW.Psir(T, rhovec); }, T, rhovec, 1);
     auto pfromderiv = rho * R * T - Psir + dPsirdrho0 + dPsirdrho1;
 
     auto t31 = std::chrono::steady_clock::now();
@@ -58,12 +58,12 @@ void test_vdwMix() {
     auto t2 = std::chrono::steady_clock::now();
 
     auto Psir = vdW.Psir(T, rhovec);
-    auto dPsirdrho0 = rhovec[0] * deriv2([&vdW, rhotot](const auto& T, const auto& rhovec) { return vdW.Psir(T, rhovec); }, T, rhovec);
-    auto dPsirdrho1 = rhovec[1] * deriv3([&vdW, rhotot](const auto& T, const auto& rhovec) { return vdW.Psir(T, rhovec); }, T, rhovec);
+    auto dPsirdrho0 = rhovec[0] * derivrhoi([&vdW, rhotot](const auto& T, const auto& rhovec) { return vdW.Psir(T, rhovec); }, T, rhovec, 0);
+    auto dPsirdrho1 = rhovec[1] * derivrhoi([&vdW, rhotot](const auto& T, const auto& rhovec) { return vdW.Psir(T, rhovec); }, T, rhovec, 1);
     auto pfromderiv = rho*R*T - Psir + dPsirdrho0 + dPsirdrho1;
     {
-        auto term0 = rhovec[0] * deriv2([&vdW, rhotot](const auto& T, const auto& rhovec) { return vdW.alphar(T, rhovec); }, T, rhovec);
-        auto term1 = rhovec[1] * deriv3([&vdW, rhotot](const auto& T, const auto& rhovec) { return vdW.alphar(T, rhovec); }, T, rhovec);
+        auto term0 = rhovec[0] * derivrhoi([&vdW, rhotot](const auto& T, const auto& rhovec) { return vdW.alphar(T, rhovec); }, T, rhovec, 0);
+        auto term1 = rhovec[1] * derivrhoi([&vdW, rhotot](const auto& T, const auto& rhovec) { return vdW.alphar(T, rhovec); }, T, rhovec, 1);
         auto pr = (term0 + term1)*rhotot*R*T;
         auto pfromderiv2 = rho*R*T + pr;
         auto err2 = pfromderiv / pfromderiv2 - 1;