text:
//          Copyright Rein Halbersma 2010-2020.
// 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)

#include <core/rules/precedence.hpp>            // is_consistent
#include <dctl/core/rules/czech.hpp>            // czech
#include <dctl/core/state/piece.hpp>            // king, pawn
#include <dctl/core/rules/type_traits.hpp>      // is_backward_pawn_jump, king_range_category, long_ranged_tag, is_trivial_precedence, jumpsep, is_absolute_modality, equal_to, less
#include <boost/test/unit_test.hpp>             // BOOST_AUTO_TEST_SUITE, BOOST_AUTO_TEST_CASE, BOOST_AUTO_TEST_SUITE_END
#include <vector>                               // vector

using namespace dctl::core;

BOOST_AUTO_TEST_SUITE(RulesCzech)

using T = czech;

BOOST_AUTO_TEST_CASE(RuleTraits)
{
        static_assert(!is_backward_pawn_jump_v<T>);
        static_assert(is_long_ranged_king_v<T>);

        static_assert(jumpsep_v<T> == ':');

        static_assert(!is_trivial_precedence_v<T>);
        static_assert(is_modality_precedence_v<T>);

        struct Action
        {
                using rules_type [[maybe_unused]] = czech;
                piece m_with;
                constexpr auto is_with_king() const noexcept { return m_with == piece::king; }
        };

        auto const moves = std::vector<Action>
        {
                { piece::pawn },
                { piece::king }
        };

        BOOST_CHECK(xxx_precedence::is_consistent(moves));
}

BOOST_AUTO_TEST_SUITE_END()

ext: cpp


text:
// vim:ts=2:sw=2:expandtab:autoindent:filetype=cpp:
/*
The MIT License

Copyright (c) 2008, 2009 Flusspferd contributors (see "CONTRIBUTORS" or
                                       http://flusspferd.org/contributors.txt)

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.
*/

#ifndef FLUSSPFERD_IO_FILESYSTEM_BASE_HPP
#define FLUSSPFERD_IO_FILESYSTEM_BASE_HPP

#include "../string.hpp"
#include "../object.hpp"
#include <boost/filesystem.hpp>

namespace flusspferd {

void load_filesystem_base_module(object container);

namespace io {

namespace fs_base {

  object raw_open(char const* name, value mode, value permissions);

  string canonical(std::string const &path);
  boost::filesystem::path canonicalize(boost::filesystem::path in);

  void move(std::string const &source, std::string const &target);
  void remove(std::string const &target);
  void touch(std::string const &path, object mtime);

  void make_directory(std::string const &target);
  void remove_directory(std::string const &target);

  string working_directory();
  void change_working_directory(std::string const &path);

  string owner(std::string const &path);
  void change_owner(std::string const &path, std::string owner);

#ifndef WIN32
  string read_link(std::string const &link);
  void link(std::string const &source, std::string const &target);
  void hard_link(std::string const &source, std::string const &target);
#endif

  bool exists(std::string const &path);
  bool is_file(std::string const &path);
  bool is_directory(std::string const &path);
  bool is_link(std::string const &path);
  bool is_readable(std::string const &path);
  bool is_writeable(std::string const &path);

  // Return a double since 1gb file limit on size would suck
  double size(std::string const &file);
  value last_modified(std::string const &path);
  bool same(std::string const &source, std::string const &target);

  array list(std::string const &path);
  object iterate(std::string const &path);
}

}}

#endif

ext: hpp

max_stars_repo_name: ruediger/flusspferd

max_issues_repo_name: ngoyal/flusspferd

max_forks_repo_name: ngoyal/flusspferd


text:
#pragma once

#include <functional>
#include <boost/algorithm/string.hpp>

#include <supermarx/message/product_summary.hpp>

#include <karl/util/levenshtein.hpp>
#include <karl/util/hungarian_fast.hpp>

namespace supermarx
{

class similarity
{
public:
	struct valuation
	{
		static constexpr size_t N = 3;

		std::array<float, N> data;

		valuation(decltype(data) _data)
		: data(_data)
		{}

		float collapse() const
		{
			static constexpr std::array<float, N> weights = {
				0.6f, 0.2f, 0.2f
			};

			float similarity = 0.0f;

			for(size_t i = 0; i < N; ++i)
				similarity += weights[i] * data[i];

			return similarity;
		}
	};

private:
	template<typename T>
	static inline float max_set(std::vector<T> const& xs, std::function<float(T const&)> f)
	{
		float result = std::numeric_limits<float>::min();

		for(auto const& x : xs)
			result = std::max(result, f(x));

		return result;
	}

	static inline float textual_compare(std::string const& x, std::string const& y)
	{
		std::vector<std::string> xs, ys;
		boost::split(xs, x, boost::is_any_of(" "));
		boost::split(ys, y, boost::is_any_of(" "));

		if(xs.size() > ys.size())
			std::swap(xs, ys);

		typedef float sim_t;
		similarity_matrix<sim_t> sm(ys.size(), xs.size());
		for(size_t yi = 0; yi < ys.size(); ++yi)
		{
			for(size_t xi = 0; xi < xs.size(); ++xi)
			{
				std::string const& ye = ys[yi];
				std::string const& xe = xs[xi];

				size_t distance_yx(levenshtein(ye, xe));
				assert(distance_yx >= 0 && distance_yx <= std::numeric_limits<sim_t>::max());

				sm(yi, xi) = static_cast<sim_t>(std::max(ye.size(), xe.size()) - distance_yx) / (sim_t)std::max(ye.size(), xe.size());
			}
		}

		hungarian_fast<decltype(sm)::similarity> h(sm);
		auto matching = h.produce();

		float similarity = 0.0f;
		for(auto matching_e : matching)
		{
			size_t yi = matching_e.first;
			size_t xi = matching_e.second;

			similarity += (float)sm(yi, xi);
		}

		float sim_min = std::min(ys.size(), xs.size()); // Due to std::swap ys.size() will always be bigger, Hungarian will thus always yield xs.size() elements
		float sim_max = std::max(ys.size(), xs.size());

		return 0.9f * similarity / sim_min + 0.1f * similarity / sim_max;
	}

	static inline float numeric_compare(float x, float y)
	{
		float result = 1.0f - std::abs(x - y) / std::max(x, y);
		assert(result >= 0.0f && result <= 1.0f);

		return result;
	}

public:
	similarity() = delete;

	static inline valuation exec(message::product_summary const& x, message::product_summary const& y)
	{
		return valuation({
			textual_compare(boost::to_lower_copy(x.name), boost::to_lower_copy(y.name)),
			numeric_compare(x.orig_price, y.orig_price),
			(x.volume_measure == y.volume_measure && x.volume == y.volume) ? 1.0f : 0.0f
		});
	}
};

}

max_stars_repo_name: Supermarx/karl

max_issues_repo_name: Supermarx/karl

max_forks_repo_name: Supermarx/karl


text: #include <boost/local_function/aux_/preprocessor/traits/decl_/index.hpp>

max_stars_repo_name: miathedev/BoostForArduino

max_issues_repo_name: miathedev/BoostForArduino

max_forks_repo_name: miathedev/BoostForArduino


text:
/// Copyright Barzer LLC 2012
/// Code is property Barzer for authorized use only
///
#include <zurch_python.h>
#include <zurch_tokenizer.h>
#include <zurch_classifier.h>

#include <barzer_universe.h>
#include <barzer_shell.h>
#include <boost/python.hpp>
#include <barzer_python.h>
#include <ay/ay_cmdproc.h>
#include <ay_translit_ru.h>
#include <boost/python/list.hpp>
#include <boost/python/dict.hpp>
#include <boost/python/object.hpp>
#include <boost/python/stl_iterator.hpp>

#include <barzer_bzspell.h>
#include <barzer_parse.h>
#include <barzer_server_response.h>
#include <barzer_el_chain.h>
#include <barzer_emitter.h>
#include <boost/variant.hpp>

#include <util/pybarzer.h>
#include <ay/ay_python.h>

#include <autotester/barzer_at_comparators.h>

using boost::python::stl_input_iterator ;
using namespace boost::python;
using namespace boost;

namespace zurch {
struct PythonTokenizer {
    ZurchTokenizer tokenizer;
    const ZurchPython& zp;
    PythonTokenizer(const ZurchPython& z): zp(z){}
    list tokenize( const std::string & s )
    {
        list retList;
        ZurchTokenVec tokVec;
        tokenizer.tokenize( tokVec, s.c_str(), s.length() );
        for( ZurchTokenVec::const_iterator i = tokVec.begin(); i!= tokVec.end(); ++i ) {
            retList.append(i->str);
        }
        return retList;
    }

};
struct PythonNormalizer {
    const ZurchPython& zp;
    ZurchTokenizer tokenizer;
    ZurchWordNormalizer normalizer;

    PythonNormalizer(const ZurchPython& z): zp(z),
    normalizer(zp.barzerPython->guaranteeUniverse().getBZSpell())
    {}
    /// must delete the returned pointer
    char* internal_tokenize( ZurchTokenVec& tokVec, const std::string& sorig )
    {
        size_t lc_sz = sorig.length();
        int lang = barzer::Lang::getLangNoUniverse( sorig.c_str(), sorig.length() );

        char * lc= new char[ lc_sz+1 ] ;

        // memcpy( lc, sorig.c_str(), sorig.length() );
        barzer::Lang::lowLevelNormalization( lc, lc_sz, sorig.c_str(), sorig.length() );
        lc[ lc_sz ] =0;

        barzer::Lang::stringToLower( lc, lc_sz, lang );

        tokenizer.tokenize( tokVec, lc, lc_sz );
        return lc;
    }

    list tokenize( const std::string& s )
    {
        list retList;
        ZurchTokenVec tokVec;
        tokenizer.tokenize( tokVec, s.c_str(), s.length() );
        for( ZurchTokenVec::const_iterator i = tokVec.begin(); i!= tokVec.end(); ++i ) {
            retList.append(i->str);
        }
        return retList;
    }
    void internal_normalize( std::vector< std::string > & vec, const std::string& s )
    {
        ZurchTokenVec tokVec;
        /// this is needed so that the char* inside tokVec are valid
        std::auto_ptr<char>( internal_tokenize( tokVec, s ) );
        std::string norm;
        ZurchWordNormalizer::NormalizerEnvironment normEnv;
        for( ZurchTokenVec::const_iterator i = tokVec.begin(); i!= tokVec.end(); ++i ) {
            norm.clear();
            if( i->str.length() > 3 ) {
                normalizer.normalize(norm,i->str.c_str(),normEnv);
                vec.push_back( norm );
            } else
                vec.push_back(i->str);
        }
    }
    list normalize( const std::string& sorig )
    {

        list retList;
        ZurchTokenVec tokVec;
        /// this is needed so that the char* inside tokVec are valid
        std::auto_ptr<char>( internal_tokenize( tokVec, sorig ) );
        std::string norm;
        ZurchWordNormalizer::NormalizerEnvironment normEnv;
        for( ZurchTokenVec::const_iterator i = tokVec.begin(); i!= tokVec.end(); ++i ) {
            norm.clear();
            if( i->str.length() > 3 ) {
                normalizer.normalize(norm,i->str.c_str(),normEnv);
                retList.append(norm);
            } else
                retList.append(i->str);

        }
        return retList;
    }
};

struct PythonClassifier {
    ZurchPython& zp;

    PythonNormalizer normalizer;
    ZurchTrainerAndClassifier* cfier;

    /// standalone tokenizer
    list tokenize( const std::string& s )
        { return normalizer.tokenize(s); }
    /// standalone normalizer
    list normalize( const std::string& s )
        { return normalizer.normalize(s); }

