Repository: automl/random_forest_run (language: python)

Issue:
Predictions in Log-Cost-Space
I hope that I correctly remember our discussion yesterday about the predictions in log-cost space since I forgot my notes in my office. @frank-hutter if anything is wrong, please correct me.

@sfalkner Frank explained yesterday how he implemented the prediction in log(cost) space in SMAC and I don't know whether this is right now possible with the new RF. I hope you can please help us here. 
* Train the RF using log(cost) values
* to get a marginalized prediction over instances
  1. compute a marginalized prediction for each tree using exp(log(cost)) of all values in the leafs -> one prediction for each tree in the original cost space
  1. mean and variance over all log(pred_t) for each t (so, again in the log-space)

How can we compute this with the RF? Is it possible using the python interface? Would it be inefficient to do it in Python? Can it be done within C++?


Discussion / hints:
IIRC, you marginalize over instances yourself right now anyway. You can use the `all_leaf_values` method to get the actual values stored in the corresponding leaf of each tree and iterate over that to compute anything you like. I don't quite understand why you want to have the final mean and variance prediction in log-space again, though.
I don't think you would gain a lot doing this in C++. If that is a major use case and you want to do a lot of predictions with that, it would be more efficient to handle the transforms during fitting such that the marginalization is fast. That would require some C++ coding, but could be done if that turns out to improve your model quality and too slow in python.
My concern still is the constant change from log to non-log space and how that affects the RF predictions. You still train it on log data, so you assume a log-normal distribution, but you want to marginalize the 'normal' values of which you will take the log again. I don't know if that's what you actually want...
I forgot to mention, I updated the `tests/pyrfr_example.py` file to show how the `all_leaf_values` method is used.
Thank you for the thorough explanation. Right now, my goal is to have reimplementation of the old SMAC. At some point, we should discuss with Frank whether this is the best way to do it or how we can evaluate alternatives.

Fix patch:
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 0ccd960..cf5bf2d 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -4,7 +4,7 @@ project (RFR CXX)
 # add a version number
 set (RFR_VERSION_MAJOR 0)
 set (RFR_VERSION_MINOR 8)
-set (RFR_VERSION_RELEASE 3)
+set (RFR_VERSION_RELEASE 4)
 
 # Print debug information about boost
 set(Boost_DEBUG 1)
@@ -88,9 +88,9 @@ add_subdirectory(tests)
 
 
 IF(CMAKE_COMPILER_IS_GNUCXX)
-  ADD_DEFINITIONS("-Wall -O0 -g")
+  ADD_DEFINITIONS("-Wall -O3 -g")
 ELSEIF(${CMAKE_CXX_COMPILER_ID} MATCHES "Clang")
-  ADD_DEFINITIONS("-Wall -fPIC -pg -std=c++11")
+  ADD_DEFINITIONS("-Wall -fPIC -pg -O3 -std=c++11")
 ELSE()
   MESSAGE(FATAL_ERROR "CMakeLists.txt has not been tested/written for your compiler.")
 ENDIF()
diff --git a/include/rfr/forests/regression_forest.hpp b/include/rfr/forests/regression_forest.hpp
index 78b60a1..0c2744a 100644
--- a/include/rfr/forests/regression_forest.hpp
+++ b/include/rfr/forests/regression_forest.hpp
@@ -24,453 +24,538 @@
 #include <sstream>
 
 
-
 #include "rfr/trees/tree_options.hpp"
 #include "rfr/forests/forest_options.hpp"
 #include "rfr/util.hpp"
 
