text:
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*  This file is part of the library KASKADE 7                               */
/*    see http://www.zib.de/projects/kaskade7-finite-element-toolbox         */
/*                                                                           */
/*  Copyright (C) 2002-2009 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

#include <complex>
#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <cmath>

#include <boost/timer.hpp>

#include "dune/common/stdstreams.hh"
#include "dune/grid/sgrid.hh"
#include "dune/grid/uggrid.hh"
#include "dune/grid/geometrygrid.hh"
#include "dune/grid/common/gridinfo.hh"
#include "dune/istl/solvers.hh"

#include "fem/assemble.hh"
#include "fem/embedded_errorest.hh"
#include "fem/istlinterface.hh"
#include "fem/functional_aux.hh"
#include "fem/hierarchicspace.hh"
#include "fem/lagrangespace.hh"
#include "linalg/trivialpreconditioner.hh"
#include "linalg/direct.hh"
#include "linalg/triplet.hh"
#include "linalg/iluprecond.hh"
#include "linalg/additiveschwarz.hh"
#include "linalg/hyprecond.hh"
#include "io/vtk.hh"
#include "io/amira.hh"

#include "utilities/kaskopt.hh"

#include "ht.hh"

 namespace Dune
  {

    namespace Capabilities
    {

      template< class Grid >
      struct hasHierarchicIndexSet
      {
        static const bool v = false;
      };

      template< class Grid >
      struct hasHierarchicIndexSet< const Grid >
      {
        static const bool v = hasHierarchicIndexSet< Grid >::v;
      };

    }

  }

class SquareToCircle
  : public Dune::AnalyticalCoordFunction< double, 2, 2, SquareToCircle >
  {
    typedef SquareToCircle This;
    typedef Dune::AnalyticalCoordFunction< double, 2, 2, This > Base;

  public:
    typedef Base::DomainVector DomainVector;
    typedef Base::RangeVector RangeVector;

    void evaluate ( const DomainVector &x, RangeVector &y ) const
    {
      double enorm = sqrt(x[0]*x[0]+x[1]*x[1]);

      if(enorm > 0.00001)
      {

        double radius = sqrt(2.0)*std::max(std::fabs(x[0]),std::fabs(x[1]));

      double scaling = radius/enorm;

      y[ 0 ] = x[ 0 ]*scaling;
      y[ 1 ] = x[ 1 ]*scaling;
      } else
      {
      y[ 0 ] = x[ 0 ];
      y[ 1 ] = x[ 1 ];
      }

    }
  };

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

	//   two-dimensional space: dim=2
	int const dim=2;
	typedef Dune::UGGrid<dim> Grid;

        typedef Dune::GeometryGrid<Grid,SquareToCircle> GGrid;

	Dune::GridFactory<Grid> factory;

	// vertex coordinates v[0], v[1]
	Dune::FieldVector<double,dim> v;
	v[0]=-1; v[1]=-1; factory.insertVertex(v);
	v[0]=1; v[1]=-1; factory.insertVertex(v);
	v[0]=1; v[1]=1; factory.insertVertex(v);
	v[0]=-1; v[1]=1; factory.insertVertex(v);
	v[0]=0; v[1]=0; factory.insertVertex(v);
	// triangle defined by 3 vertex indices
	std::vector<unsigned int> vid(3);
	Dune::GeometryType gt(Dune::GeometryType::simplex,2);
	vid[0]=0; vid[1]=1; vid[2]=4; factory.insertElement(gt,vid);
	vid[0]=1; vid[1]=2; vid[2]=4; factory.insertElement(gt,vid);
	vid[0]=2; vid[1]=3; vid[2]=4; factory.insertElement(gt,vid);
	vid[0]=3; vid[1]=0; vid[2]=4; factory.insertElement(gt,vid);
	std::auto_ptr<Grid> grid( factory.createGrid() ) ;

	// the coarse grid will be refined three times
	// some information on the refined mesh
	std::cout << "Grid: " << grid->size(0) << " triangles, " << std::endl;
	std::cout << "      " << grid->size(1) << " edges, " << std::endl;
	std::cout << "      " << grid->size(2) << " points" << std::endl;
	// a gridmanager is constructed
	// as connector between geometric and algebraic information

        SquareToCircle def;

        std::auto_ptr<GGrid> ggrid(new GGrid(*grid,def));

	ggrid->globalRefine(6);

	GridManager<GGrid> gridManager(ggrid);

//StartSnippet2
	// construction of finite element space for the scalar solution T
	typedef GGrid::LeafGridView GridView;
	typedef FEFunctionSpace<ContinuousLagrangeMapper<double,GridView> > H1Space;
	H1Space temperatureSpace(gridManager,gridManager.grid().leafView(),
							 1);
	typedef boost::fusion::vector<H1Space const*> Spaces;
	Spaces spaces(&temperatureSpace);
	// VariableDescription<int spaceId, int components, int Id>
	// spaceId: number of associated FEFunctionSpace
	// components: number of components in this variable
	// Id: number of this variable
	typedef boost::fusion::vector<VariableDescription<0,1,0> >
			VariableDescriptions;
	std::string varNames[1] = { "T" };
	typedef VariableSetDescription<Spaces,VariableDescriptions> VariableSet;
	VariableSet variableSet(spaces,varNames);
//StopSnippet2

//StartSnippet3
	typedef HeatFunctional<double,VariableSet> Functional;
	Functional F;
	typedef VariationalFunctionalAssembler<LinearizationAt<Functional> > GOP;
	typedef GOP::TestVariableRepresentation<>::type Rhs;
	GOP gop(gridManager.signals,spaces);
	VariableSet::VariableSet x(variableSet);
//StopSnippet3

//StartSnippet4

	typedef GOP::AnsatzVariableRepresentation<>::type Sol;
	Sol solution(GOP::AnsatzVariableRepresentation<>::init(gop));
	solution = 0;
	gop.assemble(linearization(F,x));

	Rhs rhs(GOP::TestVariableRepresentation<>::rhs(gop));
	AssembledGalerkinOperator<GOP,0,1,0,1> A(gop, false);
	AssembledGalerkinOperator<GOP,0,1,0,1>::matrix_type tri(A.getmat());
//	A.getmat().print();

//StartSnippet5
    boost::timer directTimer;
    directInverseOperator(A,DirectType::UMFPACK,MatrixProperties::GENERAL).applyscaleadd(-1.0,rhs,solution);
    std::cout << "direct solve: " << directTimer.elapsed() << "s\n";
    x.data = solution.data;
//StopSnippet5

//StartSnippet6
	// output of solution in VTK format for visualization,
	// the data are written as ascii stream into file temperature.vtu,
	// possible is also binary
 	IoOptions options;
 	options.outputType = IoOptions::ascii;

    typedef GGrid::LeafGridView LeafGridView;
    LeafGridView leafView = gridManager.grid().leafView();
    writeVTKFile(leafView,variableSet,x,"temperature",options);

  }

max_stars_repo_name: chenzongxiong/streambox

max_issues_repo_name: chenzongxiong/streambox

max_forks_repo_name: chenzongxiong/streambox


text:
#include <iostream>
#include <time.h>
#include <pcl/io/pcd_io.h>
#include <pcl/point_types.h>
#include <pcl/registration/icp.h>
#include <pcl/registration/ia_ransac.h>
#include <pcl/point_types.h>
#include <pcl/features/normal_3d.h>
#include <pcl/point_types.h>
#include <pcl/features/pfh.h>
#include <pcl/filters/passthrough.h>
#include <pcl/visualization/cloud_viewer.h>
#include <limits>
#include <fstream>
#include <vector>
#include <Eigen/Core>
#include "pcl/point_cloud.h"
#include "pcl/kdtree/kdtree_flann.h"
#include "pcl/filters/passthrough.h"
#include "pcl/filters/voxel_grid.h"
#include "pcl/features/fpfh.h"

using namespace pcl;
using namespace std;

namespace pc{

void voxelFilter(pcl::PointCloud<pcl::PointXYZ>::Ptr &cloud_in, pcl::PointCloud<pcl::PointXYZ>::Ptr &cloud_out, float gridsize){
	VoxelGrid<PointXYZ> vox_grid;
	vox_grid.setLeafSize(gridsize, gridsize, gridsize);
	vox_grid.setInputCloud(cloud_in);
	vox_grid.filter(*cloud_out);
	return;
}

void passFilter(pcl::PointCloud<pcl::PointXYZ>::Ptr pc, double filter_limit){
	pcl::PassThrough<pcl::PointXYZ> pass;
	pass.setInputCloud(pc);
	pass.setFilterFieldName("x");
	pass.setFilterLimits(0, filter_limit);
	pass.setFilterFieldName("y");
	pass.setFilterLimits(0, filter_limit);
	pass.setFilterFieldName("z");
	pass.setFilterLimits(0, filter_limit);
	pass.filter(*pc);
	return;
}

}

ext: cpp

max_stars_repo_name: Ewenwan/RVAF

max_issues_repo_name: Ewenwan/RVAF


text:
#include <iostream>
#include <cmath>
#include <string>
#include <random>
#include <chrono>
#include <fstream>
#include <tuple>
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <unsupported/Eigen/SparseExtra>
#include "block_iterative_solvers.hpp"
#include "linear_algebra_addon.hpp"
using namespace std;
using namespace Eigen;
using namespace std::chrono;

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

    double tol= 1e-8;
    int itermax= 1000;

    std::srand((unsigned int) time(0));
    int L=atoi(argv[1]); // dimension of matrix A

    MatrixXcd A;

    switch (argc)
    {
    case 2:
        {
        // load A from a file in matrix market format .mtx
        cout << "please enter filename for A (.mtx): " ;
        string matfname;
        cin >> matfname;
        SparseMatrix<complex<double>> mat;
        //Matrix collections is available in   http://www.cise.ufl.edu/research/sparse/matrices/
        if(!loadMarket(mat,matfname)){
            cerr << "can't open file" << matfname << endl;
            exit(EXIT_FAILURE);
        }
        cout << "A was loaded from " <<  matfname << endl;
        A= MatrixXcd(mat);

        break;
        }
    case 3:
        {
        // Set A as a random matrix
        int n=atoi(argv[2]); // number of RHS vector
        A= MatrixXcd::Random(n,n);
        A= A.transpose()*A; // complex-symmetric
        cout << "A is a " << n << " x " << n << " square random matrix" << endl;
        break;
        }
    default:
        {

        break;
        }
    }

    MatrixXcd B= MatrixXcd::Random(A.rows(),L);

    {
        time_point<steady_clock> start= steady_clock::now();
        cout << "solving AX=B using bl_cocg_rq ..." << endl;
            MatrixXcd X= bl_cocg_rq(A, B, tol, itermax); // only applicable to complex-symmetric matrix A
        cout << " bl_cocg_rq relative error: " << (A*X-B).norm()/B.norm() << endl << endl;
        auto end= steady_clock::now();
        cout << "bl_cocg_rq computation time= " << duration_cast<seconds>((end - start)).count()  << " sec" << endl << endl<< endl;
    }

    {
        time_point<steady_clock> start= steady_clock::now();
        cout << "solving AX=B using bl_bicg_rq ..." << endl;
            MatrixXcd X= bl_bicg_rq(A, B, tol, itermax); // applicable to general matrix A but needs A.adjoint()
        cout << " bl_bicg_rq relative error: " << (A*X-B).norm()/B.norm() << endl << endl;
        auto end= steady_clock::now();
        cout << "bl_bicg_rq computation time= " << duration_cast<seconds>((end - start)).count() << " sec" << endl << endl<< endl;
    }

    {
        time_point<steady_clock> start= steady_clock::now();
        cout << "solving AX=B using bl_bicr_rq ..." << endl;
            MatrixXcd X= bl_bicr_rq(A, B, tol, itermax); // applicable to general matrix A but needs A.adjoint()
        cout << " bl_bicr_rq relative error: " << (A*X-B).norm()/B.norm() << endl << endl;
        auto end= steady_clock::now();
        cout << "bl_bicr_rq computation time= " << duration_cast<seconds>((end - start)).count() << " sec" << endl << endl<< endl;
    }

    {
        time_point<steady_clock> start= steady_clock::now();
        cout << "solving AX=B using bl_bicgstab_rq ..." << endl;
            MatrixXcd X= bl_bicgstab_rq(A, B, tol, itermax); // only applicable to complex-symmetric matrix A
        cout << " bl_bicgstab_rq relative error: " << (A*X-B).norm()/B.norm() << endl << endl;
        auto end= steady_clock::now();
        cout << "bl_bicgstab_rq computation time= " << duration_cast<seconds>((end - start)).count() << " sec" << endl << endl<< endl;
    }

    {
        time_point<steady_clock> start= steady_clock::now();
        cout << "solving AX=B using bl_bicgstab ..." << endl;
            MatrixXcd X= bl_bicgstab(A, B, tol, itermax); // only applicable to complex-symmetric matrix A
        cout << " bl_bicgstab relative error: " << (A*X-B).norm()/B.norm() << endl << endl;
        auto end= steady_clock::now();
        cout << "bl_bicgstab computation time= " << duration_cast<seconds>((end - start)).count() << " sec" << endl << endl<< endl;
    }

    {
        time_point<steady_clock> start= steady_clock::now();
        cout << "solving AX=B using bl_bicggr ..." << endl;
            MatrixXcd X= bl_bicggr(A, B, tol, itermax); // only applicable to complex-symmetric matrix A
        cout << " bl_bicggr relative error: " << (A*X-B).norm()/B.norm() << endl << endl;
        auto end= steady_clock::now();
        cout << "bl_bicggr computation time= " << duration_cast<seconds>((end - start)).count() << " sec" << endl << endl<< endl;
    }

}

max_stars_repo_name: nmoteki/block-Krylov-linear-solvers

max_issues_repo_name: nmoteki/block-Krylov-linear-solvers

max_forks_repo_name: nmoteki/block-Krylov-linear-solvers


text: #include <boost/units/systems/si/plane_angle.hpp>

ext: hpp

max_stars_repo_name: miathedev/BoostForArduino

max_issues_repo_name: miathedev/BoostForArduino

max_forks_repo_name: miathedev/BoostForArduino


text:
#include "fbstab/fbstab_dense.h"

#include <Eigen/Dense>
#include <memory>
#include <stdexcept>

#include "fbstab/components/dense_cholesky_solver.h"
#include "fbstab/components/dense_data.h"
#include "fbstab/components/full_feasibility.h"
#include "fbstab/components/full_residual.h"
#include "fbstab/components/full_variable.h"
#include "fbstab/fbstab_algorithm.h"
#include "tools/output_stream.h"
#include "tools/utilities.h"

namespace fbstab {

FBstabDense::FBstabDense(int nz, int nl, int nv) {
  if (nz <= 0 || nv <= 0 || nl < 0) {
    throw std::runtime_error(
        "In FBstabDense::FBstabDense: nz and nv must be positive, nl "
        "nonnegative");
  }
  nz_ = nz;
  nl_ = nl;
  nv_ = nv;

  x1_ = tools::make_unique<FullVariable>(nz_, nl_, nv_);
  x2_ = tools::make_unique<FullVariable>(nz_, nl_, nv_);
  x3_ = tools::make_unique<FullVariable>(nz_, nl_, nv_);
  x4_ = tools::make_unique<FullVariable>(nz_, nl_, nv_);
  r1_ = tools::make_unique<FullResidual>(nz_, nl_, nv_);
  r2_ = tools::make_unique<FullResidual>(nz_, nl_, nv_);
  feasibility_checker_ = tools::make_unique<FullFeasibility>(nz_, nl_, nv_);
  linear_solver_ = tools::make_unique<DenseCholeskySolver>(nz_, nl_, nv_);

  algorithm_ = tools::make_unique<FBstabAlgoDense>(
      x1_.get(), x2_.get(), x3_.get(), x4_.get(), r1_.get(), r2_.get(),
      linear_solver_.get(), feasibility_checker_.get());

  opts_ = DefaultOptions();
}

void FBstabDense::UpdateOptions(const Options& options) {
  // No need to validate since there are no additional options and the algorithm
  // will check the algorithmic ones.
  algorithm_->UpdateParameters(&options);
}

FBstabDense::Options FBstabDense::DefaultOptions() {
  Options opts;
  opts.DefaultParameters();
  return opts;
}

FBstabDense::Options FBstabDense::ReliableOptions() {
  Options opts;
  opts.ReliableParameters();
  return opts;
}

FBstabDense::ProblemData::ProblemData(int nz, int nl, int nv) {
  H.resize(nz, nz);
  G.resize(nl, nz);
  A.resize(nv, nz);
  f.resize(nz);
  h.resize(nl);
  b.resize(nv);
}

FBstabDense::ProblemDataRef::ProblemDataRef(
    const Eigen::Map<Eigen::MatrixXd>* H_,
    const Eigen::Map<Eigen::VectorXd>* f_,
    const Eigen::Map<Eigen::MatrixXd>* G_,
    const Eigen::Map<Eigen::VectorXd>* h_,
    const Eigen::Map<Eigen::MatrixXd>* A_,
    const Eigen::Map<Eigen::VectorXd>* b_)
    : H(H_->data(), H_->rows(), H_->cols()),
      G(G_->data(), G_->rows(), G_->cols()),
      A(A_->data(), A_->rows(), A_->cols()),
      f(f_->data(), f_->size()),
      h(h_->data(), h_->size()),
      b(b_->data(), b_->size()) {}

FBstabDense::Variable::Variable(int nz, int nl, int nv) {
  z = Eigen::VectorXd::Zero(nz);
  l = Eigen::VectorXd::Zero(nl);
  v = Eigen::VectorXd::Zero(nv);
  y = Eigen::VectorXd::Zero(nv);
}

FBstabDense::VariableRef::VariableRef(Eigen::Map<Eigen::VectorXd>* z_,
                                      Eigen::Map<Eigen::VectorXd>* l_,
                                      Eigen::Map<Eigen::VectorXd>* v_,
                                      Eigen::Map<Eigen::VectorXd>* y_)
    : z(z_->data(), z_->size()),
      l(l_->data(), l_->size()),
      v(v_->data(), v_->size()),
      y(y_->data(), y_->size()) {}

void FBstabDense::VariableRef::fill(double a) {
  z.fill(a);
  l.fill(a);
  v.fill(a);
  y.fill(a);
}

// Explicit instantiation.
template class FBstabAlgorithm<FullVariable, FullResidual, DenseCholeskySolver,
                               FullFeasibility>;

}  // namespace fbstab

ext: cc

max_stars_repo_name: tcunis/fbstab

max_issues_repo_name: tcunis/fbstab

max_forks_repo_name: tcunis/fbstab


text:
#include "cpputil/StringHelper.h"

#include <boost/locale/encoding.hpp>

#include <stdarg.h>

namespace cpputil
{

void StringHelper::Split(const std::string& src, const std::string& mid, std::vector<std::string>& dst)
{
	char* cstr = new char[src.length() + 1];
	std::strcpy(cstr, src.c_str());
	char* p = std::strtok(cstr, mid.c_str());
	while (p != 0)
	{
		dst.push_back(std::string(p));
		p = strtok(nullptr, mid.c_str());
	}
	delete[] cstr;
}

void StringHelper::ReplaceAll(std::string& str, const std::string& from, const std::string& to)
{
	if (from.empty()) {
		return;
	}
	size_t start_pos = 0;
	while ((start_pos = str.find(from, start_pos)) != std::string::npos) {
		str.replace(start_pos, from.length(), to);
		start_pos += to.length();
	}
}

std::string StringHelper::UTF8ToGBK(const char* str)
{
	return boost::locale::conv::from_utf(str, "GBK");
}

std::string StringHelper::GBKToUTF8(const char* str)
{
	return boost::locale::conv::to_utf<char>(str, "GBK");
}

std::string StringHelper::Format(const std::string fmt, ...)
{
	int size = ((int)fmt.size()) * 2 + 50;   // Use a rubric appropriate for your code
	std::string str;
	va_list ap;
	while (1) {     // Maximum two passes on a POSIX system...
		str.resize(size);
		va_start(ap, fmt);
		int n = vsnprintf((char *)str.data(), size, fmt.c_str(), ap);
		va_end(ap);
		if (n > -1 && n < size) {  // Everything worked
			str.resize(n);
			return str;
		}
		if (n > -1)  // Needed size returned
			size = n + 1;   // For null char
		else
			size *= 2;      // Guess at a larger size (OS specific)
	}
	return str;
}

std::string StringHelper::Trim(const std::string& str)
{
	size_t first = str.find_first_not_of(' ');
	if (std::string::npos == first) {
		return str;
	}
	size_t last = str.find_last_not_of(' ');
	return str.substr(first, (last - first + 1));
}

}

max_stars_repo_name: xzrunner/cpputil

max_issues_repo_name: xzrunner/cpputil

max_forks_repo_name: xzrunner/cpputil


text:
//
// Copyright (c) 2011-2012 Krzysztof Jusiak (krzysztof at jusiak dot net)
//
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
//
#ifndef QFSM_BACK_AUX_DETAIL_RESULTOF_HPP
#define QFSM_BACK_AUX_DETAIL_RESULTOF_HPP

#include <boost/preprocessor/punctuation/comma_if.hpp>
#include <boost/preprocessor/repetition/repeat.hpp>
#include <boost/preprocessor/repetition/enum_params.hpp>

#ifndef RESULT_OF_LIMIT_SIZE
# define RESULT_OF_LIMIT_SIZE 5
#endif

namespace QFsm
{
namespace Back
{
namespace Aux
{
namespace Detail
{

/**
 * helper class to identify result type without decltype
 *
 * @code
 * QFsm::Back::Aux::Detail::ResultOf<int (*)(int)>::type result_type; // -> int
 * @endcode
 */
