text:
/* Copyright 2003-2013 Joaquin M Lopez Munoz.
 * Distributed under the Boost Software License, Version 1.0.
 * (See accompanying file LICENSE_1_0.txt or copy at
 * http://www.boost.org/LICENSE_1_0.txt)
 *
 * See http://www.boost.org/libs/multi_index for library home page.
 */

#ifndef BOOST_MULTI_INDEX_DETAIL_NODE_TYPE_HPP
#define BOOST_MULTI_INDEX_DETAIL_NODE_TYPE_HPP

#if defined(_MSC_VER)
#pragma once
#endif

#include <boost/config.hpp> /* keep it first to prevent nasty warns in MSVC */
#include <boost/detail/workaround.hpp>
#include <boost/mpl/bind.hpp>
#include <boost/mpl/reverse_iter_fold.hpp>
#include <boost/mpl/deref.hpp>
#include <boost/multi_index_container_fwd.hpp>
#include <boost/multi_index/detail/header_holder.hpp>
#include <boost/multi_index/detail/index_node_base.hpp>
#include <boost/multi_index/detail/is_index_list.hpp>
#include <boost/static_assert.hpp>

namespace lslboost{

namespace multi_index{

namespace detail{

/* MPL machinery to construct the internal node type associated to an
 * index list.
 */

struct index_node_applier
{
  template<typename IndexSpecifierIterator,typename Super>
  struct apply
  {
    typedef typename mpl::deref<IndexSpecifierIterator>::type index_specifier;
    typedef typename index_specifier::
      BOOST_NESTED_TEMPLATE node_class<Super>::type type;
  };
};

template<typename Value,typename IndexSpecifierList,typename Allocator>
struct multi_index_node_type
{
  BOOST_STATIC_ASSERT(detail::is_index_list<IndexSpecifierList>::value);

  typedef typename mpl::reverse_iter_fold<
    IndexSpecifierList,
    index_node_base<Value,Allocator>,
    mpl::bind2<index_node_applier,mpl::_2,mpl::_1>
  >::type type;
};

} /* namespace multi_index::detail */

} /* namespace multi_index */

} /* namespace lslboost */

#endif

ext: hpp

max_stars_repo_name: samuelpowell/liblsl

max_issues_repo_name: samuelpowell/liblsl

max_forks_repo_name: samuelpowell/liblsl


text:
//
//   Copyright 2014 QuarksLab
//
//   Licensed under the Apache License, Version 2.0 (the "License");
//   you may not use this file except in compliance with the License.
//   You may obtain a copy of the License at
//
//       http://www.apache.org/licenses/LICENSE-2.0
//
//   Unless required by applicable law or agreed to in writing, software
//   distributed under the License is distributed on an "AS IS" BASIS,
//   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
//   See the License for the specific language governing permissions and
//   limitations under the License.
//

#ifndef PE_DECODER_HPP_INCLUDED
#define PE_DECODER_HPP_INCLUDED

#include "binmap/collectors/pe.hpp"
#include "binmap/collectors/pe/data.hpp"
#include "binmap/metadata.hpp"

#include <boost/filesystem/path.hpp>

#include <vector>
#include <map>

/* \brief decode the PE format
 * uses a template helper class to distinguish 32/64 bits versions, aka PE and
 * PE+
*/

struct PEDecoder {
  virtual ~PEDecoder() {}

  // Tell if the module is a PE module.
  virtual bool is_compatible(void) const = 0;

  // get sections from PE
  virtual const std::vector<PeSectionHeader> &sections(void) const = 0;

  // get all imported modules for a given module
  virtual void
  get_imports(boost::filesystem::path const &module_path,
              std::set<boost::filesystem::path> &imports) const = 0;

  virtual bool get_exports(boost::filesystem::path const &module_path,
                           std::vector<std::string> &exports) const = 0;
  virtual bool get_delay_imports(boost::filesystem::path const &module_path,
				std::vector<std::string> &imported_symbols) const =0;

  virtual bool get_imported_symbols(boost::filesystem::path const &module_path,
                           	    std::vector<std::string> &imported_symbols) const = 0;

  virtual void extract_hardening_features(MetadataInfo&) const = 0;

  static uint16_t machine_type(std::ifstream &file);
};

/* Used to create PEDecoder instances from a given PE file path.*/
PEDecoder *PeDecoderFactory(boost::filesystem::path const &path,
                            bool full_parsing = false);
/* Used to create PEDecoder instances from a given PE file path */
PEDecoder *PeDecoderFactory(std::ifstream &file, bool full_parsing = false);

template <typename _Bits> class pe_decoder : public PEDecoder {
public:
  typedef PeNtHeadersTraits<_Bits> PeNtHeaders;
  typedef PeImageLoadConfigDirectoryTraits<_Bits> PeImageLoadConfigDirectory;

  pe_decoder(std::ifstream &file, bool full_parsing = false);

  bool is_compatible(void) const;
  uint16_t machine_type(void) const;
  bool convert_rva_to_offset(uint32_t rva, uint32_t &offset) const;
  bool convert_rva_to_offset64(uint64_t rva, uint64_t &offset) const;
  void get_imports(boost::filesystem::path const &module_path,
                   std::set<boost::filesystem::path> &imports) const;
  bool get_exports(boost::filesystem::path const &module_path,
                   std::vector<std::string> &exports) const;
  bool get_imported_symbols(boost::filesystem::path const &module_path,
                  	    std::vector<std::string> &imported_symbols) const;

  bool get_delay_imports(boost::filesystem::path const &module_path,
		   std::vector<std::string> &imported_symbols) const;
  void extract_hardening_features(MetadataInfo& mi) const;
  PeNtHeaders const nt_headers(void) const;
  std::vector<PeSectionHeader> const &sections(void) const;
  bool
  parse_manifest(std::map<std::string, boost::filesystem::path> &assembly_maps);
  char *get_section(PeDataDirectory::image_directory_entry_t entry) const;

private:
  PeData<_Bits> _pe_data;
  std::ifstream &_file;
  bool _is_compatible;
  uint16_t _machine_type;
  std::map<std::string, boost::filesystem::path> _assembly_maps;
  bool find_module_path(boost::filesystem::path const &module_path,
                        std::string &imports) const;
};

#endif


text:
#include <unistd.h>
#include <stdexcept>
#include <memory>
using namespace std;

#include <boost/scoped_ptr.hpp>
using namespace boost;

#include "calpontselectexecutionplan.h"
using namespace execplan;
#include "bytestream.h"
#include "messagequeue.h"
using namespace messageqcpp;

namespace qfe
{

MessageQueueClient* sendCSEP(CalpontSelectExecutionPlan* csep)
{
    scoped_ptr<CalpontSelectExecutionPlan> cleaner(csep);

    ByteStream bs;

    MessageQueueClient* mqc = 0;

    mqc = new MessageQueueClient("ExeMgr1");
    auto_ptr<MessageQueueClient> smqc(mqc);

    bs.reset();
    ByteStream::quadbyte wantTuples = 4;
    bs << wantTuples;
    mqc->write(bs);

    bs.reset();
    csep->serialize(bs);
    mqc->write(bs);

    SBS sbs;
    sbs = mqc->read();
    *sbs >> wantTuples;
    //cerr << "got flag: " << wantTuples << endl;
    string msg;
    sbs = mqc->read();
    *sbs >> msg;
    //cerr << "got msg: " << msg << endl;

    if (wantTuples != 0)
        throw runtime_error(msg);

    smqc.release();
    return mqc;
}

}

ext: cpp

max_stars_repo_name: zettadb/zettalib

max_issues_repo_name: zettadb/zettalib

max_forks_repo_name: zettadb/zettalib


text:
#include <string>
#include <vector>

// boost.Python
#include <boost/python.hpp>

#include "libpokercore/pokercore.h"

#ifndef POKERENGINE_H
#define POKERENGINE_H

typedef std::vector<double> CustomList;

class PokerEngine {
	public:
		PokerEngine();
		~PokerEngine();
		// 2hands & board
		CustomList EV(std::string hand1, std::string hand2, std::string board, double runs);
		std::vector<double> EV(std::vector<std::string>, std::string board, double runs);
		boost::python::list EV(boost::python::list hands, std::string board, double runs);

};

#endif

max_stars_repo_name: PhillCli/PokerEngine

max_issues_repo_name: PhillCli/PokerEngine

max_forks_repo_name: PhillCli/PokerEngine


text:
#include "stored_queries/StoredWWProbabilityQueryHandler.h"
#include <gis/Box.h>
#include <gis/OGR.h>
#include <macgyver/TimeFormatter.h>
#include <newbase/NFmiEnumConverter.h>
#include <smartmet/macgyver/Exception.h>

#include <boost/algorithm/string/replace.hpp>
#include <boost/range/adaptor/map.hpp>
#include <boost/range/algorithm/copy.hpp>
#include <iomanip>
#include <locale>

namespace
{
const char* P_FIND_NEAREST_VALID = "findNearestValid";
const char* P_LOCALE = "locale";
const char* P_MISSING_TEXT = "missingText";
const char* P_PRODUCER = "producer";
const char* P_ORIGIN_TIME = "originTime";
const char* P_CRS = "crs";
const char* P_ICAO_CODE = "icaoCode";

std::string double2string(double d, unsigned int precision)
{
  try
  {
    std::stringstream ss;

    ss << std::fixed << std::setprecision(precision) << d;

    return ss.str();
  }
  catch (...)
  {
    throw Fmi::Exception::Trace(BCP, "Operation failed!");
  }
}

#if 0
std::string bbox2string(const SmartMet::Spine::BoundingBox& bbox, OGRSpatialReference& targetSRS)
{
  try
  {
    std::stringstream ss;

    unsigned int precision(targetSRS.IsGeographic() ? 6 : 1);

    if (targetSRS.EPSGTreatsAsLatLong())
    {
      ss << std::fixed << std::setprecision(precision) << bbox.yMin << "," << bbox.xMin << " "
         << bbox.yMax << "," << bbox.xMin << " " << bbox.yMax << "," << bbox.xMax << " "
         << bbox.yMin << "," << bbox.xMax << " " << bbox.yMin << "," << bbox.xMin;
    }
    else
    {
      ss << std::fixed << std::setprecision(precision) << bbox.xMin << "," << bbox.yMin << " "
         << bbox.xMax << "," << bbox.yMin << " " << bbox.xMax << "," << bbox.yMax << " "
         << bbox.xMin << "," << bbox.yMax << " " << bbox.xMin << "," << bbox.yMin;
    }

    return ss.str();
  }
  catch (...)
  {
    throw Fmi::Exception::Trace(BCP, "Operation failed!");
  }
}
#endif

FmiParameterName get_parameter(boost::shared_ptr<SmartMet::Plugin::WFS::StoredQueryConfig> config,
                               const std::string& param_name)
{
  try
  {
    uint64_t cpid = config->get_mandatory_config_param<uint64_t>(param_name);
    return static_cast<FmiParameterName>(cpid);
  }
  catch (...)
  {
    throw Fmi::Exception::Trace(BCP, "Operation failed!");
  }
}

SmartMet::Engine::Querydata::ParameterOptions get_qengine_parameter(
    const SmartMet::Plugin::WFS::ProbabilityQueryParam& queryParam,
    const SmartMet::Spine::Parameter& smartmetParam)
{
  try
  {
    std::unique_ptr<Fmi::TimeFormatter> timeFormatter(Fmi::TimeFormatter::create("iso"));

    NFmiPoint nearestpoint(kFloatMissing, kFloatMissing);
    bool nearestFlag(false);
    std::string emptystring("");

    SmartMet::Engine::Querydata::ParameterOptions qengine_param(smartmetParam,
                                                                queryParam.producer,
                                                                *(queryParam.loc),
                                                                queryParam.loc->country,
                                                                queryParam.loc->name,
                                                                *timeFormatter,
                                                                emptystring,
                                                                queryParam.language,
                                                                queryParam.outputLocale,
                                                                queryParam.tz_name,
                                                                nearestFlag,
                                                                nearestpoint,
                                                                nearestpoint);

    return qengine_param;
  }
  catch (...)
  {
    throw Fmi::Exception::Trace(BCP, "Operation failed!");
  }
}

void populate_result_vector(
    SmartMet::Spine::TimeSeries::TimeSeriesPtr qengine_result_light,
    SmartMet::Spine::TimeSeries::TimeSeriesPtr qengine_result_moderate,
    SmartMet::Spine::TimeSeries::TimeSeriesPtr qengine_result_heavy,
    SmartMet::Plugin::WFS::WinterWeatherIntensityProbabilities& result_vector)
{
  try
  {
    unsigned int result_set_size(qengine_result_light->size());

    for (unsigned int i = 0; i < result_set_size; i++)
    {
      const SmartMet::Spine::TimeSeries::TimedValue& tv_light = qengine_result_light->at(i);
      const SmartMet::Spine::TimeSeries::TimedValue& tv_moderate = qengine_result_moderate->at(i);
      const SmartMet::Spine::TimeSeries::TimedValue& tv_heavy = qengine_result_heavy->at(i);

      if (!(boost::get<double>(&tv_light.value)))
      {
        // nan-value is string
        break;
      }

      double light_prob = *(boost::get<double>(&tv_light.value));
      double moderate_prob = *(boost::get<double>(&tv_moderate.value));
      double heavy_prob = *(boost::get<double>(&tv_heavy.value));
      result_vector.push_back(SmartMet::Plugin::WFS::WinterWeatherProbability(
          tv_light.time.utc_time(), light_prob, moderate_prob, heavy_prob));
    }
  }
  catch (...)
  {
    throw Fmi::Exception::Trace(BCP, "Operation failed!");
  }
}

}  // anonymous namespace