-namespace rfr{ namespace forests{
+namespace rfr {
+    namespace forests {
+
+        typedef cereal::PortableBinaryInputArchive binary_iarch_t;
+        typedef cereal::PortableBinaryOutputArchive binary_oarch_t;
+
+        typedef cereal::JSONInputArchive ascii_iarch_t;
+        typedef cereal::JSONOutputArchive ascii_oarch_t;
+
+
+        template<typename tree_type, typename num_t = float, typename response_t = float, typename index_t = unsigned int, typename rng_type=std::default_random_engine>
+        class regression_forest {
+        protected:
+            std::vector<tree_type> the_trees;
+            index_t num_features;
+
+            std::vector<std::vector<num_t> > bootstrap_sample_weights;
 
-typedef cereal::PortableBinaryInputArchive binary_iarch_t;
-typedef cereal::PortableBinaryOutputArchive binary_oarch_t;
+            num_t oob_error = NAN;
 
-typedef cereal::JSONInputArchive ascii_iarch_t;
-typedef cereal::JSONOutputArchive ascii_oarch_t;
+            // the forest needs to remember the data types on which it was trained
+            std::vector<index_t> types;
+            std::vector<std::array<num_t, 2> > bounds;
 
 
+        public:
 
+            forest_options<num_t, response_t, index_t> options;
 
-template <typename tree_type, typename num_t = float, typename response_t = float, typename index_t = unsigned int,  typename rng_type=std::default_random_engine>
-class regression_forest{
-  protected:
-	std::vector<tree_type> the_trees;
-	index_t num_features;
 
-	std::vector<std::vector<num_t> > bootstrap_sample_weights;
-	
-	num_t oob_error = NAN;
-	
-	// the forest needs to remember the data types on which it was trained
-	std::vector<index_t> types;
-	std::vector< std::array<num_t,2> > bounds;
-	
+            /** \brief serialize function for saving forests with cerial*/
+            template<class Archive>
+            void serialize(Archive &archive) {
+              archive(options, the_trees, num_features, bootstrap_sample_weights, oob_error, types, bounds);
+            }
 
-  public:
+            regression_forest() : options() {}
 
-	forest_options<num_t, response_t, index_t> options;
+            regression_forest(forest_options<num_t, response_t, index_t> opts) : options(opts) {}
 
+            virtual ~regression_forest() {};
 
-  	/** \brief serialize function for saving forests with cerial*/
-  	template<class Archive>
-	void serialize(Archive & archive)
-	{
-		archive( options, the_trees, num_features, bootstrap_sample_weights, oob_error, types, bounds);
-	}
+            /**\brief growing the random forest for a given data set
+             *
+             * \param data a filled data container
+             * \param rng the random number generator to be used
+             */
+            virtual void fit(const rfr::data_containers::base<num_t, response_t, index_t> &data, rng_type &rng) {
 
-	regression_forest(): options()	{}
-	
-	regression_forest(forest_options<num_t, response_t, index_t> opts): options(opts){}
+              if (options.num_trees <= 0)
+                throw std::runtime_error("The number of trees has to be positive!");
 
-	virtual ~regression_forest()	{};
+              if ((!options.do_bootstrapping) && (data.num_data_points() < options.num_data_points_per_tree))
+                throw std::runtime_error(
+                        "You cannot use more data points per tree than actual data point present without bootstrapping!");
 
-	/**\brief growing the random forest for a given data set
-	 * 
-	 * \param data a filled data container
-	 * \param rng the random number generator to be used
-	 */
-	virtual void fit(const rfr::data_containers::base<num_t, response_t, index_t> &data, rng_type &rng){
 
-		if (options.num_trees <= 0)
-			throw std::runtime_error("The number of trees has to be positive!");
+              types.resize(data.num_features());
+              bounds.resize(data.num_features());
+              for (auto i = 0u; i < data.num_features(); ++i) {
+                types[i] = data.get_type_of_feature(i);
+                auto p = data.get_bounds_of_feature(i);
+                bounds[i][0] = p.first;
+                bounds[i][1] = p.second;
+              }
 
-		if ((!options.do_bootstrapping) && (data.num_data_points() < options.num_data_points_per_tree))
-			throw std::runtime_error("You cannot use more data points per tree than actual data point present without bootstrapping!");
+              the_trees.resize(options.num_trees);
 
 
-		types.resize(data.num_features());
-		bounds.resize(data.num_features());
-		for (auto i=0u; i<data.num_features(); ++i){
-			types[i] = data.get_type_of_feature(i);
-			auto p = data.get_bounds_of_feature(i);
-			bounds[i][0] = p.first;
-			bounds[i][1] = p.second;
-		}
+              std::vector<index_t> data_indices(data.num_data_points());
+              std::iota(data_indices.begin(), data_indices.end(), 0);
 
-		the_trees.resize(options.num_trees);
+              num_features = data.num_features();
 
+              // catch some stupid things that will make the forest crash when fitting
+              if (options.num_data_points_per_tree == 0)
+                throw std::runtime_error("The number of data points per tree is set to zero!");
 
-		std::vector<index_t> data_indices( data.num_data_points());
-		std::iota(data_indices.begin(), data_indices.end(), 0);
+              if (options.tree_opts.max_features == 0)
+                throw std::runtime_error("The number of features used for a split is set to zero!");
 
-		num_features = data.num_features();
-		
-		// catch some stupid things that will make the forest crash when fitting
-		if (options.num_data_points_per_tree == 0)
-			throw std::runtime_error("The number of data points per tree is set to zero!");
-		
-		if (options.tree_opts.max_features == 0)
-			throw std::runtime_error("The number of features used for a split is set to zero!");
-		
-		bootstrap_sample_weights.clear();
+              bootstrap_sample_weights.clear();
 
-		for (auto &tree : the_trees){
-            std::vector<num_t> bssf (data.num_data_points(), 0); // BootStrap Sample Frequencies
-			// prepare the data(sub)set
-			if (options.do_bootstrapping){
-                std::uniform_int_distribution<index_t> dist (0,data.num_data_points()-1);
-                for (auto i=0u; i < options.num_data_points_per_tree; ++i){
-					bssf[dist(rng)]+=1;
-				}
-			}
-			else{
-				std::shuffle(data_indices.begin(), data_indices.end(), rng);
-                for (auto i=0u; i < options.num_data_points_per_tree; ++i)
+              for (auto &tree: the_trees) {
+                std::vector<num_t> bssf(data.num_data_points(), 0); // BootStrap Sample Frequencies
+                // prepare the data(sub)set
+                if (options.do_bootstrapping) {
+                  std::uniform_int_distribution<index_t> dist(0, data.num_data_points() - 1);
+                  for (auto i = 0u; i < options.num_data_points_per_tree; ++i) {
+                    bssf[dist(rng)] += 1;
+                  }
+                } else {
+                  std::shuffle(data_indices.begin(), data_indices.end(), rng);
+                  for (auto i = 0u; i < options.num_data_points_per_tree; ++i)
                     bssf[data_indices[i]] += 1;
-			}
-			
-			tree.fit(data, options.tree_opts, bssf, rng);
-			
-			// record sample counts for later use
-			if (options.compute_oob_error)
-				bootstrap_sample_weights.push_back(bssf);
-		}
-		
-		oob_error = NAN;
-		
-		if (options.compute_oob_error){
-			
-			rfr::util::running_statistics<num_t> oob_error_stat;
-			
-			for (auto i=0u; i < data.num_data_points(); i++){
-
-				rfr::util::running_statistics<num_t> prediction_stat;
-
-				for (auto j=0u; j<the_trees.size(); j++){
-					// only consider data points that were not part of that bootstrap sample
-					if (bootstrap_sample_weights[j][i] == 0)
-						prediction_stat.push(the_trees[j].predict( data.retrieve_data_point(i)));
-				}
-				
-				// compute squared error of prediction
-
-				if (prediction_stat.number_of_points() > 0u){
-					oob_error_stat.push(std::pow(prediction_stat.mean() - data.response(i), (num_t) 2));
-				}
-			}
-			oob_error = std::sqrt(oob_error_stat.mean());
-		}
-	}
-
-
-	/* \brief combines the prediction of all trees in the forest
-	 *
-	 * Every random tree makes an individual prediction which are averaged for the forest's prediction.
-	 *
-	 * \param feature_vector a valid vector containing the features
-	 * \return response_t the predicted value
-	 */
-    response_t predict( const std::vector<num_t> &feature_vector) const{
-
-		// collect the predictions of individual trees
-		rfr::util::running_statistics<num_t> mean_stats;
-		for (auto &tree: the_trees)
-			mean_stats.push(tree.predict(feature_vector));
-		return(mean_stats.mean());
-	}
-    
-    
-   /* \brief makes a prediction for the mean and a variance estimation
-    * 
-    * Every tree returns the mean and the variance of the leaf the feature vector falls into.
-    * These are combined to the forests mean prediction (mean of the means) and a variance estimate
-    * (mean of the variance + variance of the means).
-    * 
-    * Use weighted_data = false if the weights assigned to each data point were frequencies, not importance weights.
-    * Use this if you haven't assigned any weigths, too.
-    * 
-	* \param feature_vector a valid feature vector
-	* \param weighted_data whether the data had importance weights
-	* \return std::pair<response_t, num_t> mean and variance prediction
-    */
-    std::pair<num_t, num_t> predict_mean_var( const std::vector<num_t> &feature_vector, bool weighted_data = false){
-
-		// collect the predictions of individual trees
-		rfr::util::running_statistics<num_t> mean_stats, var_stats;
-		for (auto &tree: the_trees){
-			auto stat = tree.leaf_statistic(feature_vector);
-			mean_stats.push(stat.mean());
-			if (stat.number_of_points() > 1){
-				if (weighted_data) var_stats.push(stat.variance_unbiased_importance());
-				else var_stats.push(stat.variance_unbiased_frequency());
-			} else{
-				var_stats.push(0);
-			}
-		}
-	    num_t var = mean_stats.variance_sample();
-		if (options.compute_law_of_total_variance) {
-			return std::pair<num_t, num_t> (mean_stats.mean(), std::max<num_t>(0, var + var_stats.mean()) );
-		}
-		return std::pair<num_t, num_t> (mean_stats.mean(), std::max<num_t>(0, var) );
-	}
-
-
-	/* \brief predict the mean and the variance deviation for a configuration marginalized over a given set of partial configurations
-	 * 
-	 * This function will be mostly used to predict the mean over a given set of instances, but could be used to marginalize over any discrete set of partial configurations.
-	 * 
-	 * \param features a (partial) configuration where unset values should be set to NaN
-	 * \param set_features a array containing the (partial) assignments used for the averaging. Every NaN value will be replaced by the corresponding value from features.
-	 * \param set_size number of feature vectors in set_features
-	 * 
-	 * \return std::pair<num_t, num_t> mean and variance prediction of a feature vector averaged over 
-	 */
-    /*
-	std::pair<num_t, num_t> predict_mean_var_marginalized_over_set (num_t *features, num_t* set_features, index_t set_size){
-		
-		num_t fv[num_features];
-
-		// collect the predictions of individual trees
-		rfr::util::running_statistics<num_t> mean_stats, var_stats;
-		for (auto i=0u; i < set_size; ++i){
-			// construct the actual feature vector
-			rfr::util::merge_two_vectors(features, &set_features[i*num_features], fv, num_features);
-
-			num_t m , v;
-			std::tie(m, v) = predict_mean_var(fv);
-
-			mean_stats(m);
-			var_stats(v);
-		}
-		return(std::pair<num_t, num_t> (mean_stats.mean(), std::max<num_t>(0, mean_stats.variance() + var_stats.mean()) ));
-	}
-    */
-
-	/* \brief predict the mean and the variance of the mean prediction across a set of partial features
-	 * 
-	 * A very special function to predict the mean response of a a partial assignment for a given set.
-	 * It takes the prediction of set-mean of every individual tree and combines to estimate the mean its
-	 * total variance. The predictions of two trees are considered uncorrelated
-	 * 
-	 * \param features a (partial) configuration where unset values should be set to NaN
-	 * \param set_features a 1d-array containing the (partial) assignments used for the averaging. Every NaN value will be replaced by the corresponding value from features. The array must hold set_size times the number of features entries! There is no consistency check!
-	 * \param set_size number of feature vectors in set_features
-	 * 
-	 * \return std::tuple<num_t, num_t, num_t> mean and variance of empirical mean prediction of a feature vector averaged over. The last one is the estimated variance of a sample drawn from partial assignment.
-	 */
-    /*
-	std::tuple<num_t, num_t, num_t> predict_mean_var_of_mean_response_on_set (num_t *features, num_t* set_features, index_t set_size){
-
-			num_t fv[num_features];
-
-			rfr::util::running_statistics<num_t> mean_stats, var_stats, sample_var_stats, sample_mean_stats;
-
-			for (auto &t : the_trees){
-
-					rfr::util::running_statistics<num_t> tree_mean_stats, tree_var_stats;
-
-					for (auto i=0u; i < set_size; ++i){
-
-							rfr::util::merge_two_vectors(features, &set_features[i*num_features], fv, num_features);
-
-							num_t m , v; index_t n;
-							std::tie(m, v, n) = t.predict_mean_var_N(fv);
-
-							tree_mean_stats(m); tree_var_stats(v); sample_mean_stats(m); sample_var_stats(v);
-					}
-
-					mean_stats(tree_mean_stats.mean());
-					var_stats(std::max<num_t>(0, tree_var_stats.mean()));
-					
-			}
-			
-			return(std::make_tuple(mean_stats.mean(), std::max<num_t>(0, mean_stats.variance()) + std::max<num_t>(0, var_stats.mean()/set_size), std::max<num_t>(0,sample_mean_stats.variance() + sample_var_stats.mean())));
-	}
-    */
-
-	/* \brief estimates the covariance of two feature vectors
-	 * 
-	 * 
-	 * The covariance between to input vectors contains information about the
-	 * feature space. For other models, like GPs, this is a natural quantity
-	 * (e.g., property of the kernel). Here, we try to estimate it using the
-	 * emprical covariance of the individual tree's predictions.
-	 * 
-	 * \param f1 a valid feature vector (no sanity checks are performed!)
-	 * \param f2 a second feature vector (no sanity checks are performed!)
-	 */
-    
-	num_t covariance (const std::vector<num_t> &f1, const std::vector<num_t> &f2){
-
-		rfr::util::running_covariance<num_t> run_cov_of_means;
-
-		for (auto &t: the_trees)
-			run_cov_of_means.push(t.predict(f1),t.predict(f2));
-
-		return(run_cov_of_means.covariance());
-	}
-
-
-
-	/* \brief computes the kernel of a 'Kernel Random Forest'
-	 *
-	 * Source: "Random forests and kernel methods" by Erwan Scornet
-	 * 
-	 * The covariance between to input vectors contains information about the
-	 * feature space. For other models, like GPs, this is a natural quantity
-	 * (e.g., property of the kernel). Here, we try to estimate it using the
-	 * emprical covariance of the individual tree's predictions.
-	 * 
-	 * \param f1 a valid feature vector (no sanity checks are performed!)
-	 * \param f2 a second feature vector (no sanity checks are performed!)
-	 */
-    
-	num_t kernel (const std::vector<num_t> &f1, const std::vector<num_t> &f2){
-
-		rfr::util::running_statistics<num_t> stat;
-
-		for (auto &t: the_trees){
-			auto l1 = t.find_leaf_index(f1);
-			auto l2 = t.find_leaf_index(f2);
-
-			stat.push(l1==l2);
-		}
-		return(stat.mean());
-	}
-
-	
-	std::vector<tree_type> get_all_trees() const {return the_trees;}
-
-
-	std::vector< std::vector<num_t> > all_leaf_values (const std::vector<num_t> &feature_vector) const {
-		std::vector< std::vector<num_t> > rv;
-		rv.reserve(the_trees.size());
-
-		for (auto &t: the_trees){
-			rv.push_back(t.leaf_entries(feature_vector));
-		}
-		return(rv);
-	}
-
-
-	
-	/* \brief returns the predictions of every tree marginalized over the NAN values in the feature_vector
-	 * 
-	 * TODO: more documentation over how the 'missing values' are handled
-	 * 
-	 * \param feature_vector non-specfied values (NaN) will be marginalized over according to the training data
-	 */
-	//std::vector<num_t> marginalized_mean_predictions(const std::vector<num_t> &feature_vector) const {
-	//	std::vector<num_t> rv;
-	//	rv.reserve(the_trees.size());
-	//	for (auto &t : the_trees)
-	//		rv.emplace_back(t.marginalized_mean_prediction(feature_vector));
-	//	return(rv);
-	//}
-
-
-
-	/* \brief updates the forest by adding the provided datapoint without a complete retraining
-	 * 
-	 * 
-	 * As retraining can be quite expensive, this function can be used to quickly update the forest
-	 * by finding the leafs the datapoints belong into and just inserting them. This is, of course,
-	 * not the right way to do it for many data points, but it should be a good approximation for a few.
-	 * 
-	 * \param features a valid feature vector
-	 * \param response the corresponding response value
-	 * \param weight the associated weight
-	 */
-	void pseudo_update (std::vector<num_t> features, response_t response, num_t weight){
-		for (auto &t: the_trees)
-			t.pseudo_update(features, response, weight);
-	}
-	
-	/* \brief undoing a pseudo update by removing a point
-	 * 
-	 * This function removes one point from the corresponding leaves into
-	 * which the given feature vector falls
-	 * 
-	 * \param features a valid feature vector
-	 * \param response the corresponding response value
-	 * \param weight the associated weight
-	 */
-	void pseudo_downdate(std::vector<num_t> features, response_t response, num_t weight){
-		for (auto &t: the_trees)
-			t.pseudo_downdate(features, response, weight);
-	}
-	
-	num_t out_of_bag_error(){return(oob_error);}
-
-	/* \brief writes serialized representation into a binary file
-	 * 
-	 * \param filename name of the file to store the forest in. Make sure that the directory exists!
-	 */
-	void save_to_binary_file(const std::string filename){
-		std::ofstream ofs(filename, std::ios::binary);
-		binary_oarch_t oarch(ofs);
-		serialize(oarch);
-	}
-
-	/* \brief deserialize from a binary file created by save_to_binary_file
-	 *
-	 * \param filename name of the file in which the forest is stored. 
-	 */
-	void load_from_binary_file(const std::string filename){
-		std::ifstream ifs(filename, std::ios::binary);
-		binary_iarch_t iarch(ifs);
-		serialize(iarch);
-	}
-
-	/* serialize into a string; used for Python's pickle.dump
-	 * 
-	 * \return std::string a JSON serialization of the forest
-	 */
-	std::string ascii_string_representation(){
-		std::stringstream oss;
-		{
-			ascii_oarch_t oarch(oss);
-			serialize(oarch);
-		}
-		return(oss.str());
-	}
-
-	/* \brief deserialize from string; used for Python's pickle.load
-	 * 
-	 * \return std::string a JSON serialization of the forest
-	 */
-	void load_from_ascii_string( std::string const &str){
-		std::stringstream iss;
-		iss.str(str);
-		ascii_iarch_t iarch(iss);
-		serialize(iarch);
-	}
-
-
-
-	/* \brief stores a latex document for every individual tree
-	 * 
-	 * \param filename_template a string to specify the location and the naming scheme. Note the directory is not created, so make sure it exists.
-	 * 
-	 */
-	void save_latex_representation(const std::string filename_template){
-		for (auto i = 0u; i<the_trees.size(); i++){
-			std::stringstream filename;
-			filename << filename_template<<i<<".tex";
-			the_trees[i].save_latex_representation(filename.str().c_str());
-		}
-	}
-
-	void print_info(){
-		for (auto t: the_trees){
-			t.print_info();
-		}
-	}
-
-
-	virtual unsigned int num_trees (){ return(the_trees.size());}
-	
-};
-
-
-}}//namespace rfr::forests
+                }
+
+                tree.fit(data, options.tree_opts, bssf, rng);
+
+                // record sample counts for later use
+                if (options.compute_oob_error)
+                  bootstrap_sample_weights.push_back(bssf);
+              }
+
+              oob_error = NAN;
+
+              if (options.compute_oob_error) {
+
+                rfr::util::running_statistics<num_t> oob_error_stat;
+
+                for (auto i = 0u; i < data.num_data_points(); i++) {
+
+                  rfr::util::running_statistics<num_t> prediction_stat;
+
+                  for (auto j = 0u; j < the_trees.size(); j++) {
+                    // only consider data points that were not part of that bootstrap sample
+                    if (bootstrap_sample_weights[j][i] == 0)
+                      prediction_stat.push(the_trees[j].predict(data.retrieve_data_point(i)));
+                  }
+
+                  // compute squared error of prediction
+
+                  if (prediction_stat.number_of_points() > 0u) {
+                    oob_error_stat.push(std::pow(prediction_stat.mean() - data.response(i), (num_t) 2));
+                  }
+                }
+                oob_error = std::sqrt(oob_error_stat.mean());
+              }
+            }
+
+
+            /* \brief combines the prediction of all trees in the forest
+             *
+             * Every random tree makes an individual prediction which are averaged for the forest's prediction.
+             *
+             * \param feature_vector a valid vector containing the features
+             * \return response_t the predicted value
+             */
+            response_t predict(const std::vector<num_t> &feature_vector) const {
+
+              // collect the predictions of individual trees
+              rfr::util::running_statistics<num_t> mean_stats;
+              for (auto &tree: the_trees)
+                mean_stats.push(tree.predict(feature_vector));
+              return (mean_stats.mean());
+            }
+
+
+            /* \brief makes a prediction for the mean and a variance estimation
+             *
+             * Every tree returns the mean and the variance of the leaf the feature vector falls into.
+             * These are combined to the forests mean prediction (mean of the means) and a variance estimate
+             * (mean of the variance + variance of the means).
+             *
+             * Use weighted_data = false if the weights assigned to each data point were frequencies, not importance weights.
+             * Use this if you haven't assigned any weigths, too.
+             *
+             * \param feature_vector a valid feature vector
+             * \param weighted_data whether the data had importance weights
+             * \return std::pair<response_t, num_t> mean and variance prediction
+             */
+            std::pair<num_t, num_t>
+            predict_mean_var(const std::vector<num_t> &feature_vector, bool weighted_data = false) {
+
+              // collect the predictions of individual trees
+              rfr::util::running_statistics<num_t> mean_stats, var_stats;
+              for (auto &tree: the_trees) {
+                auto stat = tree.leaf_statistic(feature_vector);
+                mean_stats.push(stat.mean());
+                if (stat.number_of_points() > 1) {
+                  if (weighted_data) var_stats.push(stat.variance_unbiased_importance());
+                  else var_stats.push(stat.variance_unbiased_frequency());
+                } else {
+                  var_stats.push(0);
+                }
+              }
+              num_t var = mean_stats.variance_sample();
+              if (options.compute_law_of_total_variance) {
+                return std::pair<num_t, num_t>(mean_stats.mean(), std::max<num_t>(0, var + var_stats.mean()));
+              }
+              return std::pair<num_t, num_t>(mean_stats.mean(), std::max<num_t>(0, var));
+            }
+
+
+            /* \brief predict the mean and the variance deviation for a configuration marginalized over a given set of partial configurations
+             *
+             * This function will be mostly used to predict the mean over a given set of instances, but could be used to marginalize over any discrete set of partial configurations.
+             *
+             * \param features a (partial) configuration where unset values should be set to NaN
+             * \param set_features a array containing the (partial) assignments used for the averaging. Every NaN value will be replaced by the corresponding value from features.
+             * \param set_size number of feature vectors in set_features
+             *
+             * \return std::pair<num_t, num_t> mean and variance prediction of a feature vector averaged over
+             */
+            /*
+            std::pair<num_t, num_t> predict_mean_var_marginalized_over_set (num_t *features, num_t* set_features, index_t set_size){
+
+                num_t fv[num_features];
+
+                // collect the predictions of individual trees
+                rfr::util::running_statistics<num_t> mean_stats, var_stats;
+                for (auto i=0u; i < set_size; ++i){
+                    // construct the actual feature vector
+                    rfr::util::merge_two_vectors(features, &set_features[i*num_features], fv, num_features);
+
+                    num_t m , v;
+                    std::tie(m, v) = predict_mean_var(fv);
+
+                    mean_stats(m);
+                    var_stats(v);
+                }
+                return(std::pair<num_t, num_t> (mean_stats.mean(), std::max<num_t>(0, mean_stats.variance() + var_stats.mean()) ));
+            }
+            */
+
+            /* \brief predict the mean and the variance of the mean prediction across a set of partial features
+             *
+             * A very special function to predict the mean response of a a partial assignment for a given set.
+             * It takes the prediction of set-mean of every individual tree and combines to estimate the mean its
+             * total variance. The predictions of two trees are considered uncorrelated
+             *
+             * \param features a (partial) configuration where unset values should be set to NaN
+             * \param set_features a 1d-array containing the (partial) assignments used for the averaging. Every NaN value will be replaced by the corresponding value from features. The array must hold set_size times the number of features entries! There is no consistency check!
+             * \param set_size number of feature vectors in set_features
+             *
+             * \return std::tuple<num_t, num_t, num_t> mean and variance of empirical mean prediction of a feature vector averaged over. The last one is the estimated variance of a sample drawn from partial assignment.
+             */
+            /*
+            std::tuple<num_t, num_t, num_t> predict_mean_var_of_mean_response_on_set (num_t *features, num_t* set_features, index_t set_size){
+
+                    num_t fv[num_features];
+
+                    rfr::util::running_statistics<num_t> mean_stats, var_stats, sample_var_stats, sample_mean_stats;
+
+                    for (auto &t : the_trees){
+
+                            rfr::util::running_statistics<num_t> tree_mean_stats, tree_var_stats;
+
+                            for (auto i=0u; i < set_size; ++i){
+
+                                    rfr::util::merge_two_vectors(features, &set_features[i*num_features], fv, num_features);
+
+                                    num_t m , v; index_t n;
+                                    std::tie(m, v, n) = t.predict_mean_var_N(fv);
+
+                                    tree_mean_stats(m); tree_var_stats(v); sample_mean_stats(m); sample_var_stats(v);
+                            }
+
+                            mean_stats(tree_mean_stats.mean());
+                            var_stats(std::max<num_t>(0, tree_var_stats.mean()));
+
+                    }
+
+                    return(std::make_tuple(mean_stats.mean(), std::max<num_t>(0, mean_stats.variance()) + std::max<num_t>(0, var_stats.mean()/set_size), std::max<num_t>(0,sample_mean_stats.variance() + sample_var_stats.mean())));
+            }
+            */
+
+            /* \brief estimates the covariance of two feature vectors
+             *
+             *
+             * The covariance between to input vectors contains information about the
+             * feature space. For other models, like GPs, this is a natural quantity
+             * (e.g., property of the kernel). Here, we try to estimate it using the
+             * emprical covariance of the individual tree's predictions.
+             *
+             * \param f1 a valid feature vector (no sanity checks are performed!)
+             * \param f2 a second feature vector (no sanity checks are performed!)
+             */
+
+            num_t covariance(const std::vector<num_t> &f1, const std::vector<num_t> &f2) {
+
+              rfr::util::running_covariance<num_t> run_cov_of_means;
+
+              for (auto &t: the_trees)
+                run_cov_of_means.push(t.predict(f1), t.predict(f2));
+
+              return (run_cov_of_means.covariance());
+            }
+
+
+            /* \brief computes the kernel of a 'Kernel Random Forest'
+             *
+             * Source: "Random forests and kernel methods" by Erwan Scornet
+             *
+             * The covariance between to input vectors contains information about the
+             * feature space. For other models, like GPs, this is a natural quantity
+             * (e.g., property of the kernel). Here, we try to estimate it using the
+             * emprical covariance of the individual tree's predictions.
+             *
+             * \param f1 a valid feature vector (no sanity checks are performed!)
+             * \param f2 a second feature vector (no sanity checks are performed!)
+             */
+
+            num_t kernel(const std::vector<num_t> &f1, const std::vector<num_t> &f2) {
+
+              rfr::util::running_statistics<num_t> stat;
+
+              for (auto &t: the_trees) {
+                auto l1 = t.find_leaf_index(f1);
+                auto l2 = t.find_leaf_index(f2);
+
+                stat.push(l1 == l2);
+              }
+              return (stat.mean());
+            }
+
+
+            std::vector<tree_type> get_all_trees() const { return the_trees; }
+
+
+            /* \brief Collects the response values for each tree in the forest for one feature vector
+             *
+             * \param feature_vector vector with the feature values
+             * \return A nested vector with the response values for each tree.
+             */
+            std::vector<std::vector<num_t> > all_leaf_values(const std::vector<num_t> &feature_vector) const {
+              std::vector<std::vector<num_t> > rv;
+              rv.reserve(the_trees.size());
+
+              for (auto &t: the_trees) {
+                rv.push_back(t.leaf_entries(feature_vector));
+              }
+              return (rv);
+            }
+
+            /* \brief Collects the response values of the trees for a set of instance feature vectors and then
+             *        for each tree takes the average over all the collected response values.
+             *
+             *        In the case of log transformation the response values are decompressed before averaging.
+             *
+             * \param feature_matrix: nested vector where each internal vector is a feature vector
+             * \param log_y: boolean which determines if there should be corrected for log transforms
+             *
+             * \return The predicted cost marginalized over all the input vectors for each tree in the forest.
+             */
+            std::vector<num_t>
+            predict_marginalized_over_instances(const std::vector<std::vector<num_t>> &feature_matrix,
+                                                const bool log_y = false) const {
+              std::vector<num_t> tree_values(the_trees.size(), 0.0);
+              int tree_id, entry_counter;
+              num_t value;
+              for (tree_id = the_trees.size() - 1; tree_id >= 0; tree_id--) {
+                entry_counter = 1;
+                for (auto &feature_vector: feature_matrix) {
+                  for (auto val: the_trees[tree_id].leaf_entries(feature_vector)) {
+                    value = log_y ? std::exp(val) : val;
+                    tree_values[tree_id] += (value - tree_values[tree_id]) / entry_counter;
+                    entry_counter++;
+                  }
+                }
+
+                if (log_y) {
+                  tree_values[tree_id] = std::log(tree_values[tree_id]);
+                }
+              }
+
+              return (tree_values);
+            }
+
+            /* \brief Collects the predictions for each tree in the forest for multiple configurations over a
+             *        set of instances.
+             *
+             *        Each configuration vector is combined with all the instance feature vectors. Based on the
+             *        response values over all these feature vectors the mean is computed.
+             *
+             *        In the case of log transformation the response values are decompressed before averaging.
+             *
+             * \param configuration_matrix: nested vector where each internal vector is a configuration vector
+             * \param feature_matrix: nested vector where each internal vector is an instance feature vector
+             * \param log_y: boolean which determines if there should be corrected for log transforms
+             *
+             * \return For each configuration the predicted cost marginalized over the instances of each tree in the forest
+             */
+            std::vector<std::vector<num_t>>
+            predict_marginalized_over_instances_batch(const std::vector<std::vector<num_t>> configuration_matrix,
+                                                      const std::vector<std::vector<num_t>> feature_matrix,
+                                                      const bool log_y = false) const {
+              int configuration_length = configuration_matrix[0].size();
+              std::vector<num_t> features(configuration_length + feature_matrix[0].size());
+              std::vector<std::vector<num_t>> predictions(configuration_matrix.size(),
+                                                          std::vector<num_t>(the_trees.size(), 0.0));
+
+              int config_id, tree_id, entry_counter;
+              num_t value;
+
+              for (config_id = configuration_matrix.size() - 1; config_id >= 0; config_id--) {
+                std::copy(configuration_matrix[config_id].begin(), configuration_matrix[config_id].end(),
+                          features.begin());
+                for (tree_id = the_trees.size() - 1; tree_id >= 0; tree_id--) {
+                  entry_counter = 1;
+                  for (auto &feature_vector: feature_matrix) {
+                    std::copy(feature_vector.begin(), feature_vector.end(), features.begin() + configuration_length);
+                    for (auto val: the_trees[tree_id].leaf_entries(features)) {
+                      value = log_y ? std::exp(val) : val;
+                      predictions[config_id][tree_id] += (value - predictions[config_id][tree_id]) / entry_counter;
+                      entry_counter++;
+                    }
+                  }
+
+                  if (log_y) {
+                    predictions[config_id][tree_id] = std::log(predictions[config_id][tree_id]);
+                  }
+                }
+              }
+
+              return (predictions);
+            }
+
+
+
+            /* \brief returns the predictions of every tree marginalized over the NAN values in the feature_vector
+             *
+             * TODO: more documentation over how the 'missing values' are handled
+             *
+             * \param feature_vector non-specfied values (NaN) will be marginalized over according to the training data
+             */
+            //std::vector<num_t> marginalized_mean_predictions(const std::vector<num_t> &feature_vector) const {
+            //	std::vector<num_t> rv;
+            //	rv.reserve(the_trees.size());
+            //	for (auto &t : the_trees)
+            //		rv.emplace_back(t.marginalized_mean_prediction(feature_vector));
+            //	return(rv);
+            //}
+
+
+
+            /* \brief updates the forest by adding the provided datapoint without a complete retraining
+             *
+             *
+             * As retraining can be quite expensive, this function can be used to quickly update the forest
+             * by finding the leafs the datapoints belong into and just inserting them. This is, of course,
+             * not the right way to do it for many data points, but it should be a good approximation for a few.
+             *
+             * \param features a valid feature vector
+             * \param response the corresponding response value
+             * \param weight the associated weight
+             */
+            void pseudo_update(std::vector<num_t> features, response_t response, num_t weight) {
+              for (auto &t: the_trees)
+                t.pseudo_update(features, response, weight);
+            }
+
+            /* \brief undoing a pseudo update by removing a point
+             *
+             * This function removes one point from the corresponding leaves into
+             * which the given feature vector falls
+             *
+             * \param features a valid feature vector
+             * \param response the corresponding response value
+             * \param weight the associated weight
+             */
+            void pseudo_downdate(std::vector<num_t> features, response_t response, num_t weight) {
+              for (auto &t: the_trees)
+                t.pseudo_downdate(features, response, weight);
+            }
+
+            num_t out_of_bag_error() { return (oob_error); }
+
+            /* \brief writes serialized representation into a binary file
+             *
+             * \param filename name of the file to store the forest in. Make sure that the directory exists!
+             */
+            void save_to_binary_file(const std::string filename) {
+              std::ofstream ofs(filename, std::ios::binary);
+              binary_oarch_t oarch(ofs);
+              serialize(oarch);
+            }
+
+            /* \brief deserialize from a binary file created by save_to_binary_file
+             *
+             * \param filename name of the file in which the forest is stored.
+             */
+            void load_from_binary_file(const std::string filename) {
+              std::ifstream ifs(filename, std::ios::binary);
+              binary_iarch_t iarch(ifs);
+              serialize(iarch);
+            }
+
+            /* serialize into a string; used for Python's pickle.dump
+             *
+             * \return std::string a JSON serialization of the forest
+             */
+            std::string ascii_string_representation() {
+              std::stringstream oss;
+              {
+                ascii_oarch_t oarch(oss);
+                serialize(oarch);
+              }
+              return (oss.str());
+            }
+
+            /* \brief deserialize from string; used for Python's pickle.load
+             *
+             * \return std::string a JSON serialization of the forest
+             */
+            void load_from_ascii_string(std::string const &str) {
+              std::stringstream iss;
+              iss.str(str);
+              ascii_iarch_t iarch(iss);
+              serialize(iarch);
+            }
+
+
+            /* \brief stores a latex document for every individual tree
+             *
+             * \param filename_template a string to specify the location and the naming scheme. Note the directory is not created, so make sure it exists.
+             *
+             */
+            void save_latex_representation(const std::string filename_template) {
+              for (auto i = 0u; i < the_trees.size(); i++) {
+                std::stringstream filename;
+                filename << filename_template << i << ".tex";
+                the_trees[i].save_latex_representation(filename.str().c_str());
+              }
+            }
+
+            void print_info() {
+              for (auto t: the_trees) {
+                t.print_info();
+              }
+            }
+
+
+            virtual unsigned int num_trees() { return (the_trees.size()); }
+
+        };
+
+
+    }
+}//namespace rfr::forests
 #endif


