text:
/**********************************************************************
*  Copyright (c) 2008-2015, Alliance for Sustainable Energy.
*  All rights reserved.
*
*  This library is free software; you can redistribute it and/or
*  modify it under the terms of the GNU Lesser General Public
*  License as published by the Free Software Foundation; either
*  version 2.1 of the License, or (at your option) any later version.
*
*  This library is distributed in the hope that it will be useful,
*  but WITHOUT ANY WARRANTY; without even the implied warranty of
*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
*  Lesser General Public License for more details.
*
*  You should have received a copy of the GNU Lesser General Public
*  License along with this library; if not, write to the Free Software
*  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
**********************************************************************/

#ifndef UTILITIES_DATA_VECTOR_HPP
#define UTILITIES_DATA_VECTOR_HPP

#include "../UtilitiesAPI.hpp"

#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/vector.hpp>

namespace openstudio {

/// Workaround to get Vector typedef, http://www.gotw.ca/gotw/079.htm
struct UTILITIES_API VectorStruct
{
  typedef boost::numeric::ublas::vector<double> VectorType;
  typedef boost::numeric::ublas::scalar_vector<double> ScalarVectorType;
};

/// Vector
typedef VectorStruct::VectorType Vector;

/// ScalarVector
typedef VectorStruct::ScalarVectorType ScalarVector;

/// Helper function to construct Vector from std::vector<double>.
UTILITIES_API Vector createVector(const std::vector<double>& values);

/// Helper function to construct Vector from std::vector<long>
UTILITIES_API Vector createVector(const std::vector<long>& values);

UTILITIES_API std::vector<double> toStandardVector(const Vector& values);

//////////////////////////////////////////////////////////////////////////
// Begin SWIG'able, copy and paste into Vector.i
//////////////////////////////////////////////////////////////////////////

/** @name Operators */
//@{

UTILITIES_API bool operator==(const Vector& lhs, const Vector& rhs);
UTILITIES_API bool operator!=(const Vector& lhs, const Vector& rhs);

//@}
/** @name Interpolation */
//@{

// The following link has hints for interpolation
// http://o2scl.sourceforge.net/o2scl/html/index.html#intp_section

/** Enum to specify the interpolation method. */
enum InterpMethod{LinearInterp, NearestInterp, HoldLastInterp, HoldNextInterp};

/** Enum to specify the extrapolation method. */
enum ExtrapMethod{NoneExtrap, NearestExtrap};

/** Data structure for holding interpolation information. */
struct UTILITIES_API InterpInfo{
  bool extrapolated; // was point out of range
  unsigned ia, ib; // indices of two nearest points
  double wa, wb; // weights of two nearest points
};

/** Linear interpolation of the function y = f(x) at point xi. Assumes that x is strictly
 *  increasing. */
UTILITIES_API InterpInfo interpInfo(const Vector& x, double xi);

/** Linear interpolation of the function y = f(x) at point xi. Assumes that x is strictly
 *  increasing */
UTILITIES_API double interp(const Vector& x, const Vector& y, double xi,
                            InterpMethod interpMethod = LinearInterp,
                            ExtrapMethod extrapMethod = NoneExtrap);

/** Linear interpolation of the function y = f(x) at points xi. Assumes that x is strictly
 *  increasing. */
UTILITIES_API Vector interp(const Vector& x, const Vector& y, const Vector& xi,
                            InterpMethod interpMethod = LinearInterp,
                            ExtrapMethod extrapMethod = NoneExtrap);

//@}
/** @name Common Methods and Vector Operations */
//@{

/** Generates a Vector of N points randomly drawn between and including a and b. */
UTILITIES_API Vector randVector(double a, double b, unsigned N);

/** Generates a Vector of N points linearly spaced between and including a and b. */
UTILITIES_API Vector linspace(double a, double b, unsigned N);

/** Generates a Vector linearly spaced points starting at a and ending before or at b with
 *  interval delta. */
UTILITIES_API Vector deltaSpace(double a, double b, double delta);

/** Generates a Vector of N points logarithmically spaced between and including base^a and
 *  base^b. */
UTILITIES_API Vector logspace(double a, double b, unsigned N, double base = 10.0);

/** Take the natural logarithm of elements of a Vector. */
UTILITIES_API Vector log(const Vector& x);

/** Take the logarithm of elements of a Vector with certain base. */
UTILITIES_API Vector log(const Vector& x, double base);

/** Compute the cumulative sum of a Vector. */
UTILITIES_API Vector cumsum(const Vector& x, double runningSum = 0.0);

/** Returns the dot product between lhs and rhs. */
UTILITIES_API double dot(const Vector& lhs, const Vector& rhs);

/** Returns the sum of vector's values. */
UTILITIES_API double sum(const Vector& vector);

/** Returns the largest element of vector. */
UTILITIES_API double maximum(const Vector& vector);

/** Returns the smallest element of vector. */
UTILITIES_API double minimum(const Vector& vector);

/** Returns the mean of vector's values */
UTILITIES_API double mean(const Vector& vector);

/** Returns the sample variance of vector's values. */
UTILITIES_API double variance(const Vector& vector);

/** Returns the standard deviation of vector's values. */
UTILITIES_API double stdDev(const Vector& vector);

/** Returns std::function pointer to sum(const Vector&). */
UTILITIES_API std::function<double (const Vector&)> sumVectorFunctor();

/** Returns std::function pointer to maximum(const Vector&). */
UTILITIES_API std::function<double (const Vector&)> maximumVectorFunctor();

/** Returns std::function pointer to minimum(const Vector&). */
UTILITIES_API std::function<double (const Vector&)> minimumVectorFunctor();

/** Returns std::function pointer to mean(const Vector&). */
UTILITIES_API std::function<double (const Vector&)> meanVectorFunctor();

/** Returns std::function pointer to variance(const Vector&). */
UTILITIES_API std::function<double (const Vector&)> varianceVectorFunctor();

/** Returns std::function pointer to stdDev(const Vector&). */
UTILITIES_API std::function<double (const Vector&)> stdDevVectorFunctor();

/** Evaluates functor(vector). For use in SWIG bindings. */
UTILITIES_API double evaluateDoubleFromVectorFunctor(
    const std::function<double (const Vector&)>& functor,
    const Vector& vector);

//@}

} // openstudio

#endif //UTILITIES_DATA_VECTOR_HPP

ext: hpp


text:
#pragma once

#include <boost/math/special_functions/sign.hpp>

DFG_ROOT_NS_BEGIN { DFG_SUB_NS(math) {

// Returns sign of a value: 1 if value is positive, 0 if value is 0, -1 if value is < 0.
// Note: Behaviour for NaN's is unspecified.
// See also: std::signbit, std::copysign (C++11)
template <class T>
int sign(const T& val)
{
	return boost::math::sign(val);
	/*
	if (val > 0)
		return 1;
	else if (val == 0)
		return 0;
	else
		return -1;
		*/
}

template <class T>
auto signBit(const T& val) -> decltype(boost::math::signbit(val))
{
	return boost::math::signbit(val);
}

// Returns value whose absolute value is from val0 and sign from val1.
template <class T0, class T1>
auto signCopied(const T0& val0, const T1& val1) -> decltype(boost::math::copysign(val0, val1))
{
	return boost::math::copysign(val0, val1);
}

// Returns val with sign changed.
template <class T>
auto signChanged(const T& val) -> decltype(boost::math::changesign(val))
{
	return boost::math::changesign(val);
}

}} // module namespace


text:
#pragma once
#include <deip/protocol/transaction.hpp>

#include "deip_object_types.hpp"

#include <boost/multi_index/hashed_index.hpp>

namespace deip {
namespace chain {

using deip::protocol::signed_transaction;

/**
 * The purpose of this object is to enable the detection of duplicate transactions. When a transaction is included
 * in a block a transaction_object is added. At the end of block processing all transaction_objects that have
 * expired can be removed from the index.
 */
class transaction_object : public object<transaction_object_type, transaction_object>
{
    transaction_object() = delete;

public:
    template <typename Constructor, typename Allocator>
    transaction_object(Constructor&& c, allocator<Allocator> a)
        : packed_trx(a)
    {
        c(*this);
    }

    id_type id;

    bip::vector<char, allocator<char>> packed_trx;
    transaction_id_type trx_id;
    time_point_sec expiration;
};

struct by_expiration;
struct by_trx_id;
typedef multi_index_container<transaction_object,
                              indexed_by<ordered_unique<tag<by_id>,
                                                        member<transaction_object,
                                                               transaction_object_id_type,
                                                               &transaction_object::id>>,
                                         hashed_unique<tag<by_trx_id>,
                                                       BOOST_MULTI_INDEX_MEMBER(
                                                           transaction_object, transaction_id_type, trx_id),
                                                       std::hash<transaction_id_type>>,
                                         ordered_non_unique<tag<by_expiration>,
                                                            member<transaction_object,
                                                                   time_point_sec,
                                                                   &transaction_object::expiration>>>,
                              allocator<transaction_object>>
    transaction_index;
}
} // deip::chain

FC_REFLECT(deip::chain::transaction_object, (id)(packed_trx)(trx_id)(expiration))
CHAINBASE_SET_INDEX_TYPE(deip::chain::transaction_object, deip::chain::transaction_index)

max_stars_repo_name: DEIPworld/deip-chain

max_issues_repo_name: DEIPworld/deip-chain

max_forks_repo_name: DEIPworld/deip-chain


text:
#include <iostream>
#include <boost/tokenizer.hpp>
#include <string>
#include <boost/program_options.hpp>
#include <LRWrapper.h>
#include <LegionMatrix.h>
#include <boost/foreach.hpp>
#include <boost/algorithm/string.hpp>
#include <sstream>
#include "SlappyParams.h"
#include "Initializer.h"
#include "BetterMapper.h"
#include <MustEpochKernelLauncher.h>
#include <LegionHelper.h>
#include "InGroupSource.h"
#include "Sweep.h"
#include "SweepMapper.h"
//#include "gnuplot_i.h"

#include "CrossGroupSrc.h"
enum Task_IDs
{
	TOP_LEVEL_TASK_ID
};

using namespace Dragon;

int setup_diags(Context ctx, HighLevelRuntime* runtime,
                 LRWrapper& diags, LRWrapper& lines,
				 int& ndiag, int nx_l, int ny_l, int nz_l);

void setup_slgg(Context ctx, HighLevelRuntime* runtime,LRWrapper lw_slgg,int nmat,int nmom,int ng)
{
	double sigS[nmat][ng];
	for(int g=0;g<ng;g++)
		for(int m=0;m<nmat;m++)
		{
			if(g == 0)
			{
				if(m == 0)
					sigS[m][g] = 0.5;
				else
					sigS[m][g] = 1.2;
			}
			else
				sigS[m][g] = sigS[m][g-1]+0.005;
		}

RegionRequirement rr(lw_slgg.lr,READ_WRITE,EXCLUSIVE,lw_slgg.lr);
		rr.add_field(0);

		printf("Mapping slgg\n");
PhysicalRegion pr = runtime->map_region(ctx,rr);

LegionMatrix<double> slgg(lw_slgg,pr,0);

for(int g2 = 0;g2<ng;g2++)
	for(int g1=0;g1<ng;g1++)
		for(int m=0;m<nmat;m++)
		{
			if(m == 0)
			{
				if(g1 == g2)
					slgg(m,0,g1,g2) = 0.2*sigS[m][g1];
				else if(g2 < g1)
					if(g1 > 0)
						slgg(m,0,g1,g2) = 0.1*sigS[m][g1]/((double)g1);
					else
						slgg(m,0,g1,g2) = 0.3*sigS[m][g1];
				else
					if(g1 < ng-1)
						slgg(m,0,g1,g2) = 0.7*sigS[m][g1]/((double)(ng-g1-1));
					else
						slgg(m,0,g1,g2) = 0.9*sigS[m][g1];
			}
			else
			{
				if(g1 == g2)
					slgg(m,0,g1,g2) = 0.5*sigS[m][g1];
				else if(g2 < g1)
					if(g1 > 0)
						slgg(m,0,g1,g2) = 0.1*sigS[m][g1]/((double)g1);
					else
						slgg(m,0,g1,g2) = 0.6*sigS[m][g1];
				else
					if(g1 < ng-1)
						slgg(m,0,g1,g2) = 0.4*sigS[m][g1]/((double)(ng-g1-1));
					else
						slgg(m,0,g1,g2) = 0.9*sigS[m][g1];
			}
		}

if(nmom > 1)
{
	for(int m=0;m<nmat;m++)
	{
		double coef1,coef2;
		if(m==0)
		{coef1=0.1;coef2=0.5;}
		else
		{coef1=0.8;coef2=0.6;}
		for(int g2=0;g2<ng;g2++)
			for(int g1=0;g1<ng;g1++)
			{
				slgg(m,1,g1,g2) = coef1*slgg(m,0,g1,g2);
				for(int n=2;n<nmom;n++)
					slgg(m,n,g1,g2) = coef2*slgg(m,n-1,g1,g2);
			}

	}
}

runtime->unmap_region(ctx,pr);

}

