Create doc target to build documentation

There were cases where a tiny negative number was returned (less than -1e-7)
Fix mst glass test that is affected with this change

.gitignore

@@ -40,4 +40,8 @@ puml/**
 .vscode/settings.json
 sample/build
 **/.DS_Store
+docs/manual
+docs/man3
+docs/man
+docs/Doxyfile
 

.gitmodules

@@ -18,3 +18,6 @@
         url = https://github.com/catchorg/Catch2.git
 	main = main
 	update = merge
+[submodule "tests/lib/Files"]
+	path = tests/lib/Files
+	url = https://github.com/rmontanana/ArffFiles

CMakeLists.txt

@@ -25,8 +25,11 @@ set(CMAKE_CXX_EXTENSIONS                  OFF)
 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 -fno-default-inline")
+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")
+if (NOT ${CMAKE_SYSTEM_NAME} MATCHES "Darwin")
+    set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -fno-default-inline")
+endif()
 
 # Options
 # -------
@@ -63,7 +66,6 @@ endif (ENABLE_CLANG_TIDY)
 # include(FetchContent)
 add_git_submodule("lib/json")
 add_git_submodule("lib/mdlp")
-add_subdirectory("lib/Files")
 
 # Subdirectories
 # --------------
@@ -87,3 +89,14 @@ install(TARGETS BayesNet
         CONFIGURATIONS Release)
 install(DIRECTORY bayesnet/ DESTINATION include/bayesnet FILES_MATCHING CONFIGURATIONS Release PATTERN "*.h")
 install(FILES ${CMAKE_BINARY_DIR}/configured_files/include/bayesnet/config.h DESTINATION include/bayesnet CONFIGURATIONS Release)
+
+# Documentation
+# -------------
+find_package(Doxygen REQUIRED dot)
+set(DOC_DIR ${CMAKE_CURRENT_SOURCE_DIR}/docs)
+set(doxyfile_in ${DOC_DIR}/Doxyfile.in)
+set(doxyfile ${DOC_DIR}/Doxyfile)
+configure_file(${doxyfile_in} ${doxyfile} @ONLY)
+doxygen_add_docs(doxygen
+    WORKING_DIRECTORY ${DOC_DIR}
+    CONFIG_FILE ${doxyfile})

Makefile

@@ -1,6 +1,6 @@
 SHELL := /bin/bash
 .DEFAULT_GOAL := help
-.PHONY: viewcoverage coverage setup help install uninstall diagrams buildr buildd test clean debug release sample updatebadge
+.PHONY: viewcoverage coverage setup help install uninstall diagrams buildr buildd test clean debug release sample updatebadge doc
 
 f_release = build_Release
 f_debug = build_Debug
@@ -145,6 +145,11 @@ updatebadge: ## Update the coverage badge in README.md
 	@env python update_coverage.py $(f_debug)/tests
 	@echo ">>> Done";
 
+doc: ## Generate documentation
+	@echo ">>> Generating documentation..."
+	@cmake --build $(f_release) -t doxygen
+	@echo ">>> Done";
+
 help: ## Show help message
 	@IFS=$$'\n' ; \
 	help_lines=(`fgrep -h "##" $(MAKEFILE_LIST) | fgrep -v fgrep | sed -e 's/\\$$//' | sed -e 's/##/:/'`); \

bayesnet/CMakeLists.txt

@@ -1,6 +1,5 @@
 include_directories(
     ${BayesNet_SOURCE_DIR}/lib/mdlp
-    ${BayesNet_SOURCE_DIR}/lib/Files
     ${BayesNet_SOURCE_DIR}/lib/folding
     ${BayesNet_SOURCE_DIR}/lib/json/include
     ${BayesNet_SOURCE_DIR}

bayesnet/classifiers/Proposal.cc

@@ -4,7 +4,6 @@
 // SPDX-License-Identifier: MIT
 // ***************************************************************
 
-#include <ArffFiles.h>
 #include "Proposal.h"
 
 namespace bayesnet {
@@ -54,8 +53,7 @@ namespace bayesnet {
                     yJoinParents[i] += to_string(pDataset.index({ idx, i }).item<int>());
                 }
             }
-            auto arff = ArffFiles();
-            auto yxv = arff.factorize(yJoinParents);
+            auto yxv = factorize(yJoinParents);
             auto xvf_ptr = Xf.index({ index }).data_ptr<float>();
             auto xvf = std::vector<mdlp::precision_t>(xvf_ptr, xvf_ptr + Xf.size(1));
             discretizers[feature]->fit(xvf, yxv);
@@ -113,4 +111,19 @@ namespace bayesnet {
         }
         return Xtd;
     }
+    std::vector<int> Proposal::factorize(const std::vector<std::string>& labels_t)
+    {
+        std::vector<int> yy;
+        yy.reserve(labels_t.size());
+        std::map<std::string, int> labelMap;
+        int i = 0;
+        for (const std::string& label : labels_t) {
+            if (labelMap.find(label) == labelMap.end()) {
+                labelMap[label] = i++;
+                bool allDigits = std::all_of(label.begin(), label.end(), ::isdigit);
+            }
+            yy.push_back(labelMap[label]);
+        }
+        return yy;
+    }
 }

bayesnet/classifiers/Proposal.h

