text:
// Copyright 2016 Damien Buhl (alias daminetreg) for Fr. Sauter AG, CH-4016 BASEL
// 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)
#ifndef MAPBOX_DETAIL_BOOST_SPIRIT_ATTRIBUTES_HPP
#define MAPBOX_DETAIL_BOOST_SPIRIT_ATTRIBUTES_HPP

#include <mapbox/variant.hpp>

#include <mapbox/traits/is_mapbox_variant.hpp>
#include <mapbox/traits/is_type_in_variant.hpp>

#include <boost/spirit/home/support/attributes.hpp>

// Defines the trait for Boost.Spirit to treat mapbox::util::variant as a variant.
namespace boost { namespace spirit { namespace traits
{
    template <typename Domain, class... Types>
    struct not_is_variant<mapbox::util::variant<Types...>, Domain>
      : mpl::false_
    {};
}
}}

namespace boost { namespace spirit { namespace traits {

    template <class... Types>
    struct variant_which< mapbox::util::variant<Types...> >
    {
        static int call(mapbox::util::variant<Types...> const& v)
        {
            return v.which();
        }
    };

}}}

namespace boost { namespace spirit { namespace traits {

    template <typename Variant, typename Expected>
    struct compute_compatible_component_variant<Variant, Expected, mpl::false_
      , typename enable_if< typename mapbox::traits::is_mapbox_variant<Variant>::type >::type>
    {
        typedef typename traits::variant_type<Variant>::type variant_type;

        // true_ if the attribute matches one of the types in the variant
        typedef mapbox::traits::is_type_in_variant<Expected, Variant> type;
        enum { value = type::value };

        // return the type in the variant the attribute is compatible with
        typedef typename
            mpl::eval_if<type, mpl::identity<Expected>, mpl::identity<unused_type> >::type
        compatible_type;

        // return whether the given type is compatible with the Expected type
        static bool is_compatible(int which)
        {
          auto idx = (type::size - 1 - type::index); // Typelist is inverted and 0-based
          return which == idx;
        }
    };

}}}

#endif

ext: hpp

max_stars_repo_name: sauter-hq/boost-spirit-mapbox-variant


text:
#include "utils/string.h"

#include <sstream>
#include <fstream>
#include <algorithm>
#include <functional>
#include <codecvt>
#include <locale>
#include <cctype>

#include <boost/algorithm/string.hpp>

namespace util{
namespace string{

bool isspace(std::string const& str){
    return std::find_if_not(str.cbegin(),str.cend(), [](auto x){return std::isspace(x);})==str.cend();
}

std::vector<std::string> split(std::string words, const char *delim){
    std::vector<std::string> tokens;
    boost::split(tokens, words, boost::is_any_of(delim));
    return tokens;
}

std::string join(std::vector<std::string> words, std::string delim){
    std::stringstream s;
    auto n = words.size();
    for(decltype(n)i=0; i!=n; ++i){
        s<<words[i];
        if(i!=n-1) s<<delim;
    }
    return s.str();
}
std::string strip(std::string const& str){
    auto beg=std::find_if_not(str.cbegin(), str.cend(), [](auto x){return std::isspace(x);});
    auto end=std::find_if_not(str.crbegin(), str.crend(), [](auto x){return std::isspace(x);});
    return str.substr(beg-str.cbegin(), end.base()-beg);
}
std::string tolower(std::string str){
    std::transform(str.begin(), str.end(), str.begin(), ::tolower);
    return str;
}

std::vector<std::string> readlines(std::string file){
    std::ifstream is{file};
    return readlines(std::move(is));
}
std::vector<std::string> readlines(std::istream &&is){
    std::vector<std::string> lines;
    std::string line;
    while(std::getline(is, line)) lines.push_back(line);
    return lines;
}

void write_whole(std::string filename, std::string const& content){
    std::ofstream t(filename);
    t << content;
}
std::string read_whole(std::string file){
    std::ifstream t(file);
    std::stringstream buffer;
    buffer << t.rdbuf();
    return buffer.str();
}

std::vector<char> pack_words(std::vector<std::string> const &words){
    std::vector<char> vec;
    for(auto const &x:words){
        std::copy(x.cbegin(),x.cend(),std::back_inserter(vec));
        vec.push_back('\0');
    }
    return vec;
}
std::vector<char> pack_words(std::vector<const char*> const &words){
    std::vector<char> vec;
    for(auto const &word:words){
        std::string x{word};
        std::copy(x.cbegin(),x.cend(),std::back_inserter(vec));
        vec.push_back('\0');
    }
    return vec;
}

std::vector<std::string> unpack_words(std::vector<char> const &concat_words){
    std::vector<std::string> words;
    auto it =concat_words.cbegin();
    auto end=concat_words.cend();
    while(it!=end){
        words.push_back(std::string{&(*it)});
        it=std::find(it, end, '\0');
        it=std::find_if_not(it, end, [](auto x){return x=='\0';});
    }
    return words;
}

std::vector<const char *> unpack_word_views(std::vector<char> const &concat_words){
    std::vector<const char *> words;
    auto it =concat_words.cbegin();
    auto end=concat_words.cend();
    while(it!=end){
        words.push_back(&(*it));
        it=std::find(it, end, '\0');
        it=std::find_if_not(it, end, [](auto x){return x=='\0';});
    }
    return words;
}

std::vector<util::cstring_span<>> unpack_tokenized_sentence(util::cstring_span<> sentence){
    util::cstring_span<> sent{sentence};
    std::vector<util::cstring_span<>> words;
    auto beg = sent.cbegin();
    auto end = sent.cend();
    auto it=beg;
    while(it!=end){
        auto word_beg=it-beg;
        it=std::find(it, end, ' ');
        auto word_end=it-beg;
        words.push_back(sent.subspan(word_beg, word_end-word_beg));
        it=std::find_if_not(it, end, [](auto x){return x==' ';});
    }
    return words;
}

std::string substring_unicode_offset(std::string str, int beg, int end){
    std::wstring_convert<std::codecvt_utf8<wchar_t>> to_utf8;
    std::wstring wstr = to_utf8.from_bytes(str);
    auto wsubstr = wstr.substr(beg, end-beg);
    auto substr = to_utf8.to_bytes(wsubstr);
    return substr;
}

}//namespace util::string
}//namespace util

ext: cpp

max_stars_repo_name: uphere-co/nlp-prototype

max_issues_repo_name: uphere-co/nlp-prototype

max_forks_repo_name: uphere-co/nlp-prototype


text:
/**
 * @file
 * This file is part of SeisSol.
 *
 * @author Carsten Uphoff (c.uphoff AT tum.de, http://www5.in.tum.de/wiki/index.php/Carsten_Uphoff,_M.Sc.)
 * @author Sebastian Wolf (wolf.sebastian AT in.tum.de, https://www5.in.tum.de/wiki/index.php/Sebastian_Wolf,_M.Sc.)
 *
 * @section LICENSE
 * Copyright (c) 2017 - 2020, SeisSol Group
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice,
 *    this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 *    this list of conditions and the following disclaimer in the documentation
 *    and/or other materials provided with the distribution.
 *
 * 3. Neither the name of the copyright holder nor the names of its
 *    contributors may be used to endorse or promote products derived from this
 *    software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 *
 * @section DESCRIPTION
 *
 **/
#include <Eigen/Eigenvalues>
#include <Kernels/precision.hpp>
#include <Initializer/typedefs.hpp>

#include <PUML/PUML.h>
#include <PUML/Downward.h>
#include <PUML/Upward.h>
#include "LtsWeights.h"

#include <Eigen/Dense>

#include <Initializer/ParameterDB.h>
#include <Parallel/MPI.h>

#include <generated_code/tensor.h>
#include <generated_code/init.h>

class FaceSorter {
private:
	std::vector<PUML::TETPUML::face_t> const& m_faces;

public:
	FaceSorter(std::vector<PUML::TETPUML::face_t> const& faces) : m_faces(faces) {}

	bool operator()(unsigned int a, unsigned int b) const {
		return m_faces[a].gid() < m_faces[b].gid();
	}
};

void seissol::initializers::time_stepping::LtsWeights::computeMaxTimesteps( PUML::TETPUML const&  mesh,
                                                                            std::vector<double> const& pWaveVel,
                                                                            std::vector<double>& timestep ) {
  std::vector<PUML::TETPUML::cell_t> const& cells = mesh.cells();
  std::vector<PUML::TETPUML::vertex_t> const& vertices = mesh.vertices();

  for (unsigned cell = 0; cell < cells.size(); ++cell) {
    // Compute insphere radius
    Eigen::Vector3d barycentre(0.,0.,0.);
    Eigen::Vector3d x[4];
    unsigned vertLids[4];
    PUML::Downward::vertices(mesh, cells[cell], vertLids);
    for (unsigned vtx = 0; vtx < 4; ++vtx) {
      for (unsigned d = 0; d < 3; ++d) {
        x[vtx](d) = vertices[ vertLids[vtx] ].coordinate()[d];
      }
    }
    Eigen::Matrix4d A;
    A << x[0](0), x[0](1), x[0](2), 1.0,
         x[1](0), x[1](1), x[1](2), 1.0,
         x[2](0), x[2](1), x[2](2), 1.0,
         x[3](0), x[3](1), x[3](2), 1.0;

    double alpha = A.determinant();
    double Nabc = ( (x[1]-x[0]).cross(x[2]-x[0]) ).norm();
    double Nabd = ( (x[1]-x[0]).cross(x[3]-x[0]) ).norm();
    double Nacd = ( (x[2]-x[0]).cross(x[3]-x[0]) ).norm();
    double Nbcd = ( (x[2]-x[1]).cross(x[3]-x[1]) ).norm();
    double insphere = std::fabs(alpha) / (Nabc + Nabd + Nacd + Nbcd);

    // Compute maximum timestep (CFL=1)
    timestep[cell] = 2.0 * insphere / (pWaveVel[cell] * (2*CONVERGENCE_ORDER-1));
  }
}

int seissol::initializers::time_stepping::LtsWeights::getCluster( double    timestep,
                                                                  double    globalMinTimestep,
                                                                  unsigned  rate ) {
  if (rate == 1) {
    return 0;
  }

  double upper;
  upper = rate * globalMinTimestep;

  int cluster = 0;
  while (upper <= timestep) {
    upper *= rate;
    ++cluster;
  }
  return cluster;
}

