Refactor library structure

bayesnet/ensembles/AODE.cc

@@ -0,0 +1,34 @@
+#include "AODE.h"
+
+namespace bayesnet {
+    AODE::AODE(bool predict_voting) : Ensemble(predict_voting)
+    {
+        validHyperparameters = { "predict_voting" };
+
+    }
+    void AODE::setHyperparameters(const nlohmann::json& hyperparameters_)
+    {
+        auto hyperparameters = hyperparameters_;
+        if (hyperparameters.contains("predict_voting")) {
+            predict_voting = hyperparameters["predict_voting"];
+            hyperparameters.erase("predict_voting");
+        }
+        if (!hyperparameters.empty()) {
+            throw std::invalid_argument("Invalid hyperparameters" + hyperparameters.dump());
+        }
+    }
+    void AODE::buildModel(const torch::Tensor& weights)
+    {
+        models.clear();
+        significanceModels.clear();
+        for (int i = 0; i < features.size(); ++i) {
+            models.push_back(std::make_unique<SPODE>(i));
+        }
+        n_models = models.size();
+        significanceModels = std::vector<double>(n_models, 1.0);
+    }
+    std::vector<std::string> AODE::graph(const std::string& title) const
+    {
+        return Ensemble::graph(title);
+    }
+}

bayesnet/ensembles/AODE.h

@@ -0,0 +1,16 @@
+#ifndef AODE_H
+#define AODE_H
+#include "bayesnet/classifiers/SPODE.h"
+#include "Ensemble.h"
+namespace bayesnet {
+    class AODE : public Ensemble {
+    public:
+        AODE(bool predict_voting = true);
+        virtual ~AODE() {};
+        void setHyperparameters(const nlohmann::json& hyperparameters) override;
+        std::vector<std::string> graph(const std::string& title = "AODE") const override;
+    protected:
+        void buildModel(const torch::Tensor& weights) override;
+    };
+}
+#endif

bayesnet/ensembles/AODELd.cc

@@ -0,0 +1,54 @@
+#include "AODELd.h"
+
+namespace bayesnet {
+    AODELd::AODELd(bool predict_voting) : Ensemble(predict_voting), Proposal(dataset, features, className)
+    {
+        validHyperparameters = { "predict_voting" };
+
+    }
+    void AODELd::setHyperparameters(const nlohmann::json& hyperparameters_)
+    {
+        auto hyperparameters = hyperparameters_;
+        if (hyperparameters.contains("predict_voting")) {
+            predict_voting = hyperparameters["predict_voting"];
+            hyperparameters.erase("predict_voting");
+        }
+        if (!hyperparameters.empty()) {
+            throw std::invalid_argument("Invalid hyperparameters" + hyperparameters.dump());
+        }
+    }
+    AODELd& AODELd::fit(torch::Tensor& X_, torch::Tensor& y_, const std::vector<std::string>& features_, const std::string& className_, map<std::string, std::vector<int>>& states_)
+    {
+        checkInput(X_, y_);
+        features = features_;
+        className = className_;
+        Xf = X_;
+        y = y_;
+        // Fills std::vectors Xv & yv with the data from tensors X_ (discretized) & y
+        states = fit_local_discretization(y);
+        // We have discretized the input data
+        // 1st we need to fit the model to build the normal TAN structure, TAN::fit initializes the base Bayesian network
+        Ensemble::fit(dataset, features, className, states);
+        return *this;
+
+    }
+    void AODELd::buildModel(const torch::Tensor& weights)
+    {
+        models.clear();
+        for (int i = 0; i < features.size(); ++i) {
+            models.push_back(std::make_unique<SPODELd>(i));
+        }
+        n_models = models.size();
+        significanceModels = std::vector<double>(n_models, 1.0);
+    }
+    void AODELd::trainModel(const torch::Tensor& weights)
+    {
+        for (const auto& model : models) {
+            model->fit(Xf, y, features, className, states);
+        }
+    }
+    std::vector<std::string> AODELd::graph(const std::string& name) const
+    {
+        return Ensemble::graph(name);
+    }
+}