namespace SmartMet
{
namespace Plugin
{
namespace WFS
{
StoredWWProbabilityQueryHandler::StoredWWProbabilityQueryHandler(
    SmartMet::Spine::Reactor* reactor,
    StoredQueryConfig::Ptr config,
    PluginImpl& pluginData,
    boost::optional<std::string> templateFileName)

    : StoredQueryParamRegistry(config),
      SupportsExtraHandlerParams(config, false),
      RequiresGeoEngine(reactor),
      RequiresQEngine(reactor),
      StoredQueryHandlerBase(reactor, config, pluginData, templateFileName),
      SupportsLocationParameters(reactor, config, SUPPORT_KEYWORDS | INCLUDE_GEOIDS),
      SupportsBoundingBox(config, pluginData.get_crs_registry(), false),
      SupportsTimeParameters(config),
      SupportsTimeZone(reactor, config)
{
  try
  {
    register_scalar_param<int64_t>(P_FIND_NEAREST_VALID);
    register_scalar_param<std::string>(P_LOCALE);
    register_scalar_param<std::string>(P_MISSING_TEXT);
    register_scalar_param<std::string>(P_PRODUCER);
    register_scalar_param<boost::posix_time::ptime>(P_ORIGIN_TIME, false);
    register_scalar_param<std::string>(P_CRS);
    register_array_param<std::string>(P_ICAO_CODE);

    itsProbabilityConfigParams.probabilityUnit =
        config->get_mandatory_config_param<std::string>("probability_params.probabilityUnit");
    itsProbabilityConfigParams.intensityLight = config->get_mandatory_config_param<std::string>(
        "probability_params.precipitationIntensityLight");
    itsProbabilityConfigParams.intensityModerate = config->get_mandatory_config_param<std::string>(
        "probability_params.precipitationIntensityModerate");
    itsProbabilityConfigParams.intensityHeavy = config->get_mandatory_config_param<std::string>(
        "probability_params.precipitationIntensityHeavy");

    std::vector<std::string> configParamIds =
        config->get_mandatory_config_array<std::string>("probability_params.param_id", 0);

    for (const std::string& paramId : configParamIds)
    {
      ProbabilityConfigParam p;

      p.idLight = get_parameter(config, "probability_params." + paramId + ".idLight");
      p.idModerate = get_parameter(config, "probability_params." + paramId + ".idModerate");
      p.idHeavy = get_parameter(config, "probability_params." + paramId + ".idHeavy");
      p.precipitationType = config->get_mandatory_config_param<std::string>(
          "probability_params." + paramId + ".precipitationType");

      itsProbabilityConfigParams.params.push_back(p);
    }
  }
  catch (...)
  {
    throw Fmi::Exception::Trace(BCP, "Operation failed!");
  }
}

StoredWWProbabilityQueryHandler::~StoredWWProbabilityQueryHandler() {}

void StoredWWProbabilityQueryHandler::parseQueryResults(
    const ProbabilityQueryResultSet& query_results,
    const SmartMet::Spine::BoundingBox& bbox,
    const AirportLocationList& airp_llist,
    const std::string& language,
    SmartMet::Spine::CRSRegistry& crsRegistry,
    const std::string& requestedCRS,
    const boost::posix_time::ptime& origintime,
    const boost::posix_time::ptime& modificationtime,
    const std::string& tz_name,
    CTPP::CDT& hash) const
{
  try
  {
    // create targer spatial reference to get precision and coordinate order
    OGRSpatialReference targetSRS;
    targetSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
    std::string targetURN("urn:ogc:def:crs:" + requestedCRS);
    targetSRS.importFromURN(targetURN.c_str());

    // if geographic projection, precision is 6 digits, otherwise 1 digit
    unsigned int precision = ((targetSRS.IsGeographic()) ? 6 : 1);
    // coordinate order
    bool latLonOrder(targetSRS.EPSGTreatsAsLatLong());

    SmartMet::Spine::BoundingBox query_bbox = bbox;

    // handle lowerCorner and upperCorner
    if (latLonOrder)
    {
      hash["bbox_lower_corner"] = (double2string(query_bbox.yMin, precision) + " " +
                                   double2string(query_bbox.xMin, precision));
      hash["bbox_upper_corner"] = (double2string(query_bbox.yMax, precision) + " " +
                                   double2string(query_bbox.xMax, precision));
      hash["axis_labels"] = "Lat Long";
    }
    else
    {
      hash["bbox_lower_corner"] = (double2string(query_bbox.xMin, precision) + " " +
                                   double2string(query_bbox.yMin, precision));
      hash["bbox_upper_corner"] = (double2string(query_bbox.xMax, precision) + " " +
                                   double2string(query_bbox.yMax, precision));
      hash["axis_labels"] = "Long Lat";
    }

    std::string proj_uri = "UNKNOWN";
    plugin_impl.get_crs_registry().get_attribute(requestedCRS, "projUri", &proj_uri);

    boost::local_time::time_zone_ptr tzp = get_time_zone(tz_name);
    std::string runtime_timestamp = format_local_time(plugin_impl.get_time_stamp(), tzp);

    hash["proj_uri"] = proj_uri;
    hash["fmi_apikey"] = QueryBase::FMI_APIKEY_SUBST;
    hash["fmi_apikey_prefix"] = QueryBase::FMI_APIKEY_PREFIX_SUBST;
    hash["hostname"] = QueryBase::HOSTNAME_SUBST;
    hash["protocol"] = QueryBase::PROTOCOL_SUBST;
    hash["language"] = language;
    hash["response_timestamp"] = runtime_timestamp;

    // wfs members
    hash["wfs_members"] = CTPP::CDT(CTPP::CDT::ARRAY_VAL);

    unsigned int wfs_member_index(0);
    // iterate locations

    for (const AirportLocation& airp_loc : airp_llist)
    {
      SmartMet::Spine::LocationPtr loc = airp_loc.loc;

      WinterWeatherTypeProbabilities ww_tprobs = query_results.at(loc);

      std::vector<std::string> wwtypes;
      // Retrieve all wwtypes
      boost::copy(ww_tprobs | boost::adaptors::map_keys, std::back_inserter(wwtypes));

      for (const std::string& wwtype : wwtypes)
      {
        // own wfs_member for each precipitation type
        CTPP::CDT& wfs_member = hash["wfs_members"][wfs_member_index];
        wfs_member_index++;

        wfs_member["precipitation_form_id"] = wwtype;
        wfs_member["phenomenon_time"] = runtime_timestamp;
        wfs_member["analysis_time"] = format_local_time(origintime, tzp);
        wfs_member["resultTime"] = format_local_time(modificationtime, tzp);
        wfs_member["designator"] = loc->iso2;
        wfs_member["name"] = loc->name;
        wfs_member["location_indicator_icao"] = airp_loc.icao_code;
        wfs_member["field_elevation"] = double2string(loc->elevation, 1);

        auto transformation = crs_registry.create_transformation("EPSG::4326", requestedCRS);

        NFmiPoint p1(loc->longitude, loc->latitude);
        NFmiPoint p2 = transformation->transform(p1);

        double longitude(p2.X());
        double latitude(p2.Y());

        std::stringstream ss_pos;
        if (latLonOrder)
          ss_pos << double2string(latitude, precision) << " "
                 << double2string(longitude, precision);
        else
          ss_pos << double2string(longitude, precision) << " "
                 << double2string(latitude, precision);

        wfs_member["position"] = ss_pos.str();

        CTPP::CDT& result = wfs_member["result"];
        result["timesteps"] = CTPP::CDT(CTPP::CDT::ARRAY_VAL);

        WinterWeatherIntensityProbabilities ww_iprobs = ww_tprobs[wwtype];

        unsigned int timestep_index_index(0);

        for (const WinterWeatherProbability ww_iprob : ww_iprobs)
        {
          std::string timestamp(format_local_time(ww_iprob.timestamp, tzp));

          CTPP::CDT& timestep_data = result["timesteps"][timestep_index_index];
          timestep_index_index++;

          timestep_data["precipitation_form_id"] = wwtype;
          timestep_data["date_time"] = timestamp;
          timestep_data["icao_code"] = airp_loc.icao_code;
          timestep_data["latitude"] = double2string(latitude, precision);
          timestep_data["longitude"] = double2string(longitude, precision);
          timestep_data["type_light"] = itsProbabilityConfigParams.intensityLight;
          timestep_data["type_moderate"] = itsProbabilityConfigParams.intensityModerate;
          timestep_data["type_heavy"] = itsProbabilityConfigParams.intensityHeavy;
          timestep_data["value_light"] = double2string(ww_iprob.light_probability, 1);
          timestep_data["value_moderate"] = double2string(ww_iprob.moderate_probability, 1);
          timestep_data["value_heavy"] = double2string(ww_iprob.heavy_probability, 1);
          timestep_data["unit"] = itsProbabilityConfigParams.probabilityUnit;
        }
      }
    }

    hash["number_matched"] = std::to_string(wfs_member_index);
    hash["number_returned"] = std::to_string(wfs_member_index);
  }
  catch (...)
  {
    throw Fmi::Exception::Trace(BCP, "Operation failed!");
  }
}

WinterWeatherIntensityProbabilities StoredWWProbabilityQueryHandler::getProbabilities(
    const ProbabilityQueryParam& queryParam) const
{
  try
  {
    WinterWeatherIntensityProbabilities ret;

    SmartMet::Engine::Querydata::ParameterOptions qengine_param_light =
        get_qengine_parameter(queryParam, queryParam.paramLight);
    SmartMet::Engine::Querydata::ParameterOptions qengine_param_moderate =
        get_qengine_parameter(queryParam, queryParam.paramModerate);
    SmartMet::Engine::Querydata::ParameterOptions qengine_param_heavy =
        get_qengine_parameter(queryParam, queryParam.paramHeavy);

    SmartMet::Spine::TimeSeries::TimeSeriesPtr qengine_result_light =
        queryParam.q->values(qengine_param_light, queryParam.tlist);
    SmartMet::Spine::TimeSeries::TimeSeriesPtr qengine_result_moderate =
        queryParam.q->values(qengine_param_moderate, queryParam.tlist);
    SmartMet::Spine::TimeSeries::TimeSeriesPtr qengine_result_heavy =
        queryParam.q->values(qengine_param_heavy, queryParam.tlist);

    unsigned int result_set_size(qengine_result_light->size());
    // all intensities should have same timesteps
    if (qengine_result_moderate->size() != result_set_size ||
        qengine_result_heavy->size() != result_set_size)
    {
      std::ostringstream msg;
      msg << "Inconsistent timesteps in LIGHT, MODERATE, HEAVY Winter Weather Condition "
             "parameters: '"
          << queryParam.paramLight.name() << "', '" << queryParam.paramModerate.name() << "', '"
          << queryParam.paramHeavy.name() << "'";
      Fmi::Exception exception(BCP, msg.str());
      exception.addParameter(WFS_EXCEPTION_CODE, WFS_OPERATION_PROCESSING_FAILED);
      throw exception;
    }

    populate_result_vector(
        qengine_result_light, qengine_result_moderate, qengine_result_heavy, ret);

    return ret;
  }
  catch (...)
  {
    throw Fmi::Exception::Trace(BCP, "Operation failed!");
  }
}

void StoredWWProbabilityQueryHandler::query(const StoredQuery& query,
                                            const std::string& language,
					    const boost::optional<std::string>& hostname,
                                            std::ostream& output) const
{
  try
  {
    // 1) query precipitation probabilities from Querydata
    // 2) generate output xml

    const auto& sq_params = query.get_param_map();

    std::string requestedCRS = sq_params.get_single<std::string>(P_CRS);

    std::string targetURN("urn:ogc:def:crs:" + requestedCRS);
    OGRSpatialReference sr;
    sr.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
    if (sr.importFromURN(targetURN.c_str()) != OGRERR_NONE)
    {
      Fmi::Exception exception(BCP, "Invalid crs '" + requestedCRS + "'!");
      exception.addParameter(WFS_EXCEPTION_CODE, WFS_OPERATION_PROCESSING_FAILED);
      throw exception;
    }

    std::string localeName(sq_params.get_single<std::string>(P_LOCALE));
    std::locale outputLocale(localeName.c_str());
    std::string tz_name(get_tz_name(sq_params));

    typedef std::pair<std::string, SmartMet::Spine::LocationPtr> LocationListItem;
    std::list<LocationListItem> llist;

    get_location_options(sq_params, language, &llist);

    if (llist.empty())
    {
      std::cout << "No locations!!!!!!!!!!!!" << std::endl;
      return;
    }

    // bounding box where all locations must reside
    SmartMet::Spine::BoundingBox bbox;
    get_bounding_box(sq_params, requestedCRS, &bbox);

    // if requested CRS and bounding box CRS are different, do transformation
    if (requestedCRS.compare(bbox.crs) != 0)
      bbox = transform_bounding_box(bbox, requestedCRS);

    std::unique_ptr<OGRGeometry> geom;
    Fmi::OGR::CoordinatePoints geoCoordinates;
    geoCoordinates.push_back(std::pair<double, double>(bbox.xMin, bbox.yMin));
    geoCoordinates.push_back(std::pair<double, double>(bbox.xMin, bbox.yMax));
    geoCoordinates.push_back(std::pair<double, double>(bbox.xMax, bbox.yMax));
    geoCoordinates.push_back(std::pair<double, double>(bbox.xMax, bbox.yMin));
    geoCoordinates.push_back(std::pair<double, double>(bbox.xMin, bbox.yMin));
    geom.reset(Fmi::OGR::constructGeometry(geoCoordinates, wkbPolygon, 4326));

    // fetch all icao codes
    std::list<LocationListItem> llist_icao;
    get_location_options(sq_params, "icao", &llist_icao);

    // requested icao codes
    std::vector<std::string> icaoCodeVector;
    sq_params.get<std::string>(P_ICAO_CODE, std::back_inserter(icaoCodeVector));
    std::set<std::string> icaoCodeSet(icaoCodeVector.begin(), icaoCodeVector.end());

    AirportLocationList airp_llist;
    for (LocationListItem item : llist)
    {
      // if location not inside bounding box, ignore it
      OGRPoint point(item.second->longitude, item.second->latitude);
      if (!point.Within(geom.get()))
        continue;

      for (LocationListItem icao_item : llist_icao)
      {
        if (item.second->geoid == icao_item.second->geoid)
        {
          // if icao code not amongst requested, ignore location
          if (icaoCodeSet.size() > 0 &&
              icaoCodeSet.find(icao_item.second->name) == icaoCodeSet.end())
            break;

          airp_llist.push_back(AirportLocation(icao_item.first, item.second));
          break;
        }
      }
    }

    SmartMet::Engine::Querydata::Producer producer = sq_params.get_single<std::string>(P_PRODUCER);
    std::string missingText = sq_params.get_single<std::string>(P_MISSING_TEXT);
    boost::optional<boost::posix_time::ptime> requested_origintime =
        sq_params.get_optional<boost::posix_time::ptime>(P_ORIGIN_TIME);

    SmartMet::Engine::Querydata::Q q;
    if (requested_origintime)
      q = q_engine->get(producer, *requested_origintime);
    else
      q = q_engine->get(producer);

    boost::posix_time::ptime origintime = q->originTime();
    boost::posix_time::ptime modificationtime = q->modificationTime();

    boost::shared_ptr<SmartMet::Spine::TimeSeriesGeneratorOptions> pTimeOptions =
        get_time_generator_options(sq_params);
    // get data in UTC
    const std::string zone = "UTC";
    auto tz = geo_engine->getTimeZones().time_zone_from_string(zone);
    SmartMet::Spine::TimeSeriesGenerator::LocalTimeList tlist =
        SmartMet::Spine::TimeSeriesGenerator::generate(*pTimeOptions, tz);

    ProbabilityQueryResultSet query_results;
    // iterate locations
    for (AirportLocation airp_loc : airp_llist)
    {
      WinterWeatherTypeProbabilities wwprobs;

      // iterate winter weather condition parameters
      for (const ProbabilityConfigParam& configParam : itsProbabilityConfigParams.params)
      {
        // parameter read from configuration file
        SmartMet::Spine::Parameter paramLight(
            "", SmartMet::Spine::Parameter::Type::Data, configParam.idLight);
        SmartMet::Spine::Parameter paramModerate(
            "", SmartMet::Spine::Parameter::Type::Data, configParam.idModerate);
        SmartMet::Spine::Parameter paramHeavy(
            "", SmartMet::Spine::Parameter::Type::Data, configParam.idHeavy);

        ProbabilityQueryParam query_param(paramLight, paramModerate, paramHeavy, q, sr);
        query_param.tlist = tlist;
        query_param.loc = airp_loc.loc;  // item.second;
        query_param.producer = producer;
        query_param.bbox = bbox;
        query_param.missingText = missingText;
        query_param.language = language;
        query_param.outputLocale = outputLocale;
        query_param.tz_name = tz_name;
        query_param.precipitationType = configParam.precipitationType;

        // fetch probabilities for one precipitation type for all intensities
        wwprobs.insert(make_pair(configParam.precipitationType, getProbabilities(query_param)));
      }
      // probabilities for one location for all precipitation types
      query_results.insert(make_pair(airp_loc.loc, wwprobs));
    }

    SmartMet::Spine::CRSRegistry& crsRegistry = plugin_impl.get_crs_registry();
    CTPP::CDT hash;

    parseQueryResults(query_results,
                      bbox,
                      airp_llist,
                      language,
                      crsRegistry,
                      requestedCRS,
                      origintime,
                      modificationtime,
                      tz_name,
                      hash);

    // Format output
    format_output(hash, output, query.get_use_debug_format());
  }
  catch (...)
  {
    throw Fmi::Exception::Trace(BCP, "Operation failed!");
  }
}

}  // namespace WFS
}  // namespace Plugin
}  // namespace SmartMet

