Fix CUDA and mdlp library issues

Reviewed-on: #30

.devcontainer/Dockerfile

@@ -1,6 +1,6 @@
 FROM mcr.microsoft.com/devcontainers/cpp:ubuntu22.04
 
-ARG REINSTALL_CMAKE_VERSION_FROM_SOURCE="3.22.2"
+ARG REINSTALL_CMAKE_VERSION_FROM_SOURCE="3.29.3"
 
 # Optionally install the cmake for vcpkg
 COPY ./reinstall-cmake.sh /tmp/
@@ -23,7 +23,7 @@ RUN add-apt-repository ppa:ubuntu-toolchain-r/test
 RUN apt-get update
 
 # Install GCC 13.1
-RUN apt-get install -y gcc-13 g++-13
+RUN apt-get install -y gcc-13 g++-13 doxygen
 
 # Install lcov 2.1
 RUN wget --quiet https://github.com/linux-test-project/lcov/releases/download/v2.1/lcov-2.1.tar.gz && \

.gitmodules

@@ -1,8 +1,3 @@
-[submodule "lib/mdlp"]
-	path = lib/mdlp
-	url = https://github.com/rmontanana/mdlp
-	main = main
-	update = merge
 [submodule "lib/json"]
 	path = lib/json
 	url = https://github.com/nlohmann/json.git

.vscode/launch.json

@@ -14,11 +14,11 @@
             "type": "lldb",
             "request": "launch",
             "name": "test",
-            "program": "${workspaceFolder}/build_debug/tests/TestBayesNet",
+            "program": "${workspaceFolder}/build_Debug/tests/TestBayesNet",
             "args": [
-                "[Node]"
+                "[Network]"
             ],