int seissol::initializers::time_stepping::LtsWeights::getBoundaryCondition( int const* boundaryCond,
                                                                            unsigned cell,
                                                                            unsigned face ) {
  int bcCurrentFace = ((boundaryCond[cell] >> (face*8)) & 0xFF);
  if (bcCurrentFace > 64) {
     bcCurrentFace = 3;
  }
  return bcCurrentFace;
}

int seissol::initializers::time_stepping::LtsWeights::ipow(int x, int y) {
  assert(y >= 0);

  if (y == 0) {
    return 1;
  }
  int result = x;
  while(--y) {
    result *= x;
  }
  return result;
}

void seissol::initializers::time_stepping::LtsWeights::computeWeights(PUML::TETPUML const& mesh) {
  logInfo(seissol::MPI::mpi.rank()) << "Computing LTS weights.";

  std::vector<PUML::TETPUML::cell_t> const& cells = mesh.cells();
  int const* boundaryCond = mesh.cellData(1);

  std::vector<double> pWaveVel;
  pWaveVel.resize(cells.size());

  seissol::initializers::ElementBarycentreGeneratorPUML queryGen(mesh);
  //up to now we only distinguish between anisotropic elastic any other isotropic material
#ifdef USE_ANISOTROPIC
  std::vector<seissol::model::AnisotropicMaterial> materials(cells.size());
  seissol::initializers::MaterialParameterDB<seissol::model::AnisotropicMaterial> parameterDB;
#else
  std::vector<seissol::model::ElasticMaterial> materials(cells.size());
  seissol::initializers::MaterialParameterDB<seissol::model::ElasticMaterial> parameterDB;
#endif
  parameterDB.setMaterialVector(&materials);
  parameterDB.evaluateModel(m_velocityModel, queryGen);
  for(unsigned cell = 0; cell < cells.size(); ++cell) {
    pWaveVel[cell] = materials[cell].getMaxWaveSpeed();
  }
  std::vector<double> timestep;
  timestep.resize(cells.size());
  computeMaxTimesteps(mesh, pWaveVel, timestep);

  double localMinTimestep = *std::min_element(timestep.begin(), timestep.end());
  double localMaxTimestep = *std::max_element(timestep.begin(), timestep.end());
  double globalMinTimestep;
  double globalMaxTimestep;
#ifdef USE_MPI
  MPI_Allreduce(&localMinTimestep, &globalMinTimestep, 1, MPI_DOUBLE, MPI_MIN, seissol::MPI::mpi.comm());
  MPI_Allreduce(&localMaxTimestep, &globalMaxTimestep, 1, MPI_DOUBLE, MPI_MAX, seissol::MPI::mpi.comm());
#else
  globalMinTimestep = localMinTimestep;
  globalMaxTimestep = localMaxTimestep;
#endif

  int* cluster = new int[cells.size()];
  for (unsigned cell = 0; cell < cells.size(); ++cell) {
    cluster[cell] = getCluster(timestep[cell], globalMinTimestep, m_rate);
  }

  int totalNumberOfReductions = enforceMaximumDifference(mesh, cluster);

  delete[] m_vertexWeights;
  //m_ncon = 2;
  m_ncon = 1;
  m_vertexWeights = new int[cells.size() * m_ncon];
  int maxCluster = getCluster(globalMaxTimestep, globalMinTimestep, m_rate);
  int drToCellRatio = 1;
  for (unsigned cell = 0; cell < cells.size(); ++cell) {
    int dynamicRupture = 0;
    for (unsigned face = 0; face < 4; ++face) {
      dynamicRupture += ( getBoundaryCondition(boundaryCond, cell, face) == 3) ? 1 : 0;
    }

    m_vertexWeights[m_ncon * cell] = (1 + drToCellRatio*dynamicRupture) * ipow(m_rate, maxCluster - cluster[cell]);
    //m_vertexWeights[m_ncon * cell + 1] = (dynamicRupture > 0) ? 1 : 0;
  }

  delete[] cluster;

  logInfo(seissol::MPI::mpi.rank()) << "Computing LTS weights. Done. " << utils::nospace << '(' << totalNumberOfReductions << " reductions.)";
}

int seissol::initializers::time_stepping::LtsWeights::enforceMaximumDifference(PUML::TETPUML const& mesh, int* cluster) {
  int totalNumberOfReductions = 0;
  int globalNumberOfReductions;
  do {
    int localNumberOfReductions = enforceMaximumDifferenceLocal(mesh, cluster);

#ifdef USE_MPI
    MPI_Allreduce(&localNumberOfReductions, &globalNumberOfReductions, 1, MPI_INT, MPI_SUM, seissol::MPI::mpi.comm());
#else
    globalNumberOfReductions = localNumberOfReductions;
#endif // USE_MPI
    totalNumberOfReductions += globalNumberOfReductions;
  } while (globalNumberOfReductions > 0);
  return totalNumberOfReductions;
}

int seissol::initializers::time_stepping::LtsWeights::enforceMaximumDifferenceLocal(PUML::TETPUML const& mesh, int* cluster, int maxDifference) {
  int numberOfReductions = 0;

  std::vector<PUML::TETPUML::cell_t> const& cells = mesh.cells();
  std::vector<PUML::TETPUML::face_t> const& faces = mesh.faces();
	int const* boundaryCond = mesh.cellData(1);

#ifdef USE_MPI
  std::unordered_map<int, std::vector<int>> rankToSharedFaces;
  std::unordered_map<int, int> localFaceIdToLocalCellId;
#endif // USE_MPI

  for (unsigned cell = 0; cell < cells.size(); ++cell) {
    int timeCluster = cluster[cell];

		unsigned int faceids[4];
		PUML::Downward::faces(mesh, cells[cell], faceids);
    for (unsigned f = 0; f < 4; ++f) {
      int difference = maxDifference;
      int boundary = getBoundaryCondition(boundaryCond, cell, f);
      // Continue for regular, dynamic rupture, and periodic boundary cells
      if (boundary == 0 || boundary == 3 || boundary == 6) {
        // We treat MPI neighbours later
        auto const& face = faces[ faceids[f] ];
        if (!face.isShared()) {
          int cellIds[2];
          PUML::Upward::cells(mesh, face, cellIds);

          int neighbourCell = (cellIds[0] == static_cast<int>(cell)) ? cellIds[1] : cellIds[0];
          int otherTimeCluster = cluster[neighbourCell];

          if (boundary == 3) {
            difference = 0;
          }

          if (timeCluster > otherTimeCluster + difference) {
            timeCluster = otherTimeCluster + difference;
            ++numberOfReductions;
          }
        }
#ifdef USE_MPI
        else {
          rankToSharedFaces[ face.shared()[0] ].push_back(faceids[f]);
          localFaceIdToLocalCellId[ faceids[f] ] = cell;
        }
#endif // USE_MPI
      }
    }
    cluster[cell] = timeCluster;
  }

#ifdef USE_MPI
  FaceSorter faceSorter(faces);
  for (auto& sharedFaces: rankToSharedFaces) {
    std::sort(sharedFaces.second.begin(), sharedFaces.second.end(), faceSorter);
  }

  auto numExchanges = rankToSharedFaces.size();
  MPI_Request* requests = new MPI_Request[2*numExchanges];
  int** ghost = new int*[numExchanges];
  int** copy = new int*[numExchanges];
  auto exchange = rankToSharedFaces.begin();
  for (unsigned ex = 0; ex < numExchanges; ++ex) {
    auto exchangeSize = exchange->second.size();
    ghost[ex] = new int[ exchangeSize ];
    copy[ex]  = new int[ exchangeSize ];

    for (unsigned n = 0; n < exchangeSize; ++n) {
      copy[ex][n] = cluster[ localFaceIdToLocalCellId[ exchange->second[n] ] ];
    }
    MPI_Isend( copy[ex], exchangeSize, MPI_INT, exchange->first, 0, seissol::MPI::mpi.comm(), &requests[ex]);
    MPI_Irecv(ghost[ex], exchangeSize, MPI_INT, exchange->first, 0, seissol::MPI::mpi.comm(), &requests[numExchanges + ex]);
    ++exchange;
  }

  MPI_Waitall(2*numExchanges, requests, MPI_STATUSES_IGNORE);

  exchange = rankToSharedFaces.begin();
  for (unsigned ex = 0; ex < numExchanges; ++ex) {
    auto exchangeSize = exchange->second.size();
    for (unsigned n = 0; n < exchangeSize; ++n) {
      int difference = maxDifference;
      int otherTimeCluster = ghost[ex][n];

      int cellIds[2];
      PUML::Upward::cells(mesh, faces[ exchange->second[n] ], cellIds);
      int cell = (cellIds[0] >= 0) ? cellIds[0] : cellIds[1];

      unsigned int faceids[4];
      PUML::Downward::faces(mesh, cells[cell], faceids);
      unsigned f = 0;
      for (; f < 4 && static_cast<int>(faceids[f]) != exchange->second[n]; ++f);
      assert(f != 4);

      int boundary = getBoundaryCondition(boundaryCond, cell, f);
      if (boundary == 3) {
        difference = 0;
      }

      if (cluster[cell] > otherTimeCluster + difference) {
        cluster[cell] = otherTimeCluster + difference;
        ++numberOfReductions;
      }
    }
    ++exchange;
  }

  for (unsigned ex = 0; ex < numExchanges; ++ex) {
    delete[] copy[ex];
    delete[] ghost[ex];
  }
  delete[] copy;
  delete[] ghost;
  delete[] requests;
#endif // USE_MPI

  return numberOfReductions;
}

max_stars_repo_name: VMPW/SeisSol

max_issues_repo_name: VMPW/SeisSol

max_forks_repo_name: VMPW/SeisSol


text:
#include <chrono>
#include <iostream>
#include <armadillo>