    /// returns list of tag probabilities for all known tags
    /// for tags, whose probability exceeds threshold
    list classify( double threshold, const std::string& s )
    {
        list retList;
        ExtractedFeatureMap featureMap;
        cfier->extractFeaturesClassification( featureMap, s.c_str());
        DocClassifier::TagIdToProbabilityMap tagProbMap;
        cfier->classifyAll( tagProbMap, featureMap );

        for( DocClassifier::TagIdToProbabilityMap::const_iterator i= tagProbMap.begin(); i!= tagProbMap.end(); ++i ) {
            if( i->second > threshold ) {
                list l;
                l.append( i->first );
                l.append( i->second);
                retList.append(l);
            }
        }
        return retList;
    }
    /// learning function
    /// accumulates tag data for a document. tagList contains tags the document is tagged with
    list accumulate( list tagList, const std::string& s )
    {
        ExtractedFeatureMap featureMap;
        cfier->extractFeaturesLearning( featureMap, s.c_str() );
        ay::python::object_extract xtr;
        for( size_t i = 0, i_end = len(tagList); i< i_end; ++i ) {
            std::string tmp;
            xtr( tmp, tagList[i] );
            uint32_t tagId = static_cast<uint32_t>(atoi(tmp.c_str()));
            cfier->accumulate(tagId,featureMap,true);
        }
        return list();
    }
    /// should be called when accumulation is complete and stats need to be synced with the accumulators
    list computeStats()
    {
        cfier->computeStats();
        return list();
    }
    void print()
    {
        PrintContext ctxt( &(zp.barzerPython->getGP().getStringPool()) );
        cfier->print( std::cerr, ctxt );
    }
    PythonClassifier(ZurchPython& z):
        zp(z),
        normalizer(zp) ,
        cfier(0)
    {
        cfier= new ZurchTrainerAndClassifier( z.barzerPython->getGP().getStringPool() );
        cfier->init(ZurchTrainerAndClassifier::EXTRACTOR_TYPE_NORMALIZING,zp.barzerPython->guaranteeUniverse().getBZSpell());
        // cfier->getTokenizer().tokenizer().addSingleCharSeparators("()-!?");
    }
    ~PythonClassifier() { delete cfier; }
};
PythonClassifier* ZurchPython::mkClassifier()
{
    return new PythonClassifier(*this);
}
PythonTokenizer*  ZurchPython::mkTokenizer() const
{
    return new PythonTokenizer(*this);
}
PythonNormalizer* ZurchPython::mkNormalizer() const
{
    return new PythonNormalizer(*this);
}

ZurchPython::~ZurchPython() {delete barzerPython; }

ZurchPython::ZurchPython():
    barzerPython(new barzer::BarzerPython)
{
    barzerPython->guaranteeUniverse();
}

void ZurchPython::initInstance()
{
    barzerPython->guaranteeUniverse();
}
void ZurchPython::init()
{

    boost::python::class_<zurch::ZurchPython>( "Zurch" )
        .def( "mkClassifier", &ZurchPython::mkClassifier, return_value_policy<manage_new_object>() )
        .def( "mkTokenizer", &ZurchPython::mkTokenizer, return_value_policy<manage_new_object>() )
        .def( "mkNormalizer", &ZurchPython::mkNormalizer, return_value_policy<manage_new_object>() );

    boost::python::class_<zurch::PythonTokenizer>( "tokenizer", no_init )
        .def( "tokenize", &zurch::PythonTokenizer::tokenize );
    boost::python::class_<zurch::PythonNormalizer>( "normalizer", no_init )
        .def( "normalize", &zurch::PythonNormalizer::normalize );

    boost::python::class_<zurch::PythonClassifier>( "classifier", no_init )
        .def( "tokenize", &zurch::PythonClassifier::tokenize )
        .def( "normalize", &zurch::PythonClassifier::normalize )
        .def( "accumulate", &zurch::PythonClassifier::accumulate )
        .def( "computeStats", &zurch::PythonClassifier::computeStats )
        .def( "show", &zurch::PythonClassifier::print )
        .def( "classify", &zurch::PythonClassifier::classify );
    // std::cerr << "pythonInit" << std::endl;
}

} // namespace zurch

max_stars_repo_name: barzerman/barzer

max_issues_repo_name: barzerman/barzer

max_forks_repo_name: barzerman/barzer


text:
#include "cpu/graph_proc.h"

#include <set>
#include <vector>
#include <numeric> //std::iota
#include <algorithm>
#include <random>
#include <map>
#include <queue>
#include <functional>

#include <Eigen/Dense>
#include <cfloat>

using std::vector;