@@ -27,6 +27,7 @@ namespace bayesnet {
         torch::Tensor y; // y discrete nx1 tensor
         map<std::string, mdlp::CPPFImdlp*> discretizers;
     private:
+        std::vector<int> factorize(const std::vector<std::string>& labels_t);
         torch::Tensor& pDataset; // (n+1)xm tensor
         std::vector<std::string>& pFeatures;
         std::string& pClassName;

bayesnet/ensembles/Boost.cc

@@ -0,0 +1,246 @@
+// ***************************************************************
+// SPDX-FileCopyrightText: Copyright 2024 Ricardo Montañana Gómez
+// SPDX-FileType: SOURCE
+// SPDX-License-Identifier: MIT
+// ***************************************************************
+#include <folding.hpp>
+#include "bayesnet/feature_selection/CFS.h"
+#include "bayesnet/feature_selection/FCBF.h"
+#include "bayesnet/feature_selection/IWSS.h"
+#include "Boost.h"
+
+namespace bayesnet {
+    Boost::Boost(bool predict_voting) : Ensemble(predict_voting)
+    {
+        validHyperparameters = { "order", "convergence", "convergence_best", "bisection", "threshold", "maxTolerance",
+        "predict_voting", "select_features", "block_update" };
+    }
+    void Boost::setHyperparameters(const nlohmann::json& hyperparameters_)
+    {
+        auto hyperparameters = hyperparameters_;
+        if (hyperparameters.contains("order")) {
+            std::vector<std::string> algos = { Orders.ASC, Orders.DESC, Orders.RAND };
+            order_algorithm = hyperparameters["order"];
+            if (std::find(algos.begin(), algos.end(), order_algorithm) == algos.end()) {
+                throw std::invalid_argument("Invalid order algorithm, valid values [" + Orders.ASC + ", " + Orders.DESC + ", " + Orders.RAND + "]");
+            }
+            hyperparameters.erase("order");
+        }
+        if (hyperparameters.contains("convergence")) {
+            convergence = hyperparameters["convergence"];
+            hyperparameters.erase("convergence");
+        }
+        if (hyperparameters.contains("convergence_best")) {
+            convergence_best = hyperparameters["convergence_best"];
+            hyperparameters.erase("convergence_best");
+        }
+        if (hyperparameters.contains("bisection")) {
+            bisection = hyperparameters["bisection"];
+            hyperparameters.erase("bisection");
+        }
+        if (hyperparameters.contains("threshold")) {
+            threshold = hyperparameters["threshold"];
+            hyperparameters.erase("threshold");
+        }
+        if (hyperparameters.contains("maxTolerance")) {
+            maxTolerance = hyperparameters["maxTolerance"];
+            if (maxTolerance < 1 || maxTolerance > 4)
+                throw std::invalid_argument("Invalid maxTolerance value, must be greater in [1, 4]");
+            hyperparameters.erase("maxTolerance");
+        }
+        if (hyperparameters.contains("predict_voting")) {
+            predict_voting = hyperparameters["predict_voting"];
+            hyperparameters.erase("predict_voting");
+        }
+        if (hyperparameters.contains("select_features")) {
+            auto selectedAlgorithm = hyperparameters["select_features"];
+            std::vector<std::string> algos = { SelectFeatures.IWSS, SelectFeatures.CFS, SelectFeatures.FCBF };
+            selectFeatures = true;
+            select_features_algorithm = selectedAlgorithm;
+            if (std::find(algos.begin(), algos.end(), selectedAlgorithm) == algos.end()) {
+                throw std::invalid_argument("Invalid selectFeatures value, valid values [" + SelectFeatures.IWSS + ", " + SelectFeatures.CFS + ", " + SelectFeatures.FCBF + "]");
+            }
+            hyperparameters.erase("select_features");
+        }
+        if (hyperparameters.contains("block_update")) {
+            block_update = hyperparameters["block_update"];
+            hyperparameters.erase("block_update");
+        }
+        Classifier::setHyperparameters(hyperparameters);
+    }
+    void Boost::buildModel(const torch::Tensor& weights)
+    {
+        // Models shall be built in trainModel
+        models.clear();
+        significanceModels.clear();
+        n_models = 0;
+        // Prepare the validation dataset
+        auto y_ = dataset.index({ -1, "..." });
+        if (convergence) {
+            // Prepare train & validation sets from train data
+            auto fold = folding::StratifiedKFold(5, y_, 271);
+            auto [train, test] = fold.getFold(0);
+            auto train_t = torch::tensor(train);
+            auto test_t = torch::tensor(test);
+            // Get train and validation sets
+            X_train = dataset.index({ torch::indexing::Slice(0, dataset.size(0) - 1), train_t });
+            y_train = dataset.index({ -1, train_t });
+            X_test = dataset.index({ torch::indexing::Slice(0, dataset.size(0) - 1), test_t });
+            y_test = dataset.index({ -1, test_t });
+            dataset = X_train;
+            m = X_train.size(1);
+            auto n_classes = states.at(className).size();
+            // Build dataset with train data
+            buildDataset(y_train);
+            metrics = Metrics(dataset, features, className, n_classes);
+        } else {
+            // Use all data to train
+            X_train = dataset.index({ torch::indexing::Slice(0, dataset.size(0) - 1), "..." });
+            y_train = y_;
+        }
+    }
+    std::vector<int> Boost::featureSelection(torch::Tensor& weights_)
+    {
+        int maxFeatures = 0;
+        if (select_features_algorithm == SelectFeatures.CFS) {
+            featureSelector = new CFS(dataset, features, className, maxFeatures, states.at(className).size(), weights_);
+        } else if (select_features_algorithm == SelectFeatures.IWSS) {
+            if (threshold < 0 || threshold >0.5) {
+                throw std::invalid_argument("Invalid threshold value for " + SelectFeatures.IWSS + " [0, 0.5]");
+            }
+            featureSelector = new IWSS(dataset, features, className, maxFeatures, states.at(className).size(), weights_, threshold);
+        } else if (select_features_algorithm == SelectFeatures.FCBF) {
+            if (threshold < 1e-7 || threshold > 1) {
+                throw std::invalid_argument("Invalid threshold value for " + SelectFeatures.FCBF + " [1e-7, 1]");
+            }
+            featureSelector = new FCBF(dataset, features, className, maxFeatures, states.at(className).size(), weights_, threshold);
+        }
+        featureSelector->fit();
+        auto featuresUsed = featureSelector->getFeatures();
+        delete featureSelector;
+        return featuresUsed;
+    }
+    std::tuple<torch::Tensor&, double, bool> Boost::update_weights(torch::Tensor& ytrain, torch::Tensor& ypred, torch::Tensor& weights)
+    {
+        bool terminate = false;
+        double alpha_t = 0;
+        auto mask_wrong = ypred != ytrain;
+        auto mask_right = ypred == ytrain;
+        auto masked_weights = weights * mask_wrong.to(weights.dtype());
+        double epsilon_t = masked_weights.sum().item<double>();
+        if (epsilon_t > 0.5) {
+            // Inverse the weights policy (plot ln(wt))
+            // "In each round of AdaBoost, there is a sanity check to ensure that the current base 
+            // learner is better than random guess" (Zhi-Hua Zhou, 2012)
+            terminate = true;
+        } else {
+            double wt = (1 - epsilon_t) / epsilon_t;
+            alpha_t = epsilon_t == 0 ? 1 : 0.5 * log(wt);
+            // Step 3.2: Update weights for next classifier
+            // Step 3.2.1: Update weights of wrong samples
+            weights += mask_wrong.to(weights.dtype()) * exp(alpha_t) * weights;
+            // Step 3.2.2: Update weights of right samples
+            weights += mask_right.to(weights.dtype()) * exp(-alpha_t) * weights;
+            // Step 3.3: Normalise the weights
+            double totalWeights = torch::sum(weights).item<double>();
+            weights = weights / totalWeights;
+        }
+        return { weights, alpha_t, terminate };
+    }
+    std::tuple<torch::Tensor&, double, bool> Boost::update_weights_block(int k, torch::Tensor& ytrain, torch::Tensor& weights)
+    {
+        /* Update Block algorithm
+            k = # of models in block
+            n_models = # of models in ensemble to make predictions
+            n_models_bak = # models saved
+            models = vector of models to make predictions
+            models_bak = models not used to make predictions
+            significances_bak = backup of significances vector
+
+            Case list
+            A) k = 1, n_models = 1		=> n = 0 , n_models = n + k
+            B) k = 1, n_models = n + 1	=> n_models = n + k
+            C) k > 1, n_models = k + 1 	=> n= 1, n_models = n + k
+            D) k > 1, n_models = k		=> n = 0, n_models = n + k
+            E) k > 1, n_models = k + n	=> n_models = n + k
+
+            A, D) n=0, k > 0, n_models == k
+            1. n_models_bak <- n_models
+            2. significances_bak <- significances
+            3. significances = vector(k, 1)
+            4. Don’t move any classifiers out of models
+            5. n_models <- k
+            6. Make prediction, compute alpha, update weights
+            7. Don’t restore any classifiers to models
+            8. significances <- significances_bak
+            9. Update last k significances
+            10. n_models <- n_models_bak
+
+            B, C, E) n > 0, k > 0, n_models == n + k
+            1. n_models_bak <- n_models
+            2. significances_bak <- significances
+            3. significances = vector(k, 1)
+            4. Move first n classifiers to models_bak
+            5. n_models <- k
+            6. Make prediction, compute alpha, update weights
+            7. Insert classifiers in models_bak to be the first n models
+            8. significances <- significances_bak
+            9. Update last k significances
+            10. n_models <- n_models_bak
+        */
+        //
+        // Make predict with only the last k models
+        //
+        std::unique_ptr<Classifier> model;
+        std::vector<std::unique_ptr<Classifier>> models_bak;
+        // 1. n_models_bak <- n_models 2. significances_bak <- significances
+        auto significance_bak = significanceModels;
+        auto n_models_bak = n_models;
+        // 3. significances = vector(k, 1)
+        significanceModels = std::vector<double>(k, 1.0);
+        // 4. Move first n classifiers to models_bak
+        // backup the first n_models - k models (if n_models == k, don't backup any)
+        for (int i = 0; i < n_models - k; ++i) {
+            model = std::move(models[0]);
+            models.erase(models.begin());
+            models_bak.push_back(std::move(model));
+        }
+        assert(models.size() == k);
+        // 5. n_models <- k
+        n_models = k;
+        // 6. Make prediction, compute alpha, update weights
+        auto ypred = predict(X_train);
+        //
+        // Update weights
+        //
+        double alpha_t;
+        bool terminate;
+        std::tie(weights, alpha_t, terminate) = update_weights(y_train, ypred, weights);
+        //
+        // Restore the models if needed
+        //
+        // 7. Insert classifiers in models_bak to be the first n models
+        // if n_models_bak == k, don't restore any, because none of them were moved
+        if (k != n_models_bak) {
+            // Insert in the same order as they were extracted
+            int bak_size = models_bak.size();
+            for (int i = 0; i < bak_size; ++i) {
+                model = std::move(models_bak[bak_size - 1 - i]);
+                models_bak.erase(models_bak.end() - 1);
+                models.insert(models.begin(), std::move(model));
+            }
+        }
+        // 8. significances <- significances_bak
+        significanceModels = significance_bak;
+        //
+        // Update the significance of the last k models
+        //
+        // 9. Update last k significances
+        for (int i = 0; i < k; ++i) {
+            significanceModels[n_models_bak - k + i] = alpha_t;
+        }
+        // 10. n_models <- n_models_bak
+        n_models = n_models_bak;
+        return { weights, alpha_t, terminate };
+    }
+}

bayesnet/ensembles/Boost.h

@@ -0,0 +1,52 @@
+// ***************************************************************
+// SPDX-FileCopyrightText: Copyright 2024 Ricardo Montañana Gómez
+// SPDX-FileType: SOURCE
+// SPDX-License-Identifier: MIT
+// ***************************************************************
+
+#ifndef BOOST_H
+#define BOOST_H
+#include <string>
+#include <tuple>
+#include <vector>
+#include <nlohmann/json.hpp>
+#include <torch/torch.h>
+#include "Ensemble.h"
+#include "bayesnet/feature_selection/FeatureSelect.h"
+namespace bayesnet {
+    const struct {
+        std::string CFS = "CFS";
+        std::string FCBF = "FCBF";
+        std::string IWSS = "IWSS";
+    }SelectFeatures;
+    const struct {
+        std::string ASC = "asc";
+        std::string DESC = "desc";
+        std::string RAND = "rand";
+    }Orders;
+    class Boost : public Ensemble {
+    public:
+        explicit Boost(bool predict_voting = false);
+        virtual ~Boost() = default;
+        void setHyperparameters(const nlohmann::json& hyperparameters_) override;
+    protected:
+        std::vector<int> featureSelection(torch::Tensor& weights_);
+        void buildModel(const torch::Tensor& weights) override;
+        std::tuple<torch::Tensor&, double, bool> update_weights(torch::Tensor& ytrain, torch::Tensor& ypred, torch::Tensor& weights);
+        std::tuple<torch::Tensor&, double, bool> update_weights_block(int k, torch::Tensor& ytrain, torch::Tensor& weights);
+        torch::Tensor X_train, y_train, X_test, y_test;
+        // Hyperparameters
+        bool bisection = true; // if true, use bisection stratety to add k models at once to the ensemble
+        int maxTolerance = 3;
+        std::string order_algorithm; // order to process the KBest features asc, desc, rand
+        bool convergence = true; //if true, stop when the model does not improve
+        bool convergence_best = false; // wether to keep the best accuracy to the moment or the last accuracy as prior accuracy
+        bool selectFeatures = false; // if true, use feature selection
+        std::string select_features_algorithm = Orders.DESC; // Selected feature selection algorithm
+        FeatureSelect* featureSelector = nullptr;
+        double threshold = -1;
+        bool block_update = false;
+
+    };
+}
+#endif

bayesnet/ensembles/BoostA2DE.cc

@@ -0,0 +1,167 @@
+// ***************************************************************
+// SPDX-FileCopyrightText: Copyright 2024 Ricardo Montañana Gómez
+// SPDX-FileType: SOURCE
+// SPDX-License-Identifier: MIT
+// ***************************************************************
+
+#include <set>
+#include <functional>
+#include <limits.h>
+#include <tuple>
+#include <folding.hpp>
+#include "bayesnet/feature_selection/CFS.h"
+#include "bayesnet/feature_selection/FCBF.h"
+#include "bayesnet/feature_selection/IWSS.h"
+#include "BoostA2DE.h"
+
+namespace bayesnet {
+
+    BoostA2DE::BoostA2DE(bool predict_voting) : Boost(predict_voting)
+    {
+    }
+    std::vector<int> BoostA2DE::initializeModels()
+    {
+        torch::Tensor weights_ = torch::full({ m }, 1.0 / m, torch::kFloat64);
+        std::vector<int> featuresSelected = featureSelection(weights_);
+        if (featuresSelected.size() < 2) {
+            notes.push_back("No features selected in initialization");
+            status = ERROR;
+            return std::vector<int>();
+        }
+        for (int i = 0; i < featuresSelected.size() - 1; i++) {
+            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_);
+                models.push_back(std::move(model));
+                significanceModels.push_back(1.0); // They will be updated later in trainModel
+                n_models++;
+            }
+        }
+        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)
+    {
+        //
+        // Logging setup
+        //
+        // loguru::set_thread_name("BoostA2DE");
+        // loguru::g_stderr_verbosity = loguru::Verbosity_OFF;
+        // loguru::add_file("boostA2DE.log", loguru::Truncate, loguru::Verbosity_MAX);
+
+        // Algorithm based on the adaboost algorithm for classification
+        // as explained in Ensemble methods (Zhi-Hua Zhou, 2012)
+        fitted = true;
+        double alpha_t = 0;
+        torch::Tensor weights_ = torch::full({ m }, 1.0 / m, torch::kFloat64);
+        bool finished = false;
+        std::vector<int> featuresUsed;
+        if (selectFeatures) {
+            featuresUsed = initializeModels();
+            auto ypred = predict(X_train);
+            std::tie(weights_, alpha_t, finished) = update_weights(y_train, ypred, weights_);
+            // Update significance of the models
+            for (int i = 0; i < n_models; ++i) {
+                significanceModels[i] = alpha_t;
+            }
+            if (finished) {
+                return;
+            }
+        }
+        int numItemsPack = 0; // The counter of the models inserted in the current pack
+        // Variables to control the accuracy finish condition
+        double priorAccuracy = 0.0;
+        double improvement = 1.0;
+        double convergence_threshold = 1e-4;
+        int tolerance = 0; // number of times the accuracy is lower than the convergence_threshold
+        // Step 0: Set the finish condition
+        // epsilon sub t > 0.5 => inverse the weights policy
+        // validation error is not decreasing
+        // run out of features
+        bool ascending = order_algorithm == Orders.ASC;
+        std::mt19937 g{ 173 };
+        std::vector<std::pair<int, int>> pairSelection;
+        while (!finished) {
+            // Step 1: Build ranking with mutual information
+            pairSelection = metrics.SelectKPairs(weights_, featuresUsed, ascending, 0); // Get all the pairs sorted
+            if (order_algorithm == Orders.RAND) {
+                std::shuffle(pairSelection.begin(), pairSelection.end(), g);
+            }
+            int k = bisection ? pow(2, tolerance) : 1;
+            int counter = 0; // The model counter of the current pack
+            // VLOG_SCOPE_F(1, "counter=%d k=%d featureSelection.size: %zu", counter, k, featureSelection.size());
+            while (counter++ < k && pairSelection.size() > 0) {
+                auto feature_pair = pairSelection[0];
+                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_);
+                alpha_t = 0.0;
+                if (!block_update) {
+                    auto ypred = model->predict(X_train);
+                    // Step 3.1: Compute the classifier amout of say
+                    std::tie(weights_, alpha_t, finished) = update_weights(y_train, ypred, weights_);
+                }
+                // Step 3.4: Store classifier and its accuracy to weigh its future vote
+                numItemsPack++;
+                models.push_back(std::move(model));
+                significanceModels.push_back(alpha_t);
+                n_models++;
+                // VLOG_SCOPE_F(2, "numItemsPack: %d n_models: %d featuresUsed: %zu", numItemsPack, n_models, featuresUsed.size());
+            }
+            if (block_update) {
+                std::tie(weights_, alpha_t, finished) = update_weights_block(k, y_train, weights_);
+            }
+            if (convergence && !finished) {
+                auto y_val_predict = predict(X_test);
+                double accuracy = (y_val_predict == y_test).sum().item<double>() / (double)y_test.size(0);
+                if (priorAccuracy == 0) {
+                    priorAccuracy = accuracy;
+                } else {
+                    improvement = accuracy - priorAccuracy;
+                }
+                if (improvement < convergence_threshold) {
+                    // VLOG_SCOPE_F(3, "  (improvement<threshold) tolerance: %d numItemsPack: %d improvement: %f prior: %f current: %f", tolerance, numItemsPack, improvement, priorAccuracy, accuracy);
+                    tolerance++;
+                } else {
+                    // VLOG_SCOPE_F(3, "* (improvement>=threshold) Reset. tolerance: %d numItemsPack: %d improvement: %f prior: %f current: %f", tolerance, numItemsPack, improvement, priorAccuracy, accuracy);
+                    tolerance = 0; // Reset the counter if the model performs better
+                    numItemsPack = 0;
+                }
+                if (convergence_best) {
+                    // Keep the best accuracy until now as the prior accuracy
+                    priorAccuracy = std::max(accuracy, priorAccuracy);
+                } else {
+                    // Keep the last accuray obtained as the prior accuracy
+                    priorAccuracy = accuracy;
+                }
+            }
+            // VLOG_SCOPE_F(1, "tolerance: %d featuresUsed.size: %zu features.size: %zu", tolerance, featuresUsed.size(), features.size());
+            finished = finished || tolerance > maxTolerance || pairSelection.size() == 0;
+        }
+        if (tolerance > maxTolerance) {
+            if (numItemsPack < n_models) {
+                notes.push_back("Convergence threshold reached & " + std::to_string(numItemsPack) + " models eliminated");
+                // VLOG_SCOPE_F(4, "Convergence threshold reached & %d models eliminated of %d", numItemsPack, n_models);
+                for (int i = 0; i < numItemsPack; ++i) {
+                    significanceModels.pop_back();
+                    models.pop_back();
+                    n_models--;
+                }
+            } else {
+                notes.push_back("Convergence threshold reached & 0 models eliminated");
+                // VLOG_SCOPE_F(4, "Convergence threshold reached & 0 models eliminated n_models=%d numItemsPack=%d", n_models, numItemsPack);
+            }
+        }
+        if (pairSelection.size() > 0) {
+            notes.push_back("Pairs not used in train: " + std::to_string(pairSelection.size()));
+            status = WARNING;
+        }
+        notes.push_back("Number of models: " + std::to_string(n_models));
+    }
+    std::vector<std::string> BoostA2DE::graph(const std::string& title) const
+    {
+        return Ensemble::graph(title);
+    }
+}