int main() {

	for (size_t i = 8; i < 15; i ++) {
		std::cout << "Generating matrices" << std::endl;
		size_t size = 1 << i;
		arma::fmat a = arma::randu<arma::fmat>(size, size);
		arma::fmat b = arma::randu<arma::fmat>(size, size);

		std::cout << "Starting armadillo" << std::endl;
		std::chrono::steady_clock::time_point start = std::chrono::steady_clock::now();

		arma::fmat c = a * b;

		size_t millis = std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::steady_clock::now() - start).count();
		size_t secs = millis / 1000;
		millis %= 1000;

		std::cout << size << ": " << secs << "." << millis << std::endl;
	}
}

max_stars_repo_name: vorner/fastmatmult

max_issues_repo_name: vorner/fastmatmult

max_forks_repo_name: vorner/fastmatmult


text:
/**
 *  Copyright (c) 2014, Facebook, Inc.
 *  All rights reserved.
 *
 *  This source code is licensed under the BSD-style license found in the
 *  LICENSE file in the root directory of this source tree. An additional grant
 *  of patent rights can be found in the PATENTS file in the same directory.
 */
#include "ProxyLogger.h"

#include <sys/stat.h>
#include <sys/types.h>
#include <unistd.h>

#include <string>

#include <boost/filesystem/operations.hpp>
#include <boost/filesystem/path.hpp>

#include <folly/DynamicConverter.h>
#include <folly/io/async/EventBase.h>
#include <folly/json.h>

#include "mcrouter/_router.h"
#include "mcrouter/async.h"
#include "mcrouter/lib/fbi/asox_timer.h"
#include "mcrouter/lib/fbi/cpp/util.h"
#include "mcrouter/lib/fbi/debug.h"
#include "mcrouter/proxy.h"
#include "mcrouter/stats.h"

namespace facebook { namespace memcache { namespace mcrouter {

namespace {

const char* kStatsSfx = "stats";
const char* kStatsStartupOptionsSfx = "startup_options";
const char* kConfigSourcesInfoFileName = "config_sources_info";

void log_stats(const asox_timer_t timer, void* arg) {
  auto logger = reinterpret_cast<ProxyLogger*>(arg);
  logger->log();
}

bool ensure_dir_exists_and_writeable(const std::string& path) {
  boost::system::error_code ec;
  boost::filesystem::create_directories(path, ec);
  if (ec) {
    return false;
  }

  struct stat st;
  if (::stat(path.c_str(), &st) != 0) {
    return false;
  }

  if ((st.st_mode & 0777) == 0777) {
    return true;
  }

  if (::chmod(path.c_str(), 0777) != 0) {
    return false;
  }

  return true;
}

/** returns "<source_library>.<source_service>.<source_flavor>" */
inline std::string get_stats_key(const proxy_t* proxy) {
  return folly::format(
    "libmcrouter.{}.{}",
    proxy->opts.service_name,
    proxy->opts.router_name).str();
}

std::string proxy_stats_file_path(proxy_t* proxy, const std::string& suffix) {
  boost::filesystem::path path(proxy->opts.stats_root);
  path /= get_stats_key(proxy) + "." + suffix;
  return path.string();
}

/**
 * Writes string to a file.
 */
void proxy_write_file(proxy_t* proxy,
                      const std::string& suffix,
                      const std::string& str) {
  try {
    // In case the dir was deleted some time after mcrouter started
    if (!ensure_dir_exists_and_writeable(proxy->opts.stats_root)) {
      return;
    }

    auto path = proxy_stats_file_path(proxy, suffix);

    async_write_file(proxy->stats_log_writer.get(), path, str);
  } catch (...) {
    // Do nothing
  }
}

/**
 * Determines the correct location and file name and writes the stats object
 * to disk in json format. The suffix is the file extension. If the stats root
 * directory is ever removed or is unwriteable, we just give up.
 */
void proxy_write_stats_file(proxy_t* proxy,
                            const std::string& suffix,
                            const folly::dynamic& stats) {
  auto statsString = folly::toPrettyJson(stats).toStdString() + "\n";
  proxy_write_file(proxy, suffix, statsString);
}

void write_stats_to_disk(proxy_t* proxy, const std::vector<stat_t>& stats) {
  try {
    std::string prefix = get_stats_key(proxy) + ".";
    folly::dynamic jstats = folly::dynamic::object;

    for (size_t i = 0; i < stats.size(); ++i) {
      if (stats[i].group & ods_stats) {
        auto key = prefix + stats[i].name.str();

        switch (stats[i].type) {
          case stat_uint64:
            jstats[key] = stats[i].data.uint64;
            break;

          case stat_int64:
            jstats[key] = stats[i].data.int64;
            break;

          case stat_double:
            jstats[key] = stats[i].data.dbl;
            break;

          default:
            continue;
        }
      }
    }

    proxy_write_stats_file(proxy, kStatsSfx, jstats);
  } catch (...) {
    // Do nothing
  }
}

void write_config_sources_info_to_disk(const proxy_t* proxy) {
  auto config_info_json = proxy->router->configApi->getConfigSourcesInfo();

  try {
    boost::filesystem::path path(proxy->opts.stats_root);
    path /= get_stats_key(proxy) + "." + kConfigSourcesInfoFileName;
    async_write_file(proxy->stats_log_writer.get(), path.string(),
                     folly::toPrettyJson(config_info_json).toStdString());
  } catch (...) {
    LOG(ERROR) << "Error occured while writing configuration info to disk";
  }
}

}  // anonymous namespace

ProxyLogger::ProxyLogger(proxy_t* proxy,
  std::unique_ptr<AdditionalLoggerIf> additionalLogger)
    : proxy_(proxy),
      additionalLogger_(std::move(additionalLogger)) {

  if (!ensure_dir_exists_and_writeable(proxy_->opts.stats_root)) {
    LOG(ERROR) << "Can't create or chmod " << proxy_->opts.stats_root <<
                  ", disabling stats logging";
  } else {
    auto json_options = folly::toDynamic(proxy_->router->getStartupOpts());
    json_options["pid"] = folly::to<std::string>(getpid());
    proxy_write_stats_file(proxy_, kStatsStartupOptionsSfx, json_options);
    auto path = proxy_stats_file_path(proxy_, kStatsStartupOptionsSfx);
    touchStatsFilepaths_.push_back(std::move(path));

    timeval_t delay = to<timeval_t>(proxy_->opts.stats_logging_interval);
    statsLoggingTimer_ = asox_add_timer(proxy_->eventBase->getLibeventBase(),
                                        delay, log_stats, this);
  }
}

ProxyLogger::~ProxyLogger() {
  if (statsLoggingTimer_ != nullptr) {
    asox_remove_timer(statsLoggingTimer_);
  }
}

void ProxyLogger::log() {
  std::vector<stat_t> stats(num_stats);
  try {
    proxy_->router->startupLock.wait();
    std::lock_guard<ShutdownLock> lg(proxy_->router->shutdownLock());
    std::lock_guard<std::mutex> guard(proxy_->stats_lock);

    prepare_stats(proxy_, stats.data());
  } catch (const shutdown_started_exception& e) {
    return;
  }

  for (int i = 0; i < num_stats; ++i) {
    if (stats[i].group & rate_stats) {
      stats[i].type = stat_double;
      stats[i].data.dbl = stats_rate_value(proxy_, i);
    }
  }

  write_stats_to_disk(proxy_, stats);
  write_config_sources_info_to_disk(proxy_);

  for (const auto& filepath : touchStatsFilepaths_) {
    touchFile(filepath);
  }

  if (additionalLogger_) {
    additionalLogger_->log(stats);
  }
}

}}}  // facebook::memcache::mcrouter

max_stars_repo_name: wicky-info/mcrouter

max_issues_repo_name: wicky-info/mcrouter

max_forks_repo_name: wicky-info/mcrouter


text:
//          Copyright Oliver Kowalke 2009.
// 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)

#ifndef BOOST_COROUTINES_DETAIL_PUSH_COROUTINE_IMPL_H
#define BOOST_COROUTINES_DETAIL_PUSH_COROUTINE_IMPL_H

#include <boost/assert.hpp>
#include <boost/config.hpp>
#include <boost/exception_ptr.hpp>
#include <boost/throw_exception.hpp>
#include <boost/utility.hpp>

#include <boost/coroutine/detail/config.hpp>
#include <boost/coroutine/detail/coroutine_context.hpp>
#include <boost/coroutine/detail/flags.hpp>
#include <boost/coroutine/detail/parameters.hpp>
#include <boost/coroutine/detail/trampoline_push.hpp>
#include <boost/coroutine/exceptions.hpp>

#ifdef BOOST_HAS_ABI_HEADERS
#include BOOST_ABI_PREFIX
#endif