namespace graph_proc {

/**
 * \brief Compile a vertex mask that can be used to erode the provided mesh (iteratively mask out vertices at surface discontinuities, leave only non-eroded vertices)
 * \param vertex_positions array with vertex positions of dimentions (N, 3)
 * \param face_indices
 * \param iteration_count
 * \param min_neighbors minimum number of neighbors considered for erosion
 * \return (N, 1) array containing the erosion mask, such that, if applied to the original vertex postions array, would leave only non-eroded vertices
 */
py::array_t<bool> get_vertex_erosion_mask(const py::array_t<float>& vertex_positions, const py::array_t<int>& face_indices,
                                          int iteration_count, int min_neighbors) {
	int vertex_count = vertex_positions.shape(0);
	int face_count = face_indices.shape(0);

	// Init output
	py::array_t<bool> non_eroded_vertices = py::array_t<bool>({vertex_count, 1});
	std::vector<bool> non_eroded_vertices_vec(vertex_count, false);

	// Init list of eroded face indices with original list
	std::vector<Eigen::Vector3i> non_erorded_face_indices_vec;
	non_erorded_face_indices_vec.reserve(face_count);
	for (int face_idx = 0; face_idx < face_count; ++face_idx) {
		Eigen::Vector3i face(*face_indices.data(face_idx, 0), *face_indices.data(face_idx, 1), *face_indices.data(face_idx, 2));
		non_erorded_face_indices_vec.push_back(face);
	}

	// Erode mesh for a total of iteration_count
	for (int i = 0; i < iteration_count; i++) {
		face_count = non_erorded_face_indices_vec.size();

		// We compute the number of neighboring vertices for each vertex.
		vector<int> vertex_neighbor_counts(vertex_count, 0);
		for (int face_index = 0; face_index < face_count; face_index++) {
			const auto& face = non_erorded_face_indices_vec[face_index];
			vertex_neighbor_counts[face[0]] += 1;
			vertex_neighbor_counts[face[1]] += 1;
			vertex_neighbor_counts[face[2]] += 1;
		}

		std::vector<Eigen::Vector3i> tmp;
		tmp.reserve(face_count);

		for (int face_index = 0; face_index < face_count; face_index++) {
			const auto& face = non_erorded_face_indices_vec[face_index];
			if (vertex_neighbor_counts[face[0]] >= min_neighbors &&
			    vertex_neighbor_counts[face[1]] >= min_neighbors &&
			    vertex_neighbor_counts[face[2]] >= min_neighbors) {
				tmp.push_back(face);
			}
		}

		// We kill the faces with border vertices.
		non_erorded_face_indices_vec.clear();
		non_erorded_face_indices_vec = std::move(tmp);
	}

	// Mark non isolated vertices as not eroded.
	face_count = non_erorded_face_indices_vec.size();

	for (int i = 0; i < face_count; i++) {
		const auto& face = non_erorded_face_indices_vec[i];
		non_eroded_vertices_vec[face[0]] = true;
		non_eroded_vertices_vec[face[1]] = true;
		non_eroded_vertices_vec[face[2]] = true;
	}

	// Store into python array
	for (int i = 0; i < vertex_count; i++) {
		*non_eroded_vertices.mutable_data(i, 0) = non_eroded_vertices_vec[i];
	}

	return non_eroded_vertices;
}

py::tuple sample_nodes(
		const py::array_t<float>& vertex_positions_in, const py::array_t<bool>& vertex_erosion_mask_in,
		float node_coverage, const bool use_only_non_eroded_indices = true,
		const bool random_shuffle = true
) {
	assert(vertex_positions_in.ndim() == 2);
	assert(vertex_positions_in.shape(1) == 3);

	float node_coverage_squared = node_coverage * node_coverage;
	int vertex_count = vertex_positions_in.shape(0);

	// create list of shuffled indices
	std::vector<int> shuffled_vertices(vertex_count);
	std::iota(std::begin(shuffled_vertices), std::end(shuffled_vertices), 0);

	if (random_shuffle) {
		std::default_random_engine re{std::random_device{}()};
		std::shuffle(std::begin(shuffled_vertices), std::end(shuffled_vertices), re);
	}

	struct NodeInformation {
		Eigen::Vector3f vertex_position;
		int vertex_index;
	};

	std::vector<NodeInformation> node_information_vector;

	for (int vertex_index : shuffled_vertices) {
		Eigen::Vector3f point(*vertex_positions_in.data(vertex_index, 0),
		                      *vertex_positions_in.data(vertex_index, 1),
		                      *vertex_positions_in.data(vertex_index, 2));

		if (use_only_non_eroded_indices && !(*vertex_erosion_mask_in.data(vertex_index))) {
			continue;
		}

		bool is_node = true;
		for (const auto& node_information : node_information_vector) {
			if ((point - node_information.vertex_position).squaredNorm() <= node_coverage_squared) {
				is_node = false;
				break;
			}
		}

		if (is_node) {
			node_information_vector.push_back({point, vertex_index});
		}
	}

	py::array_t<float> node_positions_out({static_cast<ssize_t>(node_information_vector.size()), static_cast<ssize_t>(3)});
	py::array_t<int> node_indices_out({static_cast<ssize_t>(node_information_vector.size()), static_cast<ssize_t>(1)});

	int node_index = 0;
	for (const auto& node_information : node_information_vector) {
		*node_positions_out.mutable_data(node_index, 0) = node_information.vertex_position.x();
		*node_positions_out.mutable_data(node_index, 1) = node_information.vertex_position.y();
		*node_positions_out.mutable_data(node_index, 2) = node_information.vertex_position.z();
		*node_indices_out.mutable_data(node_index, 0) = node_information.vertex_index;
		node_index++;
	}

	return py::make_tuple(node_positions_out, node_indices_out);
}

/**
 * Custom comparison operator for shortest-path priority queue.
 */
struct CustomCompare {
	bool operator()(const std::pair<int, float>& left, const std::pair<int, float>& right) {
		return left.second > right.second;
	}
};

inline float compute_anchor_weight(const Eigen::Vector3f& point_position, const Eigen::Vector3f& node_position, float node_coverage) {
	return std::exp(-(node_position - point_position).squaredNorm() / (2.f * node_coverage * node_coverage));
}

inline float compute_anchor_weight(float dist, float node_coverage) {
	return std::exp(-(dist * dist) / (2.f * node_coverage * node_coverage));
}

inline float compute_anchor_weight_square_distance(float square_distance, float node_coverage) {
	return std::exp(-(square_distance) / (2.f * node_coverage * node_coverage));
}

template<typename TVertexCheckLambda>
inline void compute_edges_shortest_path_generic(
		const py::array_t<float>& vertex_positions_in,
		const py::array_t<int>& face_indices_in,
		const py::array_t<int>& node_indices_in,
		const int max_neighbor_count,
		const float node_coverage,
		py::array_t<int>& graph_edges_out,
		py::array_t<float>& graph_edges_weights_out,
		py::array_t<float>& graph_edges_distances_out,
		py::array_t<float>& node_to_vertex_distances_out,
		const bool enforce_total_num_neighbors,
		TVertexCheckLambda&& vertex_is_valid
) {

	assert(graph_edges_out.ndim() == 2);
	assert(graph_edges_weights_out.ndim() == 2);
	assert(graph_edges_distances_out.ndim() == 2);
	assert(node_to_vertex_distances_out.ndim() == 2);

	int vertex_count = vertex_positions_in.shape(0);
	int face_count = face_indices_in.shape(0);
	int node_count = node_indices_in.shape(0);

	assert(graph_edges_out.shape(0) == node_count);
	assert(graph_edges_out.shape(1) == max_neighbor_count);
	assert(graph_edges_weights_out.shape(0) == node_count);
	assert(graph_edges_weights_out.shape(1) == max_neighbor_count);
	assert(graph_edges_distances_out.shape(0) == node_count);
	assert(graph_edges_distances_out.shape(1) == max_neighbor_count);
	assert(node_to_vertex_distances_out.shape(0) == node_count);

	float max_influence = 2.f * node_coverage;

	// Preprocess vertex neighbors.
	vector<std::set<int>> vertex_neighbors(vertex_count);
	for (int face_index = 0; face_index < face_count; face_index++) {
		for (int face_vertex_index = 0; face_vertex_index < 3; face_vertex_index++) {
			int vertex_index = *face_indices_in.data(face_index, face_vertex_index);

			for (int other_face_vertex_index = 0; other_face_vertex_index < 3; other_face_vertex_index++) {
				int other_vertex_index = *face_indices_in.data(face_index, other_face_vertex_index);

				if (vertex_index == other_vertex_index) continue;
				vertex_neighbors[vertex_index].insert(other_vertex_index);
			}
		}
	}

	// Compute inverse vertex -> node relationship.
	vector<int> vertex_to_node_map(vertex_count, -1);

	for (int node_index = 0; node_index < node_count; node_index++) {
		int vertex_index = *node_indices_in.data(node_index);
		if (vertex_index >= 0) {
			vertex_to_node_map[vertex_index] = node_index;
		}
	}

	//TODO: parallelize or remove pragma
	// #pragma omp parallel for
	for (int node_index = 0; node_index < node_count; node_index++) {
		// vertex queue
		std::priority_queue<
				std::pair<int, float>,
				vector<std::pair<int, float>>,
				CustomCompare
		> next_vertices_with_distances;

		std::set<int> visited_vertices;

		// Add node vertex as the first vertex to be visited.
		int node_vertex_index = *node_indices_in.data(node_index);
		if (node_vertex_index < 0) continue;
		next_vertices_with_distances.push(std::make_pair(node_vertex_index, 0.f));

		// Traverse all neighbors in the monotonically increasing order.
		vector<int> neighbor_node_indices;
		vector<float> neighbor_node_weights;
		vector<float> neighbor_node_distances;
		while (!next_vertices_with_distances.empty()) {
			auto next_vertex = next_vertices_with_distances.top();
			next_vertices_with_distances.pop();

			int next_vertex_index = next_vertex.first;
			float next_vertex_distance = next_vertex.second;

			// We skip the vertex, if it was already visited before.
			if (visited_vertices.find(next_vertex_index) != visited_vertices.end()) continue;

			assert(std::forward<TVertexCheckLambda>(vertex_is_valid)(next_vertex_index));

			// We check if the vertex is a node.
			int next_node_id = vertex_to_node_map[next_vertex_index];
			if (next_node_id >= 0 && next_node_id != node_index) {
				neighbor_node_indices.push_back(next_node_id);
				neighbor_node_weights.push_back(compute_anchor_weight(next_vertex_distance, node_coverage));
				neighbor_node_distances.push_back(next_vertex_distance);
				if (neighbor_node_indices.size() >= max_neighbor_count) break;
			}

			// Note down the node-vertex distance.
			*node_to_vertex_distances_out.mutable_data(node_index, next_vertex_index) = next_vertex_distance;

			// We visit the vertex, and check all his neighbors.
			// We add only valid vertices under a certain distance
			visited_vertices.insert(next_vertex_index);
			Eigen::Vector3f next_vertex_position(*vertex_positions_in.data(next_vertex_index, 0),
			                                     *vertex_positions_in.data(next_vertex_index, 1),
			                                     *vertex_positions_in.data(next_vertex_index, 2));

			const auto& next_neighbors = vertex_neighbors[next_vertex_index];
			for (int neighbor_index : next_neighbors) {

				if (!std::forward<TVertexCheckLambda>(vertex_is_valid)(neighbor_index)) {
					continue;
				}

				Eigen::Vector3f neighbor_vertex_position(*vertex_positions_in.data(neighbor_index, 0),
				                                         *vertex_positions_in.data(neighbor_index, 1),
				                                         *vertex_positions_in.data(neighbor_index, 2));
				float distance = next_vertex_distance + (next_vertex_position - neighbor_vertex_position).norm();

				if (enforce_total_num_neighbors) {
					next_vertices_with_distances.push(std::make_pair(neighbor_index, distance));
				} else {
					if (distance <= max_influence) {
						next_vertices_with_distances.push(std::make_pair(neighbor_index, distance));
					}
				}
			}
		}

		// Store the nearest neighbors.
		int neighbor_count = neighbor_node_indices.size();

		float weight_sum = 0.f;
		for (int neighbor_index = 0; neighbor_index < neighbor_count; neighbor_index++) {
			*graph_edges_out.mutable_data(node_index, neighbor_index) = neighbor_node_indices[neighbor_index];
			weight_sum += neighbor_node_weights[neighbor_index];
		}

		// Normalize weights
		if (weight_sum > 0) {
			for (int neighbor_index = 0; neighbor_index < neighbor_count; neighbor_index++) {
				*graph_edges_weights_out.mutable_data(node_index, neighbor_index) = neighbor_node_weights[neighbor_index] / weight_sum;
			}
		} else if (neighbor_count > 0) {
			for (int neighbor_index = 0; neighbor_index < neighbor_count; neighbor_index++) {
				*graph_edges_weights_out.mutable_data(node_index, neighbor_index) = neighbor_node_weights[neighbor_index] / neighbor_count;
			}
		}

		// Store edge distance.
		for (int neighbor_index = 0; neighbor_index < neighbor_count; neighbor_index++) {
			*graph_edges_distances_out.mutable_data(node_index, neighbor_index) = neighbor_node_distances[neighbor_index];
		}
	}
}

/**
 * \brief Compute the graph edges between nodes, connecting nearest nodes using geodesic distances.
 * \param vertex_positions_in
 * \param vertex_mask_in
 * \param face_indices_in
 * \param node_indices_in
 * \param max_neighbor_count
 * \param node_coverage
 * \param graph_edges_out
 * \param graph_edges_weights_out
 * \param graph_edges_distances_out
 * \param node_to_vertex_distances_out
 * \param enforce_total_num_neighbors
 */
void compute_edges_shortest_path(
		const py::array_t<float>& vertex_positions_in,
		const py::array_t<bool>& vertex_mask_in,
		const py::array_t<int>& face_indices_in,
		const py::array_t<int>& node_indices_in,
		int max_neighbor_count, const float node_coverage,
		py::array_t<int>& graph_edges_out,
		py::array_t<float>& graph_edges_weights_out,
		py::array_t<float>& graph_edges_distances_out,
		py::array_t<float>& node_to_vertex_distances_out,
		bool enforce_total_num_neighbors
) {
	compute_edges_shortest_path_generic(vertex_positions_in, face_indices_in, node_indices_in, max_neighbor_count, node_coverage,
	                                    graph_edges_out, graph_edges_weights_out, graph_edges_distances_out, node_to_vertex_distances_out,
	                                    enforce_total_num_neighbors,
	                                    [&vertex_mask_in](int vertex_index) { return *vertex_mask_in.data(vertex_index); });
}

void compute_edges_shortest_path(
		const py::array_t<float>& vertex_positions_in,
		const py::array_t<int>& face_indices_in,
		const py::array_t<int>& node_indices_in,
		int max_neighbor_count, float node_coverage,
		py::array_t<int>& graph_edges_out,
		py::array_t<float>& graph_edges_weights_out,
		py::array_t<float>& graph_edges_distances_out,
		py::array_t<float>& node_to_vertex_distances_out,
		bool enforce_total_num_neighbors
) {
	compute_edges_shortest_path_generic(vertex_positions_in, face_indices_in, node_indices_in, max_neighbor_count, node_coverage,
	                                    graph_edges_out, graph_edges_weights_out, graph_edges_distances_out, node_to_vertex_distances_out,
	                                    enforce_total_num_neighbors,
	                                    [](int vertex_index) { return true; });
}

py::tuple compute_edges_shortest_path(
		const py::array_t<float>& vertex_positions_in,
		const py::array_t<bool>& vertex_mask_in,
		const py::array_t<int>& face_indices_in,
		const py::array_t<int>& node_indices_in,
		int max_neighbor_count,
		float node_coverage,
		bool enforce_total_num_neighbors) {

	const ssize_t node_count = node_indices_in.shape(0);
	const ssize_t vertex_count = vertex_positions_in.shape(0);

	py::array_t<int> graph_edges_out({node_count, static_cast<ssize_t>(max_neighbor_count)});
	std::fill_n(graph_edges_out.mutable_data(0, 0), graph_edges_out.size(), -1);

	py::array_t<float> graph_edges_weights_out({node_count, static_cast<ssize_t>(max_neighbor_count)});
	memset(graph_edges_weights_out.mutable_data(0, 0), 0, graph_edges_weights_out.size() * sizeof(float));

	py::array_t<float> graph_edges_distances_out({node_count, static_cast<ssize_t>(max_neighbor_count)});
	memset(graph_edges_distances_out.mutable_data(0, 0), 0, graph_edges_distances_out.size() * sizeof(float));

	py::array_t<float> node_to_vertex_distances_out({node_count, vertex_count});
	std::fill_n(node_to_vertex_distances_out.mutable_data(0, 0), node_to_vertex_distances_out.size(), -1.f);

	compute_edges_shortest_path(vertex_positions_in, vertex_mask_in, face_indices_in, node_indices_in, max_neighbor_count, node_coverage,
	                            graph_edges_out, graph_edges_weights_out, graph_edges_distances_out, node_to_vertex_distances_out,
	                            enforce_total_num_neighbors);

	return py::make_tuple(graph_edges_out, graph_edges_weights_out, graph_edges_distances_out, node_to_vertex_distances_out);
}

py::tuple compute_edges_shortest_path(
		const py::array_t<float>& vertex_positions_in,
		const py::array_t<int>& face_indices_in,
		const py::array_t<int>& node_indices_in,
		int max_neighbor_count, float node_coverage,
		bool enforce_total_num_neighbors
) {
	const ssize_t node_count = node_indices_in.shape(0);
	const ssize_t vertex_count = vertex_positions_in.shape(0);

	py::array_t<int> graph_edges_out({node_count, static_cast<ssize_t>(max_neighbor_count)});
	std::fill_n(graph_edges_out.mutable_data(0, 0), graph_edges_out.size(), -1);

	py::array_t<float> graph_edges_weights_out({node_count, static_cast<ssize_t>(max_neighbor_count)});
	memset(graph_edges_weights_out.mutable_data(0, 0), 0, graph_edges_weights_out.size() * sizeof(float));

	py::array_t<float> graph_edges_distances_out({node_count, static_cast<ssize_t>(max_neighbor_count)});
	memset(graph_edges_distances_out.mutable_data(0, 0), 0, graph_edges_distances_out.size() * sizeof(float));

	py::array_t<float> node_to_vertex_distances_out({node_count, vertex_count});
	std::fill_n(node_to_vertex_distances_out.mutable_data(0, 0), node_to_vertex_distances_out.size(), -1.f);

	compute_edges_shortest_path(vertex_positions_in, face_indices_in, node_indices_in, max_neighbor_count, node_coverage,
	                            graph_edges_out, graph_edges_weights_out, graph_edges_distances_out, node_to_vertex_distances_out,
	                            enforce_total_num_neighbors);

	return py::make_tuple(graph_edges_out, graph_edges_weights_out, graph_edges_distances_out, node_to_vertex_distances_out);
}

py::array_t<int> compute_edges_euclidean(const py::array_t<float>& node_positions, int max_neighbor_count) {
	int node_count = node_positions.shape(0);

	py::array_t<int> graph_edges = py::array_t<int>({node_count, max_neighbor_count});

	// Find nearest Euclidean neighbors for each node.
	for (int source_node_index = 0; source_node_index < node_count; source_node_index++) {
		Eigen::Vector3f node_position(*node_positions.data(source_node_index, 0),
		                              *node_positions.data(source_node_index, 1),
		                              *node_positions.data(source_node_index, 2));

		// Keep only the k nearest Euclidean neighbors.
		std::list<std::pair<int, float>> nearest_nodes_with_squared_distances;

		for (int neighbor_index = 0; neighbor_index < node_count; neighbor_index++) {
			if (neighbor_index == source_node_index) continue;

			Eigen::Vector3f neighbor_position(*node_positions.data(neighbor_index, 0), *node_positions.data(neighbor_index, 1),
			                                  *node_positions.data(neighbor_index, 2));

			float squared_distance = (node_position - neighbor_position).squaredNorm();
			bool neighbor_inserted = false;
			for (auto it = nearest_nodes_with_squared_distances.begin(); it != nearest_nodes_with_squared_distances.end(); ++it) {
				// We insert the element at the first position where its distance is smaller than the other
				// element's distance, which enables us to always keep a sorted list of at most k nearest
				// neighbors.
				if (squared_distance <= it->second) {
					it = nearest_nodes_with_squared_distances.insert(it, std::make_pair(neighbor_index, squared_distance));
					neighbor_inserted = true;
					break;
				}
			}

			if (!neighbor_inserted && nearest_nodes_with_squared_distances.size() < max_neighbor_count) {
				nearest_nodes_with_squared_distances.emplace_back(std::make_pair(neighbor_index, squared_distance));
			}

			// We keep only the list of k nearest elements.
			if (neighbor_inserted && nearest_nodes_with_squared_distances.size() > max_neighbor_count) {
				nearest_nodes_with_squared_distances.pop_back();
			}
		}

		// Store nearest neighbor indices.
		int neighbor_index = 0;
		for (auto& nearest_nodes_with_squared_distance : nearest_nodes_with_squared_distances) {
			int destination_node_index = nearest_nodes_with_squared_distance.first;
			*graph_edges.mutable_data(source_node_index, neighbor_index) = destination_node_index;
			neighbor_index++;
		}

		for (neighbor_index = nearest_nodes_with_squared_distances.size(); neighbor_index < max_neighbor_count; neighbor_index++) {
			*graph_edges.mutable_data(source_node_index, neighbor_index) = -1;
		}
	}

	return graph_edges;
}

inline int traverse_neighbors(const std::vector<std::set<int>>& node_neighbors, std::vector<int>& cluster_ids, int cluster_id, int node_id) {
	if (cluster_ids[node_id] != -1) return 0;

	std::set<int> active_node_indices;

	// Initialize with current node.
	int cluster_size = 0;
	active_node_indices.insert(node_id);

	// Process until we have no active nodes anymore.
	while (!active_node_indices.empty()) {
		int active_node_id = *active_node_indices.begin();
		active_node_indices.erase(active_node_indices.begin());

		if (cluster_ids[active_node_id] == -1) {
			cluster_ids[active_node_id] = cluster_id;
			++cluster_size;
		}

		// Look if we need to process any of the neighbors
		for (const auto& n_idx : node_neighbors[active_node_id]) {
			if (cluster_ids[n_idx] == -1) {    // If it doesn't have a cluster yet
				active_node_indices.insert(n_idx);
			}
		}
	}

	return cluster_size;
}

void node_and_edge_clean_up(const py::array_t<int>& graph_edges, py::array_t<bool>& valid_nodes_mask) {
	int num_nodes = graph_edges.shape(0);
	int max_num_neighbors = graph_edges.shape(1);

	std::list<int> removed_nodes;

	while (true) {
		int num_newly_removed_nodes = 0;

		for (int node_id = 0; node_id < num_nodes; ++node_id) {

			if (*valid_nodes_mask.data(node_id, 0) == false) {
				// if node has been already removed, continue
				continue;
			}

			int num_neighbors = 0;
			for (int i = 0; i < max_num_neighbors; ++i) {

				int neighbor_id = *graph_edges.data(node_id, i);

				// if neighboring node is -1, break, since by design 'graph_edges' has
				// the shape [2, 3, 6, -1, -1, -1, -1, -1]
				if (neighbor_id == -1) {
					break;
				}

				// if neighboring node has been marked as invalid, continue
				if (std::find(removed_nodes.begin(), removed_nodes.end(), neighbor_id) != removed_nodes.end()) {
					continue;
				}

				++num_neighbors;
			}
			// TODO: should be a parameter
			if (num_neighbors <= 1) {
				// remove node
				*valid_nodes_mask.mutable_data(node_id, 0) = false;
				removed_nodes.emplace_back(node_id);
				// std::cout << "\tremoving node_id " << node_id << std::endl;
				++num_newly_removed_nodes;
			}
		}

		// std::cout << "num_newly_removed_nodes: " << num_newly_removed_nodes << std::endl;

		if (num_newly_removed_nodes == 0) {
			break;
		}
	}
}

py::tuple compute_clusters(
		const py::array_t<int>& graph_edges_in
) {

	int node_count = graph_edges_in.shape(0);
	int max_neighbor_count = graph_edges_in.shape(1);

	py::array_t<int> graph_clusters_out({static_cast<ssize_t>(node_count), static_cast<ssize_t>(1)});
	std::fill_n(graph_clusters_out.mutable_data(0, 0), graph_clusters_out.size(), -1);

	// convert graph_edges to a vector of sets
	std::vector<std::set<int>> node_neighbors(node_count);

	for (int node_id = 0; node_id < node_count; ++node_id) {
		for (int neighbor_idx = 0; neighbor_idx < max_neighbor_count; ++neighbor_idx) {

			int neighbor_id = *graph_edges_in.data(node_id, neighbor_idx);

			if (neighbor_id == -1) {
				break;
			}

			node_neighbors[node_id].insert(neighbor_id);
			node_neighbors[neighbor_id].insert(node_id);
		}
	}

	std::vector<int> cluster_ids(node_count, -1);
	std::vector<int> cluster_sizes_out;

	int cluster_id = 0;
	for (int node_id = 0; node_id < node_count; ++node_id) {
		int cluster_size = traverse_neighbors(node_neighbors, cluster_ids, cluster_id, node_id);
		if (cluster_size > 0) {
			cluster_id++;
			cluster_sizes_out.push_back(cluster_size);
		}
	}

	for (int node_id = 0; node_id < node_count; ++node_id) {
		*graph_clusters_out.mutable_data(node_id, 0) = cluster_ids[node_id];
	}

	return py::make_tuple(cluster_sizes_out, graph_clusters_out);
}

inline void find_nearest_nodes(
		const py::array_t<float>& node_to_vertex_distance,
		const py::array_t<int>& valid_nodes_mask,
		const int vertex_id,
		std::vector<int>& nearest_node_indexes,
		std::vector<float>& nearest_node_distances
) {
	int num_nodes = node_to_vertex_distance.shape(0);

	std::map<int, float> node_map;

	for (int node_index = 0; node_index < num_nodes; ++node_index) {

		// discard node if it was marked as invalid (due to not having enough neighbors)
		if (*valid_nodes_mask.data(node_index, 0) == false) {
			continue;
		}

		float dist = *node_to_vertex_distance.data(node_index, vertex_id);

		if (dist >= 0) {
			node_map.emplace(node_index, dist);
		}
	}

	// Sort the map by distance
	// Declaring the type of Predicate that accepts 2 pairs and return a bool
	typedef std::function<bool(std::pair<int, float>, std::pair<int, float>)> IntFloatPairComparator;

	// Defining a lambda function to compare two pairs. It will compare two pairs using second field
	IntFloatPairComparator compare_by_float =
			[](std::pair<int, float> node1, std::pair<int, float> node2) {
				return node1.second < node2.second;
			};

	// Declaring a set that will store the pairs using above comparison logic
	std::set<std::pair<int, float>, IntFloatPairComparator> node_set(
			node_map.begin(), node_map.end(), compare_by_float
	);

	for (auto node : node_set) {
		nearest_node_indexes.push_back(node.first);
		nearest_node_distances.push_back(node.second);

		if (nearest_node_indexes.size() == GRAPH_K) {
			break;
		}
	}
}

void compute_pixel_anchors_shortest_path(
		const py::array_t<float>& node_to_vertex_distance,
		const py::array_t<int>& valid_nodes_mask,
		const py::array_t<float>& vertices,
		const py::array_t<int>& vertex_pixels,
		py::array_t<int>& pixel_anchors,
		py::array_t<float>& pixel_weights,
		int width, int height,
		float node_coverage
) {
	// Allocate graph node ids and corresponding skinning weights.
	// Initialize with invalid anchors.
	pixel_anchors.resize({height, width, GRAPH_K}, false);
	pixel_weights.resize({height, width, GRAPH_K}, false);

	for (int y = 0; y < height; y++) {
		for (int x = 0; x < width; x++) {
			for (int k = 0; k < GRAPH_K; k++) {
				*pixel_anchors.mutable_data(y, x, k) = -1;
				*pixel_weights.mutable_data(y, x, k) = 0.f;
			}
		}
	}

	const int vertex_count = vertices.shape(0);

#pragma omp parallel for default(none) shared(vertex_pixels, pixel_anchors, pixel_weights, valid_nodes_mask, \
    node_to_vertex_distance) firstprivate(vertex_count, node_coverage)
	for (int vertex_id = 0; vertex_id < vertex_count; vertex_id++) {
		// Get corresponding pixel location
		int u = *vertex_pixels.data(vertex_id, 0);
		int v = *vertex_pixels.data(vertex_id, 1);

		// Initialize some variables
		std::vector<int> nearest_node_ids;
		std::vector<float> distances_to_nearest_nodes;
		std::vector<float> skinning_weights;

		nearest_node_ids.reserve(GRAPH_K);
		distances_to_nearest_nodes.reserve(GRAPH_K);
		skinning_weights.reserve(GRAPH_K);

		// Find closest nodes
		find_nearest_nodes(
				node_to_vertex_distance, valid_nodes_mask, vertex_id,
				nearest_node_ids, distances_to_nearest_nodes
		);

		int anchor_count = nearest_node_ids.size();

		// Compute skinning weights.
		float weight_sum{0.f};
		for (int i = 0; i < anchor_count; ++i) {
			float shortest_path_distance_to_node = distances_to_nearest_nodes[i];

			float weight = compute_anchor_weight(shortest_path_distance_to_node, node_coverage);
			weight_sum += weight;

			skinning_weights.push_back(weight);
		}

		// Normalize the skinning weights.
		if (weight_sum > 0) {
			for (int anchor_index = 0; anchor_index < anchor_count; anchor_index++) {
				skinning_weights[anchor_index] /= weight_sum;
			}
		} else if (anchor_count > 0) {
			for (int anchor_index = 0; anchor_index < anchor_count; anchor_index++) {
				skinning_weights[anchor_index] = 1.f / static_cast<float>(anchor_count);
			}
		}

		// Store the results.
		for (int i = 0; i < anchor_count; i++) {
			*pixel_anchors.mutable_data(v, u, i) = nearest_node_ids[i];
			*pixel_weights.mutable_data(v, u, i) = skinning_weights[i];
		}
	}
}