bayesnet/ensembles/AODELd.h

@@ -0,0 +1,20 @@
+#ifndef AODELD_H
+#define AODELD_H
+#include "bayesnet/classifiers/Proposal.h"
+#include "bayesnet/classifiers/SPODELd.h"
+#include "Ensemble.h"
+
+namespace bayesnet {
+    class AODELd : public Ensemble, public Proposal {
+    public:
+        AODELd(bool predict_voting = true);
+        virtual ~AODELd() = default;
+        AODELd& fit(torch::Tensor& X_, torch::Tensor& y_, const std::vector<std::string>& features_, const std::string& className_, map<std::string, std::vector<int>>& states_) override;
+        void setHyperparameters(const nlohmann::json& hyperparameters) override;
+        std::vector<std::string> graph(const std::string& name = "AODELd") const override;
+    protected:
+        void trainModel(const torch::Tensor& weights) override;
+        void buildModel(const torch::Tensor& weights) override;
+    };
+}
+#endif // !AODELD_H

bayesnet/ensembles/BoostAODE.cc

@@ -0,0 +1,296 @@
+#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"
+
+namespace bayesnet {
+    struct {
+        std::string CFS = "CFS";
+        std::string FCBF = "FCBF";
+        std::string IWSS = "IWSS";
+    }SelectFeatures;
+    struct {
+        std::string ASC = "asc";
+        std::string DESC = "desc";
+        std::string RAND = "rand";
+    }Orders;
+    BoostAODE::BoostAODE(bool predict_voting) : Ensemble(predict_voting)
+    {
+        validHyperparameters = {
+            "repeatSparent", "maxModels", "order", "convergence", "threshold",
+            "select_features", "tolerance", "predict_voting", "predict_single"
+        };
+
+    }
+    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);
+            dataset_ = torch::clone(dataset);
+            // save input dataset
+            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();
+            metrics = Metrics(dataset, features, className, n_classes);
+            // Build dataset with train data
+            buildDataset(y_train);
+        } 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("repeatSparent")) {
+            repeatSparent = hyperparameters["repeatSparent"];
+            hyperparameters.erase("repeatSparent");
+        }
+        if (hyperparameters.contains("maxModels")) {
+            maxModels = hyperparameters["maxModels"];
+            hyperparameters.erase("maxModels");
+        }
+        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("predict_single")) {
+            predict_single = hyperparameters["predict_single"];
+            hyperparameters.erase("predict_single");
+        }
+        if (hyperparameters.contains("threshold")) {
+            threshold = hyperparameters["threshold"];
+            hyperparameters.erase("threshold");
+        }
+        if (hyperparameters.contains("tolerance")) {
+            tolerance = hyperparameters["tolerance"];
+            hyperparameters.erase("tolerance");
+        }
+        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.empty()) {
+            throw std::invalid_argument("Invalid hyperparameters" + hyperparameters.dump());
+        }
+    }
+    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::unordered_set<int> BoostAODE::initializeModels()
+    {
+        std::unordered_set<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();
+        for (const int& feature : cfsFeatures) {
+            featuresUsed.insert(feature);
+            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);
+            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;
+    }
+    torch::Tensor BoostAODE::ensemble_predict(torch::Tensor& X, SPODE* model)
+    {
+        if (initialize_prob_table) {
+            initialize_prob_table = false;
+            prob_table = model->predict_proba(X) * 1.0;
+        } else {
+            prob_table += model->predict_proba(X) * 1.0;
+        }
+        // prob_table doesn't store probabilities but the sum of them
+        // to have them we need to divide by the sum of the "weights" used to 
+        // consider the results obtanined in the model's predict_proba.
+        return prob_table.argmax(1);
+    }
+    void BoostAODE::trainModel(const torch::Tensor& weights)
+    {
+        // Algorithm based on the adaboost algorithm for classification
+        // as explained in Ensemble methods (Zhi-Hua Zhou, 2012)
+        initialize_prob_table = true;
+        fitted = true;
+        double alpha_t = 0;
+        torch::Tensor weights_ = torch::full({ m }, 1.0 / m, torch::kFloat64);
+        bool exitCondition = false;
+        std::unordered_set<int> featuresUsed;
+        if (selectFeatures) {
+            featuresUsed = initializeModels();
+            auto ypred = predict(X_train);
+            std::tie(weights_, alpha_t, exitCondition) = update_weights(y_train, ypred, weights_);
+            // Update significance of the models
+            for (int i = 0; i < n_models; ++i) {
+                significanceModels[i] = alpha_t;
+            }
+            if (exitCondition) {
+                return;
+            }
+        }
+        bool resetMaxModels = false;
+        if (maxModels == 0) {
+            maxModels = .1 * n > 10 ? .1 * n : n;
+            resetMaxModels = true; // Flag to unset maxModels
+        }
+        // Variables to control the accuracy finish condition
+        double priorAccuracy = 0.0;
+        double delta = 1.0;
+        double convergence_threshold = 1e-4;
+        int count = 0; // number of times the accuracy is lower than the convergence_threshold
+        // Step 0: Set the finish condition
+        // if not repeatSparent a finish condition is run out of features
+        // n_models == maxModels
+        // epsilon sub t > 0.5 => inverse the weights policy
+        // validation error is not decreasing
+        bool ascending = order_algorithm == Orders.ASC;
+        std::mt19937 g{ 173 };
+        while (!exitCondition) {
+            // Step 1: Build ranking with mutual information
+            auto featureSelection = metrics.SelectKBestWeighted(weights_, ascending, n); // Get all the features sorted
+            if (order_algorithm == Orders.RAND) {
+                std::shuffle(featureSelection.begin(), featureSelection.end(), g);
+            }
+            auto feature = featureSelection[0];
+            if (!repeatSparent || featuresUsed.size() < featureSelection.size()) {
+                bool used = true;
+                for (const auto& feat : featureSelection) {
+                    if (std::find(featuresUsed.begin(), featuresUsed.end(), feat) != featuresUsed.end()) {
+                        continue;
+                    }
+                    used = false;
+                    feature = feat;
+                    break;
+                }
+                if (used) {
+                    exitCondition = true;
+                    continue;
+                }
+            }
+            std::unique_ptr<Classifier> model;
+            model = std::make_unique<SPODE>(feature);
+            model->fit(dataset, features, className, states, weights_);
+            torch::Tensor ypred;
+            if (predict_single) {
+                ypred = model->predict(X_train);
+            } else {
+                ypred = ensemble_predict(X_train, dynamic_cast<SPODE*>(model.get()));
+            }
+            // Step 3.1: Compute the classifier amout of say
+            std::tie(weights_, alpha_t, exitCondition) = update_weights(y_train, ypred, weights_);
+            if (exitCondition) {
+                break;
+            }
+            // Step 3.4: Store classifier and its accuracy to weigh its future vote
+            featuresUsed.insert(feature);
+            models.push_back(std::move(model));
+            significanceModels.push_back(alpha_t);
+            n_models++;
+            if (convergence) {
+                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 {
+                    delta = accuracy - priorAccuracy;
+                }
+                if (delta < convergence_threshold) {
+                    count++;
+                }
+                priorAccuracy = accuracy;
+            }
+            exitCondition = n_models >= maxModels && repeatSparent || count > tolerance;
+        }
+        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));
+        if (resetMaxModels) {
+            maxModels = 0;
+        }
+    }
+    std::vector<std::string> BoostAODE::graph(const std::string& title) const
+    {
+        return Ensemble::graph(title);
+    }
+}