namespace boost {
namespace coroutines {

struct stack_context;

namespace detail {

template <typename Arg> class push_coroutine_impl : private noncopyable {
protected:
  int flags_;
  exception_ptr except_;
  coroutine_context *caller_;
  coroutine_context *callee_;

public:
  typedef parameters<Arg> param_type;

  push_coroutine_impl(coroutine_context *caller, coroutine_context *callee,
                      bool unwind)
      : flags_(0), except_(), caller_(caller), callee_(callee) {
    if (unwind)
      flags_ |= flag_force_unwind;
  }

  virtual ~push_coroutine_impl() {}

  bool force_unwind() const BOOST_NOEXCEPT {
    return 0 != (flags_ & flag_force_unwind);
  }

  bool unwind_requested() const BOOST_NOEXCEPT {
    return 0 != (flags_ & flag_unwind_stack);
  }

  bool is_started() const BOOST_NOEXCEPT {
    return 0 != (flags_ & flag_started);
  }

  bool is_running() const BOOST_NOEXCEPT {
    return 0 != (flags_ & flag_running);
  }

  bool is_complete() const BOOST_NOEXCEPT {
    return 0 != (flags_ & flag_complete);
  }

  void unwind_stack() BOOST_NOEXCEPT {
    if (is_started() && !is_complete() && force_unwind()) {
      flags_ |= flag_unwind_stack;
      param_type to(unwind_t::force_unwind);
      caller_->jump(*callee_, &to);
      flags_ &= ~flag_unwind_stack;

      BOOST_ASSERT(is_complete());
    }
  }

  void push(Arg const &arg) {
    BOOST_ASSERT(!is_running());
    BOOST_ASSERT(!is_complete());

    flags_ |= flag_running;
    param_type to(const_cast<Arg *>(&arg), this);
    param_type *from(static_cast<param_type *>(caller_->jump(*callee_, &to)));
    flags_ &= ~flag_running;
    if (from->do_unwind)
      throw forced_unwind();
    if (except_)
      rethrow_exception(except_);
  }

  void push(BOOST_RV_REF(Arg) arg) {
    BOOST_ASSERT(!is_running());
    BOOST_ASSERT(!is_complete());

    flags_ |= flag_running;
    param_type to(const_cast<Arg *>(&arg), this);
    param_type *from(static_cast<param_type *>(caller_->jump(*callee_, &to)));
    flags_ &= ~flag_running;
    if (from->do_unwind)
      throw forced_unwind();
    if (except_)
      rethrow_exception(except_);
  }

  virtual void destroy() = 0;
};

template <typename Arg> class push_coroutine_impl<Arg &> : private noncopyable {
protected:
  int flags_;
  exception_ptr except_;
  coroutine_context *caller_;
  coroutine_context *callee_;

public:
  typedef parameters<Arg &> param_type;

  push_coroutine_impl(coroutine_context *caller, coroutine_context *callee,
                      bool unwind)
      : flags_(0), except_(), caller_(caller), callee_(callee) {
    if (unwind)
      flags_ |= flag_force_unwind;
  }

  virtual ~push_coroutine_impl() {}

  bool force_unwind() const BOOST_NOEXCEPT {
    return 0 != (flags_ & flag_force_unwind);
  }

  bool unwind_requested() const BOOST_NOEXCEPT {
    return 0 != (flags_ & flag_unwind_stack);
  }

  bool is_started() const BOOST_NOEXCEPT {
    return 0 != (flags_ & flag_started);
  }

  bool is_running() const BOOST_NOEXCEPT {
    return 0 != (flags_ & flag_running);
  }

  bool is_complete() const BOOST_NOEXCEPT {
    return 0 != (flags_ & flag_complete);
  }

  void unwind_stack() BOOST_NOEXCEPT {
    if (is_started() && !is_complete() && force_unwind()) {
      flags_ |= flag_unwind_stack;
      param_type to(unwind_t::force_unwind);
      caller_->jump(*callee_, &to);
      flags_ &= ~flag_unwind_stack;

      BOOST_ASSERT(is_complete());
    }
  }

  void push(Arg &arg) {
    BOOST_ASSERT(!is_running());
    BOOST_ASSERT(!is_complete());

    flags_ |= flag_running;
    param_type to(&arg, this);
    param_type *from(static_cast<param_type *>(caller_->jump(*callee_, &to)));
    flags_ &= ~flag_running;
    if (from->do_unwind)
      throw forced_unwind();
    if (except_)
      rethrow_exception(except_);
  }

  virtual void destroy() = 0;
};

template <> class push_coroutine_impl<void> : private noncopyable {
protected:
  int flags_;
  exception_ptr except_;
  coroutine_context *caller_;
  coroutine_context *callee_;

public:
  typedef parameters<void> param_type;

  push_coroutine_impl(coroutine_context *caller, coroutine_context *callee,
                      bool unwind)
      : flags_(0), except_(), caller_(caller), callee_(callee) {
    if (unwind)
      flags_ |= flag_force_unwind;
  }

  virtual ~push_coroutine_impl() {}

  inline bool force_unwind() const BOOST_NOEXCEPT {
    return 0 != (flags_ & flag_force_unwind);
  }

  inline bool unwind_requested() const BOOST_NOEXCEPT {
    return 0 != (flags_ & flag_unwind_stack);
  }

  inline bool is_started() const BOOST_NOEXCEPT {
    return 0 != (flags_ & flag_started);
  }

  inline bool is_running() const BOOST_NOEXCEPT {
    return 0 != (flags_ & flag_running);
  }

  inline bool is_complete() const BOOST_NOEXCEPT {
    return 0 != (flags_ & flag_complete);
  }

  inline void unwind_stack() BOOST_NOEXCEPT {
    if (is_started() && !is_complete() && force_unwind()) {
      flags_ |= flag_unwind_stack;
      param_type to(unwind_t::force_unwind);
      caller_->jump(*callee_, &to);
      flags_ &= ~flag_unwind_stack;

      BOOST_ASSERT(is_complete());
    }
  }

  inline void push() {
    BOOST_ASSERT(!is_running());
    BOOST_ASSERT(!is_complete());

    flags_ |= flag_running;
    param_type to(this);
    param_type *from(static_cast<param_type *>(caller_->jump(*callee_, &to)));
    flags_ &= ~flag_running;
    if (from->do_unwind)
      throw forced_unwind();
    if (except_)
      rethrow_exception(except_);
  }

  virtual void destroy() = 0;
};

} // namespace detail
} // namespace coroutines
} // namespace boost

#ifdef BOOST_HAS_ABI_HEADERS
#include BOOST_ABI_SUFFIX
#endif

#endif // BOOST_COROUTINES_DETAIL_PUSH_COROUTINE_IMPL_H

max_stars_repo_name: henrywarhurst/matrix

max_issues_repo_name: henrywarhurst/matrix

max_forks_repo_name: henrywarhurst/matrix


text:
#include "CommitHistory.h"

#include <iostream>

#include <boost/date_time/posix_time/posix_time.hpp>
#include <log4cxx/logger.h>
#include <yaml-cpp/yaml.h>
#include <yaml-cpp/exceptions.h>

#include "LocalYAMLStorage.h"

namespace
{
    log4cxx::LoggerPtr spLogger(log4cxx::Logger::getLogger("CommitHistory"));
}

boost::optional<std::string> YAMLSerializer(const CommitHistory& commit)
{
    YAML::Node commitNode;
    commitNode["id"] = commit.mId;
    commitNode["parent-id"] = commit.mParent;
    commitNode["time"] = boost::posix_time::to_simple_string(commit.mTime);

    YAML::Node deltas;
    for (const auto& item : commit.mItems)
    {
        YAML::Node delta;
        delta["key"] = item.first.GetKey().ToString();
        YAML::Node old = UnpackIntoYAML(item.first);
        YAML::Node current = UnpackIntoYAML(item.second);
        delta["old"] = old;
        delta["current"] = current;
        deltas.push_back(delta);
    }

    commitNode["deltas"] = deltas;

    std::stringstream out;
    out << commitNode << std::endl;
    return boost::optional<std::string>(out.str());
}

boost::optional<CommitHistory> YAMLDeserializer(const std::string& buffer)
{
    YAML::Node commitNode = YAML::Load(buffer);
    CommitHistory commit;

    commit.mId = commitNode["id"].as<CommitHistory::CommitId>();
    commit.mParent = commitNode["parent-id"].as<CommitHistory::CommitId>();

    if (commitNode["time"])
    {
        std::string timeString = commitNode["time"].as<std::string>();
        try
        {
            commit.mTime = boost::posix_time::time_from_string(timeString);
        }
        catch(const boost::bad_lexical_cast&)
        {
            //This is the exception boost throws when it cannot convert a string to a posix_time::ptime
            commit.mTime = boost::posix_time::ptime(boost::posix_time::not_a_date_time);
        }
    }
    else
    {
        commit.mTime = boost::posix_time::ptime(boost::posix_time::not_a_date_time);
    }

    for (const auto& item : commitNode["deltas"])
    {
        std::string key = item["key"].as<std::string>();
        std::string oldtype = item["old"]["type"].as<std::string>();
        ConfigValue old = (oldtype != "NULL") ? ValueFromYAML(oldtype, item["old"]["value"]) : ConfigValue();
        std::string currenttype = item["current"]["type"].as<std::string>();
        ConfigValue current = (currenttype != "NULL") ? ValueFromYAML(currenttype, item["current"]["value"]) : ConfigValue();
        commit.mItems.emplace_back(ConfigPair(key, old), ConfigPair(key, current));
    }

    return boost::optional<CommitHistory>(commit);
}

max_stars_repo_name: kcollinssibley/sysadmin

max_issues_repo_name: kcollinssibley/sysadmin

max_forks_repo_name: kcollinssibley/sysadmin


text:
#include <boost/python/module.hpp>
#include <boost/python/class.hpp>
#include <boost/python/args.hpp>

#include <scitbx/wigner/wigner3j.h>

namespace scitbx { namespace wigner {
namespace boost_python{

  struct wigner3j_fast_wrapper
  {
    typedef wigner3j_fast < double > w_t;
    static void
    wrap()
    {
      using namespace boost::python;
      class_<w_t>("wigner3j_fast", no_init)
        .def( init< int const&
                  >
             (( arg("max")
             ))
            )
        .def("compute", &w_t::compute)
      ;
    }
  };

  struct wigner3j_wrapper
  {
    typedef wigner3j < double > w_t;
    static void
    wrap()
    {
      using namespace boost::python;
      class_<w_t>("wigner3j", no_init)
        .def( init< int const&,
                    int const&,
                    int const&,
                    int const&,
                    int const&,
                    int const&
                  >
             (( arg("j1"),
                arg("j2"),
                arg("j3"),
                arg("m1"),
                arg("m2"),
                arg("m3")
             ))
            )
        .def("check", &w_t::check)
        .def("get_value", &w_t::get_value)
      ;
    }
  };

  void
  init_module()
  {
    wigner3j_fast_wrapper::wrap();
    wigner3j_wrapper::wrap();
  }

}}} // scitbx::wigner::boost_python

BOOST_PYTHON_MODULE(scitbx_wigner_ext)
{
  scitbx::wigner::boost_python::init_module();
}


text:
/*=============================================================================
  Copyright (c) 2010-2016 Bolero MURAKAMI
  https://github.com/bolero-MURAKAMI/Sprig

  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)
=============================================================================*/
#ifndef SPRIG_TMP_ITER_AT_HPP
#define SPRIG_TMP_ITER_AT_HPP

#include <sprig/config/config.hpp>

#ifdef SPRIG_USING_PRAGMA_ONCE
#	pragma once
#endif	// #ifdef SPRIG_USING_PRAGMA_ONCE

#include <boost/mpl/advance.hpp>
#include <boost/mpl/begin.hpp>