py::tuple compute_pixel_anchors_shortest_path(
		const py::array_t<float>& node_to_vertex_distance,
		const py::array_t<int>& valid_nodes_mask,
		const py::array_t<float>& vertices,
		const py::array_t<int>& vertex_pixels,
		int width, int height,
		float node_coverage) {
	py::array_t<int> pixel_anchors;
	py::array_t<float> pixel_weights;
	compute_pixel_anchors_shortest_path(node_to_vertex_distance, valid_nodes_mask, vertices, vertex_pixels,
	                                    pixel_anchors, pixel_weights, width, height, node_coverage);
	return py::make_tuple(pixel_anchors, pixel_weights);
}

/**
 * \brief For each input pixel it computes 4 nearest anchors, using Euclidean distances.
 * It also compute skinning weights for every pixel.
 * \param graph_nodes N X 3 graph node positions
 * \param point_image height x width x 3 point cloud
 * \param node_coverage the maximal distance of each point in the point cloud to the nearest node used to generate the graph
 * \param pixel_anchors [out] height x width x 4 array, where each 2d sampling with four residuals holds the indices of the nodes
 * chosen as (controlling) anchor points for the pixel at the same 2d location
 * \param pixel_weights [out] height x width x 4 array, where each 2d sampling with four residuals holds the
 * influence weights of the nodes assigned to the pixel in pixel_anchors
 */
void compute_pixel_anchors_euclidean(
		const py::array_t<float>& graph_nodes,
		const py::array_t<float>& point_image,
		const float node_coverage,
		py::array_t<int>& pixel_anchors,
		py::array_t<float>& pixel_weights
) {
	const int node_count = graph_nodes.shape(0);
	const int width = point_image.shape(1);
	const int height = point_image.shape(0);

	// Allocate graph node ids and corresponding skinning weights.
	// Initialize with invalid anchors.
	pixel_anchors.resize({height, width, GRAPH_K}, false);
	pixel_weights.resize({height, width, GRAPH_K}, false);

	for (int y = 0; y < height; y++) {
		for (int x = 0; x < width; x++) {
			for (int k = 0; k < GRAPH_K; k++) {
				*pixel_anchors.mutable_data(y, x, k) = -1;
				*pixel_weights.mutable_data(y, x, k) = 0.f;
			}
		}
	}

	// Compute anchors for every pixel.
#pragma omp parallel for default(none) shared(point_image, pixel_anchors, pixel_weights, graph_nodes), firstprivate(height, width, node_count, node_coverage)
	for (int y = 0; y < height; y++) {
		for (int x = 0; x < width; x++) {
			// Query 3d pixel position.
			Eigen::Vector3f pixel_position(*point_image.data(y, x, 0),
			                               *point_image.data(y, x, 1),
			                               *point_image.data(y, x, 2));
			if (pixel_position.z() <= 0) continue;

			// Keep only the k nearest Euclidean neighbors.
			std::list<std::pair<int, float>> nearest_nodes_with_squared_distances;

			for (int node_index = 0; node_index < node_count; node_index++) {
				Eigen::Vector3f node_position(*graph_nodes.data(node_index, 0),
				                              *graph_nodes.data(node_index, 1),
				                              *graph_nodes.data(node_index, 2));

				float squared_distance = (pixel_position - node_position).squaredNorm();
				bool nodes_inserted = false;
				for (auto it = nearest_nodes_with_squared_distances.begin(); it != nearest_nodes_with_squared_distances.end(); ++it) {
					// We insert the element at the first position where its distance is smaller than the other
					// element's distance, which enables us to always keep a sorted list of at most k nearest
					// neighbors.
					if (squared_distance <= it->second) {
						it = nearest_nodes_with_squared_distances.insert(it, std::make_pair(node_index, squared_distance));
						nodes_inserted = true;
						break;
					}
				}

				if (!nodes_inserted && nearest_nodes_with_squared_distances.size() < GRAPH_K) {
					nearest_nodes_with_squared_distances.emplace_back(std::make_pair(node_index, squared_distance));
				}

				// We keep only the list of k nearest elements.
				if (nodes_inserted && nearest_nodes_with_squared_distances.size() > GRAPH_K) {
					nearest_nodes_with_squared_distances.pop_back();
				}
			}

			// Compute skinning weights.
			std::vector<int> nearest_euclidean_node_indices;
			nearest_euclidean_node_indices.reserve(nearest_nodes_with_squared_distances.size());

			std::vector<float> skinning_weights;
			skinning_weights.reserve(nearest_nodes_with_squared_distances.size());

			float weight_sum{0.f};
			for (auto& nearest_nodes_with_squared_distance : nearest_nodes_with_squared_distances) {
				int node_index = nearest_nodes_with_squared_distance.first;

				Eigen::Vector3f node_position(*graph_nodes.data(node_index, 0),
				                              *graph_nodes.data(node_index, 1),
				                              *graph_nodes.data(node_index, 2));
				float weight = compute_anchor_weight(pixel_position, node_position, node_coverage);
				weight_sum += weight;

				nearest_euclidean_node_indices.push_back(node_index);
				skinning_weights.push_back(weight);
			}

			// Normalize the skinning weights.
			int anchor_count = static_cast<int>(nearest_euclidean_node_indices.size());

			if (weight_sum > 0) {
				for (int i = 0; i < anchor_count; i++) skinning_weights[i] /= weight_sum;
			} else if (anchor_count > 0) {
				for (int i = 0; i < anchor_count; i++) skinning_weights[i] = 1.f / static_cast<float>(anchor_count);
			}

			// Store the results.
			for (int i = 0; i < anchor_count; i++) {
				*pixel_anchors.mutable_data(y, x, i) = nearest_euclidean_node_indices[i];
				*pixel_weights.mutable_data(y, x, i) = skinning_weights[i];
			}
		}
	}
}