bayesnet/ensembles/BoostAODE.h

@@ -0,0 +1,37 @@
+#ifndef BOOSTAODE_H
+#define BOOSTAODE_H
+#include <map>
+#include "bayesnet/classifiers/SPODE.h"
+#include "bayesnet/feature_selection/FeatureSelect.h"
+#include "Ensemble.h"
+namespace bayesnet {
+    class BoostAODE : public Ensemble {
+    public:
+        BoostAODE(bool predict_voting = true);
+        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::unordered_set<int> initializeModels();
+        torch::Tensor ensemble_predict(torch::Tensor& X, SPODE* model);
+        torch::Tensor dataset_;
+        torch::Tensor X_train, y_train, X_test, y_test;
+        // Hyperparameters
+        bool repeatSparent = false; // if true, a feature can be selected more than once
+        int maxModels = 0;
+        int tolerance = 0;
+        bool predict_single = true; // wether the last model is used to predict in training or the whole ensemble
+        std::string order_algorithm; // order to process the KBest features asc, desc, rand
+        bool convergence = false; //if true, stop when the model does not improve
+        bool selectFeatures = false; // if true, use feature selection
+        std::string select_features_algorithm = "desc"; // Selected feature selection algorithm
+        bool initialize_prob_table; // if true, initialize the prob_table with the first model (used in train)
+        torch::Tensor prob_table; // Table of probabilities for ensemble predicting if predict_single is false
+        FeatureSelect* featureSelector = nullptr;
+        double threshold = -1;
+    };
+}
+#endif

