Repository: libKriging/libKriging (language: python)

Issue:
Store logLikelihood best value (ie. after fit) in Kriging object
may be using a void theta arg when calling k.logLikelihood(theta, withgrad, withHess) from wrappers?
should also be available for leaveOneOut and logMargPost


Discussion / hints:

Fix patch:
diff --git a/README.md b/README.md
index 688393c1..b9ccdd23 100644
--- a/README.md
+++ b/README.md
@@ -57,7 +57,7 @@ import pylibkriging as lk
 k_py = lk.Kriging(y, X, "gauss")
 print(k_py.summary())
 # you can also check logLikelhood using:
-#def ll(t): return k_py.logLikelihood(t,False,False)[0]
+#def ll(t): return k_py.logLikelihoodFun(t,False,False)[0]
 #t = np.arange(0,1,1/99); pyplot.figure(1); pyplot.plot(t, [ll(ti) for ti in t]); pyplot.show()
 
 x = np.arange(0, 1, 1 / 99)
@@ -122,7 +122,7 @@ library(rlibkriging)
 k_R <- Kriging(y, X, "gauss")
 print(k_R)
 # you can also check logLikelhood using:
-# ll = function(t) logLikelihood(k_R,t)$logLikelihood; plot(ll)
+# ll = function(t) logLikelihoodFun(k_R,t)$logLikelihoodFun; plot(ll)
 x <- as.matrix(seq(0, 1, , 100))
 p <- predict(k_R, x, TRUE, FALSE)
 
@@ -170,7 +170,7 @@ y = f(X);
 k_m = Kriging(y, X, "gauss");
 disp(k_m.summary());
 % you can also check logLikelhood using:
-% function llt = ll (tt) global k_m; llt=k_m.logLikelihood(tt); endfunction; t=0:(1/99):1; plot(t,arrayfun(@ll,t))
+% function llt = ll (tt) global k_m; llt=k_m.logLikelihoodFun(tt); endfunction; t=0:(1/99):1; plot(t,arrayfun(@ll,t))
 x = reshape(0:(1/99):1,100,1);
 [p_mean, p_stdev] = k_m.predict(x, true, false);
 
diff --git a/bindings/Octave/Kriging.m b/bindings/Octave/Kriging.m
index c5a35c5e..17e55034 100644
--- a/bindings/Octave/Kriging.m
+++ b/bindings/Octave/Kriging.m
@@ -39,14 +39,26 @@ function fit(obj, varargin)
             [varargout{1:nargout}] = mLibKriging("Kriging::summary", obj.ref, varargin{:});
         end
 
+        function varargout = leaveOneOutFun(obj, varargin)
+            [varargout{1:nargout}] = mLibKriging("Kriging::leaveOneOutFun", obj.ref, varargin{:});
+        end
+
         function varargout = leaveOneOut(obj, varargin)
             [varargout{1:nargout}] = mLibKriging("Kriging::leaveOneOut", obj.ref, varargin{:});
         end
 
+        function varargout = logLikelihoodFun(obj, varargin)
+            [varargout{1:nargout}] = mLibKriging("Kriging::logLikelihoodFun", obj.ref, varargin{:});
+        end
+
         function varargout = logLikelihood(obj, varargin)
             [varargout{1:nargout}] = mLibKriging("Kriging::logLikelihood", obj.ref, varargin{:});
         end
 
+        function varargout = logMargPostFun(obj, varargin)
+            [varargout{1:nargout}] = mLibKriging("Kriging::logMargPostFun", obj.ref, varargin{:});
+        end
+
         function varargout = logMargPost(obj, varargin)
             [varargout{1:nargout}] = mLibKriging("Kriging::logMargPost", obj.ref, varargin{:});
         end
diff --git a/bindings/Octave/Kriging_binding.cpp b/bindings/Octave/Kriging_binding.cpp
index 2777e0ac..d7167716 100644
--- a/bindings/Octave/Kriging_binding.cpp
+++ b/bindings/Octave/Kriging_binding.cpp
@@ -139,7 +139,7 @@ void summary(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs) {
   output.set<0>(km->summary(), "Model description");
 }
 