Test patch:
diff --git a/tests/pyrfr_unit_test_binary_regression_forest.py b/tests/pyrfr_unit_test_binary_regression_forest.py
index b53a1c4..4ed7d89 100644
--- a/tests/pyrfr_unit_test_binary_regression_forest.py
+++ b/tests/pyrfr_unit_test_binary_regression_forest.py
@@ -9,6 +9,9 @@
 
 import pyrfr.regression as reg
 
+# Helper function
+def flatten(lists):
+	return [i for l in lists for i in l]
 
 class TestBinaryRssRegressionForest(unittest.TestCase):
 
@@ -88,6 +91,43 @@ def test_contents(self):
 		self.assertAlmostEqual(trees[0].get_node(0).get_num_split_value(), -0.20652545999380345, 2)
 		self.assertAlmostEqual(trees[1].get_node(1).get_num_split_value(), -0.8857842741236657, 2)
 
+	def test_marginalized_over_instances(self):
+		self.forest.options.num_trees = 10
+		self.forest.fit(self.data, self.rng)
+
+		data = [self.data.retrieve_data_point(i) for i in range(4)]
+		leaf_entries = [self.forest.all_leaf_values(x) for x in data]
+		leaf_entries = list(zip(*leaf_entries))  # Entries per tree
+		leaf_entries = [flatten(e) for e in leaf_entries]
+		gt_means = [sum(entries) / len(entries) for entries in leaf_entries]  # Ground truth
+
+		mean_per_tree = self.forest.predict_marginalized_over_instances(data, False)
+		self.assertEqual(len(mean_per_tree), 10)
+		for mean_test, mean_gt in zip(mean_per_tree, gt_means):
+			self.assertAlmostEqual(mean_test, mean_gt, 6)
+
+	def test_marginalized_over_instances_batch(self):
+		self.forest.options.num_trees = 10
+		self.forest.fit(self.data, self.rng)
+
+		split = len(self.data.retrieve_data_point(0)) // 2
+		data_conf = [self.data.retrieve_data_point(i)[:split] for i in range(4)]
+		data_instances = [self.data.retrieve_data_point(i)[split:] for i in range(4)]
+
+		config_means = self.forest.predict_marginalized_over_instances_batch(data_conf, data_instances, False)
+		self.assertEqual(len(config_means), 4)
+		for conf_id, mean_per_tree in enumerate(config_means):
+			self.assertEqual(len(mean_per_tree), 10)
+
+			data = [data_conf[conf_id] + inst for inst in data_instances]
+			leaf_entries = [self.forest.all_leaf_values(x) for x in data]
+			leaf_entries = list(zip(*leaf_entries))  # Entries per tree
+			leaf_entries = [flatten(e) for e in leaf_entries]
+			gt_means = [sum(entries) / len(entries) for entries in leaf_entries]  # Ground truth
+
+			for mean_test, mean_gt in zip(mean_per_tree, gt_means):
+				self.assertAlmostEqual(mean_test, mean_gt, 6)
+
 	def test_pickling(self):
 		
 		the_forest = reg.binary_rss_forest()