void top_level_task(const Task *task,
                    const std::vector<PhysicalRegion> &regions,
                    Context ctx, HighLevelRuntime *runtime)
{

	runtime->disable_profiling();
	namespace po = boost::program_options;
	using namespace std;

	const InputArgs &command_args = HighLevelRuntime::get_input_args();
	SlappyParams params(command_args);

	int nx=params.nx;
	int ny=params.ny;
	int nz = params.nz;
	int nx_l(params.nx_l),ny_l(params.ny_l),nz_l(params.nz_l);
	int nCx(params.nCx), nCy(params.nCy), nCz(params.nCz);
	int ng(params.ng),nang(params.nang),nmom(params.nmom);
	int noct(params.noct);
	int nmat(params.nmat);
	int cmom(params.cmom);

	int& nDiag = params.nDiag;
	int& nPoints = params.nPoints;

	LRWrapper allArrays;
	allArrays.create(ctx,runtime,"allArrays",{1},0,LRWrapper(),1,LRWrapper(),2,LRWrapper(),3,LRWrapper());
	// Contains the by-chunk decomposition for all chunk dependent arrays
	LRWrapper chunkArrays;
	chunkArrays.create(ctx,runtime,"chunkArrays",{nCx,nCy,nCz},
	                   FIDs::Array::angFluxes,LRWrapper(),
	                   FIDs::Array::scalFluxes,LRWrapper(),
	                   FIDs::Array::savedFlux,LRWrapper(),
	                   FIDs::Array::sigmaT,LRWrapper(),
	                   FIDs::Array::sigmaS,LRWrapper(),
	                   FIDs::Array::mat,LRWrapper(),
	                   FIDs::Array::bci,LRWrapper(),
	                   FIDs::Array::bcj,LRWrapper(),
	                   FIDs::Array::bck,LRWrapper());

	// Setup all of the interior logical regions
	{
	RegionRequirement rr(chunkArrays.lr,WRITE_ONLY,EXCLUSIVE,chunkArrays.lr);
			rr.add_field(FIDs::Array::angFluxes);
			rr.add_field(FIDs::Array::scalFluxes);
			rr.add_field(FIDs::Array::savedFlux);
			rr.add_field(FIDs::Array::sigmaT);
			rr.add_field(FIDs::Array::sigmaS);
			rr.add_field(FIDs::Array::mat);
			rr.add_field(FIDs::Array::bci);
			rr.add_field(FIDs::Array::bcj);
			rr.add_field(FIDs::Array::bck);

	PhysicalRegion pr = runtime->map_region(ctx,rr);

	LegionMatrix<LRWrapper> angFluxes(chunkArrays,pr,FIDs::Array::angFluxes);
	LegionMatrix<LRWrapper> scalFluxes(chunkArrays,pr,FIDs::Array::scalFluxes);
	LegionMatrix<LRWrapper> savedFluxes(chunkArrays,pr,FIDs::Array::savedFlux);
	LegionMatrix<LRWrapper> sigmaT(chunkArrays,pr,FIDs::Array::sigmaT);
	LegionMatrix<LRWrapper> sigmaS(chunkArrays,pr,FIDs::Array::sigmaS);
	LegionMatrix<LRWrapper> mat(chunkArrays,pr,FIDs::Array::mat);
	LegionMatrix<LRWrapper> bci(chunkArrays,pr,FIDs::Array::bci);
	LegionMatrix<LRWrapper> bcj(chunkArrays,pr,FIDs::Array::bcj);
	LegionMatrix<LRWrapper> bck(chunkArrays,pr,FIDs::Array::bck);

	for(int k=0;k<nCz;k++)
		for(int j=0;j<nCy;j++)
			for(int i=0;i<nCx;i++)
			{
				{
				char* cell_id;
				asprintf(&cell_id," %i, %i, %i,",i,j,k);
				auto cn = [&](const char* name)
		{
					std::string res(name);
					return res+cell_id;
		};
				angFluxes(i,j,k) = (createLRWrapper(ctx,runtime,cn("angFluxes").c_str(),{nang,nx_l,ny_l,nz_l,ng,noct},
				                   FIDs::AngFlux::ptr_in,double(),
				                   FIDs::AngFlux::ptr_out,double()));

				scalFluxes(i,j,k) = createLRWrapper(ctx,runtime,cn("scalFluxes").c_str(),{cmom,nx_l,ny_l,nz_l,ng},
								   FIDs::ScalFlux::q2grp,double(),
								   FIDs::ScalFlux::qtot,double(),
				                   FIDs::ScalFlux::flux,double()
				                   );
				savedFluxes(i,j,k) = createLRWrapper(ctx,runtime,cn("savedFluxes").c_str(),{nx_l,ny_l,nz_l,ng},
				                   FIDs::SavedFlux::fluxpo,double(),
				                   FIDs::SavedFlux::fluxpi,double());

				sigmaT(i,j,k) = createLRWrapper(ctx,runtime,cn("sigmaT").c_str(),{nx_l,ny_l,nz_l,ng},
				                   FIDs::SigmaT::t_xs,double(),
				                   FIDs::SigmaT::a_xs,double(),
				                   FIDs::SigmaT::qi,double());

				sigmaS(i,j,k) = createLRWrapper(ctx,runtime,cn("sigmaS").c_str(),{nmom,nx_l,ny_l,nz_l,ng},
				                   FIDs::SigmaS::s_xs,double());

				mat(i,j,k) = createLRWrapper(ctx,runtime,cn("mat").c_str(),{nx_l,ny_l,nz_l},
				                   FIDs::Mat::mat,int());

				bci(i,j,k) = createLRWrapper(ctx,runtime,cn("bci").c_str(),{nang,ny_l,nz_l,ng},
				                   FIDs::BCi::ib_out,double(),1,double());
				bcj(i,j,k) = createLRWrapper(ctx,runtime,cn("bcj").c_str(),{nang,nx_l,nz_l,ng},
				                   FIDs::BCj::jb_out,double(),1,double());
				bck(i,j,k) = createLRWrapper(ctx,runtime,cn("bck").c_str(),{nang,nx_l,ny_l,ng},
				                   FIDs::BCk::kb_out,double(),1,double());
				}

			}

	runtime->unmap_region(ctx,pr);
	runtime->unmap_all_regions(ctx);

	}

	LRWrapper lw_slgg;
	lw_slgg.create(ctx,runtime,"slgg",{nmat,nmom,ng,ng},0,double());
	setup_slgg(ctx,runtime,lw_slgg,nmat,nmom,ng);

	// Setup all of the interior logical regions

	LRWrapper diags;
	LRWrapper lines;

	printf("Setting up Diags\n");
	nPoints = setup_diags(ctx,runtime,diags,lines,nDiag,nx_l,ny_l,nz_l);

	// Initialize stuff
	Initializer setup_stuff(params);

	auto init_op = genMustEpochKernel(setup_stuff,ctx,runtime,chunkArrays,lw_slgg);
//	runtime->unmap_all_regions(ctx);

	init_op.execute(ctx,runtime);
	printf("finished init op\n");

//	{
//
//	RegionRequirement rr(chunkArrays.lr,READ_ONLY,SIMULTANEOUS,chunkArrays.lr);
//			rr.add_field(FIDs::Array::sigmaT);
//			rr.add_field(FIDs::Array::sigmaS);
//
//
//
//
//	PhysicalRegion pr = runtime->map_region(ctx,rr);
//	LegionMatrix<LRWrapper> sigmaT(chunkArrays,pr,FIDs::Array::sigmaT);
//	LegionMatrix<LRWrapper> sigmaS(chunkArrays,pr,FIDs::Array::sigmaS);
//
//	float* x = (float*)malloc(nx*sizeof(float));
//	float* y = (float*)malloc(ny*sizeof(float));
//	float* qi_z = (float*)malloc(nx*ny*sizeof(float));
//
//
//	for(int jC=0;jC<nCy;jC++)
//		for(int iC=0;iC<nCx;iC++)
//		{
//			RegionRequirement rr_sigmaT(sigmaT(iC,jC,0).cast().lr,READ_ONLY,EXCLUSIVE,sigmaT(iC,jC,0).cast().lr);
//			rr_sigmaT.add_field(FIDs::SigmaT::t_xs);
//			rr_sigmaT.add_field(FIDs::SigmaT::qi);
//
//			RegionRequirement rr_sigmaS(sigmaS(iC,jC,0).cast().lr,READ_ONLY,EXCLUSIVE,sigmaS(iC,jC,0).cast().lr);
//			rr_sigmaS.add_field(FIDs::SigmaS::s_xs);
//
//			PhysicalRegion pr_sigmaT = runtime->map_region(ctx,rr_sigmaT);
//			PhysicalRegion pr_sigmaS = runtime->map_region(ctx,rr_sigmaS);
//
//			LegionMatrix<double> t_xs(sigmaT(iC,jC,0).cast(),pr_sigmaT,FIDs::SigmaT::t_xs);
//			LegionMatrix<double> qi(sigmaT(iC,jC,0).cast(),pr_sigmaT,FIDs::SigmaT::qi);
//			LegionMatrix<double> s_xs(sigmaS(iC,jC,0).cast(),pr_sigmaS,FIDs::SigmaS::s_xs);
//
//			for(int j=0;j<ny_l;j++)
//				for(int i=0;i<nx_l;i++)
//				{
//					int ix = i+nx_l*iC;
//					int iy = j+ny_l*jC;
//
//					x[ix] = ix;
//					y[iy] = iy;
////					printf("t_xs(%i, %i, 0) = %e\n",ix,iy,t_xs(i,j,0,0).cast());
//					qi_z[ix+nx*iy] = t_xs(i,j,0,0).cast();
//
//				}
//
//
//
//			runtime->unmap_region(ctx,pr_sigmaT);
//			runtime->unmap_region(ctx,pr_sigmaS);
//
//
//		}
//
//	gnuplot_ctrl* qi_plt = gnuplot_init();
//
//	gnuplot_plot_xyz(qi_plt,x,y,qi_z,nx,ny,"qi");
//	getchar();
//
//	runtime->unmap_region(ctx,pr);
//
//	runtime->unmap_all_regions(ctx);
//
//	}

	{
	RegionRequirement rr(allArrays.lr,READ_WRITE,EXCLUSIVE,allArrays.lr);
	rr.add_field(0);
	rr.add_field(1);
	rr.add_field(2);
	rr.add_field(3);

PhysicalRegion pr = runtime->map_region(ctx,rr);
	{

	LegionMatrix<LRWrapper> _allArrays(allArrays,pr,0);
	LegionMatrix<LRWrapper> _slgg(allArrays,pr,1);
	LegionMatrix<LRWrapper> _diag(allArrays,pr,2);
	LegionMatrix<LRWrapper> _lines(allArrays,pr,3);

	_allArrays(0) = chunkArrays;
	_slgg(0) = lw_slgg;
	_diag(0) = diags;
	_lines(0) = lines;

	}

	runtime->unmap_region(ctx,pr);

	}

	CrossGroupSrc outer(params);

	auto outer_op = genSingleKernel(outer, ctx, runtime, allArrays, chunkArrays, lw_slgg, diags,lines);

	printf("Executing outer op\n");
	runtime->execute_task(ctx,outer_op).get_void_result();
	printf("finished!\n");
}

static void update_mappers(Machine machine, HighLevelRuntime *rt,
						   const std::set<Processor> &local_procs)
{
	printf("Updating mappers\n");
  for (std::set<Processor>::const_iterator it = local_procs.begin();
		it != local_procs.end(); it++)
  {
//	rt->add_mapper(1,new SweepMapper(machine, rt, *it), *it);
    rt->replace_default_mapper(new BetterMapper(machine, rt, *it), *it);
  }
}

int main(int argc, char **argv)
{

  char hostname[512];
  gethostname(hostname,512);

  srand(23948723);
//
//	  dup2(fileno(err_file),STDERR_FILENO);
  HighLevelRuntime::set_top_level_task_id(TOP_LEVEL_TASK_ID);
  HighLevelRuntime::register_legion_task<top_level_task>(TOP_LEVEL_TASK_ID,
      Processor::LOC_PROC, true/*single*/, false/*index*/);

//	  fclose(err_file);

//	  std::stringstream buffer;
  setbuf(stdout,NULL);
  setvbuf(stderr,NULL,_IOLBF,1024);

  system("export GASNET_BACKTRACE=1");

//  TaskHelper::register_hybrid_variants<FieldKernelOp<Initializer>>();
//  TaskHelper::register_hybrid_variants<FieldKernelOp<CrossGroupSrc>>();
//  TaskHelper::register_hybrid_variants<FieldKernelOp<OuterSrcScat>>();
//  TaskHelper::register_hybrid_variants<FieldKernelOp<OuterConv>>();
//  TaskHelper::register_hybrid_variants<FieldKernelOp<SaveFluxes>>();
//  TaskHelper::register_hybrid_variants<FieldKernelOp<InGroupSource>>();
//  TaskHelper::register_hybrid_variants<FieldKernelOp<InnerSrcScat>>();

//  RegisterDim3Variants();

  TaskHelper::register_hybrid_variants<LegionHelper::Setter<bool>>();
  TaskHelper::register_hybrid_variants<LegionHelper::Setter<int>>();
  TaskHelper::register_hybrid_variants<LegionHelper::Setter<double>>();
  TaskHelper::register_hybrid_variants<LegionHelper::Setter<LRWrapper>>();
  TaskHelper::register_cpu_variants<LegionHelper::Getter<int>>();
  TaskHelper::register_cpu_variants<LegionHelper::Getter<bool>>();

  CrossGroupSrc::register_cpu();
  InnerSrcScat::register_cpu();

  Initializer::register_task();

  register_sweep_tasks();

//  TaskHelper::register_hybrid_variants<FieldKernelOp<Dim3Sweep<0>>>();
//  TaskHelper::register_hybrid_variants<FieldKernelOp<Dim3Sweep<1>>>();
//  TaskHelper::register_hybrid_variants<FieldKernelOp<Dim3Sweep<2>>>();
//  TaskHelper::register_hybrid_variants<FieldKernelOp<Dim3Sweep<3>>>();
//  TaskHelper::register_hybrid_variants<FieldKernelOp<Dim3Sweep<4>>>();
//  TaskHelper::register_hybrid_variants<FieldKernelOp<Dim3Sweep<5>>>();
//  TaskHelper::register_hybrid_variants<FieldKernelOp<Dim3Sweep<6>>>();
//  TaskHelper::register_hybrid_variants<FieldKernelOp<Dim3Sweep<7>>>();
//  TaskHelper::register_hybrid_variants<FieldKernelOp<Dim3Sweep<8>>>();
//  TaskHelper::register_hybrid_variants<FieldKernelOp<Dim3Sweep<9>>>();
//  TaskHelper::register_hybrid_variants<FieldKernelOp<Dim3Sweep<10>>>();
//  TaskHelper::register_hybrid_variants<FieldKernelOp<Dim3Sweep<11>>>();
//  TaskHelper::register_hybrid_variants<FieldKernelOp<Dim3Sweep<12>>>();
//  TaskHelper::register_hybrid_variants<FieldKernelOp<Dim3Sweep<13>>>();
//  TaskHelper::register_hybrid_variants<FieldKernelOp<Dim3Sweep<14>>>();
//  TaskHelper::register_hybrid_variants<FieldKernelOp<Dim3Sweep<15>>>();

	  HighLevelRuntime::set_registration_callback(update_mappers);

  return HighLevelRuntime::start(argc, argv);
}