-void leaveOneOut(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs) {
+void leaveOneOutFun(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs) {
   MxMapper input{"Input",
                  nrhs,
                  const_cast<mxArray**>(prhs),  // NOLINT(cppcoreguidelines-pro-type-const-cast)
@@ -147,12 +147,22 @@ void leaveOneOut(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs) {
   MxMapper output{"Output", nlhs, plhs, RequiresArg::Range{1, 2}};
   auto* km = input.getObject<0, Kriging>("Kriging reference");
   const bool want_grad = flag_output_compliance<2>(input, "want_grad", output, 1);
-  auto [loo, loograd] = km->leaveOneOutEval(input.get<1, arma::vec>("theta"), want_grad);
+  auto [loo, loograd] = km->leaveOneOutFun(input.get<1, arma::vec>("theta"), want_grad);
   output.set<0>(loo, "loo");                  // FIXME better name
   output.setOptional<1>(loograd, "loograd");  // FIXME better name
 }
 
-void logLikelihood(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs) {
+void leaveOneOut(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs) {
+  MxMapper input{"Input",
+                 nrhs,
+                 const_cast<mxArray**>(prhs),  // NOLINT(cppcoreguidelines-pro-type-const-cast)
+                 RequiresArg::Exactly{1}};
+  MxMapper output{"Output", nlhs, plhs, RequiresArg::Exactly{1}};
+  auto* km = input.getObject<0, Kriging>("Kriging reference");
+  output.set<0>(km->leaveOneOut(), "Model leaveOneOut");
+}
+
+void logLikelihoodFun(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs) {
   MxMapper input{"Input",
                  nrhs,
                  const_cast<mxArray**>(prhs),  // NOLINT(cppcoreguidelines-pro-type-const-cast)
@@ -161,13 +171,23 @@ void logLikelihood(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs) {
   auto* km = input.getObject<0, Kriging>("Kriging reference");
   const bool want_grad = flag_output_compliance<2>(input, "want_grad", output, 1);
   const bool want_hess = flag_output_compliance<3>(input, "want_hess", output, 2);
-  auto [ll, llgrad, llhess] = km->logLikelihoodEval(input.get<1, arma::vec>("theta"), want_grad, want_hess);
+  auto [ll, llgrad, llhess] = km->logLikelihoodFun(input.get<1, arma::vec>("theta"), want_grad, want_hess);
   output.set<0>(ll, "ll");                  // FIXME better name
   output.setOptional<1>(llgrad, "llgrad");  // FIXME better name
   output.setOptional<2>(llhess, "llhess");  // FIXME better name
 }
 
-void logMargPost(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs) {
+void logLikelihood(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs) {
+  MxMapper input{"Input",
+                 nrhs,
+                 const_cast<mxArray**>(prhs),  // NOLINT(cppcoreguidelines-pro-type-const-cast)
+                 RequiresArg::Exactly{1}};
+  MxMapper output{"Output", nlhs, plhs, RequiresArg::Exactly{1}};
+  auto* km = input.getObject<0, Kriging>("Kriging reference");
+  output.set<0>(km->logLikelihood(), "Model logLikelihood");
+}
+
+void logMargPostFun(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs) {
   MxMapper input{"Input",
                  nrhs,
                  const_cast<mxArray**>(prhs),  // NOLINT(cppcoreguidelines-pro-type-const-cast)
@@ -175,11 +195,21 @@ void logMargPost(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs) {
   MxMapper output{"Output", nlhs, plhs, RequiresArg::Range{1, 2}};
   auto* km = input.getObject<0, Kriging>("Kriging reference");
   const bool want_grad = flag_output_compliance<2>(input, "want_grad", output, 1);
-  auto [lmp, lmpgrad] = km->logMargPostEval(input.get<1, arma::vec>("theta"), want_grad);
+  auto [lmp, lmpgrad] = km->logMargPostFun(input.get<1, arma::vec>("theta"), want_grad);
   output.set<0>(lmp, "lmp");                  // FIXME better name
   output.setOptional<1>(lmpgrad, "lmpgrad");  // FIXME better name
 }
 
+void logMargPost(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs) {
+  MxMapper input{"Input",
+                 nrhs,
+                 const_cast<mxArray**>(prhs),  // NOLINT(cppcoreguidelines-pro-type-const-cast)
+                 RequiresArg::Exactly{1}};
+  MxMapper output{"Output", nlhs, plhs, RequiresArg::Exactly{1}};
+  auto* km = input.getObject<0, Kriging>("Kriging reference");
+  output.set<0>(km->logMargPost(), "Model logMargPost");
+}
+
 void kernel(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs) {
   MxMapper input{"Input",
                  nrhs,
diff --git a/bindings/Octave/Kriging_binding.hpp b/bindings/Octave/Kriging_binding.hpp
index 184f2d90..960623c3 100644
--- a/bindings/Octave/Kriging_binding.hpp
+++ b/bindings/Octave/Kriging_binding.hpp
@@ -11,6 +11,9 @@ void predict(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs);
 void simulate(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs);
 void update(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs);
 void summary(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs);
+void leaveOneOutFun(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs);
+void logLikelihoodFun(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs);
+void logMargPostFun(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs);
 void leaveOneOut(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs);
 void logLikelihood(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs);
 void logMargPost(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs);
diff --git a/bindings/Octave/NuggetKriging.m b/bindings/Octave/NuggetKriging.m
index 6b680f2b..7c3d8b36 100644
--- a/bindings/Octave/NuggetKriging.m
+++ b/bindings/Octave/NuggetKriging.m
@@ -39,9 +39,17 @@ function fit(obj, varargin)
             [varargout{1:nargout}] = mLibKriging("NuggetKriging::summary", obj.ref, varargin{:});
         end
 
+        function varargout = logLikelihoodFun(obj, varargin)
+            [varargout{1:nargout}] = mLibKriging("NuggetKriging::logLikelihoodFun", obj.ref, varargin{:});
+        end
+
         function varargout = logLikelihood(obj, varargin)
             [varargout{1:nargout}] = mLibKriging("NuggetKriging::logLikelihood", obj.ref, varargin{:});
         end
+        
+        function varargout = logMargPostFun(obj, varargin)
+            [varargout{1:nargout}] = mLibKriging("NuggetKriging::logMargPostFun", obj.ref, varargin{:});
+        end
 
         function varargout = logMargPost(obj, varargin)
             [varargout{1:nargout}] = mLibKriging("NuggetKriging::logMargPost", obj.ref, varargin{:});
diff --git a/bindings/Octave/NuggetKriging_binding.cpp b/bindings/Octave/NuggetKriging_binding.cpp
index 077c56c9..714ed95e 100644
--- a/bindings/Octave/NuggetKriging_binding.cpp
+++ b/bindings/Octave/NuggetKriging_binding.cpp
@@ -141,7 +141,7 @@ void summary(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs) {
   output.set<0>(km->summary(), "Model description");
 }
 
-void logLikelihood(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs) {
+void logLikelihoodFun(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs) {
   MxMapper input{"Input",
                  nrhs,
                  const_cast<mxArray**>(prhs),  // NOLINT(cppcoreguidelines-pro-type-const-cast)
@@ -149,12 +149,22 @@ void logLikelihood(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs) {
   MxMapper output{"Output", nlhs, plhs, RequiresArg::Range{1, 2}};
   auto* km = input.getObject<0, NuggetKriging>("NuggetKriging reference");
   const bool want_grad = flag_output_compliance<2>(input, "want_grad", output, 1);
-  auto [ll, llgrad] = km->logLikelihoodEval(input.get<1, arma::vec>("theta"), want_grad);
+  auto [ll, llgrad] = km->logLikelihoodFun(input.get<1, arma::vec>("theta"), want_grad);
   output.set<0>(ll, "ll");                  // FIXME better name
   output.setOptional<1>(llgrad, "llgrad");  // FIXME better name
 }
 
-void logMargPost(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs) {
+void logLikelihood(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs) {
+  MxMapper input{"Input",
+                 nrhs,
+                 const_cast<mxArray**>(prhs),  // NOLINT(cppcoreguidelines-pro-type-const-cast)
+                 RequiresArg::Exactly{1}};
+  MxMapper output{"Output", nlhs, plhs, RequiresArg::Exactly{1}};
+  auto* km = input.getObject<0, NuggetKriging>("NuggetKriging reference");
+  output.set<0>(km->logLikelihood(), "Model logLikelihood");
+}
+
+void logMargPostFun(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs) {
   MxMapper input{"Input",
                  nrhs,
                  const_cast<mxArray**>(prhs),  // NOLINT(cppcoreguidelines-pro-type-const-cast)
@@ -162,11 +172,21 @@ void logMargPost(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs) {
   MxMapper output{"Output", nlhs, plhs, RequiresArg::Range{1, 2}};
   auto* km = input.getObject<0, NuggetKriging>("NuggetKriging reference");
   const bool want_grad = flag_output_compliance<2>(input, "want_grad", output, 1);
-  auto [lmp, lmpgrad] = km->logMargPostEval(input.get<1, arma::vec>("theta"), want_grad);
+  auto [lmp, lmpgrad] = km->logMargPostFun(input.get<1, arma::vec>("theta"), want_grad);
   output.set<0>(lmp, "lmp");                  // FIXME better name
   output.setOptional<1>(lmpgrad, "lmpgrad");  // FIXME better name
 }
 
+void logMargPost(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs) {
+  MxMapper input{"Input",
+                 nrhs,
+                 const_cast<mxArray**>(prhs),  // NOLINT(cppcoreguidelines-pro-type-const-cast)
+                 RequiresArg::Exactly{1}};
+  MxMapper output{"Output", nlhs, plhs, RequiresArg::Exactly{1}};
+  auto* km = input.getObject<0, NuggetKriging>("NuggetKriging reference");
+  output.set<0>(km->logMargPost(), "Model logMargPost");
+}
+
 void kernel(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs) {
   MxMapper input{"Input",
                  nrhs,
diff --git a/bindings/Octave/NuggetKriging_binding.hpp b/bindings/Octave/NuggetKriging_binding.hpp
index 354a92a5..1c593051 100644
--- a/bindings/Octave/NuggetKriging_binding.hpp
+++ b/bindings/Octave/NuggetKriging_binding.hpp
@@ -11,6 +11,8 @@ void predict(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs);
 void simulate(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs);
 void update(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs);
 void summary(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs);
+void logLikelihoodFun(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs);
+void logMargPostFun(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs);
 void logLikelihood(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs);
 void logMargPost(int nlhs, mxArray** plhs, int nrhs, const mxArray** prhs);
 
diff --git a/bindings/Octave/mLibKriging.cpp b/bindings/Octave/mLibKriging.cpp
index 05fae3dc..8456bcad 100644
--- a/bindings/Octave/mLibKriging.cpp
+++ b/bindings/Octave/mLibKriging.cpp
@@ -75,6 +75,12 @@ void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) try
       return KrigingBinding::update(nlhs, plhs, nrhs - 1, prhs + 1);
     case "Kriging::summary"_hash:
       return KrigingBinding::summary(nlhs, plhs, nrhs - 1, prhs + 1);
+    case "Kriging::leaveOneOutFun"_hash:
+      return KrigingBinding::leaveOneOutFun(nlhs, plhs, nrhs - 1, prhs + 1);
+    case "Kriging::logLikelihoodFun"_hash:
+      return KrigingBinding::logLikelihoodFun(nlhs, plhs, nrhs - 1, prhs + 1);
+    case "Kriging::logMargPostFun"_hash:
+      return KrigingBinding::logMargPostFun(nlhs, plhs, nrhs - 1, prhs + 1);
     case "Kriging::leaveOneOut"_hash:
       return KrigingBinding::leaveOneOut(nlhs, plhs, nrhs - 1, prhs + 1);
     case "Kriging::logLikelihood"_hash:
@@ -137,6 +143,10 @@ void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) try
       return NuggetKrigingBinding::update(nlhs, plhs, nrhs - 1, prhs + 1);
     case "NuggetKriging::summary"_hash:
       return NuggetKrigingBinding::summary(nlhs, plhs, nrhs - 1, prhs + 1);
+    case "NuggetKriging::logLikelihoodFun"_hash:
+      return NuggetKrigingBinding::logLikelihoodFun(nlhs, plhs, nrhs - 1, prhs + 1);
+    case "NuggetKriging::logMargPostFun"_hash:
+      return NuggetKrigingBinding::logMargPostFun(nlhs, plhs, nrhs - 1, prhs + 1);
     case "NuggetKriging::logLikelihood"_hash:
       return NuggetKrigingBinding::logLikelihood(nlhs, plhs, nrhs - 1, prhs + 1);
     case "NuggetKriging::logMargPost"_hash:
diff --git a/bindings/Python/CMakeLists.txt b/bindings/Python/CMakeLists.txt
index 0d5cb639..97b97bab 100644
--- a/bindings/Python/CMakeLists.txt
+++ b/bindings/Python/CMakeLists.txt
@@ -19,7 +19,7 @@ set(CARMA_DEV_TARGET false)
 # set(BUILD_TESTS on CACHE BOOL "Enable test in carma") 
 add_subdirectory(${LIBKRIGING_SOURCE_DIR}/dependencies/carma carma)
 
-set(pylibkriging_SRCS src/_pylibkriging/pylibkirging.cpp)
+set(pylibkriging_SRCS src/_pylibkriging/pylibkriging.cpp)
 list(APPEND pylibkriging_SRCS
         src/_pylibkriging/LinearRegression_binding.cpp src/_pylibkriging/LinearRegression_binding.hpp
         src/_pylibkriging/Kriging_binding.cpp src/_pylibkriging/Kriging_binding.hpp
diff --git a/bindings/Python/src/_pylibkriging/Kriging_binding.cpp b/bindings/Python/src/_pylibkriging/Kriging_binding.cpp
index 5d155875..b9ae2cd1 100644
--- a/bindings/Python/src/_pylibkriging/Kriging_binding.cpp
+++ b/bindings/Python/src/_pylibkriging/Kriging_binding.cpp
@@ -62,17 +62,21 @@ std::string PyKriging::summary() const {
   return m_internal->summary();
 }
 
-std::tuple<double, py::array_t<double>> PyKriging::leaveOneOutEval(const py::array_t<double>& theta,
-                                                                   const bool want_grad) {
+std::tuple<double, py::array_t<double>> PyKriging::leaveOneOutFun(const py::array_t<double>& theta,
+                                                                  const bool want_grad) {
   arma::vec vec_theta = carma::arr_to_col<double>(theta);
-  auto [llo, grad] = m_internal->leaveOneOutEval(vec_theta, want_grad);
+  auto [llo, grad] = m_internal->leaveOneOutFun(vec_theta, want_grad);
   return {llo, carma::col_to_arr(grad)};
 }
 
+double PyKriging::leaveOneOut() {
+  return m_internal->leaveOneOut();
+}
+
 std::tuple<double, py::array_t<double>, py::array_t<double>>
-PyKriging::logLikelihoodEval(const py::array_t<double>& theta, const bool want_grad, const bool want_hess) {
+PyKriging::logLikelihoodFun(const py::array_t<double>& theta, const bool want_grad, const bool want_hess) {
   arma::vec vec_theta = carma::arr_to_col<double>(theta);
-  auto [llo, grad, hess] = m_internal->logLikelihoodEval(vec_theta, want_grad, want_hess);
+  auto [llo, grad, hess] = m_internal->logLikelihoodFun(vec_theta, want_grad, want_hess);
   return {
       llo,
       carma::col_to_arr(grad),
@@ -81,9 +85,17 @@ PyKriging::logLikelihoodEval(const py::array_t<double>& theta, const bool want_g
   };
 }
 
-std::tuple<double, py::array_t<double>> PyKriging::logMargPostEval(const py::array_t<double>& theta,
-                                                                   const bool want_grad) {
+double PyKriging::logLikelihood() {
+  return m_internal->logLikelihood();
+}
+
+std::tuple<double, py::array_t<double>> PyKriging::logMargPostFun(const py::array_t<double>& theta,
+                                                                  const bool want_grad) {
   arma::vec vec_theta = carma::arr_to_col<double>(theta);
-  auto [lmp, grad] = m_internal->logMargPostEval(vec_theta, want_grad);
+  auto [lmp, grad] = m_internal->logMargPostFun(vec_theta, want_grad);
   return {lmp, carma::col_to_arr(grad)};
 }
+
+double PyKriging::logMargPost() {
+  return m_internal->logMargPost();
+}
diff --git a/bindings/Python/src/_pylibkriging/Kriging_binding.hpp b/bindings/Python/src/_pylibkriging/Kriging_binding.hpp
index 6fe4bed1..863cf377 100644
--- a/bindings/Python/src/_pylibkriging/Kriging_binding.hpp
+++ b/bindings/Python/src/_pylibkriging/Kriging_binding.hpp
@@ -47,13 +47,17 @@ class PyKriging {
 
   std::string summary() const;
 
-  std::tuple<double, py::array_t<double>> leaveOneOutEval(const py::array_t<double>& theta, const bool want_grad);
+  std::tuple<double, py::array_t<double>> leaveOneOutFun(const py::array_t<double>& theta, const bool want_grad);
 
-  std::tuple<double, py::array_t<double>, py::array_t<double>> logLikelihoodEval(const py::array_t<double>& theta,
-                                                                                 const bool want_grad,
-                                                                                 const bool want_hess);
+  std::tuple<double, py::array_t<double>, py::array_t<double>> logLikelihoodFun(const py::array_t<double>& theta,
+                                                                                const bool want_grad,
+                                                                                const bool want_hess);
 
-  std::tuple<double, py::array_t<double>> logMargPostEval(const py::array_t<double>& theta, const bool want_grad);
+  std::tuple<double, py::array_t<double>> logMargPostFun(const py::array_t<double>& theta, const bool want_grad);
+
+  double logLikelihood();
+  double leaveOneOut();
+  double logMargPost();
 
  private:
   std::unique_ptr<Kriging> m_internal;
diff --git a/bindings/Python/src/_pylibkriging/NuggetKriging_binding.cpp b/bindings/Python/src/_pylibkriging/NuggetKriging_binding.cpp
index a3393553..ab3288b6 100644
--- a/bindings/Python/src/_pylibkriging/NuggetKriging_binding.cpp
+++ b/bindings/Python/src/_pylibkriging/NuggetKriging_binding.cpp
@@ -63,11 +63,11 @@ std::string PyNuggetKriging::summary() const {
   return m_internal->summary();
 }
 
-std::tuple<double, py::array_t<double>, py::array_t<double>> PyNuggetKriging::logLikelihoodEval(
+std::tuple<double, py::array_t<double>, py::array_t<double>> PyNuggetKriging::logLikelihoodFun(
     const py::array_t<double>& theta,
     const bool want_grad) {
   arma::vec vec_theta = carma::arr_to_col<double>(theta);
-  auto [llo, grad] = m_internal->logLikelihoodEval(vec_theta, want_grad);
+  auto [llo, grad] = m_internal->logLikelihoodFun(vec_theta, want_grad);
   return {
       llo,
       carma::col_to_arr(grad),
@@ -76,9 +76,17 @@ std::tuple<double, py::array_t<double>, py::array_t<double>> PyNuggetKriging::lo
   };
 }
 
-std::tuple<double, py::array_t<double>> PyNuggetKriging::logMargPostEval(const py::array_t<double>& theta,
-                                                                         const bool want_grad) {
+double PyNuggetKriging::logLikelihood() {
+  return m_internal->logLikelihood();
+}
+
+std::tuple<double, py::array_t<double>> PyNuggetKriging::logMargPostFun(const py::array_t<double>& theta,
+                                                                        const bool want_grad) {
   arma::vec vec_theta = carma::arr_to_col<double>(theta);
-  auto [lmp, grad] = m_internal->logMargPostEval(vec_theta, want_grad);
+  auto [lmp, grad] = m_internal->logMargPostFun(vec_theta, want_grad);
   return {lmp, carma::col_to_arr(grad)};
 }
+
+double PyNuggetKriging::logMargPost() {
+  return m_internal->logMargPost();
+}
diff --git a/bindings/Python/src/_pylibkriging/NuggetKriging_binding.hpp b/bindings/Python/src/_pylibkriging/NuggetKriging_binding.hpp
index 9288b7d7..82654997 100644
--- a/bindings/Python/src/_pylibkriging/NuggetKriging_binding.hpp
+++ b/bindings/Python/src/_pylibkriging/NuggetKriging_binding.hpp
@@ -47,10 +47,13 @@ class PyNuggetKriging {
 
   std::string summary() const;
 
-  std::tuple<double, py::array_t<double>, py::array_t<double>> logLikelihoodEval(const py::array_t<double>& theta,
-                                                                                 const bool want_grad);
+  std::tuple<double, py::array_t<double>, py::array_t<double>> logLikelihoodFun(const py::array_t<double>& theta,
+                                                                                const bool want_grad);
 
-  std::tuple<double, py::array_t<double>> logMargPostEval(const py::array_t<double>& theta, const bool want_grad);
+  std::tuple<double, py::array_t<double>> logMargPostFun(const py::array_t<double>& theta, const bool want_grad);
+
+  double logLikelihood();
+  double logMargPost();
 
  private:
   std::unique_ptr<NuggetKriging> m_internal;
diff --git a/bindings/Python/src/_pylibkriging/pylibkirging.cpp b/bindings/Python/src/_pylibkriging/pylibkriging.cpp
similarity index 91%
rename from bindings/Python/src/_pylibkriging/pylibkirging.cpp
rename to bindings/Python/src/_pylibkriging/pylibkriging.cpp
index aae789c0..5251c473 100644
--- a/bindings/Python/src/_pylibkriging/pylibkirging.cpp
+++ b/bindings/Python/src/_pylibkriging/pylibkriging.cpp
@@ -108,9 +108,9 @@ PYBIND11_MODULE(_pylibkriging, m) {
       .def("simulate", &PyKriging::simulate)
       .def("update", &PyKriging::update)
       .def("summary", &PyKriging::summary)
-      .def("leaveOneOut", &PyKriging::leaveOneOutEval)
-      .def("logLikelihood", &PyKriging::logLikelihoodEval)
-      .def("logMargPost", &PyKriging::logMargPostEval);
+      .def("leaveOneOutFun", &PyKriging::leaveOneOutFun)
+      .def("logLikelihoodFun", &PyKriging::logLikelihoodFun)
+      .def("logMargPostFun", &PyKriging::logMargPostFun);
 
   // Automated mapper
   py::class_<Kriging>(m, "PyKriging")
@@ -144,9 +144,12 @@ PYBIND11_MODULE(_pylibkriging, m) {
       .def("simulate", &Kriging::simulate)
       .def("update", &Kriging::update)
       .def("summary", &Kriging::summary)
-      .def("leaveOneOut", &Kriging::leaveOneOutEval)
-      .def("logLikelihood", &Kriging::logLikelihoodEval)
-      .def("logMargPost", &Kriging::logMargPostEval)
+      .def("leaveOneOutFun", &Kriging::leaveOneOutFun)
+      .def("logLikelihoodFun", &Kriging::logLikelihoodFun)
+      .def("logMargPostFun", &Kriging::logMargPostFun)
+      .def("leaveOneOut", &Kriging::leaveOneOut)
+      .def("logLikelihood", &Kriging::logLikelihood)
+      .def("logMargPost", &Kriging::logMargPost)
 
       .def("kernel", &Kriging::kernel)
       .def("optim", &Kriging::optim)
@@ -197,8 +200,8 @@ PYBIND11_MODULE(_pylibkriging, m) {
       .def("simulate", &PyNuggetKriging::simulate)
       .def("update", &PyNuggetKriging::update)
       .def("summary", &PyNuggetKriging::summary)
-      .def("logLikelihood", &PyNuggetKriging::logLikelihoodEval)
-      .def("logMargPost", &PyNuggetKriging::logMargPostEval);
+      .def("logLikelihoodFun", &PyNuggetKriging::logLikelihoodFun)
+      .def("logMargPostFun", &PyNuggetKriging::logMargPostFun);
 
   // Automated mapper
   py::class_<NuggetKriging>(m, "PyNuggetKriging")
@@ -232,8 +235,10 @@ PYBIND11_MODULE(_pylibkriging, m) {
       .def("simulate", &NuggetKriging::simulate)
       .def("update", &NuggetKriging::update)
       .def("summary", &NuggetKriging::summary)
-      .def("logLikelihood", &NuggetKriging::logLikelihoodEval)
-      .def("logMargPost", &NuggetKriging::logMargPostEval)
+      .def("logLikelihoodFun", &NuggetKriging::logLikelihoodFun)
+      .def("logMargPostFun", &NuggetKriging::logMargPostFun)
+      .def("logLikelihood", &NuggetKriging::logLikelihood)
+      .def("logMargPost", &NuggetKriging::logMargPost)
 
       .def("kernel", &NuggetKriging::kernel)
       .def("optim", &NuggetKriging::optim)
diff --git a/bindings/R/bench_perf-full.Rmd b/bindings/R/bench_perf-full.Rmd
index 96464569..74be9cd2 100644
--- a/bindings/R/bench_perf-full.Rmd
+++ b/bindings/R/bench_perf-full.Rmd
@@ -157,9 +157,9 @@ N <- 100
 ```
 
 
-### logLikelihood (& deriv.)
+### logLikelihoodFun (& deriv.)
 
-```{r logLikelihood}
+```{r logLikelihoodFun}
 times <- list(R=rep(NA, length(logn)), cpp=rep(NA, length(logn)))
 
 theta <- 0.5
@@ -184,7 +184,7 @@ for (i in 1:length(logn)) {
              theta=matrix(k@covariance@range.val,ncol=d)))
   times$cpp[i] = system.time(
     # try({ll2x <- bench_LogLik(N,r,rep(theta,ncol(X)))})
-    try({for (t in 1:N) ll2x <- rlibkriging::logLikelihood(r,rep(theta,ncol(X)))$logLikelihood[1]}) # Loop should be done inside c++, not from R...
+    try({for (t in 1:N) ll2x <- rlibkriging::logLikelihoodFun(r,rep(theta,ncol(X)))$logLikelihoodFun[1]}) # Loop should be done inside c++, not from R...
   )[1]
   
   if (abs((llx-ll2x)/llx)>1E-3) warning("LL is not identical bw C++/R")
@@ -194,7 +194,7 @@ plot_times(times,N,"logLikelihood")
 ```
 
 
-```{r logLikelihood-grad}
+```{r logLikelihoodFun-grad}
 times <- list(R=rep(NA, length(logn)), cpp=rep(NA, length(logn)))
 
 theta <- 0.5
@@ -220,13 +220,13 @@ for (i in 1:length(logn)) {
              theta=matrix(k@covariance@range.val,ncol=d)))
   times$cpp[i] = system.time(
     # try({gll2x <- bench_LogLikGrad(N,r,rep(theta,ncol(X)))})
-    try({for (t in 1:N) gll2x <- logLikelihood(r,rep(theta,ncol(X)),grad=T)$logLikelihoodGrad}) # Loop should be done inside c++, not from R...
+    try({for (t in 1:N) gll2x <- logLikelihoodFun(r,rep(theta,ncol(X)),grad=T)$logLikelihoodGrad}) # Loop should be done inside c++, not from R...
   )[1]
   
   if (max(abs((gllx-gll2x)/gllx))>1E-3) warning("grad LL is not identical bw C++/R")
 }
 
-plot_times(times,N,"logLikelihood Grad")
+plot_times(times,N,"logLikelihoodFun Grad")
 ```
 
 
@@ -258,7 +258,7 @@ for (i in 1:length(logn)) {
              parameters=list(sigma2=k@covariance@sd2,
              theta=matrix(k@covariance@range.val,ncol=d)))
   times$cpp[i] = system.time(
-    try({for (t in 1:N) loo2x <- leaveOneOut(r,rep(theta,ncol(X)))$leaveOneOut[1]}) # Loop should be done inside c++, not from R...
+    try({for (t in 1:N) loo2x <- leaveOneOutFun(r,rep(theta,ncol(X)))$leaveOneOut[1]}) # Loop should be done inside c++, not from R...
   )[1]
   
   if (abs((loox-loo2x)/loox)>1E-3) warning("LOO is not identical bw C++/R")
@@ -295,7 +295,7 @@ for (i in 1:length(logn)) {
              theta=matrix(k@covariance@range.val,ncol=d)))
   times$cpp[i] = system.time(
     # try({gll2x <- bench_LogLikGrad(N,r,rep(theta,ncol(X)))})
-    try({for (t in 1:N) gll2x <- leaveOneOut(r,rep(theta,ncol(X)),grad=T)$leaveOneOutGrad}) # Loop should be done inside c++, not from R...
+    try({for (t in 1:N) gll2x <- leaveOneOutFun(r,rep(theta,ncol(X)),grad=T)$leaveOneOutGrad}) # Loop should be done inside c++, not from R...
   )[1]
   
   if (max(abs((gllx-gll2x)/gllx))>1E-3) warning("grad LOO is not identical bw C++/R")
@@ -307,9 +307,9 @@ plot_times(times,N,"leaveOneOut Grad")
 
 ### fit
 
-#### logLikelihood / BFGS
+#### logLikelihoodFun / BFGS
 
-```{r logLikelihood-BFGS}
+```{r logLikelihoodFun-BFGS}
 times <- list(R=rep(NA, length(logn)), cpp=rep(NA, length(logn)))
 
 for (i in 1:length(logn)) {
@@ -336,15 +336,15 @@ for (i in 1:length(logn)) {
                                     parameters=list(theta=matrix(k@parinit,ncol=d))))
   )[1]
   
-  ll_cpp <- logLikelihood(r, t(as.list(r)$theta))$logLikelihood[1]
+  ll_cpp <- logLikelihoodFun(r, t(as.list(r)$theta))$logLikelihoodFun[1]
   e <- new.env()
   ll_R <- DiceKriging::logLikFun(k@covariance@range.val, k, e)
   
-  gll_cpp <- logLikelihood(r, t(as.list(r)$theta),grad=T)$logLikelihoodGrad
+  gll_cpp <- logLikelihoodFun(r, t(as.list(r)$theta),grad=T)$logLikelihoodGrad
   gll_R <- DiceKriging::logLikGrad(k@covariance@range.val, k, e)
   
   if (abs(ll_cpp - DiceKriging::logLikFun(param=as.numeric(as.list(r)$theta),model=k))/ll_cpp>.1)
-    warning("LL function is not the same bw DiceKriging/libKriging: ",logLikelihood(r,t(as.list(r)$theta))$logLikelihood[1]," vs. ",DiceKriging::logLikFun(param=as.numeric(as.list(r)$theta),model=k))  
+    warning("LL function is not the same bw DiceKriging/libKriging: ",logLikelihoodFun(r,t(as.list(r)$theta))$logLikelihoodFun[1]," vs. ",DiceKriging::logLikFun(param=as.numeric(as.list(r)$theta),model=k))  
   
   if ((ll_cpp - ll_R)/ll_R < -.01 )
     warning("libKriging LL ",ll_cpp," << DiceKriging LL ",ll_R)
@@ -353,7 +353,7 @@ for (i in 1:length(logn)) {
     warning("DiceKriging LL ",ll_R," << libKriging LL ",ll_cpp)
 }
 
-plot_times(times,N,"fit logLikelihood / BFGS")
+plot_times(times,N,"fit logLikelihoodFun / BFGS")
 ```
 
 
@@ -387,15 +387,15 @@ for (i in 1:length(logn)) {
                                     parameters=list(theta=matrix(k@parinit,ncol=d))))
   )[1]
   
-  ll_cpp <- leaveOneOut(r, t(as.list(r)$theta))$leaveOneOut[1]
+  ll_cpp <- leaveOneOutFun(r, t(as.list(r)$theta))$leaveOneOut[1]
   e <- new.env()
   ll_R <- DiceKriging::leaveOneOutFun(k@covariance@range.val, k, e)
   
-  gll_cpp <- leaveOneOut(r, t(as.list(r)$theta),grad=T)$leaveOneOutGrad
+  gll_cpp <- leaveOneOutFun(r, t(as.list(r)$theta),grad=T)$leaveOneOutGrad
   gll_R <- DiceKriging::leaveOneOutGrad(k@covariance@range.val, k, e)
   
   if (abs(ll_cpp - DiceKriging::leaveOneOutFun(param=as.numeric(as.list(r)$theta),model=k))/ll_cpp>.1)
-    warning("LOO function is not the same bw DiceKriging/libKriging: ",leaveOneOut(r,t(as.list(r)$theta))$leaveOneOut[1]," vs. ",DiceKriging::leaveOneOutFun(param=as.numeric(as.list(r)$theta),model=k))  
+    warning("LOO function is not the same bw DiceKriging/libKriging: ",leaveOneOutFun(r,t(as.list(r)$theta))$leaveOneOut[1]," vs. ",DiceKriging::leaveOneOutFun(param=as.numeric(as.list(r)$theta),model=k))  
   
   if ((ll_cpp - ll_R)/ll_R < -.01 )
     warning("libKriging LOO ",ll_cpp," << DiceKriging LOO ",ll_R)
@@ -408,9 +408,9 @@ plot_times(times,N,"fit leaveOneOut / BFGS")
 ```
 
 
-#### logLikelihood / Newton
+#### logLikelihoodFun / Newton
 
-```{r logLikelihood-Newton}
+```{r logLikelihoodFun-Newton}
 times <- list(R=rep(NA, length(logn)), cpp=rep(NA, length(logn)))
 
 for (i in 1:length(logn)) {
@@ -437,15 +437,15 @@ for (i in 1:length(logn)) {
                                     parameters=list(theta=matrix(k@parinit,ncol=d))))
   )[1]
   
-  ll_cpp <- logLikelihood(r, t(as.list(r)$theta))$logLikelihood[1]
+  ll_cpp <- logLikelihoodFun(r, t(as.list(r)$theta))$logLikelihoodFun[1]
   e <- new.env()
   ll_R <- DiceKriging::logLikFun(k@covariance@range.val, k, e)
   
-  gll_cpp <- logLikelihood(r, t(as.list(r)$theta),grad=T)$logLikelihoodGrad
+  gll_cpp <- logLikelihoodFun(r, t(as.list(r)$theta),grad=T)$logLikelihoodGrad
   gll_R <- DiceKriging::logLikGrad(k@covariance@range.val, k, e)
   
   if (abs(ll_cpp - DiceKriging::logLikFun(param=as.numeric(as.list(r)$theta),model=k))/ll_cpp>.1)
-    warning("LL function is not the same bw DiceKriging/libKriging: ",logLikelihood(r,as.list(r)$theta)," vs. ",DiceKriging::logLikFun(param=as.numeric(as.list(r)$theta),model=k))  
+    warning("LL function is not the same bw DiceKriging/libKriging: ",logLikelihoodFun(r,as.list(r)$theta)," vs. ",DiceKriging::logLikFun(param=as.numeric(as.list(r)$theta),model=k))  
   
   if ((ll_cpp - ll_R)/ll_R < -.01 )
     warning("libKriging LL ",ll_cpp," << DiceKriging LL ",ll_R)
@@ -454,7 +454,7 @@ for (i in 1:length(logn)) {
     warning("DiceKriging LL ",ll_R," << libKriging LL ",ll_cpp)
 }
 
-plot_times(times,N,"fit logLikelihood / Newton")
+plot_times(times,N,"fit logLikelihoodFun / Newton")
 ```
 
 
diff --git a/bindings/R/rlibkriging/DiceKriging.md b/bindings/R/rlibkriging/DiceKriging.md
index ab9e0df5..9334f38c 100644
--- a/bindings/R/rlibkriging/DiceKriging.md
+++ b/bindings/R/rlibkriging/DiceKriging.md
@@ -11,8 +11,8 @@
 | predict             | predict()             |              |                 |
 | simulate            | simulate()            |              |                 |
 | update              | update()              |              |                 |
-| logLik              | logLikelihood()       |              |                 |
-| loo                 | leaveOneOut()         |              |                 |
+| logLik              | logLikelihoodFun()       |              |                 |
+| loo                 | leaveOneOutFun()         |              |                 |
 |                     |                       |              |                 |
 | kernel:gauss        | (..,"gauss"           |              |                 |
 | kernel:matern3_2    |                       |              |                 |
diff --git a/bindings/R/rlibkriging/NAMESPACE b/bindings/R/rlibkriging/NAMESPACE
index e27fabab..eb36ad7e 100644
--- a/bindings/R/rlibkriging/NAMESPACE
+++ b/bindings/R/rlibkriging/NAMESPACE
@@ -5,10 +5,15 @@ S3method(as.km,NuggetKriging)
 S3method(as.list,Kriging)
 S3method(as.list,NuggetKriging)
 S3method(leaveOneOut,Kriging)
+S3method(leaveOneOutFun,Kriging)
 S3method(logLikelihood,Kriging)
 S3method(logLikelihood,NuggetKriging)
+S3method(logLikelihoodFun,Kriging)
+S3method(logLikelihoodFun,NuggetKriging)
 S3method(logMargPost,Kriging)
 S3method(logMargPost,NuggetKriging)
+S3method(logMargPostFun,Kriging)
+S3method(logMargPostFun,NuggetKriging)
 S3method(predict,Kriging)
 S3method(predict,NuggetKriging)
 S3method(print,Kriging)
@@ -29,9 +34,18 @@ export(bench_cholsym)
 export(bench_invsympd)
 export(bench_qr)
 export(bench_solvetri)
+export(get_num_nugget)
 export(leaveOneOut)
+export(leaveOneOutFun)
+export(linalg_get_num_nugget)
+export(linalg_set_num_nugget)
 export(logLikelihood)
+export(logLikelihoodFun)
 export(logMargPost)
+export(logMargPostFun)
+export(random_reset_seed)
+export(reset_seed)
+export(set_num_nugget)
 exportMethods(predict)
 exportMethods(simulate)
 exportMethods(update)
diff --git a/bindings/R/rlibkriging/R/KrigingClass.R b/bindings/R/rlibkriging/R/KrigingClass.R
index cac619f8..4dff1f8f 100644
--- a/bindings/R/rlibkriging/R/KrigingClass.R
+++ b/bindings/R/rlibkriging/R/KrigingClass.R
@@ -495,9 +495,9 @@ update.Kriging <- function(object, newy, newX, normalize = FALSE, ...) {
 #' @return The log-Likelihood computed for given
 #'     \eqn{\boldsymbol{theta}}{\theta}.
 #' 
-#' @method logLikelihood Kriging
+#' @method logLikelihoodFun Kriging
 #' @export 
-#' @aliases logLikelihood,Kriging,Kriging-method
+#' @aliases logLikelihoodFun,Kriging,Kriging-method
 #' 
 #' @examples
 #' f <- function(x) 1 - 1 / 2 * (sin(12 * x) / (1 + x) + 2 * cos(7 * x) * x^5 + 0.7)
@@ -506,12 +506,12 @@ update.Kriging <- function(object, newy, newX, normalize = FALSE, ...) {
 #' y <- f(X)
 #' r <- Kriging(y, X, kernel = "gauss")
 #' print(r)
-#' ll <- function(theta) logLikelihood(r, theta)$logLikelihood
+#' ll <- function(theta) logLikelihoodFun(r, theta)$logLikelihood
 #' t <- seq(from = 0.0001, to = 2, length.out = 101)
 #' plot(t, ll(t), type = 'l')
 #' abline(v = as.list(r)$theta, col = "blue")
 #' 
-logLikelihood.Kriging <- function(object, theta,
+logLikelihoodFun.Kriging <- function(object, theta,
                                   grad = FALSE, hess = FALSE, ...) {
   k <- kriging_model(object) 
   if (!is.matrix(theta)) theta <- matrix(theta, ncol = ncol(k$X))
@@ -524,7 +524,7 @@ logLikelihood.Kriging <- function(object, theta,
               logLikelihoodHess = array(NA, dim = c(nrow(theta), ncol(theta),
                                                     ncol(theta))))
   for (i in 1:nrow(theta)) {
-      ll <- kriging_logLikelihood(object, theta[i, ],
+      ll <- kriging_logLikelihoodFun(object, theta[i, ],
                                   grad = isTRUE(grad), hess = isTRUE(hess))
       out$logLikelihood[i] <- ll$logLikelihood
       if (isTRUE(grad)) out$logLikelihoodGrad[i, ] <- ll$logLikelihoodGrad
@@ -536,7 +536,35 @@ logLikelihood.Kriging <- function(object, theta,
   return(out)
 }
 
-## setMethod("logLikelihood", "Kriging", logLikelihood.Kriging)
+## ****************************************************************************
+#' Get Log-Likelihood of Kriging Model
+#' 
+#' @author Yann Richet \email{yann.richet@irsn.fr}
+#' 
+#' @param object An S3 Kriging object.
+#' @param ... Not used.
+#' 
+#' @return The log-Likelihood computed for fitted
+#'     \eqn{\boldsymbol{theta}}{\theta}.
+#' 
+#' @method logLikelihood Kriging
+#' @export 
+#' @aliases logLikelihood,Kriging,Kriging-method
+#' 
+#' @examples
+#' f <- function(x) 1 - 1 / 2 * (sin(12 * x) / (1 + x) + 2 * cos(7 * x) * x^5 + 0.7)
+#' set.seed(123)
+#' X <- as.matrix(runif(5))
+#' y <- f(X)
+#' r <- Kriging(y, X, kernel = "gauss")
+#' print(r)
+#' logLikelihood(r)
+#' 
+logLikelihood.Kriging <- function(object, ...) {
+  return(kriging_logLikelihood(object))
+}
+
+## setMethod("logLikelihood", "Kriging", logLikelihoodFun.Kriging)
 
 ## ****************************************************************************
 #' Compute Leave-One-Out (LOO) error for an object with S3 class
@@ -563,9 +591,9 @@ logLikelihood.Kriging <- function(object, theta,
 #' @return The leave-One-Out value computed for the given vector
 #'     \eqn{\boldsymbol{\theta}}{\theta} of correlation ranges.
 #' 
-#' @method leaveOneOut Kriging
+#' @method leaveOneOutFun Kriging
 #' @export 
-#' @aliases leaveOneOut,Kriging,Kriging-method
+#' @aliases leaveOneOutFun,Kriging,Kriging-method
 #' 
 #' @examples
 #' f <- function(x) 1 - 1 / 2 * (sin(12 * x) / (1 + x) + 2 * cos(7 * x) * x^5 + 0.7)
@@ -574,11 +602,11 @@ logLikelihood.Kriging <- function(object, theta,
 #' y <- f(X)
 #' r <- Kriging(y, X, kernel = "gauss", objective = "LOO")
 #' print(r)
-#' loo <-  function(theta) leaveOneOut(r,theta)$leaveOneOut
+#' loo <-  function(theta) leaveOneOutFun(r,theta)$leaveOneOut
 #' t <-  seq(from = 0.0001, to = 2, length.out = 101)
 #' plot(t, loo(t), type = "l")
 #' abline(v = as.list(r)$theta, col = "blue")
-leaveOneOut.Kriging <- function(object, theta, grad = FALSE, ...) {
+leaveOneOutFun.Kriging <- function(object, theta, grad = FALSE, ...) {
     k <- kriging_model(object) 
     if (!is.matrix(theta)) theta <- matrix(theta,ncol=ncol(k$X))
     if (ncol(theta) != ncol(k$X))
@@ -588,7 +616,7 @@ leaveOneOut.Kriging <- function(object, theta, grad = FALSE, ...) {
                 leaveOneOutGrad = matrix(NA, nrow = nrow(theta),
                                          ncol = ncol(theta)))
     for (i in 1:nrow(theta)) {
-        loo <- kriging_leaveOneOut(object,theta[i,], isTRUE(grad))
+        loo <- kriging_leaveOneOutFun(object,theta[i,], isTRUE(grad))
         out$leaveOneOut[i] <- loo$leaveOneOut
         if (isTRUE(grad)) out$leaveOneOutGrad[i, ] <- loo$leaveOneOutGrad
     }    
@@ -596,6 +624,34 @@ leaveOneOut.Kriging <- function(object, theta, grad = FALSE, ...) {
     return(out)
 }
 
+## ****************************************************************************
+#' Get leaveOneOut of Kriging Model
+#' 
+#' @author Yann Richet \email{yann.richet@irsn.fr}
+#' 
+#' @param object An S3 Kriging object.
+#' @param ... Not used.
+#' 
+#' @return The leaveOneOut computed for fitted
+#'     \eqn{\boldsymbol{theta}}{\theta}.
+#' 
+#' @method leaveOneOut Kriging
+#' @export 
+#' @aliases leaveOneOut,Kriging,Kriging-method
+#' 
+#' @examples
+#' f <- function(x) 1 - 1 / 2 * (sin(12 * x) / (1 + x) + 2 * cos(7 * x) * x^5 + 0.7)
+#' set.seed(123)
+#' X <- as.matrix(runif(5))
+#' y <- f(X)
+#' r <- Kriging(y, X, kernel = "gauss")
+#' print(r)
+#' leaveOneOut(r)
+#' 
+leaveOneOut.Kriging <- function(object, ...) {
+  return(kriging_leaveOneOut(object))
+}
+
 ##****************************************************************************
 #' Compute the log-marginal posterior of a kriging model, using the
 #' prior XXXY.
@@ -614,9 +670,9 @@ leaveOneOut.Kriging <- function(object, theta, grad = FALSE, ...) {
 #' @return The value of the log-marginal posterior computed for the
 #'     given vector theta.
 #' 
-#' @method logMargPost Kriging
+#' @method logMargPostFun Kriging
 #' @export 
-#' @aliases logMargPost,Kriging,Kriging-method
+#' @aliases logMargPostFun,Kriging,Kriging-method
 #'
 #' @references
 #' XXXY A reference describing the model (prior, ...)
@@ -630,11 +686,11 @@ leaveOneOut.Kriging <- function(object, theta, grad = FALSE, ...) {
 #' y <- f(X)
 #' r <- Kriging(y, X, "gauss")
 #' print(r)
-#' lmp <- function(theta) logMargPost(r, theta)$logMargPost
+#' lmp <- function(theta) logMargPostFun(r, theta)$logMargPost
 #' t <- seq(from = 0.0001, to = 2, length.out = 101)
 #' plot(t, lmp(t), type = "l")
 #' abline(v = as.list(r)$theta, col = "blue")
-logMargPost.Kriging <- function(object, theta, grad = FALSE, ...) {
+logMargPostFun.Kriging <- function(object, theta, grad = FALSE, ...) {
     k <- kriging_model(object) 
     if (!is.matrix(theta)) theta <- matrix(theta,ncol=ncol(k$X))
     if (ncol(theta) != ncol(k$X))
@@ -644,7 +700,7 @@ logMargPost.Kriging <- function(object, theta, grad = FALSE, ...) {
                 logMargPostGrad = matrix(NA, nrow = nrow(theta),
                                          ncol = ncol(theta)))
     for (i in 1:nrow(theta)) {
-        lmp <- kriging_logMargPost(object, theta[i, ], grad = isTRUE(grad))
+        lmp <- kriging_logMargPostFun(object, theta[i, ], grad = isTRUE(grad))
         out$logMargPost[i] <- lmp$logMargPost
         if (isTRUE(grad)) out$logMargPostGrad[i, ] <- lmp$logMargPostGrad
     }
@@ -652,36 +708,30 @@ logMargPost.Kriging <- function(object, theta, grad = FALSE, ...) {
     return(out)
 }
 
-##############################################################################
-
-
 ## ****************************************************************************
-#' Print the content of a \code{Kriging} object.
-#'
-#' @title Print a \code{Kriging} object
+#' Get logMargPost of Kriging Model
+#' 
 #' @author Yann Richet \email{yann.richet@irsn.fr}
-#'
-#' @param x A (S3) \code{Kriging} Object.
-#' @param ... Ignored.
-#'
-#' @return NULL
-#'
-#' @export
-#' @method print Kriging
+#' 
+#' @param object An S3 Kriging object.
+#' @param ... Not used.
+#' 
+#' @return The logMargPost computed for fitted
+#'     \eqn{\boldsymbol{theta}}{\theta}.
+#' 
+#' @method logMargPost Kriging
+#' @export 
+#' @aliases logMargPost,Kriging,Kriging-method
 #' 
 #' @examples
 #' f <- function(x) 1 - 1 / 2 * (sin(12 * x) / (1 + x) + 2 * cos(7 * x) * x^5 + 0.7)
 #' set.seed(123)
 #' X <- as.matrix(runif(5))
 #' y <- f(X)
-#' r <- Kriging(y, X, "gauss")
+#' r <- Kriging(y, X, kernel = "gauss")
 #' print(r)
-#' ## same thing
-#' r
-print.Kriging <- function(x, ...) {
-    if (length(list(...))>0) warning("Arguments ",paste0(names(list(...)),"=",list(...),collapse=",")," are ignored.")
-    k=kriging_model(x)
-    p = paste0("Kriging model:\n\n",kriging_summary(x),"\n")
-    cat(p)
-    ## return(p)
-}
\ No newline at end of file
+#' logMargPost(r)
+#' 
+logMargPost.Kriging <- function(object, ...) {
+  return(kriging_logMargPost(object))
+}
diff --git a/bindings/R/rlibkriging/R/NuggetKrigingClass.R b/bindings/R/rlibkriging/R/NuggetKrigingClass.R
index a904183a..9a79e58a 100644
--- a/bindings/R/rlibkriging/R/NuggetKrigingClass.R
+++ b/bindings/R/rlibkriging/R/NuggetKrigingClass.R
@@ -494,9 +494,9 @@ update.NuggetKriging <- function(object, newy, newX, normalize = FALSE, ...) {
 #' @return The log-Likelihood computed for given
 #'     \eqn{\boldsymbol{theta}}{\theta}.
 #' 
-#' @method logLikelihood NuggetKriging
+#' @method logLikelihoodFun NuggetKriging
 #' @export 
-#' @aliases logLikelihood,NuggetKriging,NuggetKriging-method
+#' @aliases logLikelihoodFun,NuggetKriging,NuggetKriging-method
 #' 
 #' @examples
 #' f <- function(x) 1 - 1 / 2 * (sin(12 * x) / (1 + x) + 2 * cos(7 * x) * x^5 + 0.7)
@@ -506,12 +506,12 @@ update.NuggetKriging <- function(object, newy, newX, normalize = FALSE, ...) {
 #' r <- NuggetKriging(y, X, kernel = "gauss")
 #' print(r)
 #' alpha = as.list(r)$sigma2/(as.list(r)$nugget+as.list(r)$sigma2)
-#' ll <- function(theta) logLikelihood(r, c(theta,alpha))$logLikelihood
+#' ll <- function(theta) logLikelihoodFun(r, c(theta,alpha))$logLikelihood
 #' t <- seq(from = 0.0001, to = 2, length.out = 101)
 #' plot(t, ll(t), type = 'l')
 #' abline(v = as.list(r)$theta, col = "blue")
 #' 
-logLikelihood.NuggetKriging <- function(object, theta_alpha,
+logLikelihoodFun.NuggetKriging <- function(object, theta_alpha,
                                   grad = FALSE, ...) {
   k <- nuggetkriging_model(object) 
   if (!is.matrix(theta_alpha)) theta_alpha <- matrix(theta_alpha, ncol = ncol(k$X)+1)
@@ -522,7 +522,7 @@ logLikelihood.NuggetKriging <- function(object, theta_alpha,
               logLikelihoodGrad = matrix(NA,nrow=nrow(theta_alpha),
                                          ncol = ncol(theta_alpha)))
   for (i in 1:nrow(theta_alpha)) {
-      ll <- nuggetkriging_logLikelihood(object, theta_alpha[i, ],
+      ll <- nuggetkriging_logLikelihoodFun(object, theta_alpha[i, ],
                                   grad = isTRUE(grad))
       out$logLikelihood[i] <- ll$logLikelihood
       if (isTRUE(grad)) out$logLikelihoodGrad[i, ] <- ll$logLikelihoodGrad
@@ -532,7 +532,34 @@ logLikelihood.NuggetKriging <- function(object, theta_alpha,
   return(out)
 }
 
-## setMethod("logLikelihood", "NuggetKriging", logLikelihood.NuggetKriging)
+
+## ****************************************************************************
+#' Get logLikelihood of NuggetKriging Model
+#' 
+#' @author Yann Richet \email{yann.richet@irsn.fr}
+#' 
+#' @param object An S3 NuggetKriging object.
+#' @param ... Not used.
+#' 
+#' @return The logLikelihood computed for fitted
+#'     \eqn{\boldsymbol{theta}}{\theta}.
+#' 
+#' @method logLikelihood NuggetKriging
+#' @export 
+#' @aliases logLikelihood,NuggetKriging,NuggetKriging-method
+#' 
+#' @examples
+#' f <- function(x) 1 - 1 / 2 * (sin(12 * x) / (1 + x) + 2 * cos(7 * x) * x^5 + 0.7)
+#' set.seed(123)
+#' X <- as.matrix(runif(5))
+#' y <- f(X)
+#' r <- Kriging(y, X, kernel = "gauss")
+#' print(r)
+#' logLikelihood(r)
+#' 
+logLikelihood.NuggetKriging <- function(object, ...) {
+  return(nuggetkriging_logLikelihood(object))
+}
 
 ##****************************************************************************
 #' Compute the log-marginal posterior of a kriging model, using the
@@ -552,9 +579,9 @@ logLikelihood.NuggetKriging <- function(object, theta_alpha,
 #' @return The value of the log-marginal posterior computed for the
 #'     given vector theta.
 #' 
-#' @method logMargPost NuggetKriging
+#' @method logMargPostFun NuggetKriging
 #' @export 
-#' @aliases logMargPost,NuggetKriging,NuggetKriging-method
+#' @aliases logMargPostFun,NuggetKriging,NuggetKriging-method
 #'
 #' @references
 #' XXXY A reference describing the model (prior, ...)
@@ -569,11 +596,11 @@ logLikelihood.NuggetKriging <- function(object, theta_alpha,
 #' r <- NuggetKriging(y, X, "gauss")
 #' print(r)
 #' alpha = as.list(r)$sigma2/(as.list(r)$nugget+as.list(r)$sigma2)
-#' lmp <- function(theta) logMargPost(r, c(theta,alpha))$logMargPost
+#' lmp <- function(theta) logMargPostFun(r, c(theta,alpha))$logMargPost
 #' t <- seq(from = 0.0001, to = 2, length.out = 101)
 #' plot(t, lmp(t), type = "l")
 #' abline(v = as.list(r)$theta, col = "blue")
-logMargPost.NuggetKriging <- function(object, theta_alpha, grad = FALSE, ...) {
+logMargPostFun.NuggetKriging <- function(object, theta_alpha, grad = FALSE, ...) {
     k <- nuggetkriging_model(object) 
     if (!is.matrix(theta_alpha)) theta_alpha <- matrix(theta_alpha,ncol=ncol(k$X)+1)
     if (ncol(theta_alpha) != ncol(k$X)+1)
@@ -583,7 +610,7 @@ logMargPost.NuggetKriging <- function(object, theta_alpha, grad = FALSE, ...) {
                 logMargPostGrad = matrix(NA, nrow = nrow(theta_alpha),
                                          ncol = ncol(theta_alpha)))
     for (i in 1:nrow(theta_alpha)) {
-        lmp <- nuggetkriging_logMargPost(object, theta_alpha[i, ], grad = isTRUE(grad))
+        lmp <- nuggetkriging_logMargPostFun(object, theta_alpha[i, ], grad = isTRUE(grad))
         out$logMargPost[i] <- lmp$logMargPost
         if (isTRUE(grad)) out$logMargPostGrad[i, ] <- lmp$logMargPostGrad
     }
@@ -591,6 +618,31 @@ logMargPost.NuggetKriging <- function(object, theta_alpha, grad = FALSE, ...) {
     return(out)
 }
 
-## setMethod("logMargPost", "NuggetKriging", logMargPost.NuggetKriging)
-
+## ****************************************************************************
+#' Get logMargPost of NuggetKriging Model
+#' 
+#' @author Yann Richet \email{yann.richet@irsn.fr}
+#' 
+#' @param object An S3 NuggetKriging object.
+#' @param ... Not used.
+#' 
+#' @return The logMargPost computed for fitted
+#'     \eqn{\boldsymbol{theta}}{\theta}.
+#' 
+#' @method logMargPost NuggetKriging
+#' @export 
+#' @aliases logMargPost,NuggetKriging,NuggetKriging-method
+#' 
+#' @examples
+#' f <- function(x) 1 - 1 / 2 * (sin(12 * x) / (1 + x) + 2 * cos(7 * x) * x^5 + 0.7)
+#' set.seed(123)
+#' X <- as.matrix(runif(5))
+#' y <- f(X)
+#' r <- Kriging(y, X, kernel = "gauss")
+#' print(r)
+#' logMargPost(r)
+#' 
+logMargPost.NuggetKriging <- function(object, ...) {
+  return(nuggetkriging_logMargPost(object))
+}
 
diff --git a/bindings/R/rlibkriging/R/allGenerics.R b/bindings/R/rlibkriging/R/allGenerics.R
index 94ea449a..8d4d77f6 100644
--- a/bindings/R/rlibkriging/R/allGenerics.R
+++ b/bindings/R/rlibkriging/R/allGenerics.R
@@ -9,14 +9,14 @@
 ## ****************************************************************************
 
 ## True generic
-setGeneric(name = "logLikelihood",
-           def = function(object, ...) standardGeneric("logLikelihood"))
+setGeneric(name = "logLikelihoodFun",
+           def = function(object, ...) standardGeneric("logLikelihoodFun"))
 
-setGeneric(name = "logMargPost",
-           def = function(object, ...) standardGeneric("logMargPost"))
+setGeneric(name = "logMargPostFun",
+           def = function(object, ...) standardGeneric("logMargPostFun"))
 
-setGeneric(name = "leaveOneOut",
-           def = function(object, ...) standardGeneric("leaveOneOut"))
+setGeneric(name = "leaveOneOutFun",
+           def = function(object, ...) standardGeneric("leaveOneOutFun"))
 
 ## *****************************************************************************
 ##' Coerce an object into an object with S4 class \code{"km"} from the
@@ -39,12 +39,67 @@ as.km <- function(x, ...) {
 
 ## *****************************************************************************
 ##' Compute the leave-One-Out error of a model given in \code{object},
-##' possibly at a different value of the parameters.
+##' at a different value of the parameters.
 ##'
-##' @title Compute Leave-One-Out 
+##' @title Leave-One-Out function
 ##' 
 ##' @param object An object representing a fitted model.
-##' @param ... Further arguments for methods.
+##' @param ... Further arguments of function (eg. range).
+##'
+##' @return The Leave-One-Out sum of squares.
+##' @export
+leaveOneOutFun <- function (object, ...) {
+    UseMethod("leaveOneOutFun")
+}
+
+## *****************************************************************************
+##' Compute the log-Likelihood of a model given in \code{object},
+##' at a different value of the parameters.
+##'
+##' @title Log-Likelihood function
+##'
+##' @param object An object representing a fitted model.
+##' @param ... Further arguments of function (eg. range).
+##'
+##' @return The log-likelihood.
+##' @export
+logLikelihoodFun <- function (object, ...) {
+    UseMethod("logLikelihoodFun")
+}
+
+## *****************************************************************************
+##' Compute the log-Marginal Posterior of a model given in
+##' \code{object}, at a different value of the parameters.
+##'
+##' @title log-Marginal Posterior function
+##'
+##' @param object An object representing a fitted model.
+##' @param ... Further arguments of function (eg. range).
+##'
+##' @return The log-marginal posterior.
+##' @export
+logMargPostFun <- function (object, ...) {
+    UseMethod("logMargPostFun")
+}
+
+
+## True generic
+setGeneric(name = "logLikelihood",
+           def = function(object, ...) standardGeneric("logLikelihood"))
+
+setGeneric(name = "logMargPost",
+           def = function(object, ...) standardGeneric("logMargPost"))
+
+setGeneric(name = "leaveOneOut",
+           def = function(object, ...) standardGeneric("leaveOneOut"))
+
+## *****************************************************************************
+##' Compute the leave-One-Out error of a model given in \code{object}.
+##'
+##' @title Compute Leave-One-Out
+##' 
+##' @param object An object representing a fitted model.
+##' @param ... Ignored.
 ##'
 ##' @return The Leave-One-Out sum of squares.
 ##' @export
@@ -53,13 +108,12 @@ leaveOneOut <- function (object, ...) {
 }
 
 ## *****************************************************************************
-##' Compute the log-Likelihood of a model given in \code{object},
-##' possibly at a different value of the parameters.
+##' Compute the log-Likelihood of a model given in \code{object}.
 ##'
 ##' @title Compute Log-Likelihood
 ##'
 ##' @param object An object representing a fitted model.
-##' @param ... Further arguments for methods.
+##' @param ... Ignored.
 ##'
 ##' @return The log-likelihood.
 ##' @export
@@ -69,12 +123,12 @@ logLikelihood <- function (object, ...) {
 
 ## *****************************************************************************
 ##' Compute the log-Marginal Posterior of a model given in
-##' \code{object}, possibly at a different value of the parameters.
+##' \code{object}.
 ##'
 ##' @title Compute log-Marginal Posterior
 ##'
 ##' @param object An object representing a fitted model.
-##' @param ... Further arguments for methods.
+##' @param ... Ignored.
 ##'
 ##' @return The log-marginal posterior.
 ##' @export
diff --git a/bindings/R/rlibkriging/man/get_num_nugget.Rd b/bindings/R/rlibkriging/man/get_num_nugget.Rd
new file mode 100644
index 00000000..183b02f3
--- /dev/null
+++ b/bindings/R/rlibkriging/man/get_num_nugget.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/miscUtils.R
+\name{get_num_nugget}
+\alias{get_num_nugget}
+\title{Get numerical nugget used in Cholesky decompositions}
+\usage{
+get_num_nugget(...)
+}
+\arguments{
+\item{...}{Ignored}
+}
+\value{
+numerical nugget used.
+}
+\description{
+Get numerical nugget used in Cholesky decompositions
+}
+\author{
+Yann Richet \email{yann.richet@irsn.fr}
+}
diff --git a/bindings/R/rlibkriging/man/leaveOneOut.Kriging.Rd b/bindings/R/rlibkriging/man/leaveOneOut.Kriging.Rd
index 6a1cf4cf..92e62349 100644
--- a/bindings/R/rlibkriging/man/leaveOneOut.Kriging.Rd
+++ b/bindings/R/rlibkriging/man/leaveOneOut.Kriging.Rd
@@ -3,47 +3,32 @@
 \name{leaveOneOut.Kriging}
 \alias{leaveOneOut.Kriging}
 \alias{leaveOneOut,Kriging,Kriging-method}
-\title{Leave-One-Out Method for the S3 class \code{"Kriging"}}
+\title{Get leaveOneOut of Kriging Model}
 \usage{
-\method{leaveOneOut}{Kriging}(object, theta, grad = FALSE, ...)
+\method{leaveOneOut}{Kriging}(object, ...)
 }
 \arguments{
-\item{object}{A \code{Kriging} object.}
-
-\item{theta}{A numeric vector of range parameters at which the LOO
-will be evaluated.}
-
-\item{grad}{Logical. Should the gradient (w.r.t. \code{theta}) be
-returned?}
+\item{object}{An S3 Kriging object.}
 
 \item{...}{Not used.}
 }
 \value{
-The leave-One-Out value computed for the given vector
-    \eqn{\boldsymbol{\theta}}{\theta} of correlation ranges.
+The leaveOneOut computed for fitted
+    \eqn{\boldsymbol{theta}}{\theta}.
 }
 \description{
-Compute Leave-One-Out (LOO) error for an object with S3 class
-\code{"Kriging"} representing a kriging model.
-}
-\details{
-The returned value is the sum of squares \eqn{\sum_{i=1}^n [y_i -
-\hat{y}_{i,(-i)}]^2} where \eqn{\hat{y}_{i,(-i)}} is the
-prediction of \eqn{y_i}{y[i]} based on the the observations \eqn{y_j}{y[j]}
-with \eqn{j \neq i}{j != i}.
+Get leaveOneOut of Kriging Model
 }
 \examples{
 f <- function(x) 1 - 1 / 2 * (sin(12 * x) / (1 + x) + 2 * cos(7 * x) * x^5 + 0.7)
 set.seed(123)
 X <- as.matrix(runif(5))
 y <- f(X)
-r <- Kriging(y, X, kernel = "gauss", objective = "LOO")
+r <- Kriging(y, X, kernel = "gauss")
 print(r)
-loo <-  function(theta) leaveOneOut(r,theta)$leaveOneOut
-t <-  seq(from = 0.0001, to = 2, length.out = 101)
-plot(t, loo(t), type = "l")
-abline(v = as.list(r)$theta, col = "blue")
+leaveOneOut(r)
+
 }
 \author{
-Yann Richet \email{yann.richet.irsn.fr}
+Yann Richet \email{yann.richet@irsn.fr}
 }
diff --git a/bindings/R/rlibkriging/man/leaveOneOut.Rd b/bindings/R/rlibkriging/man/leaveOneOut.Rd
index 2f66bf13..6eee30b0 100644
--- a/bindings/R/rlibkriging/man/leaveOneOut.Rd
+++ b/bindings/R/rlibkriging/man/leaveOneOut.Rd
@@ -9,12 +9,11 @@ leaveOneOut(object, ...)
 \arguments{
 \item{object}{An object representing a fitted model.}
 
-\item{...}{Further arguments for methods.}
+\item{...}{Ignored.}
 }
 \value{
 The Leave-One-Out sum of squares.
 }
 \description{
-Compute the leave-One-Out error of a model given in \code{object},
-possibly at a different value of the parameters.
+Compute the leave-One-Out error of a model given in \code{object}.
 }
diff --git a/bindings/R/rlibkriging/man/leaveOneOutFun.Kriging.Rd b/bindings/R/rlibkriging/man/leaveOneOutFun.Kriging.Rd
new file mode 100644
index 00000000..000df693
--- /dev/null
+++ b/bindings/R/rlibkriging/man/leaveOneOutFun.Kriging.Rd
@@ -0,0 +1,49 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/KrigingClass.R
+\name{leaveOneOutFun.Kriging}
+\alias{leaveOneOutFun.Kriging}
+\alias{leaveOneOutFun,Kriging,Kriging-method}
+\title{Leave-One-Out Method for the S3 class \code{"Kriging"}}
+\usage{
+\method{leaveOneOutFun}{Kriging}(object, theta, grad = FALSE, ...)
+}
+\arguments{
+\item{object}{A \code{Kriging} object.}
+
+\item{theta}{A numeric vector of range parameters at which the LOO
+will be evaluated.}
+
+\item{grad}{Logical. Should the gradient (w.r.t. \code{theta}) be
+returned?}
+
+\item{...}{Not used.}
+}
+\value{
+The leave-One-Out value computed for the given vector
+    \eqn{\boldsymbol{\theta}}{\theta} of correlation ranges.
+}
+\description{
+Compute Leave-One-Out (LOO) error for an object with S3 class
+\code{"Kriging"} representing a kriging model.
+}
+\details{
+The returned value is the sum of squares \eqn{\sum_{i=1}^n [y_i -
+\hat{y}_{i,(-i)}]^2} where \eqn{\hat{y}_{i,(-i)}} is the
+prediction of \eqn{y_i}{y[i]} based on the the observations \eqn{y_j}{y[j]}
+with \eqn{j \neq i}{j != i}.
+}
+\examples{
+f <- function(x) 1 - 1 / 2 * (sin(12 * x) / (1 + x) + 2 * cos(7 * x) * x^5 + 0.7)
+set.seed(123)
+X <- as.matrix(runif(5))
+y <- f(X)
+r <- Kriging(y, X, kernel = "gauss", objective = "LOO")
+print(r)
+loo <-  function(theta) leaveOneOutFun(r,theta)$leaveOneOut
+t <-  seq(from = 0.0001, to = 2, length.out = 101)
+plot(t, loo(t), type = "l")
+abline(v = as.list(r)$theta, col = "blue")
+}
+\author{
+Yann Richet \email{yann.richet.irsn.fr}
+}
diff --git a/bindings/R/rlibkriging/man/leaveOneOutFun.Rd b/bindings/R/rlibkriging/man/leaveOneOutFun.Rd
new file mode 100644
index 00000000..7bbfe129
--- /dev/null
+++ b/bindings/R/rlibkriging/man/leaveOneOutFun.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/allGenerics.R
+\name{leaveOneOutFun}
+\alias{leaveOneOutFun}
+\title{Leave-One-Out function}
+\usage{
+leaveOneOutFun(object, ...)
+}
+\arguments{
+\item{object}{An object representing a fitted model.}
+
+\item{...}{Further arguments of function (eg. range).}
+}
+\value{
+The Leave-One-Out sum of squares.
+}
+\description{
+Compute the leave-One-Out error of a model given in \code{object},
+at a different value of the parameters.
+}
diff --git a/bindings/R/rlibkriging/man/logLikelihood.Kriging.Rd b/bindings/R/rlibkriging/man/logLikelihood.Kriging.Rd
index 0f100b7b..907f3717 100644
--- a/bindings/R/rlibkriging/man/logLikelihood.Kriging.Rd
+++ b/bindings/R/rlibkriging/man/logLikelihood.Kriging.Rd
@@ -3,28 +3,21 @@
 \name{logLikelihood.Kriging}
 \alias{logLikelihood.Kriging}
 \alias{logLikelihood,Kriging,Kriging-method}
-\title{Compute Log-Likelihood of Kriging Model}
+\title{Get Log-Likelihood of Kriging Model}
 \usage{
-\method{logLikelihood}{Kriging}(object, theta, grad = FALSE, hess = FALSE, ...)
+\method{logLikelihood}{Kriging}(object, ...)
 }
 \arguments{
 \item{object}{An S3 Kriging object.}
 
-\item{theta}{A numeric vector of (positive) range parameters at
-which the log-likelihood will be evaluated.}
-
-\item{grad}{Logical. Should the function return the gradient?}
-
-\item{hess}{Logical. Should the function return Hessian?}
-
 \item{...}{Not used.}
 }
 \value{
-The log-Likelihood computed for given
+The log-Likelihood computed for fitted
     \eqn{\boldsymbol{theta}}{\theta}.
 }
 \description{
-Compute Log-Likelihood of Kriging Model
+Get Log-Likelihood of Kriging Model
 }
 \examples{
 f <- function(x) 1 - 1 / 2 * (sin(12 * x) / (1 + x) + 2 * cos(7 * x) * x^5 + 0.7)
@@ -33,10 +26,7 @@ X <- as.matrix(runif(5))
 y <- f(X)
 r <- Kriging(y, X, kernel = "gauss")
 print(r)
-ll <- function(theta) logLikelihood(r, theta)$logLikelihood
-t <- seq(from = 0.0001, to = 2, length.out = 101)
-plot(t, ll(t), type = 'l')
-abline(v = as.list(r)$theta, col = "blue")
+logLikelihood(r)
 
 }
 \author{
diff --git a/bindings/R/rlibkriging/man/logLikelihood.NuggetKriging.Rd b/bindings/R/rlibkriging/man/logLikelihood.NuggetKriging.Rd
index 2807d037..c55c62bd 100644
--- a/bindings/R/rlibkriging/man/logLikelihood.NuggetKriging.Rd
+++ b/bindings/R/rlibkriging/man/logLikelihood.NuggetKriging.Rd
@@ -3,40 +3,30 @@
 \name{logLikelihood.NuggetKriging}
 \alias{logLikelihood.NuggetKriging}
 \alias{logLikelihood,NuggetKriging,NuggetKriging-method}
-\title{Compute Log-Likelihood of NuggetKriging Model}
+\title{Get logLikelihood of NuggetKriging Model}
 \usage{
-\method{logLikelihood}{NuggetKriging}(object, theta, grad = FALSE, ...)
+\method{logLikelihood}{NuggetKriging}(object, ...)
 }
 \arguments{
 \item{object}{An S3 NuggetKriging object.}
 
-\item{theta}{A numeric vector of (positive) range parameters at
-which the log-likelihood will be evaluated.}
-
-\item{grad}{Logical. Should the function return the gradient?}
-
 \item{...}{Not used.}
-
-\item{hess}{Logical. Should the function return Hessian?}
 }
 \value{
-The log-Likelihood computed for given
+The logLikelihood computed for fitted
     \eqn{\boldsymbol{theta}}{\theta}.
 }
 \description{
-Compute Log-Likelihood of NuggetKriging Model
+Get logLikelihood of NuggetKriging Model
 }
 \examples{
 f <- function(x) 1 - 1 / 2 * (sin(12 * x) / (1 + x) + 2 * cos(7 * x) * x^5 + 0.7)
 set.seed(123)
 X <- as.matrix(runif(5))
-y <- f(X) + 0.01*rnorm(nrow(X))
-r <- NuggetKriging(y, X, kernel = "gauss")
+y <- f(X)
+r <- Kriging(y, X, kernel = "gauss")
 print(r)
-ll <- function(theta) logLikelihood(r, theta)$logLikelihood
-t <- seq(from = 0.0001, to = 2, length.out = 101)
-plot(t, ll(t), type = 'l')
-abline(v = as.list(r)$theta, col = "blue")
+logLikelihood(r)
 
 }
 \author{
diff --git a/bindings/R/rlibkriging/man/logLikelihood.Rd b/bindings/R/rlibkriging/man/logLikelihood.Rd
index 81383f2c..39891045 100644
--- a/bindings/R/rlibkriging/man/logLikelihood.Rd
+++ b/bindings/R/rlibkriging/man/logLikelihood.Rd
@@ -9,12 +9,11 @@ logLikelihood(object, ...)
 \arguments{
 \item{object}{An object representing a fitted model.}
 
-\item{...}{Further arguments for methods.}
+\item{...}{Ignored.}
 }
 \value{
 The log-likelihood.
 }
 \description{
-Compute the log-Likelihood of a model given in \code{object},
-possibly at a different value of the parameters.
+Compute the log-Likelihood of a model given in \code{object}.
 }
diff --git a/bindings/R/rlibkriging/man/logLikelihoodFun.Kriging.Rd b/bindings/R/rlibkriging/man/logLikelihoodFun.Kriging.Rd
new file mode 100644
index 00000000..f098cf75
--- /dev/null
+++ b/bindings/R/rlibkriging/man/logLikelihoodFun.Kriging.Rd
@@ -0,0 +1,44 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/KrigingClass.R
+\name{logLikelihoodFun.Kriging}
+\alias{logLikelihoodFun.Kriging}
+\alias{logLikelihoodFun,Kriging,Kriging-method}
+\title{Compute Log-Likelihood of Kriging Model}
+\usage{
+\method{logLikelihoodFun}{Kriging}(object, theta, grad = FALSE, hess = FALSE, ...)
+}
+\arguments{
+\item{object}{An S3 Kriging object.}
+
+\item{theta}{A numeric vector of (positive) range parameters at
+which the log-likelihood will be evaluated.}
+
+\item{grad}{Logical. Should the function return the gradient?}
+
+\item{hess}{Logical. Should the function return Hessian?}
+
+\item{...}{Not used.}
+}
+\value{
+The log-Likelihood computed for given
+    \eqn{\boldsymbol{theta}}{\theta}.
+}
+\description{
+Compute Log-Likelihood of Kriging Model
+}
+\examples{
+f <- function(x) 1 - 1 / 2 * (sin(12 * x) / (1 + x) + 2 * cos(7 * x) * x^5 + 0.7)
+set.seed(123)
+X <- as.matrix(runif(5))
+y <- f(X)
+r <- Kriging(y, X, kernel = "gauss")
+print(r)
+ll <- function(theta) logLikelihoodFun(r, theta)$logLikelihood
+t <- seq(from = 0.0001, to = 2, length.out = 101)
+plot(t, ll(t), type = 'l')
+abline(v = as.list(r)$theta, col = "blue")
+
+}
+\author{
+Yann Richet \email{yann.richet@irsn.fr}
+}
diff --git a/bindings/R/rlibkriging/man/logLikelihoodFun.NuggetKriging.Rd b/bindings/R/rlibkriging/man/logLikelihoodFun.NuggetKriging.Rd
new file mode 100644
index 00000000..0a302edb
--- /dev/null
+++ b/bindings/R/rlibkriging/man/logLikelihoodFun.NuggetKriging.Rd
@@ -0,0 +1,43 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/NuggetKrigingClass.R
+\name{logLikelihoodFun.NuggetKriging}
+\alias{logLikelihoodFun.NuggetKriging}
+\alias{logLikelihoodFun,NuggetKriging,NuggetKriging-method}
+\title{Compute Log-Likelihood of NuggetKriging Model}
+\usage{
+\method{logLikelihoodFun}{NuggetKriging}(object, theta_alpha, grad = FALSE, ...)
+}
+\arguments{
+\item{object}{An S3 NuggetKriging object.}
+
+\item{grad}{Logical. Should the function return the gradient?}
+
+\item{...}{Not used.}
+
+\item{theta}{A numeric vector of (positive) range parameters at
+which the log-likelihood will be evaluated.}
+}
+\value{
+The log-Likelihood computed for given
+    \eqn{\boldsymbol{theta}}{\theta}.
+}
+\description{
+Compute Log-Likelihood of NuggetKriging Model
+}
+\examples{
+f <- function(x) 1 - 1 / 2 * (sin(12 * x) / (1 + x) + 2 * cos(7 * x) * x^5 + 0.7)
+set.seed(123)
+X <- as.matrix(runif(5))
+y <- f(X) + 0.01*rnorm(nrow(X))
+r <- NuggetKriging(y, X, kernel = "gauss")
+print(r)
+alpha = as.list(r)$sigma2/(as.list(r)$nugget+as.list(r)$sigma2)
+ll <- function(theta) logLikelihoodFun(r, c(theta,alpha))$logLikelihood
+t <- seq(from = 0.0001, to = 2, length.out = 101)
+plot(t, ll(t), type = 'l')
+abline(v = as.list(r)$theta, col = "blue")
+
+}
+\author{
+Yann Richet \email{yann.richet@irsn.fr}
+}
diff --git a/bindings/R/rlibkriging/man/logLikelihoodFun.Rd b/bindings/R/rlibkriging/man/logLikelihoodFun.Rd
new file mode 100644
index 00000000..cc0c0074
--- /dev/null
+++ b/bindings/R/rlibkriging/man/logLikelihoodFun.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/allGenerics.R
+\name{logLikelihoodFun}
+\alias{logLikelihoodFun}
+\title{Log-Likelihood function}
+\usage{
+logLikelihoodFun(object, ...)
+}
+\arguments{
+\item{object}{An object representing a fitted model.}
+
+\item{...}{Further arguments of function (eg. range).}
+}
+\value{
+The log-likelihood.
+}
+\description{
+Compute the log-Likelihood of a model given in \code{object},
+at a different value of the parameters.
+}
diff --git a/bindings/R/rlibkriging/man/logMargPost.Kriging.Rd b/bindings/R/rlibkriging/man/logMargPost.Kriging.Rd
index 27141687..cdd3fe5a 100644
--- a/bindings/R/rlibkriging/man/logMargPost.Kriging.Rd
+++ b/bindings/R/rlibkriging/man/logMargPost.Kriging.Rd
@@ -3,47 +3,32 @@
 \name{logMargPost.Kriging}
 \alias{logMargPost.Kriging}
 \alias{logMargPost,Kriging,Kriging-method}
-\title{Compute Log-Marginal-Posterior of Kriging Model}
+\title{Get logMargPost of Kriging Model}
 \usage{
-\method{logMargPost}{Kriging}(object, theta, grad = FALSE, ...)
+\method{logMargPost}{Kriging}(object, ...)
 }
 \arguments{
-\item{object}{S3 Kriging object.}
-
-\item{theta}{Numeric vector of coorelation range parameters at
-which the function is to be evaluated.}
-
-\item{grad}{Logical. Should the function return the gradient
-(w.r.t theta)?}
+\item{object}{An S3 Kriging object.}
 
 \item{...}{Not used.}
 }
 \value{
-The value of the log-marginal posterior computed for the
-    given vector theta.
+The logMargPost computed for fitted
+    \eqn{\boldsymbol{theta}}{\theta}.
 }
 \description{
-Compute the log-marginal posterior of a kriging model, using the
-prior XXXY.
+Get logMargPost of Kriging Model
 }
 \examples{
 f <- function(x) 1 - 1 / 2 * (sin(12 * x) / (1 + x) + 2 * cos(7 * x) * x^5 + 0.7)
 set.seed(123)
 X <- as.matrix(runif(5))
 y <- f(X)
-r <- Kriging(y, X, "gauss")
+r <- Kriging(y, X, kernel = "gauss")
 print(r)
-lmp <- function(theta) logMargPost(r, theta)$logMargPost
-t <- seq(from = 0.0001, to = 2, length.out = 101)
-plot(t, lmp(t), type = "l")
-abline(v = as.list(r)$theta, col = "blue")
-}
-\references{
-XXXY A reference describing the model (prior, ...)
-}
-\seealso{
-\code{\link[RobustGaSP]{rgasp}} in the RobustGaSP package.
+logMargPost(r)
+
 }
 \author{
-Yann Richet \email{yann.richet.irsn.fr}
+Yann Richet \email{yann.richet@irsn.fr}
 }
diff --git a/bindings/R/rlibkriging/man/logMargPost.NuggetKriging.Rd b/bindings/R/rlibkriging/man/logMargPost.NuggetKriging.Rd
index 0dd0bbb0..2b96e706 100644
--- a/bindings/R/rlibkriging/man/logMargPost.NuggetKriging.Rd
+++ b/bindings/R/rlibkriging/man/logMargPost.NuggetKriging.Rd
@@ -3,47 +3,32 @@
 \name{logMargPost.NuggetKriging}
 \alias{logMargPost.NuggetKriging}
 \alias{logMargPost,NuggetKriging,NuggetKriging-method}
-\title{Compute Log-Marginal-Posterior of NuggetKriging Model}
+\title{Get logMargPost of NuggetKriging Model}
 \usage{
-\method{logMargPost}{NuggetKriging}(object, theta, grad = FALSE, ...)
+\method{logMargPost}{NuggetKriging}(object, ...)
 }
 \arguments{
-\item{object}{S3 NuggetKriging object.}
-
-\item{theta}{Numeric vector of coorelation range parameters at
-which the function is to be evaluated.}
-
-\item{grad}{Logical. Should the function return the gradient
-(w.r.t theta)?}
+\item{object}{An S3 NuggetKriging object.}
 
 \item{...}{Not used.}
 }
 \value{
-The value of the log-marginal posterior computed for the
-    given vector theta.
+The logMargPost computed for fitted
+    \eqn{\boldsymbol{theta}}{\theta}.
 }
 \description{
-Compute the log-marginal posterior of a kriging model, using the
-prior XXXY.
+Get logMargPost of NuggetKriging Model
 }
 \examples{
 f <- function(x) 1 - 1 / 2 * (sin(12 * x) / (1 + x) + 2 * cos(7 * x) * x^5 + 0.7)
 set.seed(123)
 X <- as.matrix(runif(5))
-y <- f(X) + 0.01*rnorm(nrow(X))
-r <- NuggetKriging(y, X, "gauss")
+y <- f(X)
+r <- Kriging(y, X, kernel = "gauss")
 print(r)
-lmp <- function(theta) logMargPost(r, theta)$logMargPost
-t <- seq(from = 0.0001, to = 2, length.out = 101)
-plot(t, lmp(t), type = "l")
-abline(v = as.list(r)$theta, col = "blue")
-}
-\references{
-XXXY A reference describing the model (prior, ...)
-}
-\seealso{
-\code{\link[RobustGaSP]{rgasp}} in the RobustGaSP package.
+logMargPost(r)
+
 }
 \author{
-Yann Richet \email{yann.richet.irsn.fr}
+Yann Richet \email{yann.richet@irsn.fr}
 }
diff --git a/bindings/R/rlibkriging/man/logMargPost.Rd b/bindings/R/rlibkriging/man/logMargPost.Rd
index 7a2b3622..74c3fa37 100644
--- a/bindings/R/rlibkriging/man/logMargPost.Rd
+++ b/bindings/R/rlibkriging/man/logMargPost.Rd
@@ -9,12 +9,12 @@ logMargPost(object, ...)
 \arguments{
 \item{object}{An object representing a fitted model.}
 
-\item{...}{Further arguments for methods.}
+\item{...}{Ignored.}
 }
 \value{
 The log-marginal posterior.
 }
 \description{
 Compute the log-Marginal Posterior of a model given in
-\code{object}, possibly at a different value of the parameters.
+\code{object}.
 }
diff --git a/bindings/R/rlibkriging/man/logMargPostFun.Kriging.Rd b/bindings/R/rlibkriging/man/logMargPostFun.Kriging.Rd
new file mode 100644
index 00000000..3eacb7c4
--- /dev/null
+++ b/bindings/R/rlibkriging/man/logMargPostFun.Kriging.Rd
@@ -0,0 +1,49 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/KrigingClass.R
+\name{logMargPostFun.Kriging}
+\alias{logMargPostFun.Kriging}
+\alias{logMargPostFun,Kriging,Kriging-method}
+\title{Compute Log-Marginal-Posterior of Kriging Model}
+\usage{
+\method{logMargPostFun}{Kriging}(object, theta, grad = FALSE, ...)
+}
+\arguments{
+\item{object}{S3 Kriging object.}
+
+\item{theta}{Numeric vector of coorelation range parameters at
+which the function is to be evaluated.}
+
+\item{grad}{Logical. Should the function return the gradient
+(w.r.t theta)?}
+
+\item{...}{Not used.}
+}
+\value{
+The value of the log-marginal posterior computed for the
+    given vector theta.
+}
+\description{
+Compute the log-marginal posterior of a kriging model, using the
+prior XXXY.
+}
+\examples{
+f <- function(x) 1 - 1 / 2 * (sin(12 * x) / (1 + x) + 2 * cos(7 * x) * x^5 + 0.7)
+set.seed(123)
+X <- as.matrix(runif(5))
+y <- f(X)
+r <- Kriging(y, X, "gauss")
+print(r)
+lmp <- function(theta) logMargPostFun(r, theta)$logMargPost
+t <- seq(from = 0.0001, to = 2, length.out = 101)
+plot(t, lmp(t), type = "l")
+abline(v = as.list(r)$theta, col = "blue")
+}
+\references{
+XXXY A reference describing the model (prior, ...)
+}
+\seealso{
+\code{\link[RobustGaSP]{rgasp}} in the RobustGaSP package.
+}
+\author{
+Yann Richet \email{yann.richet.irsn.fr}
+}
diff --git a/bindings/R/rlibkriging/man/logMargPostFun.NuggetKriging.Rd b/bindings/R/rlibkriging/man/logMargPostFun.NuggetKriging.Rd
new file mode 100644
index 00000000..b4a84d2d
--- /dev/null
+++ b/bindings/R/rlibkriging/man/logMargPostFun.NuggetKriging.Rd
@@ -0,0 +1,50 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/NuggetKrigingClass.R
+\name{logMargPostFun.NuggetKriging}
+\alias{logMargPostFun.NuggetKriging}
+\alias{logMargPostFun,NuggetKriging,NuggetKriging-method}
+\title{Compute Log-Marginal-Posterior of NuggetKriging Model}
+\usage{
+\method{logMargPostFun}{NuggetKriging}(object, theta_alpha, grad = FALSE, ...)
+}
+\arguments{
+\item{object}{S3 NuggetKriging object.}
+
+\item{grad}{Logical. Should the function return the gradient
+(w.r.t theta)?}
+
+\item{...}{Not used.}
+
+\item{theta}{Numeric vector of coorelation range parameters at
+which the function is to be evaluated.}
+}
+\value{
+The value of the log-marginal posterior computed for the
+    given vector theta.
+}
+\description{
+Compute the log-marginal posterior of a kriging model, using the
+prior XXXY.
+}
+\examples{
+f <- function(x) 1 - 1 / 2 * (sin(12 * x) / (1 + x) + 2 * cos(7 * x) * x^5 + 0.7)
+set.seed(123)
+X <- as.matrix(runif(5))
+y <- f(X) + 0.01*rnorm(nrow(X))
+r <- NuggetKriging(y, X, "gauss")
+print(r)
+alpha = as.list(r)$sigma2/(as.list(r)$nugget+as.list(r)$sigma2)
+lmp <- function(theta) logMargPostFun(r, c(theta,alpha))$logMargPost
+t <- seq(from = 0.0001, to = 2, length.out = 101)
+plot(t, lmp(t), type = "l")
+abline(v = as.list(r)$theta, col = "blue")
+}
+\references{
+XXXY A reference describing the model (prior, ...)
+}
+\seealso{
+\code{\link[RobustGaSP]{rgasp}} in the RobustGaSP package.
+}
+\author{
+Yann Richet \email{yann.richet.irsn.fr}
+}
diff --git a/bindings/R/rlibkriging/man/logMargPostFun.Rd b/bindings/R/rlibkriging/man/logMargPostFun.Rd
new file mode 100644
index 00000000..1a491a79
--- /dev/null
+++ b/bindings/R/rlibkriging/man/logMargPostFun.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/allGenerics.R
+\name{logMargPostFun}
+\alias{logMargPostFun}
+\title{log-Marginal Posterior function}
+\usage{
+logMargPostFun(object, ...)
+}
+\arguments{
+\item{object}{An object representing a fitted model.}
+
+\item{...}{Further arguments of function (eg. range).}
+}
+\value{
+The log-marginal posterior.
+}
+\description{
+Compute the log-Marginal Posterior of a model given in
+\code{object}, at a different value of the parameters.
+}
diff --git a/bindings/R/rlibkriging/man/reset_seed.Rd b/bindings/R/rlibkriging/man/reset_seed.Rd
new file mode 100644
index 00000000..a64298f8
--- /dev/null
+++ b/bindings/R/rlibkriging/man/reset_seed.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/miscUtils.R
+\name{reset_seed}
+\alias{reset_seed}
+\title{Set random ssed used in fit/optim}
+\usage{
+reset_seed(x, ...)
+}
+\arguments{
+\item{x}{Seed for mt19937 RNG}
+
+\item{...}{Ignored
+#'}
+}
+\description{
+Set random ssed used in fit/optim
+}
+\author{
+Yann Richet \email{yann.richet@irsn.fr}
+}
diff --git a/bindings/R/rlibkriging/man/set_num_nugget.Rd b/bindings/R/rlibkriging/man/set_num_nugget.Rd
new file mode 100644
index 00000000..b77fd948
--- /dev/null
+++ b/bindings/R/rlibkriging/man/set_num_nugget.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/miscUtils.R
+\name{set_num_nugget}
+\alias{set_num_nugget}
+\title{Setup numerical nugget used in Cholesky decompositions}
+\usage{
+set_num_nugget(x, ...)
+}
+\arguments{
+\item{x}{Numerical nugget value.}
+
+\item{...}{Ignored}
+}
+\description{
+Setup numerical nugget used in Cholesky decompositions
+}
+\author{
+Yann Richet \email{yann.richet@irsn.fr}
+}
diff --git a/bindings/R/rlibkriging/src/kriging_binding.cpp b/bindings/R/rlibkriging/src/kriging_binding.cpp
index 2cacf60c..01b6c32a 100644
--- a/bindings/R/rlibkriging/src/kriging_binding.cpp
+++ b/bindings/R/rlibkriging/src/kriging_binding.cpp
@@ -198,14 +198,14 @@ void kriging_update(Rcpp::List k, arma::vec y, arma::mat X, bool normalize = fal
 }
 
 // [[Rcpp::export]]
-Rcpp::List kriging_logLikelihood(Rcpp::List k, arma::vec theta, bool grad = false, bool hess = false) {
+Rcpp::List kriging_logLikelihoodFun(Rcpp::List k, arma::vec theta, bool grad = false, bool hess = false) {
   if (!k.inherits("Kriging"))
     Rcpp::stop("Input must be a Kriging object.");
   SEXP impl = k.attr("object");
 
   Rcpp::XPtr<Kriging> impl_ptr(impl);
 
-  std::tuple<double, arma::vec, arma::mat> ll = impl_ptr->logLikelihoodEval(theta, grad, hess);
+  std::tuple<double, arma::vec, arma::mat> ll = impl_ptr->logLikelihoodFun(theta, grad, hess);
   if (hess) {
     return Rcpp::List::create(Rcpp::Named("logLikelihood") = std::get<0>(ll),
                               Rcpp::Named("logLikelihoodGrad") = std::get<1>(ll),
@@ -219,14 +219,25 @@ Rcpp::List kriging_logLikelihood(Rcpp::List k, arma::vec theta, bool grad = fals
 }
 
 // [[Rcpp::export]]
-Rcpp::List kriging_leaveOneOut(Rcpp::List k, arma::vec theta, bool grad = false) {
+double kriging_logLikelihood(Rcpp::List k) {
   if (!k.inherits("Kriging"))
     Rcpp::stop("Input must be a Kriging object.");
   SEXP impl = k.attr("object");
 
   Rcpp::XPtr<Kriging> impl_ptr(impl);
 
-  std::tuple<double, arma::vec> loo = impl_ptr->leaveOneOutEval(theta, grad);
+  return impl_ptr->logLikelihood();
+}
+
+// [[Rcpp::export]]
+Rcpp::List kriging_leaveOneOutFun(Rcpp::List k, arma::vec theta, bool grad = false) {
+  if (!k.inherits("Kriging"))
+    Rcpp::stop("Input must be a Kriging object.");
+  SEXP impl = k.attr("object");
+
+  Rcpp::XPtr<Kriging> impl_ptr(impl);
+
+  std::tuple<double, arma::vec> loo = impl_ptr->leaveOneOutFun(theta, grad);
   if (grad) {
     return Rcpp::List::create(Rcpp::Named("leaveOneOut") = std::get<0>(loo),
                               Rcpp::Named("leaveOneOutGrad") = std::get<1>(loo));
@@ -236,14 +247,25 @@ Rcpp::List kriging_leaveOneOut(Rcpp::List k, arma::vec theta, bool grad = false)
 }
 
 // [[Rcpp::export]]
-Rcpp::List kriging_logMargPost(Rcpp::List k, arma::vec theta, bool grad = false) {
+double kriging_leaveOneOut(Rcpp::List k) {
   if (!k.inherits("Kriging"))
     Rcpp::stop("Input must be a Kriging object.");
   SEXP impl = k.attr("object");
 
   Rcpp::XPtr<Kriging> impl_ptr(impl);
 
-  std::tuple<double, arma::vec> lmp = impl_ptr->logMargPostEval(theta, grad);
+  return impl_ptr->leaveOneOut();
+}
+
+// [[Rcpp::export]]
+Rcpp::List kriging_logMargPostFun(Rcpp::List k, arma::vec theta, bool grad = false) {
+  if (!k.inherits("Kriging"))
+    Rcpp::stop("Input must be a Kriging object.");
+  SEXP impl = k.attr("object");
+
+  Rcpp::XPtr<Kriging> impl_ptr(impl);
+
+  std::tuple<double, arma::vec> lmp = impl_ptr->logMargPostFun(theta, grad);
   if (grad) {
     return Rcpp::List::create(Rcpp::Named("logMargPost") = std::get<0>(lmp),
                               Rcpp::Named("logMargPostGrad") = std::get<1>(lmp));
@@ -252,6 +274,17 @@ Rcpp::List kriging_logMargPost(Rcpp::List k, arma::vec theta, bool grad = false)
   }
 }
 
+// [[Rcpp::export]]
+double kriging_logMargPost(Rcpp::List k) {
+  if (!k.inherits("Kriging"))
+    Rcpp::stop("Input must be a Kriging object.");
+  SEXP impl = k.attr("object");
+
+  Rcpp::XPtr<Kriging> impl_ptr(impl);
+
+  return impl_ptr->logMargPost();
+}
+
 ///////////////////////////////////////////////////////////////////////////////////
 
 // [[Rcpp::export]]
diff --git a/bindings/R/rlibkriging/src/nuggetkriging_binding.cpp b/bindings/R/rlibkriging/src/nuggetkriging_binding.cpp
index ce412042..a5cd5623 100644
--- a/bindings/R/rlibkriging/src/nuggetkriging_binding.cpp
+++ b/bindings/R/rlibkriging/src/nuggetkriging_binding.cpp
@@ -220,14 +220,14 @@ void nuggetkriging_update(Rcpp::List k, arma::vec y, arma::mat X, bool normalize
 }
 
 // [[Rcpp::export]]
-Rcpp::List nuggetkriging_logLikelihood(Rcpp::List k, arma::vec theta_alpha, bool grad = false) {
+Rcpp::List nuggetkriging_logLikelihoodFun(Rcpp::List k, arma::vec theta_alpha, bool grad = false) {
   if (!k.inherits("NuggetKriging"))
     Rcpp::stop("Input must be a NuggetKriging object.");
   SEXP impl = k.attr("object");
 
   Rcpp::XPtr<NuggetKriging> impl_ptr(impl);
 
-  std::tuple<double, arma::vec> ll = impl_ptr->logLikelihoodEval(theta_alpha, grad);
+  std::tuple<double, arma::vec> ll = impl_ptr->logLikelihoodFun(theta_alpha, grad);
   if (grad) {
     return Rcpp::List::create(Rcpp::Named("logLikelihood") = std::get<0>(ll),
                               Rcpp::Named("logLikelihoodGrad") = std::get<1>(ll));
@@ -237,14 +237,25 @@ Rcpp::List nuggetkriging_logLikelihood(Rcpp::List k, arma::vec theta_alpha, bool
 }
 
 // [[Rcpp::export]]
-Rcpp::List nuggetkriging_logMargPost(Rcpp::List k, arma::vec theta, bool grad = false) {
+double nuggetkriging_logLikelihood(Rcpp::List k) {
   if (!k.inherits("NuggetKriging"))
     Rcpp::stop("Input must be a NuggetKriging object.");
   SEXP impl = k.attr("object");
 
   Rcpp::XPtr<NuggetKriging> impl_ptr(impl);
 
-  std::tuple<double, arma::vec> lmp = impl_ptr->logMargPostEval(theta, grad);
+  return impl_ptr->logLikelihood();
+}
+
+// [[Rcpp::export]]
+Rcpp::List nuggetkriging_logMargPostFun(Rcpp::List k, arma::vec theta, bool grad = false) {
+  if (!k.inherits("NuggetKriging"))
+    Rcpp::stop("Input must be a NuggetKriging object.");
+  SEXP impl = k.attr("object");
+
+  Rcpp::XPtr<NuggetKriging> impl_ptr(impl);
+
+  std::tuple<double, arma::vec> lmp = impl_ptr->logMargPostFun(theta, grad);
   if (grad) {
     return Rcpp::List::create(Rcpp::Named("logMargPost") = std::get<0>(lmp),
                               Rcpp::Named("logMargPostGrad") = std::get<1>(lmp));
@@ -253,6 +264,17 @@ Rcpp::List nuggetkriging_logMargPost(Rcpp::List k, arma::vec theta, bool grad =
   }
 }
 
+// [[Rcpp::export]]
+double nuggetkriging_logMargPost(Rcpp::List k) {
+  if (!k.inherits("NuggetKriging"))
+    Rcpp::stop("Input must be a NuggetKriging object.");
+  SEXP impl = k.attr("object");
+
+  Rcpp::XPtr<NuggetKriging> impl_ptr(impl);
+
+  return impl_ptr->logMargPost();
+}
+
 ///////////////////////////////////////////////////////////////////////////////////
 
 // [[Rcpp::export]]
diff --git a/src/lib/Bench.cpp b/src/lib/Bench.cpp
index 42cc5be0..bda06634 100644
--- a/src/lib/Bench.cpp
+++ b/src/lib/Bench.cpp
@@ -114,7 +114,7 @@ double Bench::LogLik(Kriging& ok, const arma::vec& theta) {
 
   double s = 0;
   for (int i = 0; i < n; i++) {
-    s += ok.logLikelihood(theta, nullptr, nullptr, &okm_data);
+    s += ok._logLikelihood(theta, nullptr, nullptr, &okm_data);
   }
   return s / n;
 }
@@ -133,7 +133,7 @@ arma::vec Bench::LogLikGrad(Kriging& ok, const arma::vec& theta) {
 
   arma::vec s = arma::zeros(theta.n_elem);
   for (int i = 0; i < n; i++) {
-    s += ok.logLikelihood(theta, &grad, nullptr, &okm_data);
+    s += ok._logLikelihood(theta, &grad, nullptr, &okm_data);
   }
   return s / n;
 }
diff --git a/src/lib/Kriging.cpp b/src/lib/Kriging.cpp
index 3baf2e19..43d39a9a 100644
--- a/src/lib/Kriging.cpp
+++ b/src/lib/Kriging.cpp
@@ -81,10 +81,10 @@ LIBKRIGING_EXPORT Kriging::Kriging(const arma::colvec& y,
 
 // Objective function for fit : -logLikelihood
 
-double Kriging::logLikelihood(const arma::vec& _theta,
-                              arma::vec* grad_out,
-                              arma::mat* hess_out,
-                              Kriging::OKModel* okm_data) const {
+double Kriging::_logLikelihood(const arma::vec& _theta,
+                               arma::vec* grad_out,
+                               arma::mat* hess_out,
+                               Kriging::OKModel* okm_data) const {
   // arma::cout << " theta: " << _theta << arma::endl;
   //' @ref https://github.com/cran/DiceKriging/blob/master/R/logLikFun.R
   //  model@covariance <- vect2covparam(model@covariance, param)
@@ -301,9 +301,9 @@ double Kriging::logLikelihood(const arma::vec& _theta,
   return ll;
 }
 
-LIBKRIGING_EXPORT std::tuple<double, arma::vec, arma::mat> Kriging::logLikelihoodEval(const arma::vec& _theta,
-                                                                                      const bool _grad,
-                                                                                      const bool _hess) {
+LIBKRIGING_EXPORT std::tuple<double, arma::vec, arma::mat> Kriging::logLikelihoodFun(const arma::vec& _theta,
+                                                                                     const bool _grad,
+                                                                                     const bool _hess) {
   arma::mat T;
   arma::mat M;
   arma::colvec z;
@@ -317,13 +317,13 @@ LIBKRIGING_EXPORT std::tuple<double, arma::vec, arma::mat> Kriging::logLikelihoo
   if (_grad || _hess) {
     grad = arma::vec(_theta.n_elem);
     if (!_hess) {
-      ll = logLikelihood(_theta, &grad, nullptr, &okm_data);
+      ll = _logLikelihood(_theta, &grad, nullptr, &okm_data);
     } else {
       hess = arma::mat(_theta.n_elem, _theta.n_elem);
-      ll = logLikelihood(_theta, &grad, &hess, &okm_data);
+      ll = _logLikelihood(_theta, &grad, &hess, &okm_data);
     }
   } else
-    ll = logLikelihood(_theta, nullptr, nullptr, &okm_data);
+    ll = _logLikelihood(_theta, nullptr, nullptr, &okm_data);
 
   return std::make_tuple(ll, std::move(grad), std::move(hess));
 }
@@ -339,7 +339,7 @@ arma::colvec DiagABA(const arma::mat& A, const arma::mat& B) {
   return c;
 }
 
-double Kriging::leaveOneOut(const arma::vec& _theta, arma::vec* grad_out, Kriging::OKModel* okm_data) const {
+double Kriging::_leaveOneOut(const arma::vec& _theta, arma::vec* grad_out, Kriging::OKModel* okm_data) const {
   // arma::cout << " theta: " << _theta << arma::endl;
   //' @ref https://github.com/DiceKrigingClub/DiceKriging/blob/master/R/leaveOneOutFun.R
   // model@covariance <- vect2covparam(model@covariance, param)
@@ -463,7 +463,7 @@ double Kriging::leaveOneOut(const arma::vec& _theta, arma::vec* grad_out, Krigin
   return loo;
 }
 
-LIBKRIGING_EXPORT std::tuple<double, arma::vec> Kriging::leaveOneOutEval(const arma::vec& _theta, const bool _grad) {
+LIBKRIGING_EXPORT std::tuple<double, arma::vec> Kriging::leaveOneOutFun(const arma::vec& _theta, const bool _grad) {
   arma::mat T;
   arma::mat M;
   arma::colvec z;
@@ -475,16 +475,16 @@ LIBKRIGING_EXPORT std::tuple<double, arma::vec> Kriging::leaveOneOutEval(const a
   arma::vec grad;
   if (_grad) {
     grad = arma::vec(_theta.n_elem);
-    loo = leaveOneOut(_theta, &grad, &okm_data);
+    loo = _leaveOneOut(_theta, &grad, &okm_data);
   } else
-    loo = leaveOneOut(_theta, nullptr, &okm_data);
+    loo = _leaveOneOut(_theta, nullptr, &okm_data);
 
   return std::make_tuple(loo, std::move(grad));
 }
 
 // Objective function for fit: bayesian-like approach fromm RobustGaSP
 
-double Kriging::logMargPost(const arma::vec& _theta, arma::vec* grad_out, Kriging::OKModel* okm_data) const {
+double Kriging::_logMargPost(const arma::vec& _theta, arma::vec* grad_out, Kriging::OKModel* okm_data) const {
   // arma::cout << " theta: " << _theta << arma::endl;
 
   // In RobustGaSP:
@@ -693,7 +693,7 @@ double Kriging::logMargPost(const arma::vec& _theta, arma::vec* grad_out, Krigin
   return (log_marginal_lik + log_approx_ref_prior);
 }
 
-LIBKRIGING_EXPORT std::tuple<double, arma::vec> Kriging::logMargPostEval(const arma::vec& _theta, const bool _grad) {
+LIBKRIGING_EXPORT std::tuple<double, arma::vec> Kriging::logMargPostFun(const arma::vec& _theta, const bool _grad) {
   arma::mat T;
   arma::mat M;
   arma::colvec z;
@@ -705,13 +705,25 @@ LIBKRIGING_EXPORT std::tuple<double, arma::vec> Kriging::logMargPostEval(const a
   arma::vec grad;
   if (_grad) {
     grad = arma::vec(_theta.n_elem);
-    lmp = logMargPost(_theta, &grad, &okm_data);
+    lmp = _logMargPost(_theta, &grad, &okm_data);
   } else
-    lmp = logMargPost(_theta, nullptr, &okm_data);
+    lmp = _logMargPost(_theta, nullptr, &okm_data);
 
   return std::make_tuple(lmp, std::move(grad));
 }
 
+LIBKRIGING_EXPORT double Kriging::logLikelihood() {
+  return std::get<0>(Kriging::logLikelihoodFun(m_theta, false, false));
+}
+
+LIBKRIGING_EXPORT double Kriging::leaveOneOut() {
+  return std::get<0>(Kriging::leaveOneOutFun(m_theta, false));
+}
+
+LIBKRIGING_EXPORT double Kriging::logMargPost() {
+  return std::get<0>(Kriging::logMargPostFun(m_theta, false));
+}
+
 double optim_newton(std::function<double(arma::vec& x, arma::vec* grad_out, arma::mat* hess_out)> f,
                     arma::vec& x_0,
                     const arma::vec& x_lower,
@@ -834,7 +846,7 @@ LIBKRIGING_EXPORT void Kriging::fit(const arma::colvec& y,
         [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);
+          double ll = this->_logLikelihood(_theta, grad_out, hess_out, okm_data);
           if (grad_out != nullptr)
             *grad_out = *grad_out % _theta;
           if (hess_out != nullptr)
@@ -847,7 +859,7 @@ LIBKRIGING_EXPORT void Kriging::fit(const arma::colvec& y,
         [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);
+          double loo = this->_leaveOneOut(_theta, grad_out, okm_data);
           if (grad_out != nullptr)
             *grad_out = -*grad_out % _theta;
           return loo;
@@ -860,7 +872,7 @@ LIBKRIGING_EXPORT void Kriging::fit(const arma::colvec& y,
         [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);
+          double lmp = this->_logMargPost(_theta, grad_out, okm_data);
           if (grad_out != nullptr)
             *grad_out = *grad_out % _theta;
           return -lmp;
diff --git a/src/lib/NuggetKriging.cpp b/src/lib/NuggetKriging.cpp
index cef8244a..64b6ab7b 100644
--- a/src/lib/NuggetKriging.cpp
+++ b/src/lib/NuggetKriging.cpp
@@ -80,9 +80,9 @@ LIBKRIGING_EXPORT NuggetKriging::NuggetKriging(const arma::colvec& y,
 
 // Objective function for fit : -logLikelihood
 
-double NuggetKriging::logLikelihood(const arma::vec& _theta_alpha,
-                                    arma::vec* grad_out,
-                                    NuggetKriging::OKModel* okm_data) const {
+double NuggetKriging::_logLikelihood(const arma::vec& _theta_alpha,
+                                     arma::vec* grad_out,
+                                     NuggetKriging::OKModel* okm_data) const {
   // arma::cout << " theta, alpha: " << _theta_alpha.head(m_X.n_cols) << "," << _theta_alpha.at(m_X.n_cols) <<
   // arma::endl;
   //' @ref https://github.com/cran/DiceKriging/blob/master/R/logLikFun.R
@@ -217,8 +217,8 @@ double NuggetKriging::logLikelihood(const arma::vec& _theta_alpha,
   return ll;
 }
 
-LIBKRIGING_EXPORT std::tuple<double, arma::vec> NuggetKriging::logLikelihoodEval(const arma::vec& _theta_alpha,
-                                                                                 const bool _grad) {
+LIBKRIGING_EXPORT std::tuple<double, arma::vec> NuggetKriging::logLikelihoodFun(const arma::vec& _theta_alpha,
+                                                                                const bool _grad) {
   arma::mat T;
   arma::mat M;
   arma::colvec z;
@@ -232,17 +232,17 @@ LIBKRIGING_EXPORT std::tuple<double, arma::vec> NuggetKriging::logLikelihoodEval
   arma::vec grad;
   if (_grad) {
     grad = arma::vec(_theta_alpha.n_elem);
-    ll = logLikelihood(_theta_alpha, &grad, &okm_data);
+    ll = _logLikelihood(_theta_alpha, &grad, &okm_data);
   } else
-    ll = logLikelihood(_theta_alpha, nullptr, &okm_data);
+    ll = _logLikelihood(_theta_alpha, nullptr, &okm_data);
 
   return std::make_tuple(ll, std::move(grad));
 }
 // Objective function for fit: bayesian-like approach fromm RobustGaSP
 
-double NuggetKriging::logMargPost(const arma::vec& _theta_alpha,
-                                  arma::vec* grad_out,
-                                  NuggetKriging::OKModel* okm_data) const {
+double NuggetKriging::_logMargPost(const arma::vec& _theta_alpha,
+                                   arma::vec* grad_out,
+                                   NuggetKriging::OKModel* okm_data) const {
   // arma::cout << " theta: " << _theta << arma::endl;
 
   // In RobustGaSP:
@@ -440,8 +440,8 @@ double NuggetKriging::logMargPost(const arma::vec& _theta_alpha,
   return (log_marginal_lik + log_approx_ref_prior);
 }
 
-LIBKRIGING_EXPORT std::tuple<double, arma::vec> NuggetKriging::logMargPostEval(const arma::vec& _theta,
-                                                                               const bool _grad) {
+LIBKRIGING_EXPORT std::tuple<double, arma::vec> NuggetKriging::logMargPostFun(const arma::vec& _theta,
+                                                                              const bool _grad) {
   arma::mat T;
   arma::mat M;
   arma::colvec z;
@@ -453,13 +453,29 @@ LIBKRIGING_EXPORT std::tuple<double, arma::vec> NuggetKriging::logMargPostEval(c
   arma::vec grad;
   if (_grad) {
     grad = arma::vec(_theta.n_elem);
-    lmp = logMargPost(_theta, &grad, &okm_data);
+    lmp = _logMargPost(_theta, &grad, &okm_data);
   } else
-    lmp = logMargPost(_theta, nullptr, &okm_data);
+    lmp = _logMargPost(_theta, nullptr, &okm_data);
 
   return std::make_tuple(lmp, std::move(grad));
 }
 
+LIBKRIGING_EXPORT double NuggetKriging::logLikelihood() {
+  int d = m_theta.n_elem;
+  arma::vec _theta_alpha = arma::vec(d + 1);
+  _theta_alpha.head(d) = m_theta;
+  _theta_alpha.at(d) = m_sigma2 / (m_sigma2 + m_nugget);
+  return std::get<0>(NuggetKriging::logLikelihoodFun(_theta_alpha, false));
+}
+
+LIBKRIGING_EXPORT double NuggetKriging::logMargPost() {
+  int d = m_theta.n_elem;
+  arma::vec _theta_alpha = arma::vec(d + 1);
+  _theta_alpha.head(d) = m_theta;
+  _theta_alpha.at(d) = m_sigma2 / (m_sigma2 + m_nugget);
+  return std::get<0>(NuggetKriging::logMargPostFun(_theta_alpha, false));
+}
+
 /** Fit the kriging object on (X,y):
  * @param y is n length column vector of output
  * @param X is n*d matrix of input
@@ -490,7 +506,7 @@ LIBKRIGING_EXPORT void NuggetKriging::fit(const arma::colvec& y,
       // Change variable for opt: . -> 1/exp(.)
       arma::vec _theta_alpha = 1 / arma::exp(_gamma);
       _theta_alpha[_theta_alpha.n_elem - 1] = _gamma[_theta_alpha.n_elem - 1];  // opt!
-      double ll = this->logLikelihood(_theta_alpha, grad_out, okm_data);
+      double ll = this->_logLikelihood(_theta_alpha, grad_out, okm_data);
       if (grad_out != nullptr) {
         (*grad_out).head(_theta_alpha.n_elem - 1) %= _theta_alpha.head(_theta_alpha.n_elem - 1);
         (*grad_out).at(_theta_alpha.n_elem - 1) *= -1;
@@ -504,7 +520,7 @@ LIBKRIGING_EXPORT void NuggetKriging::fit(const arma::colvec& y,
       // Change variable for opt: . -> 1/exp(.)
       arma::vec _theta_alpha = 1 / arma::exp(_gamma);
       _theta_alpha[_theta_alpha.n_elem - 1] = _gamma[_theta_alpha.n_elem - 1];  // opt!
-      double lmp = this->logMargPost(_theta_alpha, grad_out, okm_data);
+      double lmp = this->_logMargPost(_theta_alpha, grad_out, okm_data);
       if (grad_out != nullptr) {
         (*grad_out).head(_theta_alpha.n_elem - 1) %= _theta_alpha.head(_theta_alpha.n_elem - 1);
         (*grad_out).at(_theta_alpha.n_elem - 1) *= -1;
diff --git a/src/lib/include/libKriging/Kriging.hpp b/src/lib/include/libKriging/Kriging.hpp
index 1cf872ca..9421d297 100644
--- a/src/lib/include/libKriging/Kriging.hpp
+++ b/src/lib/include/libKriging/Kriging.hpp
@@ -109,12 +109,12 @@ class Kriging {
     bool is_sigma2_estim;
   };
 
-  double logLikelihood(const arma::vec& _theta,
-                       arma::vec* grad_out,
-                       arma::mat* hess_out,
-                       Kriging::OKModel* okm_data) const;
-  double leaveOneOut(const arma::vec& _theta, arma::vec* grad_out, Kriging::OKModel* okm_data) const;
-  double logMargPost(const arma::vec& _theta, arma::vec* grad_out, Kriging::OKModel* okm_data) const;
+  double _logLikelihood(const arma::vec& _theta,
+                        arma::vec* grad_out,
+                        arma::mat* hess_out,
+                        Kriging::OKModel* okm_data) const;
+  double _leaveOneOut(const arma::vec& _theta, arma::vec* grad_out, Kriging::OKModel* okm_data) const;
+  double _logMargPost(const arma::vec& _theta, arma::vec* grad_out, Kriging::OKModel* okm_data) const;
 
   // at least, just call make_dist(kernel)
   LIBKRIGING_EXPORT Kriging(const std::string& covType);
@@ -145,12 +145,16 @@ class Kriging {
                              const Parameters& parameters
                              = Parameters{-1, false, true, arma::mat(), false, true, arma::vec(), false, true});
 
-  LIBKRIGING_EXPORT std::tuple<double, arma::vec, arma::mat> logLikelihoodEval(const arma::vec& theta,
-                                                                               const bool grad,
-                                                                               const bool hess);
-  LIBKRIGING_EXPORT std::tuple<double, arma::vec> leaveOneOutEval(const arma::vec& theta, const bool grad);
+  LIBKRIGING_EXPORT std::tuple<double, arma::vec, arma::mat> logLikelihoodFun(const arma::vec& theta,
+                                                                              const bool grad,
+                                                                              const bool hess);
+  LIBKRIGING_EXPORT std::tuple<double, arma::vec> leaveOneOutFun(const arma::vec& theta, const bool grad);
 
-  LIBKRIGING_EXPORT std::tuple<double, arma::vec> logMargPostEval(const arma::vec& theta, const bool grad);
+  LIBKRIGING_EXPORT std::tuple<double, arma::vec> logMargPostFun(const arma::vec& theta, const bool grad);
+
+  LIBKRIGING_EXPORT double logLikelihood();
+  LIBKRIGING_EXPORT double leaveOneOut();
+  LIBKRIGING_EXPORT double logMargPost();
 
   /** Compute the prediction for given points X'
    * @param Xp is m*d matrix of points where to predict output
diff --git a/src/lib/include/libKriging/NuggetKriging.hpp b/src/lib/include/libKriging/NuggetKriging.hpp
index 427f6785..786726d9 100644
--- a/src/lib/include/libKriging/NuggetKriging.hpp
+++ b/src/lib/include/libKriging/NuggetKriging.hpp
@@ -135,8 +135,8 @@ class NuggetKriging {
     double var;
   };
 
-  double logLikelihood(const arma::vec& _theta, arma::vec* grad_out, NuggetKriging::OKModel* okm_data) const;
-  double logMargPost(const arma::vec& _theta, arma::vec* grad_out, NuggetKriging::OKModel* okm_data) const;
+  double _logLikelihood(const arma::vec& _theta, arma::vec* grad_out, NuggetKriging::OKModel* okm_data) const;
+  double _logMargPost(const arma::vec& _theta, arma::vec* grad_out, NuggetKriging::OKModel* okm_data) const;
 
   // at least, just call make_dist(kernel)
   LIBKRIGING_EXPORT NuggetKriging(const std::string& covType);
@@ -166,9 +166,12 @@ class NuggetKriging {
                              const std::string& objective = "LL",
                              const Parameters& parameters = Parameters{});
 
-  LIBKRIGING_EXPORT std::tuple<double, arma::vec> logLikelihoodEval(const arma::vec& theta, const bool grad);
+  LIBKRIGING_EXPORT std::tuple<double, arma::vec> logLikelihoodFun(const arma::vec& theta, const bool grad);
 
-  LIBKRIGING_EXPORT std::tuple<double, arma::vec> logMargPostEval(const arma::vec& theta, const bool grad);
+  LIBKRIGING_EXPORT std::tuple<double, arma::vec> logMargPostFun(const arma::vec& theta, const bool grad);
+
+  LIBKRIGING_EXPORT double logLikelihood();
+  LIBKRIGING_EXPORT double logMargPost();
 
   /** Compute the prediction for given points X'
    * @param Xp is m*d matrix of points where to predict output


Test patch:
diff --git a/bindings/Octave/tests/binding_consistency_test.m b/bindings/Octave/tests/binding_consistency_test.m
index 6332a4de..3000240c 100644
--- a/bindings/Octave/tests/binding_consistency_test.m
+++ b/bindings/Octave/tests/binding_consistency_test.m
@@ -25,11 +25,11 @@
     k_m = Kriging(y, X, kernel, "constant", false, "BFGS", "LL") % lk.Parameters() % FIXME parameters arg not mapped
     s = size(X);
     x = 0.3 * ones(s(2), 1);
-    [loo, loograd] = k_m.leaveOneOut(x, true);
+    [loo, loograd] = k_m.leaveOneOutFun(x, true);
     assert(relative_error(loo, loo_ref) < tolerance)
     assert(relative_error(loograd, loograd_ref) < tolerance)
 
-    [ll, llgrad] = k_m.logLikelihood(x, true, false); % flags are optional; then results are driven by the output
+    [ll, llgrad] = k_m.logLikelihoodFun(x, true, false); % flags are optional; then results are driven by the output
     assert(relative_error(ll, ll_ref) < tolerance)
     assert(relative_error(llgrad, llgrad_ref) < tolerance)
 
@@ -59,7 +59,7 @@
         s = size(X);
         x = 0.3 * ones(s(2), 1);
 
-        [ll, llgrad] = k_m.logLikelihood(x, true, false);
+        [ll, llgrad] = k_m.logLikelihoodFun(x, true, false);
         assert(relative_error(ll, ll_ref) < tolerance)
         assert(relative_error(llgrad, llgrad_ref) < tolerance)
 
diff --git a/bindings/Python/tests/binding_consistency_test.py b/bindings/Python/tests/binding_consistency_test.py
index 499c4064..b52a97b1 100644
--- a/bindings/Python/tests/binding_consistency_test.py
+++ b/bindings/Python/tests/binding_consistency_test.py
@@ -63,11 +63,11 @@ def test_data1():
     r = lk.Kriging(y, X, kernel, lk.RegressionModel.Constant, False, "BFGS", "LL", lk.KrigingParameters())
     x = 0.3 * np.ones(np.shape(X)[1])
 
-    loo, loograd = r.leaveOneOut(x, True)
+    loo, loograd = r.leaveOneOutFun(x, True)
     assert relative_error(loo, loo_ref) < tolerance
     assert relative_error(loograd, loograd_ref) < tolerance
 
-    ll, llgrad, llhess = r.logLikelihood(x, True, False)
+    ll, llgrad, llhess = r.logLikelihoodFun(x, True, False)
     assert relative_error(ll, ll_ref) < tolerance
     assert relative_error(llgrad, llgrad_ref) < tolerance
 
@@ -90,6 +90,6 @@ def test_data2(i):
 
     x = 0.3 * np.ones(np.shape(X)[1])
 
-    ll, llgrad, llhess = r.logLikelihood(x, True, False)
+    ll, llgrad, llhess = r.logLikelihoodFun(x, True, False)
     assert relative_error(ll, ll_ref) < tolerance
     assert relative_error(llgrad, llgrad_ref) < tolerance
diff --git a/bindings/R/rlibkriging/tests/testthat/bench-KrigingFit.R b/bindings/R/rlibkriging/tests/testthat/bench-KrigingFit.R
index 6923994a..31eb1c45 100644
--- a/bindings/R/rlibkriging/tests/testthat/bench-KrigingFit.R
+++ b/bindings/R/rlibkriging/tests/testthat/bench-KrigingFit.R
@@ -27,15 +27,15 @@ for (i in 1:N) {
                     )
 
 
-  ll_cpp <- logLikelihood(r, t(as.list(r)$theta))$logLikelihood
+  ll_cpp <- logLikelihoodFun(r, t(as.list(r)$theta))$logLikelihood
   e <- new.env()
   ll_R <- DiceKriging::logLikFun(k@covariance@range.val, k, e)
   
-  gll_cpp <- logLikelihood(r, t(as.list(r)$theta))$logLikelihoodGrad
+  gll_cpp <- logLikelihoodFun(r, t(as.list(r)$theta))$logLikelihoodGrad
   gll_R <- DiceKriging::logLikGrad(k@covariance@range.val, k, e)
     
   if (abs(ll_cpp - DiceKriging::logLikFun(param=as.numeric(t(as.list(r)$theta)),model=k))/ll_cpp>.1)
-    stop("LL function is not the same bw DiceKriging/libKriging: ",logLikelihood(r,t(as.list(r)$theta))," vs. ",DiceKriging::logLikFun(param=as.numeric(t(as.list(r)$theta)),model=k))  
+    stop("LL function is not the same bw DiceKriging/libKriging: ",logLikelihoodFun(r,t(as.list(r)$theta))," vs. ",DiceKriging::logLikFun(param=as.numeric(t(as.list(r)$theta)),model=k))  
   
   if ((ll_cpp - ll_R)/ll_R < -.01 )
     warning("libKriging LL ",ll_cpp," << DiceKriging LL ",ll_R)
diff --git a/bindings/R/rlibkriging/tests/testthat/bench-KrigingLogLik.R b/bindings/R/rlibkriging/tests/testthat/bench-KrigingLogLik.R
index ec6f919f..725d08f3 100644
--- a/bindings/R/rlibkriging/tests/testthat/bench-KrigingLogLik.R
+++ b/bindings/R/rlibkriging/tests/testthat/bench-KrigingLogLik.R
@@ -24,13 +24,13 @@ for (x in xx){
 r <- Kriging(y, X, "gauss")
 i <- 1
 for (x in xx){
-  times$cpp_ll[i]=system.time(for (j in 1:times.n) ll2x <- logLikelihood(r,rep(x,3)))
-  times$cpp_gll[i]=system.time(for (j in 1:times.n) gll2x <- logLikelihood(r,rep(x,3),grad=T))
+  times$cpp_ll[i]=system.time(for (j in 1:times.n) ll2x <- logLikelihoodFun(r,rep(x,3)))
+  times$cpp_gll[i]=system.time(for (j in 1:times.n) gll2x <- logLikelihoodFun(r,rep(x,3),grad=T))
   i <- i+1
 }
 
 plot(xx,Vectorize(function(x)DiceKriging::logLikFun(rep(x,3),k))(xx))
-points(xx,Vectorize(function(x)logLikelihood(r,rep(x,3))$logLikelihood)(xx),col='red')
+points(xx,Vectorize(function(x)logLikelihoodFun(r,rep(x,3))$logLikelihood)(xx),col='red')
 
 plot(xx,log(times$R_ll),ylim=c(log(min(min(times$R_ll,na.rm = T),min(times$cpp_ll,na.rm = T))),log(max(max(times$R_ll,na.rm = T),max(times$cpp_ll,na.rm = T)))),xlab="x",ylab="log(user_time (s))", panel.first=grid())
 text(20,-1,"DiceKriging")
diff --git a/bindings/R/rlibkriging/tests/testthat/bench-KrigingLogMargPost.R b/bindings/R/rlibkriging/tests/testthat/bench-KrigingLogMargPost.R
index 8a6745b8..b3e7fb0f 100644
--- a/bindings/R/rlibkriging/tests/testthat/bench-KrigingLogMargPost.R
+++ b/bindings/R/rlibkriging/tests/testthat/bench-KrigingLogMargPost.R
@@ -47,8 +47,8 @@ for (x in xx){
 r <- Kriging(y, X, "gauss")#"matern3_2")
 i <- 1
 for (x in xx){
-  times$cpp_ll[i]=system.time(for (j in 1:times.n) ll2x <- logMargPost(r,rep(x,3))$logMargPost)
-  times$cpp_gll[i]=system.time(for (j in 1:times.n) gll2x <- logMargPost(r,rep(x,3),grad=T)$logMargPostGrad)
+  times$cpp_ll[i]=system.time(for (j in 1:times.n) ll2x <- logMargPostFun(r,rep(x,3))$logMargPost)
+  times$cpp_gll[i]=system.time(for (j in 1:times.n) gll2x <- logMargPostFun(r,rep(x,3),grad=T)$logMargPostGrad)
   i <- i+1
 }
 
diff --git a/bindings/R/rlibkriging/tests/testthat/notest-RobustGaSP.R b/bindings/R/rlibkriging/tests/testthat/notest-RobustGaSP.R
index 807759b5..29f87d6c 100644
--- a/bindings/R/rlibkriging/tests/testthat/notest-RobustGaSP.R
+++ b/bindings/R/rlibkriging/tests/testthat/notest-RobustGaSP.R
@@ -25,7 +25,7 @@ lmp_rgasp = function(X, model=k) {if (!is.matrix(X)) X = matrix(X,ncol=1);
                   # print(dim(X));
                   apply(X,1,
                     function(x) {
-                      #y=-logMargPost(r,matrix(unlist(x),ncol=2))$logMargPost
+                      #y=-logMargPostFun(r,matrix(unlist(x),ncol=2))$logMargPost
                       y=RobustGaSP:::neg_log_marginal_post_approx_ref(param=exp(x),nugget=0, nugget.est=model@nugget.est, 
                                         R0=model@R0,X=model@X, zero_mean=model@zero_mean,output=model@output, 
                                         CL=model@CL, 
@@ -41,7 +41,7 @@ lmp_lk = function(X) {if (!is.matrix(X)) X = matrix(X,ncol=1);
                   # print(dim(X));
                   apply(X,1,
                     function(x) {
-                      y=-logMargPost(r,matrix(unlist(exp(-exp(x))),ncol=1))$logMargPost
+                      y=-logMargPostFun(r,matrix(unlist(exp(-exp(x))),ncol=1))$logMargPost
                       y})}
 lines(seq(0.1,0.9,,101),lmp_lk(seq(0.1,0.9,,101)),col='red')
 
@@ -145,7 +145,7 @@ lmp_rgasp = function(X) {if (!is.matrix(X)) X = matrix(X,ncol=2);
                   # print(dim(X));
                   apply(X,1,
                     function(x) {
-                      #y=-logMargPost(r,matrix(unlist(x),ncol=2))$logMargPost
+                      #y=-logMargPostFun(r,matrix(unlist(x),ncol=2))$logMargPost
                       y=RobustGaSP:::neg_log_marginal_post_approx_ref(param=exp(x),nugget=0, nugget.est=model@nugget.est, 
                                         R0=model@R0,X=model@X, zero_mean=model@zero_mean,output=model@output, 
                                         CL=model@CL, 
@@ -159,7 +159,7 @@ lmp_lk = function(X,deriv=FALSE) {if (!is.matrix(X)) X = matrix(X,ncol=2);
                   # print(dim(X));
                   apply(X,1,
                     function(x) {
-                      y=-logMargPost(r,matrix(unlist(exp(-exp(x))),ncol=2),deriv)$logMargPost
+                      y=-logMargPostFun(r,matrix(unlist(exp(-exp(x))),ncol=2),deriv)$logMargPost
                       y})}
 
 
diff --git a/bindings/R/rlibkriging/tests/testthat/test-KrigingFit.R b/bindings/R/rlibkriging/tests/testthat/test-KrigingFit.R
index 017af0cb..b55e8a53 100644
--- a/bindings/R/rlibkriging/tests/testthat/test-KrigingFit.R
+++ b/bindings/R/rlibkriging/tests/testthat/test-KrigingFit.R
@@ -13,12 +13,12 @@ r = NULL
 k = DiceKriging::km(design=X,response=y,covtype = "gauss",control = list(trace=F))
 r <- Kriging(y, X, "gauss")
 
-ll = Vectorize(function(x) logLikelihood(r,x)$logLikelihood)
+ll = Vectorize(function(x) logLikelihoodFun(r,x)$logLikelihood)
 plot(ll,xlim=c(0.000001,1))
   for (x in seq(0.000001,1,,11)){
     envx = new.env()
-    ll2x = logLikelihood(r,x)$logLikelihood
-    gll2x = logLikelihood(r,x,grad = T)$logLikelihoodGrad
+    ll2x = logLikelihoodFun(r,x)$logLikelihood
+    gll2x = logLikelihoodFun(r,x,grad = T)$logLikelihoodGrad
     arrows(x,ll2x,x+.1,ll2x+.1*gll2x,col='red')
   }
 
@@ -53,7 +53,7 @@ ll = function(X) {if (!is.matrix(X)) X = matrix(X,ncol=2);
                     function(x) {
                       # print(dim(x))
                       #print(matrix(unlist(x),ncol=2));
-                      y=-logLikelihood(r,matrix(unlist(x),ncol=2))$logLikelihood
+                      y=-logLikelihoodFun(r,matrix(unlist(x),ncol=2))$logLikelihood
                       #print(y);
                       y})}
 #DiceView::contourview(ll,xlim=c(0.01,2),ylim=c(0.01,2))
@@ -111,7 +111,7 @@ apply(X,1,
       function(x) {
         # print(dim(x))
         #print(matrix(unlist(x),ncol=2));
-        y=-logLikelihood(r,matrix(unlist(x),ncol=2))$logLikelihood
+        y=-logLikelihoodFun(r,matrix(unlist(x),ncol=2))$logLikelihood
         #print(y);
         y})}
 #DiceView::contourview(ll,xlim=c(0.01,2),ylim=c(0.01,2))
@@ -198,7 +198,7 @@ apply(X,1,
       function(x) {
         # print(dim(x))
         #print(matrix(unlist(x),ncol=2));
-        -logLikelihood(r,matrix(unlist(x),ncol=2))$logLikelihood
+        -logLikelihoodFun(r,matrix(unlist(x),ncol=2))$logLikelihood
         #print(y);
         })}
 #DiceView::contourview(ll,xlim=c(0.01,2),ylim=c(0.01,2))
diff --git a/bindings/R/rlibkriging/tests/testthat/test-KrigingLeaveOneOut.R b/bindings/R/rlibkriging/tests/testthat/test-KrigingLeaveOneOut.R
index 75351373..ed2bac00 100644
--- a/bindings/R/rlibkriging/tests/testthat/test-KrigingLeaveOneOut.R
+++ b/bindings/R/rlibkriging/tests/testthat/test-KrigingLeaveOneOut.R
@@ -23,12 +23,12 @@ for (x in seq(0.01,1,,11)){
 }
 
 r <- Kriging(y, X, kernel)
-ll2 = function(theta) leaveOneOut(r,theta)$leaveOneOut[1]
+ll2 = function(theta) leaveOneOutFun(r,theta)$leaveOneOut[1]
 # plot(Vectorize(ll2),col='red',add=T) # FIXME fails with "error: chol(): decomposition failed"
 for (x in seq(0.01,1,,11)){
   envx = new.env()
-  ll2x = leaveOneOut(r,x)$leaveOneOut[1]
-  gll2x = leaveOneOut(r,x,grad=T)$leaveOneOutGrad
+  ll2x = leaveOneOutFun(r,x)$leaveOneOut[1]
+  gll2x = leaveOneOutFun(r,x,grad=T)$leaveOneOutGrad
   arrows(x,ll2x,x+.1,ll2x+.1*gll2x,col='red')
 }
 
@@ -36,10 +36,10 @@ precision <- 1e-8  # the following tests should work with it, since the computat
 x=.5
 xenv=new.env()
 test_that(desc="leaveOneOut is the same that DiceKriging one",
-         expect_equal(leaveOneOut(r,x)$leaveOneOut[1],DiceKriging::leaveOneOutFun(x,k,xenv),tolerance = precision))
+         expect_equal(leaveOneOutFun(r,x)$leaveOneOut[1],DiceKriging::leaveOneOutFun(x,k,xenv),tolerance = precision))
 
 test_that(desc="leaveOneOut Grad is the same that DiceKriging one",
-          expect_equal(leaveOneOut(r,x,grad=T)$leaveOneOutGrad,DiceKriging::leaveOneOutGrad(x,k,xenv),tolerance= precision))
+          expect_equal(leaveOneOutFun(r,x,grad=T)$leaveOneOutGrad,DiceKriging::leaveOneOutGrad(x,k,xenv),tolerance= precision))
 }
 
 
diff --git a/bindings/R/rlibkriging/tests/testthat/test-KrigingLeaveOneOut_3d.R b/bindings/R/rlibkriging/tests/testthat/test-KrigingLeaveOneOut_3d.R
index df5d6a5f..affe2b44 100644
--- a/bindings/R/rlibkriging/tests/testthat/test-KrigingLeaveOneOut_3d.R
+++ b/bindings/R/rlibkriging/tests/testthat/test-KrigingLeaveOneOut_3d.R
@@ -15,16 +15,16 @@ k = DiceKriging::km(design=X,response=y,covtype = kernel,control = list(trace=F)
 ll = function(theta) DiceKriging::leaveOneOutFun(theta,k)
 
 r <- Kriging(y, X, kernel)
-ll2 = function(theta) leaveOneOut(r,theta)
+ll2 = function(theta) leaveOneOutFun(r,theta)
 
 precision <- 1e-8  # the following tests should work with it, since the computations are analytical
 x=runif(d)
 xenv=new.env()
 test_that(desc="leaveOneOut is the same that DiceKriging one",
-         expect_equal(leaveOneOut(r,x)$leaveOneOut[1],DiceKriging::leaveOneOutFun(x,k,xenv),tolerance = precision))
+         expect_equal(leaveOneOutFun(r,x)$leaveOneOut[1],DiceKriging::leaveOneOutFun(x,k,xenv),tolerance = precision))
 
 test_that(desc="leaveOneOut Grad is the same that DiceKriging one",
-          expect_equal(t(leaveOneOut(r,x,grad=T)$leaveOneOutGrad),DiceKriging::leaveOneOutGrad(x,k,xenv),tolerance= precision))
+          expect_equal(t(leaveOneOutFun(r,x,grad=T)$leaveOneOutGrad),DiceKriging::leaveOneOutGrad(x,k,xenv),tolerance= precision))
 }
 
 
diff --git a/bindings/R/rlibkriging/tests/testthat/test-KrigingLogLik.R b/bindings/R/rlibkriging/tests/testthat/test-KrigingLogLik.R
index f47f1e37..1fb481c2 100644
--- a/bindings/R/rlibkriging/tests/testthat/test-KrigingLogLik.R
+++ b/bindings/R/rlibkriging/tests/testthat/test-KrigingLogLik.R
@@ -25,12 +25,12 @@ for (kernel in c("exp","matern3_2","matern5_2","gauss")) {
   
   library(rlibkriging)
   r <- Kriging(y, X, kernel)
-  ll2 = function(theta) logLikelihood(r,theta)$logLikelihood
+  ll2 = function(theta) logLikelihoodFun(r,theta)$logLikelihood
   # plot(Vectorize(ll2),col='red',add=T) # FIXME fails with "error: chol(): decomposition failed"
   for (x in seq(0.01,1,,11)){
     envx = new.env()
-    ll2x = logLikelihood(r,x)$logLikelihood
-    gll2x = logLikelihood(r,x,grad = T)$logLikelihoodGrad
+    ll2x = logLikelihoodFun(r,x)$logLikelihood
+    gll2x = logLikelihoodFun(r,x,grad = T)$logLikelihoodGrad
     arrows(x,ll2x,x+.1,ll2x+.1*gll2x,col='red')
   }
   
@@ -38,10 +38,10 @@ for (kernel in c("exp","matern3_2","matern5_2","gauss")) {
   x=.5
   xenv=new.env()
   test_that(desc="logLik is the same that DiceKriging one", 
-            expect_equal(logLikelihood(r,x)$logLikelihood[1],DiceKriging::logLikFun(x,k,xenv),tolerance = precision))
+            expect_equal(logLikelihoodFun(r,x)$logLikelihood[1],DiceKriging::logLikFun(x,k,xenv),tolerance = precision))
   
   test_that(desc="logLik Grad is the same that DiceKriging one", 
-            expect_equal(logLikelihood(r,x,grad=T)$logLikelihoodGrad,DiceKriging::logLikGrad(x,k,xenv),tolerance= precision))
+            expect_equal(logLikelihoodFun(r,x,grad=T)$logLikelihoodGrad,DiceKriging::logLikGrad(x,k,xenv),tolerance= precision))
 }
 
 
@@ -67,10 +67,10 @@ for (kernel in c("matern3_2","matern5_2","gauss","exp")) {
   x=c(.2,.5,.7)
   xenv=new.env()
   test_that(desc="logLik is the same that DiceKriging one", 
-            expect_equal(logLikelihood(r,x)$logLikelihood[1],DiceKriging::logLikFun(x,k,xenv),tolerance = precision))
+            expect_equal(logLikelihoodFun(r,x)$logLikelihood[1],DiceKriging::logLikFun(x,k,xenv),tolerance = precision))
   
   test_that(desc="logLik Grad is the same that DiceKriging one", 
-            expect_equal(logLikelihood(r,x,grad=T)$logLikelihoodGrad[1,],t(DiceKriging::logLikGrad(x,k,xenv))[1,],tolerance= precision))
+            expect_equal(logLikelihoodFun(r,x,grad=T)$logLikelihoodGrad[1,],t(DiceKriging::logLikGrad(x,k,xenv))[1,],tolerance= precision))
 }
 
 
diff --git a/bindings/R/rlibkriging/tests/testthat/test-KrigingLogLikGradHess.R b/bindings/R/rlibkriging/tests/testthat/test-KrigingLogLikGradHess.R
index b8a556d4..4030c4df 100644
--- a/bindings/R/rlibkriging/tests/testthat/test-KrigingLogLikGradHess.R
+++ b/bindings/R/rlibkriging/tests/testthat/test-KrigingLogLikGradHess.R
@@ -37,9 +37,9 @@ hll = function(theta,eps=0.0001) {
 
 library(rlibkriging)
 r <- Kriging(y, X, kernel)
-ll_C = function(theta) logLikelihood(r,theta)$logLikelihood[1]
-gll_C = function(theta) t(logLikelihood(r,theta,grad=T)$logLikelihoodGrad)
-hll_C = function(theta) logLikelihood(r,theta,hess=T)$logLikelihoodHess[,,]
+ll_C = function(theta) logLikelihoodFun(r,theta)$logLikelihood[1]
+gll_C = function(theta) t(logLikelihoodFun(r,theta,grad=T)$logLikelihoodGrad)
+hll_C = function(theta) logLikelihoodFun(r,theta,hess=T)$logLikelihoodHess[,,]
 
 
 x=runif(d)
diff --git a/bindings/R/rlibkriging/tests/testthat/test-KrigingLogMargPost.R b/bindings/R/rlibkriging/tests/testthat/test-KrigingLogMargPost.R
index f63c4fc6..f9481c81 100644
--- a/bindings/R/rlibkriging/tests/testthat/test-KrigingLogMargPost.R
+++ b/bindings/R/rlibkriging/tests/testthat/test-KrigingLogMargPost.R
@@ -56,22 +56,22 @@ for (kernel in c("matern5_2","matern3_2")) {
   r <- Kriging(y, X, kernel)
   ## Should be equal:
   #lmp(1.0); lmp_deriv(1.0);
-  #logMargPost(r,1.0,grad = T)
+  #logMargPostFun(r,1.0,grad = T)
   #lmp(0.1); lmp_deriv(0.1);
-  #logMargPost(r,0.1,grad = T)
-  #ll2 = function(theta) logMargPost(r,theta)$logMargPost
+  #logMargPostFun(r,0.1,grad = T)
+  #ll2 = function(theta) logMargPostFun(r,theta)$logMargPost
   # plot(Vectorize(ll2),col='red',add=T,xlim=c(0.01,2)) # FIXME fails with "error: chol(): decomposition failed"
   for (x in seq(0.01,2,,11)){
-    ll2x = logMargPost(r,x)$logMargPost
-    gll2x = logMargPost(r,x,grad = T)$logMargPostGrad
+    ll2x = logMargPostFun(r,x)$logMargPost
+    gll2x = logMargPostFun(r,x,grad = T)$logMargPostGrad
     arrows(x,ll2x,x+.1,ll2x+.1*gll2x,col='red')
   }
   
   precision <- 1e-8  # the following tests should work with it, since the computations are analytical
   x=.5
   test_that(desc="logMargPost is the same that RobustGaSP one", 
-            expect_equal(logMargPost(r,x)$logMargPost[1],lmp(x),tolerance = precision))
+            expect_equal(logMargPostFun(r,x)$logMargPost[1],lmp(x),tolerance = precision))
   
   test_that(desc="logMargPost Grad is the same that RobustGaSP one", 
-            expect_equal(logMargPost(r,x,grad = T)$logMargPostGrad[1],lmp_deriv(x),tolerance= precision))
+            expect_equal(logMargPostFun(r,x,grad = T)$logMargPostGrad[1],lmp_deriv(x),tolerance= precision))
 }
diff --git a/bindings/R/rlibkriging/tests/testthat/test-KrigingMethods.R b/bindings/R/rlibkriging/tests/testthat/test-KrigingMethods.R
index 219e140c..9655d034 100644
--- a/bindings/R/rlibkriging/tests/testthat/test-KrigingMethods.R
+++ b/bindings/R/rlibkriging/tests/testthat/test-KrigingMethods.R
@@ -20,11 +20,11 @@ points(X)
 r <- Kriging(y, X,"gauss",parameters = list(theta=matrix(runif(40),ncol=2)))
 
 # ll = function(X) {
-#   logLikelihood(r,X,grad=F)$logLikelihood
+#   logLikelihoodFun(r,X,grad=F)$logLikelihood
 # }
 # contour(x,x,matrix(ll(as.matrix(expand.grid(x,x))),nrow=length(x)),nlevels = 30)
 # gll = function(X) {
-#   logLikelihood(r,X,grad=T)$logLikelihoodGrad
+#   logLikelihoodFun(r,X,grad=T)$logLikelihoodGrad
 # }
 # for (ix in 1:21) {
 # for (iy in 1:21) {
@@ -41,44 +41,44 @@ p = capture.output(print(r))
 
 context("logLikelihood")
 
-ll = function(Theta){apply(Theta,1,function(theta) logLikelihood(r,theta)$logLikelihood)}
+ll = function(Theta){apply(Theta,1,function(theta) logLikelihoodFun(r,theta)$logLikelihood)}
 t=seq(0.01,2,,51)
 contour(t,t,matrix(ll(as.matrix(expand.grid(t,t))),nrow=length(t)),nlevels = 30)
 points(as.list(r)$theta[1],as.list(r)$theta[2],pch=20)
-# logLikelihood(r,t(as.list(r)$theta))
+# logLikelihoodFun(r,t(as.list(r)$theta))
 
-test_that("logLikelihood returned",
-          expect_equal(names(logLikelihood(r,runif(d))),c("logLikelihood")))
-test_that("logLikelihood logLikelihoodGrad returned",
-          expect_equal(names(logLikelihood(r,runif(d),grad=T)),c("logLikelihood","logLikelihoodGrad")))
-test_that("logLikelihood logLikelihoodGrad logLikelihoodHess returned",
-          expect_equal(names(logLikelihood(r,runif(d),grad=T,hess=T)),c("logLikelihood","logLikelihoodGrad","logLikelihoodHess")))
+test_that("logLikelihoodFun returned",
+          expect_equal(names(logLikelihoodFun(r,runif(d))),c("logLikelihood")))
+test_that("logLikelihoodFun logLikelihoodGrad returned",
+          expect_equal(names(logLikelihoodFun(r,runif(d),grad=T)),c("logLikelihood","logLikelihoodGrad")))
+test_that("logLikelihoodFun logLikelihoodGrad logLikelihoodHess returned",
+          expect_equal(names(logLikelihoodFun(r,runif(d),grad=T,hess=T)),c("logLikelihood","logLikelihoodGrad","logLikelihoodHess")))
 
-test_that("logLikelihood dim",
-          expect_equal(dim(logLikelihood(r,rbind(runif(d),runif(d)))$logLikelihood),c(2,1)))
+test_that("logLikelihoodFun dim",
+          expect_equal(dim(logLikelihoodFun(r,rbind(runif(d),runif(d)))$logLikelihood),c(2,1)))
 test_that("logLikelihoodGrad dim",
-          expect_equal(dim(logLikelihood(r,rbind(runif(d),runif(d)),grad=T)$logLikelihoodGrad),c(2,d)))
+          expect_equal(dim(logLikelihoodFun(r,rbind(runif(d),runif(d)),grad=T)$logLikelihoodGrad),c(2,d)))
 test_that("logLikelihoodHess dim",
-          expect_equal(dim(logLikelihood(r,rbind(runif(d),runif(d)),grad=T,hess=T)$logLikelihoodHess),c(2,d,d)))
+          expect_equal(dim(logLikelihoodFun(r,rbind(runif(d),runif(d)),grad=T,hess=T)$logLikelihoodHess),c(2,d,d)))
 
 
 context("leaveOneOut")
 
-loo = function(Theta){apply(Theta,1,function(theta) leaveOneOut(r,theta)$leaveOneOut)}
+loo = function(Theta){apply(Theta,1,function(theta) leaveOneOutFun(r,theta)$leaveOneOut)}
 t=seq(0.01,2,,51)
 contour(t,t,matrix(loo(as.matrix(expand.grid(t,t))),nrow=length(t)),nlevels = 30)
 points(as.list(r)$theta[1],as.list(r)$theta[2])
-# leaveOneOut(r,t(as.list(r)$theta))
+# leaveOneOutFun(r,t(as.list(r)$theta))
 
 test_that("leaveOneOut returned",
-          expect_equal(names(leaveOneOut(r,runif(d))),c("leaveOneOut")))
+          expect_equal(names(leaveOneOutFun(r,runif(d))),c("leaveOneOut")))
 test_that("leaveOneOut leaveOneOutGrad returned",
-          expect_equal(names(leaveOneOut(r,runif(d),grad=T)),c("leaveOneOut","leaveOneOutGrad")))
+          expect_equal(names(leaveOneOutFun(r,runif(d),grad=T)),c("leaveOneOut","leaveOneOutGrad")))
 
 test_that("leaveOneOut dim",
-          expect_equal(dim(leaveOneOut(r,rbind(runif(d),runif(d)))$leaveOneOut),c(2,1)))
+          expect_equal(dim(leaveOneOutFun(r,rbind(runif(d),runif(d)))$leaveOneOut),c(2,1)))
 test_that("leaveOneOutGrad dim",
-          expect_equal(dim(leaveOneOut(r,rbind(runif(d),runif(d)),grad=T)$leaveOneOutGrad),c(2,d)))
+          expect_equal(dim(leaveOneOutFun(r,rbind(runif(d),runif(d)),grad=T)$leaveOneOutGrad),c(2,d)))
 
 
 context("predict")
@@ -127,11 +127,11 @@ points(X2,col='red')
 
 #r20 <- Kriging(c(y,y2), rbind(X,X2),"gauss")
 #ll = function(X) {
-#  logLikelihood(r20,X,grad=F)$logLikelihood
+#  logLikelihoodFun(r20,X,grad=F)$logLikelihood
 #}
 #contour(x,x,matrix(ll(as.matrix(expand.grid(x,x))),nrow=length(x)),nlevels = 30)
 #gll = function(X) {
-#  logLikelihood(r20,X,grad=T)$logLikelihoodGrad
+#  logLikelihoodFun(r20,X,grad=T)$logLikelihoodGrad
 #}
 #for (ix in 1:21) {
 #for (iy in 1:21) {
@@ -148,14 +148,14 @@ contour(t,t,matrix(ll(as.matrix(expand.grid(t,t))),nrow=length(t)),nlevels = 30)
 points(as.list(r)$theta[1],as.list(r)$theta[2])
 
 r2 <- Kriging(c(y,y2), rbind(X,X2),"gauss", parameters = list(theta=matrix(as.list(r)$theta,ncol=2)))
-ll2 = function(Theta){apply(Theta,1,function(theta) logLikelihood(r2,theta)$logLikelihood)}
+ll2 = function(Theta){apply(Theta,1,function(theta) logLikelihoodFun(r2,theta)$logLikelihood)}
 contour(t,t,matrix(ll2(as.matrix(expand.grid(t,t))),nrow=length(t)),nlevels = 30,add=T,col='red')
 points(as.list(r2)$theta[1],as.list(r2)$theta[2],col='red')
 
 p2 = capture.output(print(r2))
 
 update(object=r,y2,X2)
-llu = function(Theta){apply(Theta,1,function(theta) logLikelihood(r,theta)$logLikelihood)}
+llu = function(Theta){apply(Theta,1,function(theta) logLikelihoodFun(r,theta)$logLikelihood)}
 contour(t,t,matrix(llu(as.matrix(expand.grid(t,t))),nrow=length(t)),nlevels = 30,add=T,col='blue')
 points(as.list(r)$theta[1],as.list(r)$theta[2],col='blue')
 
diff --git a/bindings/R/rlibkriging/tests/testthat/test-NuggetKrigingFit.R b/bindings/R/rlibkriging/tests/testthat/test-NuggetKrigingFit.R
index 9b08ab19..6843360d 100644
--- a/bindings/R/rlibkriging/tests/testthat/test-NuggetKrigingFit.R
+++ b/bindings/R/rlibkriging/tests/testthat/test-NuggetKrigingFit.R
@@ -22,9 +22,9 @@ alpha_r = as.list(r)$sigma2/(as.list(r)$sigma2+as.list(r)$nugget)
 test_that(desc="Nugget / Fit: 1D / fit of alpha by DiceKriging is same that libKriging",
           expect_equal(alpha_k,alpha_r, tol= 1e-4))
 
-ll = Vectorize(function(x) logLikelihood(r,c(x,alpha_k))$logLikelihood)
+ll = Vectorize(function(x) logLikelihoodFun(r,c(x,alpha_k))$logLikelihood)
 plot(ll,xlim=c(0.001,1))
-#ll = Vectorize(function(x) logLikelihood(r,c(x,alpha_r))$logLikelihood)
+#ll = Vectorize(function(x) logLikelihoodFun(r,c(x,alpha_r))$logLikelihood)
 #plot(ll_,xlim=c(0.001,1))
 
 theta_ref = optimize(ll,interval=c(0.001,1),maximum=T)$maximum
@@ -64,7 +64,7 @@ ll = function(X) {if (!is.matrix(X)) X = matrix(X,ncol=2);
                   # print(dim(X));
                   apply(X,1,
                     function(x) {
-                      y=-logLikelihood(r,c(unlist(x),alpha_k))$logLikelihood
+                      y=-logLikelihoodFun(r,c(unlist(x),alpha_k))$logLikelihood
                       #print(y);
                       y})}
 #DiceView::contourview(ll,xlim=c(0.01,2),ylim=c(0.01,2))
@@ -127,7 +127,7 @@ apply(X,1,
       function(x) {
         # print(dim(x))
         #print(matrix(unlist(x),ncol=2));
-        y=-logLikelihood(r,c(unlist(x),alpha_k))$logLikelihood
+        y=-logLikelihoodFun(r,c(unlist(x),alpha_k))$logLikelihood
         #print(y);
         y})}
 #DiceView::contourview(ll,xlim=c(0.01,2),ylim=c(0.01,2))
@@ -178,7 +178,7 @@ apply(X,1,
       function(x) {
         # print(dim(x))
         #print(matrix(unlist(x),ncol=2));
-        -logLikelihood(r,c(unlist(x),alpha_k))$logLikelihood
+        -logLikelihoodFun(r,c(unlist(x),alpha_k))$logLikelihood
         #print(y);
         })}
 #DiceView::contourview(ll,xlim=c(0.01,2),ylim=c(0.01,2))
diff --git a/bindings/R/rlibkriging/tests/testthat/test-NuggetKrigingLogLik.R b/bindings/R/rlibkriging/tests/testthat/test-NuggetKrigingLogLik.R
index 2718b691..e3f20a9c 100644
--- a/bindings/R/rlibkriging/tests/testthat/test-NuggetKrigingLogLik.R
+++ b/bindings/R/rlibkriging/tests/testthat/test-NuggetKrigingLogLik.R
@@ -27,13 +27,13 @@ for (kernel in c("exp","matern3_2","matern5_2","gauss")) {
   
   library(rlibkriging)
   r <- NuggetKriging(y, X, kernel, parameters=list(nugget=0,is_nugget_estim=TRUE))
-  ll2_theta = function(theta) logLikelihood(r,c(theta,alpha0))$logLikelihood
+  ll2_theta = function(theta) logLikelihoodFun(r,c(theta,alpha0))$logLikelihood
   # second arg is alpha=1 for nugget=0
   # plot(Vectorize(ll2),col='red'), add=T) 
   for (x in seq(0.01,1,,11)){
     envx = new.env()
-    ll2x = logLikelihood(r,c(x,alpha0))$logLikelihood
-    gll2x = logLikelihood(r,c(x,alpha0),grad = T)$logLikelihoodGrad[,1]
+    ll2x = logLikelihoodFun(r,c(x,alpha0))$logLikelihood
+    gll2x = logLikelihoodFun(r,c(x,alpha0),grad = T)$logLikelihoodGrad[,1]
     arrows(x,ll2x,x+.1,ll2x+.1*gll2x,col='red')
   }
   
@@ -46,12 +46,12 @@ for (kernel in c("exp","matern3_2","matern5_2","gauss")) {
     gllx = DiceKriging::logLikGrad(c(theta0,x),k,envx)[2,]
     arrows(x,llx,x+.1,llx+.1*gllx)
   }
-  ll2_alpha = function(alpha) logLikelihood(r,c(theta0,alpha))$logLikelihood
+  ll2_alpha = function(alpha) logLikelihoodFun(r,c(theta0,alpha))$logLikelihood
   #plot(Vectorize(ll2_alpha),col='red',add=T)
   for (x in seq(0.01,1,,11)){
     envx = new.env()
-    ll2x = logLikelihood(r,c(theta0,x))$logLikelihood
-    gll2x = logLikelihood(r,c(theta0,x),grad = T)$logLikelihoodGrad[,2]
+    ll2x = logLikelihoodFun(r,c(theta0,x))$logLikelihood
+    gll2x = logLikelihoodFun(r,c(theta0,x),grad = T)$logLikelihoodGrad[,2]
     arrows(x,ll2x,x+.1,ll2x+.1*gll2x,col='red')
   }
 
@@ -60,28 +60,28 @@ for (kernel in c("exp","matern3_2","matern5_2","gauss")) {
   xenv=new.env()
   test_that(desc="logLik is the same that DiceKriging one", 
             expect_equal(
-              logLikelihood(r,c(theta0,x))$logLikelihood[1]
+              logLikelihoodFun(r,c(theta0,x))$logLikelihood[1]
               ,
               DiceKriging::logLikFun(c(theta0,x),k,xenv)
               ,tolerance = precision))
   
   test_that(desc="logLik Grad is the same that DiceKriging one", 
             expect_equal(
-              logLikelihood(r,c(theta0,x),grad=T)$logLikelihoodGrad[,2]
+              logLikelihoodFun(r,c(theta0,x),grad=T)$logLikelihoodGrad[,2]
               ,
               DiceKriging::logLikGrad(c(theta0,x),k,xenv)[2,]
               ,tolerance= precision))
   xenv=new.env()
   test_that(desc="logLik is the same that DiceKriging one", 
             expect_equal(
-              logLikelihood(r,c(x,alpha0))$logLikelihood[1]
+              logLikelihoodFun(r,c(x,alpha0))$logLikelihood[1]
               ,
               DiceKriging::logLikFun(c(x,alpha0),k,xenv)
               ,tolerance = precision))
   
   test_that(desc="logLik Grad is the same that DiceKriging one", 
             expect_equal(
-              logLikelihood(r,c(x,alpha0),grad=T)$logLikelihoodGrad[,1]
+              logLikelihoodFun(r,c(x,alpha0),grad=T)$logLikelihoodGrad[,1]
               ,
               DiceKriging::logLikGrad(c(x,alpha0),k,xenv)[1,]
               ,tolerance= precision))
@@ -112,14 +112,14 @@ for (kernel in c("matern3_2","matern5_2","gauss","exp")) {
   xenv=new.env()
   test_that(desc="logLik is the same that DiceKriging one", 
             expect_equal(
-              logLikelihood(r,x)$logLikelihood[1]
+              logLikelihoodFun(r,x)$logLikelihood[1]
               ,
               DiceKriging::logLikFun(x,k,xenv)
               ,tolerance = precision))
   
   test_that(desc="logLik Grad is the same that DiceKriging one", 
             expect_equal(
-              logLikelihood(r,x,grad=T)$logLikelihoodGrad[1,]
+              logLikelihoodFun(r,x,grad=T)$logLikelihoodGrad[1,]
               ,
               t(DiceKriging::logLikGrad(x,k,xenv))[1,]
               ,tolerance= precision))
diff --git a/bindings/R/rlibkriging/tests/testthat/test-NuggetKrigingMethods.R b/bindings/R/rlibkriging/tests/testthat/test-NuggetKrigingMethods.R
index 4b630235..f287cd53 100644
--- a/bindings/R/rlibkriging/tests/testthat/test-NuggetKrigingMethods.R
+++ b/bindings/R/rlibkriging/tests/testthat/test-NuggetKrigingMethods.R
@@ -20,11 +20,11 @@ points(X)
 r <- NuggetKriging(y, X,"gauss",parameters = list(theta=matrix(runif(40),ncol=2),nugget=0,sigma2=1))
 l= as.list(r)
 # ll = function(X) {
-#   logLikelihood(r,X,grad=F)$logLikelihood
+#   logLikelihoodFun(r,X,grad=F)$logLikelihood
 # }
 # contour(x,x,matrix(ll(as.matrix(expand.grid(x,x))),nrow=length(x)),nlevels = 30)
 # gll = function(X) {
-#   logLikelihood(r,X,grad=T)$logLikelihoodGrad
+#   logLikelihoodFun(r,X,grad=T)$logLikelihoodGrad
 # }
 # for (ix in 1:21) {
 # for (iy in 1:21) {
@@ -41,21 +41,21 @@ p = capture.output(print(r))
 
 context("logLikelihood")
 
-ll = function(Theta){apply(Theta,1,function(theta) logLikelihood(r,c(theta,l$sigma2/(l$nugget+l$sigma2)))$logLikelihood)}
+ll = function(Theta){apply(Theta,1,function(theta) logLikelihoodFun(r,c(theta,l$sigma2/(l$nugget+l$sigma2)))$logLikelihood)}
 t=seq(0.01,2,,51)
 contour(t,t,matrix(ll(as.matrix(expand.grid(t,t))),nrow=length(t)),nlevels = 30)
 points(as.list(r)$theta[1],as.list(r)$theta[2])
-# logLikelihood(r,t(as.list(r)$theta))
+# logLikelihoodFun(r,t(as.list(r)$theta))
 
-test_that("logLikelihood returned",
-          expect_equal(names(logLikelihood(r,runif(d+1))),c("logLikelihood")))
-test_that("logLikelihood logLikelihoodGrad returned",
-          expect_equal(names(logLikelihood(r,runif(d+1),grad=T)),c("logLikelihood","logLikelihoodGrad")))
+test_that("logLikelihoodFun returned",
+          expect_equal(names(logLikelihoodFun(r,runif(d+1))),c("logLikelihood")))
+test_that("logLikelihoodFun logLikelihoodGrad returned",
+          expect_equal(names(logLikelihoodFun(r,runif(d+1),grad=T)),c("logLikelihood","logLikelihoodGrad")))
 
-test_that("logLikelihood dim",
-          expect_equal(dim(logLikelihood(r,rbind(runif(d+1),runif(d+1)))$logLikelihood),c(2,1)))
+test_that("logLikelihoodFun dim",
+          expect_equal(dim(logLikelihoodFun(r,rbind(runif(d+1),runif(d+1)))$logLikelihood),c(2,1)))
 test_that("logLikelihoodGrad dim",
-          expect_equal(dim(logLikelihood(r,rbind(runif(d+1),runif(d+1)),grad=T)$logLikelihoodGrad),c(2,d+1)))
+          expect_equal(dim(logLikelihoodFun(r,rbind(runif(d+1),runif(d+1)),grad=T)$logLikelihoodGrad),c(2,d+1)))
 
 
 context("predict")
@@ -104,11 +104,11 @@ points(X2,col='red')
 
 #r20 <- Kriging(c(y,y2), rbind(X,X2),"gauss")
 #ll = function(X) {
-#  logLikelihood(r20,X,grad=F)$logLikelihood
+#  logLikelihoodFun(r20,X,grad=F)$logLikelihood
 #}
 #contour(x,x,matrix(ll(as.matrix(expand.grid(x,x))),nrow=length(x)),nlevels = 30)
 #gll = function(X) {
-#  logLikelihood(r20,X,grad=T)$logLikelihoodGrad
+#  logLikelihoodFun(r20,X,grad=T)$logLikelihoodGrad
 #}
 #for (ix in 1:21) {
 #for (iy in 1:21) {
@@ -126,14 +126,14 @@ points(as.list(r)$theta[1],as.list(r)$theta[2],pch=20)
 
 set.seed(1234)
 r2 <- NuggetKriging(c(y,y2), rbind(X,X2),"gauss", parameters = list(theta=matrix(as.list(r)$theta,ncol=2,nrow=40,byrow = T)+matrix(0.01*runif(80),ncol=2)))
-ll2 = function(Theta){apply(Theta,1,function(theta) logLikelihood(r2,c(theta,as.list(r2)$sigma2/(as.list(r2)$sigma2+as.list(r2)$nugget)))$logLikelihood)}
+ll2 = function(Theta){apply(Theta,1,function(theta) logLikelihoodFun(r2,c(theta,as.list(r2)$sigma2/(as.list(r2)$sigma2+as.list(r2)$nugget)))$logLikelihood)}
 contour(t,t,matrix(ll2(as.matrix(expand.grid(t,t))),nrow=length(t)),nlevels = 30,add=T,col='red')
 points(as.list(r2)$theta[1],as.list(r2)$theta[2],col='red',pch=20)
 
 p2 = capture.output(print(r2))
 
 update(object=r,y2,X2)
-llu = function(Theta){apply(Theta,1,function(theta) logLikelihood(r,  c(theta,as.list(r2)$sigma2/(as.list(r2)$sigma2+as.list(r2)$nugget)) )$logLikelihood)}
+llu = function(Theta){apply(Theta,1,function(theta) logLikelihoodFun(r,  c(theta,as.list(r2)$sigma2/(as.list(r2)$sigma2+as.list(r2)$nugget)) )$logLikelihood)}
 contour(t,t,matrix(llu(as.matrix(expand.grid(t,t))),nrow=length(t)),nlevels = 30,add=T,col='blue')
 points(as.list(r)$theta[1],as.list(r)$theta[2],col='blue',pch=20)
 
diff --git a/bindings/R/rlibkriging/tests/testthat/test-RobustGaSP.R b/bindings/R/rlibkriging/tests/testthat/test-RobustGaSP.R
index 18949c8f..defc8e38 100644
--- a/bindings/R/rlibkriging/tests/testthat/test-RobustGaSP.R
+++ b/bindings/R/rlibkriging/tests/testthat/test-RobustGaSP.R
@@ -26,7 +26,7 @@ lmp_rgasp = function(X, model=k) {if (!is.matrix(X)) X = matrix(X,ncol=1);
                   # print(dim(X));
                   apply(X,1,
                     function(x) {
-                      #y=-logMargPost(r,matrix(unlist(x),ncol=2))$logMargPost
+                      #y=-logMargPostFun(r,matrix(unlist(x),ncol=2))$logMargPost
                       y=RobustGaSP:::neg_log_marginal_post_approx_ref(param=(x),nugget=0, nugget.est=model@nugget.est, 
                                         R0=model@R0,X=model@X, zero_mean=model@zero_mean,output=model@output, 
                                         CL=model@CL, 
@@ -44,7 +44,7 @@ lmp_lk = function(X) {if (!is.matrix(X)) X = matrix(X,ncol=1);
                   # print(dim(X));
                   apply(X,1,
                     function(x) {
-                      y=-logMargPost(r,matrix(unlist(exp(-(x))),ncol=1))$logMargPost
+                      y=-logMargPostFun(r,matrix(unlist(exp(-(x))),ncol=1))$logMargPost
                       y})}
 lmp_lk(1)
 
@@ -75,7 +75,7 @@ dlmp_rgasp = function(X, model=k) {if (!is.matrix(X)) X = matrix(X,ncol=1);
 #                                        b=1/(length(model@output))^{1/dim(as.matrix(model@input))[2]}*(0.2+dim(as.matrix(model@output))[2])))
 
 
-                      #y=-logMargPost(r,matrix(unlist(x),ncol=2))$logMargPost
+                      #y=-logMargPostFun(r,matrix(unlist(x),ncol=2))$logMargPost
                       y=RobustGaSP:::neg_log_marginal_post_approx_ref_deriv(param=(x),nugget=0, nugget.est=model@nugget.est, 
                                         R0=model@R0,X=model@X, zero_mean=model@zero_mean,output=model@output, 
                                         CL=model@CL, 
@@ -89,7 +89,7 @@ dlmp_rgasp(1)
 dlmp_lk = function(X) {if (!is.matrix(X)) X = matrix(X,ncol=1);
                   apply(X,1,
                     function(x) {
-                      y=-logMargPost(r,matrix(unlist(exp(-(x))),ncol=1),TRUE)$logMargPostGrad
+                      y=-logMargPostFun(r,matrix(unlist(exp(-(x))),ncol=1),TRUE)$logMargPostGrad
                       y})}
 -exp(-1)*dlmp_lk(1)
 
@@ -197,7 +197,7 @@ lmp_rgasp = function(X, model=k) {if (!is.matrix(X)) X = matrix(X,ncol=2);
                   # print(dim(X));
                   apply(X,1,
                     function(x) {
-                      #y=-logMargPost(r,matrix(unlist(x),ncol=2))$logMargPost
+                      #y=-logMargPostFun(r,matrix(unlist(x),ncol=2))$logMargPost
     
 #    print(RobustGaSP:::log_marginal_lik(param=(x),nugget=0,nugget_est=model@nugget.est, 
 #                                        R0=model@R0,X=model@X, zero_mean=model@zero_mean,
@@ -226,7 +226,7 @@ points(log(k@beta_hat[1]),log(k@beta_hat[2]))
 lmp_lk = function(X) {if (!is.matrix(X)) X = matrix(X,ncol=2);
                   apply(X,1,
                     function(x) {
-                      y=-logMargPost(r,matrix(unlist(exp(-(x))),ncol=2))$logMargPost
+                      y=-logMargPostFun(r,matrix(unlist(exp(-(x))),ncol=2))$logMargPost
                       y})}
 lmp_lk(c(1,1))
 
@@ -244,7 +244,7 @@ dlmp_rgasp = function(X, model=k) {if (!is.matrix(X)) X = matrix(X,ncol=2);
                   # print(dim(X));
                   apply(X,1,
                     function(x) {
-                      #y=-logMargPost(r,matrix(unlist(x),ncol=2))$logMargPost
+                      #y=-logMargPostFun(r,matrix(unlist(x),ncol=2))$logMargPost
                       y=RobustGaSP:::neg_log_marginal_post_approx_ref_deriv(param=(x),nugget=0, nugget.est=model@nugget.est, 
                                         R0=model@R0,X=model@X, zero_mean=model@zero_mean,output=model@output, 
                                         CL=model@CL, 
@@ -258,7 +258,7 @@ dlmp_rgasp(c(1,1))
 dlmp_lk = function(X) {if (!is.matrix(X)) X = matrix(X,ncol=2);
                   apply(X,1,
                     function(x) {
-                      y=-logMargPost(r,matrix(unlist(exp(-(x))),ncol=2),TRUE)$logMargPostGrad
+                      y=-logMargPostFun(r,matrix(unlist(exp(-(x))),ncol=2),TRUE)$logMargPostGrad
                       y})}
 -exp(-1)*dlmp_lk(c(1,1))
 
diff --git a/bindings/R/rlibkriging/tests/testthat/test-asDiceKriging.R b/bindings/R/rlibkriging/tests/testthat/test-asDiceKriging.R
index 81854a79..9afa4e89 100644
--- a/bindings/R/rlibkriging/tests/testthat/test-asDiceKriging.R
+++ b/bindings/R/rlibkriging/tests/testthat/test-asDiceKriging.R
@@ -28,17 +28,17 @@ KM1 <- KM(design = X, response = y, covtype = "gauss",
           parinit = c(.15))
 
 test_that("m1.leaveOneOutFun == KM1.leaveOneOutFun",
-          expect_true(leaveOneOutFun(km1@covariance@range.val, km1) ==
-                      leaveOneOutFun(km1@covariance@range.val, KM1)))
+          expect_true( DiceKriging::leaveOneOutFun(km1@covariance@range.val, km1) ==
+                      DiceKriging::leaveOneOutFun(km1@covariance@range.val, KM1)))
 
 test_that("m1.argmax(loo) == KM1.argmax(loo)", 
           expect_equal(km1@covariance@range.val,
                        KM1@covariance@range.val,
                        tol = 0.1))
 
-plot(Vectorize(function(.t) leaveOneOutFun(param = as.numeric(.t), model = km1)))
+plot(Vectorize(function(.t) DiceKriging::leaveOneOutFun(param = as.numeric(.t), model = km1)))
 abline(v = km1@covariance@range.val)
-plot(Vectorize(function(.t) leaveOneOut(KM1@Kriging, as.numeric(.t))),
+plot(Vectorize(function(.t) leaveOneOutFun(KM1@Kriging, as.numeric(.t))),
      add = TRUE, col = 'red')
 abline(v = KM1@covariance@range.val, col = 'red')
 
@@ -89,7 +89,7 @@ ll <- function(Theta){
 }
 as_ll <- function(Theta){
     apply(Theta, 1,
-          function(theta) logLikelihood(KM1@Kriging, theta)$logLikelihood[1])
+          function(theta) logLikelihoodFun(KM1@Kriging, theta)$logLikelihood[1])
 }
 t <- seq(from = 0.01, to = 2, length.out = 51)
 ttg <- expand.grid(t, t)
@@ -178,16 +178,16 @@ suppressWarnings(KM2 <<- rlibkriging::KM(formula = formula,design = design,
 ## plot(Vectorize(function(.t) DiceKriging::logLikFun(.t,k)[1]),
 ##      xlim = c(0.000001, 1))
 ## plot(Vectorize(function(.t)
-##         rlibkriging::logLikelihood(as_k@Kriging, .t)$logLikelihood[1]),
+##         rlibkriging::logLikelihoodFun(as_k@Kriging, .t)$logLikelihood[1]),
 ##      xlim = c(0.000001, 1))
 
 x <- runif(ncol(X))
 test_that("DiceKriging::logLik == rlibkriging::logLikelihood",
           expect_equal(DiceKriging::logLikFun(x, km2)[1],
-                       rlibkriging::logLikelihood(KM2@Kriging,x)$logLikelihood[1]))
+                       rlibkriging::logLikelihoodFun(KM2@Kriging,x)$logLikelihood[1]))
 test_that("DiceKriging::leaveOneOut == rlibkriging::leaveOneOut",
           expect_equal(DiceKriging::leaveOneOutFun(x, km2)[1],
-                       rlibkriging::leaveOneOut(KM2@Kriging,x)$leaveOneOut[1]))
+                       rlibkriging::leaveOneOutFun(KM2@Kriging,x)$leaveOneOut[1]))
 
 x <- matrix(x,ncol=d)
 test_that("Consitency of 'DiceKriging' and 'rlibkriging' 'predict' methods",
@@ -266,10 +266,10 @@ test_args <-  function(formula, design, response ,covtype, estim.method ) {
     ##}
     ##abline(v=k@covariance@range.val)
     ##if (e=="MLE") {
-    ##  plot(Vectorize(function(t)rlibkriging::logLikelihood(as_k@Kriging,t)$logLikelihood[1]),
+    ##  plot(Vectorize(function(t)rlibkriging::logLikelihoodFun(as_k@Kriging,t)$logLikelihood[1]),
     ##       xlim = c(0.0001,2), add=T, col='red')
     ##} else {
-    ##  plot(Vectorize(function(t)rlibkriging::leaveOneOut(as_k@Kriging,t)$leaveOneOut[1]),
+    ##  plot(Vectorize(function(t)rlibkriging::leaveOneOutFun(as_k@Kriging,t)$leaveOneOut[1]),
     ##       xlim=c(0.0001,2),add=T,col='red')
     ##}
     ##abline(v=as_k@covariance@range.val,col='red')
@@ -277,10 +277,10 @@ test_args <-  function(formula, design, response ,covtype, estim.method ) {
   t <- runif(ncol(X))
     test_that("DiceKriging::logLikFun == rlibkriging::logLikelihood",
               expect_equal(DiceKriging::logLikFun(t, km2)[1],
-                           rlibkriging::logLikelihood(KM2@Kriging,t)$logLikelihood[1]))
+                           rlibkriging::logLikelihoodFun(KM2@Kriging,t)$logLikelihood[1]))
     test_that("DiceKriging::leaveOneOutFun == rlibkriging::leaveOneOut",
               expect_equal(DiceKriging::leaveOneOutFun(t, km2)[1],
-                           rlibkriging::leaveOneOut(KM2@Kriging, t)$leaveOneOut[1]))
+                           rlibkriging::leaveOneOutFun(KM2@Kriging, t)$leaveOneOut[1]))
     
     x <- matrix(runif(d),ncol=d)
     test_that("DiceKriging::predict == rlibkriging::predict",
diff --git a/bindings/R/rlibkriging/tests/testthat/test-binding-consistency.R b/bindings/R/rlibkriging/tests/testthat/test-binding-consistency.R
index 9322192d..7e91c8b8 100644
--- a/bindings/R/rlibkriging/tests/testthat/test-binding-consistency.R
+++ b/bindings/R/rlibkriging/tests/testthat/test-binding-consistency.R
@@ -60,8 +60,8 @@ y = as.matrix(read.table(file = filey, header = FALSE, sep = ',')) # a matrix, n
 kernel <- "gauss"
 r <- Kriging(y, X, kernel)
 x <- 0.3 # extracted from 'for (x in seq(0.01,1,,11))' 
-llo = rlibkriging::leaveOneOut(r, x, grad = TRUE)
-loglik = rlibkriging::logLikelihood(r, x, grad = TRUE)
+llo = rlibkriging::leaveOneOutFun(r, x, grad = TRUE)
+loglik = rlibkriging::logLikelihoodFun(r, x, grad = TRUE)
 
 for (n in names(llo)) {
   filer <- file.path(refpath, sprintf("%s-result-%s.csv", prefix, n))
@@ -105,7 +105,7 @@ for (i in 1:length(logn)) {
 
   kernel <- "gauss"
   r <- Kriging(y, X, kernel)
-  loglik = rlibkriging::logLikelihood(r, x, grad = TRUE) # TODO needs also Hessian ?
+  loglik = rlibkriging::logLikelihoodFun(r, x, grad = TRUE) # TODO needs also Hessian ?
   for (n in names(loglik)) {
     filer <- file.path(refpath, sprintf("%s-%i-result-%s.csv", prefix, i, n))
     if (do_write) {
diff --git a/tests/KrigingTest.cpp b/tests/KrigingTest.cpp
index 7e599a17..2aae2d28 100644
--- a/tests/KrigingTest.cpp
+++ b/tests/KrigingTest.cpp
@@ -47,7 +47,7 @@ TEST_CASE("workflow") {
     const double theta = 0.5;
     arma::vec theta_vec(X.n_cols);
     theta_vec.fill(theta);
-    return std::get<1>(ok.logLikelihoodEval(theta_vec, true, false));
+    return std::get<1>(ok.logLikelihoodFun(theta_vec, true, false));
   });
 }
 
@@ -71,7 +71,7 @@ TEST_CASE("logLikelihoodFun benchmark", "[.benchmark]") {
     theta_vec.fill(theta);
 
     BENCHMARK("Kriging::logLikelihoodFun#" + std::to_string(i)) {
-      return std::get<0>(ok.logLikelihoodEval(theta_vec, false, false));  //
+      return std::get<0>(ok.logLikelihoodFun(theta_vec, false, false));  //
     };
   });
 }
@@ -86,7 +86,7 @@ TEST_CASE("logLikelihoodGrad benchmark", "[.benchmark]") {
     theta_vec.fill(theta);
 
     BENCHMARK("Kriging::logLikelihoodGrad#" + std::to_string(i)) {
-      return std::get<1>(ok.logLikelihoodEval(theta_vec, true, false));  //
+      return std::get<1>(ok.logLikelihoodFun(theta_vec, true, false));  //
     };
   });
 }
\ No newline at end of file