diff --git a/tests/unit_test_regression_forest.cpp b/tests/unit_test_regression_forest.cpp
index 5bc7b1e..b2bcc58 100644
--- a/tests/unit_test_regression_forest.cpp
+++ b/tests/unit_test_regression_forest.cpp
@@ -16,11 +16,11 @@
 #include <cereal/archives/xml.hpp>
 #include <cereal/archives/portable_binary.hpp>
 #include <cereal/archives/binary.hpp>
+
 typedef cereal::PortableBinaryInputArchive iarch_type;
 typedef cereal::PortableBinaryOutputArchive oarch_type;
 
 
-
 typedef double num_t;
 typedef double response_t;
 typedef unsigned int index_t;
@@ -36,253 +36,306 @@ typedef rfr::nodes::temporary_node<num_t, index_t> tmp_node_type;
 typedef rfr::trees::k_ary_random_tree<2, node_type, num_t, response_t, index_t, rng_t> tree_type;
 
 
-typedef rfr::forests::regression_forest< tree_type, num_t, response_t, index_t, rng_t> forest_type;
+typedef rfr::forests::regression_forest<tree_type, num_t, response_t, index_t, rng_t> forest_type;
 
-typedef rfr::forests::quantile_regression_forest< tree_type, num_t, response_t, index_t, rng_t> qrf_type;
+typedef rfr::forests::quantile_regression_forest<tree_type, num_t, response_t, index_t, rng_t> qrf_type;
 
 
 typedef rfr::forests::fANOVA_forest<split_type, num_t, response_t, index_t, rng_t> fANOVAf_type;
 
 