namespace
{
using namespace SmartMet::Plugin::WFS;

boost::shared_ptr<SmartMet::Plugin::WFS::StoredQueryHandlerBase>
wfs_winterweather_probabilities_query_handler_create(SmartMet::Spine::Reactor* reactor,
                                                     StoredQueryConfig::Ptr config,
                                                     PluginImpl& pluginData,
                                                     boost::optional<std::string> templateFileName)
{
  try
  {
    StoredWWProbabilityQueryHandler* qh =
        new StoredWWProbabilityQueryHandler(reactor, config, pluginData, templateFileName);
    boost::shared_ptr<SmartMet::Plugin::WFS::StoredQueryHandlerBase> result(qh);
    return result;
  }
  catch (...)
  {
    throw Fmi::Exception::Trace(BCP, "Operation failed!");
  }
}
}  // namespace

SmartMet::Plugin::WFS::StoredQueryHandlerFactoryDef
    wfs_winterweather_probabilities_query_handler_factory(
        &wfs_winterweather_probabilities_query_handler_create);

max_stars_repo_name: nakkim/smartmet-plugin-wfs

max_issues_repo_name: nakkim/smartmet-plugin-wfs

max_forks_repo_name: nakkim/smartmet-plugin-wfs


text:
#ifndef _DEPTH_ESTIMATOR_H_
#define _DEPTH_ESTIMATOR_H_

#include <vector>
#include <Eigen/Dense>
#include <memory>
#include <ceres/ceres.h>

namespace depth_estimator {

  class CostFunction {
  public:
    explicit inline CostFunction(const Eigen::Vector2d& p1, const Eigen::Vector2d& p2,
                                 int idx1, int idx2,
                                 double target_dist)
      : x1_(p1[0]), y1_(p1[1]),
        x2_(p2[0]), y2_(p2[1]),
        idx1_(idx1), idx2_(idx2),
        dist_(target_dist) {}

    template <typename T>
    inline bool operator() (const T* const d, T* residual) const {
      residual[0] = T(dist_) - ceres::sqrt(ceres::pow(d[idx1_] * T(x1_) - d[idx2_] * T(x2_), 2) +
                                           ceres::pow(d[idx1_] * T(y1_) - d[idx2_] * T(y2_), 2) +
                                           ceres::pow(d[idx1_]          - d[idx2_],          2));
      return true;
    }

  private:
    const double x1_;
    const double y1_;
    const double x2_;
    const double y2_;
    const int idx1_;
    const int idx2_;
    const double dist_;
  };

  template <int N>
  class DepthEstimator {
  public:
    inline DepthEstimator(const std::vector<Eigen::Vector2d>& points_2d,
                          const std::vector<Eigen::Vector3d>& points_3d)
      : points_2d_(points_2d), points_3d_(points_3d) {}

    inline Eigen::VectorXd Estimate() {
      std::vector<double> depths(N, 1.0);

      ceres::Problem problem;

      for (int i = N - 1; i >= 1; --i) {
        for (int j = i - 1; j >= 0; --j) {
          problem.AddResidualBlock(
            new ceres::AutoDiffCostFunction<CostFunction, 1, N>(
              new CostFunction(points_2d_[i], points_2d_[j],
                               i, j,
                               (points_3d_[i] - points_3d_[j]).norm())),
            nullptr, depths.data());
        }
      }

      ceres::Solver::Options options;
      options.linear_solver_type = ceres::DENSE_QR;
      options.minimizer_progress_to_stdout = false;
      options.parameter_tolerance = 0.01;

      ceres::Solver::Summary summary;
      ceres::Solve(options, &problem, &summary);

      Eigen::VectorXd result(N);
      for (int i = 0; i < N; ++i) {
        result[i] = depths[i];
      }
      return result;
    }

  private:
    const std::vector<Eigen::Vector2d>& points_2d_;
    const std::vector<Eigen::Vector3d>& points_3d_;
  };

}

#endif


text:
/*=============================================================================
    Copyright (c) 2001-2006 Joel de Guzman
    Copyright (c) 2005-2006 Dan Marsden

    Distributed under the Boost Software License, Version 1.0. (See accompanying
    file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
==============================================================================*/
#if !defined(BOOST_FUSION_ARRAY_27122005_1035)
#define BOOST_FUSION_ARRAY_27122005_1035

#include <boost/fusion/adapted/array/array_iterator.hpp>
#include <boost/fusion/adapted/array/tag_of.hpp>
#include <boost/fusion/adapted/array/detail/is_view_impl.hpp>
#include <boost/fusion/adapted/array/detail/is_sequence_impl.hpp>
#include <boost/fusion/adapted/array/detail/category_of_impl.hpp>
#include <boost/fusion/adapted/array/detail/begin_impl.hpp>
#include <boost/fusion/adapted/array/detail/end_impl.hpp>
#include <boost/fusion/adapted/array/detail/size_impl.hpp>
#include <boost/fusion/adapted/array/detail/at_impl.hpp>
#include <boost/fusion/adapted/array/detail/value_at_impl.hpp>

#endif

max_stars_repo_name: EricBoittier/vina-carb-docker

max_issues_repo_name: EricBoittier/vina-carb-docker

max_forks_repo_name: EricBoittier/vina-carb-docker


text:
/*Copyright (c) 2021 James Gayvert

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
*/

#include <iostream>
#include <cmath>
#include <stdlib.h>
#include <vector>
//#include <boost/math/special_functions/gamma.hpp>  remove for now, keep in case we ever end up using boost
#include <vector>
#include "Shell.h"
#include "boxcap.h"
#include "../asa/asa239.h"

using namespace std;

const double PI = 3.14159265358979;
const double sqrt_pi = 1.77245385090552;
const double OOPI3 = 1.0 / (PI * PI * PI);
const double gauss_factor[] = {1.0, 0.5, 0.75, 1.875, 6.5625, 29.53125, 162.421875};

//
// box-CAP helper functions
//
double intmod_r2(double a, int k, double cl, double cu);
int xi (double z, int i);
double intgauss(double alpha, int n, int = 0, const char* = "message ");
vector<double> polmulti(vector<double>& g, int n, vector<double>& h, int m,
	       int = 0, const char* = "message");
void poltrans (vector<double>& f, double alpha, double beta, int n,
	       int = 0, const char* = "message ");
void gauss2(double ao, double alpha, double bo, double beta,
	    double& co, double& expo, double& factor, int = 0);
void print_vector(vector<double>& v, const char msg[]);

/*
     gaussian integral of the softbox potential
     V  = { (x+boxlength)^2  for x < -boxlength
          {        0         for -boxlength < x < boxlength
          { (x-boxlength)^2  for x > boxlength

This is a Fortran to Cpp translation of our old code from the late 1990ties.
Fortran by Uwe Riss, Robin Santra, and Thomas Sommerfeld
Cpp translation by Thomas Sommerfeld, Nov 2018 T. Sommerfeld

The integrals of the so-called box-CAP are computed analytically.
Do not even try to backwards-engineer this, dig out Robin's paper.

@Article{santra99,
author = "R. Santra and L. S. Cederbaum and H.-D. Meyer",
journal = "Chem. Phys. Lett.",
volume = 303,
year =   1999,
pages =  413,
title =  "Electronic decay of molecular clusters: Non-stationary states computed by standard quantum chemistry methods"}

gto1 and gto2 define two primitive Cartesian GTOs including their l value,
say, px, or dx^2.
boxlength[3] are the cutoff parameters defined above
normalized [0] is a flag 0=unnormalized GTOs, 1=normalized GTOs
verbose can be increased to create debug information

*/

double integrate_box_cap(Shell shell1, Shell shell2,
std::array<size_t,3> l1, std::array<size_t,3> l2, double boxlength[3])
{
  double sum = 0;
	for(size_t prim1 =0;prim1<shell1.num_prims;prim1++)
  {
    gto gto1 = {shell1.origin,shell1.exps[prim1],l1};
    for(size_t prim2=0;prim2<shell2.num_prims;prim2++)
    {
      gto gto2 = {shell2.origin,shell2.exps[prim2],l2};
      sum+= shell1.coeffs[prim1] * shell2.coeffs[prim2] * boxcap(gto1,gto2,boxlength,false,0);
    }
  }
  return sum;
}