namespace sprig {
	namespace tmp {
		//
		// iter_at
		//
		template<typename Sequence, typename indexT>
		struct iter_at
			: public boost::mpl::advance<
				typename boost::mpl::begin<Sequence>::type,
				indexT
			>
		{};
	}	// namespace tmp
}	// namespace sprig

#endif	// #ifndef SPRIG_TMP_ITER_AT_HPP

max_stars_repo_name: bolero-MURAKAMI/Sprig

max_issues_repo_name: bolero-MURAKAMI/Sprig

max_forks_repo_name: bolero-MURAKAMI/Sprig


text:
/**
 *    Copyright (C) 2020-present MongoDB, Inc.
 *
 *    This program is free software: you can redistribute it and/or modify
 *    it under the terms of the Server Side Public License, version 1,
 *    as published by MongoDB, Inc.
 *
 *    This program is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    Server Side Public License for more details.
 *
 *    You should have received a copy of the Server Side Public License
 *    along with this program. If not, see
 *    <http://www.mongodb.com/licensing/server-side-public-license>.
 *
 *    As a special exception, the copyright holders give permission to link the
 *    code of portions of this program with the OpenSSL library under certain
 *    conditions as described in each individual source file and distribute
 *    linked combinations including the program with the OpenSSL library. You
 *    must comply with the Server Side Public License in all respects for
 *    all of the code used other than as permitted herein. If you modify file(s)
 *    with this exception, you may extend this exception to your version of the
 *    file(s), but you are not obligated to do so. If you do not wish to do so,
 *    delete this exception statement from your version. If you delete this
 *    exception statement from all source files in the program, then also delete
 *    it in the license file.
 */

#include "mongo/platform/basic.h"

#include <boost/intrusive_ptr.hpp>
#include <string>

#include "mongo/bson/bsonobj.h"
#include "mongo/bson/unordered_fields_bsonobj_comparator.h"
#include "mongo/db/cst/c_node.h"
#include "mongo/db/cst/cst_pipeline_translation.h"
#include "mongo/db/cst/key_fieldname.h"
#include "mongo/db/cst/key_value.h"
#include "mongo/db/exec/document_value/document.h"
#include "mongo/db/namespace_string.h"
#include "mongo/db/pipeline/document_source_limit.h"
#include "mongo/db/pipeline/document_source_single_document_transformation.h"
#include "mongo/db/pipeline/document_source_skip.h"
#include "mongo/db/pipeline/expression_context_for_test.h"
#include "mongo/db/query/util/make_data_structure.h"
#include "mongo/unittest/unittest.h"