template<typename> struct ResultOf;

#define RESULT_OF_IMPL(z, n, text)\
    template<typename Ret BOOST_PP_COMMA_IF(n) BOOST_PP_ENUM_PARAMS(n, typename T)> struct ResultOf<Ret (*)(BOOST_PP_ENUM_PARAMS(n, T))> { typedef Ret type; };

#define RESULT_OF_MEMBER_IMPL(z, n, text)\
    template<typename Ret, typename Base BOOST_PP_COMMA_IF(n) BOOST_PP_ENUM_PARAMS(n, typename T)> struct ResultOf<Ret (Base::*)(BOOST_PP_ENUM_PARAMS(n, T))> { typedef Ret type; };

BOOST_PP_REPEAT(RESULT_OF_LIMIT_SIZE, RESULT_OF_IMPL, ~)
BOOST_PP_REPEAT(RESULT_OF_LIMIT_SIZE, RESULT_OF_MEMBER_IMPL, ~)

#undef RESULT_OF_IMPL
#undef RESULT_OF_MEMBER_IMPL

} // namespace Detail
} // namespace Aux
} // namespace Back
} // namespace QFsm

#endif

max_stars_repo_name: krzysztof-jusiak/qfsm

max_issues_repo_name: krzysztof-jusiak/qfsm