// All-pass filter for easy decltype and thread-friendly default arguments
const auto TruePredicate = [](int x) { return true; };

inline decltype(TruePredicate) GetTruePredicate() {
	return TruePredicate;
}

template<typename Predicate = decltype(TruePredicate)>
inline std::list<std::pair<int, float>>
get_knn_with_squared_distances_brute_force(const Eigen::Vector3f& vertex, const py::array_t<float>& graph_nodes,
                                           const int neighbor_count, Predicate filter = GetTruePredicate()) {
// Keep only the k nearest Euclidean neighbors.
	std::list<std::pair<int, float>> nearest_nodes_with_squared_distances;

	for (int node_index = 0; node_index < graph_nodes.shape(0); node_index++) {
		if (!filter(node_index))
			continue;

		Eigen::Vector3f node_position(*graph_nodes.data(node_index, 0),
		                              *graph_nodes.data(node_index, 1),
		                              *graph_nodes.data(node_index, 2));

		float squared_distance = (vertex - node_position).squaredNorm();
		bool nodes_inserted = false;
		for (auto it = nearest_nodes_with_squared_distances.begin(); it != nearest_nodes_with_squared_distances.end(); ++it) {
			// We insert the element at the first position where its distance is smaller than the other
			// element's distance, which enables us to always keep a sorted list of at most k nearest
			// neighbors.
			if (squared_distance <= it->second) {
				it = nearest_nodes_with_squared_distances.insert(it, std::make_pair(node_index, squared_distance));
				nodes_inserted = true;
				break;
			}
		}

		if (!nodes_inserted && nearest_nodes_with_squared_distances.size() < neighbor_count) {
			nearest_nodes_with_squared_distances.emplace_back(std::make_pair(node_index, squared_distance));
		}

		// We keep only the list of k nearest elements.
		if (nodes_inserted && nearest_nodes_with_squared_distances.size() > neighbor_count) {
			nearest_nodes_with_squared_distances.pop_back();
		}
	}
	return nearest_nodes_with_squared_distances;
}

py::tuple compute_pixel_anchors_euclidean(
		const py::array_t<float>& graph_nodes,
		const py::array_t<float>& point_image,
		float node_coverage
) {
	py::array_t<int> pixel_anchors;
	py::array_t<float> pixel_weights;

	compute_pixel_anchors_euclidean(graph_nodes, point_image, node_coverage, pixel_anchors, pixel_weights);

	return py::make_tuple(pixel_anchors, pixel_weights);
}

py::tuple compute_vertex_anchors_euclidean(
		const py::array_t<float>& graph_nodes,
		const py::array_t<float>& vertices,
		float node_coverage
) {
	py::array_t<int> vertex_anchors({vertices.shape(0), static_cast<ssize_t>(GRAPH_K)});
	py::array_t<float> vertex_weights({vertices.shape(0), static_cast<ssize_t>(GRAPH_K)});
	std::fill_n(vertex_anchors.mutable_data(0), vertex_anchors.size(), -1);
	std::fill_n(vertex_weights.mutable_data(0), vertex_weights.size(), 0.f);

#pragma omp parallel for default(none) shared(vertices, graph_nodes, vertex_anchors, vertex_weights), firstprivate(node_coverage)
	for (int i_vertex = 0; i_vertex < vertices.shape(0); i_vertex++) {
		// Query 3d pixel position.
		Eigen::Vector3f vertex(*vertices.data(i_vertex, 0),
		                       *vertices.data(i_vertex, 1),
		                       *vertices.data(i_vertex, 2));

		// Keep only the k nearest Euclidean neighbors.
		std::list<std::pair<int, float>> nearest_nodes_with_squared_distances =
				get_knn_with_squared_distances_brute_force(vertex, graph_nodes, GRAPH_K);

		// Compute skinning weights.
		std::vector<int> nearest_euclidean_node_indices;
		nearest_euclidean_node_indices.reserve(nearest_nodes_with_squared_distances.size());

		std::vector<float> skinning_weights;
		skinning_weights.reserve(nearest_nodes_with_squared_distances.size());

		float weight_sum{0.f};
		for (auto& nearest_nodes_with_squared_distance : nearest_nodes_with_squared_distances) {
			int node_index = nearest_nodes_with_squared_distance.first;

			Eigen::Vector3f node_position(*graph_nodes.data(node_index, 0),
			                              *graph_nodes.data(node_index, 1),
			                              *graph_nodes.data(node_index, 2));
			float weight = compute_anchor_weight(vertex, node_position, node_coverage);
			weight_sum += weight;

			nearest_euclidean_node_indices.push_back(node_index);
			skinning_weights.push_back(weight);
		}

		// Normalize the skinning weights.
		int anchor_count = static_cast<int>(nearest_euclidean_node_indices.size());

		if (weight_sum > 0) {
			for (int i = 0; i < anchor_count; i++) skinning_weights[i] /= weight_sum;
		} else if (anchor_count > 0) {
			for (int i = 0; i < anchor_count; i++) skinning_weights[i] = 1.f / static_cast<float>(anchor_count);
		}

		// Store the results.
		for (int i_anchor_node = 0; i_anchor_node < anchor_count; i_anchor_node++) {
			*vertex_anchors.mutable_data(i_vertex, i_anchor_node) = nearest_euclidean_node_indices[i_anchor_node];
			*vertex_weights.mutable_data(i_vertex, i_anchor_node) = skinning_weights[i_anchor_node];
		}
	}
	return py::make_tuple(vertex_anchors, vertex_weights);
}

/**
  * Sample graph regularly from the image, using pixel-wise connectivity
  * (connecting each pixel with at most 8 neighbors).
  */
void construct_regular_graph(
		const py::array_t<float>& point_image,
		int x_nodes, int y_nodes,
		float edge_threshold,
		float max_point_to_node_distance,
		float max_depth,
		py::array_t<float>& graph_nodes,
		py::array_t<int>& graph_edges,
		py::array_t<int>& pixel_anchors,
		py::array_t<float>& pixel_weights
) {
	int width = point_image.shape(1);
	int height = point_image.shape(0);

	float x_step = float(width - 1) / static_cast<float>(x_nodes - 1);
	float y_step = float(height - 1) / static_cast<float>(y_nodes - 1);

	// Sample graph nodes.
	// We need to maintain the mapping from all -> valid nodes ids.
	int node_count = x_nodes * y_nodes;
	std::vector<int> sampled_node_mapping(node_count, -1);

	std::vector<Eigen::Vector3f> node_positions;
	node_positions.reserve(node_count);

	int sampled_node_count = 0;
	for (int y = 0; y < y_nodes; y++) {
		for (int x = 0; x < x_nodes; x++) {
			int linear_node_index = y * x_nodes + x;

			// We use nearest neighbor interpolation for node position
			// computation.
			int x_pixel = static_cast<int>(std::round(static_cast<float>(x) * x_step));
			int y_pixel = static_cast<int>(std::round(static_cast<float>(y) * y_step));

			Eigen::Vector3f pixel_position(*point_image.data(y_pixel, x_pixel, 0),
			                               *point_image.data(y_pixel, x_pixel, 1),
			                               *point_image.data(y_pixel, x_pixel, 2));
			if (pixel_position.z() <= 0 || pixel_position.z() > max_depth) continue;

			node_positions.push_back(pixel_position);
			sampled_node_mapping[linear_node_index] = sampled_node_count;
			sampled_node_count++;
		}
	}

	// Compute graph edges using pixel-wise connectivity. Each node
	// is connected with at most 8 neighboring pixels.
	int max_neighbor_count = 8;
	float edge_threshold_squared = edge_threshold * edge_threshold;

	std::vector<int> sampled_node_edges(sampled_node_count * max_neighbor_count, -1);
	std::vector<bool> connected_nodes(sampled_node_count, false);

	int connected_node_count = 0;
	for (int y = 0; y < y_nodes; y++) {
		for (int x = 0; x < x_nodes; x++) {
			int node_index = y * x_nodes + x;
			int sampled_node_index = sampled_node_mapping[node_index];

			if (sampled_node_index >= 0) {
				Eigen::Vector3f node_position = node_positions[sampled_node_index];

				int neighbor_count = 0;
				for (int y_delta = -1; y_delta <= 1; y_delta++) {
					for (int x_delta = -1; x_delta <= 1; x_delta++) {
						int x_neighbor = x + x_delta;
						int y_neighbor = y + y_delta;
						if (x_neighbor < 0 || x_neighbor >= x_nodes || y_neighbor < 0 || y_neighbor >= y_nodes)
							continue;

						int neighbor_index = y_neighbor * x_nodes + x_neighbor;

						if (neighbor_index == node_index || neighbor_index < 0)
							continue;

						int sampled_neighbor_index = sampled_node_mapping[neighbor_index];
						if (sampled_neighbor_index >= 0) {
							Eigen::Vector3f neighbor_position = node_positions[sampled_neighbor_index];

							if ((neighbor_position - node_position).squaredNorm() <= edge_threshold_squared) {
								sampled_node_edges[sampled_node_index * max_neighbor_count + neighbor_count] = sampled_neighbor_index;
								neighbor_count++;
							}
						}
					}
				}

				for (int i = neighbor_count; i < max_neighbor_count; i++) {
					sampled_node_edges[sampled_node_index * max_neighbor_count + i] = -1;
				}

				if (neighbor_count > 0) {
					connected_nodes[sampled_node_index] = true;
					connected_node_count += 1;
				}
			}
		}
	}

	// Filter out canonical_node_positions with no edges.
	// After changing node ids the edge ids need to be changed as well.
	std::vector<int> valid_node_mapping(sampled_node_count, -1);
	{
		graph_nodes.resize({connected_node_count, 3}, false);
		graph_edges.resize({connected_node_count, max_neighbor_count}, false);

		int valid_node_index = 0;
		for (int y = 0; y < y_nodes; y++) {
			for (int x = 0; x < x_nodes; x++) {
				int node_index = y * x_nodes + x;
				int sampled_node_index = sampled_node_mapping[node_index];

				if (sampled_node_index >= 0 && connected_nodes[sampled_node_index]) {
					valid_node_mapping[sampled_node_index] = valid_node_index;

					Eigen::Vector3f node_position = node_positions[sampled_node_index];
					*graph_nodes.mutable_data(valid_node_index, 0) = node_position.x();
					*graph_nodes.mutable_data(valid_node_index, 1) = node_position.y();
					*graph_nodes.mutable_data(valid_node_index, 2) = node_position.z();

					valid_node_index++;
				}
			}
		}
	}

	for (int y = 0; y < y_nodes; y++) {
		for (int x = 0; x < x_nodes; x++) {
			int node_index = y * x_nodes + x;
			int sampled_node_index = sampled_node_mapping[node_index];

			if (sampled_node_index >= 0 && connected_nodes[sampled_node_index]) {
				int valid_node_index = valid_node_mapping[sampled_node_index];

				if (valid_node_index >= 0) {
					for (int i = 0; i < max_neighbor_count; i++) {
						int sampledNeighborId = sampled_node_edges[sampled_node_index * max_neighbor_count + i];
						if (sampledNeighborId >= 0) {
							*graph_edges.mutable_data(valid_node_index, i) = valid_node_mapping[sampledNeighborId];
						} else {
							*graph_edges.mutable_data(valid_node_index, i) = -1;
						}
					}
				}
			}
		}
	}

	// Compute pixel anchors and weights.
	pixel_anchors.resize({height, width, 4}, false);
	pixel_weights.resize({height, width, 4}, false);

	float max_point_to_node_distance_squared = max_point_to_node_distance * max_point_to_node_distance;

	for (int y = 0; y < height; y++) {
		for (int x = 0; x < width; x++) {
			// Initialize with invalid residuals.
			for (int k = 0; k < 4; k++) {
				*pixel_anchors.mutable_data(y, x, k) = -1;
				*pixel_weights.mutable_data(y, x, k) = 0.f;
			}

			// Compute 4 nearest canonical_node_positions.
			float x_node = float(x) / x_step;
			float y_node = float(y) / y_step;

			int x0 = std::floor(x_node), x1 = x0 + 1;
			int y0 = std::floor(y_node), y1 = y0 + 1;

			// Check that all neighboring nodes are valid.
			if (x0 < 0 || x1 >= x_nodes || y0 < 0 || y1 >= y_nodes)
				continue;

			int sampled_node_00 = sampled_node_mapping[y0 * x_nodes + x0];
			int sampled_node_01 = sampled_node_mapping[y1 * x_nodes + x0];
			int sampled_node_10 = sampled_node_mapping[y0 * x_nodes + x1];
			int sampled_node_11 = sampled_node_mapping[y1 * x_nodes + x1];

			if (sampled_node_00 < 0 || sampled_node_01 < 0 || sampled_node_10 < 0 || sampled_node_11 < 0)
				continue;

			int valid_node_00 = valid_node_mapping[sampled_node_00];
			int valid_node_01 = valid_node_mapping[sampled_node_01];
			int valid_node_10 = valid_node_mapping[sampled_node_10];
			int valid_node_11 = valid_node_mapping[sampled_node_11];

			if (valid_node_00 < 0 || valid_node_01 < 0 || valid_node_10 < 0 || valid_node_11 < 0)
				continue;

			// Check that all nodes are close enough to the point.
			Eigen::Vector3f pixel_position(*point_image.data(y, x, 0),
			                               *point_image.data(y, x, 1),
			                               *point_image.data(y, x, 2));
			if (pixel_position.z() <= 0 || pixel_position.z() > max_depth) continue;

			if ((pixel_position - node_positions[sampled_node_00]).squaredNorm() > max_point_to_node_distance_squared ||
			    (pixel_position - node_positions[sampled_node_01]).squaredNorm() > max_point_to_node_distance_squared ||
			    (pixel_position - node_positions[sampled_node_10]).squaredNorm() > max_point_to_node_distance_squared ||
			    (pixel_position - node_positions[sampled_node_11]).squaredNorm() > max_point_to_node_distance_squared) {
				continue;
			}

			// Compute bilinear weights.
			float dx = x_node - static_cast<float>(x0);
			float dy = y_node - static_cast<float>(y0);

			float w00 = (1 - dx) * (1 - dy);
			float w01 = (1 - dx) * dy;
			float w10 = dx * (1 - dy);
			float w11 = dx * dy;

			*pixel_anchors.mutable_data(y, x, 0) = valid_node_00;
			*pixel_weights.mutable_data(y, x, 0) = w00;
			*pixel_anchors.mutable_data(y, x, 1) = valid_node_01;
			*pixel_weights.mutable_data(y, x, 1) = w01;
			*pixel_anchors.mutable_data(y, x, 2) = valid_node_10;
			*pixel_weights.mutable_data(y, x, 2) = w10;
			*pixel_anchors.mutable_data(y, x, 3) = valid_node_11;
			*pixel_weights.mutable_data(y, x, 3) = w11;
		}
	}
}