ext: cc

max_stars_repo_name: tomdeakin/SNAP

max_issues_repo_name: tomdeakin/SNAP

max_forks_repo_name: tomdeakin/SNAP


text: #include <boost/metaparse/v1/error/digit_expected.hpp>

max_stars_repo_name: miathedev/BoostForArduino

max_issues_repo_name: miathedev/BoostForArduino

max_forks_repo_name: miathedev/BoostForArduino


text:
//
// Copyright 2017 Pixar
//
// Licensed under the Apache License, Version 2.0 (the "Apache License")
// with the following modification; you may not use this file except in
// compliance with the Apache License and the following modification to it:
// Section 6. Trademarks. is deleted and replaced with:
//
// 6. Trademarks. This License does not grant permission to use the trade
//    names, trademarks, service marks, or product names of the Licensor
//    and its affiliates, except as required to comply with Section 4(c) of
//    the License and to reproduce the content of the NOTICE file.
//
// You may obtain a copy of the Apache License at
//
//     http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the Apache License with the above modification is
// distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
// KIND, either express or implied. See the Apache License for the specific
// language governing permissions and limitations under the Apache License.
//
#include "undoInverse.h"

#include <boost/range/adaptor/reversed.hpp>

#include "undoBlock.h"
#include "undoRouter.h"

PXR_NAMESPACE_OPEN_SCOPE

void UsdQtUndoInverse::_Append(std::function<bool()> inverse) {
    _inversion.push_back(inverse);
}

void UsdQtUndoInverse::_Invert() {
    SdfChangeBlock changeBlock;
    for (const auto& inverse : boost::adaptors::reverse(_inversion)) {
        inverse();
    }
}

void UsdQtUndoInverse::Invert() {
    UsdQtUndoRouter& router = UsdQtUndoRouter::Get();
    if (router._depth != 0) {
        TF_CODING_ERROR(
            "Inversion during open edit block may result in corrupted undo "
            "stack.");
    }

    // open up an edit change block to capture the inverse of the inversion
    UsdQtUndoBlock editBlock;
    _Invert();
    _Clear();
    // adopt the edits and clear the listeners inversion tracker.
    // when the change block is popped, no notices will be sent
    // TODO: Do we want a more explicit version of this that
    // explicitly marks that we are inverting an undo/redo as
    // opposed to a new edit?
    _Adopt(router._inversion);
    router._inversion._Clear();
}

void UsdQtUndoInverse::_Clear() { _inversion.clear(); }

void UsdQtUndoInverse::_Adopt(const UsdQtUndoInverse& inversion) {
    for (const auto& inverse : inversion._inversion) {
        _inversion.push_back(inverse);
    }
}

UsdQtUndoInverse::UsdQtUndoInverse(UsdQtUndoRouter& router) {
    _Adopt(router._inversion);
    router._inversion._Clear();
}

PXR_NAMESPACE_CLOSE_SCOPE

ext: cpp

max_stars_repo_name: chunkified/usd-qt

max_issues_repo_name: chunkified/usd-qt

max_forks_repo_name: chunkified/usd-qt


text:
#include <benchmark/benchmark.h>

#include <thread>
#include <numeric>
#include <atomic>
#include <boost/format.hpp>
#include "ps/support/config.h"
#include "ps/base/frontend.h"
#include "ps/base/cards.h"
#include "ps/base/tree.h"
#include "ps/base/board_combination_iterator.h"
#include "ps/support/index_sequence.h"
#include "ps/eval/evaluator_5_card_map.h"
#include "ps/eval/evaluator_6_card_map.h"
#include "app/pretty_printer.h"
#include "ps/base/algorithm.h"
#include "ps/eval/instruction.h"
#include "ps/base/rank_hasher.h"
#include "ps/base/suit_hasher.h"
#include "ps/base/holdem_board_decl.h"
#include "ps/eval/class_cache.h"
#include "ps/eval/pass.h"
#include "ps/eval/pass_eval_hand_instr_vec.h"
#include "ps/base/rank_board_combination_iterator.h"

#include <boost/timer/timer.hpp>

#include <boost/log/trivial.hpp>

#include <Eigen/Dense>
#include <fstream>

#include "ps/support/command.h"
#include "ps/support/persistent.h"
#include "ps/eval/holdem_class_vector_cache.h"

using namespace ps;
using namespace ps::frontend;

/*
The point of this is to save the compution of mapping AA vs KK vs QQ to hard vectors,
then can re-apply the result to any matching pattern
*/
struct holdem_class_hard_from_proto
{
        using match_ty = std::tuple<holdem_class_type, rank_id, rank_id>;
        using rank_mapping_ty = std::array<int, rank_decl::max_id>;

        static std::string pattern_to_string(std::vector<match_ty> const& pattern)
        {
                std::stringstream ss;
                ss << "<";
                std::string sep = "";
                for(auto const& m : pattern)
                {
                        switch (std::get<0>(m))
                        {
                        case holdem_class_type::pocket_pair:
                                ss << "(pp," << rank_decl::get(std::get<1>(m)).to_string() << ")" << sep;
                                break;
                        case holdem_class_type::suited:
                                ss << "(suited," << rank_decl::get(std::get<1>(m)).to_string()
                                        << "," << rank_decl::get(std::get<2>(m)).to_string() << ")" << sep;
                                break;
                        case holdem_class_type::offsuit:
                                ss << "(offsuit," << rank_decl::get(std::get<1>(m)).to_string()
                                        << "," << rank_decl::get(std::get<2>(m)).to_string() << ")" << sep;
                        }
                        sep = ",";
                }
                ss << ">";
                return ss.str();
        }

        struct prototype
        {
                prototype(
                        std::vector< holdem_hand_vector > const& hv,
                        std::vector< matrix_t > const& transform)
                        : hv_(hv), transform_(transform)
                {}
                void emit(
                        holdem_class_vector const& cv,
                        rank_mapping_ty const& mapping,
                        std::vector< holdem_hand_vector >& out,
                        std::vector<matrix_t>& transform)const
                {
                        holdem_hand_vector working_hv;
                        for (size_t idx = 0; idx != hv_.size(); ++idx)
                        {
                                working_hv.clear();
                                for (holdem_id hid : hv_[idx])
                                {
                                        auto const& hand = holdem_hand_decl::get(hid);

                                        auto mapped_rank_first = mapping[hand.first().rank()];
                                        auto mapped_rank_second = mapping[hand.second().rank()];

                                        auto mapped_first = card_decl::make_id(
                                                hand.first().suit(),
                                                mapped_rank_first);

                                        auto mapped_second = card_decl::make_id(
                                                hand.second().suit(),
                                                mapped_rank_second);

                                        auto mapped_hid = holdem_hand_decl::make_id(mapped_first, mapped_second);

                                        working_hv.push_back(mapped_hid);

                                }

                                out.push_back(working_hv);
                                transform.push_back(transform_[idx]);
                        }
                }
                std::vector< holdem_hand_vector > hv_;
                std::vector< matrix_t > transform_;
        };
        void populate(holdem_class_vector const& cv, std::vector< holdem_hand_vector >& out, std::vector<matrix_t>& transform)
        {
                auto p = make_pattern(cv);

                auto const& pattern = std::get<0>(p);

                //std::cout << cv << " => " << pattern_to_string(pattern) << "\n";

                // mapping is cv rank to standard rank
                auto const& mapping = std::get<1>(p);
                auto const& inv_mapping = std::get<2>(p);

                auto iter = prototype_map_.find(pattern);

                if (iter == prototype_map_.end())
                {
                        //std::cout << "MAKING\n";

                        std::vector<frontend::range> players;
                        for (auto const& cid : cv)
                        {
                                auto const& decl = holdem_class_decl::get(cid);

                                auto mapped_cid = holdem_class_decl::make_id(
                                        decl.category(),
                                        mapping[decl.first()],
                                        mapping[decl.second()]);

                                players.push_back(frontend::range(holdem_class_decl::get(mapped_cid).to_string()));
                        }
                        auto instr_list = frontend_to_instruction_list("Group0", players);

                        computation_pass_manager mgr;
                        mgr.add_pass<pass_permutate_class>();
                        mgr.add_pass<pass_collect_class>();
                        mgr.add_pass<pass_class2cards>();
                        mgr.add_pass<pass_permutate>();
                        mgr.add_pass<pass_sort_type>();
                        mgr.add_pass<pass_collect>();

                        const size_t common_size = cv.size();
                        const int verboseicity = 0;
                        auto comp_ctx = std::make_shared< computation_context>(common_size, verboseicity);
                        mgr.execute_(comp_ctx.get(), &instr_list, nullptr);

                        std::vector<holdem_hand_vector> agg_hv;
                        std::vector<matrix_t> agg_trans;

                        for (auto const& untyped_instr : instr_list)
                        {
                                auto instr = reinterpret_cast<const card_eval_instruction*>(untyped_instr.get());

                                agg_hv.push_back(instr->get_vector());
                                agg_trans.push_back(instr->result_desc().back().transform());

                        }

                        prototype_map_[pattern] = std::make_shared<prototype>(agg_hv, agg_trans);
                }
                return prototype_map_.find(pattern)->second->emit(cv, inv_mapping, out, transform);
        }

        std::tuple<std::vector<match_ty>,rank_mapping_ty,rank_mapping_ty > make_pattern(holdem_class_vector const& cv)const
        {
                rank_mapping_ty rank_mapping;
                rank_mapping.fill(-1);
                rank_id rank_alloc(rank_decl::max_id);
                auto map_rank = [&rank_mapping,&rank_alloc](rank_id rid)mutable->rank_id
                {
                        if(rank_mapping[rid] == -1)
                        {
                                --rank_alloc;
                                const rank_id mapped_rank = rank_alloc;
                                rank_mapping[rid] = mapped_rank;
                                return mapped_rank;
                        }
                        else
                        {
                                return rank_mapping[rid];
                        }
                };

                std::vector<match_ty> pattern;
                for(holdem_class_id cid : cv)
                {
                        holdem_class_decl const& decl = holdem_class_decl::get(cid);

                        auto mapped_first_rank = map_rank(decl.first());
                        auto mapped_second_rank = map_rank(decl.second());
                        pattern.emplace_back(
                                decl.category(),
                                mapped_first_rank,
                                mapped_second_rank);
                }

                rank_mapping_ty inv_rank_mapping;
                inv_rank_mapping.fill(-1);
                for (rank_id iter = 0; iter != rank_decl::max_id; ++iter)
                {
                        if( rank_mapping[iter] != -1 )
                                inv_rank_mapping[rank_mapping[iter]] = iter;
                }

                //std::cout << "rank_mapping = " << std_vector_to_string(rank_mapping) << "\n";
                //std::cout << "inv_rank_mapping = " << std_vector_to_string(inv_rank_mapping) << "\n";
                return std::make_tuple(pattern,rank_mapping,inv_rank_mapping);

        }
private:
        std::map<
                std::vector<match_ty>,
                std::shared_ptr<prototype>
        > prototype_map_;
};

static void BM_ClassHands(benchmark::State& state)
{
        const auto cv = holdem_class_vector::parse("AA,KK,QQ");

        for (auto _ : state)
        {
                for( auto hv : cv.get_hand_vectors()){
                        benchmark::DoNotOptimize(hv);
                }
        }
}
BENCHMARK(BM_ClassHands);

static void BM_ClassHands_MakeProto_PP(benchmark::State& state)
{
        const auto cv = holdem_class_vector::parse("AA,KK,QQ");

        std::vector<holdem_hand_vector> out;
        std::vector<matrix_t> trans;

        for (auto _ : state)
        {
                holdem_class_hard_from_proto proto;
                out.clear();
                trans.clear();
                proto.populate(cv, out, trans);
        }
}
BENCHMARK(BM_ClassHands_MakeProto_PP);

static void BM_ClassHands_FromProto_PP(benchmark::State& state)
{
        const auto cv = holdem_class_vector::parse("AA,KK,QQ");

        std::vector<holdem_hand_vector> out;
        std::vector<matrix_t> trans;

        holdem_class_hard_from_proto proto;
        out.clear();
        trans.clear();
        proto.populate(cv, out, trans);

        for (auto _ : state)
        {
                out.clear();
                trans.clear();
                proto.populate(cv, out, trans);
        }
}
BENCHMARK(BM_ClassHands_FromProto_PP);

static void BM_ClassHands_MakeProto_ThreeOffsuit(benchmark::State& state)
{
        const auto cv = holdem_class_vector::parse("AKo,QJo,T9o");

        std::vector<holdem_hand_vector> out;
        std::vector<matrix_t> trans;

        for (auto _ : state)
        {
                holdem_class_hard_from_proto proto;
                out.clear();
                trans.clear();
                proto.populate(cv, out, trans);
        }
}
BENCHMARK(BM_ClassHands_MakeProto_ThreeOffsuit);