bayesnet/ensembles/BoostA2DE.h

@@ -0,0 +1,25 @@
+// ***************************************************************
+// SPDX-FileCopyrightText: Copyright 2024 Ricardo Montañana Gómez
+// SPDX-FileType: SOURCE
+// SPDX-License-Identifier: MIT
+// ***************************************************************
+
+#ifndef BOOSTA2DE_H
+#define BOOSTA2DE_H
+#include <string>
+#include <vector>
+#include "bayesnet/classifiers/SPnDE.h"
+#include "Boost.h"
+namespace bayesnet {
+    class BoostA2DE : public Boost {
+    public:
+        explicit BoostA2DE(bool predict_voting = false);
+        virtual ~BoostA2DE() = default;
+        std::vector<std::string> graph(const std::string& title = "BoostA2DE") const override;
+    protected:
+        void trainModel(const torch::Tensor& weights) override;
+    private:
+        std::vector<int> initializeModels();
+    };
+}
+#endif

bayesnet/ensembles/BoostAODE.cc

@@ -4,275 +4,40 @@
 // SPDX-License-Identifier: MIT
 // ***************************************************************
 
+#include <random> 
 #include <set>
 #include <functional>
 #include <limits.h>
 #include <tuple>
-#include <folding.hpp>
-#include "bayesnet/feature_selection/CFS.h"
-#include "bayesnet/feature_selection/FCBF.h"
-#include "bayesnet/feature_selection/IWSS.h"
 #include "BoostAODE.h"