max_forks_repo_name: krzysztof-jusiak/qfsm


text:
/*==============================================================================
    Copyright (c) 2001-2010 Joel de Guzman
    Copyright (c) 2010 Thomas Heller

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

#include <boost/phoenix/core/limits.hpp>
#include <boost/phoenix/core/call.hpp>
#include <boost/phoenix/core/expression.hpp>
#include <boost/phoenix/core/meta_grammar.hpp>

BOOST_PHOENIX_DEFINE_EXPRESSION(
    (boost)(phoenix)(do_while)
  , (meta_grammar) // Cond
    (meta_grammar) // Do
)

namespace boost { namespace phoenix
{
    struct do_while_eval
    {
        typedef void result_type;

        template <typename Cond, typename Do, typename Context>
        result_type
        operator()(Cond const& cond, Do const& do_it, Context const & ctx) const
        {
            do
                boost::phoenix::eval(do_it, ctx);
            while (boost::phoenix::eval(cond, ctx));
        }
    };

    template <typename Dummy>
    struct default_actions::when<rule::do_while, Dummy>
        : call<do_while_eval, Dummy>
    {};

    template <typename Do>
    struct do_while_gen
    {
        do_while_gen(Do const& do_it)
            : do_(do_it) {}

        template <typename Cond>
        typename expression::do_while<Cond, Do>::type const
        while_(Cond const& cond) const
        {
            return expression::do_while<Cond, Do>::make(cond, do_);
        }

        Do const& do_;
    };

    struct do_gen
    {
        template <typename Do>
        do_while_gen<Do> const
        operator[](Do const& do_) const
        {
            return do_while_gen<Do>(do_);
        }
    };

#ifndef BOOST_PHOENIX_NO_PREDEFINED_TERMINALS
    do_gen const do_ = {};
#endif

}}

#endif


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

#include "lvt_local_map.h"
#include "lvt_image_features_handler.h"
#include "lvt_image_features_struct.h"
#include "lvt_visualization.h"
#include "lvt_logging_utils.h"
#include <cmath>
#include <Eigen/SVD>

#ifdef LVT_ENABLE_VISUALIZATION

#define VIS_BEGIN_ADDING_MATCHES()               \
    if (m_vo_params.enable_visualization)        \
    {                                            \
        m_visualization->begin_adding_matches(); \
    }

#define VIS_ADD_MATCH(_mp, _ip)               \
    if (m_vo_params.enable_visualization)     \
    {                                         \
        m_visualization->add_match(_mp, _ip); \
    }

#define VIS_FINISH_ADDING_MATCHES()               \
    if (m_vo_params.enable_visualization)         \
    {                                             \
        m_visualization->finish_adding_matches(); \
    }

#else
#define VIS_BEGIN_ADDING_MATCHES() ((void)0)
#define VIS_ADD_MATCH(_mp, _ip) ((void)0)
#define VIS_FINISH_ADDING_MATCHES() ((void)0)
#endif //LVT_ENABLE_VISUALIZATION

static float min_x, max_x, min_y, max_y;

static inline bool is_point_visible(const lvt_vector3 &pt, const lvt_matrix34 &w2c, const lvt_parameters &params,
                                    lvt_vector2 &out_projected_pt)
{
    lvt_vector4 pt_h;
    pt_h << pt, 1.0;
    lvt_vector3 pt_cam = w2c * pt_h;
    if (pt_cam.z() < params.near_plane_distance || pt_cam.z() > params.far_plane_distance)
    {
        return false;
    }
    const double inv_z = 1.0 / pt_cam.z();
    const double u = params.fx * pt_cam.x() * inv_z + params.cx;
    const double v = params.fy * pt_cam.y() * inv_z + params.cy;
    if (u < min_x || u > max_x ||
        v < min_y || v > max_y)
    {
        return false;
    }
    out_projected_pt << u, v;
    return true;
}

lvt_local_map::lvt_local_map(const lvt_parameters &vo_params,
                             lvt_image_features_handler *features_handler, lvt_visualization *vis) : m_vo_params(vo_params), m_features_handler(features_handler), m_visualization(vis)
{
    // compute image bounds
    if (fabs(m_vo_params.k1) < 1e-5)
    {
        min_x = 0.0;
        max_x = m_vo_params.img_width;
        min_y = 0.0;
        max_y = m_vo_params.img_height;
    }
    else
    {
        cv::Mat kps_mat(4, 2, CV_32F);
        kps_mat.at<float>(0, 0) = 0.0;
        kps_mat.at<float>(0, 1) = 0.0;
        kps_mat.at<float>(1, 0) = m_vo_params.img_width;
        kps_mat.at<float>(1, 1) = 0.0;
        kps_mat.at<float>(2, 0) = 0.0;
        kps_mat.at<float>(2, 1) = m_vo_params.img_height;
        kps_mat.at<float>(3, 0) = m_vo_params.img_width;
        kps_mat.at<float>(3, 1) = m_vo_params.img_height;
        cv::Matx33f intrinsics_mtrx(m_vo_params.fx, 0.0, m_vo_params.cx,
                                    0.0, m_vo_params.fy, m_vo_params.cy,
                                    0.0, 0.0, 1.0);
        std::vector<float> dist;
        dist.push_back(m_vo_params.k1);
        dist.push_back(m_vo_params.k2);
        dist.push_back(m_vo_params.p1);
        dist.push_back(m_vo_params.p2);
        dist.push_back(m_vo_params.k3);
        kps_mat = kps_mat.reshape(2);
        cv::undistortPoints(kps_mat, kps_mat, cv::Mat(intrinsics_mtrx), cv::Mat(dist), cv::Mat(), intrinsics_mtrx);
        kps_mat = kps_mat.reshape(1);
        min_x = std::min(kps_mat.at<float>(0, 0), kps_mat.at<float>(2, 0));
        max_x = std::max(kps_mat.at<float>(1, 0), kps_mat.at<float>(3, 0));
        min_y = std::min(kps_mat.at<float>(0, 1), kps_mat.at<float>(1, 1));
        max_y = std::max(kps_mat.at<float>(2, 1), kps_mat.at<float>(3, 1));
    }
}

lvt_local_map::~lvt_local_map()
{
    reset();
}

void lvt_local_map::reset()
{
    m_map_points.clear();
    m_staged_points.clear();
}

int lvt_local_map::find_matches(const lvt_pose &cam_pose, lvt_image_features_struct *left_struct,
                                lvt_vector3_array *out_map_points, std::vector<int> *out_matches_left)
{
    const lvt_matrix34 cml = lvt_pose_utils::compute_world_to_camera_transform(cam_pose);
    int matches_count = 0;
    VIS_BEGIN_ADDING_MATCHES();
    std::vector<int> matches(m_map_points.size(), -2); // mark each map point with its matching index from left_struct, -2 if not visible
    lvt_vector2_array projections(m_map_points.size());
#ifdef LVT_ENABLE_MEASURMENT_RECORD
    std::vector<float> d1s(m_map_points.size());
    std::vector<float> d2s(m_map_points.size());
#endif //LVT_ENABLE_MEASURMENT_RECORD

    for (int i = 0; i < m_map_points.size(); i++)
    {
        lvt_vector2 proj_pt_left;
        if (!is_point_visible(m_map_points[i].m_position, cml, m_vo_params, proj_pt_left))
        {
            m_map_points[i].m_counter += 1;
            matches[i] = -2;
            continue;
        }
        projections[i] = proj_pt_left;
        float d1, d2;
        int match_idx_left = left_struct->find_match_index(proj_pt_left, m_map_points[i].m_descriptor, &d1, &d2);
        matches[i] = match_idx_left;
        if (match_idx_left != -1)
        {
            matches_count++;
            left_struct->mark_as_matched(match_idx_left, true);
#ifdef LVT_ENABLE_MEASURMENT_RECORD
            d1s[i] = d1;
            d2s[i] = d2;
#endif //LVT_ENABLE_MEASURMENT_RECORD
        }
    }

    if (matches_count < LVT_N_MATCHES_TH)
    {
        matches_count = 0;
        left_struct->reset_matched_marks();
        int original_tracking_radius = left_struct->get_tracking_radius();
        left_struct->set_tracking_radius(2 * original_tracking_radius);
        for (int i = 0; i < m_map_points.size(); i++)
        {
            if (matches[i] == -2)
            {
                continue;
            }
            float d1, d2;
            int match_idx_left = left_struct->find_match_index(projections[i], m_map_points[i].m_descriptor, &d1, &d2);
            matches[i] = match_idx_left;
            if (match_idx_left != -1)
            {
                matches_count++;
                left_struct->mark_as_matched(match_idx_left, true);
#ifdef LVT_ENABLE_MEASURMENT_RECORD
                d1s[i] = d1;
                d2s[i] = d2;
#endif //LVT_ENABLE_MEASURMENT_RECORD
            }
        }
        left_struct->set_tracking_radius(original_tracking_radius);
    }

    for (int i = 0; i < matches.size(); i++)
    {
        m_map_points[i].m_match_idx = matches[i];
        if (matches[i] == -2)
        {
            continue;
        }

        if (matches[i] == -1)
        {
            m_map_points[i].m_counter += 1;
            continue;
        }

        m_map_points[i].m_age += 1;
        out_map_points->push_back(m_map_points[i].m_position);
        out_matches_left->push_back(matches[i]);
        LVT_RECORD("img feature x", left_struct->get_keypoint(matches[i]).pt.x);
        LVT_RECORD("img feature y", left_struct->get_keypoint(matches[i]).pt.y);
        LVT_RECORD("age", m_map_points[i].m_age);
        LVT_RECORD("closest descriptor distance", d1s[i]);
        LVT_RECORD("second descriptor distance", d2s[i]);
        VIS_ADD_MATCH(i, matches[i]);
    }

    VIS_FINISH_ADDING_MATCHES();
    LVT_RECORD("tracked map points", matches_count);
    return matches_count;
}

void lvt_local_map::triangulate_rgbd(const lvt_pose &cam_pose, lvt_image_features_struct *img_struct, lvt_map_point_array *out_points)
{
    const float inv_fx = 1.0f / m_vo_params.fx;
    const float inv_fy = 1.0f / m_vo_params.fy;
    lvt_matrix34 cam_to_world_mtrx;
    cam_to_world_mtrx << cam_pose.get_orientation_matrix(), cam_pose.get_position();
    for (int i = 0, count = img_struct->get_features_count(); i < count; i++)
    {
        const cv::Point2f pt = img_struct->get_keypoint(i).pt;
        const float u = pt.x;
        const float v = pt.y;
        const float z = img_struct->get_keypoint_depth(i);
        const float x = (u - m_vo_params.cx) * z * inv_fx;
        const float y = (v - m_vo_params.cy) * z * inv_fy;
        lvt_vector4 pt_h;
        pt_h << x, y, z, 1.0;
        lvt_vector3 pt_w = cam_to_world_mtrx * pt_h;

        lvt_map_point mp;
        mp.m_position = pt_w;
        mp.m_descriptor = img_struct->get_descriptor(i);
        mp.m_counter = 0;
        mp.m_age = 0;
        out_points->push_back(mp);
    }
}

void lvt_local_map::triangulate(const lvt_pose &cam_pose, lvt_image_features_struct *left_struct,
                                lvt_image_features_struct *right_struct, lvt_map_point_array *out_points)
{
    std::vector<cv::DMatch> matches;
    m_features_handler->row_match(left_struct, right_struct, &matches);
    if (matches.empty())
    {
        LVT_ASSERT(0 && "No row matches for triangulation found.");
        return;
    }

    const lvt_pose cam_pose_right = lvt_pose_utils::compute_right_camera_pose(cam_pose, m_vo_params.baseline);
    const lvt_matrix34 cml = lvt_pose_utils::compute_world_to_camera_transform(cam_pose);
    const lvt_matrix34 cmr = lvt_pose_utils::compute_world_to_camera_transform(cam_pose_right);
    const double cx = m_vo_params.cx, cy = m_vo_params.cy; // , f_inv = 1.0 / m_vo_params.focal_length;
    const double inv_fx = 1.0 / m_vo_params.fx, inv_fy = 1.0 / m_vo_params.fy;
    out_points->reserve(matches.size());

    for (size_t i = 0, count = matches.size(); i < count; i++)
    {
        const cv::Point2f u1 = left_struct->get_keypoint(matches[i].queryIdx).pt;
        const cv::Point2f u2 = right_struct->get_keypoint(matches[i].trainIdx).pt;
        double u1_x = (u1.x - cx) * inv_fx;
        double u1_y = (u1.y - cy) * inv_fy;
        double u2_x = (u2.x - cx) * inv_fx;
        double u2_y = (u2.y - cy) * inv_fy;

        // Linear-LS triangulation
        Eigen::Matrix<double, 4, 4> A;
        A << u1_x * cml.row(2) - cml.row(0),
            u1_y * cml.row(2) - cml.row(1),
            u2_x * cmr.row(2) - cmr.row(0),
            u2_y * cmr.row(2) - cmr.row(1);

        lvt_vector3 world_pt = A.leftCols<3>().jacobiSvd(Eigen::ComputeFullU | Eigen::ComputeFullV).solve(-A.col(3));
        assert(std::isfinite(world_pt.x()) && std::isfinite(world_pt.y()) && std::isfinite(world_pt.z()));

        // check if the point is in viewable region by camera
        lvt_vector2 proj_pt_l, proj_pt_r;
        if (!is_point_visible(world_pt, cml, m_vo_params, proj_pt_l) || !is_point_visible(world_pt, cmr, m_vo_params, proj_pt_r))
        {
            continue;
        }

        // check reprojection error
        {
            double err_x = proj_pt_l.x() - u1.x;
            double err_y = proj_pt_l.y() - u1.y;
            if ((err_x * err_x + err_y * err_y) > LVT_REPROJECTION_TH2)
            {
                continue;
            }
        }

        {
            double err_x = proj_pt_r.x() - u2.x;
            double err_y = proj_pt_r.y() - u2.y;
            if ((err_x * err_x + err_y * err_y) > LVT_REPROJECTION_TH2)
            {
                continue;
            }
        }

        // create map point
        lvt_map_point mp;
        mp.m_position = world_pt;
        mp.m_descriptor = left_struct->get_descriptor(matches[i].queryIdx);
        mp.m_counter = 0;
        mp.m_age = 0;
        out_points->push_back(mp);
    }
}

void lvt_local_map::update_with_new_triangulation(const lvt_pose &cam_pose,
                                                  lvt_image_features_struct *left_struct, lvt_image_features_struct *right_struct, bool dont_stage)
{
    // Triangulate new map points from features that were not matched/tracked
    lvt_map_point_array new_triangulations;
    if (left_struct->is_depth_associated())
    {
        triangulate_rgbd(cam_pose, left_struct, &new_triangulations);
    }
    else
    {
        triangulate(cam_pose, left_struct, right_struct, &new_triangulations);
    }
    LVT_ASSERT(!new_triangulations.empty() && "Nothing was triangulated");
    if (dont_stage || m_vo_params.staged_threshold == 0 || get_map_size() < LVT_N_MAP_POINTS)
    {
        m_map_points.insert(m_map_points.end(), new_triangulations.begin(), new_triangulations.end());
    }
    else
    {
        m_staged_points.insert(m_staged_points.end(), new_triangulations.begin(), new_triangulations.end());
    }
}

void lvt_local_map::update_staged_map_points(const lvt_pose &cam_pose, lvt_image_features_struct *left_struct)
{
    const lvt_matrix34 cml = lvt_pose_utils::compute_world_to_camera_transform(cam_pose);
    int n_erased_points = 0;
    int n_upgraded_points = 0;
    std::set<lvt_map_point *> points_to_be_deleted;
    for (int i = 0, count = m_staged_points.size(); i < count; i++)
    {
        lvt_map_point *mp = &(m_staged_points[i]);
        const lvt_vector3 world_pt = mp->m_position;
        lvt_vector2 proj_pt_left;
        float d1, d2;
        int match_idx_left = -1;
        if (!is_point_visible(world_pt, cml, m_vo_params, proj_pt_left) ||
            (match_idx_left = left_struct->find_match_index(proj_pt_left, mp->m_descriptor, &d1, &d2)) == -1)
        {
            points_to_be_deleted.insert(mp);
            n_erased_points++;
            continue;
        }
        left_struct->mark_as_matched(match_idx_left, true);
        mp->m_counter += 1;
        if (mp->m_counter == m_vo_params.staged_threshold || get_map_size() < LVT_N_MAP_POINTS)
        {
            m_map_points.push_back(m_staged_points[i]);
            points_to_be_deleted.insert(mp);
            n_upgraded_points++;
        }
    }

    m_staged_points.erase(std::remove_if(m_staged_points.begin(), m_staged_points.end(),
                                         [&points_to_be_deleted](lvt_map_point &pt) { return points_to_be_deleted.count(&pt) != 0; }),
                          m_staged_points.end());

    LVT_LOG("Staged map points updated: " + std::to_string(n_erased_points) + " erased, " +
            std::to_string(n_upgraded_points) + " upgraded to map points, " + std::to_string(m_staged_points.size()) + " remain.");
}

void lvt_local_map::clean_untracked_points(lvt_image_features_struct *left_struct)
{
    const int th = m_vo_params.untracked_threshold;
    lvt_map_point_array cleaned_map_points;
    cleaned_map_points.reserve(m_map_points.size());
    for (int i = 0; i < m_map_points.size(); i++)
    {
        if (m_map_points[i].m_counter >= th)
        {
            if (m_map_points[i].m_match_idx >= 0)
            {
                left_struct->mark_as_matched(m_map_points[i].m_match_idx, false);
            }
        }
        else
        {
            cleaned_map_points.push_back(m_map_points[i]);
        }
    }
    cleaned_map_points.swap(m_map_points);
}

max_stars_repo_name: iscumd/lvt

max_issues_repo_name: iscumd/lvt

max_forks_repo_name: iscumd/lvt


text:
// Copyright (C) 2011 by the BEM++ Authors
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
// THE SOFTWARE.

#include "bempp/common/config_trilinos.hpp"

#ifdef WITH_TRILINOS

#include "../check_arrays_are_close.hpp"
#include "../type_template.hpp"
#include "../random_arrays.hpp"

#include "create_regular_grid.hpp"

#include "assembly/assembly_options.hpp"
#include "assembly/discrete_boundary_operator.hpp"
#include "assembly/discrete_boundary_operator_composition.hpp"
#include "assembly/discrete_sparse_boundary_operator.hpp"
#include "assembly/boundary_operator.hpp"
#include "assembly/context.hpp"
#include "assembly/identity_operator.hpp"
#include "assembly/numerical_quadrature_strategy.hpp"
#include "assembly/sparse_cholesky.hpp"

#include "bempp/common/config_ahmed.hpp"

#include "grid/grid.hpp"

#include "space/piecewise_constant_scalar_space.hpp"
#include "space/piecewise_linear_discontinuous_scalar_space.hpp"
#include "space/piecewise_linear_continuous_scalar_space.hpp"

#include "common/armadillo_fwd.hpp"
#include "common/boost_make_shared_fwd.hpp"

#include <boost/test/unit_test.hpp>
#include <boost/test/floating_point_comparison.hpp>
#include <boost/version.hpp>
#include <complex>

using namespace Bempp;

// Tests

BOOST_AUTO_TEST_SUITE(SparseCholesky)

BOOST_AUTO_TEST_CASE(sparse_cholesky_works_for_piecewise_linears)
{
    typedef double BFT;
    typedef double RT;

    int nElementsX = 1, nElementsY = 2;
    shared_ptr<Grid> grid = createRegularTriangularGrid(nElementsX, nElementsY);

    shared_ptr<Space<BFT> > space(
        new PiecewiseLinearDiscontinuousScalarSpace<BFT>(grid));

    AssemblyOptions assemblyOptions;
    assemblyOptions.setVerbosityLevel(VerbosityLevel::LOW);
    shared_ptr<NumericalQuadratureStrategy<BFT, RT> > quadStrategy(
        new NumericalQuadratureStrategy<BFT, RT>);
    shared_ptr<Context<BFT, RT> > context(
        new Context<BFT, RT>(quadStrategy, assemblyOptions));

    BoundaryOperator<BFT, RT> op = identityOperator<BFT, RT>(
        context, space, space, space);
    shared_ptr<const DiscreteBoundaryOperator<RT> > dop = op.weakForm();
    typedef DiscreteSparseBoundaryOperator<RT> SparseOp;
    shared_ptr<const SparseOp> sdop = SparseOp::castToSparse(dop);
    shared_ptr<const Epetra_CrsMatrix> mat = sdop->epetraMatrix();
    shared_ptr<Epetra_CrsMatrix> L = sparseCholesky(*mat);

    shared_ptr<SparseOp> opL = boost::make_shared<SparseOp>(L, NO_SYMMETRY);
    shared_ptr<SparseOp> opLT = boost::make_shared<SparseOp>(L, NO_SYMMETRY,
                                                             TRANSPOSE);
    typedef DiscreteBoundaryOperatorComposition<RT> OpComposition;
    shared_ptr<DiscreteBoundaryOperatorComposition<RT> > opLLT =
            boost::make_shared<OpComposition>(opL, opLT);

    arma::Mat<double> A = dop->asMatrix();
    arma::Mat<double> LLT = opLLT->asMatrix();

//    arma::Mat<double> matL = opL->asMatrix();
//    std::cout << "A:\n" << A;
//    std::cout << "L:\n" << matL;
//    std::cout << "arma chol:\n" << arma::chol(A).t();
//    std::cout << "L*L.t():\n" << (matL * matL.t());
//    std::cout << "LLT:\n" << LLT;

    BOOST_CHECK(check_arrays_are_close<double>(
                    A, LLT, 100. * std::numeric_limits<double>::epsilon()));
}

BOOST_AUTO_TEST_CASE(sparse_cholesky_works_for_piecewise_constants)
{
    typedef double BFT;
    typedef double RT;

    int nElementsX = 1, nElementsY = 2;
    shared_ptr<Grid> grid = createRegularTriangularGrid(nElementsX, nElementsY);

    shared_ptr<Space<BFT> > space(
        new PiecewiseConstantScalarSpace<BFT>(grid));

    AssemblyOptions assemblyOptions;
    assemblyOptions.setVerbosityLevel(VerbosityLevel::LOW);
    shared_ptr<NumericalQuadratureStrategy<BFT, RT> > quadStrategy(
        new NumericalQuadratureStrategy<BFT, RT>);
    shared_ptr<Context<BFT, RT> > context(
        new Context<BFT, RT>(quadStrategy, assemblyOptions));

    BoundaryOperator<BFT, RT> op = identityOperator<BFT, RT>(
        context, space, space, space);
    shared_ptr<const DiscreteBoundaryOperator<RT> > dop = op.weakForm();
    typedef DiscreteSparseBoundaryOperator<RT> SparseOp;
    shared_ptr<const SparseOp> sdop = SparseOp::castToSparse(dop);
    shared_ptr<const Epetra_CrsMatrix> mat = sdop->epetraMatrix();
    shared_ptr<Epetra_CrsMatrix> L = sparseCholesky(*mat);

    shared_ptr<SparseOp> opL = boost::make_shared<SparseOp>(L, NO_SYMMETRY);
    shared_ptr<SparseOp> opLT = boost::make_shared<SparseOp>(L, NO_SYMMETRY,
                                                             TRANSPOSE);
    typedef DiscreteBoundaryOperatorComposition<RT> OpComposition;
    shared_ptr<DiscreteBoundaryOperatorComposition<RT> > opLLT =
            boost::make_shared<OpComposition>(opL, opLT);

    arma::Mat<double> A = dop->asMatrix();
    arma::Mat<double> LLT = opLLT->asMatrix();

    BOOST_CHECK(check_arrays_are_close<double>(
                    A, LLT, 100. * std::numeric_limits<double>::epsilon()));
}

BOOST_AUTO_TEST_CASE(sparse_cholesky_throws_for_continuous_functions)
{
    typedef double BFT;
    typedef double RT;

    int nElementsX = 1, nElementsY = 2;
    shared_ptr<Grid> grid = createRegularTriangularGrid(nElementsX, nElementsY);

    shared_ptr<Space<BFT> > space(
        new PiecewiseLinearContinuousScalarSpace<BFT>(grid));

    AssemblyOptions assemblyOptions;
    assemblyOptions.setVerbosityLevel(VerbosityLevel::LOW);
    shared_ptr<NumericalQuadratureStrategy<BFT, RT> > quadStrategy(
        new NumericalQuadratureStrategy<BFT, RT>);
    shared_ptr<Context<BFT, RT> > context(
        new Context<BFT, RT>(quadStrategy, assemblyOptions));

    BoundaryOperator<BFT, RT> op = identityOperator<BFT, RT>(
        context, space, space, space);
    shared_ptr<const DiscreteBoundaryOperator<RT> > dop = op.weakForm();
    typedef DiscreteSparseBoundaryOperator<RT> SparseOp;
    shared_ptr<const SparseOp> sdop = SparseOp::castToSparse(dop);
    shared_ptr<const Epetra_CrsMatrix> mat = sdop->epetraMatrix();
    BOOST_CHECK_THROW(sparseCholesky(*mat), std::invalid_argument   );
}

BOOST_AUTO_TEST_SUITE_END()

#endif // WITH_TRILINOS

max_stars_repo_name: mdavezac/bempp

max_issues_repo_name: mdavezac/bempp

max_forks_repo_name: mdavezac/bempp


text:
#include <cctbx/math/cos_sin_table.h>
#include <iostream>
#include <boost/timer.hpp>

int main() {
  unsigned const n = 1024, p = 1024;
  cctbx::math::cos_sin_table<double> table(n);
  cctbx::math::cos_sin_exact<double> exact;
  std::complex<double> sum = 0;
  boost::timer t;
  for (unsigned i=0; i<p*n; ++i) {
    sum += table(i);
  }
  std::cout << t.elapsed() << "\t" << sum << std::endl;
  sum = 0;
  t.restart();
  for (unsigned i=0; i<p*n; ++i) {
    sum += exact(i);
  }
  std::cout << t.elapsed() << "\t" << sum << std::endl;
  std::cout << "OK\n" << std::endl;
  return 0;
}


text:
//
// Copyright (c) 2019 Vinnie Falco (vinnie.falco@gmail.com)
//
// 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)
//
// Official repository: https://github.com/CPPAlliance/url
//

#ifndef BOOST_URL_IPV6_ADDRESS_HPP
#define BOOST_URL_IPV6_ADDRESS_HPP

#include <boost/url/detail/config.hpp>
#include <boost/url/error.hpp>
#include <boost/url/const_string.hpp>
#include <boost/url/grammar/parse_tag.hpp>
#include <array>
#include <cstdint>
#include <iosfwd>

namespace boost {
namespace urls {

#ifndef BOOST_URL_DOCS
class ipv4_address;
#endif

/** An IP version 6 style address.

    Objects of this type are used to construct,
    parse, and manipulate IP version 6 addresses.

    @par BNF
    @code
    IPv6address =                            6( h16 ":" ) ls32
                /                       "::" 5( h16 ":" ) ls32
                / [               h16 ] "::" 4( h16 ":" ) ls32
                / [ *1( h16 ":" ) h16 ] "::" 3( h16 ":" ) ls32
                / [ *2( h16 ":" ) h16 ] "::" 2( h16 ":" ) ls32
                / [ *3( h16 ":" ) h16 ] "::"    h16 ":"   ls32
                / [ *4( h16 ":" ) h16 ] "::"              ls32
                / [ *5( h16 ":" ) h16 ] "::"              h16
                / [ *6( h16 ":" ) h16 ] "::"

    ls32        = ( h16 ":" h16 ) / IPv4address
                ; least-significant 32 bits of address

    h16         = 1*4HEXDIG
                ; 16 bits of address represented in hexadecimal
    @endcode

    @par Specification
    @li IP Version 6 Addressing Architecture (rfc4291)
    @li 3.2.2. Host (rfc3986)

    @see
        @ref ipv4_address,
        @ref parse_ipv6_address.
*/
class ipv6_address
{
public:
    /** The number of characters in the longest possible IPv6 string.

        The longest IPv6 address is:
        @code
        ffff:ffff:ffff:ffff:ffff:ffff:ffff:ffff
        @endcode

        @see
            @ref to_buffer.
    */
    // ffff:ffff:ffff:ffff:ffff:ffff:ffff:ffff
    // ::ffff:255.255.255.255
    // 12345678901234567890123456789012345678901234567890
    //          1         2         3        4
    static
    constexpr
    std::size_t max_str_len = 49;