static void BM_ClassHands_FromProto_ThreeOffsuit(benchmark::State& state)
{
        const auto cv = holdem_class_vector::parse("AKo,QJo,T9o");

        std::vector<holdem_hand_vector> out;
        std::vector<matrix_t> trans;

        holdem_class_hard_from_proto proto;
        out.clear();
        trans.clear();
        proto.populate(cv, out, trans);

        for (auto _ : state)
        {
                out.clear();
                trans.clear();
                proto.populate(cv, out, trans);
        }
}
BENCHMARK(BM_ClassHands_FromProto_ThreeOffsuit);

static void BM_ClassHands_Proto_2PP(benchmark::State& state)
{
        const auto cv = holdem_class_vector::parse("AA,KK");

        holdem_class_hard_from_proto proto;

        std::vector<holdem_hand_vector> out;
        std::vector<matrix_t> trans;

        out.clear();
        trans.clear();
        proto.populate(holdem_class_vector::parse("AA,KK"), out, trans);
        std::cout << "-------------- first --------\n";
        for (size_t idx=0;idx!=out.size();++idx)
        {
                std::cout << out[idx] << " - " << matrix_to_string(trans[idx]) << "\n";
        }

        out.clear();
        trans.clear();
        proto.populate(holdem_class_vector::parse("AA,KK"), out, trans);
        std::cout << "-------------- second --------\n";
         for (size_t idx=0;idx!=out.size();++idx)
        {
                std::cout << out[idx] << " - " << matrix_to_string(trans[idx]) << "\n";
        }

        out.clear();
        trans.clear();
        proto.populate(holdem_class_vector::parse("KK,AA"), out, trans);
        std::cout << "-------------- second --------\n";
         for (size_t idx=0;idx!=out.size();++idx)
        {
                std::cout << out[idx] << " - " << matrix_to_string(trans[idx]) << "\n";
        }

        out.clear();
        trans.clear();
        proto.populate(holdem_class_vector::parse("99,22"), out, trans);
        std::cout << "-------------- second --------\n";
         for (size_t idx=0;idx!=out.size();++idx)
        {
                std::cout << out[idx] << " - " << matrix_to_string(trans[idx]) << "\n";
        }

        out.clear();
        trans.clear();
        proto.populate(holdem_class_vector::parse("AKo,QTo,98o"), out, trans);
        std::cout << "-------------- second --------\n";
         for (size_t idx=0;idx!=out.size();++idx)
        {
                std::cout << out[idx] << " - " << matrix_to_string(trans[idx]) << "\n";
        }

         out.clear();
        trans.clear();
        proto.populate(holdem_class_vector::parse("A2o,QTo,98o"), out, trans);
        std::cout << "-------------- second --------\n";
         for (size_t idx=0;idx!=out.size();++idx)
        {
                std::cout << out[idx] << " - " << matrix_to_string(trans[idx]) << "\n";
        }

        std::exit(0);
        return;

        out.clear();
        trans.clear();
        proto.populate(holdem_class_vector::parse("AA,JJ"), out, trans);
        std::cout << "-------------- third --------\n";
        for (auto const& cv : out)
        {
                std::cout << cv << "\n";
        }

}

#if 0
template <class ...Args>
static void HoldemHandVector_Permutate(benchmark::State& state, Args&&... args) {

        const auto hv = holdem_hand_vector::parse(args...);

        for (auto _ : state)
        {
                auto p =  permutate_for_the_better(hv) ;
                benchmark::DoNotOptimize(p);
        }
}
BENCHMARK_CAPTURE(HoldemHandVector_Permutate, JcTc7s7dAsKd, std::string("JcTc7s7dAsKd"))->Unit(benchmark::kMillisecond);
BENCHMARK_CAPTURE(HoldemHandVector_Permutate, JcTc7s7dAdQsAsKd, std::string("JcTc7s7dAdQsAsKd"))->Unit(benchmark::kMillisecond);
BENCHMARK_CAPTURE(HoldemHandVector_Permutate, JcTc7s7dAdQsAsKd8d9d, std::string("JcTc7s7dAdQsAsKd8d9d"))->Unit(benchmark::kMillisecond);
BENCHMARK_CAPTURE(HoldemHandVector_Permutate, JcTc7s7dAdQsAsKd8d9d3c2s, std::string("JcTc7s7dAdQsAsKd8d9d3c2s"))->Unit(benchmark::kMillisecond);
#endif

#if 0

template <class ...Args>
void ClassVector_GetHandVectors(benchmark::State& state, Args&&... args) {

        const auto cv = holdem_class_vector::parse(args...);

        for (auto _ : state)
        {
                for( auto hv : cv.get_hand_vectors()){
                        benchmark::DoNotOptimize(hv);
                }
        }
}
BENCHMARK_CAPTURE(ClassVector_GetHandVectors, AA_KK, std::string("AA,KK"))->Unit(benchmark::kMillisecond);
BENCHMARK_CAPTURE(ClassVector_GetHandVectors, AA_KK_QQ, std::string("AA,KK,QQ"))->Unit(benchmark::kMillisecond);
BENCHMARK_CAPTURE(ClassVector_GetHandVectors, QQ_KK_AA, std::string("QQ,KK,AA"))->Unit(benchmark::kMillisecond);
BENCHMARK_CAPTURE(ClassVector_GetHandVectors, AA_KK_QQ_JJ, std::string("AA,KK,QQ,JJ"))->Unit(benchmark::kMillisecond);
BENCHMARK_CAPTURE(ClassVector_GetHandVectors, AA_KK_QQ_JJ_TT, std::string("AA,KK,QQ,JJ,TT"))->Unit(benchmark::kMillisecond);
BENCHMARK_CAPTURE(ClassVector_GetHandVectors, AKo_QJo, std::string("AKo,QJo"))->Unit(benchmark::kMillisecond);
BENCHMARK_CAPTURE(ClassVector_GetHandVectors, AKo_QJo_T9o, std::string("AKo,QJo,T9o"))->Unit(benchmark::kMillisecond);
BENCHMARK_CAPTURE(ClassVector_GetHandVectors, AKo_QJo_T9o_87o, std::string("AKo,QJo,T9o,87o"))->Unit(benchmark::kMillisecond);
BENCHMARK_CAPTURE(ClassVector_GetHandVectors, AKo_QJo_T9o_87o_65o, std::string("AKo,QJo,T9o,87o,65o"))->Unit(benchmark::kMillisecond);

static void Iterator_BitmaskFiveCardBoard(benchmark::State& state) {

        //constexpr size_t end_board = 0b1111100000000000000000000000000000000000000000000000ull;
        constexpr size_t end_board = 0b11111000000000000000000000000000ull;
        constexpr size_t end = end_board + 1;
        for (auto _ : state)
        {
                size_t counter = 0;
                for (size_t iter = 0; iter < end; ++iter)
                {
                        if (detail::popcount(iter) != 5)
                        {
                                continue;
                        }
                        ++counter;

                        benchmark::DoNotOptimize(iter);
                }
                std::cout << "counter = " << counter << "\n";

        }
}
static void Iterator_FiveCardBoard(benchmark::State& state) {

        for (auto _ : state)
        {
                for(board_combination_iterator iter(5),end;iter!=end ;++iter){
                        benchmark::DoNotOptimize(*iter);
                }
        }
}
BENCHMARK(Iterator_FiveCardBoard)->Unit(benchmark::kMillisecond);

static void Iterator_ThreePlayerClass(benchmark::State& state) {

        for (auto _ : state)
        {
                for(holdem_class_iterator iter(3),end;iter!=end;++iter){
                        benchmark::DoNotOptimize(*iter);
                }
        }
}
BENCHMARK(Iterator_ThreePlayerClass)->Unit(benchmark::kMillisecond);

static void Vs_ClassRange_First1000(benchmark::State& state, size_t n) {

        for (auto _ : state)
        {
                size_t counter = 0;
                for(holdem_class_iterator iter(n),end;iter!=end;++iter,++counter){
                        auto const& cv = *iter;
                        for( auto hv : cv.get_hand_vectors()){
                                benchmark::DoNotOptimize(hv);
                        }
                        if( counter == 1000 )
                                break;
                }
        }
}
BENCHMARK_CAPTURE(Vs_ClassRange_First1000, Two, 2)->Unit(benchmark::kMillisecond);
BENCHMARK_CAPTURE(Vs_ClassRange_First1000, Three, 3)->Unit(benchmark::kMillisecond);
BENCHMARK_CAPTURE(Vs_ClassRange_First1000, Four, 4)->Unit(benchmark::kMillisecond);

#endif

BENCHMARK_MAIN();


text:
#include <stdlib.h>
#include <vector>
#include <string>
#include <Eigen/Geometry>
#include <ros/ros.h>
#include <moveit_msgs/GetPlanningScene.h>
#include <urdf_model/model.h>
#include <moveit/robot_model_loader/robot_model_loader.h>
#include <moveit/planning_scene/planning_scene.h>
#include "sdf_tools/sdf.hpp"
#include "sdf_tools/SDF.h"

#ifndef SDF_BUILDER_HPP
#define SDF_BUILDER_HPP

namespace sdf_tools
{
    static const uint8_t USE_CACHED = 0x00;
    static const uint8_t USE_ONLY_OCTOMAP = 0x01;
    static const uint8_t USE_ONLY_COLLISION_OBJECTS = 0x02;
    static const uint8_t USE_FULL_PLANNING_SCENE = 0x03;

    typedef struct
    {
        uint32_t location[3];
        uint32_t closest_point[3];
        double distance_square;
        int32_t update_direction;
    } bucket_cell;

    typedef VoxelGrid::VoxelGrid<bucket_cell> DistanceField;

    double ComputeDistanceSquared(int32_t x1, int32_t y1, int32_t z1, int32_t x2, int32_t y2, int32_t z2);

    class SDF_Builder
    {
    protected:

        bool initialized_;
        bool has_cached_sdf_;
        bool has_cached_collmap_;
        bool has_planning_scene_;
        Eigen::Isometry3d origin_transform_;
        std::string frame_;
        double x_size_;
        double y_size_;
        double z_size_;
        double resolution_;
        float OOB_value_;
        SignedDistanceField cached_sdf_;
        VoxelGrid::VoxelGrid<uint8_t> cached_collmap_;
        std::shared_ptr<planning_scene::PlanningScene> planning_scene_ptr_;
        ros::NodeHandle nh_;
        ros::ServiceClient planning_scene_client_;

        SignedDistanceField UpdateSDFFromPlanningScene();

        VoxelGrid::VoxelGrid<uint8_t> UpdateCollisionMapFromPlanningScene();

        bool BuildInternalPlanningScene();

        DistanceField BuildDistanceField(std::vector<Eigen::Vector3i>& points);

        std::vector<std::vector<std::vector<std::vector<int>>>> MakeNeighborhoods();

        static int GetDirectionNumber(int dx, int dy, int dz);

        std::string GenerateSDFComputeBotURDFString();

        std::string GenerateSDFComputeBotSRDFString();

    public:

        SDF_Builder(ros::NodeHandle& nh, Eigen::Isometry3d origin_transform, std::string frame, double x_size, double y_size, double z_size, double resolution, float OOB_value, std::string planning_scene_service);

        SDF_Builder(ros::NodeHandle& nh, std::string frame, double x_size, double y_size, double z_size, double resolution, float OOB_value, std::string planning_scene_service);

        SDF_Builder();

        void UpdatePlanningSceneFromMessage(moveit_msgs::PlanningScene& planning_scene);

        SignedDistanceField UpdateSDF(uint8_t update_mode);

        SignedDistanceField GetCachedSDF();

        VoxelGrid::VoxelGrid<uint8_t> UpdateCollisionMap(uint8_t update_mode);

        VoxelGrid::VoxelGrid<uint8_t> GetCachedCollisionMap();

    };

}

#endif // SDF_BUILDER_HPP


text:
// json parser
#include <boost/property_tree/json_parser.hpp>
#include <boost/property_tree/ptree.hpp>
#include <utility>
namespace pt = boost::property_tree;

#include "../GreensTensor/GreensTensorFactory.h"
#include "../MemoryKernel/MemoryKernelFactory.h"
#include "Polarizability.h"

// Constructor without an internal bath
Polarizability::Polarizability(double omega_a, double alpha_zero,
                               std::shared_ptr<GreensTensor> greens_tensor)
    : omega_a(omega_a), alpha_zero(alpha_zero),
      greens_tensor(std::move(greens_tensor)) {
  this->mu = nullptr;
}
// Constructor with internal bath mu
Polarizability::Polarizability(double omega_a, double alpha_zero,
                               std::shared_ptr<MemoryKernel> mu,
                               std::shared_ptr<GreensTensor> greens_tensor)
    : omega_a(omega_a), alpha_zero(alpha_zero), mu(std::move(mu)),
      greens_tensor(std::move(greens_tensor)) {}

Polarizability::Polarizability(const std::string &input_file) {
  // Create a root
  pt::ptree root;

  // Load the json file in this ptree
  pt::read_json(input_file, root);

  // read parameters
  this->omega_a = root.get<double>("Polarizability.omega_a");
  this->alpha_zero = root.get<double>("Polarizability.alpha_zero");

  // read greens tensor
  this->greens_tensor = GreensTensorFactory::create(input_file);

  // Check if the category Polarizability.MemoryKernel exists
  boost::optional<pt::ptree &> kernel_given =
      root.get_child_optional("Polarizability.MemoryKernel");

  if (kernel_given) {
    this->mu =
        MemoryKernelFactory::create(input_file, "Polarizability.MemoryKernel");
  } else {
    this->mu = nullptr;
  }
}