py::tuple construct_regular_graph(
		const py::array_t<float>& point_image,
		int x_nodes, int y_nodes,
		float edge_threshold,
		float max_point_to_node_distance,
		float max_depth
) {
	py::array_t<float> graph_nodes;
	py::array_t<int> graph_edges;
	py::array_t<int> pixel_anchors;
	py::array_t<float> pixel_weights;
	construct_regular_graph(point_image, x_nodes, y_nodes, edge_threshold,
	                        max_point_to_node_distance, max_depth,
	                        graph_nodes, graph_edges, pixel_anchors, pixel_weights);
	return py::make_tuple(graph_nodes, graph_edges, pixel_anchors, pixel_weights);
}

void update_pixel_anchors(
		const std::map<int, int>& node_id_mapping,
		py::array_t<int>& pixel_anchors
) {
	int height = pixel_anchors.shape(0);
	int width = pixel_anchors.shape(1);
	int num_anchors = pixel_anchors.shape(2);

	for (int y = 0; y < height; y++) {
		for (int x = 0; x < width; x++) {

			for (int a = 0; a < num_anchors; a++) {

				int current_anchor_id = *pixel_anchors.data(y, x, a);

				if (current_anchor_id != -1) {
					int mapped_anchor_id = node_id_mapping.at(current_anchor_id);

					// update anchor only if it would actually change something
					if (mapped_anchor_id != current_anchor_id) {
						*pixel_anchors.mutable_data(y, x, a) = mapped_anchor_id;
					}
				}
			}
		}
	}
}

/**
 * \brief Find anchor_count shortest path neighbors to the provided vertex, i.e.
 * compute euclidean nearest neighbor node to the given vertex, then find n-1 nearest shortest-path neighbors of that node.
 * (To be moved) Record indices (of said neighbors) to provided anchor array & compute weights based on the corresponding shortest-path distances.
 * \param[out] vertex_anchors
 * \param[out] vertex_weights
 * \param vertex
 * \param nodes
 * \param edges
 * \param anchor_count
 * \param node_coverage
 */
inline
void djikstra_knn(int* vertex_anchors,
                  float* vertex_weights,
                  const Eigen::Vector3f& vertex,
                  const py::array_t<float>& nodes,
                  const py::array_t<int>& edges,
                  const int anchor_count,
                  const float node_coverage) {
	std::map<int, float> distance_by_shortest_path_anchor;
	typedef std::pair<int, float> index_and_distance;
	auto by_second_comparator = [](const index_and_distance& a, const index_and_distance& b) {
		return a.second > b.second;
	};
	std::priority_queue<index_and_distance, std::vector<index_and_distance>,
			decltype(by_second_comparator)> queue(by_second_comparator);
	std::vector<index_and_distance> discovered_anchors;

	while (distance_by_shortest_path_anchor.size() < anchor_count) {
		// Keep only the anchor_count nearest Euclidean neighbors.
		std::list<std::pair<int, float>> nearest_nodes_with_squared_distances =
				get_knn_with_squared_distances_brute_force(vertex, nodes, 1,
				                                           [distance_by_shortest_path_anchor](int i) {
					                                           return distance_by_shortest_path_anchor.find(i) ==
					                                                  distance_by_shortest_path_anchor.end();
				                                           });
		if (nearest_nodes_with_squared_distances.empty()) {
			break; // no node to initialize queue with, we've got no more valid anchors to consider
		}

		auto closest_euclidean_node_and_distance = nearest_nodes_with_squared_distances.front();
		queue.emplace(closest_euclidean_node_and_distance.first, sqrt(closest_euclidean_node_and_distance.second));

		while (!queue.empty()) {
			auto graph_node_and_distance = queue.top();
			queue.pop();

			int source_node_index = graph_node_and_distance.first;
			float source_path_distance = graph_node_and_distance.second;

			if (distance_by_shortest_path_anchor.find(source_node_index) != distance_by_shortest_path_anchor.end()) {
				continue;
			}

			// insert new node in the weight map & anchor set
			distance_by_shortest_path_anchor[source_node_index] = source_path_distance;
			discovered_anchors.push_back(graph_node_and_distance);
			if (distance_by_shortest_path_anchor.size() >= anchor_count)
				break;

			Eigen::Map<const Eigen::Vector3f> source_node(nodes.data(source_node_index, 0));

			for (int i_edge = 0; i_edge < GRAPH_K; i_edge++) {
				int target_node_index = *edges.data(source_node_index, i_edge);
				if (target_node_index > -1) {
					Eigen::Map<const Eigen::Vector3f> target_node(nodes.data(target_node_index, 0));
					float distance_source_to_target = (target_node - source_node).norm();
					queue.emplace(target_node_index, distance_source_to_target + source_path_distance);
				}
			}
		}
	}

	//TODO: DRY violation -- move weight normalization out to a separate function and reuse throughout file
	float weight_sum = 0;
	int valid_anchor_count = 0;
	for (const auto anchor : discovered_anchors) {
		vertex_anchors[valid_anchor_count] = anchor.first;
		float weight = compute_anchor_weight(anchor.second, node_coverage);
		vertex_weights[valid_anchor_count] = weight;
		weight_sum += weight;
		valid_anchor_count++;
	}

	// Normalize weights
	if (weight_sum > 0) {
		for (int i_anchor = 0; i_anchor < valid_anchor_count; i_anchor++) {
			vertex_weights[i_anchor] /= weight_sum;
		}
	} else if (valid_anchor_count > 0) {
		for (int i_anchor = 0; i_anchor < valid_anchor_count; i_anchor++) {
			vertex_weights[i_anchor] = 1.f / static_cast<float>(valid_anchor_count);
		}
	}
}

/**
 * \brief compute point anchors based on shortest path within the graph
 * \param vertex_to_node_distance
 * \param vertices
 * \param pixel_anchors
 * \param pixel_weights
 * \param width
 * \param height
 * \param node_coverage
 */
void compute_vertex_anchors_shortest_path(
		const py::array_t<float>& vertices,
		const py::array_t<float>& nodes,
		// edges are assumed to have been computed in a shortest-path manner,
		// i.e. what the NNRT authors called "geodesically" -- like via the
		// compute_edges_shortest_path function in the original source,
		// for each vertex - target vertices sharing an edge, ordered by shortest edge
		const py::array_t<int>& edges,
		py::array_t<int>& anchors,
		py::array_t<float>& weights,
		int anchor_count,
		float node_coverage) {
// Allocate graph node ids and corresponding skinning weights.
	int vertex_count = static_cast<int>(vertices.shape(0));
	int node_count = static_cast<int>(nodes.shape(0));

	assert(edges.shape(1) == GRAPH_K);
	assert(edges.shape(0) == node_count);

	// Initialize with invalid anchors.
	anchors.resize({vertex_count, anchor_count}, false);
	weights.resize({vertex_count, anchor_count}, false);
	std::fill_n(anchors.mutable_data(0, 0), anchors.size(), -1);
	memset(weights.mutable_data(0, 0), 0, weights.size() * sizeof(float));

#pragma omp parallel for default(none) shared(vertices, nodes, edges, anchors, weights)\
    firstprivate(vertex_count, node_coverage, anchor_count)
	for (int i_vertex = 0; i_vertex < vertex_count; i_vertex++) {
		// Query 3d pixel position.
		Eigen::Vector3f vertex(*vertices.data(i_vertex, 0),
		                       *vertices.data(i_vertex, 1),
		                       *vertices.data(i_vertex, 2));

		// Write to output directly
		auto vertex_anchors = anchors.mutable_data(i_vertex, 0);
		auto vertex_weights = weights.mutable_data(i_vertex, 0);

		djikstra_knn(vertex_anchors, vertex_weights, vertex, nodes, edges, anchor_count, node_coverage);
	} // end vertex loop
}

py::tuple compute_vertex_anchors_shortest_path(const py::array_t<float>& vertices, const py::array_t<float>& nodes, const py::array_t<int>& edges,
                                               int anchor_count, float node_coverage) {
	py::array_t<int> anchors;
	py::array_t<float> weights;
	compute_vertex_anchors_shortest_path(vertices, nodes, edges, anchors, weights, anchor_count, node_coverage);
	return py::make_tuple(anchors, weights);
}

void compute_pixel_anchors_shortest_path(const py::array_t<float>& point_image, const py::array_t<float>& nodes, const py::array_t<int>& edges,
                                         py::array_t<int>& anchors, py::array_t<float>& weights, int anchor_count,
                                         float node_coverage) {
	// Allocate graph node ids and corresponding skinning weights.
	const int image_height = static_cast<int>(point_image.shape(0));
	const int image_width = static_cast<int>(point_image.shape(1));
	const ssize_t image_pixel_count = point_image.shape(0) * point_image.shape(1);
	const int node_count = static_cast<int>(nodes.shape(0));

	assert(edges.shape(1) == GRAPH_K);
	assert(edges.shape(0) == node_count);
	assert(point_image.shape(2) == 3);
	assert(nodes.shape(1) == 3);

	// Initialize with invalid anchors.
	anchors.resize({image_height, image_width, anchor_count}, false);
	weights.resize({image_height, image_width, anchor_count}, false);
	std::fill_n(anchors.mutable_data(0, 0), anchors.size(), -1);
	memset(weights.mutable_data(0, 0), 0, weights.size() * sizeof(float));

#pragma omp parallel for default(none) shared(point_image, nodes, edges, anchors, weights)\
    firstprivate(image_height, image_width, node_coverage, anchor_count, image_pixel_count)
	for (ssize_t i_pixel = 0; i_pixel < image_pixel_count; i_pixel++) {
		int u = static_cast<int>(i_pixel % image_width);
		int v = static_cast<int>(i_pixel / image_width);

		// Query 3d pixel position.
		Eigen::Vector3f vertex(*point_image.data(v, u, 0),
		                       *point_image.data(v, u, 1),
		                       *point_image.data(v, u, 2));

		// Write to output directly
		auto vertex_anchors = anchors.mutable_data(v, u, 0);
		auto vertex_weights = weights.mutable_data(v, u, 0);

		djikstra_knn(vertex_anchors, vertex_weights, vertex, nodes, edges, anchor_count, node_coverage);
	} // end pixel loop
}

py::tuple compute_pixel_anchors_shortest_path(
		const py::array_t<float>& point_image, const py::array_t<float>& nodes,
		const py::array_t<int>& edges, int anchor_count, float node_coverage) {
	py::array_t<int> anchors;
	py::array_t<float> weights;
	compute_pixel_anchors_shortest_path(point_image, nodes, edges, anchors, weights, anchor_count, node_coverage);
	return py::make_tuple(anchors, weights);
}

} // namespace graph_proc

max_stars_repo_name: Algomorph/NeuralTracking

max_issues_repo_name: Algomorph/NeuralTracking

max_forks_repo_name: Algomorph/NeuralTracking


text:
// MIT License
//
// shrinkler-gba: Port of the Shrinkler Amiga executable cruncher for the GBA
// Copyright (c) 2021 Thomas Mathys
//
// 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 <boost/numeric/conversion/cast.hpp>
#include <boost/test/unit_test.hpp>
#include <cstddef>
#include <fstream>
#include "shrinklergba_unittest_config.hpp"
#include "test_utilities.hpp"

