Repository: libKriging/libKriging (language: python)

Issue:
Add memoization feature to minimize re-computations
Cf Kriging.cpp (1029) bool bfgs_ok = optim::lbfgs.

What parameters should be tracked ?


Discussion / hints:

Fix patch:
diff --git a/src/lib/CMakeLists.txt b/src/lib/CMakeLists.txt
index f274fdfb..5896d48d 100644
--- a/src/lib/CMakeLists.txt
+++ b/src/lib/CMakeLists.txt
@@ -10,7 +10,9 @@ add_library(Kriging
         Kriging.cpp include/libKriging/Kriging.hpp
         include/libKriging/KrigingException.hpp
         lkalloc.cpp include/libKriging/utils/lkalloc.hpp
-        Bench.cpp include/libKriging/Bench.hpp)
+        Bench.cpp include/libKriging/Bench.hpp
+        CacheFunction.cpp include/libKriging/CacheFunction.hpp include/libKriging/utils/custom_hash_function.hpp include/libKriging/utils/cache_details.hpp
+        include/libKriging/utils/LinearHashStorage.hpp)
 
 if (CXX_CLANG_TIDY)
     set_target_properties(Kriging
@@ -28,7 +30,7 @@ target_link_libraries(Kriging PUBLIC armadillo)
 #   set(CMAKE_WINDOWS_EXPORT_ALL_SYMBOLS ON)
 
 # https://cmake.org/cmake/help/v3.11/module/GenerateExportHeader.html
-include (GenerateExportHeader)
+include(GenerateExportHeader)
 GENERATE_EXPORT_HEADER(Kriging # generates the export header `lib`_EXPORTS.h automatically
         BASE_NAME LIBKRIGING
         EXPORT_FILE_NAME libKriging/libKriging_exports.h)
@@ -51,14 +53,14 @@ set_target_properties(Kriging PROPERTIES CXX_VISIBILITY_PRESET hidden)
 ## target_include_directories (Kriging PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
 target_include_directories(Kriging
         PUBLIC
-            # path of the headers after installation
-            $<INSTALL_INTERFACE:include> # <prefix>/include/libKriging
-            # path of the headers before installation
-            $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>
-            # path of the generated headers before installation
-            $<BUILD_INTERFACE:${CMAKE_CURRENT_BINARY_DIR}>
+        # path of the headers after installation
+        $<INSTALL_INTERFACE:include> # <prefix>/include/libKriging
+        # path of the headers before installation
+        $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>
+        # path of the generated headers before installation
+        $<BUILD_INTERFACE:${CMAKE_CURRENT_BINARY_DIR}>
         PRIVATE
-            ${CMAKE_CURRENT_SOURCE_DIR}
+        ${CMAKE_CURRENT_SOURCE_DIR}
         )
 
 add_subdirectory(include)
diff --git a/src/lib/CacheFunction.cpp b/src/lib/CacheFunction.cpp
new file mode 100644
index 00000000..9f828b1c
--- /dev/null
+++ b/src/lib/CacheFunction.cpp
@@ -0,0 +1,53 @@
+#include "libKriging/CacheFunction.hpp"
+#include <algorithm>
+#include <iomanip>
+#include <iostream>
+#include <numeric>
+#include <vector>
+
+/* ----------------------------------------------------------------------------------------------------------------- */
+
+CacheFunctionCommon::~CacheFunctionCommon() {
+  if (!m_cache_hit.empty())
+    std::cout << stat() << std::endl;
+}
+
+/* ----------------------------------------------------------------------------------------------------------------- */
+
+auto CacheFunctionCommon::stat() -> CacheStat {
+  std::vector<uint32_t> hits;
+  const auto cache_size = m_cache_hit.size();
+  hits.reserve(cache_size);
+  for (auto [_, e] : m_cache_hit) {
+    hits.push_back(e);
+  }
+  const auto min_hit = [&] {
+    const auto min_element = std::min_element(hits.begin(), hits.end());
+    return (min_element == hits.end()) ? 0 : *min_element;
+  }();
+  const auto max_hit = [&] {
+    const auto max_element = std::max_element(hits.begin(), hits.end());
+    return (max_element == hits.end()) ? 0 : *max_element;
+  }();
+  const auto total_hit = std::accumulate(hits.begin(), hits.end(), uint32_t{0});
+  return {min_hit, max_hit, total_hit, cache_size, m_hash_timer / 1, m_lookup_timer / 1, m_eval_timer / 1};
+}
+
+/* ----------------------------------------------------------------------------------------------------------------- */
+
+std::ostream& operator<<(std::ostream& out, const CacheFunctionCommon::CacheStat& st) {
+  std::ostream o(out.rdbuf());  // safe RAII flags restore
+  o << st.total_hit << " hits on cache with " << st.cache_size << " entries\n";
+  o.precision(3);
+  o << "hit average: " << (1.0 * st.total_hit / st.cache_size) << " in range [ " << st.min_hit << " - " << st.max_hit
+    << " ]\n";
+  o.precision(1);
+  o << "hash time   (ns): " << st.hash_time << " (" << (100.0 * st.hash_time / st.eval_time) << "%)\n";
+  o << "lookup time (ns): " << st.lookup_time << " (" << (100.0 * st.lookup_time / st.eval_time) << "%)\n";
+  o << "call time   (ns): " << st.eval_time << "\n";
+  const auto current_cost = st.eval_time + st.lookup_time + st.hash_time;
+  const auto nominal_cost = 1. * st.eval_time / st.cache_size * st.total_hit;
+  o.precision(3);
+  o << "speed up        : " << (nominal_cost / current_cost) << "\n";
+  return out;
+}
diff --git a/src/lib/Kriging.cpp b/src/lib/Kriging.cpp
index 7b88cac5..288eb901 100644
--- a/src/lib/Kriging.cpp
+++ b/src/lib/Kriging.cpp
@@ -5,8 +5,9 @@
 // clang-format on
 
 #include "libKriging/Kriging.hpp"
+#include "libKriging/CacheFunction.hpp"
 #include "libKriging/KrigingException.hpp"
-
+#include "libKriging/utils/custom_hash_function.hpp"
 #include "libKriging/utils/lk_armadillo.hpp"
 
 #include <cassert>
@@ -884,38 +885,41 @@ LIBKRIGING_EXPORT void Kriging::fit(const arma::colvec& y,
   m_optim = optim;
   m_objective = objective;
   if (objective.compare("LL") == 0) {
-    fit_ofn = [this](const arma::vec& _gamma, arma::vec* grad_out, arma::mat* hess_out, Kriging::OKModel* okm_data) {
-      // Change variable for opt: . -> 1/exp(.)
-      arma::vec _theta = 1 / arma::exp(_gamma);
-      double ll = this->logLikelihood(_theta, grad_out, hess_out, okm_data);
-      if (grad_out != nullptr)
-        *grad_out = *grad_out % _theta;
-      if (hess_out != nullptr)
-        *hess_out = -*grad_out + *hess_out % _theta;
-      return -ll;
-    };
+    fit_ofn = CacheFunction{
+        [this](const arma::vec& _gamma, arma::vec* grad_out, arma::mat* hess_out, Kriging::OKModel* okm_data) {
+          // Change variable for opt: . -> 1/exp(.)
+          arma::vec _theta = 1 / arma::exp(_gamma);
+          double ll = this->logLikelihood(_theta, grad_out, hess_out, okm_data);
+          if (grad_out != nullptr)
+            *grad_out = *grad_out % _theta;
+          if (hess_out != nullptr)
+            *hess_out = -*grad_out + *hess_out % _theta;
+          return -ll;
+        }};
 
   } else if (objective.compare("LOO") == 0) {
-    fit_ofn = [this](const arma::vec& _gamma, arma::vec* grad_out, arma::mat* hess_out, Kriging::OKModel* okm_data) {
-      // Change variable for opt: . -> 1/exp(.)
-      arma::vec _theta = 1 / arma::exp(_gamma);
-      double loo = this->leaveOneOut(_theta, grad_out, okm_data);
-      if (grad_out != nullptr)
-        *grad_out = -*grad_out % _theta;
-      return loo;
-    };
+    fit_ofn = CacheFunction{
+        [this](const arma::vec& _gamma, arma::vec* grad_out, arma::mat* hess_out, Kriging::OKModel* okm_data) {
+          // Change variable for opt: . -> 1/exp(.)
+          arma::vec _theta = 1 / arma::exp(_gamma);
+          double loo = this->leaveOneOut(_theta, grad_out, okm_data);
+          if (grad_out != nullptr)
+            *grad_out = -*grad_out % _theta;
+          return loo;
+        }};
 
   } else if (objective.compare("LMP") == 0) {
     // Our impl. of https://github.com/cran/RobustGaSP/blob/5cf21658e6a6e327be6779482b93dfee25d24592/R/rgasp.R#L303
     //@see Mengyang Gu, Xiao-jing Wang and Jim Berger, 2018, Annals of Statistics.
-    fit_ofn = [this](const arma::vec& _gamma, arma::vec* grad_out, arma::mat* hess_out, Kriging::OKModel* okm_data) {
-      // Change variable for opt: . -> 1/exp(.)
-      const arma::vec& _theta = 1 / arma::exp(_gamma);
-      double lmp = this->logMargPost(_theta, grad_out, okm_data);
-      if (grad_out != nullptr)
-        *grad_out = *grad_out % _theta;
-      return -lmp;
-    };
+    fit_ofn = CacheFunction{
+        [this](const arma::vec& _gamma, arma::vec* grad_out, arma::mat* hess_out, Kriging::OKModel* okm_data) {
+          // Change variable for opt: . -> 1/exp(.)
+          const arma::vec& _theta = 1 / arma::exp(_gamma);
+          double lmp = this->logMargPost(_theta, grad_out, okm_data);
+          if (grad_out != nullptr)
+            *grad_out = *grad_out % _theta;
+          return -lmp;
+        }};
 
   } else
     throw std::invalid_argument("Unsupported fit objective: " + objective);
@@ -1040,7 +1044,7 @@ LIBKRIGING_EXPORT void Kriging::fit(const arma::colvec& y,
 
       bool bfgs_ok = optim::lbfgs(
           gamma_tmp,
-          [&okm_data, this, fit_ofn](const arma::vec& vals_inp, arma::vec* grad_out, void*) -> double {
+          [&okm_data, this, &fit_ofn](const arma::vec& vals_inp, arma::vec* grad_out, void*) -> double {
             return fit_ofn(vals_inp, grad_out, nullptr, &okm_data);
           },
           nullptr,
@@ -1093,7 +1097,7 @@ LIBKRIGING_EXPORT void Kriging::fit(const arma::colvec& y,
       Kriging::OKModel okm_data{T, M, z, beta, parameters.estim_beta, sigma2, parameters.estim_sigma2};
 
       double min_ofn_tmp = optim_newton(
-          [&okm_data, this, fit_ofn](const arma::vec& vals_inp, arma::vec* grad_out, arma::mat* hess_out) -> double {
+          [&okm_data, this, &fit_ofn](const arma::vec& vals_inp, arma::vec* grad_out, arma::mat* hess_out) -> double {
             return fit_ofn(vals_inp, grad_out, hess_out, &okm_data);
           },
           gamma_tmp,
diff --git a/src/lib/include/libKriging/CacheFunction.hpp b/src/lib/include/libKriging/CacheFunction.hpp
new file mode 100644
index 00000000..156eaccf
--- /dev/null
+++ b/src/lib/include/libKriging/CacheFunction.hpp
@@ -0,0 +1,130 @@
+#ifndef LIBKRIGING_SRC_LIB_CACHE_HPP
+#define LIBKRIGING_SRC_LIB_CACHE_HPP
+
+#include <cassert>
+#include <iostream>
+#include <tuple>
+#include <unordered_map>
+#include "libKriging/libKriging_exports.h"
+#include "libKriging/utils/LinearHashStorage.hpp"
+#include "libKriging/utils/cache_details.hpp"
+
+#define LIBKRIGING_CACHE_ANALYSE  // Use timers to measure part of the cache processing
+#define LIBKRIGING_CACHE_VERIFY
+
+#ifdef LIBKRIGING_CACHE_ANALYSE
+#include <chrono>
+#define ANALYSE(expr) expr;
+template <typename Timer>
+inline uint64_t diffAndUpdateTimer(Timer& t_init) {
+  const auto t_end = std::chrono::high_resolution_clock::now();
+  const auto diff = (t_end - t_init).count();  // ns count
+  t_init = std::move(t_end);
+  return diff;
+}
+#else
+#define ANALYSE(expr)
+#endif
+
+/* ----------------------------------------------------------------------------------------------------------------- */
+
+class CacheFunctionCommon {
+ public:
+  struct CacheStat {
+    uint32_t min_hit;
+    uint32_t max_hit;
+    uint32_t total_hit;
+    size_t cache_size;
+    uint64_t hash_time;
+    uint64_t lookup_time;
+    uint64_t eval_time;
+  };
+
+ protected:
+  using HashKey = std::size_t;
+
+ public:
+  LIBKRIGING_EXPORT ~CacheFunctionCommon();
+  LIBKRIGING_EXPORT auto stat() -> CacheStat;
+
+ protected:
+  mutable std::unordered_map<HashKey, uint32_t> m_cache_hit;
+#ifdef LIBKRIGING_CACHE_ANALYSE
+  mutable uint64_t m_hash_timer = 0;
+  mutable uint64_t m_lookup_timer = 0;
+  mutable uint64_t m_eval_timer = 0;
+#endif
+};
+
+/* ----------------------------------------------------------------------------------------------------------------- */
+
+std::ostream& operator<<(std::ostream& o, const CacheFunctionCommon::CacheStat&);
+
+/* ----------------------------------------------------------------------------------------------------------------- */
+
+template <typename Callable, typename Signature, typename... Contexts>
+class CacheFunction {};
+
+template <typename Callable,    // true function type for performance; could be a lambda
+          typename R,           // return type
+          typename... Args,     // input parameters
+          typename... Contexts  // non-const external context to manage
+          >
+class CacheFunction<Callable, std::function<R(Args...)>, Contexts...> : public CacheFunctionCommon {
+ public:
+  using type = R(Args...);
+
+ public:
+  LIBKRIGING_EXPORT explicit CacheFunction(const Callable& callable, const Contexts&... contexts)
+      : m_callable(callable), m_context(contexts...) {}
+
+  LIBKRIGING_EXPORT auto operator()(Args... args) const -> R {
+    ANALYSE(auto t = std::chrono::high_resolution_clock::now());
+    const auto arg_key = hash_args(args...);
+    ANALYSE(m_hash_timer += diffAndUpdateTimer(t));
+    auto [finder, is_new] = m_cache.emplace(arg_key, R{});
+    ANALYSE(m_lookup_timer += diffAndUpdateTimer(t));
+    ++m_cache_hit[arg_key];
+    ANALYSE(diffAndUpdateTimer(t));
+    if (is_new) {
+      finder->second = m_callable(std::forward<Args>(args)...);
+      ANALYSE(m_eval_timer += diffAndUpdateTimer(t));
+    } else {
+#ifdef LIBKRIGING_CACHE_VERIFY
+      auto expected_result = m_callable(std::forward<Args>(args)...);
+      assert(expected_result == finder->second);
+#endif
+    }
+    return finder->second;
+  }
+
+ public:
+  auto hash_args(Args... args) const -> HashKey {
+    const auto args_hash = details::tupleHash(std::forward_as_tuple(args...));
+    if constexpr (std::tuple_size<decltype(m_context)>{} == 0) {
+      return args_hash;
+    } else {
+      const auto context_hash = details::tupleHash(m_context);
+      return details::composeHash(context_hash, args_hash);
+    }
+  }
+
+  auto inspect(Args... args) -> uint32_t {
+    const auto arg_key = hash_args(args...);
+    const auto finder = m_cache_hit.find(arg_key);
+    return (finder == m_cache_hit.end()) ? 0 : finder->second;
+  }
+
+ private:
+  Callable m_callable;
+  std::tuple<const Contexts&...> m_context;
+  //  mutable LinearHashStorage<HashKey, R> m_cache; // this struct could be simpler to optimize as circular buffer
+  mutable std::unordered_map<HashKey, R> m_cache;
+};
+
+/* ----------------------------------------------------------------------------------------------------------------- */
+
+template <typename F, typename... Contexts>
+CacheFunction(const F& f, const Contexts&...) -> CacheFunction<F, decltype(std::function{f}), Contexts...>;
+
+#endif  // LIBKRIGING_SRC_LIB_CACHE_HPP
diff --git a/src/lib/include/libKriging/utils/LinearHashStorage.hpp b/src/lib/include/libKriging/utils/LinearHashStorage.hpp
new file mode 100644
index 00000000..a7fbefc1
--- /dev/null
+++ b/src/lib/include/libKriging/utils/LinearHashStorage.hpp
@@ -0,0 +1,38 @@
+#ifndef LIBKRIGING_SRC_LIB_INCLUDE_LIBKRIGING_UTILS_LINEARHASHSTORAGE_HPP
+#define LIBKRIGING_SRC_LIB_INCLUDE_LIBKRIGING_UTILS_LINEARHASHSTORAGE_HPP
+
+#include <utility>
+#include <vector>
+
+template <typename Key, typename T>
+class LinearHashStorage {
+ public:
+  using value_type = std::pair<const Key, T>;
+  struct iterator {
+    T& second;
+    iterator* operator->() { return this; }
+  };
+
+ public:
+  LinearHashStorage(size_t reservation = 20) {
+    m_keys.reserve(reservation);
+    m_data.reserve(reservation);
+  }
+
+  std::pair<iterator, bool> emplace(const Key& key, T&& data) {
+    for (size_t pos = 0; pos < m_keys.size(); ++pos) {
+      if (key == m_keys[pos]) {
+        return {iterator{m_data[pos]}, false};
+      }
+    }
+    m_keys.emplace_back(key);
+    m_data.emplace_back(std::forward<T>(data));
+    return {iterator{m_data.back()}, true};
+  }
+
+ private:
+  std::vector<Key> m_keys;
+  std::vector<T> m_data;
+};
+
+#endif  // LIBKRIGING_SRC_LIB_INCLUDE_LIBKRIGING_UTILS_LINEARHASHSTORAGE_HPP
diff --git a/src/lib/include/libKriging/utils/cache_details.hpp b/src/lib/include/libKriging/utils/cache_details.hpp
new file mode 100644
index 00000000..f470eb2d
--- /dev/null
+++ b/src/lib/include/libKriging/utils/cache_details.hpp
@@ -0,0 +1,44 @@
+#ifndef LIBKRIGING_SRC_LIB_INCLUDE_LIBKRIGING_DETAILS_HPP
+#define LIBKRIGING_SRC_LIB_INCLUDE_LIBKRIGING_DETAILS_HPP
+
+#include <functional>
+#include <tuple>
+#include <utility>
+
+namespace details {
+
+inline std::size_t composeHash(std::size_t acc, std::size_t new_hash) {
+  return acc ^ (new_hash + 0x9e3779b9 + (acc << 6) + (acc >> 2));  // should not be symmetric => this hack
+  //                                                             // same magic number as in boost::hash
+  // https://stackoverflow.com/questions/4948780/magic-number-in-boosthash-combine/4948967#4948967
+}
+
+template <typename T>
+static std::size_t hashValue(T&& t) {
+  return std::hash<std::decay_t<T>>{}(std::forward<T>(t));
+}
+
+template <typename Tuple, std::size_t... ids>
+static std::size_t _tupleHash(Tuple&& tuple, const std::index_sequence<ids...>&&) {
+  if constexpr (sizeof...(ids) == 0) {
+    return 1;
+  } else if constexpr (sizeof...(ids) == 1) {
+    return hashValue(std::get<0>(tuple));
+  } else {
+    std::size_t result = 0;
+#pragma unroll
+    for (auto const& hash : {hashValue(std::get<ids>(tuple))...}) {
+      result = composeHash(result, hash);
+    }
+    return result;
+  }
+}
+
+template <typename Tuple>
+static std::size_t tupleHash(Tuple&& tuple) {
+  return _tupleHash(std::forward<Tuple>(tuple), std::make_index_sequence<std::tuple_size<std::decay_t<Tuple>>{}>{});
+}
+
+}  // namespace details
+
+#endif  // LIBKRIGING_SRC_LIB_INCLUDE_LIBKRIGING_DETAILS_HPP
diff --git a/src/lib/include/libKriging/utils/custom_hash_function.hpp b/src/lib/include/libKriging/utils/custom_hash_function.hpp
new file mode 100644
index 00000000..82f3d8a1
--- /dev/null
+++ b/src/lib/include/libKriging/utils/custom_hash_function.hpp
@@ -0,0 +1,76 @@
+#ifndef LIBKRIGING_SRC_LIB_INCLUDE_LIBKRIGING_HASH_FUNCTION_HPP
+#define LIBKRIGING_SRC_LIB_INCLUDE_LIBKRIGING_HASH_FUNCTION_HPP
+
+#include <armadillo>
+#include <functional>
+#include "libKriging/Kriging.hpp"
+#include "libKriging/utils/cache_details.hpp"
+
+template <>
+struct std::hash<arma::vec> {
+  std::size_t operator()(const arma::vec& s) const noexcept {
+    size_t result = 0;
+    std::hash<arma::vec::value_type> hasher{};
+    for (auto&& e : s) {
+      size_t hash = hasher(e);
+      result = details::composeHash(result, hash);
+    }
+    return result;
+  }
+};
+
+template <>
+struct std::hash<arma::mat> {
+  std::size_t operator()(const arma::mat& s) const noexcept {
+    size_t result = 0;
+    std::hash<arma::mat::value_type> hasher{};
+    for (auto&& e : s) {
+      size_t hash = hasher(e);
+      result = details::composeHash(result, hash);
+    }
+    return result;
+  }
+};
+
+template <>
+struct std::hash<arma::vec*> {
+  std::size_t operator()(const arma::vec* s) const noexcept {
+    if (s != nullptr)
+      return std::hash<arma::vec>{}(*s);
+    else
+      return 0;
+  }
+};
+
+template <>
+struct std::hash<arma::mat*> {
+  std::size_t operator()(const arma::mat* s) const noexcept {
+    if (s != nullptr)
+      return std::hash<arma::mat>{}(*s);
+    else
+      return 0;
+  }
+};
+
+template <>
+struct std::hash<Kriging::OKModel*> {
+  std::size_t operator()(const Kriging::OKModel* s) const noexcept {
+    if (s != nullptr) {
+      struct OKModel {
+        arma::mat T;
+        arma::mat M;
+        arma::colvec z;
+        arma::colvec beta;
+        bool estim_beta;
+        double sigma2;
+        bool estim_sigma2;
+      };
+
+      //      return std::hash<arma::mat>{}(*s);
+    } else {
+      return 0;
+    }
+  }
+};
+
+#endif  // LIBKRIGING_SRC_LIB_INCLUDE_LIBKRIGING_HASH_FUNCTION_HPP
diff --git a/src/lib/include/libKriging/utils/utils.hpp b/src/lib/include/libKriging/utils/utils.hpp
index 1a6e4689..b9da5a59 100644
--- a/src/lib/include/libKriging/utils/utils.hpp
+++ b/src/lib/include/libKriging/utils/utils.hpp
@@ -1,6 +1,10 @@
 #ifndef LIBKRIGING_UTILS_HPP
 #define LIBKRIGING_UTILS_HPP
 
+#include <ostringstream>
+#include <string>
+#include <utility>
+
 template <typename... Args>
 std::string asString(Args&&... args) {
   std::ostringstream oss;


Test patch:
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index ce034852..51e5c7d6 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -31,6 +31,11 @@ add_executable(KrigingTest KrigingTest.cpp)
 target_link_libraries(KrigingTest LINK_PUBLIC Kriging Catch2)
 add_dependencies(all_test_binaries KrigingTest)
 
+add_executable(CacheTest CacheFunctionTest.cpp)
+target_link_libraries(CacheTest LINK_PUBLIC Kriging Catch2)
+add_dependencies(all_test_binaries CacheTest)
+
+
 # Can be enabled if necessary
 #if (CXX_CLANG_TIDY)
 #    foreach (bin ${all_test_binaries})
@@ -57,6 +62,7 @@ endif ()
 # Test using catch2 (parse catch tests)
 catch_discover_tests(catch2_unit_test)
 catch_discover_tests(regression_unit_test)
+catch_discover_tests(CacheTest)
 
 set(AdditionalCatchParameters --benchmark-samples 10)
 catch_discover_tests(KrigingTest)
diff --git a/tests/CacheFunctionTest.cpp b/tests/CacheFunctionTest.cpp
new file mode 100644
index 00000000..e8d5f0a4
--- /dev/null
+++ b/tests/CacheFunctionTest.cpp
@@ -0,0 +1,190 @@
+#define CATCH_CONFIG_MAIN
+#include <catch2/catch.hpp>
+#include <tuple>
+#include "libKriging/CacheFunction.hpp"
+
+TEST_CASE("Cache of 0-arg function", "[core]") {
+  auto f = []() -> double { return 1; };
+
+  CacheFunction f_cached(f);
+
+  REQUIRE(f_cached.inspect() == 0);
+  REQUIRE(f_cached() == f());
+  REQUIRE(f_cached.inspect() == 1);
+  REQUIRE(f_cached() == f());
+  REQUIRE(f_cached.inspect() == 2);
+
+  auto stat = f_cached.stat();
+  REQUIRE(stat.min_hit == 2);
+  REQUIRE(stat.total_hit == 2);
+  REQUIRE(stat.cache_size == 1);
+}
+
+TEST_CASE("Cache of 1-arg function", "[core]") {
+  auto f = [](double x) -> double { return x; };
+
+  SECTION("Same arg should increase hint count") {
+    CacheFunction f_cached(f);
+
+    double x = 1;
+    REQUIRE(f_cached.inspect(x) == 0);
+    REQUIRE(f_cached(x) == f(x));
+    REQUIRE(f_cached.inspect(x) == 1);
+    REQUIRE(f_cached(x) == f(x));
+    REQUIRE(f_cached.inspect(x) == 2);
+
+    auto stat = f_cached.stat();
+    REQUIRE(stat.min_hit == 2);
+    REQUIRE(stat.total_hit == 2);
+    REQUIRE(stat.cache_size == 1);
+  }
+
+  SECTION("Different args should not increase hint count") {
+    CacheFunction f_cached(f);
+
+    double x = 1;
+    REQUIRE(f_cached.inspect(x) == 0);
+    REQUIRE(f_cached(x) == f(x));
+    REQUIRE(f_cached.inspect(x) == 1);
+    x = 2;
+    REQUIRE(f_cached.inspect(x) == 0);
+    REQUIRE(f_cached(x) == f(x));
+    REQUIRE(f_cached.inspect(x) == 1);
+
+    auto stat = f_cached.stat();
+    REQUIRE(stat.min_hit == 1);
+    REQUIRE(stat.total_hit == 2);
+    REQUIRE(stat.cache_size == 2);
+  }
+}
+
+TEST_CASE("Cache of 2-args function", "[core]") {
+  auto f = [](double x, double y) -> double { return x + y; };
+
+  SECTION("Same arg should increase hint count") {
+    CacheFunction f_cached(f);
+
+    double x = 1, y = 2;
+    REQUIRE(f_cached.inspect(x, y) == 0);
+    REQUIRE(f_cached(x, y) == f(x, y));
+    REQUIRE(f_cached.inspect(x, y) == 1);
+    REQUIRE(f_cached(x, y) == f(x, y));
+    REQUIRE(f_cached.inspect(x, y) == 2);
+
+    auto stat = f_cached.stat();
+    REQUIRE(stat.min_hit == 2);
+    REQUIRE(stat.total_hit == 2);
+    REQUIRE(stat.cache_size == 1);
+  }
+
+  SECTION("Different args should not increase hint count") {
+    CacheFunction f_cached(f);
+
+    double x = 1, y = 2;
+    REQUIRE(f_cached.inspect(x, y) == 0);
+    REQUIRE(f_cached(x, y) == f(x, y));
+    REQUIRE(f_cached.inspect(x, y) == 1);
+    x = 2;  // different order should be a different value
+    y = 1;
+    REQUIRE(f_cached.inspect(x, y) == 0);
+    REQUIRE(f_cached(x, y) == f(x, y));
+    REQUIRE(f_cached.inspect(x, y) == 1);
+
+    auto stat = f_cached.stat();
+    REQUIRE(stat.min_hit == 1);
+    REQUIRE(stat.total_hit == 2);
+    REQUIRE(stat.cache_size == 2);
+  }
+}
+
+TEST_CASE("Cache of 0-arg function with global context", "[core]") {
+  double context = 1;
+  auto f = [&context]() -> double { return context; };
+
+  SECTION("Global context change should increase hit count") {
+    CacheFunction f_cached(f, context);
+
+    REQUIRE(f_cached.inspect() == 0);
+    REQUIRE(f_cached() == f());
+    REQUIRE(f_cached.inspect() == 1);
+    context = 2;
+    REQUIRE(f_cached.inspect() == 0);
+    REQUIRE(f_cached() == f());
+    REQUIRE(f_cached.inspect() == 1);
+
+    auto stat = f_cached.stat();
+    REQUIRE(stat.max_hit == 1);
+    REQUIRE(stat.total_hit == 2);
+    REQUIRE(stat.cache_size == 2);
+  }
+}
+
+TEST_CASE("Cache of 1-arg function with global context", "[core]") {
+  double context = 1;
+  auto f = [&context](double x) -> double { return x * context; };
+
+  SECTION("Global context change should increase hint count") {
+    CacheFunction f_cached(f, context);
+
+    double x = 1;
+    REQUIRE(f_cached.inspect(x) == 0);
+    REQUIRE(f_cached(x) == f(x));
+    REQUIRE(f_cached.inspect(x) == 1);
+    context = 2;
+    REQUIRE(f_cached.inspect(x) == 0);
+    REQUIRE(f_cached(x) == f(x));
+    REQUIRE(f_cached.inspect(x) == 1);
+
+    auto stat = f_cached.stat();
+    REQUIRE(stat.max_hit == 1);
+    REQUIRE(stat.total_hit == 2);
+    REQUIRE(stat.cache_size == 2);
+  }
+}
+
+TEST_CASE("Cache of 2-args function with global context", "[core]") {
+  double context = 1;
+  auto f = [&context](double x, double y) -> double { return (x + y) * context; };
+
+  SECTION("Global context change should increase hint count") {
+    CacheFunction f_cached(f, context);
+
+    double x = 1, y = 2;
+    REQUIRE(f_cached.inspect(x, y) == 0);
+    REQUIRE(f_cached(x, y) == f(x, y));
+    REQUIRE(f_cached.inspect(x, y) == 1);
+    context = 2;
+    REQUIRE(f_cached.inspect(x, y) == 0);
+    REQUIRE(f_cached(x, y) == f(x, y));
+    REQUIRE(f_cached.inspect(x, y) == 1);
+
+    auto stat = f_cached.stat();
+    REQUIRE(stat.max_hit == 1);
+    REQUIRE(stat.total_hit == 2);
+    REQUIRE(stat.cache_size == 2);
+  }
+}
+
+TEST_CASE("Cache function behaves like std::function", "[core]") {
+  double context = 1;
+  auto f = [&context](double x, double y) -> double { return (x + y) * context; };
+  CacheFunction f_cached(f, context);
+
+  std::function<decltype(f_cached)::type> f2(f_cached);
+  auto f_cached2 = f2.target<decltype(f_cached)>();
+  REQUIRE(f_cached2 != nullptr);
+
+  double x = 1, y = 2;
+  REQUIRE(f_cached2->inspect(x, y) == 0);
+  REQUIRE(f2(x, y) == f(x, y));
+  REQUIRE(f_cached2->inspect(x, y) == 1);
+  context = 2;
+  REQUIRE(f_cached2->inspect(x, y) == 0);
+  REQUIRE(f2(x, y) == f(x, y));
+  REQUIRE(f_cached2->inspect(x, y) == 1);
+
+  auto stat = f_cached2->stat();
+  REQUIRE(stat.max_hit == 1);
+  REQUIRE(stat.total_hit == 2);
+  REQUIRE(stat.cache_size == 2);
+}
\ No newline at end of file