void Polarizability::calculate_tensor(double omega, cx_mat::fixed<3, 3> &alpha,
                                      Tensor_Options fancy_complex) const {
  // imaginary unit
  std::complex<double> I(0.0, 1.0);

  // calculate diagonal entries
  cx_mat::fixed<3, 3> diag;
  diag.zeros();
  if (mu != nullptr) {
    diag(0, 0) = diag(1, 1) = diag(2, 2) =
        (omega_a * omega_a - omega * omega) - I * omega * mu->calculate(omega);
  } else {
    diag(0, 0) = diag(1, 1) = diag(2, 2) = omega_a * omega_a - omega * omega;
  }
  // calculate integral over green's tensor with fancy R
  cx_mat::fixed<3, 3> greens_R;
  this->greens_tensor->integrate_k(omega, greens_R, RE, UNIT);

  // calculate integral over green's tensor with fancy I
  cx_mat::fixed<3, 3> greens_I;
  this->greens_tensor->integrate_k(omega, greens_I, IM, UNIT);

  // put everything together
  alpha =
      alpha_zero * omega_a * omega_a *
      inv(diag - alpha_zero * omega_a * omega_a * (greens_R + I * greens_I));

  if (fancy_complex == IM) {
    alpha = (alpha - trans(alpha)) /
            (2.0 * I); // trans is hermitean conjugation in armadillo
  } else if (fancy_complex == RE) {
    alpha = (alpha + trans(alpha)) /
            (2.0); // trans is hermitean conjugation in armadillo
  }
}

double Polarizability::integrand_omega(double omega,
                                       const uvec::fixed<2> &indices,
                                       Tensor_Options fancy_complex) const {
  cx_mat::fixed<3, 3> alpha;
  calculate_tensor(omega, alpha, fancy_complex);

  assert(imag(alpha(indices(0), indices(1))) == 0);
  double result = real(alpha(indices(0), indices(1)));

  return result;
}

double Polarizability::integrate_omega(const uvec::fixed<2> &indices,
                                       Tensor_Options fancy_complex,
                                       double omega_min, double omega_max,
                                       double relerr, double abserr) const {
  auto F = [=](double x) -> double {
    return this->integrand_omega(x, indices, fancy_complex);
  };
  return cquad(F, omega_min, omega_max, relerr, abserr);
}

void Polarizability::print_info(std::ostream &stream) const {
  stream << "# Polarizability\n#\n"
         << "# omega_a = " << omega_a << "\n"
         << "# alpha_zero = " << alpha_zero << "\n";
  if (mu != nullptr) {
    mu->print_info(stream);
  }
  greens_tensor->print_info(stream);
}

max_stars_repo_name: QuaCaTeam/quaca

max_issues_repo_name: myoelmy/quaca

max_forks_repo_name: myoelmy/quaca


text:
#include "VSCIMMapping.h"
#include "VSCIMMappingImplementations.h"

#include "VSCOBScene.h"
#include "VSCOBElement.h"
#include "VSCOBCollision.h"

#include <Ogre/Ogre.h>
#include <boost/assert.hpp>

using namespace VSC;
using namespace VSC::IM;

std::ostream& VSC::IM::operator << (std::ostream& stream, const Mapping& mapping)
{
    stream << "Mapping with type: " << StringForMappingType(mapping.getMappingType());
    stream << " offset: " << mapping.getOffset();
    stream << " scale factor: " << mapping.getScaleFactor();
    return stream;
}

Mapping::Mapping() : mOffset(0.0), mScaleFactor(1.0), mMappingType(MappingTypeNone)
{

}

Float Mapping::applyOffsetAndScaleFactor(Float rawValue)
{
    return mOffset + (rawValue * mScaleFactor);
}

void Mapping::setMappingType(MappingType mappingType)
{
    switch (mappingType)
    {
        case MappingTypeConstant:
            mImplementation = Implementation::SPtr(new ImplementationConstant);
            break;

        case MappingTypeCollisionVelocity:
            mImplementation = Implementation::SPtr(new ImplementationCollisionVelocity);
            break;

        case MappingTypeCollisionDistance:
            mImplementation = Implementation::SPtr(new ImplementationCollisionDistance);
            break;

        default:
            BOOST_ASSERT_MSG(false, "MappingType is not handled");
            mImplementation = Implementation::SPtr(new ImplementationConstant);
            mMappingType = MappingTypeConstant;
            return;
    }

    mMappingType = mappingType;

}

Float Mapping::mappedValue(Trigger trigger, TriggerPayload::SPtr payload)
{
    Float rawValue = 0.0;

    BOOST_ASSERT(mImplementation);
    if (mImplementation)
    {
        rawValue = mImplementation->mappedValue(trigger, payload);
    }

    return this->applyOffsetAndScaleFactor(rawValue);
}


text:
#include <vee/helper/strmagic.h>
#include <boost/tokenizer.hpp>

namespace vee {

namespace strmagic {

std::vector<::std::string> split(const ::std::string& string, const char* token)
{
    using char_separator = ::boost::char_separator<char>;
    using tokenizer = boost::tokenizer<char_separator>;

    // Parse raw data per lines (token: \n)
    ::std::vector<::std::string> ret;
    {
        char_separator sep(token);
        tokenizer tok(string, sep);
        for (auto& it : tok)
        {
            ret.push_back(it);
        }
    }

    return ret;
}

} // !namespace strmagic

} // !namespace vee


text:
#include "luminosityHadronic.h"

#include <fparameters/parameters.h>
#include <fmath/RungeKutta.h>
#include <fmath/interpolation.h>
#include <flosses/crossSectionInel.h>
#include <fmath/physics.h>
#include <algorithm>

#include <boost/math/special_functions/bessel.hpp>

/*double fpp(double Ep, double E)         //funcion a integrar   x=Eproton; E=Epion
{
	double L = log(Ep/1.6); //el 1.6 son TeV en erg
	double x = E/Ep;
	double Bg = 1.3+0.14*L+0.011*L*L;
	double beta = 1.0/(1.79+0.11*L+0.008*L*L);
	double kappa = 1.0/(0.801+0.049*L+0.014*L*L);

	double equis_b = pow(x,beta);
	double factor = 1-equis_b;
	double f;
	if (factor =! 0)	{

		double f1 = Bg*log(x)/x;
		double f2 = factor/(1.0+kappa*equis_b*factor);
		double f3 = 1.0/log(x) - 4.0*beta*equis_b/factor - 4.0*kappa*beta*equis_b*(1.0-2.0*equis_b)/(1.0+kappa*equis_b*factor);

		f      =  f1*pow(f2,4)*f3;
	}
	else	{
		f = 0.0;
	}
	return f;
}*/

double cHadron(double E)  //limite inferior
{
	double Kpi = 0.17;
	double thr = 0.0016; //1GeV

	return std::max(E+P2(neutralPionMass*cLight2)/(4*E),thr*Kpi); //== Ekin > Ethr
}

double dHadron()         //limite superior
{
	return 1.6e-12*pow(10.0,13.0);   //esto es un infinito
}

/*double fHadron(double Ep, double E, const Particle& p, const SpaceCoord& psc) //funcion a integrar   x=Ecreator; L=L(Ega)
{
	double Kpi = 0.17;
	double eval = p.mass*cLight2+Ep/Kpi;
	double Ekin = Ep/Kpi;
	double distCreator=0.0;
	if (Ep < p.emax() && Ep >= p.emin()) {
		distCreator = p.distribution.interpolate({ { 0, Ep } }, &psc);
	}

	double thr = 0.0016; //1GeV
	//double sigma = 30e-27*(0.95+0.06*log(Ekin/thr));

	double l = log10((protonMass*cLight2+Ep/Kpi)/1.6);
	double sigma = 1.e-27 * (34.3+1.88*l+0.25*l*l);
	double pionEmiss = sigma*distCreator*fpp(Ep,E)/Ep; //Kpi;  //sigma = crossSectionHadronicDelta(Ekin)

	double result = pionEmiss; ///sqrt(P2(Epi)-P2(neutralPionMass*cLight2));
	return result;
}*/

double heaviside(double x,double a,double b)
{
	return (a <= x && x <= b ? 1.0 : 0.0);
}

double auxf3(double dGeV, double sGeV, double gx, double GammaGeV, double isoGeV)
{
	double pGeV = protonMass*cLight2/1.6e-3;
	double piGeV = neutralPionMass*cLight2/1.6e-3;
	double gd = (sGeV + dGeV*dGeV - pGeV*pGeV)/(2.0*dGeV*sqrt(sGeV));
	double betad = sqrt(1.0-1.0/(gd*gd));
	double epi = (dGeV*dGeV+piGeV*piGeV-pGeV*pGeV)/(2.0*dGeV);
	double ppi = sqrt(epi*epi-piGeV*piGeV);
	double aux1 = 0.5/(betad*gd*ppi);   // REVISAR ESTE p_pi
	double aux2 = P2(dGeV-isoGeV)+GammaGeV*GammaGeV;
	double h = heaviside(gx*piGeV,gd*(epi-betad*ppi),gd*(epi+betad*ppi));

	return aux1*h/aux2;    // [GeV^-3]

}

double inclusiveSigma(double sGeV)
{
	double pGeV = protonMass*cLight2/1.6e-3;
	double piGeV = neutralPionMass*cLight2/1.6e-3;
	double eta=sqrt(P2(sGeV-P2(piGeV)-P2(2.0*pGeV))-4.0*P2(piGeV*2.0*pGeV))/(2.0*piGeV*sqrt(sGeV));
	double pthr = 0.78; //[GeV]
	double sigma=0.0;;

	// s = 2mp(Ek+2mp) = 2mp(E+mp)
	// E = s/2mp - mp
	// p = sqrt(E**2-m**2)
	double p = sqrt(P2(0.5*sGeV/pGeV-pGeV)-P2(pGeV));

	if (p >= pthr && p <= 0.96) {
		sigma = 0.032*eta*eta+0.04*pow(eta,6)+0.047*pow(eta,8);
	} else if (p > 0.96 && p <= 1.27) {
		sigma = 32.6*pow(p-0.8,3.21);
	} else if (p > 1.27 && p <= 8.0) {
		sigma = 5.4*pow(p-0.8,0.81);
	} else if (p > 8.0) {
		sigma = 32.0 * log(p) + 48.5 / sqrt(p) - 59.5;
	}
	return 1.0e-27 * sigma;     // [cm^2]
}
/*
double dsigma(double gx, double gr, double sGeV)
{
	double GammaGeV = 0.0575;
	double pGeV = protonMass*cLight2/1.6e-3;
	double piGeV = neutralPionMass * cLight2 / 1.6e-3;
	double isoGeV = 1.236;
	double atan1 = atan((sqrt(sGeV)-pGeV-isoGeV)/GammaGeV);
	double atan2 = atan((pGeV+piGeV-isoGeV)/GammaGeV);
	double aux1 = GammaGeV/(atan1-atan2);

	double Min = pGeV+piGeV;
	double Max = sqrt(sGeV)-pGeV;
	double integ= integSimpson(log(Min),log(Max),[&](double log_dGeV)
					{
						double dGeV = exp(log_dGeV);
						return auxf3(dGeV,sGeV,gx,GammaGeV,isoGeV)*dGeV;
					},30);
	return inclusiveSigma(sGeV)*aux1*integ*piGeV;  // [cm^2]
}*/

double dsigma(double gx, double gr, double sGeV)
{
	double GammaGeV = 0.0575;
	double pGeV = protonMass*cLight2/1.6e-3;
	double piGeV = neutralPionMass * cLight2 / 1.6e-3;
	double isoGeV = 1.236;
	double atan1 = atan((sqrt(sGeV)-pGeV-isoGeV)/GammaGeV);
	double atan2 = atan((pGeV+piGeV-isoGeV)/GammaGeV);
	double aux1 = GammaGeV/(atan1-atan2);

	double Min = pGeV+piGeV;
	double Max = sqrt(sGeV)-pGeV;

	double integ= RungeKuttaSimple(Min,Max,[&](double dGeV)
					{return auxf3(dGeV,sGeV,gx,GammaGeV,isoGeV);});

	return inclusiveSigma(sGeV)*aux1*integ*piGeV;  // [cm^2]
}

double auxf2(double gx, double gr, double epi, double normtemp, double sGeV)
{
	double piGeV = neutralPionMass*cLight2/1.6e-3;
	double g = epi/piGeV;
	double beta = sqrt(1.0-1.0/(g*g));
	double betax = sqrt(1.0-1.0/(gx*gx));
	double q = sqrt(2.0*(gr+1.0))/normtemp;

	double f1 = exp(-q*g*gx*(1.0-beta*betax))-exp(-q*g*gx*(1.0+beta*betax));

	return f1/(betax*gx) * dsigma(gx,gr,sGeV);
}

/*
double auxf(double gr, double epi, double normtemp)
{
	double piGeV = neutralPionMass*cLight2 / 1.6e-3;
	double pGeV = protonMass*cLight2/1.6e-3;
	double sGeV = 2.0*P2(pGeV)*(gr+1.0);
	double ji = (sGeV-4.0*P2(pGeV)+P2(piGeV))/(2.0*sqrt(sGeV));
	double Max = ji/piGeV;
	double integral = integSimpson(0.0,log(Max),[gr,epi,normtemp,sGeV](double log_gx)
						{
							double gx = exp(log_gx);
							return auxf2(gx,gr,epi,normtemp,sGeV)*gx;
						},30);
	double result = (gr*gr-1.0) / sqrt(2.0*(gr+1.0)) * integral;
}*/