namespace shrinklergba_unittest
{

using std::filesystem::path;
using std::vector;

std::vector<unsigned char> load_binary_file(const std::filesystem::path& filename)
{
    path full_path = SHRINKLERGBA_UNITTEST_TESTDATA_DIRECTORY / filename;

    // Open file, stop it from eating whitespace.
    std::ifstream file(full_path, std::ios::binary);
    BOOST_REQUIRE_MESSAGE(file, "Could not open " + full_path.string());
    file.unsetf(std::ios::skipws);

    // Create vector with sufficient capacity to hold entire file.
    auto filesize = std::filesystem::file_size(full_path);
    vector<unsigned char> data;
    data.reserve(boost::numeric_cast<size_t>(filesize));

    // Read entire file
    data.insert(
        data.begin(),
        std::istream_iterator<unsigned char>(file),
        std::istream_iterator<unsigned char>());
    BOOST_REQUIRE_MESSAGE(file || file.eof(), "Error reading " + full_path.string());
    file.close();

    return data;
}

const std::vector<unsigned char> make_bytevector(const char* s)
{
    return std::vector<unsigned char>(s, s + std::strlen(s));
}

}


text:
#include "annotate.hpp"

#include "CgStr.hpp"
#include "MultiTuProcessor.hpp"
#include "cgWrappers.hpp"
#include "FileIdSupport.hpp"
#include "highlight.hpp"
#include "output.hpp"
#include "xref.hpp"

#include <boost/assert.hpp>

#include <climits>
#include <iostream>

using namespace synth;

// This is needed for macro arguments: clang_equalLocations is only true if two
// locations are truly equal. That is if either of spelling-, source- or
// file-location is different, it returns false. However we are only interested
// in the file location, hence this function.
// We must work with the offset here because it is in fact often different (and
// more "correct" for our purposes) from what line:column would suggest. See
// also comment inside processToken().
static bool equalFileLocations(CXSourceLocation loc1, CXSourceLocation loc2)
{
    CXFile f1;
    unsigned off1;
    clang_getFileLocation(loc1, &f1, nullptr, nullptr, &off1);

    CXFile f2;
    unsigned off2;
    clang_getFileLocation(loc2, &f2, nullptr, nullptr, &off2);

    return off1 == off2 && clang_File_isEqual(f1, f2);
}

static unsigned getLocOffset(CXSourceLocation loc)
{
    unsigned r;
    clang_getFileLocation(loc, nullptr, nullptr, nullptr, &r);
    return r;
}

namespace {

struct FileAnnotationState {
    CXFile file;
    HighlightedFile& hlFile;
    CgTokensHandle tokens;
    std::vector<CXCursor> annotations;
    std::vector<bool> annotationBad;

    // Maps from token begin offsets to indices (in tokens and tokCurs).
    std::unordered_map<unsigned, std::size_t> locationMap;

    void populateLocationMap(CXTranslationUnit tu)
    {
        for (std::size_t i = 0; i < tokens.size(); ++i) {
            CXToken tok = tokens.tokens()[i];
            unsigned off = getLocOffset(clang_getTokenLocation(tu, tok));
            BOOST_VERIFY(locationMap.insert({off, i}).second);
        }
    }
};

using TuAnnotationMap = std::unordered_map<CXFileUniqueID, FileAnnotationState>;

struct TuState {
    TuAnnotationMap annotationMap;
    CXTranslationUnit tu;
    MultiTuProcessor& multiTuProcessor;
    bool isC;
};

struct FileState {
    TuState& tuState;
    HighlightedFile& hlFile;
    bool lnkPending;
};

} // anonymous namespace

static FileAnnotationState* lookupFileAnnotations(TuAnnotationMap& m, CXFile f)
{
    CXFileUniqueID fuid;
    if (!f || clang_getFileUniqueID(f, &fuid) != 0)
        return nullptr;
    auto it = m.find(fuid);
    return it == m.end() ? nullptr : &it->second;
}

static void processToken(FileState& state, CXToken tok, CXCursor cur)
{
    auto& markups = state.hlFile.markups;
    markups.emplace_back();
    Markup* m = &markups.back();
    CXTranslationUnit tu = clang_Cursor_getTranslationUnit(cur);
    CXSourceRange rng = clang_getTokenExtent(tu, tok);
    unsigned lineno;
    clang_getFileLocation(
        clang_getRangeStart(rng), nullptr, &lineno, nullptr, &m->beginOffset);
    m->endOffset = getLocOffset(clang_getRangeEnd(rng));
    if (m->beginOffset == m->endOffset) {
        markups.pop_back();
        return;
    }

    CgStr hsp(clang_getTokenSpelling(tu, tok));
    boost::string_ref sp = hsp.gets();
    m->attrs = getTokenAttributes(tok, cur, sp);

    CXTokenKind tk = clang_getTokenKind(tok);
    if (tk == CXToken_Comment || tk == CXToken_Literal)
        return;

    CXCursorKind k = clang_getCursorKind(cur);
    if (state.lnkPending) {
        if (tk == CXToken_Punctuation && (sp == "(" || sp == "[")) {
            // This is the "("/"[" of an operator overload and we also want
            // to highlight the closing ")"/"]".
            return;
        }

        state.lnkPending = false;
        Markup lnk = {};
        lnk.beginOffset = getLocOffset(clang_getCursorLocation(cur));
        lnk.endOffset = m->endOffset;
        assert(lnk.beginOffset < lnk.endOffset);
        markups.push_back(std::move(lnk));
        m = &markups.back();
    } else if (!equalFileLocations(
        clang_getRangeStart(rng), clang_getCursorLocation(cur))
    ) {
        // Note that there is magic in the offset with which equalFileLocations
        // works (and clang_equalLocations too); it is sometimes different from
        // what line:column would suggest. Otherwise this condition would not
        // correctly work because e.g. for a function with built-in return type,
        // the cursor location of the FunctionDecl starts at the return type
        // declaration according to line:column, but at the function name (or
        // operator keyword or dtor tilde) according to the offset.
        return;
    } else if (k == CXCursor_InclusionDirective) {
        Markup incLnk = {};
        CXSourceRange incrng = clang_getCursorExtent(cur);
        incLnk.beginOffset = getLocOffset(clang_getRangeStart(incrng));
        incLnk.endOffset = getLocOffset(clang_getRangeEnd(incrng));
        linkCursor(incLnk, cur, state.tuState.multiTuProcessor);
        if (incLnk.isRef())
            state.hlFile.markups.push_back(std::move(incLnk));
        return;
    } else if (k == CXCursor_Destructor) {
        // This is the "~" of a dtor. Include the next part in the link.
        state.lnkPending = true;
        return;
    } else if (tk == CXToken_Keyword
        && (k == CXCursor_FunctionDecl || k == CXCursor_CXXMethod)
        && sp == "operator"
    ) {
        state.lnkPending = true;
        return;
    }

    // Non-KWs: Avoids emitting duplicate names at least for
    // template <typename T> using Foo = /* .. */;
    // where both the "using" and the "Foo" were independenty linked.
    // "{" was highlighted as definition for anonymous namespaces.
    if (tk != CXToken_Keyword
        && (tk != CXToken_Punctuation || (sp != "{" && sp != ";"))
    ) {
        SymbolDeclaration* decl = nullptr;
        auto const loadDecl = [&]() {
            if (decl)
                return;
            decl = &state.tuState.multiTuProcessor.createSymbol(
                state.hlFile,
                lineno,
                m->beginOffset);

            assert(decl->fileUniqueName.empty());
            std::size_t maxIdSz = state.tuState.multiTuProcessor.maxIdSz();
            if (maxIdSz > 0) {
                std::string name = fileUniqueName(cur, state.tuState.isC);
                if (name.size() < maxIdSz) {
                    decl->fileUniqueName = std::move(name);
                    m->fileUniqueName = &decl->fileUniqueName;
                }
            }
        };

        if (clang_isDeclaration(k)) {
            m->attrs |= TokenAttributes::flagDecl;
            loadDecl();
        }

        CXCursor defcur = clang_getCursorDefinition(cur);
        if (clang_equalCursors(defcur, cur)) { // This is a definition:
            m->attrs |= TokenAttributes::flagDef;
            loadDecl();
            CgStr usr(clang_getCursorUSR(cur));
            if (!usr.empty())
                state.tuState.multiTuProcessor.registerDef(usr.get(), decl);
        }
    }

    assert(m->beginOffset < m->endOffset);
    linkCursor(*m, cur, state.tuState.multiTuProcessor);
}

static CXChildVisitResult annotateVisit(
    CXCursor c, CXCursor, CXClientData ud)
{
    auto& state = *static_cast<TuState*>(ud);
    if (state.isC) {
        CXLanguageKind lang = clang_getCursorLanguage(c);
        state.isC = lang == CXLanguage_C || lang == CXLanguage_Invalid;
    }
    CXFile f;
    unsigned ln, off;
    clang_getFileLocation(
        clang_getCursorLocation(c), &f, &ln, nullptr, &off);
    FileAnnotationState* astate = lookupFileAnnotations(state.annotationMap, f);
    if (!astate)
        return ln == 0 ? CXChildVisit_Recurse : CXChildVisit_Continue;

    auto itIdx = astate->locationMap.find(off);
    if (itIdx == astate->locationMap.end()
        || !astate->annotationBad[itIdx->second]
    ) {
        return CXChildVisit_Recurse;
    }

    CXCursor& acur = astate->annotations[itIdx->second];
    acur = c;
    return CXChildVisit_Recurse;
}

static void annotate(TuState& state, CXCursor root)
{
    clang_visitChildren(root, &annotateVisit, &state);
}

static void processFile(
    CXFile file, CXSourceLocation*, unsigned, CXClientData ud)
{
    CXFileUniqueID fuid;
    if (clang_getFileUniqueID(file, &fuid) != 0)
        return;

    auto& state = *static_cast<TuState*>(ud);
    CXTranslationUnit tu = state.tu;

    CXSourceLocation beg = clang_getLocationForOffset(tu, file, 0);
    CXSourceLocation end = clang_getLocation(tu, file, UINT_MAX, UINT_MAX);

    HighlightedFile* hlFile = state.multiTuProcessor.prepareToProcess(file);
    if (!hlFile)
        return;

    CXToken* tokens;
    unsigned numTokens;
    clang_tokenize(tu, clang_getRange(beg, end), &tokens, &numTokens);
    CgTokensHandle hToks(tokens, numTokens, tu);

    if (numTokens == 0)
        return;

    std::vector<CXCursor> annotations(numTokens);
    clang_annotateTokens(tu, tokens, numTokens, annotations.data());
    std::vector<bool> annotationBad(numTokens);

    for (std::size_t i = 0; i < numTokens; ++i) {
        CXCursor& cur = annotations[i];
        CXSourceLocation tokLoc = clang_getTokenLocation(tu, tokens[i]);
        if (!equalFileLocations(clang_getCursorLocation(cur), tokLoc)) {
            CXCursor c2 = clang_getCursor(tu, tokLoc);
            CXSourceLocation loc2 = clang_getCursorLocation(c2);
            if (equalFileLocations(loc2, tokLoc))
                cur = c2;
            else
                annotationBad[i] = true;
        }
    }

    CgSourceRangesHandle skippedRanges(clang_getSkippedRanges(tu, file));
    for (unsigned i = 0; i < skippedRanges->count; ++i) {
        CXSourceRange rng = skippedRanges->ranges[i];
        std::pair<unsigned, unsigned> skippedRng;
        clang_getFileLocation(
            clang_getRangeStart(rng),
            nullptr,
            &skippedRng.first,
            nullptr,
            nullptr);
        clang_getFileLocation(
            clang_getRangeEnd(rng),
            nullptr,
            &skippedRng.second,
            nullptr,
            nullptr);
        hlFile->disabledLines.push_back(skippedRng);
    }
    std::sort(hlFile->disabledLines.begin(), hlFile->disabledLines.end());

    FileAnnotationState fstate {
        std::move(file),
        *hlFile,
        std::move(hToks),
        std::move(annotations),
        std::move(annotationBad),
        std::unordered_map<unsigned, std::size_t>()
    };
    auto kv = std::make_pair(std::move(fuid), std::move(fstate));
    auto insRes = state.annotationMap.insert(std::move(kv));
    assert(insRes.second);
    insRes.first->second.populateLocationMap(tu);

}

static void writeHlTokens(TuState& state)
{
    for (auto& fAnnotationsEntry : state.annotationMap) {
        FileAnnotationState& fAnnotations = fAnnotationsEntry.second;
        FileState fstate {state, fAnnotations.hlFile, /*lnkPending=*/false};
        for (std::size_t i = 0; i < fAnnotations.tokens.size(); ++i) {
            CXToken tok = fAnnotations.tokens.tokens()[i];
            CXCursor cur = fAnnotations.annotations[i];
            processToken(fstate, tok, cur);
        }
        fAnnotations.hlFile.markups.shrink_to_fit();
    }
}