bayesnet/ensembles/Ensemble.cc

@@ -0,0 +1,216 @@
+#include "Ensemble.h"
+
+namespace bayesnet {
+
+    Ensemble::Ensemble(bool predict_voting) : Classifier(Network()), n_models(0), predict_voting(predict_voting)
+    {
+
+    };
+    const std::string ENSEMBLE_NOT_FITTED = "Ensemble has not been fitted";
+    void Ensemble::trainModel(const torch::Tensor& weights)
+    {
+        n_models = models.size();
+        for (auto i = 0; i < n_models; ++i) {
+            // fit with std::vectors
+            models[i]->fit(dataset, features, className, states);
+        }
+    }
+    std::vector<int> Ensemble::compute_arg_max(std::vector<std::vector<double>>& X)
+    {
+        std::vector<int> y_pred;
+        for (auto i = 0; i < X.size(); ++i) {
+            auto max = std::max_element(X[i].begin(), X[i].end());
+            y_pred.push_back(std::distance(X[i].begin(), max));
+        }
+        return y_pred;
+    }
+    torch::Tensor Ensemble::compute_arg_max(torch::Tensor& X)
+    {
+        auto y_pred = torch::argmax(X, 1);
+        return y_pred;
+    }
+    torch::Tensor Ensemble::voting(torch::Tensor& votes)
+    {
+        // Convert m x n_models tensor to a m x n_class_states with voting probabilities
+        auto y_pred_ = votes.accessor<int, 2>();
+        std::vector<int> y_pred_final;
+        int numClasses = states.at(className).size();
+        // votes is m x n_models with the prediction of every model for each sample
+        auto result = torch::zeros({ votes.size(0), numClasses }, torch::kFloat32);
+        auto sum = std::reduce(significanceModels.begin(), significanceModels.end());
+        for (int i = 0; i < votes.size(0); ++i) {
+            // n_votes store in each index (value of class) the significance added by each model
+            // i.e. n_votes[0] contains how much value has the value 0 of class. That value is generated by the models predictions
+            std::vector<double> n_votes(numClasses, 0.0);
+            for (int j = 0; j < n_models; ++j) {
+                n_votes[y_pred_[i][j]] += significanceModels.at(j);
+            }
+            result[i] = torch::tensor(n_votes);
+        }
+        // To only do one division and gain precision
+        result /= sum;
+        return result;
+    }
+    std::vector<std::vector<double>> Ensemble::predict_proba(std::vector<std::vector<int>>& X)
+    {
+        if (!fitted) {
+            throw std::logic_error(ENSEMBLE_NOT_FITTED);
+        }
+        return predict_voting ? predict_average_voting(X) : predict_average_proba(X);
+    }
+    torch::Tensor Ensemble::predict_proba(torch::Tensor& X)
+    {
+        if (!fitted) {
+            throw std::logic_error(ENSEMBLE_NOT_FITTED);
+        }
+        return predict_voting ? predict_average_voting(X) : predict_average_proba(X);
+    }
+    std::vector<int> Ensemble::predict(std::vector<std::vector<int>>& X)
+    {
+        auto res = predict_proba(X);
+        return compute_arg_max(res);
+    }
+    torch::Tensor Ensemble::predict(torch::Tensor& X)
+    {
+        auto res = predict_proba(X);
+        return compute_arg_max(res);
+    }
+    torch::Tensor Ensemble::predict_average_proba(torch::Tensor& X)
+    {
+        auto n_states = models[0]->getClassNumStates();
+        torch::Tensor y_pred = torch::zeros({ X.size(1), n_states }, torch::kFloat32);
+        auto threads{ std::vector<std::thread>() };
+        std::mutex mtx;
+        for (auto i = 0; i < n_models; ++i) {
+            threads.push_back(std::thread([&, i]() {
+                auto ypredict = models[i]->predict_proba(X);
+                std::lock_guard<std::mutex> lock(mtx);
+                y_pred += ypredict * significanceModels[i];
+                }));
+        }
+        for (auto& thread : threads) {
+            thread.join();
+        }
+        auto sum = std::reduce(significanceModels.begin(), significanceModels.end());
+        y_pred /= sum;
+        return y_pred;
+    }
+    std::vector<std::vector<double>> Ensemble::predict_average_proba(std::vector<std::vector<int>>& X)
+    {
+        auto n_states = models[0]->getClassNumStates();
+        std::vector<std::vector<double>> y_pred(X[0].size(), std::vector<double>(n_states, 0.0));
+        auto threads{ std::vector<std::thread>() };
+        std::mutex mtx;
+        for (auto i = 0; i < n_models; ++i) {
+            threads.push_back(std::thread([&, i]() {
+                auto ypredict = models[i]->predict_proba(X);
+                assert(ypredict.size() == y_pred.size());
+                assert(ypredict[0].size() == y_pred[0].size());
+                std::lock_guard<std::mutex> lock(mtx);
+                // Multiply each prediction by the significance of the model and then add it to the final prediction
+                for (auto j = 0; j < ypredict.size(); ++j) {
+                    std::transform(y_pred[j].begin(), y_pred[j].end(), ypredict[j].begin(), y_pred[j].begin(),
+                        [significanceModels = significanceModels[i]](double x, double y) { return x + y * significanceModels; });
+                }
+                }));
+        }
+        for (auto& thread : threads) {
+            thread.join();
+        }
+        auto sum = std::reduce(significanceModels.begin(), significanceModels.end());
+        //Divide each element of the prediction by the sum of the significances
+        for (auto j = 0; j < y_pred.size(); ++j) {
+            std::transform(y_pred[j].begin(), y_pred[j].end(), y_pred[j].begin(), [sum](double x) { return x / sum; });
+        }
+        return y_pred;
+    }
+    std::vector<std::vector<double>> Ensemble::predict_average_voting(std::vector<std::vector<int>>& X)
+    {
+        torch::Tensor Xt = bayesnet::vectorToTensor(X, false);
+        auto y_pred = predict_average_voting(Xt);
+        std::vector<std::vector<double>> result = tensorToVectorDouble(y_pred);
+        return result;
+    }
+    torch::Tensor Ensemble::predict_average_voting(torch::Tensor& X)
+    {
+        // Build a m x n_models tensor with the predictions of each model
+        torch::Tensor y_pred = torch::zeros({ X.size(1), n_models }, torch::kInt32);
+        auto threads{ std::vector<std::thread>() };
+        std::mutex mtx;
+        for (auto i = 0; i < n_models; ++i) {
+            threads.push_back(std::thread([&, i]() {
+                auto ypredict = models[i]->predict(X);
+                std::lock_guard<std::mutex> lock(mtx);
+                y_pred.index_put_({ "...", i }, ypredict);
+                }));
+        }
+        for (auto& thread : threads) {
+            thread.join();
+        }
+        return voting(y_pred);
+    }
+    float Ensemble::score(torch::Tensor& X, torch::Tensor& y)
+    {
+        auto y_pred = predict(X);
+        int correct = 0;
+        for (int i = 0; i < y_pred.size(0); ++i) {
+            if (y_pred[i].item<int>() == y[i].item<int>()) {
+                correct++;
+            }
+        }
+        return (double)correct / y_pred.size(0);
+    }
+    float Ensemble::score(std::vector<std::vector<int>>& X, std::vector<int>& y)
+    {
+        auto y_pred = predict(X);
+        int correct = 0;
+        for (int i = 0; i < y_pred.size(); ++i) {
+            if (y_pred[i] == y[i]) {
+                correct++;
+            }
+        }
+        return (double)correct / y_pred.size();
+    }
+    std::vector<std::string> Ensemble::show() const
+    {
+        auto result = std::vector<std::string>();
+        for (auto i = 0; i < n_models; ++i) {
+            auto res = models[i]->show();
+            result.insert(result.end(), res.begin(), res.end());
+        }
+        return result;
+    }
+    std::vector<std::string> Ensemble::graph(const std::string& title) const
+    {
+        auto result = std::vector<std::string>();
+        for (auto i = 0; i < n_models; ++i) {
+            auto res = models[i]->graph(title + "_" + std::to_string(i));
+            result.insert(result.end(), res.begin(), res.end());
+        }
+        return result;
+    }
+    int Ensemble::getNumberOfNodes() const
+    {
+        int nodes = 0;
+        for (auto i = 0; i < n_models; ++i) {
+            nodes += models[i]->getNumberOfNodes();
+        }
+        return nodes;
+    }
+    int Ensemble::getNumberOfEdges() const
+    {
+        int edges = 0;
+        for (auto i = 0; i < n_models; ++i) {
+            edges += models[i]->getNumberOfEdges();
+        }
+        return edges;
+    }
+    int Ensemble::getNumberOfStates() const
+    {
+        int nstates = 0;
+        for (auto i = 0; i < n_models; ++i) {
+            nstates += models[i]->getNumberOfStates();
+        }
+        return nstates;
+    }
+}