    /** The type used to represent an address as an array of bytes.

        Octets are stored in network byte order.
    */
    using bytes_type = std::array<
        unsigned char, 16>;

    /** Constructor.

        Default constructed objects represent
        the unspecified address.

        @li 2.5.2. The Unspecified Address

        @see
            @ref is_unspecified
    */
    ipv6_address() = default;

    /** Constructor.
    */
    ipv6_address(
        ipv6_address const&) = default;

    /** Copy Assignment
    */
    ipv6_address&
    operator=(
        ipv6_address const&) = default;

    /** Construct from an array of bytes.

        This function constructs an address
        from the array in `bytes`, which will
        be interpreted in big-endian.

        @param bytes The value to construct from.
    */
    BOOST_URL_DECL
    ipv6_address(
        bytes_type const& bytes) noexcept;

    /** Construct from an IPv4 address.

        This function constructs an IPv6 address
        from the IPv4 address `addr`. The resulting
        address is an IPv4-Mapped IPv6 Address.

        @param addr The address to construct from.

        @par Specification
        @li 2.5.5.2. IPv4-Mapped IPv6 Address (rfc4291)
    */
    BOOST_URL_DECL
    ipv6_address(
        ipv4_address const& addr) noexcept;

    /** Construct from a string.

        This function constructs an address from
        the string `s`, which must contain a valid
        IPv6 address string or else an exception
        is thrown.

        @note For a non-throwing parse function,
        use @ref parse_ipv6_address.

        @throw std::invalid_argument parse error.

        @param s The string to parse.

        @par Specification
        @li 3.2.2. Host (rfc3986)

        @see
            @ref parse_ipv6_address.
    */
    BOOST_URL_DECL
    ipv6_address(
        string_view s);