double boxcap(gto gto1, gto gto2, double boxlength[3],
	      int normalized, int verbose) {

  int debug = verbose;

  double cen[3] = {0, 0, 0};
  double alpha = gto1.a;
  double beta  = gto2.a;
  double d0, delta, fakt;

  double ovl[3] = {0,0,0};
  double dip[3] = {0,0,0};
  for (int k = 0; k < 3; ++k) {
    int m = gto1.l[k];
    int n = gto2.l[k];
    int nm = n+m;
    double ao = gto1.r[k];
    double bo = gto2.r[k];
    gauss2(ao, alpha,  bo, beta, d0, delta, fakt, verbose);

    vector<double> g (m+1,0);
    vector<double> h (n+1,0);
    g[m] = 1;
    h[n] = 1;
    poltrans(g, ao, d0, m, debug, " g = ");
    poltrans(h, bo, d0, n, debug, " h = ");
    vector<double> f = polmulti(g, m, h, n, debug, " f = ");

    double cl = cen[k] - d0 - boxlength[k];
    double cu = cen[k] - d0 + boxlength[k];

    double temp1 = 0;
    double temp2 = 0;
    for (int i = 0; i <= nm; ++i)
    {
      temp1 += f[i] * intgauss( delta, i, debug, "  I = ");
      temp2 += f[i] * intmod_r2(delta, i, cl, cu);
    }
    ovl[k] = fakt * temp1;
    dip[k] = fakt * temp2;
  }
  if (verbose > 1) {
    cout << "Ovl: " << ovl[0] << ", " <<  ovl[1] << ", " <<  ovl[2] << endl;
    cout << "Dip: " << dip[0] << ", " <<  dip[1] << ", " <<  dip[2] << endl;
  }

  double cap = dip[0]*ovl[1]*ovl[2] + ovl[0]*dip[1]*ovl[2] + ovl[0]*ovl[1]*dip[2];
  if (normalized == 0)
    return cap;
  else {
    alpha *= 2;
    beta *= 2;
    double ovl[3] = {0, 0, 0};
    for (int k=0; k < 3; ++k) {
      int m = 2 * gto1.l[k];
      int n = 2 * gto2.l[k];
      ovl[k] = intgauss(alpha,m) * intgauss(beta,n);
    }
    if (verbose > 1)
      cout << "<0>: " << ovl[0] << ", " <<  ovl[1] << ", " <<  ovl[2] << endl;
    return cap / sqrt(ovl[0]*ovl[1]*ovl[2]);
  }
}

//
//
// integrate(x**k * f(x) * exp(-a*x**2), x, -oo, oo)
//
// where
//
//        (x + cl)**2 if x < cl
// f(x) =    0        if cl <= x <= cu
//        (x - cu)**2 if x > cu
//
// This integral can be expressed in terms of the incomplete gamma-function.
// see:
//  @Article{santra99,
//  author = "R. Santra and L. S. Cederbaum and H.-D. Meyer",
//  journal = "Chem. Phys. Lett.",
//  volume = 303,
//  year =   1999,
//  pages =  413,
//  title =  "Electronic decay of molecular clusters:
//            Non-stationary states computed by standard quantum chemistry methods"}
//
//
double intmod_r2(double a, int k, double cl, double cu) {
  int sg = -1;
  if ((k % 2) == 0)
    sg = 1;
  double ql = a*cl*cl;
  double qu = a*cu*cu;
  double result = 0;
  int ifault = 0;
  for (int i = 0; i < 3; ++i) {
    int j = k + i;
    double par = 0.5*(j + 1);
    int ex = 2 - i;
    double dl = pow(cl,ex);
    double du = pow(-cu,ex);
    double term = tgamma(par)*(sg*dl + du);
    double norm_factor = tgamma(par);
    term -= sg*dl*gammad(ql,par,&ifault)*norm_factor*xi(-cl,j);
    term -=    du*gammad(qu,par,&ifault)*norm_factor*xi( cu,j);
    //term -= sg*dl*boost::math::tgamma_lower(par,ql)*xi(-cl,j);  func is not normalized w/ boost
    //term -=    du*boost::math::tgamma_lower(par,qu)*xi( cu,j);  func is not normalized w/ boost
    term /= pow(a,par);
    if (i == 1)
      term *= 2;
    result += term;
  }
  return 0.5*result;
}

//
//  computes a sign needed in intmod_r2()
//
int xi (double z, int i) {
  if (z >= 0)
    return 1;
  else {
    if ((i%2) == 0)
      return -1;
    else
      return 1;
  }
}

//
//  returns int dx x^n exp(-alpha*x^2)
//
double intgauss(double alpha, int n, int debug, const char* msg) {
  double integral = 0.0;
  if (n%2 == 0) {
    int k = n/2;
    double a = pow(alpha,k);
    integral = gauss_factor[k] * sqrt(PI / alpha) / a;
  }
  if (debug > 3)
    cout << msg << integral << endl;
  return integral;
}

//
//   {sum} a_i x^i  *  {sum} b_i x^i  ->  {sum} c_i x^i
//
vector<double> polmulti(vector<double>& g, int n, vector<double>& h, int m,
			int verbose, const char msg[]) {
  vector<double> f (n+m+1,0);
  for (int i = 0; i <= n; ++i) {
    for (int j = 0; j <= m; ++j) {
      f[i+j] += g[i] * h[j];
    }
  }
  if (verbose > 2) {
    print_vector(f, msg);
  }
  return f;
}

//
// {sum} a_i ( x - a )^i   ->   {sum} b_i ( x - b )^i
//
void poltrans (vector<double> &f, double alpha, double beta, int n,
	       int verbose, const char msg[]) {
  double diff = beta - alpha;
  int k = 1;
  for (int j = 0; j <= n; ++j) {
    double w = 0;
    for (int i = n; i >= j; --i) {
      w = w*diff + f[i];
      f[i] *= (i-j);
    }
    f[j] = w / k;
    k = (j+1)*k;
  }
  if (verbose > 2)
    print_vector(f, msg);
}

void gauss2(double ao, double alpha, double bo, double beta,
	    double& co, double& expo, double& factor, int verbose) {
  expo = alpha + beta;
  co = (ao*alpha + bo*beta) / expo;
  factor = alpha * beta * (bo-ao) * (bo-ao) / expo;
  factor = exp(-factor);
  if (verbose > 2)
    cout << "gs2:" << co << "  " << expo << "  "  << factor << endl;
}

void print_vector(vector<double> &v, const char msg[]) {
  int n = static_cast<int>(v.size());
  cout << msg;
  for (int i = 0; i < n; ++i) {
    cout << v[i];
    if ((i+1)%5 == 0)
      cout << endl;
    else
      cout << ", ";
  }
  if (n%5 != 0)
    cout << endl;
}

max_stars_repo_name: SoubhikM/opencap

max_issues_repo_name: SoubhikM/opencap

max_forks_repo_name: SoubhikM/opencap


text:
/*

Part of Aalto University Game Tools. See LICENSE.txt for licensing info.

*/

#ifndef GENERIC_DENSITY_FOREST_LINEAR_H
#define GENERIC_DENSITY_FOREST_LINEAR_H

#include <Eigen/Eigen>
#include <memory>
#include <list>
#include <deque>
#include <omp.h>
#include <algorithm>
#include "ProbUtils.hpp"
#include "MathUtils.h"

#include <exception>

class ZeroKeyVectorException : std::exception
{
public:
	const char* what() const throw() // my call to the std exception class function (doesn't nessasarily have to be virtual).
	{
		return "The key vector does not match normal."; // my error message
	}
};

typedef float generic_density_forest_scalar;
typedef Eigen::Matrix<generic_density_forest_scalar, Eigen::Dynamic, Eigen::Dynamic> generic_density_forest_matrix;
typedef Eigen::Matrix<generic_density_forest_scalar, Eigen::Dynamic, 1> generic_density_forest_vector;

template<typename first_type>
static bool second_is_smaller(const std::pair<first_type, generic_density_forest_scalar>& first, const std::pair<first_type, generic_density_forest_scalar>& second) {
	return first.second < second.second;
}

template<typename data_type>
static bool is_in_array(const data_type elements[], int end_idx, const data_type element) {
	for (int i = 0; i < end_idx; i++) {
		if (elements[i] == element) {
			return true;
		}
	}

	return false;
}

enum Covariance_computation_mode
{
	AXIS_ALIGNED, FULL
};

enum Split_Optimization_Mode
{
	soINFO_GAIN, soKMEANS
};

template<typename DataType>
class GenericDensityNode
{

	typedef GenericDensityNode<DataType>* NodePointer;
	typedef std::unique_ptr<GenericDensityNode<DataType> > UniqueNodePointer;
private:

	generic_density_forest_vector left_mean_;
	generic_density_forest_vector right_mean_;

	generic_density_forest_matrix left_cov_;
	generic_density_forest_matrix right_cov_;

public:
	static const int MAX_SAMPLES_FOR_INFO_GAIN_MEASURING = 256;

	generic_density_forest_vector separating_normal_;
	generic_density_forest_scalar separating_bias_;
	std::pair<generic_density_forest_vector, generic_density_forest_matrix> gaussian_model_;

	//Parameters for regressor use.
	generic_density_forest_matrix regressor_mean_;
	generic_density_forest_vector regressand_mean_;
	generic_density_forest_matrix regressand_cov_cholesky_;
	generic_density_forest_matrix regression_;

	Covariance_computation_mode covariance_computation_mode_;

	NodePointer parent_;
	UniqueNodePointer left_child_;
	UniqueNodePointer right_child_;

	Split_Optimization_Mode split_mode;
	std::list<std::unique_ptr<DataType> > data_items_;

	const generic_density_forest_vector* (*get_key_vector_)(const DataType& datum);

	void copy_data(GenericDensityNode<DataType>& other) {
		separating_bias_ = other.separating_bias_;
		separating_normal_ = other.separating_normal_;
		gaussian_model_ = other.gaussian_model_;
		regressor_mean_ = other.regressor_mean_;
		regressand_mean_ = other.regressand_mean_;
		regressand_cov_cholesky_ = other.regressand_cov_cholesky_;
		regression_ = other.regression_;
		get_key_vector_ = other.get_key_vector_;
		covariance_computation_mode_ = other.covariance_computation_mode_;
		split_mode = other.split_mode;

		std::list<std::unique_ptr<DataType> >::iterator iter = other.data_items_.begin();

		while (iter != other.data_items_.end()) {
			data_items_.push_back(std::unique_ptr<DataType>(new DataType));
			*(data_items_.back()) = **iter;
			iter++;
		}

	}

	GenericDensityNode() {
		parent_ = nullptr;
		left_child_ = UniqueNodePointer(nullptr);
		right_child_ = UniqueNodePointer(nullptr);
		gaussian_model_ = std::make_pair<generic_density_forest_vector, generic_density_forest_matrix>(generic_density_forest_vector::Zero(0), generic_density_forest_matrix::Zero(0, 0));
		separating_bias_ = std::numeric_limits<generic_density_forest_scalar>::quiet_NaN();
		separating_normal_ = generic_density_forest_vector::Zero(0);
		regressor_mean_ = generic_density_forest_vector::Zero(0);
		regressand_mean_ = generic_density_forest_vector::Zero(0);
		regressand_cov_cholesky_ = generic_density_forest_matrix::Zero(0, 0);
		regression_ = generic_density_forest_matrix::Zero(0, 0);
		covariance_computation_mode_ = AXIS_ALIGNED;
		split_mode = soINFO_GAIN;
		get_key_vector_ = nullptr;
	}

	GenericDensityNode(GenericDensityNode<DataType>& other) {
		parent_ = nullptr;
		left_child_ = UniqueNodePointer(nullptr);
		right_child_ = UniqueNodePointer(nullptr);
		copy_data(other);
	}

	GenericDensityNode<DataType> operator=(GenericDensityNode<DataType>& other) {
		parent_ = nullptr;
		left_child_ = UniqueNodePointer(nullptr);
		right_child_ = UniqueNodePointer(nullptr);
		copy_data(other);

		return *this;
	}

	void build_empty_tree(int depth, int data_dimension, const generic_density_forest_vector* (*key_vector_function)(const DataType&)) {

		get_key_vector_ = key_vector_function;
		separating_bias_ *= 0;
		separating_normal_ = generic_density_forest_vector::Zero(data_dimension);
		left_mean_ = generic_density_forest_vector::Zero(data_dimension);
		right_mean_ = generic_density_forest_vector::Zero(data_dimension);

		if (covariance_computation_mode_ == Covariance_computation_mode::FULL) {
			left_cov_ = generic_density_forest_matrix::Zero(data_dimension, data_dimension);
			right_cov_ = generic_density_forest_matrix::Zero(data_dimension, data_dimension);
		}

		if (covariance_computation_mode_ == Covariance_computation_mode::AXIS_ALIGNED) {
			left_cov_ = generic_density_forest_matrix::Zero(data_dimension, 1);
			right_cov_ = generic_density_forest_matrix::Zero(data_dimension, 1);
		}

		gaussian_model_.first = generic_density_forest_vector::Zero(data_dimension);
		gaussian_model_.second = generic_density_forest_matrix::Zero(data_dimension, data_dimension);

		if (depth < 0) {
			return;
		}

		UniqueNodePointer tmp = UniqueNodePointer(new GenericDensityNode<DataType>());
		tmp->get_key_vector_ = get_key_vector_;
		setLeftChild(tmp);
		left_child_->build_empty_tree(depth - 1, data_dimension, key_vector_function);
		tmp = UniqueNodePointer(new GenericDensityNode<DataType>());
		tmp->get_key_vector_ = get_key_vector_;
		setRightChild(tmp);
		right_child_->build_empty_tree(depth - 1, data_dimension, key_vector_function);

	}