bayesnet/ensembles/Ensemble.h

@@ -0,0 +1,46 @@
+#ifndef ENSEMBLE_H
+#define ENSEMBLE_H
+#include <torch/torch.h>
+#include "bayesnet/utils/BayesMetrics.h"
+#include "bayesnet/utils/bayesnetUtils.h"
+#include "bayesnet/classifiers/Classifier.h"
+
+namespace bayesnet {
+    class Ensemble : public Classifier {
+    public:
+        Ensemble(bool predict_voting = true);
+        virtual ~Ensemble() = default;
+        torch::Tensor predict(torch::Tensor& X) override;
+        std::vector<int> predict(std::vector<std::vector<int>>& X) override;
+        torch::Tensor predict_proba(torch::Tensor& X) override;
+        std::vector<std::vector<double>> predict_proba(std::vector<std::vector<int>>& X) override;
+        float score(torch::Tensor& X, torch::Tensor& y) override;
+        float score(std::vector<std::vector<int>>& X, std::vector<int>& y) override;
+        int getNumberOfNodes() const override;
+        int getNumberOfEdges() const override;
+        int getNumberOfStates() const override;
+        std::vector<std::string> show() const override;
+        std::vector<std::string> graph(const std::string& title) const override;
+        std::vector<std::string> topological_order()  override
+        {
+            return std::vector<std::string>();
+        }
+        void dump_cpt() const override
+        {
+        }
+    protected:
+        torch::Tensor predict_average_voting(torch::Tensor& X);
+        std::vector<std::vector<double>> predict_average_voting(std::vector<std::vector<int>>& X);
+        torch::Tensor predict_average_proba(torch::Tensor& X);
+        std::vector<std::vector<double>> predict_average_proba(std::vector<std::vector<int>>& X);
+        torch::Tensor compute_arg_max(torch::Tensor& X);
+        std::vector<int> compute_arg_max(std::vector<std::vector<double>>& X);
+        torch::Tensor voting(torch::Tensor& votes);
+        unsigned n_models;
+        std::vector<std::unique_ptr<Classifier>> models;
+        std::vector<double> significanceModels;
+        void trainModel(const torch::Tensor& weights) override;
+        bool predict_voting;
+    };
+}
+#endif