    /** Return the address as bytes, in network byte order
    */
    bytes_type
    to_bytes() const noexcept
    {
        return addr_;
    }

    /** Return the address as a string.

        The returned string will not
        contain surrounding square brackets.

        @par Exception Safety
        Strong guarantee.
        Calls to allocate may throw.

        @param a An optional allocator the returned
        string will use. If this parameter is omitted,
        the default allocator is used.

        @return A @ref const_string using the
        specified allocator.

        @par Specification
        @li
            2.2. Text Representation of Addresses (rfc4291)
    */
    template<class Allocator =
        std::allocator<char>>
    const_string
    to_string(Allocator const& a = {}) const
    {
        char buf[max_str_len];
        auto const n = print_impl(buf);
        return const_string(
            string_view(buf, n), a);
    }

    /** Write a dotted decimal string representing the address to a buffer

        The resulting buffer is not null-terminated.

        @throw std::length_error `dest_size < ipv6_address::max_str_len`

        @return The formatted string

        @param dest The buffer in which to write,
        which must have at least `dest_size` space.

        @param dest_size The size of the output buffer.
    */
    BOOST_URL_DECL
    string_view
    to_buffer(
        char* dest,
        std::size_t dest_size) const;

    /** Return true if the address is unspecified

        The address 0:0:0:0:0:0:0:0 is called the
        unspecified address. It indicates the
        absence of an address.

        @par Specification
        @li
            2.5.2. The Unspecified Address (rfc4291)
    */
    BOOST_URL_DECL
    bool
    is_unspecified() const noexcept;

