Add selectKParis method

Reviewed-on: #28

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,172 @@
+// ***************************************************************
+// 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 };
+        // while (!finished) {
+        //     // Step 1: Build ranking with mutual information
+        //     auto pairSelection = metrics.SelectKBestWeighted(weights_, ascending, n); // Get all the features sorted
+        //     if (order_algorithm == Orders.RAND) {
+        //         std::shuffle(featureSelection.begin(), featureSelection.end(), g);
+        //     }
+        //     // Remove used features
+        //     featureSelection.erase(remove_if(begin(featureSelection), end(featureSelection), [&](auto x)
+        //         { return std::find(begin(featuresUsed), end(featuresUsed), x) != end(featuresUsed);}),
+        //         end(featureSelection)
+        //     );
+        //     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 && featureSelection.size() > 0) {
+        //         auto feature = featureSelection[0];
+        //         featureSelection.erase(featureSelection.begin());
+        //         std::unique_ptr<Classifier> model;
+        //         model = std::make_unique<SPODE>(feature);
+        //         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++;
+        //         featuresUsed.push_back(feature);
+        //         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 || 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);
+        //         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 (featuresUsed.size() != features.size()) {
+        //     notes.push_back("Used features in train: " + std::to_string(featuresUsed.size()) + " of " + std::to_string(features.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,44 @@ namespace bayesnet {
         }
         samples.index_put_({ -1, "..." }, torch::tensor(labels, torch::kInt32));
     }
+    std::vector<std::pair<int, int>> Metrics::SelectKPairs(const torch::Tensor& weights, bool ascending, unsigned k)
+    {
+        // Return the K Best features 
+        auto n = features.size();
+        if (k == 0) {
+            k = n;
+        }
+        // compute scores
+        scoresKPairs.clear();
+        pairsKBest.clear();
+        auto label = samples.index({ -1, "..." });
+        // for (int i = 0; i < n; ++i) {
+        //     for (int j = i + 1; j < n; ++j) {
+        //         scoresKBest.push_back(mutualInformation(samples.index({ i, "..." }), samples.index({ j, "..." }), weights));
+        //         featuresKBest.push_back(i);
+        //         featuresKBest.push_back(j);
+        //     }
+        // }
+        // // sort & reduce scores and features
+        // if (ascending) {
+        //     sort(featuresKBest.begin(), featuresKBest.end(), [&](int i, int j)
+        //         { return scoresKBest[i] < scoresKBest[j]; });
+        //     sort(scoresKBest.begin(), scoresKBest.end(), std::less<double>());
+        //     if (k < n) {
+        //         for (int i = 0; i < n - k; ++i) {
+        //             featuresKBest.erase(featuresKBest.begin());
+        //             scoresKBest.erase(scoresKBest.begin());
+        //         }
+        //     }
+        // } else {
+        //     sort(featuresKBest.begin(), featuresKBest.end(), [&](int i, int j)
+        //         { return scoresKBest[i] > scoresKBest[j]; });
+        //     sort(scoresKBest.begin(), scoresKBest.end(), std::greater<double>());
+        //     featuresKBest.resize(k);
+        //     scoresKBest.resize(k);
+        // }
+        return pairsKBest;
+    }
     std::vector<int> Metrics::SelectKBestWeighted(const torch::Tensor& weights, bool ascending, unsigned k)
     {
         // Return the K Best features 
@@ -187,7 +225,7 @@ namespace bayesnet {
             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) {

bayesnet/utils/BayesMetrics.h

@@ -16,6 +16,7 @@ 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, bool ascending = false, unsigned k = 0);
         std::vector<double> getScoresKBest() 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);
@@ -41,7 +42,7 @@ namespace bayesnet {
             }
             return result;
         }
-        template <class T>
+            template <class T>
         T pop_first(std::vector<T>& v)
         {
             T temp = v[0];
@@ -52,6 +53,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::map<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)

tests/CMakeLists.txt

@@ -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})
+        TestUtils.cc TestBayesEnsemble.cc TestModulesVersions.cc TestBoostA2DE.cc ${BayesNet_SOURCES})
     target_link_libraries(TestBayesNet PUBLIC "${TORCH_LIBRARIES}" ArffFiles 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

@@ -100,6 +100,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);

tests/TestBoostA2DE.cc

@@ -0,0 +1,226 @@
+// ***************************************************************
+// 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);
+    bayesnet::Metrics metrics(raw.dataset, raw.features, raw.className, raw.classNumStates);
+    auto expected = std::map<std::pair<int, int>, double>{
+        { { 0, 1 }, 0.0 },
+        { { 0, 2 }, 0.287696 },
+        { { 0, 3 }, 0.403749 },
+        { { 1, 2 }, 1.17112 },
+        { { 1, 3 }, 1.31852 },
+        { { 2, 3 }, 0.210068 },
+    };
+    for (int i = 0; i < raw.features.size() - 1; ++i) {
+        for (int j = i + 1; j < raw.features.size(); ++j) {
+            double result = metrics.conditionalMutualInformation(raw.dataset.index({ i, "..." }), raw.dataset.index({ j, "..." }), raw.yt, raw.weights);
+            // REQUIRE(result == Catch::Approx(expected.at({ i, j })).epsilon(raw.epsilon));
+            auto clf = bayesnet::SPnDE({ i, j });
+            clf.fit(raw.Xv, raw.yv, raw.features, raw.className, raw.states);
+            auto score = clf.score(raw.Xt, raw.yt);
+            std::cout << " i " << i << " j " << j << " cmi "
+                << std::setw(8) << std::setprecision(6) << fixed << result
+                << " score = " << score << std::endl;
+        }
+    }
+}
+// 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

@@ -104,7 +104,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 +120,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 +131,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);
     }