namespace mongo {
namespace {
using namespace std::string_literals;

auto getExpCtx() {
    auto nss = NamespaceString{"db", "coll"};
    return boost::intrusive_ptr<ExpressionContextForTest>{new ExpressionContextForTest(nss)};
}

auto makePipelineContainingProjectStageWithLiteral(CNode&& literal) {
    return CNode{CNode::ArrayChildren{CNode{CNode::ObjectChildren{
        {KeyFieldname::projectInclusion,
         CNode{CNode::ObjectChildren{
             {ProjectionPath{makeVector<std::string>("a")},
              CNode{CNode::ObjectChildren{{KeyFieldname::literal, std::move(literal)}}}}}}}}}}};
}

TEST(CstLiteralsTest, TranslatesDouble) {
    auto cst = makePipelineContainingProjectStageWithLiteral(CNode{UserDouble{5e-324}});
    // DocumenSourceSingleDoucmentTransformation reorders fields so we need to be insensitive.
    ASSERT(UnorderedFieldsBSONObjComparator{}.evaluate(
        BSON("_id" << true << "a" << BSON("$const" << 5e-324)) ==
        dynamic_cast<DocumentSourceSingleDocumentTransformation&>(
            **cst_pipeline_translation::translatePipeline(cst, getExpCtx())->getSources().begin())
            .getTransformer()
            .serializeTransformation(boost::none)
            .toBson()));
}

TEST(CstLiteralsTest, TranslatesString) {
    auto cst = makePipelineContainingProjectStageWithLiteral(CNode{UserString{"soup can"}});
    // DocumenSourceSingleDoucmentTransformation reorders fields so we need to be insensitive.
    ASSERT(UnorderedFieldsBSONObjComparator{}.evaluate(
        BSON("_id" << true << "a"
                   << BSON("$const"
                           << "soup can")) ==
        dynamic_cast<DocumentSourceSingleDocumentTransformation&>(
            **cst_pipeline_translation::translatePipeline(cst, getExpCtx())->getSources().begin())
            .getTransformer()
            .serializeTransformation(boost::none)
            .toBson()));
}

TEST(CstLiteralsTest, TranslatesBinary) {
    auto cst =
        makePipelineContainingProjectStageWithLiteral(CNode{UserBinary{"a\0b", 3, BinDataGeneral}});
    // DocumenSourceSingleDoucmentTransformation reorders fields so we need to be insensitive.
    ASSERT(UnorderedFieldsBSONObjComparator{}.evaluate(
        BSON("_id" << true << "a" << BSON("$const" << BSONBinData("a\0b", 3, BinDataGeneral))) ==
        dynamic_cast<DocumentSourceSingleDocumentTransformation&>(
            **cst_pipeline_translation::translatePipeline(cst, getExpCtx())->getSources().begin())
            .getTransformer()
            .serializeTransformation(boost::none)
            .toBson()));
}

TEST(CstLiteralsTest, TranslatesUndefined) {
    auto cst = makePipelineContainingProjectStageWithLiteral(CNode{UserUndefined{}});
    // DocumenSourceSingleDoucmentTransformation reorders fields so we need to be insensitive.
    ASSERT(UnorderedFieldsBSONObjComparator{}.evaluate(
        BSON("_id" << true << "a" << BSON("$const" << BSONUndefined)) ==
        dynamic_cast<DocumentSourceSingleDocumentTransformation&>(
            **cst_pipeline_translation::translatePipeline(cst, getExpCtx())->getSources().begin())
            .getTransformer()
            .serializeTransformation(boost::none)
            .toBson()));
}

TEST(CstLiteralsTest, TranslatesObjectId) {
    auto cst = makePipelineContainingProjectStageWithLiteral(
        CNode{UserObjectId{"01234567890123456789aaaa"}});
    // DocumenSourceSingleDoucmentTransformation reorders fields so we need to be insensitive.
    ASSERT(UnorderedFieldsBSONObjComparator{}.evaluate(
        BSON("_id" << true << "a" << BSON("$const" << OID("01234567890123456789aaaa"))) ==
        dynamic_cast<DocumentSourceSingleDocumentTransformation&>(
            **cst_pipeline_translation::translatePipeline(cst, getExpCtx())->getSources().begin())
            .getTransformer()
            .serializeTransformation(boost::none)
            .toBson()));
}

TEST(CstLiteralsTest, TranslatesBoolean) {
    auto cst = makePipelineContainingProjectStageWithLiteral(CNode{UserBoolean{false}});
    // DocumenSourceSingleDoucmentTransformation reorders fields so we need to be insensitive.
    ASSERT(UnorderedFieldsBSONObjComparator{}.evaluate(
        BSON("_id" << true << "a" << BSON("$const" << false)) ==
        dynamic_cast<DocumentSourceSingleDocumentTransformation&>(
            **cst_pipeline_translation::translatePipeline(cst, getExpCtx())->getSources().begin())
            .getTransformer()
            .serializeTransformation(boost::none)
            .toBson()));
}

TEST(CstLiteralsTest, TranslatesDate) {
    auto cst = makePipelineContainingProjectStageWithLiteral(
        CNode{UserDate{Date_t::fromMillisSinceEpoch(424242)}});
    // DocumenSourceSingleDoucmentTransformation reorders fields so we need to be insensitive.
    ASSERT(UnorderedFieldsBSONObjComparator{}.evaluate(
        BSON("_id" << true << "a" << BSON("$const" << Date_t::fromMillisSinceEpoch(424242))) ==
        dynamic_cast<DocumentSourceSingleDocumentTransformation&>(
            **cst_pipeline_translation::translatePipeline(cst, getExpCtx())->getSources().begin())
            .getTransformer()
            .serializeTransformation(boost::none)
            .toBson()));
}

TEST(CstLiteralsTest, TranslatesNull) {
    auto cst = makePipelineContainingProjectStageWithLiteral(CNode{UserNull{}});
    // DocumenSourceSingleDoucmentTransformation reorders fields so we need to be insensitive.
    ASSERT(UnorderedFieldsBSONObjComparator{}.evaluate(
        BSON("_id" << true << "a" << BSON("$const" << BSONNULL)) ==
        dynamic_cast<DocumentSourceSingleDocumentTransformation&>(
            **cst_pipeline_translation::translatePipeline(cst, getExpCtx())->getSources().begin())
            .getTransformer()
            .serializeTransformation(boost::none)
            .toBson()));
}

TEST(CstLiteralsTest, TranslatesRegex) {
    auto cst = makePipelineContainingProjectStageWithLiteral(CNode{UserRegex{".*", "i"}});
    // DocumenSourceSingleDoucmentTransformation reorders fields so we need to be insensitive.
    ASSERT(UnorderedFieldsBSONObjComparator{}.evaluate(
        BSON("_id" << true << "a" << BSON("$const" << BSONRegEx(".*", "i"))) ==
        dynamic_cast<DocumentSourceSingleDocumentTransformation&>(
            **cst_pipeline_translation::translatePipeline(cst, getExpCtx())->getSources().begin())
            .getTransformer()
            .serializeTransformation(boost::none)
            .toBson()));
}

TEST(CstLiteralsTest, TranslatesDBPointer) {
    auto cst = makePipelineContainingProjectStageWithLiteral(
        CNode{UserDBPointer{"db.c", OID("010203040506070809101112")}});
    // DocumenSourceSingleDoucmentTransformation reorders fields so we need to be insensitive.
    ASSERT(UnorderedFieldsBSONObjComparator{}.evaluate(
        BSON("_id" << true << "a"
                   << BSON("$const" << BSONDBRef("db.c", OID("010203040506070809101112")))) ==
        dynamic_cast<DocumentSourceSingleDocumentTransformation&>(
            **cst_pipeline_translation::translatePipeline(cst, getExpCtx())->getSources().begin())
            .getTransformer()
            .serializeTransformation(boost::none)
            .toBson()));
}

TEST(CstLiteralsTest, TranslatesJavascript) {
    auto cst = makePipelineContainingProjectStageWithLiteral(CNode{UserJavascript{"5 === 5"}});
    // DocumenSourceSingleDoucmentTransformation reorders fields so we need to be insensitive.
    ASSERT(UnorderedFieldsBSONObjComparator{}.evaluate(
        BSON("_id" << true << "a" << BSON("$const" << BSONCode("5 === 5"))) ==
        dynamic_cast<DocumentSourceSingleDocumentTransformation&>(
            **cst_pipeline_translation::translatePipeline(cst, getExpCtx())->getSources().begin())
            .getTransformer()
            .serializeTransformation(boost::none)
            .toBson()));
}

TEST(CstLiteralsTest, TranslatesSymbol) {
    auto cst = makePipelineContainingProjectStageWithLiteral(CNode{UserSymbol{"foo"}});
    // DocumenSourceSingleDoucmentTransformation reorders fields so we need to be insensitive.
    ASSERT(UnorderedFieldsBSONObjComparator{}.evaluate(
        BSON("_id" << true << "a" << BSON("$const" << BSONSymbol("foo"))) ==
        dynamic_cast<DocumentSourceSingleDocumentTransformation&>(
            **cst_pipeline_translation::translatePipeline(cst, getExpCtx())->getSources().begin())
            .getTransformer()
            .serializeTransformation(boost::none)
            .toBson()));
}

TEST(CstLiteralsTest, TranslatesJavascriptWithScope) {
    auto cst = makePipelineContainingProjectStageWithLiteral(
        CNode{UserJavascriptWithScope{"6 === 6", BSONObj{}}});
    // DocumenSourceSingleDoucmentTransformation reorders fields so we need to be insensitive.
    ASSERT(UnorderedFieldsBSONObjComparator{}.evaluate(
        BSON("_id" << true << "a" << BSON("$const" << BSONCodeWScope("6 === 6", BSONObj()))) ==
        dynamic_cast<DocumentSourceSingleDocumentTransformation&>(
            **cst_pipeline_translation::translatePipeline(cst, getExpCtx())->getSources().begin())
            .getTransformer()
            .serializeTransformation(boost::none)
            .toBson()));
}

TEST(CstLiteralsTest, TranslatesInt) {
    auto cst = makePipelineContainingProjectStageWithLiteral(CNode{UserInt{777}});
    // DocumenSourceSingleDoucmentTransformation reorders fields so we need to be insensitive.
    ASSERT(UnorderedFieldsBSONObjComparator{}.evaluate(
        BSON("_id" << true << "a" << BSON("$const" << 777)) ==
        dynamic_cast<DocumentSourceSingleDocumentTransformation&>(
            **cst_pipeline_translation::translatePipeline(cst, getExpCtx())->getSources().begin())
            .getTransformer()
            .serializeTransformation(boost::none)
            .toBson()));
}

TEST(CstLiteralsTest, TranslatesTimestamp) {
    auto cst = makePipelineContainingProjectStageWithLiteral(CNode{UserTimestamp{4102444800, 1}});
    // DocumenSourceSingleDoucmentTransformation reorders fields so we need to be insensitive.
    ASSERT(UnorderedFieldsBSONObjComparator{}.evaluate(
        BSON("_id" << true << "a" << BSON("$const" << Timestamp(4102444800, 1))) ==
        dynamic_cast<DocumentSourceSingleDocumentTransformation&>(
            **cst_pipeline_translation::translatePipeline(cst, getExpCtx())->getSources().begin())
            .getTransformer()
            .serializeTransformation(boost::none)
            .toBson()));
}

TEST(CstLiteralsTest, TranslatesLong) {
    auto cst = makePipelineContainingProjectStageWithLiteral(CNode{UserLong{777777777777777777ll}});
    // DocumenSourceSingleDoucmentTransformation reorders fields so we need to be insensitive.
    ASSERT(UnorderedFieldsBSONObjComparator{}.evaluate(
        BSON("_id" << true << "a" << BSON("$const" << 777777777777777777ll)) ==
        dynamic_cast<DocumentSourceSingleDocumentTransformation&>(
            **cst_pipeline_translation::translatePipeline(cst, getExpCtx())->getSources().begin())
            .getTransformer()
            .serializeTransformation(boost::none)
            .toBson()));
}

TEST(CstLiteralsTest, TranslatesDecimal) {
    auto cst = makePipelineContainingProjectStageWithLiteral(
        CNode{UserDecimal{Decimal128::kLargestNegative}});
    // DocumenSourceSingleDoucmentTransformation reorders fields so we need to be insensitive.
    ASSERT(UnorderedFieldsBSONObjComparator{}.evaluate(
        BSON("_id" << true << "a" << BSON("$const" << Decimal128::kLargestNegative)) ==
        dynamic_cast<DocumentSourceSingleDocumentTransformation&>(
            **cst_pipeline_translation::translatePipeline(cst, getExpCtx())->getSources().begin())
            .getTransformer()
            .serializeTransformation(boost::none)
            .toBson()));
}

TEST(CstLiteralsTest, TranslatesMinKey) {
    auto cst = makePipelineContainingProjectStageWithLiteral(CNode{UserMinKey{}});
    // DocumenSourceSingleDoucmentTransformation reorders fields so we need to be insensitive.
    ASSERT(UnorderedFieldsBSONObjComparator{}.evaluate(
        BSON("_id" << true << "a" << BSON("$const" << MINKEY)) ==
        dynamic_cast<DocumentSourceSingleDocumentTransformation&>(
            **cst_pipeline_translation::translatePipeline(cst, getExpCtx())->getSources().begin())
            .getTransformer()
            .serializeTransformation(boost::none)
            .toBson()));
}

TEST(CstLiteralsTest, TranslatesMaxKey) {
    auto cst = makePipelineContainingProjectStageWithLiteral(CNode{UserMaxKey{}});
    // DocumenSourceSingleDoucmentTransformation reorders fields so we need to be insensitive.
    ASSERT(UnorderedFieldsBSONObjComparator{}.evaluate(
        BSON("_id" << true << "a" << BSON("$const" << MAXKEY)) ==
        dynamic_cast<DocumentSourceSingleDocumentTransformation&>(
            **cst_pipeline_translation::translatePipeline(cst, getExpCtx())->getSources().begin())
            .getTransformer()
            .serializeTransformation(boost::none)
            .toBson()));
}

TEST(CstLiteralsTest, TranslatesArray) {
    auto cst = makePipelineContainingProjectStageWithLiteral(CNode{});
    // DocumenSourceSingleDoucmentTransformation reorders fields so we need to be insensitive.
    ASSERT(UnorderedFieldsBSONObjComparator{}.evaluate(
        BSON("_id" << true << "a" << BSON("$const" << BSONArray())) ==
        dynamic_cast<DocumentSourceSingleDocumentTransformation&>(
            **cst_pipeline_translation::translatePipeline(cst, getExpCtx())->getSources().begin())
            .getTransformer()
            .serializeTransformation(boost::none)
            .toBson()));
}

TEST(CstLiteralsTest, TranslatesObject) {
    auto cst = makePipelineContainingProjectStageWithLiteral(CNode{CNode::ObjectChildren{}});
    // DocumenSourceSingleDoucmentTransformation reorders fields so we need to be insensitive.
    ASSERT(UnorderedFieldsBSONObjComparator{}.evaluate(
        BSON("_id" << true << "a" << BSON("$const" << BSONObj())) ==
        dynamic_cast<DocumentSourceSingleDocumentTransformation&>(
            **cst_pipeline_translation::translatePipeline(cst, getExpCtx())->getSources().begin())
            .getTransformer()
            .serializeTransformation(boost::none)
            .toBson()));
}

TEST(CstLiteralsTest, TranslatesNestedLiteral) {
    auto cst = makePipelineContainingProjectStageWithLiteral(CNode{CNode::ObjectChildren{
        {UserFieldname{"a"}, CNode{UserMaxKey{}}},
        {UserFieldname{"b"},
         CNode{CNode::ObjectChildren{{UserFieldname{"1"}, CNode{UserDecimal{1.0}}},
                                     {UserFieldname{"2"}, CNode{UserLong{2ll}}}}}},
        {UserFieldname{"c"},
         CNode{CNode::ArrayChildren{CNode{UserString{"foo"}}, CNode{UserSymbol{"bar"}}}}}}});
    // DocumenSourceSingleDoucmentTransformation reorders fields so we need to be insensitive.
    ASSERT(UnorderedFieldsBSONObjComparator{}.evaluate(
        BSON("_id" << true << "a"
                   << BSON("$const" << BSON("a" << MAXKEY << "b"
                                                << BSON("1" << Decimal128(1.0) << "2" << 2ll) << "c"
                                                << BSON_ARRAY("foo" << BSONSymbol("bar"))))) ==
        dynamic_cast<DocumentSourceSingleDocumentTransformation&>(
            **cst_pipeline_translation::translatePipeline(cst, getExpCtx())->getSources().begin())
            .getTransformer()
            .serializeTransformation(boost::none)
            .toBson()));
}

}  // namespace
}  // namespace mongo

max_stars_repo_name: benety/mongo

max_issues_repo_name: benety/mongo

max_forks_repo_name: benety/mongo


text:
/*
 * Software License Agreement (BSD License)
 *
 *  Point Cloud Library (PCL) - www.pointclouds.org
 *  Copyright (c) 2010-2011, Willow Garage, Inc.
 *
 *  All rights reserved.
 *
 *  Redistribution and use in source and binary forms, with or without
 *  modification, are permitted provided that the following conditions
 *  are met:
 *
 *   * Redistributions of source code must retain the above copyright
 *     notice, this list of conditions and the following disclaimer.
 *   * Redistributions in binary form must reproduce the above
 *     copyright notice, this list of conditions and the following
 *     disclaimer in the documentation and/or other materials provided
 *     with the distribution.
 *   * Neither the name of the copyright holder(s) nor the names of its
 *     contributors may be used to endorse or promote products derived
 *     from this software without specific prior written permission.
 *
 *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 *  "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 *  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
 *  FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
 *  COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
 *  INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
 *  BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 *  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 *  CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 *  LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
 *  ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 *  POSSIBILITY OF SUCH DAMAGE.
 *
 */