+data_container_type load_diabetes_data() {
+  data_container_type data(10);
 
-data_container_type load_diabetes_data(){
-	data_container_type data(10);
-	
-    std::string feature_file, response_file;
-    
-    feature_file  = std::string(boost::unit_test::framework::master_test_suite().argv[1]) + "diabetes_features.csv";
-    response_file = std::string(boost::unit_test::framework::master_test_suite().argv[1]) + "diabetes_responses.csv";
+  std::string feature_file, response_file;
 
-    data.import_csv_files(feature_file, response_file);
-    return(data);
+  feature_file = std::string(boost::unit_test::framework::master_test_suite().argv[1]) + "diabetes_features.csv";
+  response_file = std::string(boost::unit_test::framework::master_test_suite().argv[1]) + "diabetes_responses.csv";
+
+  data.import_csv_files(feature_file, response_file);
+  return (data);
 }
 
-BOOST_AUTO_TEST_CASE( regression_forest_serialize_test ){
-    
-    
-    auto data = load_diabetes_data();
+BOOST_AUTO_TEST_CASE(regression_forest_serialize_test) {
 
-	rfr::trees::tree_options<num_t, response_t, index_t> tree_opts;
-	tree_opts.min_samples_to_split = 2;
-	tree_opts.min_samples_in_leaf = 1;
-	tree_opts.max_features = data.num_features()*3/4;
-	tree_opts.max_depth=5;
-	tree_opts.max_num_nodes=2048;
 
-	
-	rfr::forests::forest_options<num_t, response_t, index_t> forest_opts(tree_opts);
+  auto data = load_diabetes_data();
 
-	forest_opts.num_data_points_per_tree = data.num_data_points()/2;
-	forest_opts.num_trees = 16;
-	forest_opts.do_bootstrapping = true;
-	forest_opts.compute_oob_error= true;
-	
-	forest_type the_forest(forest_opts);
-	
-	rng_t rng;
+  rfr::trees::tree_options<num_t, response_t, index_t> tree_opts;
+  tree_opts.min_samples_to_split = 2;
+  tree_opts.min_samples_in_leaf = 1;
+  tree_opts.max_features = data.num_features() * 3 / 4;
+  tree_opts.max_depth = 5;
+  tree_opts.max_num_nodes = 2048;
 
-	BOOST_REQUIRE(std::isnan(the_forest.out_of_bag_error()));
 
-	the_forest.fit(data, rng);
+  rfr::forests::forest_options<num_t, response_t, index_t> forest_opts(tree_opts);
 
-	BOOST_REQUIRE(!std::isnan(the_forest.out_of_bag_error()));
+  forest_opts.num_data_points_per_tree = data.num_data_points() / 2;
+  forest_opts.num_trees = 16;
+  forest_opts.do_bootstrapping = true;
+  forest_opts.compute_oob_error = true;
 
-	BOOST_REQUIRE_EQUAL(the_forest.get_all_trees().size(), 16);
+  forest_type the_forest(forest_opts);
 
-    auto tmp = the_forest.predict(data.retrieve_data_point(5));
+  rng_t rng;
 
-	std::ostringstream oss;
-	
-	{
-		oarch_type oarchive(oss);
-		oarchive(the_forest);
-	}
+  BOOST_REQUIRE(std::isnan(the_forest.out_of_bag_error()));
 
-	
-	the_forest.save_to_binary_file("regression_forest_test.bin");
-	
-	forest_type the_forest2;
+  the_forest.fit(data, rng);
 
-	auto s = the_forest.ascii_string_representation();
+  BOOST_REQUIRE(!std::isnan(the_forest.out_of_bag_error()));
 
-	the_forest2.load_from_ascii_string(s);
+  BOOST_REQUIRE_EQUAL(the_forest.get_all_trees().size(), 16);
 
-	
-	forest_type the_forest3;
-	the_forest3.load_from_binary_file("regression_forest_test.bin");
+  auto tmp = the_forest.predict(data.retrieve_data_point(5));
 
+  std::ostringstream oss;
 
-	for (auto i=0u; i < data.num_data_points(); ++i){
-		auto v1 = the_forest.predict(data.retrieve_data_point(i));
-		auto v2 = the_forest2.predict(data.retrieve_data_point(i));
-		auto v3 = the_forest3.predict(data.retrieve_data_point(i));
+  {
+    oarch_type oarchive(oss);
+    oarchive(the_forest);
+  }
 
-		BOOST_REQUIRE_EQUAL(v1,v2);
-		BOOST_REQUIRE_EQUAL(v1,v3);
-	}
 
+  the_forest.save_to_binary_file("regression_forest_test.bin");
+
+  forest_type the_forest2;
+
+  auto s = the_forest.ascii_string_representation();
+
+  the_forest2.load_from_ascii_string(s);
+
+
+  forest_type the_forest3;
+  the_forest3.load_from_binary_file("regression_forest_test.bin");
+
+
+  for (auto i = 0u; i < data.num_data_points(); ++i) {
+    auto v1 = the_forest.predict(data.retrieve_data_point(i));
+    auto v2 = the_forest2.predict(data.retrieve_data_point(i));
+    auto v3 = the_forest3.predict(data.retrieve_data_point(i));
+
+    BOOST_REQUIRE_EQUAL(v1, v2);
+    BOOST_REQUIRE_EQUAL(v1, v3);
+  }
 
-	
 
 }
 
 
+BOOST_AUTO_TEST_CASE(regression_forest_update_downdate_tests) {
 
+  double unique_value = 42.424242;
 
-BOOST_AUTO_TEST_CASE( regression_forest_update_downdate_tests ){
-	
-	double unique_value = 42.424242;
-	
-	auto data = load_diabetes_data();
+  auto data = load_diabetes_data();
 
-	rfr::trees::tree_options<num_t, response_t, index_t> tree_opts;
-	tree_opts.min_samples_to_split = 2;
-	tree_opts.min_samples_in_leaf = 1;
-	tree_opts.max_features = 10;
+  rfr::trees::tree_options<num_t, response_t, index_t> tree_opts;
+  tree_opts.min_samples_to_split = 2;
+  tree_opts.min_samples_in_leaf = 1;
+  tree_opts.max_features = 10;
 
-	rfr::forests::forest_options<num_t, response_t, index_t> forest_opts(tree_opts);
+  rfr::forests::forest_options<num_t, response_t, index_t> forest_opts(tree_opts);
 
-	forest_opts.num_data_points_per_tree = data.num_data_points();
-	forest_opts.num_trees = 10;
-	forest_opts.do_bootstrapping = true;
+  forest_opts.num_data_points_per_tree = data.num_data_points();
+  forest_opts.num_trees = 10;
+  forest_opts.do_bootstrapping = true;
 
 
-	forest_type the_forest(forest_opts);
-	
-	rng_t rng;
+  forest_type the_forest(forest_opts);
 
-	//fit forest
-	the_forest.fit(data, rng);
+  rng_t rng;
 
-	// get reference leaf values for one configuration
-	std::vector<std::vector< num_t> > before = the_forest.all_leaf_values(data.retrieve_data_point(0));
+  //fit forest
+  the_forest.fit(data, rng);
 
-	// update forest with that configuration and a unique response value
-	the_forest.pseudo_update(data.retrieve_data_point(0), unique_value, 1);
-	
-	// get new leaf values
-	std::vector<std::vector< num_t> > after_update = the_forest.all_leaf_values(data.retrieve_data_point(0));
-	
-	// compare them to ensure the data point has been added correctly
-	
-	for (auto i=0u; i < before.size(); ++i){
-		BOOST_REQUIRE(before[i].size() == after_update[i].size() -1);
-		BOOST_CHECK_EQUAL_COLLECTIONS( before[i].begin(), before[i].end(), after_update[i].begin(), --after_update[i].end());
-		BOOST_REQUIRE( after_update[i].back() == unique_value);
-	}
-	
-	// downdate the tree
-	the_forest.pseudo_downdate(data.retrieve_data_point(0), unique_value, 1);
-	
-	// get new leaf values
-	std::vector<std::vector< num_t> > after_downdate = the_forest.all_leaf_values(data.retrieve_data_point(0));
-	
-	// compare them to ensure the last data point has been removed
-	for (auto i=0u; i < before.size(); ++i)
-		BOOST_CHECK_EQUAL_COLLECTIONS( before[i].begin(), before[i].end(), after_downdate[i].begin(), after_downdate[i].end());
-	
-	auto m = the_forest.predict(data.retrieve_data_point(0));
-}
+  // get reference leaf values for one configuration
+  std::vector<std::vector<num_t> > before = the_forest.all_leaf_values(data.retrieve_data_point(0));
 
+  // update forest with that configuration and a unique response value
+  the_forest.pseudo_update(data.retrieve_data_point(0), unique_value, 1);
 
-BOOST_AUTO_TEST_CASE( regression_forest_exceptions_tests ){
-    
-    auto data = load_diabetes_data();
+  // get new leaf values
+  std::vector<std::vector<num_t> > after_update = the_forest.all_leaf_values(data.retrieve_data_point(0));
 
-	rfr::trees::tree_options<num_t, response_t, index_t> tree_opts;
-	tree_opts.min_samples_to_split = 2;
-	tree_opts.min_samples_in_leaf = 1;
-	tree_opts.max_features = data.num_data_points()*3/4;
+  // compare them to ensure the data point has been added correctly
 
-	rfr::forests::forest_options<num_t, response_t, index_t> forest_opts(tree_opts);
+  for (auto i = 0u; i < before.size(); ++i) {
+    BOOST_REQUIRE(before[i].size() == after_update[i].size() - 1);
+    BOOST_CHECK_EQUAL_COLLECTIONS(before[i].begin(), before[i].end(), after_update[i].begin(), --after_update[i].end());
+    BOOST_REQUIRE(after_update[i].back() == unique_value);
+  }
 
-	forest_opts.num_data_points_per_tree = data.num_data_points();
-	forest_opts.num_trees = 8;
-	forest_opts.do_bootstrapping = false;
-	forest_opts.compute_oob_error= true;
-	
-	forest_type the_forest(forest_opts);
-	
-	rng_t rng;
+  // downdate the tree
+  the_forest.pseudo_downdate(data.retrieve_data_point(0), unique_value, 1);
 
-	// no trees in the forest
-	the_forest.options.num_trees = 0;
-	BOOST_REQUIRE_THROW(the_forest.fit(data, rng), std::runtime_error);
+  // get new leaf values
+  std::vector<std::vector<num_t> > after_downdate = the_forest.all_leaf_values(data.retrieve_data_point(0));
 
-	// no datapoints in any trees
-	the_forest.options.num_trees = 8;
-	the_forest.options.num_data_points_per_tree = 0;
-	BOOST_REQUIRE_THROW(the_forest.fit(data, rng), std::runtime_error);
+  // compare them to ensure the last data point has been removed
+  for (auto i = 0u; i < before.size(); ++i)
+    BOOST_CHECK_EQUAL_COLLECTIONS(before[i].begin(), before[i].end(), after_downdate[i].begin(),
+                                  after_downdate[i].end());
 
-	// not enough datapoints and no bootstrapping
-	the_forest.options.num_data_points_per_tree = 32768;
-	BOOST_REQUIRE_THROW(the_forest.fit(data, rng), std::runtime_error);
+  auto m = the_forest.predict(data.retrieve_data_point(0));
 
-	// a successful training without bootstrapping for covarage
-	the_forest.options.num_data_points_per_tree = 128;
-	the_forest.fit(data, rng);
+}
 
-	// no features to split
-	the_forest.options.num_data_points_per_tree = 128;
-	the_forest.options.tree_opts.max_features = 0;
-	BOOST_REQUIRE_THROW(the_forest.fit(data, rng), std::runtime_error);
+BOOST_AUTO_TEST_CASE(regression_forest_predict_marginalized) {
+  auto data = load_diabetes_data();
+
+  rfr::trees::tree_options<num_t, response_t, index_t> tree_opts;
+  tree_opts.min_samples_to_split = 2;
+  tree_opts.min_samples_in_leaf = 1;
+  tree_opts.max_features = 10;
+
+  rfr::forests::forest_options<num_t, response_t, index_t> forest_opts(tree_opts);
+
+  forest_opts.num_data_points_per_tree = data.num_data_points();
+  forest_opts.num_trees = 10;
+  forest_opts.do_bootstrapping = true;
+
+
+  forest_type the_forest(forest_opts);
+
+  rng_t rng;
+
+  // fit forest
+  the_forest.fit(data, rng);
+
+  // compute marginal cost over set of feature vectors
+  const int n_samples = 4;
+  std::vector<std::vector<num_t>> feature_matrix(n_samples);
+  for (int i = 0; i < n_samples; i++) {
+    feature_matrix.push_back(data.retrieve_data_point(i));
+  }
+
+  std::vector<num_t> mean_per_tree = the_forest.predict_marginalized_over_instances(feature_matrix, false);
+  BOOST_CHECK_EQUAL(mean_per_tree.size(), 10);
+
+
+  // predict_marginalized_over_instances_batch
+  int vectorsize = feature_matrix[0].size();
+  int splitpoint = vectorsize / 2;
+  int i, j;
+  std::vector<std::vector<num_t>> config_matrix(n_samples, std::vector<num_t>(splitpoint));
+  std::vector<std::vector<num_t>> feature_matrix2(n_samples, std::vector<num_t>(vectorsize - splitpoint));
+  for (i = 0; i < n_samples; i++) {
+    for (j = 0; j < vectorsize; j++) {
+      if (j < splitpoint) {
+        config_matrix[i][j] = feature_matrix[i][j];
+      } else {
+        feature_matrix2[i][j - splitpoint] = feature_matrix[i][j];
+      }
+    }
+  }
+
+  std::vector<std::vector<num_t>> config_means = the_forest.predict_marginalized_over_instances_batch(config_matrix,
+                                                                                                      feature_matrix2,
+                                                                                                      false);
+  BOOST_CHECK_EQUAL(config_means.size(), n_samples);
+  for (auto &mean_per_tree: config_means) {
+    BOOST_CHECK_EQUAL(mean_per_tree.size(), 10);
+  }
 }
 
 
+BOOST_AUTO_TEST_CASE(regression_forest_exceptions_tests) {
 
+  auto data = load_diabetes_data();
 
-BOOST_AUTO_TEST_CASE( quantile_regression_forest_test ){
-	
-	auto data = load_diabetes_data();
+  rfr::trees::tree_options<num_t, response_t, index_t> tree_opts;
+  tree_opts.min_samples_to_split = 2;
+  tree_opts.min_samples_in_leaf = 1;
+  tree_opts.max_features = data.num_data_points() * 3 / 4;
 
-	rng_t rng;
+  rfr::forests::forest_options<num_t, response_t, index_t> forest_opts(tree_opts);
 
-	rfr::trees::tree_options<num_t, response_t, index_t> tree_opts;
-	tree_opts.min_samples_to_split = 2;
-	tree_opts.min_samples_in_leaf = 1;
-	tree_opts.max_features = 10;
+  forest_opts.num_data_points_per_tree = data.num_data_points();
+  forest_opts.num_trees = 8;
+  forest_opts.do_bootstrapping = false;
+  forest_opts.compute_oob_error = true;
 
-	// don't split anything
+  forest_type the_forest(forest_opts);
 
-	tree_opts.max_num_nodes = 1;
-	tree_opts.max_depth = 0;
+  rng_t rng;
 
-	rfr::forests::forest_options<num_t, response_t, index_t> forest_opts(tree_opts);
+  // no trees in the forest
+  the_forest.options.num_trees = 0;
+  BOOST_REQUIRE_THROW(the_forest.fit(data, rng), std::runtime_error);
 
-	forest_opts.num_data_points_per_tree = data.num_data_points();
-	forest_opts.num_trees = 1;
-	forest_opts.do_bootstrapping = false;
+  // no datapoints in any trees
+  the_forest.options.num_trees = 8;
+  the_forest.options.num_data_points_per_tree = 0;
+  BOOST_REQUIRE_THROW(the_forest.fit(data, rng), std::runtime_error);
 
+  // not enough datapoints and no bootstrapping
+  the_forest.options.num_data_points_per_tree = 32768;
+  BOOST_REQUIRE_THROW(the_forest.fit(data, rng), std::runtime_error);
 
-	// just to test the default constructor
-	qrf_type sudowoodo;
-	BOOST_REQUIRE_THROW(sudowoodo.fit(data, rng), std::runtime_error);
+  // a successful training without bootstrapping for covarage
+  the_forest.options.num_data_points_per_tree = 128;
+  the_forest.fit(data, rng);
 
-	qrf_type the_forest(forest_opts);
+  // no features to split
+  the_forest.options.num_data_points_per_tree = 128;
+  the_forest.options.tree_opts.max_features = 0;
+  BOOST_REQUIRE_THROW(the_forest.fit(data, rng), std::runtime_error);
+}
 
-	//fit forest
-	the_forest.fit(data, rng);
 
-	auto mew = data.retrieve_data_point(0);
-	std::vector<num_t> mew2 = {0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1};
+BOOST_AUTO_TEST_CASE(quantile_regression_forest_test) {
+
+  auto data = load_diabetes_data();
+
+  rng_t rng;
+
+  rfr::trees::tree_options<num_t, response_t, index_t> tree_opts;
+  tree_opts.min_samples_to_split = 2;
+  tree_opts.min_samples_in_leaf = 1;
+  tree_opts.max_features = 10;
+
+  // don't split anything
+
+  tree_opts.max_num_nodes = 1;
+  tree_opts.max_depth = 0;
+
+  rfr::forests::forest_options<num_t, response_t, index_t> forest_opts(tree_opts);
+
+  forest_opts.num_data_points_per_tree = data.num_data_points();
+  forest_opts.num_trees = 1;
+  forest_opts.do_bootstrapping = false;
+
+
+  // just to test the default constructor
+  qrf_type sudowoodo;
+  BOOST_REQUIRE_THROW(sudowoodo.fit(data, rng), std::runtime_error);
+
+  qrf_type the_forest(forest_opts);
+
+  //fit forest
+  the_forest.fit(data, rng);
+
+  auto mew = data.retrieve_data_point(0);
+  std::vector<num_t> mew2 = {0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1};
 
-	auto qv1 = the_forest.predict_quantiles(mew , mew2);
+  auto qv1 = the_forest.predict_quantiles(mew, mew2);
 
-	BOOST_REQUIRE_EQUAL(qv1.size(), mew2.size());
+  BOOST_REQUIRE_EQUAL(qv1.size(), mew2.size());
 
-	// check that shuffling doesn't affect the outcome
-	std::shuffle(mew2.begin(), mew2.end(), rng);
-	auto qv2 = the_forest.predict_quantiles(mew , mew2);
-	BOOST_CHECK_EQUAL_COLLECTIONS ( qv1.begin(), qv1.end(), qv2.begin(), qv2.end());
+  // check that shuffling doesn't affect the outcome
+  std::shuffle(mew2.begin(), mew2.end(), rng);
+  auto qv2 = the_forest.predict_quantiles(mew, mew2);
+  BOOST_CHECK_EQUAL_COLLECTIONS (qv1.begin(), qv1.end(), qv2.begin(), qv2.end());
 
-	//
-	std::vector<num_t> qv_numpy_ref = {  25.,   60.,   77.,   94.,  115.,  141.,  168.,  197.,  233., 268.,  346.};
+  //
+  std::vector<num_t> qv_numpy_ref = {25., 60., 77., 94., 115., 141., 168., 197., 233., 268., 346.};
 
-	BOOST_CHECK_EQUAL_COLLECTIONS ( qv1.begin(), qv1.end(), qv_numpy_ref.begin(), qv_numpy_ref.end());
+  BOOST_CHECK_EQUAL_COLLECTIONS (qv1.begin(), qv1.end(), qv_numpy_ref.begin(), qv_numpy_ref.end());
 
 
-	BOOST_REQUIRE_THROW(the_forest.predict_quantiles(mew, {0,-0.5,1}) ,std::runtime_error);
-	BOOST_REQUIRE_THROW(the_forest.predict_quantiles(mew, {1.1,0.5}) ,std::runtime_error);
+  BOOST_REQUIRE_THROW(the_forest.predict_quantiles(mew, {0, -0.5, 1}), std::runtime_error);
+  BOOST_REQUIRE_THROW(the_forest.predict_quantiles(mew, {1.1, 0.5}), std::runtime_error);
 }
 
 /* not interesting right now!