    /** Return true if the address is a loopback address

        The unicast address 0:0:0:0:0:0:0:1 is called
        the loopback address. It may be used by a node
        to send an IPv6 packet to itself.

        @par Specification
        @li
            2.5.3. The Loopback Address (rfc4291)
    */
    BOOST_URL_DECL
    bool
    is_loopback() const noexcept;

    /** Return true if the address is a mapped IPv4 address

        This address type is used to represent the
        addresses of IPv4 nodes as IPv6 addresses.

        @par Specification
        @li
            2.5.5.2. IPv4-Mapped IPv6 Address (rfc4291)
    */
    BOOST_URL_DECL
    bool
    is_v4_mapped() const noexcept;

    /** Return true if two addresses are equal
    */
    BOOST_URL_DECL
    friend
    bool
    operator==(
        ipv6_address const& a1,
        ipv6_address const& a2) noexcept;

    /** Return true if two addresses are not equal
    */
    friend
    bool
    operator!=(
        ipv6_address const& a1,
        ipv6_address const& a2) noexcept
    {
        return !( a1 == a2 );
    }

    /** Return an address object that represents the loopback address

        The unicast address 0:0:0:0:0:0:0:1 is called
        the loopback address. It may be used by a node
        to send an IPv6 packet to itself.

        @par Specification
        @li
            2.5.3. The Loopback Address (rfc4291)
    */
    BOOST_URL_DECL
    static
    ipv6_address
    loopback() noexcept;

    /** Parse a string containing an IPv6 address.

        This function overload is used with
        @ref grammar::parse.

        @return true if the parse was successful.

        @param it A pointer to the beginning of
        the string to be parsed. Upon return, this
        will point to one past the last character
        visited.

        @param end A pointer to one past the last
        character in the string to be parsed.

        @param ec Set to the error, if any occurred.

        @param t Set to the result upon success.
    */
    BOOST_URL_DECL
    friend
    void
    tag_invoke(
        grammar::parse_tag const&,
        char const*& it,
        char const* const end,
        error_code& ec,
        ipv6_address& t) noexcept;

private:
    BOOST_URL_DECL
    std::size_t
    print_impl(
        char* dest) const noexcept;

    bytes_type addr_{};
};

//------------------------------------------------

/** Format the address to an output stream

    This function writes the address to an
    output stream using standard notation.

    @return The output stream, for chaining.

    @param os The output stream to write to.

    @param addr The address to write.
*/
BOOST_URL_DECL
std::ostream&
operator<<(
    std::ostream& os,
    ipv6_address const& addr);

/** Parse a string containing an IPv6 address.

    This function attempts to parse the string
    as an IPv6 address and returns a result
    containing the address upon success, or
    an error code if the string does not contain
    a valid IPv6 address.

    @par Exception Safety
    Throws nothing.

    @return A result containing the address.

    @param s The string to parse.
*/
BOOST_URL_DECL
result<ipv6_address>
parse_ipv6_address(
    string_view s) noexcept;

} // urls
} // boost

#endif

max_stars_repo_name: vgrigoriu/url

max_issues_repo_name: vgrigoriu/url

max_forks_repo_name: vgrigoriu/url


text:
#define BOOST_TEST_MODULE TestBackupCatalog
#include <boost/test/unit_test.hpp>
#include <common.hxx>
#include <BackupCatalog.hxx>

using namespace pgbckctl;