#ifndef PCL_VOXEL_GRID_COVARIANCE_IMPL_OMP_H_
#define PCL_VOXEL_GRID_COVARIANCE_IMPL_OMP_H_

#include "voxel_grid_covariance_omp.h"
#include <Eigen/Cholesky>
#include <Eigen/Dense>
#include <pcl/common/common.h>
#include <pcl/filters/boost.h>

////////////////////////////////////////////////////////////////////////////////
template <typename PointT>
void pclomp::VoxelGridCovariance<PointT>::applyFilter(PointCloud &output) {
  voxel_centroids_leaf_indices_.clear();

  // Has the input dataset been set already?
  if (!input_) {
    PCL_WARN("[pcl::%s::applyFilter] No input dataset given!\n",
             getClassName().c_str());
    output.width = output.height = 0;
    output.points.clear();
    return;
  }

  // Copy the header (and thus the frame_id) + allocate enough space for points
  output.height = 1;      // downsampling breaks the organized structure
  output.is_dense = true; // we filter out invalid points
  output.points.clear();

  Eigen::Vector4f min_p, max_p;
  // Get the minimum and maximum dimensions
  if (!filter_field_name_
           .empty()) // If we don't want to process the entire cloud...
    pcl::getMinMax3D<PointT>(input_, filter_field_name_,
                             static_cast<float>(filter_limit_min_),
                             static_cast<float>(filter_limit_max_), min_p,
                             max_p, filter_limit_negative_);
  else
    pcl::getMinMax3D<PointT>(*input_, min_p, max_p);

  // Check that the leaf size is not too small, given the size of the data
  int64_t dx =
      static_cast<int64_t>((max_p[0] - min_p[0]) * inverse_leaf_size_[0]) + 1;
  int64_t dy =
      static_cast<int64_t>((max_p[1] - min_p[1]) * inverse_leaf_size_[1]) + 1;
  int64_t dz =
      static_cast<int64_t>((max_p[2] - min_p[2]) * inverse_leaf_size_[2]) + 1;

  if ((dx * dy * dz) > std::numeric_limits<int32_t>::max()) {
    PCL_WARN("[pcl::%s::applyFilter] Leaf size is too small for the input "
             "dataset. Integer indices would overflow.",
             getClassName().c_str());
    output.clear();
    return;
  }

  // Compute the minimum and maximum bounding box values
  min_b_[0] = static_cast<int>(floor(min_p[0] * inverse_leaf_size_[0]));
  max_b_[0] = static_cast<int>(floor(max_p[0] * inverse_leaf_size_[0]));
  min_b_[1] = static_cast<int>(floor(min_p[1] * inverse_leaf_size_[1]));
  max_b_[1] = static_cast<int>(floor(max_p[1] * inverse_leaf_size_[1]));
  min_b_[2] = static_cast<int>(floor(min_p[2] * inverse_leaf_size_[2]));
  max_b_[2] = static_cast<int>(floor(max_p[2] * inverse_leaf_size_[2]));

  // Compute the number of divisions needed along all axis
  div_b_ = max_b_ - min_b_ + Eigen::Vector4i::Ones();
  div_b_[3] = 0;

  // Clear the leaves
  leaves_.clear();
  //  leaves_.reserve(8192);

  // Set up the division multiplier
  divb_mul_ = Eigen::Vector4i(1, div_b_[0], div_b_[0] * div_b_[1], 0);

  int centroid_size = 4;

  if (downsample_all_data_)
    centroid_size = boost::mpl::size<FieldList>::value;

  // ---[ RGB special case
  std::vector<pcl::PCLPointField> fields;
  int rgba_index = -1;
  rgba_index = pcl::getFieldIndex(*input_, "rgb", fields);
  if (rgba_index == -1)
    rgba_index = pcl::getFieldIndex(*input_, "rgba", fields);
  if (rgba_index >= 0) {
    rgba_index = fields[rgba_index].offset;
    centroid_size += 3;
  }

  // If we don't want to process the entire cloud, but rather filter points far
  // away from the viewpoint first...
  if (!filter_field_name_.empty()) {
    // Get the distance field index
    std::vector<pcl::PCLPointField> fields;
    int distance_idx = pcl::getFieldIndex(*input_, filter_field_name_, fields);
    if (distance_idx == -1)
      PCL_WARN(
          "[pcl::%s::applyFilter] Invalid filter field name. Index is %d.\n",
          getClassName().c_str(), distance_idx);

    // First pass: go over all points and insert them into the right leaf
    for (size_t cp = 0; cp < input_->points.size(); ++cp) {
      if (!input_->is_dense)
        // Check if the point is invalid
        if (!pcl_isfinite(input_->points[cp].x) ||
            !pcl_isfinite(input_->points[cp].y) ||
            !pcl_isfinite(input_->points[cp].z))
          continue;

      // Get the distance value
      const uint8_t *pt_data =
          reinterpret_cast<const uint8_t *>(&input_->points[cp]);
      float distance_value = 0;
      memcpy(&distance_value, pt_data + fields[distance_idx].offset,
             sizeof(float));

      if (filter_limit_negative_) {
        // Use a threshold for cutting out points which inside the interval
        if ((distance_value < filter_limit_max_) &&
            (distance_value > filter_limit_min_))
          continue;
      } else {
        // Use a threshold for cutting out points which are too close/far away
        if ((distance_value > filter_limit_max_) ||
            (distance_value < filter_limit_min_))
          continue;
      }

      int ijk0 =
          static_cast<int>(floor(input_->points[cp].x * inverse_leaf_size_[0]) -
                           static_cast<float>(min_b_[0]));
      int ijk1 =
          static_cast<int>(floor(input_->points[cp].y * inverse_leaf_size_[1]) -
                           static_cast<float>(min_b_[1]));
      int ijk2 =
          static_cast<int>(floor(input_->points[cp].z * inverse_leaf_size_[2]) -
                           static_cast<float>(min_b_[2]));

      // Compute the centroid leaf index
      int idx = ijk0 * divb_mul_[0] + ijk1 * divb_mul_[1] + ijk2 * divb_mul_[2];

      Leaf &leaf = leaves_[idx];
      if (leaf.nr_points == 0) {
        leaf.centroid.resize(centroid_size);
        leaf.centroid.setZero();
      }

      Eigen::Vector3d pt3d(input_->points[cp].x, input_->points[cp].y,
                           input_->points[cp].z);
      // Accumulate point sum for centroid calculation
      leaf.mean_ += pt3d;
      // Accumulate x*xT for single pass covariance calculation
      leaf.cov_ += pt3d * pt3d.transpose();

      // Do we need to process all the fields?
      if (!downsample_all_data_) {
        Eigen::Vector4f pt(input_->points[cp].x, input_->points[cp].y,
                           input_->points[cp].z, 0);
        leaf.centroid.template head<4>() += pt;
      } else {
        // Copy all the fields
        Eigen::VectorXf centroid = Eigen::VectorXf::Zero(centroid_size);
        // ---[ RGB special case
        if (rgba_index >= 0) {
          // fill r/g/b data
          int rgb;
          memcpy(&rgb,
                 reinterpret_cast<const char *>(&input_->points[cp]) +
                     rgba_index,
                 sizeof(int));
          centroid[centroid_size - 3] =
              static_cast<float>((rgb >> 16) & 0x0000ff);
          centroid[centroid_size - 2] =
              static_cast<float>((rgb >> 8) & 0x0000ff);
          centroid[centroid_size - 1] = static_cast<float>((rgb)&0x0000ff);
        }
        pcl::for_each_type<FieldList>(
            pcl::NdCopyPointEigenFunctor<PointT>(input_->points[cp], centroid));
        leaf.centroid += centroid;
      }
      ++leaf.nr_points;
    }
  }
  // No distance filtering, process all data
  else {
    // First pass: go over all points and insert them into the right leaf
    for (size_t cp = 0; cp < input_->points.size(); ++cp) {
      if (!input_->is_dense)
        // Check if the point is invalid
        if (!pcl_isfinite(input_->points[cp].x) ||
            !pcl_isfinite(input_->points[cp].y) ||
            !pcl_isfinite(input_->points[cp].z))
          continue;

      int ijk0 =
          static_cast<int>(floor(input_->points[cp].x * inverse_leaf_size_[0]) -
                           static_cast<float>(min_b_[0]));
      int ijk1 =
          static_cast<int>(floor(input_->points[cp].y * inverse_leaf_size_[1]) -
                           static_cast<float>(min_b_[1]));
      int ijk2 =
          static_cast<int>(floor(input_->points[cp].z * inverse_leaf_size_[2]) -
                           static_cast<float>(min_b_[2]));

      // Compute the centroid leaf index
      int idx = ijk0 * divb_mul_[0] + ijk1 * divb_mul_[1] + ijk2 * divb_mul_[2];

      // int idx = (((input_->points[cp].getArray4fmap () *
      // inverse_leaf_size_).template cast<int> ()).matrix () - min_b_).dot
      // (divb_mul_);
      Leaf &leaf = leaves_[idx];
      if (leaf.nr_points == 0) {
        leaf.centroid.resize(centroid_size);
        leaf.centroid.setZero();
      }

      Eigen::Vector3d pt3d(input_->points[cp].x, input_->points[cp].y,
                           input_->points[cp].z);
      // Accumulate point sum for centroid calculation
      leaf.mean_ += pt3d;
      // Accumulate x*xT for single pass covariance calculation
      leaf.cov_ += pt3d * pt3d.transpose();

      // Do we need to process all the fields?
      if (!downsample_all_data_) {
        Eigen::Vector4f pt(input_->points[cp].x, input_->points[cp].y,
                           input_->points[cp].z, 0);
        leaf.centroid.template head<4>() += pt;
      } else {
        // Copy all the fields
        Eigen::VectorXf centroid = Eigen::VectorXf::Zero(centroid_size);
        // ---[ RGB special case
        if (rgba_index >= 0) {
          // Fill r/g/b data, assuming that the order is BGRA
          int rgb;
          memcpy(&rgb,
                 reinterpret_cast<const char *>(&input_->points[cp]) +
                     rgba_index,
                 sizeof(int));
          centroid[centroid_size - 3] =
              static_cast<float>((rgb >> 16) & 0x0000ff);
          centroid[centroid_size - 2] =
              static_cast<float>((rgb >> 8) & 0x0000ff);
          centroid[centroid_size - 1] = static_cast<float>((rgb)&0x0000ff);
        }
        pcl::for_each_type<FieldList>(
            pcl::NdCopyPointEigenFunctor<PointT>(input_->points[cp], centroid));
        leaf.centroid += centroid;
      }
      ++leaf.nr_points;
    }
  }

  // Second pass: go over all leaves and compute centroids and covariance
  // matrices
  output.points.reserve(leaves_.size());
  if (searchable_)
    voxel_centroids_leaf_indices_.reserve(leaves_.size());
  int cp = 0;
  if (save_leaf_layout_)
    leaf_layout_.resize(div_b_[0] * div_b_[1] * div_b_[2], -1);

  // Eigen values and vectors calculated to prevent near singluar matrices
  Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> eigensolver;
  Eigen::Matrix3d eigen_val;
  Eigen::Vector3d pt_sum;

  // Eigen values less than a threshold of max eigen value are inflated to a set
  // fraction of the max eigen value.
  double min_covar_eigvalue;

  for (auto it = leaves_.begin(); it != leaves_.end(); ++it) {

    // Normalize the centroid
    Leaf &leaf = it->second;

    // Normalize the centroid
    leaf.centroid /= static_cast<float>(leaf.nr_points);
    // Point sum used for single pass covariance calculation
    pt_sum = leaf.mean_;
    // Normalize mean
    leaf.mean_ /= leaf.nr_points;

    // If the voxel contains sufficient points, its covariance is calculated and
    // is added to the voxel centroids and output clouds. Points with less than
    // the minimum points will have a can not be accuratly approximated using a
    // normal distribution.
    if (leaf.nr_points >= min_points_per_voxel_) {
      if (save_leaf_layout_)
        leaf_layout_[it->first] = cp++;

      output.push_back(PointT());

      // Do we need to process all the fields?
      if (!downsample_all_data_) {
        output.points.back().x = leaf.centroid[0];
        output.points.back().y = leaf.centroid[1];
        output.points.back().z = leaf.centroid[2];
      } else {
        pcl::for_each_type<FieldList>(
            pcl::NdCopyEigenPointFunctor<PointT>(leaf.centroid, output.back()));
        // ---[ RGB special case
        if (rgba_index >= 0) {
          // pack r/g/b into rgb
          float r = leaf.centroid[centroid_size - 3],
                g = leaf.centroid[centroid_size - 2],
                b = leaf.centroid[centroid_size - 1];
          int rgb = (static_cast<int>(r)) << 16 | (static_cast<int>(g)) << 8 |
                    (static_cast<int>(b));
          memcpy(reinterpret_cast<char *>(&output.points.back()) + rgba_index,
                 &rgb, sizeof(float));
        }
      }

      // Stores the voxel indice for fast access searching
      if (searchable_)
        voxel_centroids_leaf_indices_.push_back(static_cast<int>(it->first));

      // Single pass covariance calculation
      leaf.cov_ =
          (leaf.cov_ - 2 * (pt_sum * leaf.mean_.transpose())) / leaf.nr_points +
          leaf.mean_ * leaf.mean_.transpose();
      leaf.cov_ *= (leaf.nr_points - 1.0) / leaf.nr_points;

      // Normalize Eigen Val such that max no more than 100x min.
      eigensolver.compute(leaf.cov_);
      eigen_val = eigensolver.eigenvalues().asDiagonal();
      leaf.evecs_ = eigensolver.eigenvectors();

      if (eigen_val(0, 0) < 0 || eigen_val(1, 1) < 0 || eigen_val(2, 2) <= 0) {
        leaf.nr_points = -1;
        continue;
      }

      // Avoids matrices near singularities (eq 6.11)[Magnusson 2009]

      min_covar_eigvalue = min_covar_eigvalue_mult_ * eigen_val(2, 2);
      if (eigen_val(0, 0) < min_covar_eigvalue) {
        eigen_val(0, 0) = min_covar_eigvalue;

        if (eigen_val(1, 1) < min_covar_eigvalue) {
          eigen_val(1, 1) = min_covar_eigvalue;
        }

        leaf.cov_ = leaf.evecs_ * eigen_val * leaf.evecs_.inverse();
      }
      leaf.evals_ = eigen_val.diagonal();

      leaf.icov_ = leaf.cov_.inverse();
      if (leaf.icov_.maxCoeff() == std::numeric_limits<float>::infinity() ||
          leaf.icov_.minCoeff() == -std::numeric_limits<float>::infinity()) {
        leaf.nr_points = -1;
      }
    }
  }

  output.width = static_cast<uint32_t>(output.points.size());
}