-#include "lib/log/loguru.cpp"
 
 namespace bayesnet {
 
-    BoostAODE::BoostAODE(bool predict_voting) : Ensemble(predict_voting)
+    BoostAODE::BoostAODE(bool predict_voting) : Boost(predict_voting)
     {
-        validHyperparameters = {
-            "maxModels", "bisection", "order", "convergence", "convergence_best", "threshold",
-            "select_features", "maxTolerance", "predict_voting", "block_update"
-        };
-
-    }
-    void BoostAODE::buildModel(const torch::Tensor& weights)
-    {
-        // Models shall be built in trainModel
-        models.clear();
-        significanceModels.clear();
-        n_models = 0;
-        // Prepare the validation dataset
-        auto y_ = dataset.index({ -1, "..." });
-        if (convergence) {
-            // Prepare train & validation sets from train data
-            auto fold = folding::StratifiedKFold(5, y_, 271);
-            auto [train, test] = fold.getFold(0);
-            auto train_t = torch::tensor(train);
-            auto test_t = torch::tensor(test);
-            // Get train and validation sets
-            X_train = dataset.index({ torch::indexing::Slice(0, dataset.size(0) - 1), train_t });
-            y_train = dataset.index({ -1, train_t });
-            X_test = dataset.index({ torch::indexing::Slice(0, dataset.size(0) - 1), test_t });
-            y_test = dataset.index({ -1, test_t });
-            dataset = X_train;
-            m = X_train.size(1);
-            auto n_classes = states.at(className).size();
-            // Build dataset with train data
-            buildDataset(y_train);
-            metrics = Metrics(dataset, features, className, n_classes);
-        } else {
-            // Use all data to train
-            X_train = dataset.index({ torch::indexing::Slice(0, dataset.size(0) - 1), "..." });
-            y_train = y_;
-        }
-    }
-    void BoostAODE::setHyperparameters(const nlohmann::json& hyperparameters_)
-    {
-        auto hyperparameters = hyperparameters_;
-        if (hyperparameters.contains("order")) {
-            std::vector<std::string> algos = { Orders.ASC, Orders.DESC, Orders.RAND };
-            order_algorithm = hyperparameters["order"];
-            if (std::find(algos.begin(), algos.end(), order_algorithm) == algos.end()) {
-                throw std::invalid_argument("Invalid order algorithm, valid values [" + Orders.ASC + ", " + Orders.DESC + ", " + Orders.RAND + "]");
-            }
-            hyperparameters.erase("order");
-        }
-        if (hyperparameters.contains("convergence")) {
-            convergence = hyperparameters["convergence"];
-            hyperparameters.erase("convergence");
-        }
-        if (hyperparameters.contains("convergence_best")) {
-            convergence_best = hyperparameters["convergence_best"];
-            hyperparameters.erase("convergence_best");
-        }
-        if (hyperparameters.contains("bisection")) {
-            bisection = hyperparameters["bisection"];
-            hyperparameters.erase("bisection");
-        }
-        if (hyperparameters.contains("threshold")) {
-            threshold = hyperparameters["threshold"];
-            hyperparameters.erase("threshold");
-        }
-        if (hyperparameters.contains("maxTolerance")) {
-            maxTolerance = hyperparameters["maxTolerance"];
-            if (maxTolerance < 1 || maxTolerance > 4)
-                throw std::invalid_argument("Invalid maxTolerance value, must be greater in [1, 4]");
-            hyperparameters.erase("maxTolerance");
-        }
-        if (hyperparameters.contains("predict_voting")) {
-            predict_voting = hyperparameters["predict_voting"];
-            hyperparameters.erase("predict_voting");
-        }
-        if (hyperparameters.contains("select_features")) {
-            auto selectedAlgorithm = hyperparameters["select_features"];
-            std::vector<std::string> algos = { SelectFeatures.IWSS, SelectFeatures.CFS, SelectFeatures.FCBF };
-            selectFeatures = true;
-            select_features_algorithm = selectedAlgorithm;
-            if (std::find(algos.begin(), algos.end(), selectedAlgorithm) == algos.end()) {
-                throw std::invalid_argument("Invalid selectFeatures value, valid values [" + SelectFeatures.IWSS + ", " + SelectFeatures.CFS + ", " + SelectFeatures.FCBF + "]");
-            }
-            hyperparameters.erase("select_features");
-        }
-        if (hyperparameters.contains("block_update")) {
-            block_update = hyperparameters["block_update"];
-            hyperparameters.erase("block_update");
-        }
-        Classifier::setHyperparameters(hyperparameters);
-    }
-    std::tuple<torch::Tensor&, double, bool> update_weights(torch::Tensor& ytrain, torch::Tensor& ypred, torch::Tensor& weights)
-    {
-        bool terminate = false;
-        double alpha_t = 0;
-        auto mask_wrong = ypred != ytrain;
-        auto mask_right = ypred == ytrain;
-        auto masked_weights = weights * mask_wrong.to(weights.dtype());
-        double epsilon_t = masked_weights.sum().item<double>();
-        if (epsilon_t > 0.5) {
-            // Inverse the weights policy (plot ln(wt))
-            // "In each round of AdaBoost, there is a sanity check to ensure that the current base 
-            // learner is better than random guess" (Zhi-Hua Zhou, 2012)
-            terminate = true;
-        } else {
-            double wt = (1 - epsilon_t) / epsilon_t;
-            alpha_t = epsilon_t == 0 ? 1 : 0.5 * log(wt);
-            // Step 3.2: Update weights for next classifier
-            // Step 3.2.1: Update weights of wrong samples
-            weights += mask_wrong.to(weights.dtype()) * exp(alpha_t) * weights;
-            // Step 3.2.2: Update weights of right samples
-            weights += mask_right.to(weights.dtype()) * exp(-alpha_t) * weights;
-            // Step 3.3: Normalise the weights
-            double totalWeights = torch::sum(weights).item<double>();
-            weights = weights / totalWeights;
-        }
-        return { weights, alpha_t, terminate };
-    }
-    std::tuple<torch::Tensor&, double, bool> BoostAODE::update_weights_block(int k, torch::Tensor& ytrain, torch::Tensor& weights)
-    {
-        /* Update Block algorithm
-            k = # of models in block
-            n_models = # of models in ensemble to make predictions
-            n_models_bak = # models saved
-            models = vector of models to make predictions
-            models_bak = models not used to make predictions
-            significances_bak = backup of significances vector
-
-            Case list
-            A) k = 1, n_models = 1		=> n = 0 , n_models = n + k
-            B) k = 1, n_models = n + 1	=> n_models = n + k
-            C) k > 1, n_models = k + 1 	=> n= 1, n_models = n + k
-            D) k > 1, n_models = k		=> n = 0, n_models = n + k
-            E) k > 1, n_models = k + n	=> n_models = n + k
-
-            A, D) n=0, k > 0, n_models == k
-            1. n_models_bak <- n_models
-            2. significances_bak <- significances
-            3. significances = vector(k, 1)
-            4. Don’t move any classifiers out of models
-            5. n_models <- k
-            6. Make prediction, compute alpha, update weights
-            7. Don’t restore any classifiers to models
-            8. significances <- significances_bak
-            9. Update last k significances
-            10. n_models <- n_models_bak
-
-            B, C, E) n > 0, k > 0, n_models == n + k
-            1. n_models_bak <- n_models
-            2. significances_bak <- significances
-            3. significances = vector(k, 1)
-            4. Move first n classifiers to models_bak
-            5. n_models <- k
-            6. Make prediction, compute alpha, update weights
-            7. Insert classifiers in models_bak to be the first n models
-            8. significances <- significances_bak
-            9. Update last k significances
-            10. n_models <- n_models_bak
-        */
-        //
-        // Make predict with only the last k models
-        //
-        std::unique_ptr<Classifier> model;
-        std::vector<std::unique_ptr<Classifier>> models_bak;
-        // 1. n_models_bak <- n_models 2. significances_bak <- significances
-        auto significance_bak = significanceModels;
-        auto n_models_bak = n_models;
-        // 3. significances = vector(k, 1)
-        significanceModels = std::vector<double>(k, 1.0);
-        // 4. Move first n classifiers to models_bak
-        // backup the first n_models - k models (if n_models == k, don't backup any)
-        for (int i = 0; i < n_models - k; ++i) {
-            model = std::move(models[0]);
-            models.erase(models.begin());
-            models_bak.push_back(std::move(model));
-        }
-        assert(models.size() == k);
-        // 5. n_models <- k
-        n_models = k;
-        // 6. Make prediction, compute alpha, update weights
-        auto ypred = predict(X_train);
-        //
-        // Update weights
-        //
-        double alpha_t;
-        bool terminate;
-        std::tie(weights, alpha_t, terminate) = update_weights(y_train, ypred, weights);
-        //
-        // Restore the models if needed
-        //
-        // 7. Insert classifiers in models_bak to be the first n models
-        // if n_models_bak == k, don't restore any, because none of them were moved
-        if (k != n_models_bak) {
-            // Insert in the same order as they were extracted
-            int bak_size = models_bak.size();
-            for (int i = 0; i < bak_size; ++i) {
-                model = std::move(models_bak[bak_size - 1 - i]);
-                models_bak.erase(models_bak.end() - 1);
-                models.insert(models.begin(), std::move(model));
-            }
-        }
-        // 8. significances <- significances_bak
-        significanceModels = significance_bak;
-        //
-        // Update the significance of the last k models
-        //
-        // 9. Update last k significances
-        for (int i = 0; i < k; ++i) {
-            significanceModels[n_models_bak - k + i] = alpha_t;
-        }
-        // 10. n_models <- n_models_bak
-        n_models = n_models_bak;
-        return { weights, alpha_t, terminate };
     }
     std::vector<int> BoostAODE::initializeModels()
     {
-        std::vector<int> featuresUsed;
         torch::Tensor weights_ = torch::full({ m }, 1.0 / m, torch::kFloat64);
-        int maxFeatures = 0;
-        if (select_features_algorithm == SelectFeatures.CFS) {
-            featureSelector = new CFS(dataset, features, className, maxFeatures, states.at(className).size(), weights_);
-        } else if (select_features_algorithm == SelectFeatures.IWSS) {
-            if (threshold < 0 || threshold >0.5) {
-                throw std::invalid_argument("Invalid threshold value for " + SelectFeatures.IWSS + " [0, 0.5]");
-            }
-            featureSelector = new IWSS(dataset, features, className, maxFeatures, states.at(className).size(), weights_, threshold);
-        } else if (select_features_algorithm == SelectFeatures.FCBF) {
-            if (threshold < 1e-7 || threshold > 1) {
-                throw std::invalid_argument("Invalid threshold value for " + SelectFeatures.FCBF + " [1e-7, 1]");
-            }
-            featureSelector = new FCBF(dataset, features, className, maxFeatures, states.at(className).size(), weights_, threshold);
-        }
-        featureSelector->fit();
-        auto cfsFeatures = featureSelector->getFeatures();
-        auto scores = featureSelector->getScores();
-        for (const int& feature : cfsFeatures) {
-            featuresUsed.push_back(feature);
+        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_);
             models.push_back(std::move(model));
             significanceModels.push_back(1.0); // They will be updated later in trainModel
             n_models++;
         }
-        notes.push_back("Used features in initialization: " + std::to_string(featuresUsed.size()) + " of " + std::to_string(features.size()) + " with " + select_features_algorithm);
-        delete featureSelector;
-        return featuresUsed;
+        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)
     {
         //
         // Logging setup
         //
-        loguru::set_thread_name("BoostAODE");
-        loguru::g_stderr_verbosity = loguru::Verbosity_OFF;
-        loguru::add_file("boostAODE.log", loguru::Truncate, loguru::Verbosity_MAX);
+        // loguru::set_thread_name("BoostAODE");
+        // loguru::g_stderr_verbosity = loguru::Verbosity_OFF;
+        // loguru::add_file("boostAODE.log", loguru::Truncate, loguru::Verbosity_MAX);
 
         // Algorithm based on the adaboost algorithm for classification
         // as explained in Ensemble methods (Zhi-Hua Zhou, 2012)
@@ -318,7 +83,7 @@ namespace bayesnet {
             );
             int k = bisection ? pow(2, tolerance) : 1;
             int counter = 0; // The model counter of the current pack
-            VLOG_SCOPE_F(1, "counter=%d k=%d featureSelection.size: %zu", counter, k, featureSelection.size());
+            // VLOG_SCOPE_F(1, "counter=%d k=%d featureSelection.size: %zu", counter, k, featureSelection.size());
             while (counter++ < k && featureSelection.size() > 0) {
                 auto feature = featureSelection[0];
                 featureSelection.erase(featureSelection.begin());
@@ -337,7 +102,7 @@ namespace bayesnet {
                 models.push_back(std::move(model));
                 significanceModels.push_back(alpha_t);
                 n_models++;
-                VLOG_SCOPE_F(2, "numItemsPack: %d n_models: %d featuresUsed: %zu", numItemsPack, n_models, featuresUsed.size());
+                // VLOG_SCOPE_F(2, "numItemsPack: %d n_models: %d featuresUsed: %zu", numItemsPack, n_models, featuresUsed.size());
             }
             if (block_update) {
                 std::tie(weights_, alpha_t, finished) = update_weights_block(k, y_train, weights_);
@@ -351,10 +116,10 @@ namespace bayesnet {
                     improvement = accuracy - priorAccuracy;
                 }
                 if (improvement < convergence_threshold) {
-                    VLOG_SCOPE_F(3, "  (improvement<threshold) tolerance: %d numItemsPack: %d improvement: %f prior: %f current: %f", tolerance, numItemsPack, improvement, priorAccuracy, accuracy);
+                    // VLOG_SCOPE_F(3, "  (improvement<threshold) tolerance: %d numItemsPack: %d improvement: %f prior: %f current: %f", tolerance, numItemsPack, improvement, priorAccuracy, accuracy);
                     tolerance++;
                 } else {
-                    VLOG_SCOPE_F(3, "* (improvement>=threshold) Reset. tolerance: %d numItemsPack: %d improvement: %f prior: %f current: %f", tolerance, numItemsPack, improvement, priorAccuracy, accuracy);
+                    // VLOG_SCOPE_F(3, "* (improvement>=threshold) Reset. tolerance: %d numItemsPack: %d improvement: %f prior: %f current: %f", tolerance, numItemsPack, improvement, priorAccuracy, accuracy);
                     tolerance = 0; // Reset the counter if the model performs better
                     numItemsPack = 0;
                 }
@@ -366,13 +131,13 @@ namespace bayesnet {
                     priorAccuracy = accuracy;
                 }
             }
-            VLOG_SCOPE_F(1, "tolerance: %d featuresUsed.size: %zu features.size: %zu", tolerance, featuresUsed.size(), features.size());
+            // VLOG_SCOPE_F(1, "tolerance: %d featuresUsed.size: %zu features.size: %zu", tolerance, featuresUsed.size(), features.size());
             finished = finished || tolerance > maxTolerance || featuresUsed.size() == features.size();
         }
         if (tolerance > maxTolerance) {
             if (numItemsPack < n_models) {
                 notes.push_back("Convergence threshold reached & " + std::to_string(numItemsPack) + " models eliminated");
-                VLOG_SCOPE_F(4, "Convergence threshold reached & %d models eliminated of %d", numItemsPack, n_models);
+                // VLOG_SCOPE_F(4, "Convergence threshold reached & %d models eliminated of %d", numItemsPack, n_models);
                 for (int i = 0; i < numItemsPack; ++i) {
                     significanceModels.pop_back();
                     models.pop_back();
@@ -380,7 +145,7 @@ namespace bayesnet {
                 }
             } else {
                 notes.push_back("Convergence threshold reached & 0 models eliminated");
-                VLOG_SCOPE_F(4, "Convergence threshold reached & 0 models eliminated n_models=%d numItemsPack=%d", n_models, numItemsPack);
+                // VLOG_SCOPE_F(4, "Convergence threshold reached & 0 models eliminated n_models=%d numItemsPack=%d", n_models, numItemsPack);
             }
         }
         if (featuresUsed.size() != features.size()) {

bayesnet/ensembles/BoostAODE.h

@@ -6,45 +6,21 @@
 
 #ifndef BOOSTAODE_H
 #define BOOSTAODE_H
-#include <map>
+#include <string>
+#include <vector>
 #include "bayesnet/classifiers/SPODE.h"
-#include "bayesnet/feature_selection/FeatureSelect.h"
-#include "Ensemble.h"
+#include "Boost.h"
+
 namespace bayesnet {
-    const struct {
-        std::string CFS = "CFS";
-        std::string FCBF = "FCBF";
-        std::string IWSS = "IWSS";
-    }SelectFeatures;
-    const struct {
-        std::string ASC = "asc";
-        std::string DESC = "desc";
-        std::string RAND = "rand";
-    }Orders;
-    class BoostAODE : public Ensemble {
+    class BoostAODE : public Boost {
     public:
         explicit BoostAODE(bool predict_voting = false);
         virtual ~BoostAODE() = default;
         std::vector<std::string> graph(const std::string& title = "BoostAODE") const override;
-        void setHyperparameters(const nlohmann::json& hyperparameters_) override;
     protected:
-        void buildModel(const torch::Tensor& weights) override;
         void trainModel(const torch::Tensor& weights) override;
     private:
-        std::tuple<torch::Tensor&, double, bool> update_weights_block(int k, torch::Tensor& ytrain, torch::Tensor& weights);
         std::vector<int> initializeModels();
-        torch::Tensor X_train, y_train, X_test, y_test;
-        // Hyperparameters
-        bool bisection = true; // if true, use bisection stratety to add k models at once to the ensemble
-        int maxTolerance = 3;
-        std::string order_algorithm; // order to process the KBest features asc, desc, rand
-        bool convergence = true; //if true, stop when the model does not improve
-        bool convergence_best = false; // wether to keep the best accuracy to the moment or the last accuracy as prior accuracy
-        bool selectFeatures = false; // if true, use feature selection
-        std::string select_features_algorithm = Orders.DESC; // Selected feature selection algorithm
-        FeatureSelect* featureSelector = nullptr;
-        double threshold = -1;
-        bool block_update = false;
     };
 }
 #endif

bayesnet/utils/BayesMetrics.cc

@@ -30,6 +30,53 @@ namespace bayesnet {
         }
         samples.index_put_({ -1, "..." }, torch::tensor(labels, torch::kInt32));
     }
+    std::vector<std::pair<int, int>> Metrics::SelectKPairs(const torch::Tensor& weights, std::vector<int>& featuresExcluded, bool ascending, unsigned k)
+    {
+        // Return the K Best features 
+        auto n = features.size();
+        // compute scores
+        scoresKPairs.clear();
+        pairsKBest.clear();
+        auto labels = samples.index({ -1, "..." });
+        for (int i = 0; i < n - 1; ++i) {
+            if (std::find(featuresExcluded.begin(), featuresExcluded.end(), i) != featuresExcluded.end()) {
+                continue;
+            }
+            for (int j = i + 1; j < n; ++j) {
+                if (std::find(featuresExcluded.begin(), featuresExcluded.end(), j) != featuresExcluded.end()) {
+                    continue;
+                }
+                auto key = std::make_pair(i, j);
+                auto value = conditionalMutualInformation(samples.index({ i, "..." }), samples.index({ j, "..." }), labels, weights);
+                scoresKPairs.push_back({ key, value });
+            }
+        }
+        // sort scores
+        if (ascending) {
+            sort(scoresKPairs.begin(), scoresKPairs.end(), [](auto& a, auto& b)
+                { return a.second < b.second; });
+
+        } else {
+            sort(scoresKPairs.begin(), scoresKPairs.end(), [](auto& a, auto& b)
+                { return a.second > b.second; });
+        }
+        for (auto& [pairs, score] : scoresKPairs) {
+            pairsKBest.push_back(pairs);
+        }
+        if (k != 0 && k < pairsKBest.size()) {
+            if (ascending) {
+                int limit = pairsKBest.size() - k;
+                for (int i = 0; i < limit; i++) {
+                    pairsKBest.erase(pairsKBest.begin());
+                    scoresKPairs.erase(scoresKPairs.begin());
+                }
+            } else {
+                pairsKBest.resize(k);
+                scoresKPairs.resize(k);
+            }
+        }
+        return pairsKBest;
+    }
     std::vector<int> Metrics::SelectKBestWeighted(const torch::Tensor& weights, bool ascending, unsigned k)
     {
         // Return the K Best features 
@@ -69,7 +116,10 @@ namespace bayesnet {
     {
         return scoresKBest;
     }
-
+    std::vector<std::pair<std::pair<int, int>, double>> Metrics::getScoresKPairs() const
+    {
+        return scoresKPairs;
+    }
     torch::Tensor Metrics::conditionalEdge(const torch::Tensor& weights)
     {
         auto result = std::vector<double>();
@@ -148,24 +198,20 @@ namespace bayesnet {
         }
         return entropyValue;
     }
-    // H(Y|X,C) = sum_{x in X, c in C} p(x,c) H(Y|X=x,C=c)
+    // H(X|Y,C) = sum_{y in Y, c in C} p(x,c) H(X|Y=y,C=c)
     double Metrics::conditionalEntropy(const torch::Tensor& firstFeature, const torch::Tensor& secondFeature, const torch::Tensor& labels, const torch::Tensor& weights)
     {
         // Ensure the tensors are of the same length
         assert(firstFeature.size(0) == secondFeature.size(0) && firstFeature.size(0) == labels.size(0) && firstFeature.size(0) == weights.size(0));
-
         // Convert tensors to vectors for easier processing
         auto firstFeatureData = firstFeature.accessor<int, 1>();
         auto secondFeatureData = secondFeature.accessor<int, 1>();
         auto labelsData = labels.accessor<int, 1>();
         auto weightsData = weights.accessor<double, 1>();
-
         int numSamples = firstFeature.size(0);
-
         // Maps for joint and marginal probabilities
         std::map<std::tuple<int, int, int>, double> jointCount;
         std::map<std::tuple<int, int>, double> marginalCount;
-
         // Compute joint and marginal counts
         for (int i = 0; i < numSamples; ++i) {
             auto keyJoint = std::make_tuple(firstFeatureData[i], labelsData[i], secondFeatureData[i]);
@@ -174,34 +220,29 @@ namespace bayesnet {
             jointCount[keyJoint] += weightsData[i];
             marginalCount[keyMarginal] += weightsData[i];
         }
-
         // Total weight sum
         double totalWeight = torch::sum(weights).item<double>();
         if (totalWeight == 0)
             return 0;
-
         // Compute the conditional entropy
         double conditionalEntropy = 0.0;
-
         for (const auto& [keyJoint, jointFreq] : jointCount) {
             auto [x, c, y] = keyJoint;
             auto keyMarginal = std::make_tuple(x, c);
-
-            double p_xc = marginalCount[keyMarginal] / totalWeight;
+            //double p_xc = marginalCount[keyMarginal] / totalWeight;
             double p_y_given_xc = jointFreq / marginalCount[keyMarginal];
-
             if (p_y_given_xc > 0) {
                 conditionalEntropy -= (jointFreq / totalWeight) * std::log(p_y_given_xc);
             }
         }
         return conditionalEntropy;
     }
-    // I(X;Y) = H(Y) - H(Y|X)
+    // I(X;Y) = H(Y) - H(Y|X) ; I(X;Y) >= 0
     double Metrics::mutualInformation(const torch::Tensor& firstFeature, const torch::Tensor& secondFeature, const torch::Tensor& weights)
     {
-        return entropy(firstFeature, weights) - conditionalEntropy(firstFeature, secondFeature, weights);
+        return std::max(entropy(firstFeature, weights) - conditionalEntropy(firstFeature, secondFeature, weights), 0.0);
     }
-    // I(X;Y|C) = H(Y|C) - H(Y|X,C)
+    // I(X;Y|C) = H(X|C) - H(X|Y,C) >= 0
     double Metrics::conditionalMutualInformation(const torch::Tensor& firstFeature, const torch::Tensor& secondFeature, const torch::Tensor& labels, const torch::Tensor& weights)
     {
         return std::max(conditionalEntropy(firstFeature, labels, weights) - conditionalEntropy(firstFeature, secondFeature, labels, weights), 0.0);

bayesnet/utils/BayesMetrics.h

@@ -16,7 +16,9 @@ namespace bayesnet {
         Metrics(const torch::Tensor& samples, const std::vector<std::string>& features, const std::string& className, const int classNumStates);
         Metrics(const std::vector<std::vector<int>>& vsamples, const std::vector<int>& labels, const std::vector<std::string>& features, const std::string& className, const int classNumStates);
         std::vector<int> SelectKBestWeighted(const torch::Tensor& weights, bool ascending = false, unsigned k = 0);
+        std::vector<std::pair<int, int>> SelectKPairs(const torch::Tensor& weights, std::vector<int>& featuresExcluded, bool ascending = false, unsigned k = 0);
         std::vector<double> getScoresKBest() const;
+        std::vector<std::pair<std::pair<int, int>, double>> getScoresKPairs() const;
         double mutualInformation(const torch::Tensor& firstFeature, const torch::Tensor& secondFeature, const torch::Tensor& weights);
         double conditionalMutualInformation(const torch::Tensor& firstFeature, const torch::Tensor& secondFeature, const torch::Tensor& labels, const torch::Tensor& weights);
         torch::Tensor conditionalEdge(const torch::Tensor& weights);
@@ -33,7 +35,7 @@ namespace bayesnet {
         std::vector<std::pair<T, T>> doCombinations(const std::vector<T>& source)
         {
             std::vector<std::pair<T, T>> result;
-            for (int i = 0; i < source.size(); ++i) {
+            for (int i = 0; i < source.size() - 1; ++i) {
                 T temp = source[i];
                 for (int j = i + 1; j < source.size(); ++j) {
                     result.push_back({ temp, source[j] });
@@ -52,6 +54,8 @@ namespace bayesnet {
         int classNumStates = 0;
         std::vector<double> scoresKBest;
         std::vector<int> featuresKBest; // sorted indices of the features
+        std::vector<std::pair<int, int>> pairsKBest; // sorted indices of the pairs
+        std::vector<std::pair<std::pair<int, int>, double>> scoresKPairs;
         double conditionalEntropy(const torch::Tensor& firstFeature, const torch::Tensor& secondFeature, const torch::Tensor& weights);
     };
 }

docs/BoostAODE.md

@@ -27,4 +27,4 @@ The hyperparameters defined in the algorithm are:
 
 ## Operation
 
-### [Algorithm](./algorithm.md)
+### [Base Algorithm](./algorithm.md)

lib/Files/ArffFiles.cc

@@ -1,171 +0,0 @@
-#include "ArffFiles.h"
-#include <fstream>
-#include <sstream>
-#include <map>
-#include <iostream>
-
-std::vector<std::string> ArffFiles::getLines() const
-{
-    return lines;
-}
-
-unsigned long int ArffFiles::getSize() const
-{
-    return lines.size();
-}
-
-std::vector<std::pair<std::string, std::string>> ArffFiles::getAttributes() const
-{
-    return attributes;
-}
-
-std::string ArffFiles::getClassName() const
-{
-    return className;
-}
-
-std::string ArffFiles::getClassType() const
-{
-    return classType;
-}
-
-std::vector<std::vector<float>>& ArffFiles::getX()
-{
-    return X;
-}
-
-std::vector<int>& ArffFiles::getY()
-{
-    return y;
-}
-
-void ArffFiles::loadCommon(const std::string& fileName)
-{
-    std::ifstream file(fileName);
-    if (!file.is_open()) {
-        throw std::invalid_argument("Unable to open file");
-    }
-    std::string line;
-    std::string keyword;
-    std::string attribute;
-    std::string type;
-    std::string type_w;
-    // Read file
-    while (getline(file, line)) {
-        if (line.empty() || line[0] == '%' || line == "\r" || line == " ") {
-            // Skip comments and empty lines
-            continue;
-        }
-        if (line.find("@attribute") != std::string::npos || line.find("@ATTRIBUTE") != std::string::npos) {
-            // Read attributes
-            std::stringstream ss(line);
-            ss >> keyword >> attribute;
-            type = "";
-            while (ss >> type_w)
-                type += type_w + " ";
-            attributes.emplace_back(trim(attribute), trim(type));
-            continue;
-        }
-        if (line[0] == '@') {
-            continue;
-        }
-        // Read data
-        lines.push_back(line);
-    }
-    file.close();
-    if (attributes.empty())
-        throw std::invalid_argument("No attributes found");
-}
-
-void ArffFiles::load(const std::string& fileName, bool classLast)
-{
-    int labelIndex;
-    loadCommon(fileName);
-    if (classLast) {
-        className = std::get<0>(attributes.back());
-        classType = std::get<1>(attributes.back());
-        attributes.pop_back();
-        labelIndex = static_cast<int>(attributes.size());
-    } else {
-        className = std::get<0>(attributes.front());
-        classType = std::get<1>(attributes.front());
-        attributes.erase(attributes.begin());
-        labelIndex = 0;
-    }
-    generateDataset(labelIndex);
-}
-void ArffFiles::load(const std::string& fileName, const std::string& name)
-{
-    int labelIndex;
-    loadCommon(fileName);
-    bool found = false;
-    for (int i = 0; i < attributes.size(); ++i) {
-        if (attributes[i].first == name) {
-            className = std::get<0>(attributes[i]);
-            classType = std::get<1>(attributes[i]);
-            attributes.erase(attributes.begin() + i);
-            labelIndex = i;
-            found = true;
-            break;
-        }
-    }
-    if (!found) {
-        throw std::invalid_argument("Class name not found");
-    }
-    generateDataset(labelIndex);
-}
-
-void ArffFiles::generateDataset(int labelIndex)
-{
-    X = std::vector<std::vector<float>>(attributes.size(), std::vector<float>(lines.size()));
-    auto yy = std::vector<std::string>(lines.size(), "");
-    auto removeLines = std::vector<int>(); // Lines with missing values
-    size_t numLine = 0;
-    for (numLine = 0; numLine < lines.size(); numLine++) {
-        std::stringstream ss(lines[numLine]);
-        std::string value;
-        int pos = 0;
-        int xIndex = 0;
-        while (getline(ss, value, ',')) {
-            if (pos++ == labelIndex) {
-                yy[numLine] = value;
-            } else {
-                if (value == "?") {
-                    X[xIndex++][numLine] = -1;
-                    removeLines.push_back(numLine);
-                } else
-                    X[xIndex++][numLine] = stof(value);
-            }
-        }
-    }
-    for (auto i : removeLines) {
-        yy.erase(yy.begin() + i);
-        for (auto& x : X) {
-            x.erase(x.begin() + i);
-        }
-    }
-    y = factorize(yy);
-}
-
-std::string ArffFiles::trim(const std::string& source)
-{
-    std::string s(source);
-    s.erase(0, s.find_first_not_of(" '\n\r\t"));
-    s.erase(s.find_last_not_of(" '\n\r\t") + 1);
-    return s;
-}
-
-std::vector<int> ArffFiles::factorize(const std::vector<std::string>& labels_t)
-{
-    std::vector<int> yy;
-    yy.reserve(labels_t.size());
-    std::map<std::string, int> labelMap;
-    int i = 0;
-    for (const std::string& label : labels_t) {
-        if (labelMap.find(label) == labelMap.end()) {
-            labelMap[label] = i++;
-        }
-        yy.push_back(labelMap[label]);
-    }
-    return yy;
-}

lib/Files/ArffFiles.h

@@ -1,33 +0,0 @@
-#ifndef ARFFFILES_H
-#define ARFFFILES_H
-
-#include <string>
-#include <vector>
-
-class ArffFiles {
-public:
-    ArffFiles() = default;
-    void load(const std::string&, bool = true);
-    void load(const std::string&, const std::string&);
-    std::vector<std::string> getLines() const;
-    unsigned long int getSize() const;
-    std::string getClassName() const;
-    std::string getClassType() const;
-    static std::string trim(const std::string&);
-    std::vector<std::vector<float>>& getX();
-    std::vector<int>& getY();
-    std::vector<std::pair<std::string, std::string>> getAttributes() const;
-    static std::vector<int> factorize(const std::vector<std::string>& labels_t);
-protected:
-    std::vector<std::string> lines;
-    std::vector<std::pair<std::string, std::string>> attributes;
-    std::string className;
-    std::string classType;
-    std::vector<std::vector<float>> X;
-    std::vector<int> y;
-    int maxLines = 0;
-    void generateDataset(int);
-    void loadCommon(const std::string&);
-};
-
-#endif

lib/Files/CMakeLists.txt

@@ -1 +0,0 @@
-add_library(ArffFiles ArffFiles.cc)

sample/CMakeLists.txt

@@ -14,7 +14,6 @@ include_directories(
     /usr/local/include
 )
 
-add_subdirectory(lib/Files)
 add_subdirectory(lib/mdlp)
 add_executable(bayesnet_sample sample.cc) 
-target_link_libraries(bayesnet_sample ArffFiles mdlp "${TORCH_LIBRARIES}" "${BayesNet}")
+target_link_libraries(bayesnet_sample mdlp "${TORCH_LIBRARIES}" "${BayesNet}")

sample/lib/Files/ArffFiles.cc

@@ -1,174 +0,0 @@
-// ***************************************************************
-// SPDX-FileCopyrightText: Copyright 2024 Ricardo Montañana Gómez
-// SPDX-FileType: SOURCE
-// SPDX-License-Identifier: MIT
-// ***************************************************************
-
-#include "ArffFiles.h"
-#include <fstream>
-#include <sstream>
-#include <map>
-#include <iostream>
-
-ArffFiles::ArffFiles() = default;
-
-std::vector<std::string> ArffFiles::getLines() const
-{
-    return lines;
-}
-
-unsigned long int ArffFiles::getSize() const
-{
-    return lines.size();
-}
-
-std::vector<std::pair<std::string, std::string>> ArffFiles::getAttributes() const
-{
-    return attributes;
-}
-
-std::string ArffFiles::getClassName() const
-{
-    return className;
-}
-
-std::string ArffFiles::getClassType() const
-{
-    return classType;
-}
-
-std::vector<std::vector<float>>& ArffFiles::getX()
-{
-    return X;
-}
-
-std::vector<int>& ArffFiles::getY()
-{
-    return y;
-}
-
-void ArffFiles::loadCommon(std::string fileName)
-{
-    std::ifstream file(fileName);
-    if (!file.is_open()) {
-        throw std::invalid_argument("Unable to open file");
-    }
-    std::string line;
-    std::string keyword;
-    std::string attribute;
-    std::string type;
-    std::string type_w;
-    while (getline(file, line)) {
-        if (line.empty() || line[0] == '%' || line == "\r" || line == " ") {
-            continue;
-        }
-        if (line.find("@attribute") != std::string::npos || line.find("@ATTRIBUTE") != std::string::npos) {
-            std::stringstream ss(line);
-            ss >> keyword >> attribute;
-            type = "";
-            while (ss >> type_w)
-                type += type_w + " ";
-            attributes.emplace_back(trim(attribute), trim(type));
-            continue;
-        }
-        if (line[0] == '@') {
-            continue;
-        }
-        lines.push_back(line);
-    }
-    file.close();
-    if (attributes.empty())
-        throw std::invalid_argument("No attributes found");
-}
-
-void ArffFiles::load(const std::string& fileName, bool classLast)
-{
-    int labelIndex;
-    loadCommon(fileName);
-    if (classLast) {
-        className = std::get<0>(attributes.back());
-        classType = std::get<1>(attributes.back());
-        attributes.pop_back();
-        labelIndex = static_cast<int>(attributes.size());
-    } else {
-        className = std::get<0>(attributes.front());
-        classType = std::get<1>(attributes.front());
-        attributes.erase(attributes.begin());
-        labelIndex = 0;
-    }
-    generateDataset(labelIndex);
-}
-void ArffFiles::load(const std::string& fileName, const std::string& name)
-{
-    int labelIndex;
-    loadCommon(fileName);
-    bool found = false;
-    for (int i = 0; i < attributes.size(); ++i) {
-        if (attributes[i].first == name) {
-            className = std::get<0>(attributes[i]);
-            classType = std::get<1>(attributes[i]);
-            attributes.erase(attributes.begin() + i);
-            labelIndex = i;
-            found = true;
-            break;
-        }
-    }
-    if (!found) {
-        throw std::invalid_argument("Class name not found");
-    }
-    generateDataset(labelIndex);
-}
-
-void ArffFiles::generateDataset(int labelIndex)
-{
-    X = std::vector<std::vector<float>>(attributes.size(), std::vector<float>(lines.size()));
-    auto yy = std::vector<std::string>(lines.size(), "");
-    auto removeLines = std::vector<int>(); // Lines with missing values
-    for (size_t i = 0; i < lines.size(); i++) {
-        std::stringstream ss(lines[i]);
-        std::string value;
-        int pos = 0;
-        int xIndex = 0;
-        while (getline(ss, value, ',')) {
-            if (pos++ == labelIndex) {
-                yy[i] = value;
-            } else {
-                if (value == "?") {
-                    X[xIndex++][i] = -1;
-                    removeLines.push_back(i);
-                } else
-                    X[xIndex++][i] = stof(value);
-            }
-        }
-    }
-    for (auto i : removeLines) {
-        yy.erase(yy.begin() + i);
-        for (auto& x : X) {
-            x.erase(x.begin() + i);
-        }
-    }
-    y = factorize(yy);
-}
-
-std::string ArffFiles::trim(const std::string& source)
-{
-    std::string s(source);
-    s.erase(0, s.find_first_not_of(" '\n\r\t"));
-    s.erase(s.find_last_not_of(" '\n\r\t") + 1);
-    return s;
-}
-
-std::vector<int> ArffFiles::factorize(const std::vector<std::string>& labels_t)
-{
-    std::vector<int> yy;
-    yy.reserve(labels_t.size());
-    std::map<std::string, int> labelMap;
-    int i = 0;
-    for (const std::string& label : labels_t) {
-        if (labelMap.find(label) == labelMap.end()) {
-            labelMap[label] = i++;
-        }
-        yy.push_back(labelMap[label]);
-    }
-    return yy;
-}

sample/lib/Files/ArffFiles.h

@@ -1,38 +0,0 @@
-// ***************************************************************
-// SPDX-FileCopyrightText: Copyright 2024 Ricardo Montañana Gómez
-// SPDX-FileType: SOURCE
-// SPDX-License-Identifier: MIT
-// ***************************************************************
-
-#ifndef ARFFFILES_H
-#define ARFFFILES_H
-
-#include <string>
-#include <vector>
-
-class ArffFiles {
-private:
-    std::vector<std::string> lines;
-    std::vector<std::pair<std::string, std::string>> attributes;
-    std::string className;
-    std::string classType;
-    std::vector<std::vector<float>> X;
-    std::vector<int> y;
-    void generateDataset(int);
-    void loadCommon(std::string);
-public:
-    ArffFiles();
-    void load(const std::string&, bool = true);
-    void load(const std::string&, const std::string&);
-    std::vector<std::string> getLines() const;
-    unsigned long int getSize() const;
-    std::string getClassName() const;
-    std::string getClassType() const;
-    static std::string trim(const std::string&);
-    std::vector<std::vector<float>>& getX();
-    std::vector<int>& getY();
-    std::vector<std::pair<std::string, std::string>> getAttributes() const;
-    static std::vector<int> factorize(const std::vector<std::string>& labels_t);
-};
-
-#endif

sample/lib/Files/CMakeLists.txt

@@ -1 +0,0 @@
-add_library(ArffFiles ArffFiles.cc)

sample/sample.cc

@@ -4,7 +4,7 @@
 // SPDX-License-Identifier: MIT
 // ***************************************************************
 
-#include <ArffFiles.h>
+#include <ArffFiles.hpp>
 #include <CPPFImdlp.h>
 #include <bayesnet/ensembles/BoostAODE.h>
 

tests/CMakeLists.txt

@@ -1,6 +1,6 @@
 if(ENABLE_TESTING)
     include_directories(
-        ${BayesNet_SOURCE_DIR}/lib/Files
+        ${BayesNet_SOURCE_DIR}/tests/lib/Files
         ${BayesNet_SOURCE_DIR}/lib/folding
         ${BayesNet_SOURCE_DIR}/lib/mdlp
         ${BayesNet_SOURCE_DIR}/lib/json/include
@@ -10,17 +10,18 @@ if(ENABLE_TESTING)
     file(GLOB_RECURSE BayesNet_SOURCES "${BayesNet_SOURCE_DIR}/bayesnet/*.cc")
     add_executable(TestBayesNet TestBayesNetwork.cc TestBayesNode.cc TestBayesClassifier.cc 
         TestBayesModels.cc TestBayesMetrics.cc TestFeatureSelection.cc TestBoostAODE.cc TestA2DE.cc 
-        TestUtils.cc TestBayesEnsemble.cc TestModulesVersions.cc ${BayesNet_SOURCES})
-    target_link_libraries(TestBayesNet PUBLIC "${TORCH_LIBRARIES}" ArffFiles mdlp PRIVATE Catch2::Catch2WithMain)
+        TestUtils.cc TestBayesEnsemble.cc TestModulesVersions.cc TestBoostA2DE.cc ${BayesNet_SOURCES})
+    target_link_libraries(TestBayesNet PUBLIC "${TORCH_LIBRARIES}" mdlp PRIVATE Catch2::Catch2WithMain)
     add_test(NAME BayesNetworkTest COMMAND TestBayesNet)
-    add_test(NAME Network COMMAND TestBayesNet "[Network]")
-    add_test(NAME Node COMMAND TestBayesNet "[Node]")
-    add_test(NAME Metrics COMMAND TestBayesNet "[Metrics]")
-    add_test(NAME FeatureSelection COMMAND TestBayesNet "[FeatureSelection]")
+    add_test(NAME A2DE COMMAND TestBayesNet "[A2DE]")
+    add_test(NAME BoostA2DE COMMAND TestBayesNet "[BoostA2DE]")
+    add_test(NAME BoostAODE COMMAND TestBayesNet "[BoostAODE]")
     add_test(NAME Classifier COMMAND TestBayesNet "[Classifier]")
     add_test(NAME Ensemble COMMAND TestBayesNet "[Ensemble]")
+    add_test(NAME FeatureSelection COMMAND TestBayesNet "[FeatureSelection]")
+    add_test(NAME Metrics COMMAND TestBayesNet "[Metrics]")
     add_test(NAME Models COMMAND TestBayesNet "[Models]")
-    add_test(NAME BoostAODE COMMAND TestBayesNet "[BoostAODE]")
-    add_test(NAME A2DE COMMAND TestBayesNet "[A2DE]")
     add_test(NAME Modules COMMAND TestBayesNet "[Modules]")
+    add_test(NAME Network COMMAND TestBayesNet "[Network]")
+    add_test(NAME Node COMMAND TestBayesNet "[Node]")
 endif(ENABLE_TESTING)

tests/TestBayesMetrics.cc

@@ -11,7 +11,6 @@
 #include "TestUtils.h"
 #include "Timer.h"
 
-
 TEST_CASE("Metrics Test", "[Metrics]")
 {
     std::string file_name = GENERATE("glass", "iris", "ecoli", "diabetes");
@@ -28,8 +27,8 @@ TEST_CASE("Metrics Test", "[Metrics]")
         {"diabetes", 0.0345470614}
     };
     map<pair<std::string, int>, std::vector<pair<int, int>>> resultsMST = {
-        { {"glass", 0}, { {0, 6}, {0, 5}, {0, 3}, {5, 1}, {5, 8}, {5, 4}, {6, 2}, {6, 7} } },
-        { {"glass", 1}, { {1, 5}, {5, 0}, {5, 8}, {5, 4}, {0, 6}, {0, 3}, {6, 2}, {6, 7} } },
+        { {"glass", 0}, { {0, 6}, {0, 5}, {0, 3}, {3, 4}, {5, 1}, {5, 8}, {6, 2}, {6, 7} } },
+        { {"glass", 1}, { {1, 5}, {5, 0}, {5, 8}, {0, 6}, {0, 3}, {3, 4}, {6, 2}, {6, 7} } },
         { {"iris", 0}, { {0, 1}, {0, 2}, {1, 3} } },
         { {"iris", 1}, { {1, 0}, {1, 3}, {0, 2} } },
         { {"ecoli", 0}, { {0, 1}, {0, 2}, {1, 5}, {1, 3}, {5, 6}, {5, 4} } },
@@ -100,6 +99,25 @@ TEST_CASE("Entropy Test", "[Metrics]")
     REQUIRE(result == Catch::Approx(0.693147180559945).epsilon(raw.epsilon));
 }
 TEST_CASE("Conditional Entropy", "[Metrics]")
+{
+    auto raw = RawDatasets("iris", true);
+    bayesnet::Metrics metrics(raw.dataset, raw.features, raw.className, raw.classNumStates);
+    auto expected = std::map<std::pair<int, int>, double>{
+        { { 0, 1 }, 1.32674 },
+        { { 0, 2 }, 0.236253 },
+        { { 0, 3 }, 0.1202 },
+        { { 1, 2 }, 0.252551 },
+        { { 1, 3 }, 0.10515 },
+        { { 2, 3 }, 0.108323 },
+    };
+    for (int i = 0; i < raw.features.size() - 1; ++i) {
+        for (int j = i + 1; j < raw.features.size(); ++j) {
+            double result = metrics.conditionalEntropy(raw.dataset.index({ i, "..." }), raw.dataset.index({ j, "..." }), raw.yt, raw.weights);
+            REQUIRE(result == Catch::Approx(expected.at({ i, j })).epsilon(raw.epsilon));
+        }
+    }
+}
+TEST_CASE("Conditional Mutual Information", "[Metrics]")
 {
     auto raw = RawDatasets("iris", true);
     bayesnet::Metrics metrics(raw.dataset, raw.features, raw.className, raw.classNumStates);
@@ -118,3 +136,132 @@ TEST_CASE("Conditional Entropy", "[Metrics]")
         }
     }
 }
+TEST_CASE("Select K Pairs descending", "[Metrics]")
+{
+    auto raw = RawDatasets("iris", true);
+    bayesnet::Metrics metrics(raw.dataset, raw.features, raw.className, raw.classNumStates);
+    std::vector<int> empty;
+    auto results = metrics.SelectKPairs(raw.weights, empty, false);
+    auto expected = std::vector<std::pair<std::pair<int, int>, double>>{
+        { { 1, 3 }, 1.31852 },
+        { { 1, 2 }, 1.17112 },
+        { { 0, 3 }, 0.403749 },
+        { { 0, 2 }, 0.287696 },
+        { { 2, 3 }, 0.210068 },
+        { { 0, 1 }, 0.0 },
+    };
+    auto scores = metrics.getScoresKPairs();
+    for (int i = 0; i < results.size(); ++i) {
+        auto result = results[i];
+        auto expect = expected[i];
+        auto score = scores[i];
+        REQUIRE(result.first == expect.first.first);
+        REQUIRE(result.second == expect.first.second);
+        REQUIRE(score.first.first == expect.first.first);
+        REQUIRE(score.first.second == expect.first.second);
+        REQUIRE(score.second == Catch::Approx(expect.second).epsilon(raw.epsilon));
+    }
+    REQUIRE(results.size() == 6);
+    REQUIRE(scores.size() == 6);
+}
+TEST_CASE("Select K Pairs ascending", "[Metrics]")
+{
+    auto raw = RawDatasets("iris", true);
+    bayesnet::Metrics metrics(raw.dataset, raw.features, raw.className, raw.classNumStates);
+    std::vector<int> empty;
+    auto results = metrics.SelectKPairs(raw.weights, empty, true);
+    auto expected = std::vector<std::pair<std::pair<int, int>, double>>{
+        { { 0, 1 }, 0.0 },
+        { { 2, 3 }, 0.210068 },
+        { { 0, 2 }, 0.287696 },
+        { { 0, 3 }, 0.403749 },
+        { { 1, 2 }, 1.17112 },
+        { { 1, 3 }, 1.31852 },
+    };
+    auto scores = metrics.getScoresKPairs();
+    for (int i = 0; i < results.size(); ++i) {
+        auto result = results[i];
+        auto expect = expected[i];
+        auto score = scores[i];
+        REQUIRE(result.first == expect.first.first);
+        REQUIRE(result.second == expect.first.second);
+        REQUIRE(score.first.first == expect.first.first);
+        REQUIRE(score.first.second == expect.first.second);
+        REQUIRE(score.second == Catch::Approx(expect.second).epsilon(raw.epsilon));
+    }
+    REQUIRE(results.size() == 6);
+    REQUIRE(scores.size() == 6);
+}
+TEST_CASE("Select K Pairs with features excluded", "[Metrics]")
+{
+    auto raw = RawDatasets("iris", true);
+    bayesnet::Metrics metrics(raw.dataset, raw.features, raw.className, raw.classNumStates);
+    std::vector<int> excluded = { 0, 3 };
+    auto results = metrics.SelectKPairs(raw.weights, excluded, true);
+    auto expected = std::vector<std::pair<std::pair<int, int>, double>>{
+        { { 1, 2 }, 1.17112 },
+    };
+    auto scores = metrics.getScoresKPairs();
+    for (int i = 0; i < results.size(); ++i) {
+        auto result = results[i];
+        auto expect = expected[i];
+        auto score = scores[i];
+        REQUIRE(result.first == expect.first.first);
+        REQUIRE(result.second == expect.first.second);
+        REQUIRE(score.first.first == expect.first.first);
+        REQUIRE(score.first.second == expect.first.second);
+        REQUIRE(score.second == Catch::Approx(expect.second).epsilon(raw.epsilon));
+    }
+    REQUIRE(results.size() == 1);
+    REQUIRE(scores.size() == 1);
+}
+TEST_CASE("Select K Pairs with number of pairs descending", "[Metrics]")
+{
+    auto raw = RawDatasets("iris", true);
+    bayesnet::Metrics metrics(raw.dataset, raw.features, raw.className, raw.classNumStates);
+    std::vector<int> empty;
+    auto results = metrics.SelectKPairs(raw.weights, empty, false, 3);
+    auto expected = std::vector<std::pair<std::pair<int, int>, double>>{
+        { { 1, 3 }, 1.31852 },
+        { { 1, 2 }, 1.17112 },
+        { { 0, 3 }, 0.403749 }
+    };
+    auto scores = metrics.getScoresKPairs();
+    REQUIRE(results.size() == 3);
+    REQUIRE(scores.size() == 3);
+    for (int i = 0; i < results.size(); ++i) {
+        auto result = results[i];
+        auto expect = expected[i];
+        auto score = scores[i];
+        REQUIRE(result.first == expect.first.first);
+        REQUIRE(result.second == expect.first.second);
+        REQUIRE(score.first.first == expect.first.first);
+        REQUIRE(score.first.second == expect.first.second);
+        REQUIRE(score.second == Catch::Approx(expect.second).epsilon(raw.epsilon));
+    }
+}
+TEST_CASE("Select K Pairs with number of pairs ascending", "[Metrics]")
+{
+    auto raw = RawDatasets("iris", true);
+    bayesnet::Metrics metrics(raw.dataset, raw.features, raw.className, raw.classNumStates);
+    std::vector<int> empty;
+    auto results = metrics.SelectKPairs(raw.weights, empty, true, 3);
+    auto expected = std::vector<std::pair<std::pair<int, int>, double>>{
+        { { 0, 3 }, 0.403749 },
+        { { 1, 2 }, 1.17112 },
+        { { 1, 3 }, 1.31852 }
+    };
+    auto scores = metrics.getScoresKPairs();
+    REQUIRE(results.size() == 3);
+    REQUIRE(scores.size() == 3);
+    for (int i = 0; i < results.size(); ++i) {
+        auto result = results[i];
+        auto expect = expected[i];
+        auto score = scores[i];
+        REQUIRE(result.first == expect.first.first);
+        REQUIRE(result.second == expect.first.second);
+        REQUIRE(score.first.first == expect.first.first);
+        REQUIRE(score.first.second == expect.first.second);
+        REQUIRE(score.second == Catch::Approx(expect.second).epsilon(raw.epsilon));
+    }
+}

tests/TestBayesModels.cc

@@ -56,14 +56,14 @@ TEST_CASE("Test Bayesian Classifiers score & version", "[Models]")
             auto raw = RawDatasets(file_name, discretize);
             clf->fit(raw.Xt, raw.yt, raw.features, raw.className, raw.states);
             auto score = clf->score(raw.Xt, raw.yt);
-            INFO("Classifier: " + name + " File: " + file_name);
+            INFO("Classifier: " << name << " File: " << file_name);
             REQUIRE(score == Catch::Approx(scores[{file_name, name}]).epsilon(raw.epsilon));
             REQUIRE(clf->getStatus() == bayesnet::NORMAL);
         }
     }
     SECTION("Library check version")
     {
-        INFO("Checking version of " + name + " classifier");
+        INFO("Checking version of " << name << " classifier");
         REQUIRE(clf->getVersion() == ACTUAL_VERSION);
     }
     delete clf;

tests/TestBoostA2DE.cc

@@ -0,0 +1,215 @@
+// ***************************************************************
+// SPDX-FileCopyrightText: Copyright 2024 Ricardo Montañana Gómez
+// SPDX-FileType: SOURCE
+// SPDX-License-Identifier: MIT
+// ***************************************************************
+
+#include <type_traits>
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/catch_approx.hpp>
+#include <catch2/generators/catch_generators.hpp>
+#include "bayesnet/utils/BayesMetrics.h"
+#include "bayesnet/ensembles/BoostA2DE.h"
+#include "TestUtils.h"
+
+
+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);
+    REQUIRE(clf.getNumberOfNodes() == 342);
+    REQUIRE(clf.getNumberOfEdges() == 684);
+    REQUIRE(clf.getNotes().size() == 3);
+    REQUIRE(clf.getNotes()[0] == "Convergence threshold reached & 15 models eliminated");
+    REQUIRE(clf.getNotes()[1] == "Pairs not used in train: 20");
+    REQUIRE(clf.getNotes()[2] == "Number of models: 38");
+    auto score = clf.score(raw.Xv, raw.yv);
+    REQUIRE(score == Catch::Approx(0.919271).epsilon(raw.epsilon));
+}
+// 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);
+//     REQUIRE(clf.getNumberOfNodes() == 90);
+//     REQUIRE(clf.getNumberOfEdges() == 153);
+//     REQUIRE(clf.getNotes().size() == 2);
+//     REQUIRE(clf.getNotes()[0] == "Used features in initialization: 4 of 9 with IWSS");
+//     REQUIRE(clf.getNotes()[1] == "Number of models: 9");
+// }
+// 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);
+//     REQUIRE(clf.getNumberOfNodes() == 90);
+//     REQUIRE(clf.getNumberOfEdges() == 153);
+//     REQUIRE(clf.getNotes().size() == 2);
+//     REQUIRE(clf.getNotes()[0] == "Used features in initialization: 4 of 9 with FCBF");
+//     REQUIRE(clf.getNotes()[1] == "Number of models: 9");
+// }
+// TEST_CASE("Test used features in train note and score", "[BoostAODE]")
+// {
+//     auto raw = RawDatasets("diabetes", true);
+//     auto clf = bayesnet::BoostAODE(true);
+//     clf.setHyperparameters({
+//         {"order", "asc"},
+//         {"convergence", true},
+//         {"select_features","CFS"},
+//         });
+//     clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states);
+//     REQUIRE(clf.getNumberOfNodes() == 72);
+//     REQUIRE(clf.getNumberOfEdges() == 120);
+//     REQUIRE(clf.getNotes().size() == 2);
+//     REQUIRE(clf.getNotes()[0] == "Used features in initialization: 6 of 8 with CFS");
+//     REQUIRE(clf.getNotes()[1] == "Number of models: 8");
+//     auto score = clf.score(raw.Xv, raw.yv);
+//     auto scoret = clf.score(raw.Xt, raw.yt);
+//     REQUIRE(score == Catch::Approx(0.809895813).epsilon(raw.epsilon));
+//     REQUIRE(scoret == Catch::Approx(0.809895813).epsilon(raw.epsilon));
+// }
+// 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);
+//     auto score_proba = clf.score(raw.Xv, raw.yv);
+//     auto pred_proba = clf.predict_proba(raw.Xv);
+//     clf.setHyperparameters({
+//         {"predict_voting",true},
+//         });
+//     auto score_voting = clf.score(raw.Xv, raw.yv);
+//     auto pred_voting = clf.predict_proba(raw.Xv);
+//     REQUIRE(score_proba == Catch::Approx(0.97333).epsilon(raw.epsilon));
+//     REQUIRE(score_voting == Catch::Approx(0.98).epsilon(raw.epsilon));
+//     REQUIRE(pred_voting[83][2] == Catch::Approx(1.0).epsilon(raw.epsilon));
+//     REQUIRE(pred_proba[83][2] == Catch::Approx(0.86121525).epsilon(raw.epsilon));
+//     REQUIRE(clf.dump_cpt() == "");
+//     REQUIRE(clf.topological_order() == std::vector<std::string>());
+// }
+// TEST_CASE("Order asc, desc & random", "[BoostAODE]")
+// {
+//     auto raw = RawDatasets("glass", true);
+//     std::map<std::string, double> scores{
+//         {"asc", 0.83645f }, { "desc", 0.84579f }, { "rand", 0.84112 }
+//     };
+//     for (const std::string& order : { "asc", "desc", "rand" }) {
+//         auto clf = bayesnet::BoostAODE();
+//         clf.setHyperparameters({
+//             {"order", order},
+//             {"bisection", false},
+//             {"maxTolerance", 1},
+//             {"convergence", false},
+//             });
+//         clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states);
+//         auto score = clf.score(raw.Xv, raw.yv);
+//         auto scoret = clf.score(raw.Xt, raw.yt);
+//         INFO("BoostAODE order: " + order);
+//         REQUIRE(score == Catch::Approx(scores[order]).epsilon(raw.epsilon));
+//         REQUIRE(scoret == Catch::Approx(scores[order]).epsilon(raw.epsilon));
+//     }
+// }
+// TEST_CASE("Oddities", "[BoostAODE]")
+// {
+//     auto clf = bayesnet::BoostAODE();
+//     auto raw = RawDatasets("iris", true);
+//     auto bad_hyper = nlohmann::json{
+//         { { "order", "duck" } },
+//         { { "select_features", "duck" } },
+//         { { "maxTolerance", 0 } },
+//         { { "maxTolerance", 5 } },
+//     };
+//     for (const auto& hyper : bad_hyper.items()) {
+//         INFO("BoostAODE hyper: " + hyper.value().dump());
+//         REQUIRE_THROWS_AS(clf.setHyperparameters(hyper.value()), std::invalid_argument);
+//     }
+//     REQUIRE_THROWS_AS(clf.setHyperparameters({ {"maxTolerance", 0 } }), std::invalid_argument);
+//     auto bad_hyper_fit = nlohmann::json{
+//         { { "select_features","IWSS" }, { "threshold", -0.01 } },
+//         { { "select_features","IWSS" }, { "threshold", 0.51 } },
+//         { { "select_features","FCBF" }, { "threshold", 1e-8 } },
+//         { { "select_features","FCBF" }, { "threshold", 1.01 } },
+//     };
+//     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);
+//     }
+// }
+
+// TEST_CASE("Bisection Best", "[BoostAODE]")
+// {
+//     auto clf = bayesnet::BoostAODE();
+//     auto raw = RawDatasets("kdd_JapaneseVowels", true, 1200, true, false);
+//     clf.setHyperparameters({
+//         {"bisection", true},
+//         {"maxTolerance", 3},
+//         {"convergence", true},
+//         {"block_update", false},
+//         {"convergence_best", false},
+//         });
+//     clf.fit(raw.X_train, raw.y_train, raw.features, raw.className, raw.states);
+//     REQUIRE(clf.getNumberOfNodes() == 210);
+//     REQUIRE(clf.getNumberOfEdges() == 378);
+//     REQUIRE(clf.getNotes().size() == 1);
+//     REQUIRE(clf.getNotes().at(0) == "Number of models: 14");
+//     auto score = clf.score(raw.X_test, raw.y_test);
+//     auto scoret = clf.score(raw.X_test, raw.y_test);
+//     REQUIRE(score == Catch::Approx(0.991666675f).epsilon(raw.epsilon));
+//     REQUIRE(scoret == Catch::Approx(0.991666675f).epsilon(raw.epsilon));
+// }
+// TEST_CASE("Bisection Best vs Last", "[BoostAODE]")
+// {
+//     auto raw = RawDatasets("kdd_JapaneseVowels", true, 1500, true, false);
+//     auto clf = bayesnet::BoostAODE(true);
+//     auto hyperparameters = nlohmann::json{
+//         {"bisection", true},
+//         {"maxTolerance", 3},
+//         {"convergence", true},
+//         {"convergence_best", true},
+//     };
+//     clf.setHyperparameters(hyperparameters);
+//     clf.fit(raw.X_train, raw.y_train, raw.features, raw.className, raw.states);
+//     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);
+//     auto score_last = clf.score(raw.X_test, raw.y_test);
+//     REQUIRE(score_last == Catch::Approx(0.976666689f).epsilon(raw.epsilon));
+// }
+
+// TEST_CASE("Block Update", "[BoostAODE]")
+// {
+//     auto clf = bayesnet::BoostAODE();
+//     auto raw = RawDatasets("mfeat-factors", true, 500);
+//     clf.setHyperparameters({
+//         {"bisection", true},
+//         {"block_update", true},
+//         {"maxTolerance", 3},
+//         {"convergence", true},
+//         });
+//     clf.fit(raw.X_train, raw.y_train, raw.features, raw.className, raw.states);
+//     REQUIRE(clf.getNumberOfNodes() == 868);
+//     REQUIRE(clf.getNumberOfEdges() == 1724);
+//     REQUIRE(clf.getNotes().size() == 3);
+//     REQUIRE(clf.getNotes()[0] == "Convergence threshold reached & 15 models eliminated");
+//     REQUIRE(clf.getNotes()[1] == "Used features in train: 19 of 216");
+//     REQUIRE(clf.getNotes()[2] == "Number of models: 4");
+//     auto score = clf.score(raw.X_test, raw.y_test);
+//     auto scoret = clf.score(raw.X_test, raw.y_test);
+//     REQUIRE(score == Catch::Approx(0.99f).epsilon(raw.epsilon));
+//     REQUIRE(scoret == Catch::Approx(0.99f).epsilon(raw.epsilon));
+//     //
+//     // std::cout << "Number of nodes " << clf.getNumberOfNodes() << std::endl;
+//     // std::cout << "Number of edges " << clf.getNumberOfEdges() << std::endl;
+//     // std::cout << "Notes size " << clf.getNotes().size() << std::endl;
+//     // for (auto note : clf.getNotes()) {
+//     //     std::cout << note << std::endl;
+//     // }
+//     // std::cout << "Score " << score << std::endl;
+// }

tests/TestBoostAODE.cc

@@ -8,6 +8,7 @@
 #include <catch2/catch_test_macros.hpp>
 #include <catch2/catch_approx.hpp>
 #include <catch2/generators/catch_generators.hpp>
+#include <catch2/matchers/catch_matchers.hpp>
 #include "bayesnet/ensembles/BoostAODE.h"
 #include "TestUtils.h"
 
@@ -104,7 +105,7 @@ TEST_CASE("Order asc, desc & random", "[BoostAODE]")
         clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states);
         auto score = clf.score(raw.Xv, raw.yv);
         auto scoret = clf.score(raw.Xt, raw.yt);
-        INFO("BoostAODE order: " + order);
+        INFO("BoostAODE order: " << order);
         REQUIRE(score == Catch::Approx(scores[order]).epsilon(raw.epsilon));
         REQUIRE(scoret == Catch::Approx(scores[order]).epsilon(raw.epsilon));
     }
@@ -120,7 +121,7 @@ TEST_CASE("Oddities", "[BoostAODE]")
         { { "maxTolerance", 5 } },
     };
     for (const auto& hyper : bad_hyper.items()) {
-        INFO("BoostAODE hyper: " + hyper.value().dump());
+        INFO("BoostAODE hyper: " << hyper.value().dump());
         REQUIRE_THROWS_AS(clf.setHyperparameters(hyper.value()), std::invalid_argument);
     }
     REQUIRE_THROWS_AS(clf.setHyperparameters({ {"maxTolerance", 0 } }), std::invalid_argument);
@@ -131,7 +132,7 @@ TEST_CASE("Oddities", "[BoostAODE]")
         { { "select_features","FCBF" }, { "threshold", 1.01 } },
     };
     for (const auto& hyper : bad_hyper_fit.items()) {
-        INFO("BoostAODE hyper: " + hyper.value().dump());
+        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);
     }

tests/TestModulesVersions.cc

@@ -18,7 +18,8 @@
 std::map<std::string, std::string> modules = {
     { "mdlp", "1.1.2" },
     { "Folding", "1.1.0" },
-    { "json", "3.11" }
+    { "json", "3.11" },
+    { "ArffFiles", "1.0.0" }
 };
 
 TEST_CASE("MDLP", "[Modules]")
@@ -35,3 +36,8 @@ TEST_CASE("NLOHMANN_JSON", "[Modules]")
 {
     REQUIRE(JSON_VERSION == modules["json"]);
 }
+TEST_CASE("ArffFiles", "[Modules]")
+{
+    auto handler = ArffFiles();
+    REQUIRE(handler.version() == modules["ArffFiles"]);
+}

tests/TestUtils.h

@@ -11,7 +11,7 @@
 #include <vector>
 #include <map>
 #include <tuple>
-#include <ArffFiles.h>
+#include <ArffFiles.hpp>
 #include <CPPFImdlp.h>
 #include <folding.hpp>