BOOST_AUTO_TEST_CASE(TestBackupCatalogSetup)
{

  std::shared_ptr<BackupCatalog> catalog = nullptr;

  /* 1 should not throw */
  BOOST_REQUIRE_NO_THROW( BackupCatalog() );

  /* 2 Should throw */
  BOOST_CHECK_THROW( BackupCatalog ("/tmp/foobar"),
                     std::exception);

  /* 3 Create valid SQLite database handle */
  BOOST_REQUIRE_NO_THROW( catalog
                          = std::make_shared<BackupCatalog>(".pg_backup_ctl.sqlite") );

  /* 4 Open backup catalog for read/write */
  BOOST_REQUIRE_NO_THROW( catalog->open_rw() );

  /* 6 Close backup catalog database */
  BOOST_REQUIRE_NO_THROW( catalog->close() );

}

BOOST_AUTO_TEST_CASE(TestBackupCatalogCreateArchive)
{

  std::shared_ptr<BackupCatalog> catalog = nullptr;

  /* 1 should not throw */
  BOOST_REQUIRE_NO_THROW( catalog
                          = std::make_shared<BackupCatalog>(".pg_backup_ctl.sqlite") );

  /* 2 Open backup catalog for read/write */
  BOOST_REQUIRE_NO_THROW( catalog->open_rw() );

  /* 3 Backup catalog should be available */
  BOOST_CHECK( catalog->available() );

  /*
   * 4 Create archive "test"
   *
   * NOTE: BackupCatalog::createArchive doesn't require
   *                      pushAffectedAttributes() anymore,
   *                      since we always need the same number
   *                      of attributes here.
   */
  {
    std::shared_ptr<CatalogDescr> desc = std::make_shared<CatalogDescr>();
    std::shared_ptr<CatalogDescr> check_desc;

    BOOST_REQUIRE_NO_THROW( catalog->startTransaction() );

    /* CatalogDescr::archive_name required, should throw */
    BOOST_CHECK_THROW( catalog->createArchive(desc),
                       CCatalogIssue );

    desc->archive_name = "test";

    /* CatalogDescr::directory required, should throw */
    BOOST_CHECK_THROW( catalog->createArchive(desc),
                       CCatalogIssue );

    desc->directory = "/tmp";

    /*
     * Set connection type. this is important since BackupCatalog::createArchive()
     * insists this to be set! Thus, this should throw.
     */
    BOOST_CHECK_THROW( catalog->createArchive(desc),
                       CCatalogIssue );
    desc->coninfo->type = ConnectionDescr::CONNECTION_TYPE_BASEBACKUP;

    /* Should succeed now */
    /* Force compression off, though not required */
    desc->compression = false;
    BOOST_REQUIRE_NO_THROW( catalog->createArchive(desc) );

    /* Check whether archive exists */
    BOOST_REQUIRE_NO_THROW( check_desc = catalog->existsByName("test") );
    BOOST_CHECK( check_desc->id != -1 );
    BOOST_TEST( check_desc->id > -1 );
    BOOST_CHECK_EQUAL( check_desc->archive_name, "test" );
    BOOST_CHECK_EQUAL( check_desc->directory, "/tmp" );
    BOOST_CHECK_EQUAL( check_desc->compression, false );

    /* Create a basebackup streaming connection */

    /* This must fail, since no connection DSN/attributes currently set */
    BOOST_CHECK_THROW( catalog->createCatalogConnection(check_desc->coninfo),
                       CCatalogIssue );

    /* Assign a DSN string and attributes, this should succeed now */
    check_desc->coninfo->pushAffectedAttribute(SQL_CON_DSN_ATTNO);
    check_desc->coninfo->pushAffectedAttribute(SQL_CON_ARCHIVE_ID_ATTNO);
    check_desc->coninfo->pushAffectedAttribute(SQL_CON_TYPE_ATTNO);

    check_desc->coninfo->dsn = "host=bar.server.name dbname=foo user=test";
    check_desc->coninfo->archive_id = check_desc->id;
    check_desc->coninfo->type = ConnectionDescr::CONNECTION_TYPE_BASEBACKUP;

    BOOST_REQUIRE_NO_THROW( catalog->createCatalogConnection(check_desc->coninfo) );

    /* Drop archive and recheck */
    BOOST_REQUIRE_NO_THROW( catalog->dropArchive("test") );

    BOOST_REQUIRE_NO_THROW( check_desc = catalog->existsByName("test") );
    BOOST_CHECK( check_desc->id == -1 );

    BOOST_REQUIRE_NO_THROW( catalog->commitTransaction() );
  }

  /* 5 Close works */
  BOOST_REQUIRE_NO_THROW( catalog->close() );

  /* 5 Catalog not available anymore */
  BOOST_TEST( !catalog->available() );

}

ext: cxx


text:
//==============================================================================
//         Copyright 2003 - 2012   LASMEA UMR 6602 CNRS/Univ. Clermont II
//         Copyright 2009 - 2012   LRI    UMR 8623 CNRS/Univ Paris Sud XI
//
//          Distributed under the Boost Software License, Version 1.0.
//                 See accompanying file LICENSE.txt or copy at
//                     http://www.boost.org/LICENSE_1_0.txt
//==============================================================================
#ifndef BOOST_SIMD_SDK_SIMD_EXTENSIONS_X86_SSE2_HPP_INCLUDED
#define BOOST_SIMD_SDK_SIMD_EXTENSIONS_X86_SSE2_HPP_INCLUDED

#if defined(__SSE2__) || defined(_M_X64) || _M_IX86_FP >= 2
#  ifndef BOOST_SIMD_HAS_SSE2_SUPPORT
#    define BOOST_SIMD_HAS_SSE2_SUPPORT
#  endif
#elif defined(BOOST_SIMD_HAS_SSE2_SUPPORT) && !defined(_MSC_VER)
#  undef BOOST_SIMD_HAS_SSE2_SUPPORT
#endif

#if defined(BOOST_SIMD_HAS_SSE2_SUPPORT) && !defined(BOOST_SIMD_HAS_SSE_SUPPORT)
#  define BOOST_SIMD_HAS_SSE_SUPPORT
#endif

#if !defined(BOOST_SIMD_DETECTED) && defined(BOOST_SIMD_HAS_SSE2_SUPPORT)

////////////////////////////////////////////////////////////////////////////////
// SSE2 extensions flags
////////////////////////////////////////////////////////////////////////////////
#define BOOST_SIMD_DETECTED
#define BOOST_SIMD_SSE2
#define BOOST_SIMD_SSE_FAMILY
#define BOOST_SIMD_STRING             "SSE2"
#define BOOST_SIMD_STRING_LIST        "SSE2"
#define BOOST_SIMD_BYTES              16
#define BOOST_SIMD_BITS               128
#define BOOST_SIMD_CARDINALS          (2)(4)(8)(16)
#define BOOST_SIMD_TAG_SEQ            (::boost::simd::tag::sse_)
#define BOOST_SIMD_DEFAULT_EXTENSION  ::boost::simd::tag::sse_
#define BOOST_SIMD_DEFAULT_SITE       ::boost::simd::tag::sse2_

#include <emmintrin.h>

#include <boost/simd/sdk/simd/extensions/meta/sse.hpp>

#endif
#endif


text:
#include "util.h"

#include <stdint.h>
#include <vector>

#include "chainparams.h"
#include "key.h"
#include "validation.h"
#include "pubkey.h"
#include "txdb.h"
#include "txmempool.h"
#include "sigma.h"

#include "test/fixtures.h"
#include "test/testutil.h"

#include "wallet/db.h"
#include "wallet/wallet.h"
#include "wallet/walletexcept.h"

#include <boost/filesystem.hpp>
#include <boost/test/unit_test.hpp>
#include <boost/thread.hpp>

using namespace std;

BOOST_FIXTURE_TEST_SUITE(sigma_mintspend, ZerocoinTestingSetup200)