int synth::processTu(
    CXIndex cidx,
    MultiTuProcessor& multiTuProcessor,
    char const* const* args,
    int nargs)
{
    CXTranslationUnit tu = nullptr;
    CXErrorCode err = clang_parseTranslationUnit2FullArgv(
        cidx,
        /*source_filename:*/ nullptr, // Included in commandline.
        args,
        nargs,
        /*unsaved_files:*/ nullptr,
        /*num_unsaved_files:*/ 0,
        CXTranslationUnit_DetailedPreprocessingRecord,
        &tu);
    CgTuHandle htu(tu);
    if (err != CXError_Success) {
        std::cerr << "Failed parsing translation unit (code "
                  << static_cast<int>(err)
                  << ")\n";
        std::cerr << "  args:";
        for (int i = 0; i < nargs; ++i)
            std::cerr << ' ' << args[i];
        std::cerr << '\n';
        return err + 10;
    }

    TuState state {TuAnnotationMap(), tu, multiTuProcessor, /*isC=*/ true};
    clang_getInclusions(tu, &processFile, &state);
    annotate(state, clang_getTranslationUnitCursor(tu));
    writeHlTokens(state);

    return EXIT_SUCCESS;
}


text:
//==================================================================================================
/*!

  Copyright 2015 NumScale SAS
  Copyright 2015 J.T. Lapreste

  Distributed under the Boost Software License, Version 1.0.
  (See accompanying file LICENSE.md or copy at http://boost.org/LICENSE_1_0.txt)
*/
//==================================================================================================
#include <boost/simd/function/scalar/asinpi.hpp>
#include <simd_test.hpp>
#include <boost/simd/constant/inf.hpp>
#include <boost/simd/constant/minf.hpp>
#include <boost/simd/constant/nan.hpp>
#include <boost/simd/constant/one.hpp>
#include <boost/simd/constant/mone.hpp>
#include <boost/simd/constant/zero.hpp>
#include <boost/simd/constant/mzero.hpp>

STF_CASE_TPL (" asinpi",  STF_IEEE_TYPES)
{
  namespace bs = boost::simd;
  namespace bd = boost::dispatch;
  using bs::asinpi;

  using r_t = decltype(asinpi(T()));

  // return type conformity test
  STF_TYPE_IS(r_t, T);

  // specific values tests
#ifndef BOOST_SIMD_NO_INVALIDS
  STF_ULP_EQUAL(asinpi(bs::Nan<T>()), bs::Nan<r_t>(), 0.5);
#endif
  STF_ULP_EQUAL(asinpi(bs::Half<T>()), T(1)/6, 0.5);
  STF_ULP_EQUAL(asinpi(bs::Mhalf<T>()), -T(1)/6, 0.5);
  STF_ULP_EQUAL(asinpi(bs::Mone<T>()), -0.5, 0.5);
  STF_ULP_EQUAL(asinpi(bs::One<T>()), 0.5, 0.5);
  STF_ULP_EQUAL(asinpi(bs::Zero<T>()), bs::Zero<r_t>(), 0.5);
}


text:
/*
 * The MIT License (MIT)
 *
 * Copyright (c) 2018 Sylko Olzscher
 *
 */

/**	@file service.hpp
 *
 *	Windows NT service template. Here we are in windows land.
 *	This code is shurely not portable.
 */
#pragma once

#include <boost/asio.hpp>

#include <boost/noncopyable.hpp>
#include <boost/bind.hpp>
#include <boost/function.hpp>

typedef std::function<void(DWORD)>	ctrl_handler_f;

/**
 *	The following template interface must be provided:
 *
 *	@li typedef service_impl_type and start_parameter_type
 *	@li Constructor with ( DWORD argc, char* argv[] )
 *	@li method: stop()
 *	@li method: run();
 *	@li method: reload();
 *	@li method: process_XML( start_parameter_type& );
 *	@li method: start_parameter_type::get_service_name()
 */
template < typename SRV_T >
class service
	: private boost::noncopyable
{

	/**
	 * Broker between windows API calls and implementation class.
	 */
	struct broker
	{
		static service*	this_service;

		static void WINAPI service_main(DWORD argc, char* argv[])
		{
			this_service->service_main(argc, argv);
			this_service->update_status();
		}

		static void WINAPI control_handler_(DWORD code)
		{
			this_service->control_handler(code);
		}
	};

public:
	/**
	 *	Constructor
	 *
	 *	The following events could be accepted
	 *	@li SERVICE_ACCEPT_STOP                    0x00000001
	 *	@li SERVICE_ACCEPT_PAUSE_CONTINUE          0x00000002
	 *	@li SERVICE_ACCEPT_SHUTDOWN                0x00000004
	 *	@li SERVICE_ACCEPT_PARAMCHANGE             0x00000008
	 *	@li SERVICE_ACCEPT_NETBINDCHANGE           0x00000010
	 *	@li SERVICE_ACCEPT_HARDWAREPROFILECHANGE   0x00000020
	 *	@li SERVICE_ACCEPT_POWEREVENT              0x00000040
	 *	@li SERVICE_ACCEPT_SESSIONCHANGE           0x00000080
	 *
	 */
	service(SRV_T&& srv, const std::string& service_name)
		: service_(std::move(srv))
		, handle_(NULL)
		, status_()
		, service_name_()
		, ctrl_handler_cb()
		, shutdown_(false)
	{
		BOOST_ASSERT_MSG(service_name.size() < sizeof(service_name_), "service name to long");
		BOOST_ASSERT_MSG(!service_name.empty(), "no service name");
		BOOST_ASSERT_MSG(service_name.at(0) != '"', "quoted service name");

		::memset(service_name_, 0, sizeof(service_name_));

		//	copy max. 256 characters
		::strcpy_s(service_name_, sizeof(service_name_), service_name.c_str());

		//status_.dwServiceType = SERVICE_WIN32_SHARE_PROCESS;
		status_.dwServiceType = SERVICE_WIN32_OWN_PROCESS;
		status_.dwCurrentState = SERVICE_STOPPED;
		status_.dwControlsAccepted = SERVICE_ACCEPT_STOP;
		status_.dwWin32ExitCode = NO_ERROR;
		status_.dwServiceSpecificExitCode = 0;
		status_.dwCheckPoint = 0;

		//	The estimated time required for a pending start, stop, pause, or continue operation, in milliseconds.
		status_.dwWaitHint = 12 * 1000;	//	12 sec

		BOOST_ASSERT_MSG(broker::this_service == nullptr, "wrong service broker");
		broker::this_service = this;

		std::stringstream ss;
		ss
			<< "service base class("
			<< service_name_
			<< ") successful constructed"
			<< std::endl
			;
		::OutputDebugString(ss.str().c_str());
	}

	/**
	 *	This parameters are not from the command line.
	 *
	 *	@param argc argument count
	 *	@param argv argument vector
	 */
	void service_main(DWORD argc, char* argv[])
	{
#ifdef _DEBUG
		//	now we have 30 sec. to startup
		//::DebugBreak();
#endif

		//
		//	create service shell on stack
		//

		if (register_control_handler(std::bind(& typename SRV_T::control_handler, &service_, std::placeholders::_1)))
		{
			::OutputDebugString("update status: SERVICE-START-PENDING\n");
			update_status(SERVICE_START_PENDING);

			::OutputDebugString("update status: SERVICE-RUNNING\n");
			update_status(SERVICE_RUNNING);

			//	startup server
			service_.run(false);

			while (!shutdown_)
			{
				//	Wait for all threads in the pool to exit.
				::Sleep(500);	//	1/2 second
			}
			::OutputDebugString("update status: SERVICE-STOP-PENDING\n");
			update_status(SERVICE_STOP_PENDING);
		}

		::Sleep(1500);	//	1+1/2 second more to calm down
		update_status(SERVICE_STOPPED);
		::OutputDebugString("update status: SERVICE-STOPPED\n");
	}

	/**
	*	run as service. Starts the dispatcher.
	*	This calls the service main() function and the
	*	service is running.
	*
	*	@return If the function succeeds, the return value is nonzero.
	*/
	DWORD run()
	{
		//	start dispatcher
		SERVICE_TABLE_ENTRY service_table[] =
		{
			{ service_name_, (LPSERVICE_MAIN_FUNCTION)&broker::service_main },
			{ NULL, NULL }
		};

		std::stringstream ss;
		ss
			<< "StartServiceCtrlDispatcher("
			<< service_name_
			<< ")"
			<< std::endl
			;
		::OutputDebugString(ss.str().c_str());

		if (broker::this_service != this)
		{
			::OutputDebugString("***Error: wrong service broker");
		}
		BOOST_ASSERT_MSG(broker::this_service == this, "wrong service broker");

		ss.str("");
		ss
			<< "main function address of "
			<< service_name_
			<< ": "
			<< (void*)&broker::service_main
			<< std::endl
			;
		::OutputDebugString(ss.str().c_str());

		if (0 == ::StartServiceCtrlDispatcher(service_table))
		{
			status_.dwWin32ExitCode = ::GetLastError();
		}
		return status_.dwWin32ExitCode;
	}

private:

	void control_handler(DWORD code)
	{
		std::stringstream ss;
		ss
			<< "service "
			<< service_name_
			<< " received control code: "
			<< code
			<< std::endl
			;
		::OutputDebugString(ss.str().c_str());

		//	Just switch on the control code sent to
		//	us and call the relavent virtual function
		switch (code)
		{
		case SERVICE_CONTROL_STOP:
			update_status(SERVICE_STOP_PENDING);
			try
			{
				shutdown_ = true;
				if (ctrl_handler_cb)
				{
					ctrl_handler_cb(SERVICE_CONTROL_STOP);
				}
			}
			catch (std::exception const& ex)
			{
				ss.str("");
				ss
					<< "service "
					<< service_name_
					<< " got an exception: "
					<< ex.what()
					<< std::endl
					;

				::OutputDebugString(ss.str().c_str());

			}
			catch (...)
			{
				ss.str("");
				ss
					<< "service "
					<< service_name_
					<< " got an exception "
					<< std::endl
					;

				::OutputDebugString(ss.str().c_str());
			}

#ifdef _SRV_DEBUG
			::OutputDebugString("NT Service: on_stop(SERVICE_STOPPED)\n");
#endif
			update_status(SERVICE_STOPPED);
			break;
		case SERVICE_CONTROL_PAUSE:
			//on_pause();
			break;
		case SERVICE_CONTROL_CONTINUE:
			//on_continue();
			break;
		case SERVICE_CONTROL_INTERROGATE:
			//on_interrogate();
			break;
		case SERVICE_CONTROL_SHUTDOWN:
			ctrl_handler_cb(SERVICE_CONTROL_SHUTDOWN);
			break;
		default:
			ctrl_handler_cb(code);
			break;
		}

		//	Any request from the SCM will be acked by this service
		//update_status();
	}

	/**
	 *	Update status in SCM database.
	 */
	bool update_status()
	{
		if (handle_ != NULL)
		{
			return 0 == ::SetServiceStatus(handle_, &status_);
		}
		return false;
	}

	bool update_status(DWORD status)
	{
		status_.dwCurrentState = status;
		return update_status();
	}

	bool register_control_handler(ctrl_handler_f cb)
	{

		//::DebugBreak();
		std::stringstream ss;
		ss
			<< "RegisterServiceCtrlHandler("
			<< service_name_
			<< ")"
			<< std::endl
			;
		::OutputDebugString(ss.str().c_str());

		if (broker::this_service != this)
		{
			::OutputDebugString("***Error: wrong service broker");
		}
		BOOST_ASSERT_MSG(broker::this_service == this, "wrong service broker");

		ss.str("");
		ss
			<< "control handler address of "
			<< service_name_
			<< ": "
			<< (void*)&broker::control_handler_
			<< std::endl
			;
		::OutputDebugString(ss.str().c_str());

		handle_ = ::RegisterServiceCtrlHandler(service_name_, (LPHANDLER_FUNCTION)&broker::control_handler_);
		if (handle_ == 0)
		{
			status_.dwWin32ExitCode = ::GetLastError();

			ss.str("");
			ss
				<< "RegisterServiceCtrlHandler("
				<< status_.dwWin32ExitCode
				<< ") failed"
				<< std::endl
				;
			::OutputDebugString(ss.str().c_str());

			return false;
		}

		ctrl_handler_cb.swap(cb);

		ss.str("");
		ss
			<< "RegisterServiceCtrlHandler("
			<< service_name_
			<< ") successful registered"
			<< std::endl
			;
		::OutputDebugString(ss.str().c_str());
		return true;
	}

private:
	SRV_T	service_;

	SERVICE_STATUS_HANDLE	handle_;
	SERVICE_STATUS			status_;
	char					service_name_[256];
	ctrl_handler_f			ctrl_handler_cb;
	bool					shutdown_;

	static	broker	broker_;
};

//
//	Initialize static members of this template
//
template < typename SRV_T >
service< SRV_T>*
service< SRV_T>::broker::this_service = nullptr;

max_stars_repo_name: solosTec/cyng

max_issues_repo_name: solosTec/cyng

max_forks_repo_name: solosTec/cyng