double auxf(double gr, double epi, double normtemp)
{
	double piGeV = neutralPionMass*cLight2 / 1.6e-3;
	double pGeV = protonMass*cLight2/1.6e-3;
	double sGeV = 2.0*P2(pGeV)*(gr+1.0);
	double ji = (sGeV-4.0*P2(pGeV)+P2(piGeV))/(2.0*sqrt(sGeV));
	double Max = ji/piGeV;
	double integral = RungeKuttaSimple(1.0,Max,[gr,epi,normtemp,sGeV](double gx)
						{return auxf2(gx,gr,epi,normtemp,sGeV);});
	//double integral = qMidPointLog(1.0,Max,[gr,epi,normtemp,sGeV](double gx)
	//					{return auxf2(gx,gr,epi,normtemp,sGeV);},1.0e-2);
	double result = (gr*gr-1.0) / sqrt(2.0*(gr+1.0)) * integral;
}

/*
double fPion(double epi, double density, double temp)
{
	double normtemp = boltzmann* temp / (protonMass*cLight2);
	double k2 = boost::math::cyl_bessel_k(2, 1.0/normtemp);
	double piGeV = neutralPionMass*cLight2/1.6e-3;
	double constant = cLight*density*density/(4.0*piGeV*normtemp*k2*k2);
	return constant * integSimpson(0.0,log(1.0e3),[epi,normtemp](double log_gr)
						{
							double gr = exp(log_gr);
							return auxf(gr,epi,normtemp)*gr;
						},30);	//  [cm-3 s-1 GeV-1]
}*/

double fPion(double epi, double density, double temp)
{
	double normtemp = boltzmann* temp / (protonMass*cLight2);
	double k2 = boost::math::cyl_bessel_k(2, 1.0/normtemp);
	double piGeV = neutralPionMass*cLight2/1.6e-3;
	double constant = cLight*density*density/(4.0*piGeV*normtemp*k2*k2);

	return constant * RungeKuttaSimple(1.0,1.0e3,[&](double gr)
						{return auxf(gr,epi,normtemp);});         //  [cm-3 s-1 GeV-1]
	//return constant * qImpropLog(1.0,1.0e3,[epi,normtemp](double gr)
	//					{return auxf(gr,epi,normtemp);},1.0e-2);
}

double fHadron(double epi, double density, double temp)
{
	double piGeV = neutralPionMass*cLight2/1.6e-3;
	double qpi = fPion(epi,density,temp);
	return qpi / sqrt(P2(epi)-P2(piGeV));   // [cm-3 s-1 GeV-2]
}

/*
double luminosityHadronic(double E,
	const double density, double temp)
{
	double Kpi = 0.17;
	double thr = 0.0016; //1GeV

	//double Max  = dHadron();   //esto es un infinito
	double Min  = cHadron(E);
	//double Min = E+P2(0.5*neutralPionMass*cLight2)/E;
	Min = Min / 1.6e-3;
	double Max = 1.0e3;
	double integral = 2.0*integSimpson(log(Min),log(Max),[density,temp](double logepi)
						{
							double epi = exp(logepi);
							return fHadron(epi,density,temp)*epi;
						},30); // [cm-3 s-1 GeV-1]

	double jpp = integral * E*planck * 0.25/pi; // [erg s^-1 Hz^-1 cm^-3]
	jpp = jpp / (1.6e-3);
	return jpp;
}
*/

double luminosityHadronic(double E,
	const double density, double temp)
{
	double Kpi = 0.17;
	double thr = 0.0016; //1GeV

	//double Max  = dHadron();   //esto es un infinito
	double Min  = cHadron(E);
	//double Min = E+P2(0.5*neutralPionMass*cLight2)/E;
	Min = Min / 1.6e-3;
	double Max = 1.0e3;  // [en GeV]

	//double integral = 2.0*cLight*density*RungeKuttaSimple(Min, Max,
	//	[&E,&p,&psc](double x) {return fHadron(x, E, p, psc); });    //integra entre Emin y Emax

	double integral = 2.0*RungeKuttaSimple(Min,Max,[&](double epi)
							{return fHadron(epi,density,temp);}); // [cm-3 s-1 GeV-1]
	//double integral = 2.0*qMidPointLog(Min,Max,[density,temp](double epi)
	//						{return fHadron(epi,density,temp);},1.0e-2);
	//double integral = 2.0*qromoLog(Min,Max,[density,temp](double epi){return fHadron(epi,density,temp);});
	/*double integral = 2.0*cLight*density*RungeKutta(p.emin(), p.emax(),
		[E](double u){
			return cHadron(u,E);
		},
		[E](double u){
			return dHadron(u,E);
		},
		[&p,&psc](double u, double t){
			return fHadron(u,t,p, psc);
		}); */

	double jpp = integral * E*planck * 0.25/pi; // [erg s^-1 Hz^-1 cm^-3]
	jpp = jpp / (1.6e-3);
	return jpp;
}

/*
double fPP(double x, double E, Particle& creator )         //funcion a integrar   x=Eproton; E=Epion
{
	//DataInjection* data = (DataInjection*)voiddata;
	//double E = data->E;    //esta E corresponde a la energia del foton emitido; E=Ega
	//double mass   = data->mass;
	//Vector& Ncreator = data->Ncreator;
	//Vector& Ecreator = data->Ecreator;
	//const double mass   = particle.mass;
	//const Vector& Ncreator = creator.distribution.values;
	//const Vector& Ecreator = creator.eDim()->values;
	double L      = log(x/1.6); //el 1.6 son TeV en erg
	double ap     = 3.67+0.83*L+0.075*P2(L);
	double se     = crossSectionHadronic(x);
	double Bp     = ap+0.25;
	double r      = 2.6*pow(ap,-0.5);
	double alpha  = 0.98*pow(ap,-0.5);
	double equis  = E/x;
	double factor = 1-pow(equis,alpha);
	double f;
	if (factor =! 0)	{
		f      = 4*alpha*Bp*pow(equis,(alpha-1))*pow((factor/(1+r*pow(equis,alpha)*factor)),4)
      	         *(1/factor+r*(1-2*pow(equis,alpha))/(1+r*pow(equis,alpha)*factor))*
     	         pow((1-chargedPionMass*cLight2/(equis*x)),0.5);
	}
	else	{
		f = 0.0;
	}*/


text:
// Copyright (C) Ben Pope 2014.
// 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 <boost/graph/directed_graph.hpp>
#include <boost/graph/undirected_graph.hpp>

template<typename Graph>
void test_member_swap()
{
  Graph lhs, rhs;
  lhs.swap(rhs);
}

int main()
{
    test_member_swap<boost::adjacency_list<> >();
    test_member_swap<boost::directed_graph<> >();
    test_member_swap<boost::undirected_graph<> >();
}


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

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

#include <boost/fusion/support/config.hpp>
#include <boost/fusion/view/joint_view.hpp>

namespace boost { namespace fusion {

    namespace result_of
    {
        template<typename LhSequence, typename RhSequence>
        struct join
        {
            typedef joint_view<LhSequence, RhSequence> type;
        };
    }

    template<typename LhSequence, typename RhSequence>
    BOOST_FUSION_GPU_ENABLED
    inline typename result_of::join<LhSequence const, RhSequence const>::type
    join(LhSequence const& lhs, RhSequence const& rhs)
    {
        return typename result_of::join<LhSequence const, RhSequence const>::type(
            lhs, rhs);
    }
}}

#endif


text:
#include <pcl/point_types.h>
#include <pcl/features/normal_3d.h>
#include <pcl/search/kdtree.h>
#include <pcl/filters/extract_indices.h>
#include <pcl/octree/octree_pointcloud.h>

#include <pcl/features/vfh.h>
#include <pcl/io/pcd_io.h>

#include <Eigen/Dense>

#include <pcl/features/integral_image_normal.h>

#include "minipcl.h"

// set ksearch and radius to < 0 to disable
void mpcl_compute_normals(const pcl::PointCloud<pcl::PointXYZ>& cloud,
                          int ksearch,
                          double searchRadius,
                          pcl::PointCloud<pcl::Normal> &out)
{
    pcl::search::KdTree<pcl::PointXYZ>::Ptr tree (new pcl::search::KdTree<pcl::PointXYZ> ());
    pcl::NormalEstimation<pcl::PointXYZ, pcl::Normal> ne;

    ne.setSearchMethod (tree);
    ne.setInputCloud (cloud.makeShared());
    if (ksearch >= 0)
        ne.setKSearch (ksearch);
    if (searchRadius >= 0.0)
        ne.setRadiusSearch (searchRadius);
    ne.compute (out);
}

void mpcl_compute_normals_PointXYZI(const pcl::PointCloud<pcl::PointXYZI>& cloud,
                          int ksearch,
                          double searchRadius,
                          pcl::PointCloud<pcl::Normal> &out)
{
    pcl::search::KdTree<pcl::PointXYZI>::Ptr tree (new pcl::search::KdTree<pcl::PointXYZI> ());
    pcl::NormalEstimation<pcl::PointXYZI, pcl::Normal> ne;

    ne.setSearchMethod (tree);
    ne.setInputCloud (cloud.makeShared());
    if (ksearch >= 0)
        ne.setKSearch (ksearch);
    if (searchRadius >= 0.0)
        ne.setRadiusSearch (searchRadius);
    ne.compute (out);
}

void mpcl_compute_normals_PointXYZRGB(const pcl::PointCloud<pcl::PointXYZRGB>& cloud,
                          int ksearch,
                          double searchRadius,
                          pcl::PointCloud<pcl::Normal> &out)
{
    pcl::search::KdTree<pcl::PointXYZRGB>::Ptr tree (new pcl::search::KdTree<pcl::PointXYZRGB> ());
    pcl::NormalEstimation<pcl::PointXYZRGB, pcl::Normal> ne;

    ne.setSearchMethod (tree);
    ne.setInputCloud (cloud.makeShared());
    if (ksearch >= 0)
        ne.setKSearch (ksearch);
    if (searchRadius >= 0.0)
        ne.setRadiusSearch (searchRadius);
    ne.compute (out);
}

void mpcl_compute_normals_PointXYZRGBA(const pcl::PointCloud<pcl::PointXYZRGBA>& cloud,
                          int ksearch,
                          double searchRadius,
                          pcl::PointCloud<pcl::Normal> &out)
{
    pcl::search::KdTree<pcl::PointXYZRGBA>::Ptr tree (new pcl::search::KdTree<pcl::PointXYZRGBA> ());
    pcl::NormalEstimation<pcl::PointXYZRGBA, pcl::Normal> ne;

    ne.setSearchMethod (tree);
    ne.setInputCloud (cloud.makeShared());
    if (ksearch >= 0)
        ne.setKSearch (ksearch);
    if (searchRadius >= 0.0)
        ne.setRadiusSearch (searchRadius);
    ne.compute (out);
}

// set ksearch and radius to < 0 to disable
void mpcl_sacnormal_set_axis(pcl::SACSegmentationFromNormals<pcl::PointXYZ, pcl::Normal> &sac,
                             double ax, double ay, double az)
{
    Eigen::Vector3f vect(ax,ay,az);
    sac.setAxis(vect);
}

void mpcl_sacnormal_set_axis_PointXYZI(pcl::SACSegmentationFromNormals<pcl::PointXYZI, pcl::Normal> &sac,
                             double ax, double ay, double az)
{
    Eigen::Vector3f vect(ax,ay,az);
    sac.setAxis(vect);
}

void mpcl_sacnormal_set_axis_PointXYZRGB(pcl::SACSegmentationFromNormals<pcl::PointXYZRGB, pcl::Normal> &sac,
                             double ax, double ay, double az)
{
    Eigen::Vector3f vect(ax,ay,az);
    sac.setAxis(vect);
}

void mpcl_sacnormal_set_axis_PointXYZRGBA(pcl::SACSegmentationFromNormals<pcl::PointXYZRGBA, pcl::Normal> &sac,
                             double ax, double ay, double az)
{
    Eigen::Vector3f vect(ax,ay,az);
    sac.setAxis(vect);
}

//
void mpcl_extract(pcl::PointCloud<pcl::PointXYZ>::Ptr &incloud,
                  pcl::PointCloud<pcl::PointXYZ> *outcloud,
                  pcl::PointIndices *indices,
                  bool negative)
{
    pcl::PointIndices::Ptr indicesptr (indices);
    pcl::ExtractIndices<pcl::PointXYZ> ext;
    ext.setInputCloud(incloud);
    ext.setIndices(indicesptr);
    ext.setNegative(negative);
    ext.filter(*outcloud);
}

void mpcl_extract_PointXYZI(pcl::PointCloud<pcl::PointXYZI>::Ptr &incloud,
                   pcl::PointCloud<pcl::PointXYZI> *outcloud,
                   pcl::PointIndices *indices,
                   bool negative)
{
    pcl::PointIndices::Ptr indicesptr (indices);
    pcl::ExtractIndices<pcl::PointXYZI> ext;
    ext.setInputCloud(incloud);
    ext.setIndices(indicesptr);
    ext.setNegative(negative);
    ext.filter(*outcloud);
}

void mpcl_extract_PointXYZRGB(pcl::PointCloud<pcl::PointXYZRGB>::Ptr &incloud,
                   pcl::PointCloud<pcl::PointXYZRGB> *outcloud,
                   pcl::PointIndices *indices,
                   bool negative)
{
    pcl::PointIndices::Ptr indicesptr (indices);
    pcl::ExtractIndices<pcl::PointXYZRGB> ext;
    ext.setInputCloud(incloud);
    ext.setIndices(indicesptr);
    ext.setNegative(negative);
    ext.filter(*outcloud);
}

