BOOST_FORCEINLINE
 void del(A0& a0, A1 & a1, const char& orient, boost::mpl::long_<2> const&) const
 {
   BOOST_AUTO_TPL(r1, boost::proto::child_c<1>(a1));
   BOOST_AUTO_TPL(q1, boost::proto::child_c<0>(a1));
   do_it(a0,orient,q1,r1);
 }

    static inline void apply(Polygon const& polygon_in, ring_identifier id,
                Polygon& polygon_out, MarkMap const& mark_map)
    {
        typedef typename geometry::ring_type<Polygon>::type ring_type;

        typedef range_remove_marked<ring_type, MarkMap> per_range;
        id.ring_index = -1;
        per_range::apply(exterior_ring(polygon_in), id, exterior_ring(polygon_out), mark_map);


        typename interior_return_type<Polygon const>::type rings_in
                    = interior_rings(polygon_in);
        typename interior_return_type<Polygon>::type rings_out
                    = interior_rings(polygon_out);

        rings_out.resize(boost::size(interior_rings(polygon_in)));
        BOOST_AUTO_TPL(out, boost::begin(rings_out));

        for (BOOST_AUTO_TPL(it, boost::begin(rings_in));
            it != boost::end(rings_in);
            ++it, ++out)
        {
            id.ring_index++;
            per_range::apply(*it, id, *out, mark_map);
        }
    }

    static inline void apply(Polygon1 const& source, Polygon2& destination)
    {
        // Clearing managed per ring, and in the resizing of interior rings

        per_ring::apply(geometry::exterior_ring(source),
                        geometry::exterior_ring(destination));

        // Container should be resizeable
        traits::resize
        <
        typename pdalboost::remove_reference
        <
        typename traits::interior_mutable_type<Polygon2>::type
        >::type
        >::apply(interior_rings(destination), num_interior_rings(source));

        typename interior_return_type<Polygon1 const>::type rings_source
            = interior_rings(source);
        typename interior_return_type<Polygon2>::type rings_dest
            = interior_rings(destination);

        BOOST_AUTO_TPL(it_source, pdalboost::begin(rings_source));
        BOOST_AUTO_TPL(it_dest, pdalboost::begin(rings_dest));

        for ( ; it_source != pdalboost::end(rings_source); ++it_source, ++it_dest)
        {
            per_ring::apply(*it_source, *it_dest);
        }
    }

 static void evenx2(A0& out, std::size_t n)
 {
   value_type nn = n;
   BOOST_AUTO_TPL(p, nt2::_(value_type(1), nn));
   BOOST_AUTO_TPL(j, nt2::trunc(nt2::mod(p,nt2::Four<value_type>())*nt2::Half<value_type>()));
   BOOST_AUTO_TPL(k, nt2::sx(nt2::tag::is_equal_(), j, nt2::colvect(j)));
   BOOST_AUTO_TPL(p1, colvect(nt2::_(value_type(1), nn, nt2::sqr(nn))));
   out =  nt2::sx(nt2::tag::plus_(),p1,nt2::minusone(p));
   out(k) = nt2::oneplus(sqr(nn)) - out(k);
 }

 static void compute(A0& out,const C& c, const R& r)
 {
   size_t p = numel(r);
   size_t m = numel(c);
   BOOST_AUTO_TPL(idx, nt2::_(ptrdiff_t(p), ptrdiff_t(-1), ptrdiff_t(2)));
   BOOST_AUTO_TPL(ridx, nt2::colvect(r(idx)));
   BOOST_AUTO_TPL(x, catv(ridx, c)); //build vector of user data
   BOOST_AUTO_TPL(v1, nt2::fliplr(nt2::cif(m, p, meta::as_<ptrdiff_t>())));
   BOOST_AUTO_TPL(v2, nt2::ric(m, p, meta::as_<ptrdiff_t>()));
   out(nt2::_) = x(v1+v2);
 }

    BOOST_FORCEINLINE
      static void matrix_compute(const size_t & n, const size_t & k, A1 & a1)
    {
      typedef typename boost::proto::result_of::child_c<A1&,0>::type       a_t;
      a_t & a = boost::proto::child_c<0>(a1);
      BOOST_AUTO_TPL(c, nt2::cumprod(nt2::_(value_t(2), value_t(n+1))));
      BOOST_AUTO_TPL(d, nt2::cumprod(nt2::_(value_t(n+1), value_t(2*n)))*c(n-1));

      a = nt2::hankel(c, d);
      if (k) a = nt2::rec(a);
    }

 static inline A0 finalize(const A0& a0, const A0& y)
 {
 #ifdef BOOST_SIMD_NO_NANS
   BOOST_AUTO_TPL(test, nt2::is_ltz(a0));
 #else
   BOOST_AUTO_TPL(test, nt2::logical_or(nt2::is_ltz(a0), nt2::is_nan(a0)));
 #endif
   A0 y1 = nt2::if_nan_else(test, y);
 #ifndef BOOST_SIMD_NO_INFINITIES
   y1 = if_else(nt2::is_equal(a0, nt2::Inf<A0>()), a0, y1);
 #endif
   return if_else(is_eqz(a0), nt2::Minf<A0>(), y1);
 }

 value_t compute(const  FUNC &f, const X&ranges, const input_t & vol,  input_t fact)
 {
     size_t n = size(ranges, 1);
     BOOST_AUTO_TPL(r1, ranges(nt2::_, begin_));
     BOOST_AUTO_TPL(r1ex, nt2::expand_to(r1,   nt2::of_size(n, nbpts_)));
     BOOST_AUTO_TPL(r2ex, nt2::expand_to(fact, nt2::of_size(n, nbpts_)));
     rtab_t rnd = nt2::quasi(n, nbpts_, nt2::meta::as_<real_t>());
     BOOST_AUTO_TPL(x, nt2::fma(rnd, r2ex, r1ex));
     BOOST_AUTO_TPL(z, f(x));
     if (compute_err_)
     {
         real_t fnbpt = real_t(nbpts_);
         err_ += nt2::pow(nt2::log(fnbpt), real_t(n))*nt2::Three<real_t>()*vol*nt2::real(nt2::globalstdev(z))/fnbpt;
     }
     return vol*nt2::globalmean(z);
 }

 static inline void apply(MultiGeometry& multi, Functor& f)
 {
     for(BOOST_AUTO_TPL(it, boost::begin(multi)); it != boost::end(multi); ++it)
     {
         Policy::apply(*it, f);
     }
 }

 BOOST_FORCEINLINE static void conf_bounds(const A0& a0, A1& a1,
                                           const value_type& alpha )
 {
   nt2::container::table<value_type> pcov =  boost::proto::child_c<3>(a0);
   const In0& x  = boost::proto::child_c<0>(a0);
   const In1& mu = boost::proto::child_c<1>(a0);
   const In2& sigma = boost::proto::child_c<2>(a0);
   BOOST_AUTO_TPL(z, (x-mu)/sigma);
   BOOST_AUTO_TPL(zvar, pcov(1,1) + (Two<value_type>()*pcov(1,2) + pcov(2,2)*z)*z);
   BOOST_ASSERT_MSG(nt2::globalall(nt2::is_gez(zvar)), "Covariance matrix must be positive");
   value_type normz = -nt2::norminv(alpha*nt2::Half<value_type>());
   BOOST_AUTO_TPL(halfwidth, normz*nt2::sqrt(zvar)/sigma);
   boost::proto::child_c<0>(a1) = Half<value_type>()*nt2::erfc(-Sqrt_2o_2<value_type>()*z);
   boost::proto::child_c<1>(a1) = Half<value_type>()*nt2::erfc(-Sqrt_2o_2<value_type>()*(z-halfwidth));
   boost::proto::child_c<2>(a1) = Half<value_type>()*nt2::erfc(-Sqrt_2o_2<value_type>()*(z+halfwidth));
 }

 void compute( const FUNC& f, const X & x, const o_t & o)
 {
   real_t tmptol = tol_;
   init(o, x);
   BOOST_AUTO_TPL(dif, nt2::abs(nt2::diff(nt2::rowvect(wpts_))));
   real_t tol1= tol_/nt2::globalasum1(dif);
   size_t l = numel(dif);
   res_.resize(extent(wpts_));
   res_(1) = nt2::Zero<value_t>();
   tol_ = tol1*dif(1);
   res_(2) = compute<true>(f, wpts_(1), wpts_(2));
   for(size_t i=2; i < l; ++i)
   {
     tol_ = tol1*dif(i);
     res_(i+1) = res_(i)+compute<false>(f, wpts_(i), wpts_(i+1));
   }
   if (l >= 2)
   {
     tol_ = tol1*dif(l);
     res_(l+1) = res_(l)+compute<true>(f, wpts_(l), wpts_(l+1));
   }
   tol_ = tmptol;
   if (!o.return_waypoints)
   {
     res_(begin_) = res_(end_);
     res_.resize(nt2::of_size(1, 1));
   }
 }

    static inline void apply(Range& range, Functor& f)
    {
        typedef typename add_const_if_c
            <
                is_const<Range>::value,
                typename point_type<Range>::type
            >::type point_type;

        BOOST_AUTO_TPL(it, boost::begin(range));
        BOOST_AUTO_TPL(previous, it++);
        while(it != boost::end(range))
        {
            model::referring_segment<point_type> s(*previous, *it);
            f(s);
            previous = it++;
        }
    }

    BOOST_FORCEINLINE typename Target::type
    operator()(Pos const& p, Size const& sz, Target const&) const
    {
      typedef typename Target::type                                         type;
      typedef typename meta::as_index<type>::type                            i_t;
      typedef typename meta::call<tag::enumerate_(Pos,meta::as_<i_t>)>::type p_t;
      typedef typename nt2::result_of::as_subscript<Size,p_t>::type          s_t;
      typedef typename s_t::value_type                                       v_t;

      // Retrieve 2D position from the linear index
      s_t const pos = as_subscript(sz,nt2::enumerate<i_t>(p));
      // Return a 1 where it belongs
      BOOST_AUTO_TPL(ii, nt2::splat<v_t>(i_));
      BOOST_AUTO_TPL(jj, nt2::splat<v_t>(j_));
      return nt2::if_else(nt2::logical_and(nt2::eq(pos[0],ii),
                                           nt2::eq(pos[1],jj)),
                          nt2::One<type>(), nt2::Zero<type>());
    }

 BOOST_FORCEINLINE static void conf_bounds(const A0& a0, A1& a1,
                                           const value_type& alpha )
 {
   typedef typename boost::proto::result_of::child_c<A0&,1>::type          In1;
   value_type pcov =  boost::proto::child_c<2>(a0);
   BOOST_ASSERT_MSG(pcov >= Zero<value_type>(), "Variance pcov must be non-negative");
   const In0& q  = boost::proto::child_c<0>(a0);
   const In1& mu = boost::proto::child_c<1>(a0);
   BOOST_ASSERT_MSG(nt2::globalall(nt2::is_gtz(mu)), "mu parameter must be positive");
   BOOST_AUTO_TPL(x, (-nt2::log1p(-q)*mu));
   //      nt2::container::table<value_type> logx = nt2::log(q*mu);
   value_type normz = -nt2::norminv(alpha*nt2::Half<value_type>());
   BOOST_AUTO_TPL(halfwidth, normz * nt2::sqrt(pcov/nt2::sqr(mu)));
   BOOST_AUTO_TPL(exp_halfwidth, exp(halfwidth));
   boost::proto::child_c<0>(a1) = x;
   boost::proto::child_c<2>(a1) = x*exp_halfwidth;
   boost::proto::child_c<1>(a1) = x/exp_halfwidth;
 }

 static inline ReturnType sum_interior_rings(Rings const& rings, Strategy const& strategy)
 {
     ReturnType sum = ReturnType();
     for (BOOST_AUTO_TPL(it, boost::begin(rings)); it != boost::end(rings); ++it)
     {
         sum += Policy::apply(*it, strategy);
     }
     return sum;
 }

 static inline Functor apply(
                 typename add_const_if_c<IsConst, MultiGeometry>::type& multi,
                 Functor f)
 {
     for(BOOST_AUTO_TPL(it, boost::begin(multi)); it != boost::end(multi); ++it)
     {
         f = Policy::apply(*it, f);
     }
     return f;
 }

 result_type operator()(A0& out, const A1& in) const
 {
   BOOST_AUTO_TPL(x,boost::proto::child_c<0>(in));
   size_t k = boost::proto::child_c<1>(in);
   size_t m = boost::proto::child_c<2>(in);
   size_t n = numel(x);
   BOOST_AUTO_TPL(a, nt2::expand(nt2::from_diag(x), nt2::of_size(m, n)));
   tab_t aa;
   if (!k)
   {
     tab_t q = nt2::qmult(n, nt2::meta::as_<value_t>());
     aa = nt2::mtimes(nt2::mtimes(q, a),nt2::transpose(q));//Not exploiting symmetry here.
   }
   else
   {
     aa = a;
   }
   do
   {
     tab_t diaga = nt2::diag_of(aa);
     nt2::table<ptrdiff_t> y = nt2::find(nt2::lt(diaga, nt2::One<value_t>()));
     nt2::table<ptrdiff_t> z = nt2::find(nt2::gt(diaga, nt2::One<value_t>()));
     if(nt2::isempty(y) || nt2::isempty(z)) break;
     size_t i = y(nt2::iceil(nt2::rand(nt2::meta::as_<value_t>())*nt2::numel(y)));
     size_t j = z(nt2::iceil(nt2::rand(nt2::meta::as_<value_t>())*nt2::numel(z)));
     if (i > j) std::swap(i, j);
     value_t aij = aa(i, j);
     value_t alpha = nt2::sqrt(nt2::sqr(aij)- nt2::minusone(diaga(i))*nt2::minusone(diaga(j)));
     value_t t1 = (aij + nt2::signnz(aij)*alpha)/nt2::minusone(diaga(j));
     value_t t2 =  nt2::oneminus(diaga(i))/(nt2::oneminus(diaga(j))*t1);
     value_t t = (nt2::rand(nt2::meta::as_<value_t>()) >  Half<value_t>()) ? t1 : t2;//Choose randomly from the two roots.
     value_t c =  nt2::rsqrt(nt2::oneplus(nt2::sqr(t))) ;
     value_t s =  t*c ;
     BOOST_AUTO_TPL( ij, nt2::cath(i, j));
     BOOST_AUTO_TPL(rot, nt2::catv(nt2::cath(c, s), nt2::cath(-s, c)));
     aa(nt2::_,ij)  =   nt2::mtimes(aa(nt2::_,ij), rot) ;
     aa(ij, nt2::_) =   nt2::mtimes(nt2::transpose(rot), aa(ij, nt2::_));
     aa(i, i) = nt2::One<value_t>();  //Ensure (i,i) element is exactly 1.
   } while (true);
   out =  aa;
   return out;
 }

 result_type operator()(A0& out, const A1& in) const
 {
   BOOST_AUTO_TPL(x,boost::proto::child_c<0>(in));
   std::size_t k = boost::proto::child_c<1>(in);
   std::size_t m = boost::proto::child_c<2>(in);
   std::size_t n = numel(x);
   BOOST_AUTO_TPL(a, nt2::expand(nt2::from_diag(x), nt2::of_size(m, n)));
   table<value_t> aa;
   if (!k)
   {
     // A = U*A*V where U and V are two random orthogonal matrices.
     aa = nt2::transpose(nt2::qmult(nt2::transpose(qmult(a))));
   }
   else
   {
     aa = a;
   }
   table<value_t> suma = nt2::sum(nt2::sqr(aa));
   do
   {
     nt2::table<ptrdiff_t> y = nt2::find(nt2::lt(suma, nt2::One<value_t>()));
     nt2::table<ptrdiff_t> z = nt2::find(nt2::gt(suma, nt2::One<value_t>()));
     if(nt2::isempty(y) || nt2::isempty(z)) break;
     std::size_t i = y(nt2::iceil(nt2::rand(nt2::meta::as_<value_t>())*nt2::numel(y)));
     std::size_t j = z(nt2::iceil(nt2::rand(nt2::meta::as_<value_t>())*nt2::numel(z)));
     if (i > j) std::swap(i, j);
     value_t aij = nt2::dot(aa(nt2::_, i), aa(nt2::_, j));
     value_t alpha = nt2::sqrt(nt2::sqr(aij)- nt2::minusone(suma(i))*nt2::minusone(suma(j)));
     value_t t1 = (aij + nt2::signnz(aij)*alpha)/nt2::minusone(suma(j));
     value_t t2 =  nt2::oneminus(suma(i))/(nt2::oneminus(suma(j))*t1);
     value_t t = (nt2::rand(nt2::meta::as_<value_t>()) >  Half<value_t>()) ? t1 : t2;
     value_t c =  nt2::rsqrt(nt2::oneplus(nt2::sqr(t))) ;
     value_t s =  t*c ;
     BOOST_AUTO_TPL( ij, nt2::cath(i, j));
     aa(nt2::_,ij) =   nt2::mtimes(aa(nt2::_,ij), nt2::catv(nt2::cath(c, s), nt2::cath(-s, c))) ;
     suma(j) +=  nt2::minusone(suma(i));
     suma(i) = nt2::One<value_t>();
   } while (true);
   out =  aa;
   return out;
 }

    static inline bool apply(Polygon1 const& poly1, Polygon2& poly2,
                Strategy const& strategy)
    {
        typedef typename ring_type<Polygon1>::type ring1_type;
        typedef typename ring_type<Polygon2>::type ring2_type;
        typedef typename point_type<Polygon2>::type point2_type;

        geometry::clear(poly2);

        if (!transform_range_out<point2_type>(exterior_ring(poly1),
                    std::back_inserter(exterior_ring(poly2)), strategy))
        {
            return false;
        }

        // Note: here a resizeable container is assumed.
        traits::resize
            <
                typename boost::remove_reference
                <
                    typename traits::interior_mutable_type<Polygon2>::type
                >::type
            >::apply(interior_rings(poly2), num_interior_rings(poly1));

        typename interior_return_type<Polygon1 const>::type rings1
                    = interior_rings(poly1);
        typename interior_return_type<Polygon2>::type rings2
                    = interior_rings(poly2);
        BOOST_AUTO_TPL(it1, boost::begin(rings1));
        BOOST_AUTO_TPL(it2, boost::begin(rings2));
        for ( ; it1 != boost::end(interior_rings(poly1)); ++it1, ++it2)
        {
            if (!transform_range_out<point2_type>(*it1,
                std::back_inserter(*it2), strategy))
            {
                return false;
            }
        }

        return true;
    }

    static inline void apply(std::basic_ostream<Char, Traits>& os,
        Polygon const& polygon, std::string const& style, int size)
    {
        typedef typename geometry::ring_type<Polygon>::type ring_type;
        typedef typename boost::range_iterator<ring_type const>::type iterator_type;

        bool first = true;
        os << "<g fill-rule=\"evenodd\"><path d=\"";

        ring_type const& ring = geometry::exterior_ring(polygon);
        for (iterator_type it = boost::begin(ring);
            it != boost::end(ring);
            ++it, first = false)
        {
            os << (first ? "M" : " L") << " "
                << geometry::get<0>(*it)
                << ","
                << geometry::get<1>(*it);
        }

        // Inner rings:
        {
            typename interior_return_type<Polygon const>::type rings
                        = interior_rings(polygon);
            for (BOOST_AUTO_TPL(rit, boost::begin(rings));
                rit != boost::end(rings); ++rit)
            {
                first = true;
                for (BOOST_AUTO_TPL(it, boost::begin(*rit)); it != boost::end(*rit);
                    ++it, first = false)
                {
                    os << (first ? "M" : " L") << " "
                        << geometry::get<0>(*it)
                        << ","
                        << geometry::get<1>(*it);
                }
            }
        }
        os << " z \" style=\"" << style << "\"/></g>";

    }