	void setLeftChild(UniqueNodePointer& newLeftChild) {

		if (newLeftChild.get()) {
			if (left_child_.get()) {
				left_child_->parent_ = nullptr;
			}
			left_child_ = std::move(newLeftChild);
			left_child_->parent_ = this;
		}
		else {
			left_child_ = UniqueNodePointer(nullptr);
		}

	}

	void setRightChild(UniqueNodePointer& newRightChild) {

		if (newRightChild.get()) {
			if (right_child_.get()) {
				right_child_->parent_ = nullptr;
			}
			right_child_ = std::move(newRightChild);
			right_child_->parent_ = this;
		}
		else {
			right_child_ = UniqueNodePointer(nullptr);
		}
	}

	bool is_empty_leaf(void) {
		if (!hasChildren()) {

			if (data_items_.size() == 0) {
				return true;
			}

		}
		return false;
	}

	NodePointer get_leaf(const generic_density_forest_vector &key_vector) {
		NodePointer node = this;
		NodePointer prev_node = node;
		while (node) {
			prev_node = node;

			if (node->left_or_right(key_vector)) {
				node = node->right_child_.get();
			}
			else {
				node = node->left_child_.get();
			}

		}
		return prev_node;
	}

	NodePointer get_leaf(const DataType& datum) {
		const generic_density_forest_vector &key_vector = *(get_key_vector_(datum));
		return get_leaf(key_vector);
	}

	void form_regression_parameters(int regressor_dim) {

		int regressand_dim = gaussian_model_.first.size() - regressor_dim;

		assert(regressand_dim > 0);

		regressor_mean_ = gaussian_model_.first.tail(regressor_dim);
		regressand_mean_ = gaussian_model_.first.head(regressand_dim);

		generic_density_forest_matrix cov_xy = gaussian_model_.second.topRightCorner(regressand_dim, regressor_dim);
		generic_density_forest_matrix cov_yy = gaussian_model_.second.bottomRightCorner(regressor_dim, regressor_dim);
		generic_density_forest_matrix cov_yy_inv = invert_positive_definite(cov_yy);

		regression_ = cov_xy*cov_yy_inv;

		regressand_cov_cholesky_ = gaussian_model_.second.topLeftCorner(regressand_dim, regressand_dim) - cov_xy*cov_yy_inv*cov_xy.transpose();
		regressand_cov_cholesky_ = make_invertable_positive_definite(regressand_cov_cholesky_).llt().matrixL();

	}

	bool hasChildren() {
		if (right_child_.get() && left_child_.get()) {
			return true;
		}
		else {
			return false;
		}
	}

	generic_density_forest_scalar information_gain(std::list<DataType* >& left_data, std::list<DataType* >& right_data) {

		if (left_data.size() == 0 || right_data.size() == 0) {
			return -std::numeric_limits<generic_density_forest_scalar>::infinity();
		}

		const generic_density_forest_vector& tmp_key_vector = *get_key_vector_(**left_data.begin());
		int data_dim = tmp_key_vector.size();

		if (left_cov_.rows() != data_dim || right_cov_.rows() != data_dim) {

			if (covariance_computation_mode_ == Covariance_computation_mode::FULL) {
				left_cov_ = generic_density_forest_matrix::Zero(data_dim, data_dim);
				right_cov_ = generic_density_forest_matrix::Zero(data_dim, data_dim);
			}

			if (covariance_computation_mode_ == Covariance_computation_mode::AXIS_ALIGNED) {
				left_cov_ = generic_density_forest_matrix::Zero(data_dim, 1);
				right_cov_ = generic_density_forest_matrix::Zero(data_dim, 1);
			}
		}
		else {
			left_cov_ *= 0;
			right_cov_ *= 0;
		}

		if (left_mean_.size() != data_dim || right_mean_.size() != data_dim) {
			left_mean_ = generic_density_forest_vector::Zero(data_dim);
			right_mean_ = generic_density_forest_vector::Zero(data_dim);
		}
		else {
			left_mean_ *= 0;
			right_mean_ *= 0;
		}

		int left_size = 0;
		int right_size = 0;

		for (DataType* datum : left_data) {
			const generic_density_forest_vector& sample_key_vector = *get_key_vector_(*datum);
			left_mean_ += sample_key_vector;
			left_size++;
		}
		left_mean_ /= (generic_density_forest_scalar)left_size;

		for (DataType* datum : right_data) {
			const generic_density_forest_vector& sample_key_vector = *get_key_vector_(*datum);
			right_mean_ += sample_key_vector;
			right_size++;
		}
		right_mean_ /= (generic_density_forest_scalar)right_size;

		generic_density_forest_vector diff;

		for (DataType* datum : left_data) {
			const generic_density_forest_vector& sample_key_vector = *get_key_vector_(*datum);
			diff = sample_key_vector - left_mean_;
			if (covariance_computation_mode_ == Covariance_computation_mode::FULL) {
				left_cov_ += diff*diff.transpose();
			}

			if (covariance_computation_mode_ == Covariance_computation_mode::AXIS_ALIGNED) {
				left_cov_ += diff.cwiseAbs2();
			}
		}
		left_cov_ /= (generic_density_forest_scalar)left_size;

		for (DataType* datum : right_data) {
			const generic_density_forest_vector& sample_key_vector = *get_key_vector_(*datum);
			diff = sample_key_vector - right_mean_;

			if (covariance_computation_mode_ == Covariance_computation_mode::FULL) {
				right_cov_ += diff*diff.transpose();
			}

			if (covariance_computation_mode_ == Covariance_computation_mode::AXIS_ALIGNED) {
				right_cov_ += diff.cwiseAbs2();
			}
		}
		right_cov_ /= (generic_density_forest_scalar)right_size;

		int joint_size = left_size + right_size;

		generic_density_forest_scalar info_gain = (generic_density_forest_scalar)0;
		if (covariance_computation_mode_ == Covariance_computation_mode::FULL) {
			info_gain -= (generic_density_forest_scalar)left_size / (generic_density_forest_scalar)joint_size*log_determinant_positive_definite(left_cov_);
			assert(finiteNumber(info_gain));
			info_gain -= (generic_density_forest_scalar)right_size / (generic_density_forest_scalar)joint_size*log_determinant_positive_definite(right_cov_);
			assert(finiteNumber(info_gain));
		}

		if (covariance_computation_mode_ == Covariance_computation_mode::AXIS_ALIGNED) {
			for (int i = 0; i < left_cov_.size(); i++) {
				left_cov_(i) = std::max(std::numeric_limits<generic_density_forest_scalar>::min(), left_cov_(i));
			}
			info_gain -= (generic_density_forest_scalar)left_size / (generic_density_forest_scalar)joint_size*(left_cov_.array().log().sum());
			//assert(finiteNumber(info_gain));
			for (int i = 0; i < right_cov_.size(); i++) {
				right_cov_(i) = std::max(std::numeric_limits<generic_density_forest_scalar>::min(), right_cov_(i));
			}
			info_gain -= (generic_density_forest_scalar)right_size / (generic_density_forest_scalar)joint_size*(right_cov_.array().log().sum());
			//assert(finiteNumber(info_gain));
		}

		return info_gain;

	}

	//generic_density_forest_scalar information_gain_maximum_separation(std::vector<DataType*> left_data, std::vector<DataType*> right_data){

	//	generic_density_forest_scalar margin = std::numeric_limits<generic_density_forest_scalar>::infinity();
	//	generic_density_forest_scalar current_margin = std::numeric_limits<generic_density_forest_scalar>::infinity();

	//	for (DataType* datum : left_data){
	//		current_margin = std::abs(separating_bias_ + separating_normal_.dot(*get_key_vector_(*datum)));
	//		margin = std::min(margin,current_margin);
	//	}

	//	for (DataType* datum : right_data){
	//		current_margin = std::abs(separating_bias_ + separating_normal_.dot(*get_key_vector_(*datum)));
	//		margin = std::min(margin,current_margin);
	//	}

	//	return margin*left_data.size()*right_data.size();

	//}

	////Obsolete
	//std::vector<std::pair<DataType*,generic_density_forest_scalar> > project_data(std::vector<DataType*>& data, generic_density_forest_vector& start, generic_density_forest_vector& direction){
	//	std::vector<std::pair<DataType*,generic_density_forest_scalar> > data_with_projection;
	//	data_with_projection.reserve(data.size());
	//	for (int i = 0; i < (int)data.size(); i++){
	//		const generic_density_forest_vector* location = get_key_vector_(*(data[i]));
	//		generic_density_forest_scalar alpha = projection_of_point_to_line(*(location),start,direction);
	//		data_with_projection.push_back(std::make_pair(data[i],alpha));
	//	}
	//	return data_with_projection;
	//}

	////Obsolete
	////It is assumed that the <gaussian_model_> is computed.
	//void get_separating_hyper_plane_perpendicular_to_first_principal_axis(std::vector<DataType*>& data){

	//	int data_dim = (*get_key_vector_(*data[0])).size();

	//	if (data.size() < 2){
	//		separating_normal_ = generic_density_forest_vector::Zero(data_dim);
	//		separating_bias_ = 0;
	//		return;
	//	}

	//	separating_normal_ = get_principal_axes(gaussian_model_.second).col(0);

	//	left_cov_ = generic_density_forest_matrix::Zero(data_dim,1);
	//	right_cov_ = generic_density_forest_matrix::Zero(data_dim,1);
	//	left_mean_ = generic_density_forest_vector::Zero(data_dim);
	//	right_mean_ = generic_density_forest_vector::Zero(data_dim);

	//	//Project the data to the first principal axis of the data.
	//	std::vector<std::pair<DataType*,generic_density_forest_scalar> > data_with_projection = project_data(data,gaussian_model_.first,separating_normal_);

	//	//Sort
	//	std::sort(data_with_projection.begin(),data_with_projection.end(),second_is_smaller<DataType*>);

	//	int left_size = 0;
	//	int right_size = data.size();
	//	generic_density_forest_scalar joint_size = (generic_density_forest_scalar)(left_size + right_size);
	//	right_mean_ = gaussian_model_.first;
	//	right_cov_ = gaussian_model_.second.diagonal();

	//	generic_density_forest_scalar max_info_gain = -std::numeric_limits<generic_density_forest_scalar>::infinity();
	//	generic_density_forest_scalar current_info_gain = -std::numeric_limits<generic_density_forest_scalar>::infinity();
	//	int max_gain_index = 0; // The index of the first element that should be left on the

	//	generic_density_forest_vector left_square_sum = left_mean_*0;
	//	right_mean_ *= (generic_density_forest_scalar)right_size;
	//	generic_density_forest_vector right_square_sum = right_mean_*0;
	//	for (int i = 0; i < (int)data.size();i++){
	//		right_square_sum += get_key_vector_(*data[i])->cwiseAbs2();
	//	}

	//	//NB! left_mean_ and right_mean_ are actually used to store the sum instead of mean in this loop!!!
	//	for (int i = 0; i < (int)data_with_projection.size()-1; i++){
	//		current_info_gain = 0;

	//		const generic_density_forest_vector* location = get_key_vector_(*(data_with_projection[i].first));
	//		left_square_sum += location->cwiseAbs2();
	//		left_mean_ += *location;
	//		right_square_sum -= location->cwiseAbs2();
	//		right_mean_ -= *location;

	//		left_size++;
	//		right_size--;

	//		left_cov_ = left_square_sum - (left_mean_.cwiseAbs2())/(generic_density_forest_scalar)left_size;
	//		left_cov_ /= (generic_density_forest_scalar)left_size;
	//		right_cov_ = right_square_sum - (right_mean_.cwiseAbs2())/(generic_density_forest_scalar)right_size;
	//		right_cov_ /= (generic_density_forest_scalar)right_size;

	//		for (int i = 0; i < left_cov_.size(); i++){
	//			left_cov_(i) = std::max(std::numeric_limits<generic_density_forest_scalar>::min(),left_cov_(i));
	//		}
	//		current_info_gain -= (generic_density_forest_scalar)left_size/(generic_density_forest_scalar)joint_size*(left_cov_.array().log().sum());
	//		for (int i = 0; i < right_cov_.size(); i++){
	//			right_cov_(i) = std::max(std::numeric_limits<generic_density_forest_scalar>::min(),right_cov_(i));
	//		}
	//		current_info_gain -= (generic_density_forest_scalar)right_size/(generic_density_forest_scalar)joint_size*(right_cov_.array().log().sum());

	//		if (current_info_gain > max_info_gain){
	//			max_info_gain = current_info_gain;
	//			max_gain_index = i;
	//		}

	//	}

	//	const generic_density_forest_vector* location1 = get_key_vector_(*(data_with_projection[max_gain_index].first));
	//	const generic_density_forest_vector* location2 = get_key_vector_(*(data_with_projection[max_gain_index+1].first));
	//	generic_density_forest_scalar bias1 = -location1->dot(separating_normal_);
	//	generic_density_forest_scalar bias2 = -location2->dot(separating_normal_);