void mpcl_extract_PointXYZRGBA(pcl::PointCloud<pcl::PointXYZRGBA>::Ptr &incloud,
                   pcl::PointCloud<pcl::PointXYZRGBA> *outcloud,
                   pcl::PointIndices *indices,
                   bool negative)
{
    pcl::PointIndices::Ptr indicesptr (indices);
    pcl::ExtractIndices<pcl::PointXYZRGBA> ext;
    ext.setInputCloud(incloud);
    ext.setIndices(indicesptr);
    ext.setNegative(negative);
    ext.filter(*outcloud);
}

// EuclideanClusterExtraction
// Octree
// void mpcl_deleteVoxelAtPoint(pcl::octree::OctreePointCloud<pcl::PointXYZ>& inOctree, pcl::PointXYZ incloud)
// {
//     inOctree.deleteVoxelAtPoint(incloud);
// }
//
// void mpcl_deleteVoxelAtPoint(pcl::octree::OctreePointCloud<pcl::PointXYZI>& inOctree, pcl::PointXYZI incloud)
// {
//     inOctree.deleteVoxelAtPoint(incloud);
// }
//
// void mpcl_deleteVoxelAtPoint(pcl::octree::OctreePointCloud<pcl::PointXYZRGB>& inOctree, pcl::PointXYZRGB incloud)
// {
//     inOctree.deleteVoxelAtPoint(incloud);
// }
//
// void mpcl_deleteVoxelAtPoint(pcl::octree::OctreePointCloud<pcl::PointXYZRGBA>& inOctree, pcl::PointXYZRGBA incloud)
// {
//     inOctree.deleteVoxelAtPoint(incloud);
// }
//
//
// int mpcl_getOccupiedVoxelCenters(pcl::octree::OctreePointCloud<pcl::PointXYZ>& inOctree, pcl::PointCloud<pcl::PointXYZRGBA>::Ptr &incloud, vector<pcl::PointXYZ, Eigen::aligned_allocator<pcl::PointXYZ> > alignPoint)
// {
//     return inOctree.getOccupiedVoxelCenters(alignPoint);
// }
//
// int mpcl_getOccupiedVoxelCenters(pcl::octree::OctreePointCloud<pcl::PointXYZI>& inOctree, pcl::PointCloud<pcl::PointXYZRGBA>::Ptr &incloud, vector<pcl::PointXYZI, Eigen::aligned_allocator<pcl::PointXYZI> > alignPoint)
// {
//     return inOctree.getOccupiedVoxelCenters(alignPoint);
// }
//
// int mpcl_getOccupiedVoxelCenters(pcl::octree::OctreePointCloud<pcl::PointXYZRGB>& inOctree, pcl::PointCloud<pcl::PointXYZRGB>::Ptr &incloud, vector<pcl::PointXYZRGB, Eigen::aligned_allocator<pcl::PointXYZRGB> > alignPoint)
// {
//     return inOctree.getOccupiedVoxelCenters(alignPoint);
// }
//
// int mpcl_getOccupiedVoxelCenters(pcl::octree::OctreePointCloud<pcl::PointXYZRGBA>& inOctree, pcl::PointCloud<pcl::PointXYZRGBA>::Ptr &incloud, vector<pcl::PointXYZRGBA, Eigen::aligned_allocator<pcl::PointXYZRGBA> > alignPoint)
// {
//     return inOctree.getOccupiedVoxelCenters(alignPoint);
// }
//

void mpcl_extract_VFH(pcl::PointCloud<pcl::PointXYZ>::Ptr cloud)
{
    // http://virtuemarket-lab.blogspot.jp/2015/03/viewpoint-feature-histogram.html
    // pcl::PointCloud<pcl::VFHSignature308>::Ptr Extract_VFH(pcl::PointCloud<pcl::PointXYZ>::Ptr cloud)
    pcl::PointCloud<pcl::Normal>::Ptr cloud_normals (new pcl::PointCloud<pcl::Normal> ());
    pcl::VFHEstimation<pcl::PointXYZ, pcl::Normal, pcl::VFHSignature308> vfh;
    pcl::PointCloud<pcl::VFHSignature308>::Ptr vfhs (new pcl::PointCloud<pcl::VFHSignature308> ());

    // cloud_normals = surface_normals(cloud);
    vfh.setInputCloud (cloud);
    vfh.setInputNormals (cloud_normals);

    pcl::search::KdTree<pcl::PointXYZ>::Ptr tree (new pcl::search::KdTree<pcl::PointXYZ> ());
    vfh.setSearchMethod (tree);

    vfh.compute (*vfhs);
    // return vfhs;
}

/*
// pcl1.6
#include <pcl/keypoints/harris_keypoint3D.h>
// use 1.7�`
#include <pcl/keypoints/harris_3d.h>

// HarrisKeypoint3D
// NG
// outcloud set pcl::PointXYZI
void mpcl_extract_HarrisKeypoint3D(pcl::PointCloud<pcl::PointXYZ>::Ptr &incloud,
                                   pcl::PointCloud<pcl::PointXYZ> *outcloud)
{
    // pcl::HarrisKeypoint3D<pcl::PointXYZ, pcl::PointXYZ> detector;
    pcl::HarrisKeypoint3D<pcl::PointXYZ, pcl::PointXYZI> detector;

    detector.setInputCloud(incloud);

    detector.setNonMaxSupression (true);
    detector.setRadius (0.01);
    // detector.setRadiusSearch (100);
    // detector.setIndices(indicesptr);

    // NG
    // detector.compute(*outcloud);

    // OK
    pcl::PointCloud<pcl::PointXYZI>::Ptr keypoints(new pcl::PointCloud<pcl::PointXYZI>());
    detector.compute(*keypoints);
}

// HarrisKeypoint3D
void mpcl_extract_HarrisKeypoint3D(pcl::PointCloud<pcl::PointXYZ>::Ptr &incloud,
                                   pcl::PointCloud<pcl::PointXYZI> *outcloud)
{
    pcl::HarrisKeypoint3D<pcl::PointXYZ, pcl::PointXYZI> detector;

    detector.setInputCloud(incloud);

    detector.setNonMaxSupression (true);
    detector.setRadius (0.01);
    //detector.setRadiusSearch (100);

    // detector.setIndices(indicesptr);
    // detector.compute(*outcloud);
    detector.compute(*outcloud);
}
*/

// features
// integral_image_normal.h
void mpcl_features_NormalEstimationMethod_AVERAGE_3D_GRADIENT(pcl::IntegralImageNormalEstimation<pcl::PointXYZ, pcl::Normal> ne)
{
    ne.setNormalEstimationMethod (ne.AVERAGE_3D_GRADIENT);
}

void mpcl_features_NormalEstimationMethod_COVARIANCE_MATRIX(pcl::IntegralImageNormalEstimation<pcl::PointXYZ, pcl::Normal> ne)
{
    ne.setNormalEstimationMethod (ne.COVARIANCE_MATRIX);
}

void mpcl_features_NormalEstimationMethod_AVERAGE_DEPTH_CHANGE(pcl::IntegralImageNormalEstimation<pcl::PointXYZ, pcl::Normal> ne)
{
    ne.setNormalEstimationMethod (ne.AVERAGE_DEPTH_CHANGE);
}

void mpcl_features_NormalEstimationMethod_SIMPLE_3D_GRADIENT(pcl::IntegralImageNormalEstimation<pcl::PointXYZ, pcl::Normal> ne)
{
    ne.setNormalEstimationMethod (ne.SIMPLE_3D_GRADIENT);
}

void mpcl_features_NormalEstimationMethod_compute(pcl::IntegralImageNormalEstimation<pcl::PointXYZ, pcl::Normal> ne, pcl::PointCloud<pcl::Normal> &out)
{
    // NG : out Variant Function end error
    printf("compute start.\n");
    // �Q�ƃJ�E���g��NG?(�֐����������ɃG���[)
    ne.compute (out);
    // pcl 1.7.2 error
    // printf("out = %p.\n", out);
    printf("compute end.\n");
}


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

//// Solver
// Build AIC matrices and solve iteratively the "body" and the "field"
// Body module: compute impermeability condition and solve linear system of equations (corresponding to the
// Linear Potential Equation)
// Field module: compute through AIC and FD the variables in the field, as well as the field sources (source term in
// the Full Potential Equation)
//
// I/O:
// - numC: list of numerical parameters (structure)
// - symY: defines symmetry about Y axis
// - sRef: reference surface of the full wing
// - alpha: freestream angle of attack
// - vInf: freestream velocity vector
// - Minf: freestream Mach number
// - bPan: (network of) body panels (structure)
// - wPan: (network of) wake panels (structure)
// - fPan: field panels (structure)
// - sp: sub-panels (structure)
// - cL: lift coefficient
// - cD: drag coefficient
//
// Output:
// Outputs:
// - 0, if function succeeded
// - 1, if function failed

#include <iostream>
#include <Eigen/Dense>
#include "solver.h"
#include "id_subpanel.h"
#include "build_AIC.h"
#include "solve_body.h"
#include "solve_field.h"
#include "compute_sVars.h"

#define ANSI_COLOR_RED     "\x1b[1;31m"
#define ANSI_COLOR_GREEN   "\x1b[1;32m"
#define ANSI_COLOR_YELLOW  "\x1b[1;33m"
#define ANSI_COLOR_BLUE    "\x1b[1;34m"
#define ANSI_COLOR_MAGENTA "\x1b[1;35m"
#define ANSI_COLOR_CYAN    "\x1b[1;36m"
#define ANSI_COLOR_WHITE   "\e[1;37m"
#define ANSI_COLOR_RESET   "\x1b[0m"

#define NDIM 3

using namespace std;
using namespace Eigen;

int solver(Numerical_CST &numC, bool symY, double sRef, double alpha, Vector3d &vInf, double Minf,
           Network &bPan, Network &wPan, Field &fPan, Subpanel &sp, double &cL, double &cD) {

    //// Begin solver
    cout << ANSI_COLOR_BLUE;
    cout << "*********************" << endl;
    cout << "*Beginning solver...*" << endl;
    cout << "*********************";
    cout << ANSI_COLOR_RESET << endl;

    //// Initialization
    // AIC matrices
    Body_AIC b2bAIC = {}, b2fAIC = {}; // body to body and body to field
    Field2field_AIC f2fAIC = {}; // field to field
    Field2body_AIC f2bAIC = {}; // field to field and field to body
    Subpanel_AIC spAIC = {}; // body to field (sub-panel)

    // Singularities
    bPan.tau.resize(bPan.nP);
    bPan.mu.resize(bPan.nP);
    fPan.sigma = VectorXd::Zero(fPan.nF);
    // Flow variables // TODO check if initialization is useful
    fPan.phi = VectorXd::Zero(fPan.nF);
    fPan.M = VectorXd::Zero(fPan.nF);
    fPan.U = MatrixX3d::Zero(fPan.nF, NDIM);
    fPan.rho = VectorXd::Zero(fPan.nF);
    fPan.dRho = MatrixX3d::Zero(fPan.nF, NDIM);
    fPan.a = VectorXd::Zero(fPan.nF);
    // Numerics
    fPan.epsilon = VectorXd::Zero(fPan.nF);

    // Temporary variables
    int itCnt = 0; // Global iteration counter
    double deltaSigma0 = 0; // Initial sigma change
    VectorXd sigmaTmp, deltaSigma; // To store sigma and delta sigma during iteration
    deltaSigma.resize(fPan.nF);
    sigmaTmp.resize(fPan.nF);
    VectorXd RHS; // Right hand side
    RHS.resize(bPan.nP);

    MatrixX3d vSigma; // Velocity induced by field sources on body
    vSigma = MatrixXd::Zero(bPan.nP, NDIM);

    //// Identify sub-panels
    id_subpanel(bPan, fPan, sp, spAIC);

    //// Build AIC matrices
    build_AIC(symY, bPan, wPan, fPan, b2bAIC, b2fAIC, f2fAIC, f2bAIC, sp, spAIC);

    //// Solver
    // Field Panel Method
    if (Minf != 0) {
        cout << "|¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯|" << endl;
        cout << "|Compressible computation.|" << endl;
        cout << " ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ " << endl << endl;
        do {
            // Store new field sources
            sigmaTmp = fPan.sigma;
            // Panel prediction and boundary condition
            solve_body(vInf, RHS, vSigma, bPan, b2bAIC);
            // Field correction
            solve_field(Minf, vInf, bPan, fPan, sp, b2fAIC, f2fAIC, spAIC);
            // Source induced velocity (to recompute B.C.)
            vSigma.col(0) = f2bAIC.Cu * fPan.sigma;
            vSigma.col(1) = f2bAIC.Cv * fPan.sigma;
            vSigma.col(2) = f2bAIC.Cw * fPan.sigma;
            // Stop criterion
            for (int i = 0; i < fPan.nF; ++i)
                deltaSigma(i) = abs(fPan.sigma(i) - sigmaTmp(i));
            if (!itCnt)
                deltaSigma0 = deltaSigma.norm();
            if (isnan(deltaSigma.norm())) {
                break;
            }
            cout << ANSI_COLOR_CYAN;
            cout << "Relative source change at iteration " << itCnt << ": " << log10(deltaSigma.norm()/deltaSigma0) << endl;
            cout << "FPE global residual at iteration " << itCnt << ": " << log10(fPan.epsilon.norm());
            cout << ANSI_COLOR_RESET << endl << endl;
            itCnt++;
        } while(log10(deltaSigma.norm()/deltaSigma0) > -numC.RRED);
        //TODO Consider using true residual (div(U) - sigma -> 0). Currently impossible since accuracy.
        if (isnan(deltaSigma.norm())) {
            cout << ANSI_COLOR_RED;
            cout << "∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨" << endl;
            cout << ">> Process diverged at iteration #" << itCnt + 1 << "!" << endl;
            cout << "∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧";
            cout << ANSI_COLOR_RESET << endl << endl;
        }
        else {
            cout << ANSI_COLOR_GREEN;
            cout << "∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨" << endl;
            cout << ">> Process converged in " << itCnt << " iteration(s)!" << endl;
            cout << "∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧∧";
            cout << ANSI_COLOR_RESET << endl << endl;
        }
    }
    // Panel Method
    else {
        cout << "|¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯|" << endl;
        cout << "|Incompressible computation.|" << endl;
        cout << " ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ " << endl << endl;
        solve_body(vInf, RHS, vSigma, bPan, b2bAIC);
    }

    // Surface velocity and pressure computation
    compute_sVars(symY, sRef, alpha, Minf, vInf, vSigma, bPan, cL, cD);

    //// End solver
    cout << ANSI_COLOR_BLUE;
    cout << "********************" << endl;
    cout << "*Solver successful!*" << endl;
    cout << "********************";
    cout << ANSI_COLOR_RESET << endl;
    return 0;
}