-            "cwd": "${workspaceFolder}/build_debug/tests"
+            "cwd": "${workspaceFolder}/build_Debug/tests"
         },
         {
             "name": "(gdb) Launch",

CHANGELOG.md

@@ -19,9 +19,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Add tests to check the correct version of the mdlp, folding and json libraries.
 - Library documentation generated with Doxygen.
 - Link to documentation in the README.md.
+- Three types of smoothing the Bayesian Network OLD_LAPLACE, LAPLACE and CESTNIK.
 
 ### Internal
 
+- Ignore CUDA language in CMake CodeCoverage module.
+- Remove mdlp library from submodules and add it as a dependency.
 - Create library ShuffleArffFile to limit the number of samples with a parameter and shuffle them.
 - Refactor catch2 library location to test/lib
 - Refactor loadDataset function in tests.

CMakeLists.txt

@@ -1,7 +1,7 @@
 cmake_minimum_required(VERSION 3.20)
 
 project(BayesNet
-  VERSION 1.0.5.1
+  VERSION 1.0.6
   DESCRIPTION "Bayesian Network and basic classifiers Library."
   HOMEPAGE_URL "https://github.com/rmontanana/bayesnet"
   LANGUAGES CXX
@@ -26,7 +26,7 @@ set(CMAKE_EXPORT_COMPILE_COMMANDS          ON)
 set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${TORCH_CXX_FLAGS}")
 SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -pthread")
 set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -fprofile-arcs -ftest-coverage -fno-elide-constructors")
-set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -O3")
+set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -Ofast")
 if (NOT ${CMAKE_SYSTEM_NAME} MATCHES "Darwin")
     set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -fno-default-inline")
 endif()
@@ -49,11 +49,12 @@ if (CMAKE_BUILD_TYPE STREQUAL "Debug")
   set(CODE_COVERAGE ON)
 endif (CMAKE_BUILD_TYPE STREQUAL "Debug")
 
-
+get_property(LANGUAGES GLOBAL PROPERTY ENABLED_LANGUAGES)
+message(STATUS "Languages=${LANGUAGES}")
 if (CODE_COVERAGE)
   enable_testing()
   include(CodeCoverage)
-  MESSAGE("Code coverage enabled")
+  MESSAGE(STATUS "Code coverage enabled")
   SET(GCC_COVERAGE_LINK_FLAGS " ${GCC_COVERAGE_LINK_FLAGS} -lgcov --coverage")
 endif (CODE_COVERAGE)
 
@@ -63,9 +64,16 @@ endif (ENABLE_CLANG_TIDY)
 
 # External libraries - dependencies of BayesNet
 # ---------------------------------------------
+
+find_library(FImdlp NAMES libfimdlp.a REQUIRED)
+find_path(FImdlp_INCLUDE_DIRS REQUIRED NAMES fimdlp)
+
+message(STATUS "FImdlp=${FImdlp}")
+message(STATUS "FImdlp_INCLUDE_DIRS=${FImdlp_INCLUDE_DIRS}")
+
+
 # include(FetchContent)
 add_git_submodule("lib/json")
-add_git_submodule("lib/mdlp")
 
 # Subdirectories
 # --------------
@@ -75,7 +83,7 @@ add_subdirectory(bayesnet)
 # Testing
 # -------
 if (ENABLE_TESTING)
-MESSAGE("Testing enabled")
+MESSAGE(STATUS "Testing enabled")
   add_subdirectory(tests/lib/catch2)
   include(CTest)
   add_subdirectory(tests)
@@ -93,6 +101,7 @@ install(FILES ${CMAKE_BINARY_DIR}/configured_files/include/bayesnet/config.h DES
 # Documentation
 # -------------
 find_package(Doxygen)
+if (Doxygen_FOUND)
   set(DOC_DIR ${CMAKE_CURRENT_SOURCE_DIR}/docs)
   set(doxyfile_in ${DOC_DIR}/Doxyfile.in)
   set(doxyfile ${DOC_DIR}/Doxyfile)
@@ -100,3 +109,6 @@ configure_file(${doxyfile_in} ${doxyfile} @ONLY)
   doxygen_add_docs(doxygen
       WORKING_DIRECTORY ${DOC_DIR}
     CONFIG_FILE ${doxyfile})
+else (Doxygen_FOUND) 
+  MESSAGE("* Doxygen not found")
+endif (Doxygen_FOUND)

Makefile

@@ -12,7 +12,6 @@ plantuml = plantuml
 lcov = lcov
 genhtml = genhtml
 dot = dot
-n_procs = -j 16
 docsrcdir = docs/manual
 mansrcdir = docs/man3
 mandestdir = /usr/local/share/man
@@ -59,10 +58,10 @@ diagrams: ## Create an UML class diagram & depnendency of the project (diagrams/
 	@$(dot) -Tsvg $(f_debug)/dependency.dot.BayesNet -o $(f_diagrams)/dependency.svg
 
 buildd: ## Build the debug targets
-	cmake --build $(f_debug) -t $(app_targets) $(n_procs)
+	cmake --build $(f_debug) -t $(app_targets) --parallel $(CMAKE_BUILD_PARALLEL_LEVEL)
 
 buildr: ## Build the release targets
-	cmake --build $(f_release) -t $(app_targets) $(n_procs)
+	cmake --build $(f_release) -t $(app_targets) --parallel $(CMAKE_BUILD_PARALLEL_LEVEL)
 
 clean: ## Clean the tests info
 	@echo ">>> Cleaning Debug BayesNet tests...";
@@ -106,7 +105,7 @@ opt = ""
 test: ## Run tests (opt="-s") to verbose output the tests, (opt="-c='Test Maximum Spanning Tree'") to run only that section
 	@echo ">>> Running BayesNet tests...";
 	@$(MAKE) clean
-	@cmake --build $(f_debug) -t $(test_targets) $(n_procs)
+	@cmake --build $(f_debug) -t $(test_targets) --parallel $(CMAKE_BUILD_PARALLEL_LEVEL)
 	@for t in $(test_targets); do \
 		echo ">>> Running $$t...";\
 		if [ -f $(f_debug)/tests/$$t ]; then \
@@ -119,7 +118,7 @@ test: ## Run tests (opt="-s") to verbose output the tests, (opt="-c='Test Maximu
 
 coverage: ## Run tests and generate coverage report (build/index.html)
 	@echo ">>> Building tests with coverage..."
-	@which $(lcov) || (echo ">>> Please install lcov"; exit 1)
+	@which $(lcov) || (echo ">>ease install lcov"; exit 1)
 	@if [ ! -f $(f_debug)/tests/coverage.info ] ; then $(MAKE) test ; fi
 	@echo ">>> Building report..."
 	@cd $(f_debug)/tests; \

README.md

@@ -7,7 +7,7 @@
 [![Security Rating](https://sonarcloud.io/api/project_badges/measure?project=rmontanana_BayesNet&metric=security_rating)](https://sonarcloud.io/summary/new_code?id=rmontanana_BayesNet)
 [![Reliability Rating](https://sonarcloud.io/api/project_badges/measure?project=rmontanana_BayesNet&metric=reliability_rating)](https://sonarcloud.io/summary/new_code?id=rmontanana_BayesNet)
 ![Gitea Last Commit](https://img.shields.io/gitea/last-commit/rmontanana/bayesnet?gitea_url=https://gitea.rmontanana.es:3000&logo=gitea)
-[![Coverage Badge](https://img.shields.io/badge/Coverage-97,3%25-green)](html/index.html)
+[![Coverage Badge](https://img.shields.io/badge/Coverage-97,1%25-green)](html/index.html)
 
 Bayesian Network Classifiers using libtorch from scratch
 

bayesnet/BaseClassifier.h

@@ -8,16 +8,18 @@
 #include <vector>
 #include <torch/torch.h>
 #include <nlohmann/json.hpp>
+#include "bayesnet/network/Network.h"
+
 namespace bayesnet {
     enum status_t { NORMAL, WARNING, ERROR };
     class BaseClassifier {
     public:
         // X is nxm std::vector, y is nx1 std::vector
-        virtual BaseClassifier& fit(std::vector<std::vector<int>>& X, std::vector<int>& y, const std::vector<std::string>& features, const std::string& className, std::map<std::string, std::vector<int>>& states) = 0;
+        virtual BaseClassifier& fit(std::vector<std::vector<int>>& X, std::vector<int>& y, const std::vector<std::string>& features, const std::string& className, std::map<std::string, std::vector<int>>& states, const Smoothing_t smoothing) = 0;
         // X is nxm tensor, y is nx1 tensor
-        virtual BaseClassifier& fit(torch::Tensor& X, torch::Tensor& y, const std::vector<std::string>& features, const std::string& className, std::map<std::string, std::vector<int>>& states) = 0;
-        virtual BaseClassifier& fit(torch::Tensor& dataset, const std::vector<std::string>& features, const std::string& className, std::map<std::string, std::vector<int>>& states) = 0;
-        virtual BaseClassifier& fit(torch::Tensor& dataset, const std::vector<std::string>& features, const std::string& className, std::map<std::string, std::vector<int>>& states, const torch::Tensor& weights) = 0;
+        virtual BaseClassifier& fit(torch::Tensor& X, torch::Tensor& y, const std::vector<std::string>& features, const std::string& className, std::map<std::string, std::vector<int>>& states, const Smoothing_t smoothing) = 0;
+        virtual BaseClassifier& fit(torch::Tensor& dataset, const std::vector<std::string>& features, const std::string& className, std::map<std::string, std::vector<int>>& states, const Smoothing_t smoothing) = 0;
+        virtual BaseClassifier& fit(torch::Tensor& dataset, const std::vector<std::string>& features, const std::string& className, std::map<std::string, std::vector<int>>& states, const torch::Tensor& weights, const Smoothing_t smoothing) = 0;
         virtual ~BaseClassifier() = default;
         torch::Tensor virtual predict(torch::Tensor& X) = 0;
         std::vector<int> virtual predict(std::vector<std::vector<int >>& X) = 0;
@@ -39,7 +41,7 @@ namespace bayesnet {
         virtual void setHyperparameters(const nlohmann::json& hyperparameters) = 0;
         std::vector<std::string>& getValidHyperparameters() { return validHyperparameters; }
     protected:
-        virtual void trainModel(const torch::Tensor& weights) = 0;
+        virtual void trainModel(const torch::Tensor& weights, const Smoothing_t smoothing) = 0;
         std::vector<std::string> validHyperparameters;
     };
 }

bayesnet/CMakeLists.txt

@@ -1,12 +1,13 @@
 include_directories(
-    ${BayesNet_SOURCE_DIR}/lib/mdlp
+    ${BayesNet_SOURCE_DIR}/lib/mdlp/src
     ${BayesNet_SOURCE_DIR}/lib/folding
     ${BayesNet_SOURCE_DIR}/lib/json/include
     ${BayesNet_SOURCE_DIR}
     ${CMAKE_BINARY_DIR}/configured_files/include
+    ${FImdlp_INCLUDE_DIRS}
 )
 
 file(GLOB_RECURSE Sources "*.cc")
 
 add_library(BayesNet ${Sources})
-target_link_libraries(BayesNet mdlp "${TORCH_LIBRARIES}")
+target_link_libraries(BayesNet ${FImdlp} "${TORCH_LIBRARIES}")

bayesnet/classifiers/Classifier.cc

@@ -11,7 +11,7 @@
 namespace bayesnet {
     Classifier::Classifier(Network model) : model(model), m(0), n(0), metrics(Metrics()), fitted(false) {}
     const std::string CLASSIFIER_NOT_FITTED = "Classifier has not been fitted";
-    Classifier& Classifier::build(const std::vector<std::string>& features, const std::string& className, std::map<std::string, std::vector<int>>& states, const torch::Tensor& weights)
+    Classifier& Classifier::build(const std::vector<std::string>& features, const std::string& className, std::map<std::string, std::vector<int>>& states, const torch::Tensor& weights, const Smoothing_t smoothing)
     {
         this->features = features;
         this->className = className;
@@ -23,7 +23,7 @@ namespace bayesnet {
         metrics = Metrics(dataset, features, className, n_classes);
         model.initialize();
         buildModel(weights);
-        trainModel(weights);
+        trainModel(weights, smoothing);
         fitted = true;
         return *this;
     }
@@ -41,20 +41,20 @@ namespace bayesnet {
             throw std::runtime_error(oss.str());
         }
     }
-    void Classifier::trainModel(const torch::Tensor& weights)
+    void Classifier::trainModel(const torch::Tensor& weights, Smoothing_t smoothing)
     {
-        model.fit(dataset, weights, features, className, states);
+        model.fit(dataset, weights, features, className, states, smoothing);
     }
     // X is nxm where n is the number of features and m the number of samples
-    Classifier& Classifier::fit(torch::Tensor& X, torch::Tensor& y, const std::vector<std::string>& features, const std::string& className, std::map<std::string, std::vector<int>>& states)
+    Classifier& Classifier::fit(torch::Tensor& X, torch::Tensor& y, const std::vector<std::string>& features, const std::string& className, std::map<std::string, std::vector<int>>& states, const Smoothing_t smoothing)
     {
         dataset = X;
         buildDataset(y);
         const torch::Tensor weights = torch::full({ dataset.size(1) }, 1.0 / dataset.size(1), torch::kDouble);
-        return build(features, className, states, weights);
+        return build(features, className, states, weights, smoothing);
     }
     // X is nxm where n is the number of features and m the number of samples
-    Classifier& Classifier::fit(std::vector<std::vector<int>>& X, std::vector<int>& y, const std::vector<std::string>& features, const std::string& className, std::map<std::string, std::vector<int>>& states)
+    Classifier& Classifier::fit(std::vector<std::vector<int>>& X, std::vector<int>& y, const std::vector<std::string>& features, const std::string& className, std::map<std::string, std::vector<int>>& states, const Smoothing_t smoothing)
     {
         dataset = torch::zeros({ static_cast<int>(X.size()), static_cast<int>(X[0].size()) }, torch::kInt32);
         for (int i = 0; i < X.size(); ++i) {
@@ -63,18 +63,18 @@ namespace bayesnet {
         auto ytmp = torch::tensor(y, torch::kInt32);
         buildDataset(ytmp);
         const torch::Tensor weights = torch::full({ dataset.size(1) }, 1.0 / dataset.size(1), torch::kDouble);
-        return build(features, className, states, weights);
+        return build(features, className, states, weights, smoothing);
     }
-    Classifier& Classifier::fit(torch::Tensor& dataset, const std::vector<std::string>& features, const std::string& className, std::map<std::string, std::vector<int>>& states)
+    Classifier& Classifier::fit(torch::Tensor& dataset, const std::vector<std::string>& features, const std::string& className, std::map<std::string, std::vector<int>>& states, const Smoothing_t smoothing)
     {
         this->dataset = dataset;
         const torch::Tensor weights = torch::full({ dataset.size(1) }, 1.0 / dataset.size(1), torch::kDouble);
-        return build(features, className, states, weights);
+        return build(features, className, states, weights, smoothing);
     }
-    Classifier& Classifier::fit(torch::Tensor& dataset, const std::vector<std::string>& features, const std::string& className, std::map<std::string, std::vector<int>>& states, const torch::Tensor& weights)
+    Classifier& Classifier::fit(torch::Tensor& dataset, const std::vector<std::string>& features, const std::string& className, std::map<std::string, std::vector<int>>& states, const torch::Tensor& weights, const Smoothing_t smoothing)
     {
         this->dataset = dataset;
-        return build(features, className, states, weights);
+        return build(features, className, states, weights, smoothing);
     }
     void Classifier::checkFitParameters()
     {

bayesnet/classifiers/Classifier.h

@@ -8,7 +8,6 @@
 #define CLASSIFIER_H
 #include <torch/torch.h>
 #include "bayesnet/utils/BayesMetrics.h"
-#include "bayesnet/network/Network.h"
 #include "bayesnet/BaseClassifier.h"
 
 namespace bayesnet {
@@ -16,10 +15,10 @@ namespace bayesnet {
     public:
         Classifier(Network model);
         virtual ~Classifier() = default;
-        Classifier& fit(std::vector<std::vector<int>>& X, std::vector<int>& y, const std::vector<std::string>& features, const std::string& className, std::map<std::string, std::vector<int>>& states) override;
-        Classifier& fit(torch::Tensor& X, torch::Tensor& y, const std::vector<std::string>& features, const std::string& className, std::map<std::string, std::vector<int>>& states) override;
-        Classifier& fit(torch::Tensor& dataset, const std::vector<std::string>& features, const std::string& className, std::map<std::string, std::vector<int>>& states) override;
-        Classifier& fit(torch::Tensor& dataset, const std::vector<std::string>& features, const std::string& className, std::map<std::string, std::vector<int>>& states, const torch::Tensor& weights) override;
+        Classifier& fit(std::vector<std::vector<int>>& X, std::vector<int>& y, const std::vector<std::string>& features, const std::string& className, std::map<std::string, std::vector<int>>& states, const Smoothing_t smoothing) override;
+        Classifier& fit(torch::Tensor& X, torch::Tensor& y, const std::vector<std::string>& features, const std::string& className, std::map<std::string, std::vector<int>>& states, const Smoothing_t smoothing) override;
+        Classifier& fit(torch::Tensor& dataset, const std::vector<std::string>& features, const std::string& className, std::map<std::string, std::vector<int>>& states, const Smoothing_t smoothing) override;
+        Classifier& fit(torch::Tensor& dataset, const std::vector<std::string>& features, const std::string& className, std::map<std::string, std::vector<int>>& states, const torch::Tensor& weights, const Smoothing_t smoothing) override;
         void addNodes();
         int getNumberOfNodes() const override;
         int getNumberOfEdges() const override;
@@ -51,10 +50,10 @@ namespace bayesnet {
         std::vector<std::string> notes; // Used to store messages occurred during the fit process
         void checkFitParameters();
         virtual void buildModel(const torch::Tensor& weights) = 0;
-        void trainModel(const torch::Tensor& weights) override;
+        void trainModel(const torch::Tensor& weights, const Smoothing_t smoothing) override;
         void buildDataset(torch::Tensor& y);
     private:
-        Classifier& build(const std::vector<std::string>& features, const std::string& className, std::map<std::string, std::vector<int>>& states, const torch::Tensor& weights);
+        Classifier& build(const std::vector<std::string>& features, const std::string& className, std::map<std::string, std::vector<int>>& states, const torch::Tensor& weights, const Smoothing_t smoothing);
     };
 }
 #endif

bayesnet/classifiers/KDBLd.cc

@@ -8,7 +8,7 @@
 
 namespace bayesnet {
     KDBLd::KDBLd(int k) : KDB(k), Proposal(dataset, features, className) {}
-    KDBLd& KDBLd::fit(torch::Tensor& X_, torch::Tensor& y_, const std::vector<std::string>& features_, const std::string& className_, map<std::string, std::vector<int>>& states_)
+    KDBLd& KDBLd::fit(torch::Tensor& X_, torch::Tensor& y_, const std::vector<std::string>& features_, const std::string& className_, map<std::string, std::vector<int>>& states_, const Smoothing_t smoothing)
     {
         checkInput(X_, y_);
         features = features_;
@@ -19,7 +19,7 @@ namespace bayesnet {
         states = fit_local_discretization(y);
         // We have discretized the input data
         // 1st we need to fit the model to build the normal KDB structure, KDB::fit initializes the base Bayesian network
-        KDB::fit(dataset, features, className, states);
+        KDB::fit(dataset, features, className, states, smoothing);
         states = localDiscretizationProposal(states, model);
         return *this;
     }

bayesnet/classifiers/KDBLd.h

@@ -15,7 +15,7 @@ namespace bayesnet {
     public:
         explicit KDBLd(int k);
         virtual ~KDBLd() = default;
-        KDBLd& fit(torch::Tensor& X, torch::Tensor& y, const std::vector<std::string>& features, const std::string& className, map<std::string, std::vector<int>>& states) override;
+        KDBLd& fit(torch::Tensor& X, torch::Tensor& y, const std::vector<std::string>& features, const std::string& className, map<std::string, std::vector<int>>& states, const Smoothing_t smoothing) override;
         std::vector<std::string> graph(const std::string& name = "KDB") const override;
         torch::Tensor predict(torch::Tensor& X) override;
         static inline std::string version() { return "0.0.1"; };

bayesnet/classifiers/Proposal.cc

@@ -70,7 +70,7 @@ namespace bayesnet {
                 states[pFeatures[index]] = xStates;
             }
             const torch::Tensor weights = torch::full({ pDataset.size(1) }, 1.0 / pDataset.size(1), torch::kDouble);
-            model.fit(pDataset, weights, pFeatures, pClassName, states);
+            model.fit(pDataset, weights, pFeatures, pClassName, states, Smoothing_t::ORIGINAL);
         }
         return states;
     }

bayesnet/classifiers/Proposal.h

@@ -9,7 +9,7 @@
 #include <string>
 #include <map>
 #include <torch/torch.h>
-#include <CPPFImdlp.h>
+#include <fimdlp/CPPFImdlp.h>
 #include "bayesnet/network/Network.h"
 #include "Classifier.h"
 

bayesnet/classifiers/SPODELd.cc

@@ -8,25 +8,25 @@
 
 namespace bayesnet {
     SPODELd::SPODELd(int root) : SPODE(root), Proposal(dataset, features, className) {}
-    SPODELd& SPODELd::fit(torch::Tensor& X_, torch::Tensor& y_, const std::vector<std::string>& features_, const std::string& className_, map<std::string, std::vector<int>>& states_)
+    SPODELd& SPODELd::fit(torch::Tensor& X_, torch::Tensor& y_, const std::vector<std::string>& features_, const std::string& className_, map<std::string, std::vector<int>>& states_, const Smoothing_t smoothing)
     {
         checkInput(X_, y_);
         Xf = X_;
         y = y_;
-        return commonFit(features_, className_, states_);
+        return commonFit(features_, className_, states_, smoothing);
     }
 
-    SPODELd& SPODELd::fit(torch::Tensor& dataset, const std::vector<std::string>& features_, const std::string& className_, map<std::string, std::vector<int>>& states_)
+    SPODELd& SPODELd::fit(torch::Tensor& dataset, const std::vector<std::string>& features_, const std::string& className_, map<std::string, std::vector<int>>& states_, const Smoothing_t smoothing)
     {
         if (!torch::is_floating_point(dataset)) {
             throw std::runtime_error("Dataset must be a floating point tensor");
         }
         Xf = dataset.index({ torch::indexing::Slice(0, dataset.size(0) - 1), "..." }).clone();
         y = dataset.index({ -1, "..." }).clone().to(torch::kInt32);
-        return commonFit(features_, className_, states_);
+        return commonFit(features_, className_, states_, smoothing);
     }
 
-    SPODELd& SPODELd::commonFit(const std::vector<std::string>& features_, const std::string& className_, map<std::string, std::vector<int>>& states_)
+    SPODELd& SPODELd::commonFit(const std::vector<std::string>& features_, const std::string& className_, map<std::string, std::vector<int>>& states_, const Smoothing_t smoothing)
     {
         features = features_;
         className = className_;
@@ -34,7 +34,7 @@ namespace bayesnet {
         states = fit_local_discretization(y);
         // We have discretized the input data
         // 1st we need to fit the model to build the normal SPODE structure, SPODE::fit initializes the base Bayesian network
-        SPODE::fit(dataset, features, className, states);
+        SPODE::fit(dataset, features, className, states, smoothing);
         states = localDiscretizationProposal(states, model);
         return *this;
     }

bayesnet/classifiers/SPODELd.h

@@ -14,10 +14,10 @@ namespace bayesnet {
     public:
         explicit SPODELd(int root);
         virtual ~SPODELd() = default;
-        SPODELd& fit(torch::Tensor& X, torch::Tensor& y, const std::vector<std::string>& features, const std::string& className, map<std::string, std::vector<int>>& states) override;
-        SPODELd& fit(torch::Tensor& dataset, const std::vector<std::string>& features, const std::string& className, map<std::string, std::vector<int>>& states) override;
-        SPODELd& commonFit(const std::vector<std::string>& features, const std::string& className, map<std::string, std::vector<int>>& states);
-        std::vector<std::string> graph(const std::string& name = "SPODE") const override;
+        SPODELd& fit(torch::Tensor& X, torch::Tensor& y, const std::vector<std::string>& features, const std::string& className, map<std::string, std::vector<int>>& states, const Smoothing_t smoothing) override;
+        SPODELd& fit(torch::Tensor& dataset, const std::vector<std::string>& features, const std::string& className, map<std::string, std::vector<int>>& states, const Smoothing_t smoothing) override;
+        SPODELd& commonFit(const std::vector<std::string>& features, const std::string& className, map<std::string, std::vector<int>>& states, const Smoothing_t smoothing);
+        std::vector<std::string> graph(const std::string& name = "SPODELd") const override;
         torch::Tensor predict(torch::Tensor& X) override;
         static inline std::string version() { return "0.0.1"; };
     };

bayesnet/classifiers/TANLd.cc

@@ -8,7 +8,7 @@
 
 namespace bayesnet {
     TANLd::TANLd() : TAN(), Proposal(dataset, features, className) {}
-    TANLd& TANLd::fit(torch::Tensor& X_, torch::Tensor& y_, const std::vector<std::string>& features_, const std::string& className_, map<std::string, std::vector<int>>& states_)
+    TANLd& TANLd::fit(torch::Tensor& X_, torch::Tensor& y_, const std::vector<std::string>& features_, const std::string& className_, map<std::string, std::vector<int>>& states_, const Smoothing_t smoothing)
     {
         checkInput(X_, y_);
         features = features_;
@@ -19,7 +19,7 @@ namespace bayesnet {
         states = fit_local_discretization(y);
         // We have discretized the input data
         // 1st we need to fit the model to build the normal TAN structure, TAN::fit initializes the base Bayesian network
-        TAN::fit(dataset, features, className, states);
+        TAN::fit(dataset, features, className, states, smoothing);
         states = localDiscretizationProposal(states, model);
         return *this;
 

bayesnet/classifiers/TANLd.h

@@ -15,10 +15,9 @@ namespace bayesnet {
     public:
         TANLd();
         virtual ~TANLd() = default;
-        TANLd& fit(torch::Tensor& X, torch::Tensor& y, const std::vector<std::string>& features, const std::string& className, map<std::string, std::vector<int>>& states) override;
-        std::vector<std::string> graph(const std::string& name = "TAN") const override;
+        TANLd& fit(torch::Tensor& X, torch::Tensor& y, const std::vector<std::string>& features, const std::string& className, map<std::string, std::vector<int>>& states, const Smoothing_t smoothing) override;
+        std::vector<std::string> graph(const std::string& name = "TANLd") const override;
         torch::Tensor predict(torch::Tensor& X) override;
-        static inline std::string version() { return "0.0.1"; };
     };
 }
 #endif // !TANLD_H

bayesnet/ensembles/AODELd.cc

@@ -10,7 +10,7 @@ namespace bayesnet {
     AODELd::AODELd(bool predict_voting) : Ensemble(predict_voting), Proposal(dataset, features, className)
     {
     }
-    AODELd& AODELd::fit(torch::Tensor& X_, torch::Tensor& y_, const std::vector<std::string>& features_, const std::string& className_, map<std::string, std::vector<int>>& states_)
+    AODELd& AODELd::fit(torch::Tensor& X_, torch::Tensor& y_, const std::vector<std::string>& features_, const std::string& className_, map<std::string, std::vector<int>>& states_, const Smoothing_t smoothing)
     {
         checkInput(X_, y_);
         features = features_;
@@ -21,7 +21,7 @@ namespace bayesnet {
         states = fit_local_discretization(y);
         // We have discretized the input data
         // 1st we need to fit the model to build the normal TAN structure, TAN::fit initializes the base Bayesian network
-        Ensemble::fit(dataset, features, className, states);
+        Ensemble::fit(dataset, features, className, states, smoothing);
         return *this;
 
     }
@@ -34,10 +34,10 @@ namespace bayesnet {
         n_models = models.size();
         significanceModels = std::vector<double>(n_models, 1.0);
     }
-    void AODELd::trainModel(const torch::Tensor& weights)
+    void AODELd::trainModel(const torch::Tensor& weights, const Smoothing_t smoothing)
     {
         for (const auto& model : models) {
-            model->fit(Xf, y, features, className, states);
+            model->fit(Xf, y, features, className, states, smoothing);
         }
     }
     std::vector<std::string> AODELd::graph(const std::string& name) const

bayesnet/ensembles/AODELd.h

@@ -15,10 +15,10 @@ namespace bayesnet {
     public:
         AODELd(bool predict_voting = true);
         virtual ~AODELd() = default;
-        AODELd& fit(torch::Tensor& X_, torch::Tensor& y_, const std::vector<std::string>& features_, const std::string& className_, map<std::string, std::vector<int>>& states_) override;
+        AODELd& fit(torch::Tensor& X_, torch::Tensor& y_, const std::vector<std::string>& features_, const std::string& className_, map<std::string, std::vector<int>>& states_, const Smoothing_t smoothing) override;
         std::vector<std::string> graph(const std::string& name = "AODELd") const override;
     protected:
-        void trainModel(const torch::Tensor& weights) override;
+        void trainModel(const torch::Tensor& weights, const Smoothing_t smoothing) override;
         void buildModel(const torch::Tensor& weights) override;
     };
 }

bayesnet/ensembles/BoostA2DE.cc

@@ -19,7 +19,7 @@ namespace bayesnet {
     BoostA2DE::BoostA2DE(bool predict_voting) : Boost(predict_voting)
     {
     }
-    std::vector<int> BoostA2DE::initializeModels()
+    std::vector<int> BoostA2DE::initializeModels(const Smoothing_t smoothing)
     {
         torch::Tensor weights_ = torch::full({ m }, 1.0 / m, torch::kFloat64);
         std::vector<int> featuresSelected = featureSelection(weights_);
@@ -32,7 +32,7 @@ namespace bayesnet {
             for (int j = i + 1; j < featuresSelected.size(); j++) {
                 auto parents = { featuresSelected[i], featuresSelected[j] };
                 std::unique_ptr<Classifier> model = std::make_unique<SPnDE>(parents);
-                model->fit(dataset, features, className, states, weights_);
+                model->fit(dataset, features, className, states, weights_, smoothing);
                 models.push_back(std::move(model));
                 significanceModels.push_back(1.0); // They will be updated later in trainModel
                 n_models++;
@@ -41,7 +41,7 @@ namespace bayesnet {
         notes.push_back("Used features in initialization: " + std::to_string(featuresSelected.size()) + " of " + std::to_string(features.size()) + " with " + select_features_algorithm);
         return featuresSelected;
     }
-    void BoostA2DE::trainModel(const torch::Tensor& weights)
+    void BoostA2DE::trainModel(const torch::Tensor& weights, const Smoothing_t smoothing)
     {
         //
         // Logging setup
@@ -58,7 +58,7 @@ namespace bayesnet {
         bool finished = false;
         std::vector<int> featuresUsed;
         if (selectFeatures) {
-            featuresUsed = initializeModels();
+            featuresUsed = initializeModels(smoothing);
             auto ypred = predict(X_train);
             std::tie(weights_, alpha_t, finished) = update_weights(y_train, ypred, weights_);
             // Update significance of the models
@@ -96,7 +96,7 @@ namespace bayesnet {
                 pairSelection.erase(pairSelection.begin());
                 std::unique_ptr<Classifier> model;
                 model = std::make_unique<SPnDE>(std::vector<int>({ feature_pair.first, feature_pair.second }));
-                model->fit(dataset, features, className, states, weights_);
+                model->fit(dataset, features, className, states, weights_, smoothing);
                 alpha_t = 0.0;
                 if (!block_update) {
                     auto ypred = model->predict(X_train);

bayesnet/ensembles/BoostA2DE.h

@@ -17,9 +17,9 @@ namespace bayesnet {
         virtual ~BoostA2DE() = default;
         std::vector<std::string> graph(const std::string& title = "BoostA2DE") const override;
     protected:
-        void trainModel(const torch::Tensor& weights) override;
+        void trainModel(const torch::Tensor& weights, const Smoothing_t smoothing) override;
     private:
-        std::vector<int> initializeModels();
+        std::vector<int> initializeModels(const Smoothing_t smoothing);
     };
 }
 #endif

bayesnet/ensembles/BoostAODE.cc

@@ -16,13 +16,13 @@ namespace bayesnet {
     BoostAODE::BoostAODE(bool predict_voting) : Boost(predict_voting)
     {
     }
-    std::vector<int> BoostAODE::initializeModels()
+    std::vector<int> BoostAODE::initializeModels(const Smoothing_t smoothing)
     {
         torch::Tensor weights_ = torch::full({ m }, 1.0 / m, torch::kFloat64);
         std::vector<int> featuresSelected = featureSelection(weights_);
         for (const int& feature : featuresSelected) {
             std::unique_ptr<Classifier> model = std::make_unique<SPODE>(feature);
-            model->fit(dataset, features, className, states, weights_);
+            model->fit(dataset, features, className, states, weights_, smoothing);
             models.push_back(std::move(model));
             significanceModels.push_back(1.0); // They will be updated later in trainModel
             n_models++;
@@ -30,7 +30,7 @@ namespace bayesnet {
         notes.push_back("Used features in initialization: " + std::to_string(featuresSelected.size()) + " of " + std::to_string(features.size()) + " with " + select_features_algorithm);
         return featuresSelected;
     }
-    void BoostAODE::trainModel(const torch::Tensor& weights)
+    void BoostAODE::trainModel(const torch::Tensor& weights, const Smoothing_t smoothing)
     {
         //
         // Logging setup
@@ -47,7 +47,7 @@ namespace bayesnet {
         bool finished = false;
         std::vector<int> featuresUsed;
         if (selectFeatures) {
-            featuresUsed = initializeModels();
+            featuresUsed = initializeModels(smoothing);
             auto ypred = predict(X_train);
             std::tie(weights_, alpha_t, finished) = update_weights(y_train, ypred, weights_);
             // Update significance of the models
@@ -89,7 +89,7 @@ namespace bayesnet {
                 featureSelection.erase(featureSelection.begin());
                 std::unique_ptr<Classifier> model;
                 model = std::make_unique<SPODE>(feature);
-                model->fit(dataset, features, className, states, weights_);
+                model->fit(dataset, features, className, states, weights_, smoothing);
                 alpha_t = 0.0;
                 if (!block_update) {
                     auto ypred = model->predict(X_train);

bayesnet/ensembles/BoostAODE.h

@@ -18,9 +18,9 @@ namespace bayesnet {
         virtual ~BoostAODE() = default;
         std::vector<std::string> graph(const std::string& title = "BoostAODE") const override;
     protected:
-        void trainModel(const torch::Tensor& weights) override;
+        void trainModel(const torch::Tensor& weights, const Smoothing_t smoothing) override;
     private:
-        std::vector<int> initializeModels();
+        std::vector<int> initializeModels(const Smoothing_t smoothing);
     };
 }
 #endif

bayesnet/ensembles/Ensemble.cc

@@ -3,22 +3,21 @@
 // SPDX-FileType: SOURCE
 // SPDX-License-Identifier: MIT
 // ***************************************************************
-
 #include "Ensemble.h"
+#include "bayesnet/utils/CountingSemaphore.h"
 
 namespace bayesnet {
 
     Ensemble::Ensemble(bool predict_voting) : Classifier(Network()), n_models(0), predict_voting(predict_voting)
     {
-
     };
     const std::string ENSEMBLE_NOT_FITTED = "Ensemble has not been fitted";
-    void Ensemble::trainModel(const torch::Tensor& weights)
+    void Ensemble::trainModel(const torch::Tensor& weights, const Smoothing_t smoothing)
     {
         n_models = models.size();
         for (auto i = 0; i < n_models; ++i) {
             // fit with std::vectors
-            models[i]->fit(dataset, features, className, states);
+            models[i]->fit(dataset, features, className, states, smoothing);
         }
     }
     std::vector<int> Ensemble::compute_arg_max(std::vector<std::vector<double>>& X)
@@ -85,17 +84,9 @@ namespace bayesnet {
     {
         auto n_states = models[0]->getClassNumStates();
         torch::Tensor y_pred = torch::zeros({ X.size(1), n_states }, torch::kFloat32);
-        auto threads{ std::vector<std::thread>() };
-        std::mutex mtx;
         for (auto i = 0; i < n_models; ++i) {
-            threads.push_back(std::thread([&, i]() {
             auto ypredict = models[i]->predict_proba(X);
-                std::lock_guard<std::mutex> lock(mtx);
             y_pred += ypredict * significanceModels[i];
-                }));
-        }
-        for (auto& thread : threads) {
-            thread.join();
         }
         auto sum = std::reduce(significanceModels.begin(), significanceModels.end());
         y_pred /= sum;
@@ -105,23 +96,15 @@ namespace bayesnet {
     {
         auto n_states = models[0]->getClassNumStates();
         std::vector<std::vector<double>> y_pred(X[0].size(), std::vector<double>(n_states, 0.0));
-        auto threads{ std::vector<std::thread>() };
-        std::mutex mtx;
         for (auto i = 0; i < n_models; ++i) {
-            threads.push_back(std::thread([&, i]() {
             auto ypredict = models[i]->predict_proba(X);
             assert(ypredict.size() == y_pred.size());
             assert(ypredict[0].size() == y_pred[0].size());
-                std::lock_guard<std::mutex> lock(mtx);
             // Multiply each prediction by the significance of the model and then add it to the final prediction
             for (auto j = 0; j < ypredict.size(); ++j) {
                 std::transform(y_pred[j].begin(), y_pred[j].end(), ypredict[j].begin(), y_pred[j].begin(),
                     [significanceModels = significanceModels[i]](double x, double y) { return x + y * significanceModels; });
             }
-                }));
-        }
-        for (auto& thread : threads) {
-            thread.join();
         }
         auto sum = std::reduce(significanceModels.begin(), significanceModels.end());
         //Divide each element of the prediction by the sum of the significances
@@ -141,17 +124,9 @@ namespace bayesnet {
     {
         // Build a m x n_models tensor with the predictions of each model
         torch::Tensor y_pred = torch::zeros({ X.size(1), n_models }, torch::kInt32);
-        auto threads{ std::vector<std::thread>() };
-        std::mutex mtx;
         for (auto i = 0; i < n_models; ++i) {
-            threads.push_back(std::thread([&, i]() {
             auto ypredict = models[i]->predict(X);
-                std::lock_guard<std::mutex> lock(mtx);
             y_pred.index_put_({ "...", i }, ypredict);
-                }));
-        }
-        for (auto& thread : threads) {
-            thread.join();
         }
         return voting(y_pred);
     }

bayesnet/ensembles/Ensemble.h

@@ -46,7 +46,7 @@ namespace bayesnet {
         unsigned n_models;
         std::vector<std::unique_ptr<Classifier>> models;
         std::vector<double> significanceModels;
-        void trainModel(const torch::Tensor& weights) override;
+        void trainModel(const torch::Tensor& weights, const Smoothing_t smoothing) override;
         bool predict_voting;
     };
 }

bayesnet/network/Network.cc

@@ -5,20 +5,20 @@
 // ***************************************************************
 
 #include <thread>
-#include <mutex>
 #include <sstream>
+#include <numeric>
+#include <algorithm>
 #include "Network.h"
 #include "bayesnet/utils/bayesnetUtils.h"
+#include "bayesnet/utils/CountingSemaphore.h"
+#include <pthread.h>
+#include <fstream>
 namespace bayesnet {
-    Network::Network() : fitted{ false }, maxThreads{ 0.95 }, classNumStates{ 0 }, laplaceSmoothing{ 0 }
+    Network::Network() : fitted{ false }, classNumStates{ 0 }
     {
     }
-    Network::Network(float maxT) : fitted{ false }, maxThreads{ maxT }, classNumStates{ 0 }, laplaceSmoothing{ 0 }
-    {
-
-    }
-    Network::Network(const Network& other) : laplaceSmoothing(other.laplaceSmoothing), features(other.features), className(other.className), classNumStates(other.getClassNumStates()),
-        maxThreads(other.getMaxThreads()), fitted(other.fitted), samples(other.samples)
+    Network::Network(const Network& other) : features(other.features), className(other.className), classNumStates(other.getClassNumStates()),
+        fitted(other.fitted), samples(other.samples)
     {
         if (samples.defined())
             samples = samples.clone();
@@ -35,16 +35,15 @@ namespace bayesnet {
         nodes.clear();
         samples = torch::Tensor();
     }
-    float Network::getMaxThreads() const
-    {
-        return maxThreads;
-    }
     torch::Tensor& Network::getSamples()
     {
         return samples;
     }
     void Network::addNode(const std::string& name)
     {
+        if (fitted) {
+            throw std::invalid_argument("Cannot add node to a fitted network. Initialize first.");
+        }
         if (name == "") {
             throw std::invalid_argument("Node name cannot be empty");
         }
@@ -94,12 +93,21 @@ namespace bayesnet {
     }
     void Network::addEdge(const std::string& parent, const std::string& child)
     {
+        if (fitted) {
+            throw std::invalid_argument("Cannot add edge to a fitted network. Initialize first.");
+        }
         if (nodes.find(parent) == nodes.end()) {
             throw std::invalid_argument("Parent node " + parent + " does not exist");
         }
         if (nodes.find(child) == nodes.end()) {
             throw std::invalid_argument("Child node " + child + " does not exist");
         }
+        // Check if the edge is already in the graph
+        for (auto& node : nodes[parent]->getChildren()) {
+            if (node->getName() == child) {
+                throw std::invalid_argument("Edge " + parent + " -> " + child + " already exists");
+            }
+        }
         // Temporarily add edge to check for cycles
         nodes[parent]->addChild(nodes[child].get());
         nodes[child]->addParent(nodes[parent].get());
@@ -155,7 +163,7 @@ namespace bayesnet {
         classNumStates = nodes.at(className)->getNumStates();
     }
     // X comes in nxm, where n is the number of features and m the number of samples
-    void Network::fit(const torch::Tensor& X, const torch::Tensor& y, const torch::Tensor& weights, const std::vector<std::string>& featureNames, const std::string& className, const std::map<std::string, std::vector<int>>& states)
+    void Network::fit(const torch::Tensor& X, const torch::Tensor& y, const torch::Tensor& weights, const std::vector<std::string>& featureNames, const std::string& className, const std::map<std::string, std::vector<int>>& states, const Smoothing_t smoothing)
     {
         checkFitData(X.size(1), X.size(0), y.size(0), featureNames, className, states, weights);
         this->className = className;
@@ -164,17 +172,17 @@ namespace bayesnet {
         for (int i = 0; i < featureNames.size(); ++i) {
             auto row_feature = X.index({ i, "..." });
         }
-        completeFit(states, weights);
+        completeFit(states, weights, smoothing);
     }
-    void Network::fit(const torch::Tensor& samples, const torch::Tensor& weights, const std::vector<std::string>& featureNames, const std::string& className, const std::map<std::string, std::vector<int>>& states)
+    void Network::fit(const torch::Tensor& samples, const torch::Tensor& weights, const std::vector<std::string>& featureNames, const std::string& className, const std::map<std::string, std::vector<int>>& states, const Smoothing_t smoothing)
     {
         checkFitData(samples.size(1), samples.size(0) - 1, samples.size(1), featureNames, className, states, weights);
         this->className = className;
         this->samples = samples;
-        completeFit(states, weights);
+        completeFit(states, weights, smoothing);
     }
     // input_data comes in nxm, where n is the number of features and m the number of samples
-    void Network::fit(const std::vector<std::vector<int>>& input_data, const std::vector<int>& labels, const std::vector<double>& weights_, const std::vector<std::string>& featureNames, const std::string& className, const std::map<std::string, std::vector<int>>& states)
+    void Network::fit(const std::vector<std::vector<int>>& input_data, const std::vector<int>& labels, const std::vector<double>& weights_, const std::vector<std::string>& featureNames, const std::string& className, const std::map<std::string, std::vector<int>>& states, const Smoothing_t smoothing)
     {
         const torch::Tensor weights = torch::tensor(weights_, torch::kFloat64);
         checkFitData(input_data[0].size(), input_data.size(), labels.size(), featureNames, className, states, weights);
@@ -185,21 +193,57 @@ namespace bayesnet {
             samples.index_put_({ i, "..." }, torch::tensor(input_data[i], torch::kInt32));
         }
         samples.index_put_({ -1, "..." }, torch::tensor(labels, torch::kInt32));
-        completeFit(states, weights);
+        completeFit(states, weights, smoothing);
     }
-    void Network::completeFit(const std::map<std::string, std::vector<int>>& states, const torch::Tensor& weights)
+    void Network::completeFit(const std::map<std::string, std::vector<int>>& states, const torch::Tensor& weights, const Smoothing_t smoothing)
     {
         setStates(states);
-        laplaceSmoothing = 1.0 / samples.size(1); // To use in CPT computation
         std::vector<std::thread> threads;
+        auto& semaphore = CountingSemaphore::getInstance();
+        const double n_samples = static_cast<double>(samples.size(1));
+        auto worker = [&](std::pair<const std::string, std::unique_ptr<Node>>& node, int i) {
+            std::string threadName = "FitWorker-" + std::to_string(i);
+#if defined(__linux__)
+            pthread_setname_np(pthread_self(), threadName.c_str());
+#else
+            pthread_setname_np(threadName.c_str());
+#endif
+            double numStates = static_cast<double>(node.second->getNumStates());
+            double smoothing_factor = 0.0;
+            switch (smoothing) {
+                case Smoothing_t::ORIGINAL:
+                    smoothing_factor = 1.0 / n_samples;
+                    break;
+                case Smoothing_t::LAPLACE:
+                    smoothing_factor = 1.0;
+                    break;
+                case Smoothing_t::CESTNIK:
+                    smoothing_factor = 1 / numStates;
+                    break;
+                default:
+                    throw std::invalid_argument("Smoothing method not recognized " + std::to_string(static_cast<int>(smoothing)));
+            }
+            node.second->computeCPT(samples, features, smoothing_factor, weights);
+            semaphore.release();
+            };
+        int i = 0;
         for (auto& node : nodes) {
-            threads.emplace_back([this, &node, &weights]() {
-                node.second->computeCPT(samples, features, laplaceSmoothing, weights);
-                });
+            semaphore.acquire();
+            threads.emplace_back(worker, std::ref(node), i++);
         }
         for (auto& thread : threads) {
             thread.join();
         }
+        // std::fstream file;
+        // file.open("cpt.txt", std::fstream::out | std::fstream::app);
+        // file << std::string(80, '*') << std::endl;
+        // for (const auto& item : graph("Test")) {
+        //     file << item << std::endl;
+        // }
+        // file << std::string(80, '-') << std::endl;
+        // file << dump_cpt() << std::endl;
+        // file << std::string(80, '=') << std::endl;
+        // file.close();
         fitted = true;
     }
     torch::Tensor Network::predict_tensor(const torch::Tensor& samples, const bool proba)
@@ -207,15 +251,39 @@ namespace bayesnet {
         if (!fitted) {
             throw std::logic_error("You must call fit() before calling predict()");
         }
+        // Ensure the sample size is equal to the number of features
+        if (samples.size(0) != features.size() - 1) {
+            throw std::invalid_argument("(T) Sample size (" + std::to_string(samples.size(0)) +
+                ") does not match the number of features (" + std::to_string(features.size() - 1) + ")");
+        }
         torch::Tensor result;
+        std::vector<std::thread> threads;
+        std::mutex mtx;
+        auto& semaphore = CountingSemaphore::getInstance();
         result = torch::zeros({ samples.size(1), classNumStates }, torch::kFloat64);
-        for (int i = 0; i < samples.size(1); ++i) {
-            const torch::Tensor sample = samples.index({ "...", i });
+        auto worker = [&](const torch::Tensor& sample, int i) {
+            std::string threadName = "PredictWorker-" + std::to_string(i);
+#if defined(__linux__)
+            pthread_setname_np(pthread_self(), threadName.c_str());
+#else
+            pthread_setname_np(threadName.c_str());
+#endif
             auto psample = predict_sample(sample);
             auto temp = torch::tensor(psample, torch::kFloat64);
-            //            result.index_put_({ i, "..." }, torch::tensor(predict_sample(sample), torch::kFloat64));
+            {
+                std::lock_guard<std::mutex> lock(mtx);
                 result.index_put_({ i, "..." }, temp);
             }
+            semaphore.release();
+            };
+        for (int i = 0; i < samples.size(1); ++i) {
+            semaphore.acquire();
+            const torch::Tensor sample = samples.index({ "...", i });
+            threads.emplace_back(worker, sample, i);
+        }
+        for (auto& thread : threads) {
+            thread.join();
+        }
         if (proba)
             return result;
         return result.argmax(1);
@@ -239,18 +307,38 @@ namespace bayesnet {
         if (!fitted) {
             throw std::logic_error("You must call fit() before calling predict()");
         }
-        std::vector<int> predictions;
+        // Ensure the sample size is equal to the number of features
+        if (tsamples.size() != features.size() - 1) {
+            throw std::invalid_argument("(V) Sample size (" + std::to_string(tsamples.size()) +
+                ") does not match the number of features (" + std::to_string(features.size() - 1) + ")");
+        }
+        std::vector<int> predictions(tsamples[0].size(), 0);
         std::vector<int> sample;
+        std::vector<std::thread> threads;
+        auto& semaphore = CountingSemaphore::getInstance();
+        auto worker = [&](const std::vector<int>& sample, const int row, int& prediction) {
+            std::string threadName = "(V)PWorker-" + std::to_string(row);
+#if defined(__linux__)
+            pthread_setname_np(pthread_self(), threadName.c_str());
+#else
+            pthread_setname_np(threadName.c_str());
+#endif
+            auto classProbabilities = predict_sample(sample);
+            auto maxElem = max_element(classProbabilities.begin(), classProbabilities.end());
+            int predictedClass = distance(classProbabilities.begin(), maxElem);
+            prediction = predictedClass;
+            semaphore.release();
+            };
         for (int row = 0; row < tsamples[0].size(); ++row) {
             sample.clear();
             for (int col = 0; col < tsamples.size(); ++col) {
                 sample.push_back(tsamples[col][row]);
             }
-            std::vector<double> classProbabilities = predict_sample(sample);
-            // Find the class with the maximum posterior probability
-            auto maxElem = max_element(classProbabilities.begin(), classProbabilities.end());
-            int predictedClass = distance(classProbabilities.begin(), maxElem);
-            predictions.push_back(predictedClass);
+            semaphore.acquire();
+            threads.emplace_back(worker, sample, row, std::ref(predictions[row]));
+        }
+        for (auto& thread : threads) {
+            thread.join();
         }
         return predictions;
     }
@@ -261,14 +349,36 @@ namespace bayesnet {
         if (!fitted) {
             throw std::logic_error("You must call fit() before calling predict_proba()");
         }
-        std::vector<std::vector<double>> predictions;
+        // Ensure the sample size is equal to the number of features
+        if (tsamples.size() != features.size() - 1) {
+            throw std::invalid_argument("(V) Sample size (" + std::to_string(tsamples.size()) +
+                ") does not match the number of features (" + std::to_string(features.size() - 1) + ")");
+        }
+        std::vector<std::vector<double>> predictions(tsamples[0].size(), std::vector<double>(classNumStates, 0.0));
         std::vector<int> sample;
+        std::vector<std::thread> threads;
+        auto& semaphore = CountingSemaphore::getInstance();
+        auto worker = [&](const std::vector<int>& sample, int row, std::vector<double>& predictions) {
+            std::string threadName = "(V)PWorker-" + std::to_string(row);
+#if defined(__linux__)
+            pthread_setname_np(pthread_self(), threadName.c_str());
+#else
+            pthread_setname_np(threadName.c_str());
+#endif
+            std::vector<double> classProbabilities = predict_sample(sample);
+            predictions = classProbabilities;
+            semaphore.release();
+            };
         for (int row = 0; row < tsamples[0].size(); ++row) {
             sample.clear();
             for (int col = 0; col < tsamples.size(); ++col) {
                 sample.push_back(tsamples[col][row]);
             }
-            predictions.push_back(predict_sample(sample));
+            semaphore.acquire();
+            threads.emplace_back(worker, sample, row, std::ref(predictions[row]));
+        }
+        for (auto& thread : threads) {
+            thread.join();
         }
         return predictions;
     }
@@ -286,11 +396,6 @@ namespace bayesnet {
     // Return 1xn std::vector of probabilities
     std::vector<double> Network::predict_sample(const std::vector<int>& sample)
     {
-        // Ensure the sample size is equal to the number of features
-        if (sample.size() != features.size() - 1) {
-            throw std::invalid_argument("Sample size (" + std::to_string(sample.size()) +
-                ") does not match the number of features (" + std::to_string(features.size() - 1) + ")");
-        }
         std::map<std::string, int> evidence;
         for (int i = 0; i < sample.size(); ++i) {
             evidence[features[i]] = sample[i];
@@ -300,44 +405,26 @@ namespace bayesnet {
     // Return 1xn std::vector of probabilities
     std::vector<double> Network::predict_sample(const torch::Tensor& sample)
     {
-        // Ensure the sample size is equal to the number of features
-        if (sample.size(0) != features.size() - 1) {
-            throw std::invalid_argument("Sample size (" + std::to_string(sample.size(0)) +
-                ") does not match the number of features (" + std::to_string(features.size() - 1) + ")");
-        }
         std::map<std::string, int> evidence;
         for (int i = 0; i < sample.size(0); ++i) {
             evidence[features[i]] = sample[i].item<int>();
         }
         return exactInference(evidence);
     }
-    double Network::computeFactor(std::map<std::string, int>& completeEvidence)
-    {
-        double result = 1.0;
-        for (auto& node : getNodes()) {
-            result *= node.second->getFactorValue(completeEvidence);
-        }
-        return result;
-    }
     std::vector<double> Network::exactInference(std::map<std::string, int>& evidence)
     {
         std::vector<double> result(classNumStates, 0.0);
-        std::vector<std::thread> threads;
-        std::mutex mtx;
-        for (int i = 0; i < classNumStates; ++i) {
-            threads.emplace_back([this, &result, &evidence, i, &mtx]() {
         auto completeEvidence = std::map<std::string, int>(evidence);
+        for (int i = 0; i < classNumStates; ++i) {
             completeEvidence[getClassName()] = i;
-                double factor = computeFactor(completeEvidence);
-                std::lock_guard<std::mutex> lock(mtx);
-                result[i] = factor;
-                });
+            double partial = 1.0;
+            for (auto& node : getNodes()) {
+                partial *= node.second->getFactorValue(completeEvidence);
             }
-        for (auto& thread : threads) {
-            thread.join();
+            result[i] = partial;
         }
         // Normalize result
-        double sum = accumulate(result.begin(), result.end(), 0.0);
+        double sum = std::accumulate(result.begin(), result.end(), 0.0);
         transform(result.begin(), result.end(), result.begin(), [sum](const double& value) { return value / sum; });
         return result;
     }

bayesnet/network/Network.h

@@ -12,14 +12,18 @@
 #include "Node.h"
 
 namespace bayesnet {
+    enum class Smoothing_t {
+        NONE = -1,
+        ORIGINAL = 0,
+        LAPLACE,
+        CESTNIK
+    };
     class Network {
     public:
         Network();
-        explicit Network(float);
         explicit Network(const Network&);
         ~Network() = default;
         torch::Tensor& getSamples();
-        float getMaxThreads() const;
         void addNode(const std::string&);
         void addEdge(const std::string&, const std::string&);
         std::map<std::string, std::unique_ptr<Node>>& getNodes();
@@ -32,9 +36,9 @@ namespace bayesnet {
         /*
         Notice: Nodes have to be inserted in the same order as they are in the dataset, i.e., first node is first column and so on.
         */
-        void fit(const std::vector<std::vector<int>>& input_data, const std::vector<int>& labels, const std::vector<double>& weights, const std::vector<std::string>& featureNames, const std::string& className, const std::map<std::string, std::vector<int>>& states);
-        void fit(const torch::Tensor& X, const torch::Tensor& y, const torch::Tensor& weights, const std::vector<std::string>& featureNames, const std::string& className, const std::map<std::string, std::vector<int>>& states);
-        void fit(const torch::Tensor& samples, const torch::Tensor& weights, const std::vector<std::string>& featureNames, const std::string& className, const std::map<std::string, std::vector<int>>& states);
+        void fit(const std::vector<std::vector<int>>& input_data, const std::vector<int>& labels, const std::vector<double>& weights, const std::vector<std::string>& featureNames, const std::string& className, const std::map<std::string, std::vector<int>>& states, const Smoothing_t smoothing);
+        void fit(const torch::Tensor& X, const torch::Tensor& y, const torch::Tensor& weights, const std::vector<std::string>& featureNames, const std::string& className, const std::map<std::string, std::vector<int>>& states, const Smoothing_t smoothing);
+        void fit(const torch::Tensor& samples, const torch::Tensor& weights, const std::vector<std::string>& featureNames, const std::string& className, const std::map<std::string, std::vector<int>>& states, const Smoothing_t smoothing);
         std::vector<int> predict(const std::vector<std::vector<int>>&); // Return mx1 std::vector of predictions
         torch::Tensor predict(const torch::Tensor&); // Return mx1 tensor of predictions
         torch::Tensor predict_tensor(const torch::Tensor& samples, const bool proba);
@@ -50,19 +54,16 @@ namespace bayesnet {
     private:
         std::map<std::string, std::unique_ptr<Node>> nodes;
         bool fitted;
-        float maxThreads = 0.95;
         int classNumStates;
         std::vector<std::string> features; // Including classname
         std::string className;
-        double laplaceSmoothing;
         torch::Tensor samples; // n+1xm tensor used to fit the model
         bool isCyclic(const std::string&, std::unordered_set<std::string>&, std::unordered_set<std::string>&);
         std::vector<double> predict_sample(const std::vector<int>&);
         std::vector<double> predict_sample(const torch::Tensor&);
         std::vector<double> exactInference(std::map<std::string, int>&);
-        double computeFactor(std::map<std::string, int>&);
-        void completeFit(const std::map<std::string, std::vector<int>>& states, const torch::Tensor& weights);
-        void checkFitData(int n_features, int n_samples, int n_samples_y, const std::vector<std::string>& featureNames, const std::string& className, const std::map<std::string, std::vector<int>>& states, const torch::Tensor& weights);
+        void completeFit(const std::map<std::string, std::vector<int>>& states, const torch::Tensor& weights, const Smoothing_t smoothing);
+        void checkFitData(int n_samples, int n_features, int n_samples_y, const std::vector<std::string>& featureNames, const std::string& className, const std::map<std::string, std::vector<int>>& states, const torch::Tensor& weights);
         void setStates(const std::map<std::string, std::vector<int>>&);
     };
 }

bayesnet/network/Node.cc

@@ -90,51 +90,54 @@ namespace bayesnet {
         }
         return result;
     }
-    void Node::computeCPT(const torch::Tensor& dataset, const std::vector<std::string>& features, const double laplaceSmoothing, const torch::Tensor& weights)
+    void Node::computeCPT(const torch::Tensor& dataset, const std::vector<std::string>& features, const double smoothing, const torch::Tensor& weights)
     {
         dimensions.clear();
         // Get dimensions of the CPT
         dimensions.push_back(numStates);
         transform(parents.begin(), parents.end(), back_inserter(dimensions), [](const auto& parent) { return parent->getNumStates(); });
         // Create a tensor of zeros with the dimensions of the CPT
-        cpTable = torch::zeros(dimensions, torch::kFloat) + laplaceSmoothing;
+        cpTable = torch::zeros(dimensions, torch::kDouble) + smoothing;
         // Fill table with counts
         auto pos = find(features.begin(), features.end(), name);
         if (pos == features.end()) {
             throw std::logic_error("Feature " + name + " not found in dataset");
         }
         int name_index = pos - features.begin();
-        for (int n_sample = 0; n_sample < dataset.size(1); ++n_sample) {
         c10::List<c10::optional<at::Tensor>> coordinates;
-            coordinates.push_back(dataset.index({ name_index, n_sample }));
+        for (int n_sample = 0; n_sample < dataset.size(1); ++n_sample) {
+            coordinates.clear();
+            auto sample = dataset.index({ "...", n_sample });
+            coordinates.push_back(sample[name_index]);
             for (auto parent : parents) {
                 pos = find(features.begin(), features.end(), parent->getName());
                 if (pos == features.end()) {
                     throw std::logic_error("Feature parent " + parent->getName() + " not found in dataset");
                 }
                 int parent_index = pos - features.begin();
-                coordinates.push_back(dataset.index({ parent_index, n_sample }));
+                coordinates.push_back(sample[parent_index]);
             }
             // Increment the count of the corresponding coordinate
-            cpTable.index_put_({ coordinates }, cpTable.index({ coordinates }) + weights.index({ n_sample }).item<double>());
+            cpTable.index_put_({ coordinates }, weights.index({ n_sample }), true);
         }
         // Normalize the counts
+        // Divide each row by the sum of the row
         cpTable = cpTable / cpTable.sum(0);
     }
-    float Node::getFactorValue(std::map<std::string, int>& evidence)
+    double Node::getFactorValue(std::map<std::string, int>& evidence)
     {
         c10::List<c10::optional<at::Tensor>> coordinates;
         // following predetermined order of indices in the cpTable (see Node.h)
         coordinates.push_back(at::tensor(evidence[name]));
         transform(parents.begin(), parents.end(), std::back_inserter(coordinates), [&evidence](const auto& parent) { return at::tensor(evidence[parent->getName()]); });
-        return cpTable.index({ coordinates }).item<float>();
+        return cpTable.index({ coordinates }).item<double>();
     }
     std::vector<std::string> Node::graph(const std::string& className)
     {
         auto output = std::vector<std::string>();
         auto suffix = name == className ? ", fontcolor=red, fillcolor=lightblue, style=filled " : "";
-        output.push_back(name + " [shape=circle" + suffix + "] \n");
-        transform(children.begin(), children.end(), back_inserter(output), [this](const auto& child) { return name + " -> " + child->getName(); });
+        output.push_back("\"" + name + "\" [shape=circle" + suffix + "] \n");
+        transform(children.begin(), children.end(), back_inserter(output), [this](const auto& child) { return "\"" + name + "\" -> \"" + child->getName() + "\""; });
         return output;
     }
 }

bayesnet/network/Node.h

@@ -23,12 +23,12 @@ namespace bayesnet {
         std::vector<Node*>& getParents();
         std::vector<Node*>& getChildren();
         torch::Tensor& getCPT();
-        void computeCPT(const torch::Tensor& dataset, const std::vector<std::string>& features, const double laplaceSmoothing, const torch::Tensor& weights);
+        void computeCPT(const torch::Tensor& dataset, const std::vector<std::string>& features, const double smoothing, const torch::Tensor& weights);
         int getNumStates() const;
         void setNumStates(int);
         unsigned minFill();
         std::vector<std::string> graph(const std::string& clasName); // Returns a std::vector of std::strings representing the graph in graphviz format
-        float getFactorValue(std::map<std::string, int>&);
+        double getFactorValue(std::map<std::string, int>&);
     private:
         std::string name;
         std::vector<Node*> parents;

bayesnet/utils/CountingSemaphore.h

@@ -0,0 +1,46 @@
+#ifndef COUNTING_SEMAPHORE_H
+#define COUNTING_SEMAPHORE_H
+#include <mutex>
+#include <condition_variable>
+#include <algorithm>
+#include <thread>
+#include <mutex>
+#include <condition_variable>
+#include <thread>
+
+class CountingSemaphore {
+public:
+    static CountingSemaphore& getInstance()
+    {
+        static CountingSemaphore instance;
+        return instance;
+    }
+    // Delete copy constructor and assignment operator
+    CountingSemaphore(const CountingSemaphore&) = delete;
+    CountingSemaphore& operator=(const CountingSemaphore&) = delete;
+    void acquire()
+    {
+        std::unique_lock<std::mutex> lock(mtx_);
+        cv_.wait(lock, [this]() { return count_ > 0; });
+        --count_;
+    }
+    void release()
+    {
+        std::lock_guard<std::mutex> lock(mtx_);
+        ++count_;
+        if (count_ <= max_count_) {
+            cv_.notify_one();
+        }
+    }
+private:
+    CountingSemaphore()
+        : max_count_(std::max(1u, static_cast<uint>(0.95 * std::thread::hardware_concurrency()))),
+        count_(max_count_)
+    {
+    }
+    std::mutex mtx_;
+    std::condition_variable cv_;
+    const uint max_count_;
+    uint count_;
+};
+#endif

cmake/modules/CodeCoverage.cmake

@@ -137,7 +137,7 @@
 
 include(CMakeParseArguments)
 
-option(CODE_COVERAGE_VERBOSE "Verbose information" FALSE)
+option(CODE_COVERAGE_VERBOSE "Verbose information" TRUE)
 
 # Check prereqs
 find_program( GCOV_PATH gcov )
@@ -160,8 +160,12 @@ foreach(LANG ${LANGUAGES})
     endif()
   elseif(NOT "${CMAKE_${LANG}_COMPILER_ID}" MATCHES "GNU"
          AND NOT "${CMAKE_${LANG}_COMPILER_ID}" MATCHES "(LLVM)?[Ff]lang")
+         if ("${LANG}" MATCHES "CUDA")
+             message(STATUS "Ignoring CUDA")
+        else()
             message(FATAL_ERROR "Compiler is not GNU or Flang! Aborting...")
         endif()
+  endif()
 endforeach()
 
 set(COVERAGE_COMPILER_FLAGS "-g --coverage"

sample/CMakeLists.txt

@@ -5,15 +5,22 @@ project(bayesnet_sample)
 set(CMAKE_CXX_STANDARD 17)
 
 find_package(Torch REQUIRED)
-find_library(BayesNet NAMES BayesNet.a libBayesNet.a REQUIRED)
+find_library(BayesNet NAMES libBayesNet BayesNet libBayesNet.a  REQUIRED)
+find_path(Bayesnet_INCLUDE_DIRS REQUIRED NAMES bayesnet)
+find_library(FImdlp NAMES libfimdlp.a  PATHS REQUIRED)
+find_path(FImdlp_INCLUDE_DIRS REQUIRED NAMES fimdlp)
+
+message(STATUS "FImdlp=${FImdlp}")
+message(STATUS "FImdlp_INCLUDE_DIRS=${FImdlp_INCLUDE_DIRS}")
+message(STATUS "BayesNet=${BayesNet}")
+message(STATUS "Bayesnet_INCLUDE_DIRS=${Bayesnet_INCLUDE_DIRS}")
 
 include_directories(
-    lib/Files
-    lib/mdlp
+    ../tests/lib/Files
     lib/json/include
     /usr/local/include
+    ${FImdlp_INCLUDE_DIRS}
 )
 
-add_subdirectory(lib/mdlp)
 add_executable(bayesnet_sample sample.cc) 
-target_link_libraries(bayesnet_sample mdlp "${TORCH_LIBRARIES}" "${BayesNet}")
+target_link_libraries(bayesnet_sample ${FImdlp} "${TORCH_LIBRARIES}" "${BayesNet}")

sample/lib/mdlp/CMakeLists.txt

@@ -1,11 +0,0 @@
-cmake_minimum_required(VERSION 3.20)
-project(mdlp)
-
-if (POLICY CMP0135)
-    cmake_policy(SET CMP0135 NEW)
-endif ()
-
-set(CMAKE_CXX_STANDARD 11)
-
-add_library(mdlp CPPFImdlp.cpp Metrics.cpp)
-

sample/lib/mdlp/CPPFImdlp.cpp

@@ -1,222 +0,0 @@
-#include <numeric>
-#include <algorithm>
-#include <set>
-#include <cmath>
-#include "CPPFImdlp.h"
-#include "Metrics.h"
-
-namespace mdlp {
-
-    CPPFImdlp::CPPFImdlp(size_t min_length_, int max_depth_, float proposed) : min_length(min_length_),
-        max_depth(max_depth_),
-        proposed_cuts(proposed)
-    {
-    }
-
-    CPPFImdlp::CPPFImdlp() = default;
-
-    CPPFImdlp::~CPPFImdlp() = default;
-
-    size_t CPPFImdlp::compute_max_num_cut_points() const
-    {
-        // Set the actual maximum number of cut points as a number or as a percentage of the number of samples
-        if (proposed_cuts == 0) {
-            return numeric_limits<size_t>::max();
-        }
-        if (proposed_cuts < 0 || proposed_cuts > static_cast<float>(X.size())) {
-            throw invalid_argument("wrong proposed num_cuts value");
-        }
-        if (proposed_cuts < 1)
-            return static_cast<size_t>(round(static_cast<float>(X.size()) * proposed_cuts));
-        return static_cast<size_t>(proposed_cuts);
-    }
-
-    void CPPFImdlp::fit(samples_t& X_, labels_t& y_)
-    {
-        X = X_;
-        y = y_;
-        num_cut_points = compute_max_num_cut_points();
-        depth = 0;
-        discretizedData.clear();
-        cutPoints.clear();
-        if (X.size() != y.size()) {
-            throw invalid_argument("X and y must have the same size");
-        }
-        if (X.empty() || y.empty()) {
-            throw invalid_argument("X and y must have at least one element");
-        }
-        if (min_length < 3) {
-            throw invalid_argument("min_length must be greater than 2");
-        }
-        if (max_depth < 1) {
-            throw invalid_argument("max_depth must be greater than 0");
-        }
-        indices = sortIndices(X_, y_);
-        metrics.setData(y, indices);
-        computeCutPoints(0, X.size(), 1);
-        sort(cutPoints.begin(), cutPoints.end());
-        if (num_cut_points > 0) {
-            // Select the best (with lower entropy) cut points
-            while (cutPoints.size() > num_cut_points) {
-                resizeCutPoints();
-            }
-        }
-    }
-
-    pair<precision_t, size_t> CPPFImdlp::valueCutPoint(size_t start, size_t cut, size_t end)
-    {
-        size_t n;
-        size_t m;
-        size_t idxPrev = cut - 1 >= start ? cut - 1 : cut;
-        size_t idxNext = cut + 1 < end ? cut + 1 : cut;
-        bool backWall; // true if duplicates reach beginning of the interval
-        precision_t previous;
-        precision_t actual;
-        precision_t next;
-        previous = X[indices[idxPrev]];
-        actual = X[indices[cut]];
-        next = X[indices[idxNext]];
-        // definition 2 of the paper => X[t-1] < X[t]
-        // get the first equal value of X in the interval
-        while (idxPrev > start && actual == previous) {
-            previous = X[indices[--idxPrev]];
-        }
-        backWall = idxPrev == start && actual == previous;
-        // get the last equal value of X in the interval
-        while (idxNext < end - 1 && actual == next) {
-            next = X[indices[++idxNext]];
-        }
-        // # of duplicates before cutpoint
-        n = cut - 1 - idxPrev;
-        // # of duplicates after cutpoint
-        m = idxNext - cut - 1;
-        // Decide which values to use
-        cut = cut + (backWall ? m + 1 : -n);
-        actual = X[indices[cut]];
-        return { (actual + previous) / 2, cut };
-    }
-
-    void CPPFImdlp::computeCutPoints(size_t start, size_t end, int depth_)
-    {
-        size_t cut;
-        pair<precision_t, size_t> result;
-        // Check if the interval length and the depth are Ok
-        if (end - start < min_length || depth_ > max_depth)
-            return;
-        depth = depth_ > depth ? depth_ : depth;
-        cut = getCandidate(start, end);
-        if (cut == numeric_limits<size_t>::max())
-            return;
-        if (mdlp(start, cut, end)) {
-            result = valueCutPoint(start, cut, end);
-            cut = result.second;
-            cutPoints.push_back(result.first);
-            computeCutPoints(start, cut, depth_ + 1);
-            computeCutPoints(cut, end, depth_ + 1);
-        }
-    }
-
-    size_t CPPFImdlp::getCandidate(size_t start, size_t end)
-    {
-        /* Definition 1: A binary discretization for A is determined by selecting the cut point TA for which
-        E(A, TA; S) is minimal amongst all the candidate cut points. */
-        size_t candidate = numeric_limits<size_t>::max();
-        size_t elements = end - start;
-        bool sameValues = true;
-        precision_t entropy_left;
-        precision_t entropy_right;
-        precision_t minEntropy;
-        // Check if all the values of the variable in the interval are the same
-        for (size_t idx = start + 1; idx < end; idx++) {
-            if (X[indices[idx]] != X[indices[start]]) {
-                sameValues = false;
-                break;
-            }
-        }
-        if (sameValues)
-            return candidate;
-        minEntropy = metrics.entropy(start, end);
-        for (size_t idx = start + 1; idx < end; idx++) {
-            // Cutpoints are always on boundaries (definition 2)
-            if (y[indices[idx]] == y[indices[idx - 1]])
-                continue;
-            entropy_left = precision_t(idx - start) / static_cast<precision_t>(elements) * metrics.entropy(start, idx);
-            entropy_right = precision_t(end - idx) / static_cast<precision_t>(elements) * metrics.entropy(idx, end);
-            if (entropy_left + entropy_right < minEntropy) {
-                minEntropy = entropy_left + entropy_right;
-                candidate = idx;
-            }
-        }
-        return candidate;
-    }
-
-    bool CPPFImdlp::mdlp(size_t start, size_t cut, size_t end)
-    {
-        int k;
-        int k1;
-        int k2;
-        precision_t ig;
-        precision_t delta;
-        precision_t ent;
-        precision_t ent1;
-        precision_t ent2;
-        auto N = precision_t(end - start);
-        k = metrics.computeNumClasses(start, end);
-        k1 = metrics.computeNumClasses(start, cut);
-        k2 = metrics.computeNumClasses(cut, end);
-        ent = metrics.entropy(start, end);
-        ent1 = metrics.entropy(start, cut);
-        ent2 = metrics.entropy(cut, end);
-        ig = metrics.informationGain(start, cut, end);
-        delta = static_cast<precision_t>(log2(pow(3, precision_t(k)) - 2) -
-            (precision_t(k) * ent - precision_t(k1) * ent1 - precision_t(k2) * ent2));
-        precision_t term = 1 / N * (log2(N - 1) + delta);
-        return ig > term;
-    }
-
-    // Argsort from https://stackoverflow.com/questions/1577475/c-sorting-and-keeping-track-of-indexes
-    indices_t CPPFImdlp::sortIndices(samples_t& X_, labels_t& y_)
-    {
-        indices_t idx(X_.size());
-        iota(idx.begin(), idx.end(), 0);
-        stable_sort(idx.begin(), idx.end(), [&X_, &y_](size_t i1, size_t i2) {
-            if (X_[i1] == X_[i2])
-                return y_[i1] < y_[i2];
-            else
-                return X_[i1] < X_[i2];
-            });
-        return idx;
-    }
-
-    void CPPFImdlp::resizeCutPoints()
-    {
-        //Compute entropy of each of the whole cutpoint set and discards the biggest value
-        precision_t maxEntropy = 0;
-        precision_t entropy;
-        size_t maxEntropyIdx = 0;
-        size_t begin = 0;
-        size_t end;
-        for (size_t idx = 0; idx < cutPoints.size(); idx++) {
-            end = begin;
-            while (X[indices[end]] < cutPoints[idx] && end < X.size())
-                end++;
-            entropy = metrics.entropy(begin, end);
-            if (entropy > maxEntropy) {
-                maxEntropy = entropy;
-                maxEntropyIdx = idx;
-            }
-            begin = end;
-        }
-        cutPoints.erase(cutPoints.begin() + static_cast<long>(maxEntropyIdx));
-    }
-    labels_t& CPPFImdlp::transform(const samples_t& data)
-    {
-        discretizedData.clear();
-        discretizedData.reserve(data.size());
-        for (const precision_t& item : data) {
-            auto upper = upper_bound(cutPoints.begin(), cutPoints.end(), item);
-            discretizedData.push_back(upper - cutPoints.begin());
-        }
-        return discretizedData;
-    }
-}

sample/lib/mdlp/CPPFImdlp.h

@@ -1,51 +0,0 @@
-// ***************************************************************
-// SPDX-FileCopyrightText: Copyright 2024 Ricardo Montañana Gómez
-// SPDX-FileType: SOURCE
-// SPDX-License-Identifier: MIT
-// ***************************************************************
-
-#ifndef CPPFIMDLP_H
-#define CPPFIMDLP_H
-
-#include "typesFImdlp.h"
-#include "Metrics.h"
-#include <limits>
-#include <utility>
-#include <string>
-
-namespace mdlp {
-    class CPPFImdlp {
-    protected:
-        size_t min_length = 3;
-        int depth = 0;
-        int max_depth = numeric_limits<int>::max();
-        float proposed_cuts = 0;
-        indices_t indices = indices_t();
-        samples_t X = samples_t();
-        labels_t y = labels_t();
-        Metrics metrics = Metrics(y, indices);
-        cutPoints_t cutPoints;
-        size_t num_cut_points = numeric_limits<size_t>::max();
-        labels_t discretizedData = labels_t();
-
-        static indices_t sortIndices(samples_t&, labels_t&);
-
-        void computeCutPoints(size_t, size_t, int);
-        void resizeCutPoints();
-        bool mdlp(size_t, size_t, size_t);
-        size_t getCandidate(size_t, size_t);
-        size_t compute_max_num_cut_points() const;
-        pair<precision_t, size_t> valueCutPoint(size_t, size_t, size_t);
-
-    public:
-        CPPFImdlp();
-        CPPFImdlp(size_t, int, float);
-        ~CPPFImdlp();
-        void fit(samples_t&, labels_t&);
-        inline cutPoints_t getCutPoints() const { return cutPoints; };
-        labels_t& transform(const samples_t&);
-        inline int get_depth() const { return depth; };
-        static inline string version() { return "1.1.2"; };
-    };
-}
-#endif

sample/lib/mdlp/LICENSE

@@ -1,21 +0,0 @@
-MIT License
-
-Copyright (c) 2022 Ricardo Montañana Gómez
-
-Permission is hereby granted, free of charge, to any person obtaining a copy
-of this software and associated documentation files (the "Software"), to deal
-in the Software without restriction, including without limitation the rights
-to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-copies of the Software, and to permit persons to whom the Software is
-furnished to do so, subject to the following conditions:
-
-The above copyright notice and this permission notice shall be included in all
-copies or substantial portions of the Software.
-
-THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
-SOFTWARE.

sample/lib/mdlp/Metrics.cpp

@@ -1,78 +0,0 @@
-#include "Metrics.h"
-#include <set>
-#include <cmath>
-
-using namespace std;
-namespace mdlp {
-    Metrics::Metrics(labels_t& y_, indices_t& indices_): y(y_), indices(indices_),
-        numClasses(computeNumClasses(0, indices.size()))
-    {
-    }
-
-    int Metrics::computeNumClasses(size_t start, size_t end)
-    {
-        set<int> nClasses;
-        for (auto i = start; i < end; ++i) {
-            nClasses.insert(y[indices[i]]);
-        }
-        return static_cast<int>(nClasses.size());
-    }
-
-    void Metrics::setData(const labels_t& y_, const indices_t& indices_)
-    {
-        indices = indices_;
-        y = y_;
-        numClasses = computeNumClasses(0, indices.size());
-        entropyCache.clear();
-        igCache.clear();
-    }
-
-    precision_t Metrics::entropy(size_t start, size_t end)
-    {
-        precision_t p;
-        precision_t ventropy = 0;
-        int nElements = 0;
-        labels_t counts(numClasses + 1, 0);
-        if (end - start < 2)
-            return 0;
-        if (entropyCache.find({ start, end }) != entropyCache.end()) {
-            return entropyCache[{start, end}];
-        }
-        for (auto i = &indices[start]; i != &indices[end]; ++i) {
-            counts[y[*i]]++;
-            nElements++;
-        }
-        for (auto count : counts) {
-            if (count > 0) {
-                p = static_cast<precision_t>(count) / static_cast<precision_t>(nElements);
-                ventropy -= p * log2(p);
-            }
-        }
-        entropyCache[{start, end}] = ventropy;
-        return ventropy;
-    }
-
-    precision_t Metrics::informationGain(size_t start, size_t cut, size_t end)
-    {
-        precision_t iGain;
-        precision_t entropyInterval;
-        precision_t entropyLeft;
-        precision_t entropyRight;
-        size_t nElementsLeft = cut - start;
-        size_t nElementsRight = end - cut;
-        size_t nElements = end - start;
-        if (igCache.find(make_tuple(start, cut, end)) != igCache.end()) {
-            return igCache[make_tuple(start, cut, end)];
-        }
-        entropyInterval = entropy(start, end);
-        entropyLeft = entropy(start, cut);
-        entropyRight = entropy(cut, end);
-        iGain = entropyInterval -
-            (static_cast<precision_t>(nElementsLeft) * entropyLeft +
-                static_cast<precision_t>(nElementsRight) * entropyRight) /
-            static_cast<precision_t>(nElements);
-        igCache[make_tuple(start, cut, end)] = iGain;
-        return iGain;
-    }
-
-}

sample/lib/mdlp/Metrics.h

@@ -1,28 +0,0 @@
-// ***************************************************************
-// SPDX-FileCopyrightText: Copyright 2024 Ricardo Montañana Gómez
-// SPDX-FileType: SOURCE
-// SPDX-License-Identifier: MIT
-// ***************************************************************
-
-#ifndef CCMETRICS_H
-#define CCMETRICS_H
-
-#include "typesFImdlp.h"
-
-namespace mdlp {
-    class Metrics {
-    protected:
-        labels_t& y;
-        indices_t& indices;
-        int numClasses;
-        cacheEnt_t entropyCache = cacheEnt_t();
-        cacheIg_t igCache = cacheIg_t();
-    public:
-        Metrics(labels_t&, indices_t&);
-        void setData(const labels_t&, const indices_t&);
-        int computeNumClasses(size_t, size_t);
-        precision_t entropy(size_t, size_t);
-        precision_t informationGain(size_t, size_t, size_t);
-    };
-}
-#endif

sample/lib/mdlp/README.md

@@ -1,41 +0,0 @@
-[![Build](https://github.com/rmontanana/mdlp/actions/workflows/build.yml/badge.svg)](https://github.com/rmontanana/mdlp/actions/workflows/build.yml)
-[![Quality Gate Status](https://sonarcloud.io/api/project_badges/measure?project=rmontanana_mdlp&metric=alert_status)](https://sonarcloud.io/summary/new_code?id=rmontanana_mdlp)
-[![Reliability Rating](https://sonarcloud.io/api/project_badges/measure?project=rmontanana_mdlp&metric=reliability_rating)](https://sonarcloud.io/summary/new_code?id=rmontanana_mdlp)
-
-# mdlp
-
-Discretization algorithm based on the paper by Fayyad & Irani [Multi-Interval Discretization of Continuous-Valued Attributes for Classification Learning](https://www.ijcai.org/Proceedings/93-2/Papers/022.pdf)
-
-The implementation tries to mitigate the problem of different label values with the same value of the variable:
-
-- Sorts the values of the variable using the label values as a tie-breaker
-- Once found a valid candidate for the split, it checks if the previous value is the same as actual one, and tries to get previous one, or next if the former is not possible.
-
-Other features:
-
-- Intervals with the same value of the variable are not taken into account for cutpoints.
-- Intervals have to have more than two examples to be evaluated.
-
-The algorithm returns the cut points for the variable.
-
-## Sample
-
-To run the sample, just execute the following commands:
-
-```bash
-cd sample
-cmake -B build
-cd build
-make
-./sample -f iris -m 2
-./sample -h
-```
-
-## Test
-
-To run the tests and see coverage (llvm & gcovr have to be installed), execute the following commands:
-
-```bash
-cd tests
-./test
-```

sample/lib/mdlp/typesFImdlp.h

@@ -1,24 +0,0 @@
-// ***************************************************************
-// SPDX-FileCopyrightText: Copyright 2024 Ricardo Montañana Gómez
-// SPDX-FileType: SOURCE
-// SPDX-License-Identifier: MIT
-// ***************************************************************
-
-#ifndef TYPES_H
-#define TYPES_H
-
-#include <vector>
-#include <map>
-#include <stdexcept>
-
-using namespace std;
-namespace mdlp {
-    typedef float precision_t;
-    typedef vector<precision_t> samples_t;
-    typedef vector<int> labels_t;
-    typedef vector<size_t> indices_t;
-    typedef vector<precision_t> cutPoints_t;
-    typedef map<pair<int, int>, precision_t> cacheEnt_t;
-    typedef map<tuple<int, int, int>, precision_t> cacheIg_t;
-}
-#endif

sample/sample.cc

@@ -5,7 +5,7 @@
 // ***************************************************************
 
 #include <ArffFiles.hpp>
-#include <CPPFImdlp.h>
+#include <fimdlp/CPPFImdlp.h>
 #include <bayesnet/ensembles/BoostAODE.h>
 
 std::vector<mdlp::labels_t> discretizeDataset(std::vector<mdlp::samples_t>& X, mdlp::labels_t& y)
@@ -60,9 +60,9 @@ int main(int argc, char* argv[])
     auto clf = bayesnet::BoostAODE(false); // false for not using voting in predict
     std::cout << "Library version: " << clf.getVersion() << std::endl;
     tie(X, y, features, className, states) = loadDataset(file_name, true);
-    clf.fit(X, y, features, className, states);
+    clf.fit(X, y, features, className, states, bayesnet::Smoothing_t::LAPLACE);
     auto score = clf.score(X, y);
-    std::cout << "File: " << file_name << " score: " << score << std::endl;
+    std::cout << "File: " << file_name << " Model: BoostAODE score: " << score << std::endl;
     return 0;
 }
 

tests/CMakeLists.txt

@@ -2,10 +2,11 @@ if(ENABLE_TESTING)
     include_directories(
         ${BayesNet_SOURCE_DIR}/tests/lib/Files
         ${BayesNet_SOURCE_DIR}/lib/folding
-        ${BayesNet_SOURCE_DIR}/lib/mdlp
+        ${BayesNet_SOURCE_DIR}/lib/mdlp/src
         ${BayesNet_SOURCE_DIR}/lib/json/include
         ${BayesNet_SOURCE_DIR}
         ${CMAKE_BINARY_DIR}/configured_files/include
+        ${FImdlp_INCLUDE_DIRS}
     )
     file(GLOB_RECURSE BayesNet_SOURCES "${BayesNet_SOURCE_DIR}/bayesnet/*.cc")
     add_executable(TestBayesNet TestBayesNetwork.cc TestBayesNode.cc TestBayesClassifier.cc 

tests/TestA2DE.cc

@@ -16,7 +16,7 @@ TEST_CASE("Fit and Score", "[A2DE]")
 {
     auto raw = RawDatasets("glass", true);
     auto clf = bayesnet::A2DE();
-    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states);
+    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states, raw.smoothing);
     REQUIRE(clf.score(raw.Xv, raw.yv) == Catch::Approx(0.831776).epsilon(raw.epsilon));
     REQUIRE(clf.getNumberOfNodes() == 360);
     REQUIRE(clf.getNumberOfEdges() == 756);
@@ -30,18 +30,18 @@ TEST_CASE("Test score with predict_voting", "[A2DE]")
        {"predict_voting", true},
     };
     clf.setHyperparameters(hyperparameters);
-    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states);
+    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states, raw.smoothing);
     REQUIRE(clf.score(raw.Xv, raw.yv) == Catch::Approx(0.82243).epsilon(raw.epsilon));
     hyperparameters["predict_voting"] = false;
     clf.setHyperparameters(hyperparameters);
-    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states);
+    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states, raw.smoothing);
     REQUIRE(clf.score(raw.Xv, raw.yv) == Catch::Approx(0.83178).epsilon(raw.epsilon));
 }
 TEST_CASE("Test graph", "[A2DE]")
 {
     auto raw = RawDatasets("iris", true);
     auto clf = bayesnet::A2DE();
-    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states);
+    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states, raw.smoothing);
     auto graph = clf.graph();
     REQUIRE(graph.size() == 78);
     REQUIRE(graph[0] == "digraph BayesNet {\nlabel=<BayesNet A2DE_0>\nfontsize=30\nfontcolor=blue\nlabelloc=t\nlayout=circo\n");

tests/TestBayesClassifier.cc

@@ -18,38 +18,38 @@ TEST_CASE("Test Cannot build dataset with wrong data vector", "[Classifier]")
     auto model = bayesnet::TAN();
     auto raw = RawDatasets("iris", true);
     raw.yv.pop_back();
-    REQUIRE_THROWS_AS(model.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states), std::runtime_error);
-    REQUIRE_THROWS_WITH(model.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states), "* Error in X and y dimensions *\nX dimensions: [4, 150]\ny dimensions: [149]");
+    REQUIRE_THROWS_AS(model.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states, raw.smoothing), std::runtime_error);
+    REQUIRE_THROWS_WITH(model.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states, raw.smoothing), "* Error in X and y dimensions *\nX dimensions: [4, 150]\ny dimensions: [149]");
 }
 TEST_CASE("Test Cannot build dataset with wrong data tensor", "[Classifier]")
 {
     auto model = bayesnet::TAN();
     auto raw = RawDatasets("iris", true);
     auto yshort = torch::zeros({ 149 }, torch::kInt32);
-    REQUIRE_THROWS_AS(model.fit(raw.Xt, yshort, raw.features, raw.className, raw.states), std::runtime_error);
-    REQUIRE_THROWS_WITH(model.fit(raw.Xt, yshort, raw.features, raw.className, raw.states), "* Error in X and y dimensions *\nX dimensions: [4, 150]\ny dimensions: [149]");
+    REQUIRE_THROWS_AS(model.fit(raw.Xt, yshort, raw.features, raw.className, raw.states, raw.smoothing), std::runtime_error);
+    REQUIRE_THROWS_WITH(model.fit(raw.Xt, yshort, raw.features, raw.className, raw.states, raw.smoothing), "* Error in X and y dimensions *\nX dimensions: [4, 150]\ny dimensions: [149]");
 }
 TEST_CASE("Invalid data type", "[Classifier]")
 {
     auto model = bayesnet::TAN();
     auto raw = RawDatasets("iris", false);
-    REQUIRE_THROWS_AS(model.fit(raw.Xt, raw.yt, raw.features, raw.className, raw.states), std::invalid_argument);
-    REQUIRE_THROWS_WITH(model.fit(raw.Xt, raw.yt, raw.features, raw.className, raw.states), "dataset (X, y) must be of type Integer");
+    REQUIRE_THROWS_AS(model.fit(raw.Xt, raw.yt, raw.features, raw.className, raw.states, raw.smoothing), std::invalid_argument);
+    REQUIRE_THROWS_WITH(model.fit(raw.Xt, raw.yt, raw.features, raw.className, raw.states, raw.smoothing), "dataset (X, y) must be of type Integer");
 }
 TEST_CASE("Invalid number of features", "[Classifier]")
 {
     auto model = bayesnet::TAN();
     auto raw = RawDatasets("iris", true);
     auto Xt = torch::cat({ raw.Xt, torch::zeros({ 1, 150 }, torch::kInt32) }, 0);
-    REQUIRE_THROWS_AS(model.fit(Xt, raw.yt, raw.features, raw.className, raw.states), std::invalid_argument);
-    REQUIRE_THROWS_WITH(model.fit(Xt, raw.yt, raw.features, raw.className, raw.states), "Classifier: X 5 and features 4 must have the same number of features");
+    REQUIRE_THROWS_AS(model.fit(Xt, raw.yt, raw.features, raw.className, raw.states, raw.smoothing), std::invalid_argument);
+    REQUIRE_THROWS_WITH(model.fit(Xt, raw.yt, raw.features, raw.className, raw.states, raw.smoothing), "Classifier: X 5 and features 4 must have the same number of features");
 }
 TEST_CASE("Invalid class name", "[Classifier]")
 {
     auto model = bayesnet::TAN();
     auto raw = RawDatasets("iris", true);
-    REQUIRE_THROWS_AS(model.fit(raw.Xt, raw.yt, raw.features, "duck", raw.states), std::invalid_argument);
-    REQUIRE_THROWS_WITH(model.fit(raw.Xt, raw.yt, raw.features, "duck", raw.states), "class name not found in states");
+    REQUIRE_THROWS_AS(model.fit(raw.Xt, raw.yt, raw.features, "duck", raw.states, raw.smoothing), std::invalid_argument);
+    REQUIRE_THROWS_WITH(model.fit(raw.Xt, raw.yt, raw.features, "duck", raw.states, raw.smoothing), "class name not found in states");
 }
 TEST_CASE("Invalid feature name", "[Classifier]")
 {
@@ -57,8 +57,8 @@ TEST_CASE("Invalid feature name", "[Classifier]")
     auto raw = RawDatasets("iris", true);
     auto statest = raw.states;
     statest.erase("petallength");
-    REQUIRE_THROWS_AS(model.fit(raw.Xt, raw.yt, raw.features, raw.className, statest), std::invalid_argument);
-    REQUIRE_THROWS_WITH(model.fit(raw.Xt, raw.yt, raw.features, raw.className, statest), "feature [petallength] not found in states");
+    REQUIRE_THROWS_AS(model.fit(raw.Xt, raw.yt, raw.features, raw.className, statest, raw.smoothing), std::invalid_argument);
+    REQUIRE_THROWS_WITH(model.fit(raw.Xt, raw.yt, raw.features, raw.className, statest, raw.smoothing), "feature [petallength] not found in states");
 }
 TEST_CASE("Invalid hyperparameter", "[Classifier]")
 {
@@ -71,7 +71,7 @@ TEST_CASE("Topological order", "[Classifier]")
 {
     auto model = bayesnet::TAN();
     auto raw = RawDatasets("iris", true);
-    model.fit(raw.Xt, raw.yt, raw.features, raw.className, raw.states);
+    model.fit(raw.Xt, raw.yt, raw.features, raw.className, raw.states, raw.smoothing);
     auto order = model.topological_order();
     REQUIRE(order.size() == 4);
     REQUIRE(order[0] == "petallength");
@@ -83,7 +83,7 @@ TEST_CASE("Dump_cpt", "[Classifier]")
 {
     auto model = bayesnet::TAN();
     auto raw = RawDatasets("iris", true);
-    model.fit(raw.Xt, raw.yt, raw.features, raw.className, raw.states);
+    model.fit(raw.Xt, raw.yt, raw.features, raw.className, raw.states, raw.smoothing);
     auto cpt = model.dump_cpt();
     REQUIRE(cpt.size() == 1713);
 }
@@ -111,7 +111,7 @@ TEST_CASE("KDB Graph", "[Classifier]")
 {
     auto model = bayesnet::KDB(2);
     auto raw = RawDatasets("iris", true);
-    model.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states);
+    model.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states, raw.smoothing);
     auto graph = model.graph();
     REQUIRE(graph.size() == 15);
 }
@@ -119,7 +119,7 @@ TEST_CASE("KDBLd Graph", "[Classifier]")
 {
     auto model = bayesnet::KDBLd(2);
     auto raw = RawDatasets("iris", false);
-    model.fit(raw.Xt, raw.yt, raw.features, raw.className, raw.states);
+    model.fit(raw.Xt, raw.yt, raw.features, raw.className, raw.states, raw.smoothing);
     auto graph = model.graph();
     REQUIRE(graph.size() == 15);
 }

tests/TestBayesEnsemble.cc

@@ -18,7 +18,7 @@ TEST_CASE("Topological Order", "[Ensemble]")
 {
     auto raw = RawDatasets("glass", true);
     auto clf = bayesnet::BoostAODE();
-    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states);
+    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states, raw.smoothing);
     auto order = clf.topological_order();
     REQUIRE(order.size() == 0);
 }
@@ -26,7 +26,7 @@ TEST_CASE("Dump CPT", "[Ensemble]")
 {
     auto raw = RawDatasets("glass", true);
     auto clf = bayesnet::BoostAODE();
-    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states);
+    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states, raw.smoothing);
     auto dump = clf.dump_cpt();
     REQUIRE(dump == "");
 }
@@ -34,7 +34,7 @@ TEST_CASE("Number of States", "[Ensemble]")
 {
     auto clf = bayesnet::BoostAODE();
     auto raw = RawDatasets("iris", true);
-    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states);
+    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states, raw.smoothing);
     REQUIRE(clf.getNumberOfStates() == 76);
 }
 TEST_CASE("Show", "[Ensemble]")
@@ -46,7 +46,7 @@ TEST_CASE("Show", "[Ensemble]")
             {"maxTolerance", 1},
             {"convergence", false},
         });
-    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states);
+    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states, raw.smoothing);
     std::vector<std::string> expected = {
         "class -> sepallength, sepalwidth, petallength, petalwidth, ",
         "petallength -> sepallength, sepalwidth, petalwidth, ",
@@ -78,16 +78,16 @@ TEST_CASE("Graph", "[Ensemble]")
 {
     auto clf = bayesnet::BoostAODE();
     auto raw = RawDatasets("iris", true);
-    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states);
+    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states, raw.smoothing);
     auto graph = clf.graph();
     REQUIRE(graph.size() == 56);
     auto clf2 = bayesnet::AODE();
-    clf2.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states);
+    clf2.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states, raw.smoothing);
     graph = clf2.graph();
     REQUIRE(graph.size() == 56);
     raw = RawDatasets("glass", false);
     auto clf3 = bayesnet::AODELd();
-    clf3.fit(raw.Xt, raw.yt, raw.features, raw.className, raw.states);
+    clf3.fit(raw.Xt, raw.yt, raw.features, raw.className, raw.states, raw.smoothing);
     graph = clf3.graph();
     REQUIRE(graph.size() == 261);
 }

tests/TestBayesModels.cc

@@ -20,7 +20,7 @@
 #include "bayesnet/ensembles/BoostAODE.h"
 #include "TestUtils.h"
 
-const std::string ACTUAL_VERSION = "1.0.5.1";
+const std::string ACTUAL_VERSION = "1.0.6";
 
 TEST_CASE("Test Bayesian Classifiers score & version", "[Models]")
 {
@@ -54,7 +54,7 @@ TEST_CASE("Test Bayesian Classifiers score & version", "[Models]")
             auto clf = models[name];
             auto discretize = name.substr(name.length() - 2) != "Ld";
             auto raw = RawDatasets(file_name, discretize);
-            clf->fit(raw.Xt, raw.yt, raw.features, raw.className, raw.states);
+            clf->fit(raw.Xt, raw.yt, raw.features, raw.className, raw.states, raw.smoothing);
             auto score = clf->score(raw.Xt, raw.yt);
             INFO("Classifier: " << name << " File: " << file_name);
             REQUIRE(score == Catch::Approx(scores[{file_name, name}]).epsilon(raw.epsilon));
@@ -81,7 +81,7 @@ TEST_CASE("Models features & Graph", "[Models]")
     {
         auto raw = RawDatasets("iris", true);
         auto clf = bayesnet::TAN();
-        clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states);
+        clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states, raw.smoothing);
         REQUIRE(clf.getNumberOfNodes() == 5);
         REQUIRE(clf.getNumberOfEdges() == 7);
         REQUIRE(clf.getNumberOfStates() == 19);
@@ -93,7 +93,7 @@ TEST_CASE("Models features & Graph", "[Models]")
     {
         auto clf = bayesnet::TANLd();
         auto raw = RawDatasets("iris", false);
-        clf.fit(raw.Xt, raw.yt, raw.features, raw.className, raw.states);
+        clf.fit(raw.Xt, raw.yt, raw.features, raw.className, raw.states, raw.smoothing);
         REQUIRE(clf.getNumberOfNodes() == 5);
         REQUIRE(clf.getNumberOfEdges() == 7);
         REQUIRE(clf.getNumberOfStates() == 19);
@@ -106,7 +106,7 @@ TEST_CASE("Get num features & num edges", "[Models]")
 {
     auto raw = RawDatasets("iris", true);
     auto clf = bayesnet::KDB(2);
-    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states);
+    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states, raw.smoothing);
     REQUIRE(clf.getNumberOfNodes() == 5);
     REQUIRE(clf.getNumberOfEdges() == 8);
 }
@@ -166,7 +166,7 @@ TEST_CASE("Model predict_proba", "[Models]")
     SECTION("Test " + model + " predict_proba")
     {
         auto clf = models[model];
-        clf->fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states);
+        clf->fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states, raw.smoothing);
         auto y_pred_proba = clf->predict_proba(raw.Xv);
         auto yt_pred_proba = clf->predict_proba(raw.Xt);
         auto y_pred = clf->predict(raw.Xv);
@@ -203,7 +203,7 @@ TEST_CASE("AODE voting-proba", "[Models]")
 {
     auto raw = RawDatasets("glass", true);
     auto clf = bayesnet::AODE(false);
-    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states);
+    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states, raw.smoothing);
     auto score_proba = clf.score(raw.Xv, raw.yv);
     auto pred_proba = clf.predict_proba(raw.Xv);
     clf.setHyperparameters({
@@ -222,9 +222,9 @@ TEST_CASE("SPODELd dataset", "[Models]")
     auto raw = RawDatasets("iris", false);
     auto clf = bayesnet::SPODELd(0);
     // raw.dataset.to(torch::kFloat32);
-    clf.fit(raw.dataset, raw.features, raw.className, raw.states);
+    clf.fit(raw.dataset, raw.features, raw.className, raw.states, raw.smoothing);
     auto score = clf.score(raw.Xt, raw.yt);
-    clf.fit(raw.Xt, raw.yt, raw.features, raw.className, raw.states);
+    clf.fit(raw.Xt, raw.yt, raw.features, raw.className, raw.states, raw.smoothing);
     auto scoret = clf.score(raw.Xt, raw.yt);
     REQUIRE(score == Catch::Approx(0.97333f).epsilon(raw.epsilon));
     REQUIRE(scoret == Catch::Approx(0.97333f).epsilon(raw.epsilon));
@@ -233,13 +233,13 @@ TEST_CASE("KDB with hyperparameters", "[Models]")
 {
     auto raw = RawDatasets("glass", true);
     auto clf = bayesnet::KDB(2);
-    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states);
+    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states, raw.smoothing);
     auto score = clf.score(raw.Xv, raw.yv);
     clf.setHyperparameters({
         {"k", 3},
         {"theta", 0.7},
         });
-    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states);
+    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states, raw.smoothing);
     auto scoret = clf.score(raw.Xv, raw.yv);
     REQUIRE(score == Catch::Approx(0.827103).epsilon(raw.epsilon));
     REQUIRE(scoret == Catch::Approx(0.761682).epsilon(raw.epsilon));
@@ -248,7 +248,7 @@ TEST_CASE("Incorrect type of data for SPODELd", "[Models]")
 {
     auto raw = RawDatasets("iris", true);
     auto clf = bayesnet::SPODELd(0);
-    REQUIRE_THROWS_AS(clf.fit(raw.dataset, raw.features, raw.className, raw.states), std::runtime_error);
+    REQUIRE_THROWS_AS(clf.fit(raw.dataset, raw.features, raw.className, raw.states, raw.smoothing), std::runtime_error);
 }
 TEST_CASE("Predict, predict_proba & score without fitting", "[Models]")
 {

tests/TestBayesNetwork.cc

@@ -15,6 +15,7 @@
 #include "bayesnet/network/Node.h"
 #include "bayesnet/utils/bayesnetUtils.h"
 
+const double threshold = 1e-4;
 void buildModel(bayesnet::Network& net, const std::vector<std::string>& features, const std::string& className)
 {
     std::vector<pair<int, int>> network = { {0, 1}, {0, 2}, {1, 3} };
@@ -29,13 +30,11 @@ void buildModel(bayesnet::Network& net, const std::vector<std::string>& features
         net.addEdge(className, feature);
     }
 }
-
 TEST_CASE("Test Bayesian Network", "[Network]")
 {
 
     auto raw = RawDatasets("iris", true);
     auto net = bayesnet::Network();
-    double threshold = 1e-4;
 
     SECTION("Test get features")
     {
@@ -115,9 +114,9 @@ TEST_CASE("Test Bayesian Network", "[Network]")
             REQUIRE(children == children3);
         }
         // Fit networks
-        net.fit(raw.Xv, raw.yv, raw.weightsv, raw.features, raw.className, raw.states);
-        net2.fit(raw.dataset, raw.weights, raw.features, raw.className, raw.states);
-        net3.fit(raw.Xt, raw.yt, raw.weights, raw.features, raw.className, raw.states);
+        net.fit(raw.Xv, raw.yv, raw.weightsv, raw.features, raw.className, raw.states, raw.smoothing);
+        net2.fit(raw.dataset, raw.weights, raw.features, raw.className, raw.states, raw.smoothing);
+        net3.fit(raw.Xt, raw.yt, raw.weights, raw.features, raw.className, raw.states, raw.smoothing);
         REQUIRE(net.getStates() == net2.getStates());
         REQUIRE(net.getStates() == net3.getStates());
         REQUIRE(net.getFeatures() == net2.getFeatures());
@@ -150,6 +149,7 @@ TEST_CASE("Test Bayesian Network", "[Network]")
     }
     SECTION("Test show")
     {
+        INFO("Test show");
         net.addNode("A");
         net.addNode("B");
         net.addNode("C");
@@ -163,6 +163,7 @@ TEST_CASE("Test Bayesian Network", "[Network]")
     }
     SECTION("Test topological_sort")
     {
+        INFO("Test topological sort");
         net.addNode("A");
         net.addNode("B");
         net.addNode("C");
@@ -176,6 +177,7 @@ TEST_CASE("Test Bayesian Network", "[Network]")
     }
     SECTION("Test graph")
     {
+        INFO("Test graph");
         net.addNode("A");
         net.addNode("B");
         net.addNode("C");
@@ -193,8 +195,9 @@ TEST_CASE("Test Bayesian Network", "[Network]")
     }
     SECTION("Test predict")
     {
+        INFO("Test predict");
         buildModel(net, raw.features, raw.className);
-        net.fit(raw.Xv, raw.yv, raw.weightsv, raw.features, raw.className, raw.states);
+        net.fit(raw.Xv, raw.yv, raw.weightsv, raw.features, raw.className, raw.states, raw.smoothing);
         std::vector<std::vector<int>> test = { {1, 2, 0, 1, 1}, {0, 1, 2, 0, 1}, {0, 0, 0, 0, 1}, {2, 2, 2, 2, 1} };
         std::vector<int> y_test = { 2, 2, 0, 2, 1 };
         auto y_pred = net.predict(test);
@@ -202,8 +205,9 @@ TEST_CASE("Test Bayesian Network", "[Network]")
     }
     SECTION("Test predict_proba")
     {
+        INFO("Test predict_proba");
         buildModel(net, raw.features, raw.className);
-        net.fit(raw.Xv, raw.yv, raw.weightsv, raw.features, raw.className, raw.states);
+        net.fit(raw.Xv, raw.yv, raw.weightsv, raw.features, raw.className, raw.states, raw.smoothing);
         std::vector<std::vector<int>> test = { {1, 2, 0, 1, 1}, {0, 1, 2, 0, 1}, {0, 0, 0, 0, 1}, {2, 2, 2, 2, 1} };
         std::vector<std::vector<double>> y_test = {
             {0.450237, 0.0866621, 0.463101},
@@ -223,15 +227,17 @@ TEST_CASE("Test Bayesian Network", "[Network]")
     }
     SECTION("Test score")
     {
+        INFO("Test score");
         buildModel(net, raw.features, raw.className);
-        net.fit(raw.Xv, raw.yv, raw.weightsv, raw.features, raw.className, raw.states);
+        net.fit(raw.Xv, raw.yv, raw.weightsv, raw.features, raw.className, raw.states, raw.smoothing);
         auto score = net.score(raw.Xv, raw.yv);
         REQUIRE(score == Catch::Approx(0.97333333).margin(threshold));
     }
     SECTION("Copy constructor")
     {
+        INFO("Test copy constructor");
         buildModel(net, raw.features, raw.className);
-        net.fit(raw.Xv, raw.yv, raw.weightsv, raw.features, raw.className, raw.states);
+        net.fit(raw.Xv, raw.yv, raw.weightsv, raw.features, raw.className, raw.states, raw.smoothing);
         auto net2 = bayesnet::Network(net);
         REQUIRE(net.getFeatures() == net2.getFeatures());
         REQUIRE(net.getEdges() == net2.getEdges());
@@ -253,6 +259,7 @@ TEST_CASE("Test Bayesian Network", "[Network]")
     }
     SECTION("Test oddities")
     {
+        INFO("Test oddities");
         buildModel(net, raw.features, raw.className);
         // predict without fitting
         std::vector<std::vector<int>> test = { {1, 2, 0, 1, 1}, {0, 1, 2, 0, 1}, {0, 0, 0, 0, 1}, {2, 2, 2, 2, 1} };
@@ -268,27 +275,27 @@ TEST_CASE("Test Bayesian Network", "[Network]")
         // predict with wrong data
         auto netx = bayesnet::Network();
         buildModel(netx, raw.features, raw.className);
-        netx.fit(raw.Xv, raw.yv, raw.weightsv, raw.features, raw.className, raw.states);
+        netx.fit(raw.Xv, raw.yv, raw.weightsv, raw.features, raw.className, raw.states, raw.smoothing);
         std::vector<std::vector<int>> test2 = { {1, 2, 0, 1, 1}, {0, 1, 2, 0, 1}, {0, 0, 0, 0, 1} };
         auto test_tensor2 = bayesnet::vectorToTensor(test2, false);
-        REQUIRE_THROWS_AS(netx.predict(test2), std::logic_error);
-        REQUIRE_THROWS_WITH(netx.predict(test2), "Sample size (3) does not match the number of features (4)");
-        REQUIRE_THROWS_AS(netx.predict(test_tensor2), std::logic_error);
-        REQUIRE_THROWS_WITH(netx.predict(test_tensor2), "Sample size (3) does not match the number of features (4)");
+        REQUIRE_THROWS_AS(netx.predict(test2), std::invalid_argument);
+        REQUIRE_THROWS_WITH(netx.predict(test2), "(V) Sample size (3) does not match the number of features (4)");
+        REQUIRE_THROWS_AS(netx.predict(test_tensor2), std::invalid_argument);
+        REQUIRE_THROWS_WITH(netx.predict(test_tensor2), "(T) Sample size (3) does not match the number of features (4)");
         // fit with wrong data
         // Weights
         auto net2 = bayesnet::Network();
-        REQUIRE_THROWS_AS(net2.fit(raw.Xv, raw.yv, std::vector<double>(), raw.features, raw.className, raw.states), std::invalid_argument);
+        REQUIRE_THROWS_AS(net2.fit(raw.Xv, raw.yv, std::vector<double>(), raw.features, raw.className, raw.states, raw.smoothing), std::invalid_argument);
         std::string invalid_weights = "Weights (0) must have the same number of elements as samples (150) in Network::fit";
-        REQUIRE_THROWS_WITH(net2.fit(raw.Xv, raw.yv, std::vector<double>(), raw.features, raw.className, raw.states), invalid_weights);
+        REQUIRE_THROWS_WITH(net2.fit(raw.Xv, raw.yv, std::vector<double>(), raw.features, raw.className, raw.states, raw.smoothing), invalid_weights);
         // X & y
         std::string invalid_labels = "X and y must have the same number of samples in Network::fit (150 != 0)";
-        REQUIRE_THROWS_AS(net2.fit(raw.Xv, std::vector<int>(), raw.weightsv, raw.features, raw.className, raw.states), std::invalid_argument);
-        REQUIRE_THROWS_WITH(net2.fit(raw.Xv, std::vector<int>(), raw.weightsv, raw.features, raw.className, raw.states), invalid_labels);
+        REQUIRE_THROWS_AS(net2.fit(raw.Xv, std::vector<int>(), raw.weightsv, raw.features, raw.className, raw.states, raw.smoothing), std::invalid_argument);
+        REQUIRE_THROWS_WITH(net2.fit(raw.Xv, std::vector<int>(), raw.weightsv, raw.features, raw.className, raw.states, raw.smoothing), invalid_labels);
         // Features
         std::string invalid_features = "X and features must have the same number of features in Network::fit (4 != 0)";
-        REQUIRE_THROWS_AS(net2.fit(raw.Xv, raw.yv, raw.weightsv, std::vector<std::string>(), raw.className, raw.states), std::invalid_argument);
-        REQUIRE_THROWS_WITH(net2.fit(raw.Xv, raw.yv, raw.weightsv, std::vector<std::string>(), raw.className, raw.states), invalid_features);
+        REQUIRE_THROWS_AS(net2.fit(raw.Xv, raw.yv, raw.weightsv, std::vector<std::string>(), raw.className, raw.states, raw.smoothing), std::invalid_argument);
+        REQUIRE_THROWS_WITH(net2.fit(raw.Xv, raw.yv, raw.weightsv, std::vector<std::string>(), raw.className, raw.states, raw.smoothing), invalid_features);
         // Different number of features
         auto net3 = bayesnet::Network();
         auto test2y = { 1, 2, 3, 4, 5 };
@@ -296,23 +303,23 @@ TEST_CASE("Test Bayesian Network", "[Network]")
         auto features3 = raw.features;
         features3.pop_back();
         std::string invalid_features2 = "X and local features must have the same number of features in Network::fit (3 != 4)";
-        REQUIRE_THROWS_AS(net3.fit(test2, test2y, std::vector<double>(5, 0), features3, raw.className, raw.states), std::invalid_argument);
-        REQUIRE_THROWS_WITH(net3.fit(test2, test2y, std::vector<double>(5, 0), features3, raw.className, raw.states), invalid_features2);
+        REQUIRE_THROWS_AS(net3.fit(test2, test2y, std::vector<double>(5, 0), features3, raw.className, raw.states, raw.smoothing), std::invalid_argument);
+        REQUIRE_THROWS_WITH(net3.fit(test2, test2y, std::vector<double>(5, 0), features3, raw.className, raw.states, raw.smoothing), invalid_features2);
         // Uninitialized network
         std::string network_invalid = "The network has not been initialized. You must call addNode() before calling fit()";
-        REQUIRE_THROWS_AS(net2.fit(raw.Xv, raw.yv, raw.weightsv, raw.features, "duck", raw.states), std::invalid_argument);
-        REQUIRE_THROWS_WITH(net2.fit(raw.Xv, raw.yv, raw.weightsv, raw.features, "duck", raw.states), network_invalid);
+        REQUIRE_THROWS_AS(net2.fit(raw.Xv, raw.yv, raw.weightsv, raw.features, "duck", raw.states, raw.smoothing), std::invalid_argument);
+        REQUIRE_THROWS_WITH(net2.fit(raw.Xv, raw.yv, raw.weightsv, raw.features, "duck", raw.states, raw.smoothing), network_invalid);
         // Classname
         std::string invalid_classname = "Class Name not found in Network::features";
-        REQUIRE_THROWS_AS(net.fit(raw.Xv, raw.yv, raw.weightsv, raw.features, "duck", raw.states), std::invalid_argument);
-        REQUIRE_THROWS_WITH(net.fit(raw.Xv, raw.yv, raw.weightsv, raw.features, "duck", raw.states), invalid_classname);
+        REQUIRE_THROWS_AS(net.fit(raw.Xv, raw.yv, raw.weightsv, raw.features, "duck", raw.states, raw.smoothing), std::invalid_argument);
+        REQUIRE_THROWS_WITH(net.fit(raw.Xv, raw.yv, raw.weightsv, raw.features, "duck", raw.states, raw.smoothing), invalid_classname);
         // Invalid feature
         auto features2 = raw.features;
         features2.pop_back();
         features2.push_back("duck");
         std::string invalid_feature = "Feature duck not found in Network::features";
-        REQUIRE_THROWS_AS(net.fit(raw.Xv, raw.yv, raw.weightsv, features2, raw.className, raw.states), std::invalid_argument);
-        REQUIRE_THROWS_WITH(net.fit(raw.Xv, raw.yv, raw.weightsv, features2, raw.className, raw.states), invalid_feature);
+        REQUIRE_THROWS_AS(net.fit(raw.Xv, raw.yv, raw.weightsv, features2, raw.className, raw.states, raw.smoothing), std::invalid_argument);
+        REQUIRE_THROWS_WITH(net.fit(raw.Xv, raw.yv, raw.weightsv, features2, raw.className, raw.states, raw.smoothing), invalid_feature);
         // Add twice the same node name to the network => Nothing should happen
         net.addNode("A");
         net.addNode("A");
@@ -320,8 +327,8 @@ TEST_CASE("Test Bayesian Network", "[Network]")
         auto net4 = bayesnet::Network();
         buildModel(net4, raw.features, raw.className);
         std::string invalid_state = "Feature sepallength not found in states";
-        REQUIRE_THROWS_AS(net4.fit(raw.Xv, raw.yv, raw.weightsv, raw.features, raw.className, std::map<std::string, std::vector<int>>()), std::invalid_argument);
-        REQUIRE_THROWS_WITH(net4.fit(raw.Xv, raw.yv, raw.weightsv, raw.features, raw.className, std::map<std::string, std::vector<int>>()), invalid_state);
+        REQUIRE_THROWS_AS(net4.fit(raw.Xv, raw.yv, raw.weightsv, raw.features, raw.className, std::map<std::string, std::vector<int>>(), raw.smoothing), std::invalid_argument);
+        REQUIRE_THROWS_WITH(net4.fit(raw.Xv, raw.yv, raw.weightsv, raw.features, raw.className, std::map<std::string, std::vector<int>>(), raw.smoothing), invalid_state);
     }
 
 }
@@ -342,15 +349,6 @@ TEST_CASE("Cicle in Network", "[Network]")
     REQUIRE_THROWS_AS(net.addEdge("C", "A"), std::invalid_argument);
     REQUIRE_THROWS_WITH(net.addEdge("C", "A"), "Adding this edge forms a cycle in the graph.");
 }
-TEST_CASE("Test max threads constructor", "[Network]")
-{
-    auto net = bayesnet::Network();
-    REQUIRE(net.getMaxThreads() == 0.95f);
-    auto net2 = bayesnet::Network(4);
-    REQUIRE(net2.getMaxThreads() == 4);
-    auto net3 = bayesnet::Network(1.75);
-    REQUIRE(net3.getMaxThreads() == 1.75);
-}
 TEST_CASE("Edges troubles", "[Network]")
 {
     auto net = bayesnet::Network();
@@ -360,13 +358,16 @@ TEST_CASE("Edges troubles", "[Network]")
     REQUIRE_THROWS_WITH(net.addEdge("A", "C"), "Child node C does not exist");
     REQUIRE_THROWS_AS(net.addEdge("C", "A"), std::invalid_argument);
     REQUIRE_THROWS_WITH(net.addEdge("C", "A"), "Parent node C does not exist");
+    net.addEdge("A", "B");
+    REQUIRE_THROWS_AS(net.addEdge("A", "B"), std::invalid_argument);
+    REQUIRE_THROWS_WITH(net.addEdge("A", "B"), "Edge A -> B already exists");
 }
 TEST_CASE("Dump CPT", "[Network]")
 {
     auto net = bayesnet::Network();
     auto raw = RawDatasets("iris", true);
     buildModel(net, raw.features, raw.className);
-    net.fit(raw.Xv, raw.yv, raw.weightsv, raw.features, raw.className, raw.states);
+    net.fit(raw.Xv, raw.yv, raw.weightsv, raw.features, raw.className, raw.states, raw.smoothing);
     auto res = net.dump_cpt();
     std::string expected = R"(* class: (3) : [3]
  0.3333
@@ -459,3 +460,108 @@ TEST_CASE("Dump CPT", "[Network]")
     REQUIRE(res == expected);
 }
 
+TEST_CASE("Test Smoothing A", "[Network]")
+{
+    /*
+    Tomando m = 1 Pa = 0.5
+    Si estoy calculando P(A | C), con C en{ 0,1,2 } y tengo :
+    AC = { 11, 12, 11, 10, 10, 12, 10, 01, 00, 02 }
+    Entonces:
+    P(A = 1 | C = 0) = (3 + 1 / 2 * 1) / (4 + 1) = 3.5 / 5
+    P(A = 0 | C = 0) = (1 + 1 / 2 * 1) / (4 + 1) = 1.5 / 5
+    Donde m aquí es el número de veces de C = 0 que es la que condiciona y la a priori vuelve a ser sobre A que es sobre las que estaríamos calculando esas marginales.
+    P(A = 1 | C = 1) = (2 + 1 / 2 * 1) / (3 + 1) = 2.5 / 4
+    P(A = 0 | C = 1) = (1 + 1 / 2 * 1) / (3 + 1) = 1.5 / 4
+    P(A = 1 | C = 2) = (2 + 1 / 2 * 1) / (3 + 1) = 2.5 / 5
+    P(A = 0 | C = 2) = (1 + 1 / 2 * 1) / (3 + 1) = 1.5 / 5
+    En realidad es parecido a Laplace, que en este caso p.e.con C = 0 sería
+    P(A = 1 | C = 0) = (3 + 1) / (4 + 2) = 4 / 6
+    P(A = 0 | C = 0) = (1 + 1) / (4 + 2) = 2 / 6
+    */
+    auto net = bayesnet::Network();
+    net.addNode("A");
+    net.addNode("C");
+    net.addEdge("C", "A");
+    std::vector<int> C = { 1, 2, 1, 0, 0, 2, 0, 1, 0, 2 };
+    std::vector<std::vector<int>> A = { { 1, 1, 1, 1, 1, 1, 1, 0, 0, 0 } };
+    std::map<std::string, std::vector<int>> states = { { "A", {0, 1} }, { "C", {0, 1, 2} } };
+    auto weights = std::vector<double>(C.size(), 1);
+    //
+    // Laplace
+    //
+    net.fit(A, C, weights, { "A" }, "C", states, bayesnet::Smoothing_t::LAPLACE);
+    auto cpt_c_laplace = net.getNodes().at("C")->getCPT();
+    REQUIRE(cpt_c_laplace.size(0) == 3);
+    auto laplace_c = std::vector<float>({ 0.3846, 0.3077, 0.3077 });
+    for (int i = 0; i < laplace_c.size(); ++i) {
+        REQUIRE(cpt_c_laplace.index({ i }).item<float>() == Catch::Approx(laplace_c[i]).margin(threshold));
+    }
+    auto cpt_a_laplace = net.getNodes().at("A")->getCPT();
+    REQUIRE(cpt_a_laplace.size(0) == 2);
+    REQUIRE(cpt_a_laplace.size(1) == 3);
+    auto laplace_a = std::vector<std::vector<float>>({ {0.3333, 0.4000,0.4000}, {0.6667,  0.6000,  0.6000} });
+    for (int i = 0; i < 2; ++i) {
+        for (int j = 0; j < 3; ++j) {
+            REQUIRE(cpt_a_laplace.index({ i, j }).item<float>() == Catch::Approx(laplace_a[i][j]).margin(threshold));
+        }
+    }
+    //
+    // Cestnik
+    //
+    net.fit(A, C, weights, { "A" }, "C", states, bayesnet::Smoothing_t::CESTNIK);
+    auto cpt_c_cestnik = net.getNodes().at("C")->getCPT();
+    REQUIRE(cpt_c_cestnik.size(0) == 3);
+    auto cestnik_c = std::vector<float>({ 0.3939, 0.3030, 0.3030 });
+    for (int i = 0; i < laplace_c.size(); ++i) {
+        REQUIRE(cpt_c_cestnik.index({ i }).item<float>() == Catch::Approx(cestnik_c[i]).margin(threshold));
+    }
+    auto cpt_a_cestnik = net.getNodes().at("A")->getCPT();
+    REQUIRE(cpt_a_cestnik.size(0) == 2);
+    REQUIRE(cpt_a_cestnik.size(1) == 3);
+    auto cestnik_a = std::vector<std::vector<float>>({ {0.3000, 0.3750, 0.3750}, {0.7000, 0.6250, 0.6250} });
+    for (int i = 0; i < 2; ++i) {
+        for (int j = 0; j < 3; ++j) {
+            REQUIRE(cpt_a_cestnik.index({ i, j }).item<float>() == Catch::Approx(cestnik_a[i][j]).margin(threshold));
+        }
+    }
+}
+TEST_CASE("Test Smoothing B", "[Network]")
+{
+    auto net = bayesnet::Network();
+    net.addNode("X");
+    net.addNode("Y");
+    net.addNode("Z");
+    net.addNode("C");
+    net.addEdge("C", "X");
+    net.addEdge("C", "Y");
+    net.addEdge("C", "Z");
+    net.addEdge("Y", "Z");
+    std::vector<int> C = { 0, 1, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1 };
+    std::vector<std::vector<int>> Data = {
+        { 0,0,1,1,0,1,0,1,0,1,0,0,0,1,0,1,0,0},
+        { 1,2,0,2,2,2,1,0,0,1,1,1,0,1,2,1,0,2},
+        { 2,1,3,3,2,0,0,1,3,2,1,2,2,3,0,0,1,2}
+    };
+    std::map<std::string, std::vector<int>> states = {
+        { "X", {0, 1} },
+        { "Y", {0, 1, 2} },
+        { "Z", {0, 1, 2, 3} },
+        { "C", {0, 1} }
+    };
+    auto weights = std::vector<double>(C.size(), 1);
+    // Simple
+    std::cout << "LAPLACE\n";
+    net.fit(Data, C, weights, { "X", "Y", "Z" }, "C", states, bayesnet::Smoothing_t::LAPLACE);
+    std::cout << net.dump_cpt();
+    std::cout << "Predict proba of {0, 1, 2} y {1, 2, 3} = " << net.predict_proba({ {0, 1}, {1, 2}, {2, 3} }) << std::endl;
+    std::cout << "ORIGINAL\n";
+    net.fit(Data, C, weights, { "X", "Y", "Z" }, "C", states, bayesnet::Smoothing_t::ORIGINAL);
+    std::cout << net.dump_cpt();
+    std::cout << "Predict proba of {0, 1, 2} y {1, 2, 3} = " << net.predict_proba({ {0, 1}, {1, 2}, {2, 3} }) << std::endl;
+    std::cout << "CESTNIK\n";
+    net.fit(Data, C, weights, { "X", "Y", "Z" }, "C", states, bayesnet::Smoothing_t::CESTNIK);
+    std::cout << net.dump_cpt();
+    std::cout << "Predict proba of {0, 1, 2} y {1, 2, 3} = " << net.predict_proba({ {0, 1}, {1, 2}, {2, 3} }) << std::endl;
+
+
+}

tests/TestBoostA2DE.cc

@@ -17,7 +17,7 @@ TEST_CASE("Build basic model", "[BoostA2DE]")
 {
     auto raw = RawDatasets("diabetes", true);
     auto clf = bayesnet::BoostA2DE();
-    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states);
+    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states, raw.smoothing);
     REQUIRE(clf.getNumberOfNodes() == 342);
     REQUIRE(clf.getNumberOfEdges() == 684);
     REQUIRE(clf.getNotes().size() == 3);
@@ -32,7 +32,7 @@ TEST_CASE("Build basic model", "[BoostA2DE]")
 //     auto raw = RawDatasets("glass", true);
 //     auto clf = bayesnet::BoostAODE();
 //     clf.setHyperparameters({ {"select_features", "IWSS"}, {"threshold", 0.5 } });
-//     clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states);
+//     clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states, raw.smoothing);
 //     REQUIRE(clf.getNumberOfNodes() == 90);
 //     REQUIRE(clf.getNumberOfEdges() == 153);
 //     REQUIRE(clf.getNotes().size() == 2);
@@ -44,7 +44,7 @@ TEST_CASE("Build basic model", "[BoostA2DE]")
 //     auto raw = RawDatasets("glass", true);
 //     auto clf = bayesnet::BoostAODE();
 //     clf.setHyperparameters({ {"select_features", "FCBF"}, {"threshold", 1e-7 } });
-//     clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states);
+//     clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states, raw.smoothing);
 //     REQUIRE(clf.getNumberOfNodes() == 90);
 //     REQUIRE(clf.getNumberOfEdges() == 153);
 //     REQUIRE(clf.getNotes().size() == 2);
@@ -60,7 +60,7 @@ TEST_CASE("Build basic model", "[BoostA2DE]")
 //         {"convergence", true},
 //         {"select_features","CFS"},
 //         });
-//     clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states);
+//     clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states, raw.smoothing);
 //     REQUIRE(clf.getNumberOfNodes() == 72);
 //     REQUIRE(clf.getNumberOfEdges() == 120);
 //     REQUIRE(clf.getNotes().size() == 2);
@@ -75,7 +75,7 @@ TEST_CASE("Build basic model", "[BoostA2DE]")
 // {
 //     auto raw = RawDatasets("iris", true);
 //     auto clf = bayesnet::BoostAODE(false);
-//     clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states);
+//     clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states, raw.smoothing);
 //     auto score_proba = clf.score(raw.Xv, raw.yv);
 //     auto pred_proba = clf.predict_proba(raw.Xv);
 //     clf.setHyperparameters({
@@ -104,7 +104,7 @@ TEST_CASE("Build basic model", "[BoostA2DE]")
 //             {"maxTolerance", 1},
 //             {"convergence", false},
 //             });
-//         clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states);
+//         clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states, raw.smoothing);
 //         auto score = clf.score(raw.Xv, raw.yv);
 //         auto scoret = clf.score(raw.Xt, raw.yt);
 //         INFO("BoostAODE order: " + order);
@@ -136,7 +136,7 @@ TEST_CASE("Build basic model", "[BoostA2DE]")
 //     for (const auto& hyper : bad_hyper_fit.items()) {
 //         INFO("BoostAODE hyper: " + hyper.value().dump());
 //         clf.setHyperparameters(hyper.value());
-//         REQUIRE_THROWS_AS(clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states), std::invalid_argument);
+//         REQUIRE_THROWS_AS(clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states, raw.smoothing, std::invalid_argument);
 //     }
 // }
 
@@ -151,7 +151,7 @@ TEST_CASE("Build basic model", "[BoostA2DE]")
 //         {"block_update", false},
 //         {"convergence_best", false},
 //         });
-//     clf.fit(raw.X_train, raw.y_train, raw.features, raw.className, raw.states);
+//     clf.fit(raw.X_train, raw.y_train, raw.features, raw.className, raw.states, raw.smoothing);
 //     REQUIRE(clf.getNumberOfNodes() == 210);
 //     REQUIRE(clf.getNumberOfEdges() == 378);
 //     REQUIRE(clf.getNotes().size() == 1);
@@ -172,13 +172,13 @@ TEST_CASE("Build basic model", "[BoostA2DE]")
 //         {"convergence_best", true},
 //     };
 //     clf.setHyperparameters(hyperparameters);
-//     clf.fit(raw.X_train, raw.y_train, raw.features, raw.className, raw.states);
+//     clf.fit(raw.X_train, raw.y_train, raw.features, raw.className, raw.states, raw.smoothing);
 //     auto score_best = clf.score(raw.X_test, raw.y_test);
 //     REQUIRE(score_best == Catch::Approx(0.980000019f).epsilon(raw.epsilon));
 //     // Now we will set the hyperparameter to use the last accuracy
 //     hyperparameters["convergence_best"] = false;
 //     clf.setHyperparameters(hyperparameters);
-//     clf.fit(raw.X_train, raw.y_train, raw.features, raw.className, raw.states);
+//     clf.fit(raw.X_train, raw.y_train, raw.features, raw.className, raw.states, raw.smoothing);
 //     auto score_last = clf.score(raw.X_test, raw.y_test);
 //     REQUIRE(score_last == Catch::Approx(0.976666689f).epsilon(raw.epsilon));
 // }
@@ -193,7 +193,7 @@ TEST_CASE("Build basic model", "[BoostA2DE]")
 //         {"maxTolerance", 3},
 //         {"convergence", true},
 //         });
-//     clf.fit(raw.X_train, raw.y_train, raw.features, raw.className, raw.states);
+//     clf.fit(raw.X_train, raw.y_train, raw.features, raw.className, raw.states, raw.smoothing);
 //     REQUIRE(clf.getNumberOfNodes() == 868);
 //     REQUIRE(clf.getNumberOfEdges() == 1724);
 //     REQUIRE(clf.getNotes().size() == 3);

tests/TestBoostAODE.cc

@@ -18,7 +18,7 @@ TEST_CASE("Feature_select CFS", "[BoostAODE]")
     auto raw = RawDatasets("glass", true);
     auto clf = bayesnet::BoostAODE();
     clf.setHyperparameters({ {"select_features", "CFS"} });
-    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states);
+    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states, raw.smoothing);
     REQUIRE(clf.getNumberOfNodes() == 90);
     REQUIRE(clf.getNumberOfEdges() == 153);
     REQUIRE(clf.getNotes().size() == 2);
@@ -30,7 +30,7 @@ TEST_CASE("Feature_select IWSS", "[BoostAODE]")
     auto raw = RawDatasets("glass", true);
     auto clf = bayesnet::BoostAODE();
     clf.setHyperparameters({ {"select_features", "IWSS"}, {"threshold", 0.5 } });
-    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states);
+    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states, raw.smoothing);
     REQUIRE(clf.getNumberOfNodes() == 90);
     REQUIRE(clf.getNumberOfEdges() == 153);
     REQUIRE(clf.getNotes().size() == 2);
@@ -42,7 +42,7 @@ TEST_CASE("Feature_select FCBF", "[BoostAODE]")
     auto raw = RawDatasets("glass", true);
     auto clf = bayesnet::BoostAODE();
     clf.setHyperparameters({ {"select_features", "FCBF"}, {"threshold", 1e-7 } });
-    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states);
+    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states, raw.smoothing);
     REQUIRE(clf.getNumberOfNodes() == 90);
     REQUIRE(clf.getNumberOfEdges() == 153);
     REQUIRE(clf.getNotes().size() == 2);
@@ -58,7 +58,7 @@ TEST_CASE("Test used features in train note and score", "[BoostAODE]")
         {"convergence", true},
         {"select_features","CFS"},
         });
-    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states);
+    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states, raw.smoothing);
     REQUIRE(clf.getNumberOfNodes() == 72);
     REQUIRE(clf.getNumberOfEdges() == 120);
     REQUIRE(clf.getNotes().size() == 2);
@@ -73,7 +73,7 @@ TEST_CASE("Voting vs proba", "[BoostAODE]")
 {
     auto raw = RawDatasets("iris", true);
     auto clf = bayesnet::BoostAODE(false);
-    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states);
+    clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states, raw.smoothing);
     auto score_proba = clf.score(raw.Xv, raw.yv);
     auto pred_proba = clf.predict_proba(raw.Xv);
     clf.setHyperparameters({
@@ -102,7 +102,7 @@ TEST_CASE("Order asc, desc & random", "[BoostAODE]")
             {"maxTolerance", 1},
             {"convergence", false},
             });
-        clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states);
+        clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states, raw.smoothing);
         auto score = clf.score(raw.Xv, raw.yv);
         auto scoret = clf.score(raw.Xt, raw.yt);
         INFO("BoostAODE order: " << order);
@@ -134,7 +134,7 @@ TEST_CASE("Oddities", "[BoostAODE]")
     for (const auto& hyper : bad_hyper_fit.items()) {
         INFO("BoostAODE hyper: " << hyper.value().dump());
         clf.setHyperparameters(hyper.value());
-        REQUIRE_THROWS_AS(clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states), std::invalid_argument);
+        REQUIRE_THROWS_AS(clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states, raw.smoothing), std::invalid_argument);
     }
 }
 
@@ -149,7 +149,7 @@ TEST_CASE("Bisection Best", "[BoostAODE]")
         {"block_update", false},
         {"convergence_best", false},
         });
-    clf.fit(raw.X_train, raw.y_train, raw.features, raw.className, raw.states);
+    clf.fit(raw.X_train, raw.y_train, raw.features, raw.className, raw.states, raw.smoothing);
     REQUIRE(clf.getNumberOfNodes() == 210);
     REQUIRE(clf.getNumberOfEdges() == 378);
     REQUIRE(clf.getNotes().size() == 1);
@@ -170,13 +170,13 @@ TEST_CASE("Bisection Best vs Last", "[BoostAODE]")
         {"convergence_best", true},
     };
     clf.setHyperparameters(hyperparameters);
-    clf.fit(raw.X_train, raw.y_train, raw.features, raw.className, raw.states);
+    clf.fit(raw.X_train, raw.y_train, raw.features, raw.className, raw.states, raw.smoothing);
     auto score_best = clf.score(raw.X_test, raw.y_test);
     REQUIRE(score_best == Catch::Approx(0.980000019f).epsilon(raw.epsilon));
     // Now we will set the hyperparameter to use the last accuracy
     hyperparameters["convergence_best"] = false;
     clf.setHyperparameters(hyperparameters);
-    clf.fit(raw.X_train, raw.y_train, raw.features, raw.className, raw.states);
+    clf.fit(raw.X_train, raw.y_train, raw.features, raw.className, raw.states, raw.smoothing);
     auto score_last = clf.score(raw.X_test, raw.y_test);
     REQUIRE(score_last == Catch::Approx(0.976666689f).epsilon(raw.epsilon));
 }
@@ -191,7 +191,7 @@ TEST_CASE("Block Update", "[BoostAODE]")
         {"maxTolerance", 3},
         {"convergence", true},
         });
-    clf.fit(raw.X_train, raw.y_train, raw.features, raw.className, raw.states);
+    clf.fit(raw.X_train, raw.y_train, raw.features, raw.className, raw.states, raw.smoothing);
     REQUIRE(clf.getNumberOfNodes() == 868);
     REQUIRE(clf.getNumberOfEdges() == 1724);
     REQUIRE(clf.getNotes().size() == 3);

tests/TestModulesVersions.cc

@@ -7,7 +7,7 @@
 #include <catch2/catch_test_macros.hpp>
 #include <catch2/matchers/catch_matchers.hpp>
 #include <string>
-#include <CPPFImdlp.h>
+#include <fimdlp/CPPFImdlp.h>
 #include <folding.hpp>
 #include <nlohmann/json.hpp>
 #define TO_STR2(x) #x
@@ -16,10 +16,10 @@
 #include "TestUtils.h"
 
 std::map<std::string, std::string> modules = {
-    { "mdlp", "1.1.2" },
+    { "mdlp", "2.0.0" },
     { "Folding", "1.1.0" },
     { "json", "3.11" },
-    { "ArffFiles", "1.0.0" }
+    { "ArffFiles", "1.1.0" }
 };
 
 TEST_CASE("MDLP", "[Modules]")

tests/TestUtils.h

@@ -12,8 +12,9 @@
 #include <map>
 #include <tuple>
 #include <ArffFiles.hpp>
-#include <CPPFImdlp.h>
+#include <fimdlp/CPPFImdlp.h>
 #include <folding.hpp>
+#include <bayesnet/network/Network.h>
 
 
 class RawDatasets {
@@ -32,6 +33,7 @@ public:
     bool discretize;
     int num_samples = 0;
     bool shuffle = false;
+    bayesnet::Smoothing_t smoothing = bayesnet::Smoothing_t::ORIGINAL;
 private:
     std::string to_string()
     {