	//	separating_bias_ = bias1 + (bias2-bias1)*sampleUniform<generic_density_forest_scalar>();

	//}

	std::vector<DataType*> compute_gaussian_models(void) {

		std::vector<DataType*> data;

		if (!hasChildren()) {
			for (const std::unique_ptr<DataType>& tmp : data_items_) {
				data.push_back(tmp.get());
			}
			assert(data.size() > 0);
		}
		else {
			std::vector<DataType*> tmp_data = left_child_->compute_gaussian_models();
			for (DataType* datum : tmp_data) {
				data.push_back(datum);
			}
			assert(data.size() > 0);

			tmp_data = right_child_->compute_gaussian_models();
			for (DataType* datum : tmp_data) {
				data.push_back(datum);
			}
			assert(data.size() > 0);
		}

		assert(data.size() > 0);

		//Build a gaussian model of all the incoming data
		const generic_density_forest_vector& tmp_key_vector = *get_key_vector_(*(data[0]));
		gaussian_model_.first = tmp_key_vector * 0;
		gaussian_model_.second = generic_density_forest_matrix::Zero(gaussian_model_.first.size(), gaussian_model_.first.size());

		for (size_t i = 0; i < data.size(); i++) {
			const generic_density_forest_vector& sample_key_vector = *get_key_vector_(*(data[i]));
			gaussian_model_.first += sample_key_vector;
		}
		gaussian_model_.first /= (generic_density_forest_scalar)data.size();

		generic_density_forest_vector diff;
		for (size_t i = 0; i < data.size(); i++) {
			const generic_density_forest_vector& sample_key_vector = *get_key_vector_(*(data[i]));
			diff = sample_key_vector - gaussian_model_.first;
			gaussian_model_.second += diff*diff.transpose();
		}
		gaussian_model_.second /= (generic_density_forest_scalar)data.size();

		return data;

	}

	void build_tree(int tries, int minimum_data_in_leaf, std::list<std::unique_ptr<DataType> > data) {

		data_items_.clear();

		UniqueNodePointer tmp = UniqueNodePointer(nullptr);
		setLeftChild(tmp);
		tmp = UniqueNodePointer(nullptr);
		setRightChild(tmp);

		//If the split breaks too few data points away, this branch is ready.
		if (data.size() <= (unsigned)minimum_data_in_leaf) {

			data_items_ = std::move(data);

			return;

		}

		if (split_mode == soKMEANS)
		{
			//split plane through kmeans with 2 clusters and downsampling of data
			const generic_density_forest_vector& vect_tmp = *get_key_vector_(**data.begin());
			int keyDim = vect_tmp.size();
			generic_density_forest_vector mean[2], sum[2];
			float w[2];
			for (int i = 0; i < 2; i++)
			{
				mean[i] = generic_density_forest_vector(keyDim);
				sum[i] = generic_density_forest_vector(keyDim);
				mean[i].setZero();
				sum[i].setZero();
				w[i] = 0;
			}
			int nSamples = std::min(tries, (int)data.size());

			int counter = 0;
			//first init clusters randomly
			for (std::unique_ptr<DataType>& datum : data)
			{
				const generic_density_forest_vector &key = *get_key_vector_(*datum);
				int clusterIdx;
				if (w[0] == 0)
					clusterIdx = 0;
				else if (w[1] == 0)
					clusterIdx = 1;
				else
					clusterIdx = AaltoGames::rand01();
				sum[clusterIdx] += key;
				w[clusterIdx]++;
				counter++;
				if (counter == nSamples) {
					break;
				}
			}
			for (int i = 0; i < 2; i++)
			{
				if (w[i] != 0)
					mean[i] = sum[i] / w[i];
				sum[i].setZero();
				w[i] = 0;
			}

			//now do the k-means iteration
			const int kMeansIter = 3;
			for (int iter = 0; iter < kMeansIter; iter++)
			{
				counter = 0;
				for (std::unique_ptr<DataType>& datum : data)
				{
					const generic_density_forest_vector &key = *get_key_vector_(*datum);
					int clusterIdx;
					float dist0 = (key - mean[0]).norm();
					float dist1 = (key - mean[1]).norm();
					if (dist0 < dist1)
						clusterIdx = 0;
					else if (dist1 < dist0)
						clusterIdx = 1;
					else
					{
						if (w[0] == 0)
							clusterIdx = 0;
						else if (w[1] == 0)
							clusterIdx = 1;
						else
							clusterIdx = AaltoGames::rand01();
					}
					sum[clusterIdx] += key;
					w[clusterIdx]++;
					counter++;
					if (counter == nSamples) {
						break;
					}
				}
				for (int i = 0; i < 2; i++)
				{
					if (w[i] != 0)
						mean[i] = sum[i] / w[i];
					sum[i].setZero();
					w[i] = 0;
				}
			}
			std::pair<generic_density_forest_scalar, generic_density_forest_vector> hyper_plane;
			hyper_plane.second = (mean[1] - mean[0]).normalized();
			separating_normal_ = hyper_plane.second;
			hyper_plane.first = -separating_normal_.dot(((generic_density_forest_scalar)0.5)*(mean[0] + mean[1]));
			separating_bias_ = hyper_plane.first;
		}
		else
		{
			std::pair<generic_density_forest_scalar, generic_density_forest_vector> hyper_plane = get_separating_hyperplane(data);
			separating_bias_ = hyper_plane.first;
			separating_normal_ = hyper_plane.second;

			std::list<DataType* > left_samples;
			std::list<DataType* > right_samples;

			if (data.size() > MAX_SAMPLES_FOR_INFO_GAIN_MEASURING)
				divide_data_sampled(data, left_samples, right_samples, MAX_SAMPLES_FOR_INFO_GAIN_MEASURING);
			else
				divide_data_info(data, left_samples, right_samples);

			generic_density_forest_scalar max_information_gain = information_gain(left_samples, right_samples);
			//generic_density_forest_scalar max_information_gain = information_gain_maximum_separation(separated_data.first,separated_data.second);
			for (int i = 0; i < tries; i++) {

				std::pair<generic_density_forest_scalar, generic_density_forest_vector> tmp_plane = get_separating_hyperplane(data);
				separating_bias_ = tmp_plane.first;
				separating_normal_ = tmp_plane.second;
				if (data.size() > MAX_SAMPLES_FOR_INFO_GAIN_MEASURING)
					divide_data_sampled(data, left_samples, right_samples, MAX_SAMPLES_FOR_INFO_GAIN_MEASURING);
				else
					divide_data_info(data, left_samples, right_samples);

				generic_density_forest_scalar tmp_information_gain = information_gain(left_samples, right_samples);
				//generic_density_forest_scalar tmp_information_gain = information_gain_maximum_separation(separated_data.first,separated_data.second);

				if (tmp_information_gain > max_information_gain) {
					max_information_gain = tmp_information_gain;
					hyper_plane = tmp_plane;
				}

			}

			separating_bias_ = hyper_plane.first;
			separating_normal_ = hyper_plane.second;

			//std::pair<std::vector<DataType*>,std::vector<DataType*> > separated_data;
			//get_separating_hyper_plane_perpendicular_to_first_principal_axis(data);

		}

		std::list<std::unique_ptr<DataType> > left_data;
		std::list<std::unique_ptr<DataType> > right_data;
		divide_data(data, left_data, right_data);

		//If the split breaks too few data points away, this branch is ready.
		if (left_data.size() == 0 || right_data.size() == 0) {

			data_items_.clear();
			for (std::unique_ptr<DataType>& datum : left_data) {
				std::unique_ptr<DataType> tmp_ptr = std::unique_ptr<DataType>(nullptr);
				datum.swap(tmp_ptr);
				data_items_.push_back(std::move(tmp_ptr));
			}

			for (std::unique_ptr<DataType>& datum : right_data) {
				std::unique_ptr<DataType> tmp_ptr = std::unique_ptr<DataType>(nullptr);
				datum.swap(tmp_ptr);
				data_items_.push_back(std::move(tmp_ptr));
			}

			return;

		}

		tmp = UniqueNodePointer(new GenericDensityNode<DataType>());
		tmp->get_key_vector_ = get_key_vector_;
		setLeftChild(tmp);
		left_child_->build_tree(tries, minimum_data_in_leaf, std::move(left_data));

		tmp = UniqueNodePointer(new GenericDensityNode<DataType>());
		tmp->get_key_vector_ = get_key_vector_;
		setRightChild(tmp);
		right_child_->build_tree(tries, minimum_data_in_leaf, std::move(right_data));

	}

	void add_sample(int tries, int minimum_data_in_leaf, int maximum_data_in_leaf, std::unique_ptr<DataType> datum) {

		data_items_.push_back(std::move(datum));

		if ((int)data_items_.size() > maximum_data_in_leaf) {
			std::list<std::unique_ptr<DataType> > new_data_items = std::move(data_items_);
			data_items_.clear();
			build_tree(tries, minimum_data_in_leaf, std::move(new_data_items));
		}

	}

	////Obsolete
	//void fill_tree(int tries, int minimum_data_in_leaf, std::vector<DataType*>& data){

	//	data_items_.clear();

	//	//Build a gaussian model of all the incoming data
	//	gaussian_model_.first = *get_key_vector_(*(data[0])) * 0;
	//	gaussian_model_.second = generic_density_forest_matrix::Zero(gaussian_model_.first.size(),gaussian_model_.first.size());

	//	for (size_t i = 0; i < data.size(); i++){
	//		gaussian_model_.first += *get_key_vector_(*(data[i]));
	//	}
	//	gaussian_model_.first /= (generic_density_forest_scalar)data.size();

	//	generic_density_forest_vector diff;
	//	for (size_t i = 0; i < data.size(); i++){
	//		diff = *get_key_vector_(*(data[i])) - gaussian_model_.first;
	//		gaussian_model_.second += diff*diff.transpose();
	//	}
	//	gaussian_model_.second /= (generic_density_forest_scalar)data.size();

	//	//If the split breaks too few data points away, this branch is ready.
	//	if (data.size() < (unsigned)minimum_data_in_leaf){

	//		left_child_.reset();
	//		right_child_.reset();

	//		data_items_.clear();
	//		data_item_pointers_.clear();
	//		if (data_storage_mode_ == Storage_mode::COPY){
	//			data_items_.reserve(data.size());
	//			for (size_t i = 0; i < data.size(); i++){
	//				data_items_.push_back(std::move(*data[i]));
	//			}
	//		}

	//		if (data_storage_mode_ == Storage_mode::POINTER){
	//			data_item_pointers_.reserve(data.size());
	//			for (size_t i = 0; i < data.size(); i++){
	//				data_item_pointers_.push_back(data[i]);
	//			}
	//		}

	//		return;

	//	}

	//	std::pair<generic_density_forest_scalar,generic_density_forest_vector> hyper_plane = get_separating_hyperplane(data);
	//	separating_bias_ = hyper_plane.first;
	//	separating_normal_ = hyper_plane.second;

	//	std::pair<std::vector<DataType*>,std::vector<DataType*> > separated_data;
	//	divide_data(data,separated_data.first,separated_data.second);

	//	generic_density_forest_scalar max_information_gain = information_gain(separated_data.first,separated_data.second);
	//	//generic_density_forest_scalar max_information_gain = information_gain_maximum_separation(separated_data.first,separated_data.second);
	//	for (int i = 0; i < tries; i++){

	//		std::pair<generic_density_forest_scalar,generic_density_forest_vector> tmp_plane = get_separating_hyperplane(data);
	//		separating_bias_ = tmp_plane.first;
	//		separating_normal_ = tmp_plane.second;

	//		divide_data(data,separated_data.first,separated_data.second);

	//		generic_density_forest_scalar tmp_information_gain = information_gain(separated_data.first,separated_data.second);
	//		//generic_density_forest_scalar tmp_information_gain = information_gain_maximum_separation(separated_data.first,separated_data.second);

	//		if (tmp_information_gain > max_information_gain){
	//			max_information_gain = tmp_information_gain;
	//			hyper_plane = tmp_plane;
	//		}

	//	}

	//	separating_bias_ = hyper_plane.first;
	//	separating_normal_ = hyper_plane.second;

	//	divide_data(data,separated_data.first,separated_data.second);

	//	left_child_->fill_tree(tries,minimum_data_in_leaf,separated_data.first);
	//	right_child_->fill_tree(tries,minimum_data_in_leaf,separated_data.second);

	//}