max_stars_repo_name: acrovato/aero

max_issues_repo_name: acrovato/aero

max_forks_repo_name: acrovato/aero


text:
// Boost.Geometry (aka GGL, Generic Geometry Library)

// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
// Copyright (c) 2013 Adam Wulkiewicz, Lodz, Poland.

// This file was modified by Oracle on 2013, 2014.
// Modifications copyright (c) 2013, 2014 Oracle and/or its affiliates.

// Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
// (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.

// Use, modification and distribution is subject to 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)

// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle

#ifndef BOOST_GEOMETRY_STRATEGY_AGNOSTIC_POINT_IN_POLY_WINDING_HPP
#define BOOST_GEOMETRY_STRATEGY_AGNOSTIC_POINT_IN_POLY_WINDING_HPP

#include <boost/core/ignore_unused.hpp>

#include <boost/geometry/util/math.hpp>
#include <boost/geometry/util/select_calculation_type.hpp>

#include <boost/geometry/strategies/side.hpp>
#include <boost/geometry/strategies/covered_by.hpp>
#include <boost/geometry/strategies/within.hpp>

namespace geofeatures_boost {} namespace boost = geofeatures_boost; namespace geofeatures_boost { namespace geometry
{

namespace strategy { namespace within
{

// Fix for https://svn.boost.org/trac/boost/ticket/9628
// For floating point coordinates, the <1> coordinate of a point is compared
// with the segment's points using some EPS. If the coordinates are "equal"
// the sides are calculated. Therefore we can treat a segment as a long areal
// geometry having some width. There is a small ~triangular area somewhere
// between the segment's effective area and a segment's line used in sides
// calculation where the segment is on the one side of the line but on the
// other side of a segment (due to the width).
// For the s1 of a segment going NE the real side is RIGHT but the point may
// be detected as LEFT, like this:
//                     RIGHT
//                 ___----->
//                  ^      O Pt  __ __
//                 EPS     __ __
//                  v__ __ BUT DETECTED AS LEFT OF THIS LINE
//             _____7
//       _____/
// _____/
template <typename CSTag>
struct winding_side_equal
{
    typedef typename strategy::side::services::default_strategy
        <
            CSTag
        >::type strategy_side_type;

    template <size_t D, typename Point, typename PointOfSegment>
    static inline int apply(Point const& point,
                            PointOfSegment const& se,
                            int count)
    {
        // Create a vertical segment intersecting the original segment's endpoint
        // equal to the point, with the derived direction (UP/DOWN).
        // Set only the 2 first coordinates, the other ones are ignored
        PointOfSegment ss1, ss2;
        set<1-D>(ss1, get<1-D>(se));
        set<1-D>(ss2, get<1-D>(se));
        if (count > 0) // UP
        {
            set<D>(ss1, 0);
            set<D>(ss2, 1);
        }
        else // DOWN
        {
            set<D>(ss1, 1);
            set<D>(ss2, 0);
        }
        // Check the side using this vertical segment
        return strategy_side_type::apply(ss1, ss2, point);
    }
};

// The optimization for cartesian
template <>
struct winding_side_equal<cartesian_tag>
{
    template <size_t D, typename Point, typename PointOfSegment>
    static inline int apply(Point const& point,
                            PointOfSegment const& se,
                            int count)
    {
        return math::equals(get<1-D>(point), get<1-D>(se)) ?
                0 :
                get<1-D>(point) < get<1-D>(se) ?
                    // assuming count is equal to 1 or -1
                    count : // ( count > 0 ? 1 : -1) :
                    -count; // ( count > 0 ? -1 : 1) ;
    }
};

template <typename CSTag>
struct winding_side_between
{
    typedef typename strategy::side::services::default_strategy
        <
            CSTag
        >::type strategy_side_type;

    template <size_t D, typename Point, typename PointOfSegment>
    static inline int apply(Point const& point,
                            PointOfSegment const& s1, PointOfSegment const& s2,
                            int count)
    {
        // Create a vertical segment intersecting the original segment's endpoint
        // equal to the point, with the derived direction (UP/DOWN).
        // Set only the 2 first coordinates, the other ones are ignored
        PointOfSegment ss1, ss2;
        set<1-D>(ss1, get<1-D>(s1));
        set<1-D>(ss2, get<1-D>(s1));

        if (count > 0) // UP
        {
            set<D>(ss1, 0);
            set<D>(ss2, 1);
        }
        else // DOWN
        {
            set<D>(ss1, 1);
            set<D>(ss2, 0);
        }

        int const seg_side = strategy_side_type::apply(ss1, ss2, s2);

        if (seg_side != 0) // segment not vertical
        {
            if (strategy_side_type::apply(ss1, ss2, point) == -seg_side) // point on the opposite side than s2
            {
                return -seg_side;
            }
            else
            {
                set<1-D>(ss1, get<1-D>(s2));
                set<1-D>(ss2, get<1-D>(s2));

                if (strategy_side_type::apply(ss1, ss2, point) == seg_side) // point behind s2
                {
                    return seg_side;
                }
            }
        }

        // segment is vertical or point is between p1 and p2
        return strategy_side_type::apply(s1, s2, point);
    }
};

// The specialization for cartesian
template <>
struct winding_side_between<cartesian_tag>
{
    typedef strategy::side::services::default_strategy
        <
            cartesian_tag
        >::type strategy_side_type;

    template <size_t D, typename Point, typename PointOfSegment>
    static inline int apply(Point const& point,
                            PointOfSegment const& s1, PointOfSegment const& s2,
                            int /*count*/)
    {
        return strategy_side_type::apply(s1, s2, point);
    }
};

/*!
\brief Within detection using winding rule
\ingroup strategies
\tparam Point \tparam_point
\tparam PointOfSegment \tparam_segment_point
\tparam CalculationType \tparam_calculation
\author Barend Gehrels
\note The implementation is inspired by terralib http://www.terralib.org (LGPL)
\note but totally revised afterwards, especially for cases on segments
\note Only dependant on "side", -> agnostic, suitable for spherical/latlong

\qbk{
[heading See also]
[link geometry.reference.algorithms.within.within_3_with_strategy within (with strategy)]
}
 */
template
<
    typename Point,
    typename PointOfSegment = Point,
    typename CalculationType = void
>
class winding
{
    typedef typename select_calculation_type
        <
            Point,
            PointOfSegment,
            CalculationType
        >::type calculation_type;

    typedef typename strategy::side::services::default_strategy
        <
            typename cs_tag<Point>::type
        >::type strategy_side_type;

    /*! subclass to keep state */
    class counter
    {
        int m_count;
        bool m_touches;

        inline int code() const
        {
            return m_touches ? 0 : m_count == 0 ? -1 : 1;
        }

    public :
        friend class winding;

        inline counter()
            : m_count(0)
            , m_touches(false)
        {}

    };

    template <size_t D>
    static inline int check_touch(Point const& point,
                PointOfSegment const& seg1, PointOfSegment const& seg2,
                counter& state)
    {
        calculation_type const p = get<D>(point);
        calculation_type const s1 = get<D>(seg1);
        calculation_type const s2 = get<D>(seg2);
        if ((s1 <= p && s2 >= p) || (s2 <= p && s1 >= p))
        {
            state.m_touches = true;
        }
        return 0;
    }

    template <size_t D>
    static inline int check_segment(Point const& point,
                PointOfSegment const& seg1, PointOfSegment const& seg2,
                counter& state, bool& eq1, bool& eq2)
    {
        calculation_type const p = get<D>(point);
        calculation_type const s1 = get<D>(seg1);
        calculation_type const s2 = get<D>(seg2);

        // Check if one of segment endpoints is at same level of point
        eq1 = math::equals(s1, p);
        eq2 = math::equals(s2, p);

        if (eq1 && eq2)
        {
            // Both equal p -> segment is horizontal (or vertical for D=0)
            // The only thing which has to be done is check if point is ON segment
            return check_touch<1 - D>(point, seg1, seg2, state);
        }

        return
              eq1 ? (s2 > p ?  1 : -1)  // Point on level s1, UP/DOWN depending on s2
            : eq2 ? (s1 > p ? -1 :  1)  // idem
            : s1 < p && s2 > p ?  2     // Point between s1 -> s2 --> UP
            : s2 < p && s1 > p ? -2     // Point between s2 -> s1 --> DOWN
            : 0;
    }

public :

    // Typedefs and static methods to fulfill the concept
    typedef Point point_type;
    typedef PointOfSegment segment_point_type;
    typedef counter state_type;

    static inline bool apply(Point const& point,
                PointOfSegment const& s1, PointOfSegment const& s2,
                counter& state)
    {
        typedef typename cs_tag<Point>::type cs_t;

        bool eq1 = false;
        bool eq2 = false;
        geofeatures_boost::ignore_unused(eq2);

        int count = check_segment<1>(point, s1, s2, state, eq1, eq2);
        if (count != 0)
        {
            int side = 0;
            if (count == 1 || count == -1)
            {
                side = winding_side_equal<cs_t>
                            ::template apply<1>(point, eq1 ? s1 : s2, count);
            }
            else // count == 2 || count == -2
            {
                side = winding_side_between<cs_t>
                            ::template apply<1>(point, s1, s2, count);
            }

            if (side == 0)
            {
                // Point is lying on segment
                state.m_touches = true;
                state.m_count = 0;
                return false;
            }

            // Side is NEG for right, POS for left.
            // The count is -2 for down, 2 for up (or -1/1)
            // Side positive thus means UP and LEFTSIDE or DOWN and RIGHTSIDE
            // See accompagnying figure (TODO)
            if (side * count > 0)
            {
                state.m_count += count;
            }
        }
        return ! state.m_touches;
    }

    static inline int result(counter const& state)
    {
        return state.code();
    }
};

#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS

namespace services
{

// Register using "areal_tag" for ring, polygon, multi-polygon
template <typename AnyTag, typename Point, typename Geometry>
struct default_strategy<point_tag, AnyTag, point_tag, areal_tag, cartesian_tag, cartesian_tag, Point, Geometry>
{
    typedef winding<Point, typename geometry::point_type<Geometry>::type> type;
};

template <typename AnyTag, typename Point, typename Geometry>
struct default_strategy<point_tag, AnyTag, point_tag, areal_tag, spherical_tag, spherical_tag, Point, Geometry>
{
    typedef winding<Point, typename geometry::point_type<Geometry>::type> type;
};

// TODO: use linear_tag and pointlike_tag the same way how areal_tag is used

template <typename Point, typename Geometry, typename AnyTag>
struct default_strategy<point_tag, AnyTag, point_tag, AnyTag, cartesian_tag, cartesian_tag, Point, Geometry>
{
    typedef winding<Point, typename geometry::point_type<Geometry>::type> type;
};

template <typename Point, typename Geometry, typename AnyTag>
struct default_strategy<point_tag, AnyTag, point_tag, AnyTag, spherical_tag, spherical_tag, Point, Geometry>
{
    typedef winding<Point, typename geometry::point_type<Geometry>::type> type;
};

} // namespace services

#endif

}} // namespace strategy::within

#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
namespace strategy { namespace covered_by { namespace services
{

// Register using "areal_tag" for ring, polygon, multi-polygon
template <typename AnyTag, typename Point, typename Geometry>
struct default_strategy<point_tag, AnyTag, point_tag, areal_tag, cartesian_tag, cartesian_tag, Point, Geometry>
{
    typedef strategy::within::winding<Point, typename geometry::point_type<Geometry>::type> type;
};

template <typename AnyTag, typename Point, typename Geometry>
struct default_strategy<point_tag, AnyTag, point_tag, areal_tag, spherical_tag, spherical_tag, Point, Geometry>
{
    typedef strategy::within::winding<Point, typename geometry::point_type<Geometry>::type> type;
};

// TODO: use linear_tag and pointlike_tag the same way how areal_tag is used

template <typename Point, typename Geometry, typename AnyTag>
struct default_strategy<point_tag, AnyTag, point_tag, AnyTag, cartesian_tag, cartesian_tag, Point, Geometry>
{
    typedef strategy::within::winding<Point, typename geometry::point_type<Geometry>::type> type;
};

template <typename Point, typename Geometry, typename AnyTag>
struct default_strategy<point_tag, AnyTag, point_tag, AnyTag, spherical_tag, spherical_tag, Point, Geometry>
{
    typedef strategy::within::winding<Point, typename geometry::point_type<Geometry>::type> type;
};

}}} // namespace strategy::covered_by::services
#endif

}} // namespace geofeatures_boost::geometry

#endif // BOOST_GEOMETRY_STRATEGY_AGNOSTIC_POINT_IN_POLY_WINDING_HPP

max_stars_repo_name: xarvey/Yuuuuuge

max_issues_repo_name: xarvey/Yuuuuuge

max_forks_repo_name: xarvey/Yuuuuuge