////////////////////////////////////////////////////////////////////////////////
template <typename PointT>
int pclomp::VoxelGridCovariance<PointT>::getNeighborhoodAtPoint(
    const Eigen::MatrixXi &relative_coordinates, const PointT &reference_point,
    std::vector<LeafConstPtr> &neighbors) const {
  neighbors.clear();

  // Find displacement coordinates
  Eigen::Vector4i ijk(
      static_cast<int>(floor(reference_point.x / leaf_size_[0])),
      static_cast<int>(floor(reference_point.y / leaf_size_[1])),
      static_cast<int>(floor(reference_point.z / leaf_size_[2])), 0);
  Eigen::Array4i diff2min = min_b_ - ijk;
  Eigen::Array4i diff2max = max_b_ - ijk;
  neighbors.reserve(relative_coordinates.cols());

  // Check each neighbor to see if it is occupied and contains sufficient points
  // Slower than radius search because needs to check 26 indices
  for (int ni = 0; ni < relative_coordinates.cols(); ni++) {
    Eigen::Vector4i displacement =
        (Eigen::Vector4i() << relative_coordinates.col(ni), 0).finished();
    // Checking if the specified cell is in the grid
    if ((diff2min <= displacement.array()).all() &&
        (diff2max >= displacement.array()).all()) {
      auto leaf_iter =
          leaves_.find(((ijk + displacement - min_b_).dot(divb_mul_)));
      if (leaf_iter != leaves_.end() &&
          leaf_iter->second.nr_points >= min_points_per_voxel_) {
        LeafConstPtr leaf = &(leaf_iter->second);
        neighbors.push_back(leaf);
      }
    }
  }

  return (static_cast<int>(neighbors.size()));
}

////////////////////////////////////////////////////////////////////////////////
template <typename PointT>
int pclomp::VoxelGridCovariance<PointT>::getNeighborhoodAtPoint(
    const PointT &reference_point, std::vector<LeafConstPtr> &neighbors) const {
  neighbors.clear();

  // Find displacement coordinates
  Eigen::MatrixXi relative_coordinates = pcl::getAllNeighborCellIndices();
  return getNeighborhoodAtPoint(relative_coordinates, reference_point,
                                neighbors);
}

////////////////////////////////////////////////////////////////////////////////
template <typename PointT>
int pclomp::VoxelGridCovariance<PointT>::getNeighborhoodAtPoint7(
    const PointT &reference_point, std::vector<LeafConstPtr> &neighbors) const {
  neighbors.clear();

  Eigen::MatrixXi relative_coordinates(3, 7);
  relative_coordinates.setZero();
  relative_coordinates(0, 1) = 1;
  relative_coordinates(0, 2) = -1;
  relative_coordinates(1, 3) = 1;
  relative_coordinates(1, 4) = -1;
  relative_coordinates(2, 5) = 1;
  relative_coordinates(2, 6) = -1;

  return getNeighborhoodAtPoint(relative_coordinates, reference_point,
                                neighbors);
}

//////////////////////////////////////////////////////////////////////////////////////////
template <typename PointT>
int pclomp::VoxelGridCovariance<PointT>::getNeighborhoodAtPoint1(
    const PointT &reference_point, std::vector<LeafConstPtr> &neighbors) const {
  neighbors.clear();
  return getNeighborhoodAtPoint(Eigen::MatrixXi::Zero(3, 1), reference_point,
                                neighbors);
}

////////////////////////////////////////////////////////////////////////////////
template <typename PointT>
void pclomp::VoxelGridCovariance<PointT>::getDisplayCloud(
    pcl::PointCloud<pcl::PointXYZ> &cell_cloud) {
  cell_cloud.clear();

  int pnt_per_cell = 1000;
  boost::mt19937 rng;
  boost::normal_distribution<> nd(0.0, leaf_size_.head(3).norm());
  boost::variate_generator<boost::mt19937 &, boost::normal_distribution<>>
      var_nor(rng, nd);

  Eigen::LLT<Eigen::Matrix3d> llt_of_cov;
  Eigen::Matrix3d cholesky_decomp;
  Eigen::Vector3d cell_mean;
  Eigen::Vector3d rand_point;
  Eigen::Vector3d dist_point;

  // Generate points for each occupied voxel with sufficient points.
  for (auto it = leaves_.begin(); it != leaves_.end(); ++it) {
    Leaf &leaf = it->second;

    if (leaf.nr_points >= min_points_per_voxel_) {
      cell_mean = leaf.mean_;
      llt_of_cov.compute(leaf.cov_);
      cholesky_decomp = llt_of_cov.matrixL();

      // Random points generated by sampling the normal distribution given by
      // voxel mean and covariance matrix
      for (int i = 0; i < pnt_per_cell; i++) {
        rand_point = Eigen::Vector3d(var_nor(), var_nor(), var_nor());
        dist_point = cell_mean + cholesky_decomp * rand_point;
        cell_cloud.push_back(pcl::PointXYZ(static_cast<float>(dist_point(0)),
                                           static_cast<float>(dist_point(1)),
                                           static_cast<float>(dist_point(2))));
      }
    }
  }
}

#define PCL_INSTANTIATE_VoxelGridCovariance(T)                                 \
  template class PCL_EXPORTS pcl::VoxelGridCovariance<T>;

#endif // PCL_VOXEL_GRID_COVARIANCE_IMPL_H_