	std::pair<generic_density_forest_scalar, generic_density_forest_vector> get_separating_hyperplane(const std::list<std::unique_ptr<DataType> >& data) {

		if (data.size() < 1) {
			throw ZeroKeyVectorException();
			return std::make_pair(0, generic_density_forest_vector::Zero(0));
		}

		const generic_density_forest_vector& tmp_key_vect = *get_key_vector_(**data.begin());
		generic_density_forest_vector normal = generic_density_forest_vector::Random(tmp_key_vect.size());
		normal.normalize();

		int tmp_idx_1 = rand() % data.size();
		int tmp_idx_2 = rand() % data.size();
		while (tmp_idx_2 == tmp_idx_1) {
			tmp_idx_2 = rand() % data.size();
		}

		DataType* element1 = nullptr;
		DataType* element2 = nullptr;

		//if (data.size() <= 1){
		//	std::cout << "Too few elements.";
		//}

		int counter = 0;
		for (const std::unique_ptr<DataType>& datum : data) {

			if (counter == tmp_idx_1) {
				element1 = datum.get();
			}

			if (counter == tmp_idx_2) {
				element2 = datum.get();
			}

			if (element1 && element2) {
				break;
			}
			counter++;
		}

		//generic_density_forest_vector first_random_element = -*get_key_vector_(*(data[tmp_idx_1]));
		//generic_density_forest_vector second_random_element = -*get_key_vector_(*(data[tmp_idx_2]));

		//std::cout << element1 << " " << element2 << std::endl;

		const generic_density_forest_vector key_vect1 = *get_key_vector_(*element1);
		const generic_density_forest_vector key_vect2 = *get_key_vector_(*element2);
		generic_density_forest_scalar bias_1 = normal.dot(-key_vect1);
		generic_density_forest_scalar bias_2 = normal.dot(-key_vect2);

		generic_density_forest_scalar bias = bias_1 + (bias_2 - bias_1)*sampleUniform<generic_density_forest_scalar>();

		return std::make_pair(bias, normal);

	}

	//True signals right, false signals left.
	bool left_or_right(const generic_density_forest_vector& datum) {

		generic_density_forest_scalar decision = 0;

		if (separating_normal_.size() == 0) {
			if (AaltoGames::rand01() == 0) {
				return false;
			}
			else {
				return true;
			}
		}

		if (datum.size() != separating_normal_.size()) {
			throw ZeroKeyVectorException();
		}

		decision = separating_normal_.dot(datum) + separating_bias_;

		if (decision > 0) {
			return true;
		}
		else if (decision < 0) {
			return false;
		}
		else {
			if (AaltoGames::rand01() == 0) {
				return false;
			}
			else {
				return true;
			}
		}

	}

	//The data will be split to left and right data, which are stored to <left_data> and <right_data>.
	void divide_data_info(std::list<std::unique_ptr<DataType> >& data, std::list<DataType* >& left_data, std::list<DataType* >& right_data) {

		left_data.clear();;
		right_data.clear();

		for (const std::unique_ptr<DataType>& datum : data) {

			DataType* datum_ptr = datum.get();
			const generic_density_forest_vector& key_vector = *get_key_vector_(*datum_ptr);

			if (left_or_right(key_vector)) {
				right_data.push_back(datum_ptr);
			}
			else {
				left_data.push_back(datum_ptr);
			}
		}

	}

	//The data will be split to left and right data, which are stored to <left_data> and <right_data>.
	void divide_data(std::list<std::unique_ptr<DataType> >& data, std::list<std::unique_ptr<DataType> >& left_data, std::list<std::unique_ptr<DataType> >& right_data) {

		left_data.clear();;
		right_data.clear();

		while (data.size() > 0) {
			std::unique_ptr<DataType> datum = std::unique_ptr<DataType>(nullptr);
			datum.swap(data.back());
			data.pop_back();

			const generic_density_forest_vector& key_vect = *get_key_vector_(*datum);

			if (left_or_right(key_vect)) {
				right_data.push_back(std::move(datum));
			}
			else {
				left_data.push_back(std::move(datum));
			}
		}

		data.clear();

	}

	//The data will be split to left and right data, which are stored to <left_data> and <right_data>.
	void divide_data_sampled(std::list<std::unique_ptr<DataType> >& data, std::list<DataType* >& left_data, std::list<DataType* >& right_data, int nSamples) {

		left_data.clear();;
		right_data.clear();

		if (data.size() < 1) {
			return;
		}

		while ((int)(right_data.size() + left_data.size()) < nSamples) {

			int idx = AaltoGames::randInt(0, (int)data.size() - 1);

			std::list<std::unique_ptr<DataType> >::iterator iter = data.begin();

			int counter = 0;
			while (counter < idx && iter != data.end()) {
				iter++;
				counter++;
			}

			if (iter == data.end()) {
				continue;
			}

			std::unique_ptr<DataType>& datum = *iter;

			//try{
			const generic_density_forest_vector& key_vect = *get_key_vector_(*datum);

			if (left_or_right(key_vect)) {
				right_data.push_back(datum.get());
			}
			else {
				left_data.push_back(datum.get());
			}
			//}
			//catch (std::exception e){
			//	std::cout << data.size() << " " << left_data.size() << " " << right_data.size() << std::endl;
			//	std::cout << counter << std::endl;
			//	std::cout << datum << std::endl;
			//	std::cout << datum.get()->state_ << std::endl;
			//	std::cout << datum.get()->control_ << std::endl;
			//	std::cout << datum.get()->future_state_ << std::endl;
			//	std::cout << datum.get()->key_vector_ << std::endl;
			//}

		}

	}

};

template<typename DataType>
class GenericDensityTree {

	typedef GenericDensityNode<DataType>* NodePointer;
	typedef GenericDensityTree* TreePointer;
	typedef std::unique_ptr<GenericDensityNode<DataType> > UniqueNodePointer;

public:

	UniqueNodePointer root_;
	//std::list<DensityNode > mNodes; //This is just a container to hold the nodes.

	GenericDensityTree() {
		//mNodes.clear();
		root_ = UniqueNodePointer(new GenericDensityNode<DataType>());
	}

	GenericDensityTree(GenericDensityTree& otherTree) {

		root_ = UniqueNodePointer(new GenericDensityNode<DataType>(*(otherTree.root_)));

		std::vector<GenericDensityNode<DataType>* > nodesOther;
		nodesOther.push_back(&(*(otherTree.root_)));

		std::vector<GenericDensityNode<DataType>* > nodes;
		nodes.push_back(&(*root_));

		while (nodes.size() > 0) {
			NodePointer nodeOther = nodesOther.back();
			nodesOther.pop_back();

			NodePointer node = nodes.back();
			nodes.pop_back();

			UniqueNodePointer newLeft = UniqueNodePointer(nullptr);
			if (nodeOther->left_child_.get()) {
				newLeft = UniqueNodePointer(new GenericDensityNode<DataType>(*(nodeOther->left_child_.get())));
			}
			UniqueNodePointer newRight = UniqueNodePointer(nullptr);
			if (nodeOther->right_child_.get()) {
				newRight = UniqueNodePointer(new GenericDensityNode<DataType>(*(nodeOther->right_child_.get())));
			}

			node->setLeftChild(newLeft);
			node->setRightChild(newRight);

			if (nodeOther->left_child_.get()) {
				nodesOther.push_back(nodeOther->left_child_.get());
				nodes.push_back(node->left_child_.get());
			}

			if (nodeOther->right_child_.get()) {
				nodesOther.push_back(nodeOther->right_child_.get());
				nodes.push_back(node->right_child_.get());
			}

		}

	}

	GenericDensityTree operator=(GenericDensityTree& otherTree) {

		if (&otherTree == this) {
			return *this;
		}

		root_ = UniqueNodePointer(new GenericDensityNode<DataType>(*(otherTree.root_)));

		std::vector<GenericDensityNode<DataType>* > nodesOther;
		nodesOther.push_back(&(*(otherTree.root_)));

		std::vector<GenericDensityNode<DataType>* > nodes;
		nodes.push_back(&(*root_));

		while (nodes.size() > 0) {
			NodePointer nodeOther = nodesOther.back();
			nodesOther.pop_back();

			NodePointer node = nodes.back();
			nodes.pop_back();

			UniqueNodePointer newLeft = UniqueNodePointer(nullptr);
			if (nodeOther->left_child_.get()) {
				newLeft = UniqueNodePointer(new GenericDensityNode<DataType>(*(nodeOther->left_child_)));
			}
			UniqueNodePointer newRight = UniqueNodePointer(nullptr);
			if (nodeOther->right_child_.get()) {
				newRight = UniqueNodePointer(new GenericDensityNode<DataType>(*(nodeOther->right_child_)));
			}

			node->setLeftChild(newLeft);
			node->setRightChild(newRight);

			if (nodeOther->left_child_.get()) {
				nodesOther.push_back(nodeOther->left_child_.get());
				nodes.push_back(node->left_child_.get());
			}

			if (nodeOther->right_child_.get()) {
				nodesOther.push_back(nodeOther->right_child_.get());
				nodes.push_back(node->right_child_.get());
			}

		}

		return *this;
	}

	DataType* get_random() {

		NodePointer node = root_.get();
		NodePointer next = root_.get();

		while (next) {
			node = next;
			if (rand() % 2 == 0) {
				next = node->left_child_.get();
			}
			else {
				next = node->right_child_.get();
			}
		}

		DataType* nearest = nullptr;

		if (node->data_items_.size() > 0) {

			int rand_idx = rand() % node->data_items_.size();
			std::list<std::unique_ptr<DataType> >::iterator iter = node->data_items_.begin();

			while (rand_idx > 0 && iter != node->data_items_.end()) {
				rand_idx--;
				iter++;
			}

			if (iter == node->data_items_.end()) {
				iter--;
			}

			nearest = iter->get();

		}

		return nearest;

	}

	DataType* get_approximate_nearest(const DataType& datum) {

		NodePointer node = root_->get_leaf(datum);

		generic_density_forest_scalar nearest_dist = std::numeric_limits<generic_density_forest_scalar>::infinity();
		generic_density_forest_scalar current_dist = std::numeric_limits<generic_density_forest_scalar>::infinity();

		const generic_density_forest_vector& key_vector = *(node->get_key_vector_(datum));

		DataType* nearest = nullptr;

		for (const std::unique_ptr<DataType>& tmp : node->data_items_) {
			const generic_density_forest_vector& data_key_vector = *(node->get_key_vector_(*tmp));
			current_dist = (key_vector - data_key_vector).norm();

			if (current_dist < nearest_dist) {
				nearest_dist = current_dist;
				nearest = tmp.get();
			}
		}

		return nearest;

	}

	std::vector<DataType*> get_neighborhood(const DataType& datum) {

		NodePointer node = root_->get_leaf(datum);
		std::vector<DataType*> neighbors;

		for (const std::unique_ptr<DataType>& tmp : node->data_items_) {

			neighbors.push_back(tmp.get());

		}

		return neighbors;

	}

	int get_neighborhood(const DataType& datum, DataType* data_ptr_array[], int max_end_index) {

		NodePointer node = root_->get_leaf(datum);
		int end_idx = 0;

		if (node->data_storage_mode_ == Storage_mode::COPY) {
			for (const std::unique_ptr<DataType>& tmp : node->data_items_) {
				data_ptr_array[end_idx] = tmp.get();
				end_idx++;
				if (max_end_index == end_idx) {
					break;
				}
			}
		}

		return end_idx;

	}

	void rebuild_tree(int tries, int minimum_data_in_leaf, int remember_instances = -1) {
		std::list<std::unique_ptr<DataType> > data;

		std::vector<NodePointer> nodes_to_visit;
		nodes_to_visit.push_back(root_.get());

		while (nodes_to_visit.size() > 0) {
			NodePointer node = nodes_to_visit.back();
			nodes_to_visit.pop_back();
			if (node) {
				nodes_to_visit.push_back(node->left_child_.get());
				nodes_to_visit.push_back(node->right_child_.get());

				for (std::unique_ptr<DataType>& tmp : node->data_items_) {
					data.push_back(std::move(tmp));
				}
			}

		}

		if (remember_instances > 0) {
			while (remember_instances < (int)data.size()) {

				std::list<std::unique_ptr<DataType> >::iterator iter = data.begin();
				int rand_idx = rand() % data.size();

				while (rand_idx > 0) {
					rand_idx--;

					iter++;

					if (iter == data.end()) {
						break;
					}

				}
				if (iter == data.end()) {
					iter--;
				}

				data.erase(iter);

			}
		}

		const generic_density_forest_vector* (*key_vector_fun)(const DataType& datum) = root_->get_key_vector_;

		UniqueNodePointer tmp = UniqueNodePointer(new GenericDensityNode<DataType>());
		tmp->get_key_vector_ = key_vector_fun;
		root_.swap(tmp);

		//std::cout << "Tree has nodes: " << data.size() << std::endl;

		root_->build_tree(tries, minimum_data_in_leaf, std::move(data));

	}

	//Assuming that data_ptr_array has room for k pointers.
	void get_up_to_k_nearest(const DataType& datum, DataType* data_ptr_array[], int k) {
		generic_density_forest_vector key_vector = *(root_->get_key_vector_(datum));
		get_up_to_k_nearest(key_vector, data_ptr_array, k);
	}