/*
* 1. Create one denomination pair and check it can't be spend till 2 conf of mint
* 2. Make one more mint of denom pair and check it can't be spend till 2 conf
* 3. Create two spend transactions using same mint
* 4. Double spend with previous spend in last block
*/
BOOST_AUTO_TEST_CASE(sigma_mintspend_test)
{
    sigma::CSigmaState *sigmaState = sigma::CSigmaState::GetState();
    string denomination;
    std::vector<string> denominations = {"0.05", "0.1", "0.5", "1", "10", "25", "100"};

    // foreach denom from denominations
    for(auto denomination : denominations)
    {
        string stringError;
        // Make sure that transactions get to mempool
        pwalletMain->SetBroadcastTransactions(true);

        sigma::CoinDenomination denom;
        BOOST_CHECK_MESSAGE(StringToDenomination(denomination, denom), "Unable to convert denomination string to value.");

        // Verify Mint is successful
        std::vector<sigma::PrivateCoin> privCoins;
        const auto& sigmaParams = sigma::Params::get_default();
        privCoins.push_back(sigma::PrivateCoin(sigmaParams, denom));
        vector<CHDMint> vDMints;
        auto vecSend = CWallet::CreateSigmaMintRecipients(privCoins, vDMints);
        {
            CWalletTx wtx;
            stringError = pwalletMain->MintAndStoreSigma(vecSend, privCoins, vDMints, wtx);
        }
        BOOST_CHECK_MESSAGE(stringError == "", "Create Mint Failed");

        // Verify Mint gets in the mempool
        BOOST_CHECK_MESSAGE(mempool.size() == 1, "Mint was not added to mempool");

        int previousHeight = chainActive.Height();
        CBlock b = CreateAndProcessBlock(scriptPubKey);
        BOOST_CHECK_MESSAGE(previousHeight + 1 == chainActive.Height(), "Block not added to chain");

        previousHeight = chainActive.Height();

        // Generate address
        CPubKey newKey;
        BOOST_CHECK_MESSAGE(pwalletMain->GetKeyFromPool(newKey), "Fail to get new address");

        const CBitcoinAddress randomAddr(newKey.GetID());

        CAmount nValue;
        DenominationToInteger(denom, nValue);

        std::vector<CRecipient> recipients = {
                {GetScriptForDestination(randomAddr.Get()), nValue, true},
        };

        // Add 1 more blocks and verify that Mint can not be spent until 2 blocks verification
        {
            CWalletTx wtx;
            BOOST_CHECK_THROW(pwalletMain->SpendSigma(recipients, wtx), WalletError); //this must throw as 2 blocks have not passed yet,
        }

        b = CreateAndProcessBlock(scriptPubKey);

        BOOST_CHECK_MESSAGE(previousHeight + 1 == chainActive.Height(), "Block not added to chain");

        {
            CWalletTx wtx;
            BOOST_CHECK_THROW(pwalletMain->SpendSigma(recipients, wtx), WalletError); //this must throw as it has to have at least two mint coins with at least 2 confirmation
        }

        privCoins.clear();
        privCoins.push_back(sigma::PrivateCoin(sigmaParams, denom));
        vDMints.clear();
        vecSend = CWallet::CreateSigmaMintRecipients(privCoins, vDMints);
        {
            CWalletTx wtx;
            stringError = pwalletMain->MintAndStoreSigma(vecSend, privCoins, vDMints, wtx);
        }
        BOOST_CHECK_MESSAGE(stringError == "", "Create Mint Failed");

        BOOST_CHECK_MESSAGE(mempool.size() == 1, "Mint was not added to mempool");

        previousHeight = chainActive.Height();
        b = CreateAndProcessBlock(scriptPubKey);
        BOOST_CHECK_MESSAGE(previousHeight + 1 == chainActive.Height(), "Block not added to chain");

        previousHeight = chainActive.Height();
        //Add 1 more blocks and verify that Mint can not be spent until 2 blocks verification
        {
            CWalletTx wtx;
            BOOST_CHECK_THROW(pwalletMain->SpendSigma(recipients, wtx), WalletError); //this must throw as 2 blocks have not passed yet,
        }
        b = CreateAndProcessBlock(scriptPubKey);

        BOOST_CHECK_MESSAGE(previousHeight + 1 == chainActive.Height(), "Block not added to chain");

        std::vector<CSigmaEntry> coins;
        CWalletTx wtx;
        {
            std::vector <CHDMint> changes;
            CAmount fee;
            bool fChangeAddedToFee;
            wtx = pwalletMain->CreateSigmaSpendTransaction(recipients, fee, coins, changes, fChangeAddedToFee);
            pwalletMain->CommitSigmaTransaction(wtx, coins, changes);
        }
        //Set mints unused, and try to spend again
        for(auto mint : coins)
            pwalletMain->zwallet->GetTracker().SetPubcoinNotUsed(primitives::GetPubCoinValueHash(mint.value));

        wtx.Init(NULL);
        //try double spend
        BOOST_CHECK_NO_THROW(pwalletMain->SpendSigma(recipients, wtx));
        //Verify spend got into mempool
        BOOST_CHECK_MESSAGE(mempool.size() == 1, "Spend was not added to mempool");

        b = CreateBlock(scriptPubKey);
        previousHeight = chainActive.Height();
        BOOST_CHECK_MESSAGE(ProcessBlock(b), "ProcessBlock failed although valid spend inside");
        BOOST_CHECK_MESSAGE(previousHeight + 1 == chainActive.Height(), "Block not added to chain");

        BOOST_CHECK_MESSAGE(mempool.size() == 0, "Mempool not cleared");

        //roll back mints used
        for(auto mint : coins)
            pwalletMain->zwallet->GetTracker().SetPubcoinUsed(primitives::GetPubCoinValueHash(mint.value), uint256());

        coins.clear();
        wtx.Init(NULL);
        {
            std::vector <CHDMint> changes;
            CAmount fee;
            bool fChangeAddedToFee;
            wtx = pwalletMain->CreateSigmaSpendTransaction(recipients, fee, coins, changes, fChangeAddedToFee);
            pwalletMain->CommitSigmaTransaction(wtx, coins, changes);
        }

        //Verify spend got into mempool
        BOOST_CHECK_MESSAGE(mempool.size() == 1, "Spend was not added to mempool");

        b = CreateBlock(scriptPubKey);
        previousHeight = chainActive.Height();
        BOOST_CHECK_MESSAGE(ProcessBlock(b), "ProcessBlock failed although valid spend inside");
        BOOST_CHECK_MESSAGE(previousHeight + 1 == chainActive.Height(), "Block not added to chain");

        BOOST_CHECK_MESSAGE(mempool.size() == 0, "Mempool not cleared");

        //Set mints unused, and try to spend again
        for(auto mint : coins)
            pwalletMain->zwallet->GetTracker().SetPubcoinNotUsed(primitives::GetPubCoinValueHash(mint.value));

        //try double spend
        BOOST_CHECK_NO_THROW(pwalletMain->SpendSigma(recipients, wtx));

        BOOST_CHECK_MESSAGE(mempool.size() == 0, "Mempool not empty although mempool should reject double spend");

        //Temporary disable usedCoinSerials check to force double spend in mempool
        auto tempSerials = sigmaState->containers.usedCoinSerials;
        sigmaState->containers.usedCoinSerials.clear();

        {
            //Set mints unused, and try to spend again
            for(auto mint : coins)
                pwalletMain->zwallet->GetTracker().SetPubcoinNotUsed(primitives::GetPubCoinValueHash(mint.value));
            CWalletTx wtx;
            pwalletMain->SpendSigma(recipients, wtx);
            BOOST_CHECK_MESSAGE(mempool.size() == 1, "Spend was not added to mempool");
        }

        sigmaState->containers.usedCoinSerials = tempSerials;

        BOOST_CHECK_EXCEPTION(CreateBlock(scriptPubKey), std::runtime_error, no_check);
        BOOST_CHECK_MESSAGE(mempool.size() == 1, "Mempool not set");
        tempSerials = sigmaState->containers.usedCoinSerials;
        sigmaState->containers.usedCoinSerials.clear();
        CreateBlock(scriptPubKey);

        sigmaState->containers.usedCoinSerials = tempSerials;

        mempool.clear();
        previousHeight = chainActive.Height();
        BOOST_CHECK_MESSAGE(ProcessBlock(b), "ProcessBlock failed");
        //This test confirms that a block containing a double spend is rejected and not added in the chain
        BOOST_CHECK_MESSAGE(previousHeight == chainActive.Height(), "Double spend - Block added to chain even though same spend in previous block");

        mempool.clear();
        sigmaState->Reset();
    }
}
BOOST_AUTO_TEST_SUITE_END()

max_stars_repo_name: KuroGuo/firo

max_issues_repo_name: KuroGuo/firo

max_forks_repo_name: KuroGuo/firo


text:
#include <iostream>

extern "C"
{
	#include "lua.h"
	#include "lauxlib.h"
	#include "lualib.h"
}

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

#include "directory_iterator.hpp"

const boost::filesystem::path&
identity(const boost::filesystem::path& x)
{
	return x;
}

void bind_filesystem(lua_State* L)
{
	using namespace luabind;
	using namespace boost::filesystem;

	namespace fs = boost::filesystem;

	module(L, "filesystem")
	[
		class_<fs::path>("path")
			.def(constructor<>())
			.def(constructor<const char*>())
			.def("string", &fs::path::string)
			.def("native_file_string", &fs::path::native_file_string)
			.def("native_directory_string", &fs::path::native_directory_string)
			.def("root_path", &fs::path::root_path)
			.def("root_name", &fs::path::root_name)
			.def("root_directory", &fs::path::root_directory)
			.def("relative_path", &fs::path::relative_path)
			.def("leaf", &fs::path::leaf)
			.def("branch_path", &fs::path::branch_path)

			.def("empty", &fs::path::empty)
			.def("is_complete", &fs::path::is_complete)
			.def("is_directory", &fs::is_directory)
			.def("is_empty", &fs::is_empty)
			.def("has_root_path", &fs::path::has_root_path)
			.def("has_root_name", &fs::path::has_root_name)
			.def("has_root_directory", &fs::path::has_root_directory)
			.def("has_relative_path", &fs::path::has_relative_path)
			.def("has_leaf", &fs::path::has_leaf)
			.def("has_branch_path", &fs::path::has_branch_path)

			.def(const_self / const_self)
			.def(other<const char*>() / const_self)
			.def(const_self / other<const char*>())

//			.property("contents", &identity, return_directory_iterator)
			,

		def("exists", &fs::exists),
		def("is_directory", &fs::is_directory),
		def("is_empty", &fs::is_empty),
		def("create_directory", &fs::create_directory),
		def("remove", &fs::remove),
		def("remove_all", &fs::remove_all),
		def("rename", &fs::rename),
		def("copy_file", &fs::copy_file),

		def("initial_path", &fs::initial_path),
		def("current_path", &fs::current_path),
		def("complete", &fs::complete),
		def("system_complete", &fs::system_complete)
	];
}

int main(int argc, const char* argv[])
{
	lua_State* L = luaL_newstate();
	luaopen_base(L);
	luaopen_string(L);
	luaopen_table(L);
	luaopen_math(L);
	luaopen_io(L);
	luaopen_debug(L);

	if (argc < 2)
	{
		std::cout << "This example will bind parts of the boost.Filesystem library\n"
			"and then run the given lua file.\n\n"
			"usage: filesystem filename [args]\n";
		return 1;
	}

	using namespace luabind;

	open(L);

	bind_filesystem(L);

	object args = newtable(L);

	for (int i = 0; i < argc; ++i)
	{
		args[i + 1] = argv[i];
	}

	args["n"] = argc;
	globals(L)["args"] = args;

	luaL_dofile(L, argv[1]);
}