	void get_up_to_k_nearest(const generic_density_forest_vector &key_vector, DataType* data_ptr_array[], int k, bool findInSibling = false) {

		NodePointer node = root_->get_leaf(key_vector);
		if (findInSibling)
		{
			if (node->parent_ != nullptr)
			{
				NodePointer left = node->parent_->left_child_.get();
				NodePointer right = node->parent_->right_child_.get();
				node = left == node ? right : left;
				node = node->get_leaf(key_vector);
			}
		}

		auto first_is_closer = [&](DataType* datum1, DataType* datum2) {
			if (!datum1 && !datum2) {
				return false;
			}
			if (datum1 && !datum2) {
				return true;
			}
			if (!datum1 && datum2) {
				return false;
			}

			const generic_density_forest_vector& key_vect1 = *(root_->get_key_vector_(*datum1));
			const generic_density_forest_vector& key_vect2 = *(root_->get_key_vector_(*datum2));

			float dist_to_1 = (key_vector - key_vect1).norm();
			float dist_to_2 = (key_vector - key_vect2).norm();
			return dist_to_1 < dist_to_2;
		};

		for (const std::unique_ptr<DataType>& tmp : node->data_items_) {
			DataType* ptr = tmp.get();

			bool is_in_array_already = false;
			for (int i = 0; i < k; i++) {
				if (data_ptr_array[i]) {
					if (*data_ptr_array[i] == *ptr) {
						is_in_array_already = true;
						break;
					}
				}
			}
			if (is_in_array_already) {
				continue;
			}

			if (first_is_closer(ptr, data_ptr_array[k - 1])) {
				data_ptr_array[k - 1] = ptr;
				std::sort(data_ptr_array, &(data_ptr_array[k]), first_is_closer);
			}

		}

	}

	//DataType* get_approximate_nearest(const DataType& datum){

	//	NodePointer node = get_leaf(datum);

	//	return &(node->data_items_[rand()%node->data_items_.size()]);

	//}

	void remove_empty_leaf(NodePointer node) {

		if (!(node->is_empty_leaf())) {
			return;
		}

		//Node is root
		if (!(node->parent_)) {
			return;
		}

		NodePointer grand_parent = nullptr;
		NodePointer parent = node->parent_;
		if (parent) {
			grand_parent = parent->parent_;
		}
		UniqueNodePointer sibling = nullptr;
		if (parent->left_child_.get() == node) {
			sibling = std::move(parent->right_child_);
		}
		else {
			sibling = std::move(parent->left_child_);
		}

		//The sibling should become the new root
		if (!grand_parent) {
			sibling->parent_ = nullptr;
			root_ = std::move(sibling);
			return;
		}

		assert(root_->hasChildren());

		//The sibling replaces parent
		if (grand_parent->left_child_.get() == parent) {
			grand_parent->setLeftChild(sibling);
			assert(grand_parent->hasChildren());
		}
		else {
			grand_parent->setRightChild(sibling);
			assert(grand_parent->hasChildren());
		}

	}

	std::vector<DataType*> find_data_item_in_tree_exhaustive_search(const DataType& datum) {

		std::vector<NodePointer> nodes;
		nodes.push_back(root_.get());

		std::vector<DataType*> data_items;

		while (nodes.size() > 0) {

			NodePointer node = nodes[nodes.size() - 1];
			nodes.pop_back();

			if (node) {
				NodePointer tmp_ptr = node->left_child_.get();
				if (tmp_ptr) {
					nodes.push_back(tmp_ptr);
				}

				tmp_ptr = node->right_child_.get();
				if (tmp_ptr) {
					nodes.push_back(tmp_ptr);
				}

			}

			std::list<std::unique_ptr<DataType> >::iterator iter = node->data_items_.begin();

			while (iter != node->data_items_.end()) {

				if ((**iter) == datum) {
					data_items.push_back((*iter).get());
				}

				iter++;

			}

		}

		return data_items;

	}

	void remove_data_point(const DataType& datum) {

		std::vector<NodePointer> nodes;
		nodes.push_back(root_.get());

		std::vector<NodePointer> leaves;

		while (nodes.size() > 0) {

			NodePointer node = nodes[nodes.size() - 1];
			nodes.pop_back();

			if (node) {
				bool has_child = false;
				NodePointer tmp_ptr = node->left_child_.get();
				if (tmp_ptr) {
					nodes.push_back(tmp_ptr);
					has_child = true;
				}

				tmp_ptr = node->right_child_.get();
				if (tmp_ptr) {
					nodes.push_back(tmp_ptr);
					has_child = true;
				}

				if (!has_child) {
					leaves.push_back(node);
				}

			}

			std::list<std::unique_ptr<DataType> >::iterator iter = node->data_items_.begin();

			while (iter != node->data_items_.end()) {

				if ((**iter) == datum) {
					iter = node->data_items_.erase(iter);
				}
				else {
					iter++;
				}
			}

		}

		for (NodePointer ptr : leaves) {
			remove_empty_leaf(ptr);
		}

	}

	void build_tree_copying(int tries, int minimum_data_in_leaf, std::vector<DataType >& data) {

		std::list<std::unique_ptr<DataType> > tmp;

		for (unsigned i = 0; i < data.size(); i++) {
			tmp.push_back(std::unique_ptr<DataType>(new DataType));
			*(tmp.back()) = data[i];
		}

		root_->build_tree(tries, minimum_data_in_leaf, std::move(tmp));

	}

	void add_sample(int tries, int minimum_data_in_leaf, int maximum_data_in_leaf, DataType datum) {

		NodePointer node = root_->get_leaf(datum);

		std::unique_ptr<DataType> datum_ptr = std::unique_ptr<DataType>(new DataType());
		*datum_ptr = datum;

		node->add_sample(tries, minimum_data_in_leaf, maximum_data_in_leaf, std::move(datum_ptr));

	}

	generic_density_forest_vector sample_conditioned(const generic_density_forest_vector& regressor) {

		generic_density_forest_vector key = root_->gaussian_model_.first*std::numeric_limits<generic_density_forest_scalar>::quiet_NaN();
		key.tail(regressor.size()) = regressor;

		NodePointer node = root_.get();

		while (node->hasChildren()) {

			key.head(node->regressand_mean_.size()) = node->regressand_mean_ + node->regression_*(regressor - node->regressor_mean_);
			key.head(node->regressand_mean_.size()) += node->regressand_cov_cholesky_*BoxMuller<generic_density_forest_scalar>(node->regressand_mean_.size());

			if (node->left_or_right(key)) {
				node = node->right_child_.get();
			}
			else {
				node = node->left_child_.get();
			}

		}

		return key;

	}

	void form_regression(int regressor_dim) {

		root_->compute_gaussian_models();

		std::vector<NodePointer> nodes;
		nodes.push_back(root_.get());

		while (nodes.size() > 0) {

			NodePointer node = nodes[nodes.size() - 1];
			nodes.pop_back();

			if (node) {
				node->form_regression_parameters(regressor_dim);
				nodes.push_back(node->left_child_.get());
				nodes.push_back(node->right_child_.get());
			}

		}

	}

};

template<typename DataType>
class GenericDensityForest {

	typedef GenericDensityNode<DataType>* NodePointer;
	typedef GenericDensityTree<DataType>* TreePointer;

private:

public:

	std::vector<GenericDensityTree<DataType> > forest_;

	GenericDensityForest() {
		forest_.push_back(GenericDensityTree<DataType>());
	}

	GenericDensityForest(unsigned int numberOfTrees) {
		for (unsigned int treeNum = 0; treeNum < numberOfTrees; treeNum++) {
			forest_.push_back(GenericDensityTree<DataType>());
		}
	}

	GenericDensityForest(GenericDensityForest& other) {
		forest_ = other.forest_;
	}

	GenericDensityForest<DataType> operator=(const GenericDensityForest<DataType>& other) {
		forest_ = other.forest_;
		return *this;
	}

	void set_key_vector_function(const generic_density_forest_vector* (*key_vector_function)(const DataType&)) {
		for (size_t i = 0; i < forest_.size(); i++) {
			forest_[i].root_->get_key_vector_ = key_vector_function;
		}
	}

	void set_split_optimization_mode(Split_Optimization_Mode mode) {
		for (size_t i = 0; i < forest_.size(); i++) {
			forest_[i].root_->split_mode = mode;
		}
	}

	void add_sample(int tries, int minimum_data_in_leaf, int maximum_data_in_leaf, DataType datum) {
		for (size_t i = 0; i < forest_.size(); i++) {
			forest_[i].add_sample(tries, minimum_data_in_leaf, maximum_data_in_leaf, datum);
		}
	}

	DataType* get_approximate_nearest(const DataType& datum) {

		std::vector<DataType*> nearest_ones(forest_.size(), nullptr);

		//#pragma omp parallel for
		for (size_t i = 0; i < (int)forest_.size(); i++) {
			nearest_ones[i] = forest_[i].get_approximate_nearest(datum);
		}

		NodePointer random_node = forest_[0].root_.get();
		const generic_density_forest_vector& key_vector = *(random_node->get_key_vector_(datum));

		int nearest_idx = 0;
		generic_density_forest_scalar nearest_dist = std::numeric_limits<generic_density_forest_scalar>::infinity();
		generic_density_forest_scalar current_dist = std::numeric_limits<generic_density_forest_scalar>::infinity();

		for (size_t i = 0; i < nearest_ones.size(); i++) {

			if (nearest_ones[i]) {

				const generic_density_forest_vector& comparison_key_vector = *(random_node->get_key_vector_(*(nearest_ones[i])));
				current_dist = (key_vector - comparison_key_vector).norm();

				if (current_dist < nearest_dist) {
					nearest_dist = current_dist;
					nearest_idx = i;
				}

			}

		}

		return nearest_ones[nearest_idx];

	}

	std::vector<DataType*> get_neighborhood(const DataType& datum) {

		std::vector<DataType*> nearest_ones;

		for (size_t i = 0; i < (int)forest_.size(); i++) {
			std::vector<DataType*> nearest_ones_tree;
			nearest_ones_tree = forest_[i].get_neighborhood(datum);
			for (size_t j = 0; j < nearest_ones_tree.size(); j++) {
				nearest_ones.push_back(nearest_ones_tree[j]);
			}
		}

		return nearest_ones;

	}

	//Assuming that data_ptr_array has room for k pointers.
	int get_up_to_k_nearest(const DataType& datum, DataType* data_ptr_array[], int k) {
		const generic_density_forest_vector& key_vector = *forest_[0].root_->get_key_vector_(datum);
		return get_up_to_k_nearest(key_vector, data_ptr_array, k);
	}

	int get_up_to_k_nearest(const generic_density_forest_vector &key_vector, DataType* data_ptr_array[], int k) {

		for (int i = 0; i < k; i++) {
			data_ptr_array[i] = nullptr;
		}

		if (k == 0) {
			return 0;
		}

		for (size_t i = 0; i < (int)forest_.size(); i++) {
			forest_[i].get_up_to_k_nearest(key_vector, data_ptr_array, k);
			forest_[i].get_up_to_k_nearest(key_vector, data_ptr_array, k, true);  //for robustness against splits that leave a child with only few samples, also search the sibling
		}

		int found_count = 0;

		for (int i = 0; i < k; i++) {
			if (data_ptr_array[i]) {
				found_count++;
			}
		}

		return found_count;

	}

	int get_neighborhood(const DataType& datum, DataType* nearest_ones[], int max_amount) {

		int end_idx = 0;

		for (size_t i = 0; i < (int)forest_.size(); i++) {
			end_idx += forest_[i].get_neighborhood(datum, &nearest_ones[end_idx], max_amount - end_idx);
		}

		return end_idx;

	}

	//DataType* get_approximate_nearest(const DataType& datum){

	//	return forest_[rand()%forest_.size()].get_approximate_nearest(datum);

	//}

	void build_forest(int tries, int minimum_data_in_leaf, std::vector<std::unique_ptr<DataType> > data) {
		for (int i = 0; i < (int)forest_.size(); i++) {
			GenericDensityTree<DataType>& tree = forest_[i];
			tree.root_->build_tree(tries, minimum_data_in_leaf, data);
		}
	}

	void build_forest_copying(int tries, int minimum_data_in_leaf, std::vector<DataType >& data) {
		for (int i = 0; i < (int)forest_.size(); i++) {
			GenericDensityTree<DataType>& tree = forest_[i];
			tree.build_tree_copying(tries, minimum_data_in_leaf, data);
		}
	}

	generic_density_forest_vector sample_conditioned(const generic_density_forest_vector& regressor) {

		return forest_[rand() % forest_.size()].sample_conditioned(regressor);

	}

	void form_regressions(int regressor_dim) {
		for (int i = 0; i < (int)forest_.size(); i++) {
			GenericDensityTree<DataType>& tree = forest_[i];
			tree.form_regression(regressor_dim);
		}

	}

	void claim_ownership_of_forest(GenericDensityForest<DataType>& other_forest) {
		forest_ = std::move(other_forest.forest_);
	}

	void claim_ownership_of_tree(GenericDensityTree<DataType>& tree, size_t number_of_trees_in_this_forest) {
		if (forest_.size() < number_of_trees_in_this_forest) {
			forest_.push_back(std::move(tree));
			return;
		}

		if (forest_.size() > number_of_trees_in_this_forest) {
			forest_.pop_back();
		}

		int swap_idx = rand() % forest_.size();
		//forest_[swap_idx].root_.reset();
		//forest_[swap_idx].root_ = std::move(tree.root_);
		forest_[swap_idx].root_.swap(tree.root_);
	}

	void remove_data_point(DataType datum) {
		for (int i = 0; i < (int)forest_.size(); i++) {
			GenericDensityTree<DataType>& tree = forest_[i];
			tree.remove_data_point(datum);
		}
	}

	DataType* get_random() {
		return forest_[rand() % forest_.size()].get_random();
	}

	bool forest_ready_to_use(void) {

		if (forest_.size() == 0) {
			return false;
		}

		if (forest_[0].root_->separating_normal_.size() == 0 && forest_[0].root_->data_items_.size() == 0) {
			return false;
		}

		return true;

	}

};

#endif