Boost-Commit :

Date view	Thread view	Subject view	Author view

Subject: [Boost-commit] svn:boost r51174 - in sandbox/improved_fast_gauss_transform: . boost boost/math boost/math/ifgt boost/math/ifgt/detail boost/utility libs libs/math libs/math/ifgt libs/math/ifgt/doc libs/math/ifgt/src libs/math/ifgt/src/example
From: erwann.rogard_at_[hidden]
Date: 2009-02-10 04:15:21

Next message: pbristow_at_[hidden]: "[Boost-commit] svn:boost r51175 - sandbox/SOC/2007/visualization/libs/svg_plot/test"
Previous message: erwann.rogard_at_[hidden]: "[Boost-commit] svn:boost r51173 - sandbox/miscellanea_algorithm/boost/algorithm"

Author: e_r
Date: 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
New Revision: 51174
URL: http://svn.boost.org/trac/boost/changeset/51174

Log:
Adding package improved_fast_gauss_transform
Added:
   sandbox/improved_fast_gauss_transform/
   sandbox/improved_fast_gauss_transform/boost/
   sandbox/improved_fast_gauss_transform/boost/math/
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/accumulator_base.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/benchmark.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/call_gauss_transform.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/call_nadaraya_watson_estimate.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/call_rozenblatt_parzen_estimate.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/call_wrapper.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/cluster.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/cluster_call_center.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/cluster_evaluator.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/clusters_evaluator.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/coefficients.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/coefficients_evaluator.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/cutoff_radius_none.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/cutoff_radius_rydg.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/detail/
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/detail/divider.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/detail/product_pdf.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/evaluator_base.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/evaluator_common.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/exact_accumulator.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/exact_contribution.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/exact_contributions_evaluator.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/exact_evaluator.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/fast_accumulator.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/fast_evaluator.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/find_nearest_cluster.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/for_each_accumulate.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/for_each_accumulate_rest_weights.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/for_each_evaluate.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/for_each_rozenblatt_parzen_estimate.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/include.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/normal_kernel_properties.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/optimal_bandwidth.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/optimal_parameter_given_max_degree.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/rest_weights_wrapper.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/tag.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/traits.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/truncation_degree_constant.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/truncation_properties.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/math/ifgt/zscore.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/utility/
   sandbox/improved_fast_gauss_transform/boost/utility/dont_care.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/boost/utility/nested_type.hpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/libs/
   sandbox/improved_fast_gauss_transform/libs/math/
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/doc/
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/doc/readme.txt (contents, props changed)
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/benchmark_exact.cpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/benchmark_exact.h (contents, props changed)
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/benchmark_fast.cpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/benchmark_fast.h (contents, props changed)
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/cluster.cpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/cluster.h (contents, props changed)
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/cluster_evaluator.cpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/cluster_evaluator.h (contents, props changed)
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/clusters_evaluator.cpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/clusters_evaluator.h (contents, props changed)
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/coefficients.cpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/coefficients.h (contents, props changed)
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/coefficients_evaluator.cpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/coefficients_evaluator.h (contents, props changed)
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/exact_gauss_transform.cpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/exact_gauss_transform.h (contents, props changed)
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/fast_accumulator.cpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/fast_accumulator.h (contents, props changed)
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/fast_evaluator.cpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/fast_evaluator.h (contents, props changed)
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/find_nearest_cluster.cpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/find_nearest_cluster.h (contents, props changed)
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/for_each_accumulate.cpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/for_each_accumulate.h (contents, props changed)
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/for_each_evaluate.cpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/for_each_evaluate.h (contents, props changed)
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/main.cpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/truncation_degree_constant.cpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/truncation_degree_constant.h (contents, props changed)
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/truncation_properties.cpp (contents, props changed)
   sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/truncation_properties.h (contents, props changed)

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/accumulator_base.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/accumulator_base.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,23 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/accumulator_base.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_ACCUMULATOR_BASE_HPP_ER_2009
+#define BOOST_MATH_IFGT_ACCUMULATOR_BASE_HPP_ER_2009
+
+namespace boost{namespace math{namespace ifgt{
+
+ template<typename Derived>
+ struct accumulator_base{};
+
+ // For now these two are just informative
+ template<typename Derived>
+ struct exact_accumulator_base : accumulator_base<Derived>{};
+
+ template<typename Derived>
+ struct fast_accumulator_base : accumulator_base<Derived>{};
+
+}}}
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/benchmark.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/benchmark.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,321 @@
+///////////////////////////////////////////////////////////////////////////////
+// benchmark.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_BENCHMARK_HPP_ER_2009
+#define BOOST_MATH_IFGT_BENCHMARK_HPP_ER_2009
+#include <ostream>
+#include <string>
+#include <vector>
+#include <algorithm>
+#include <functional>
+#include <numeric>
+#include <iterator>
+#include <boost/timer.hpp>
+#include <boost/static_assert.hpp>
+#include <boost/assign/std/vector.hpp>
+#include <boost/math/distributions/normal.hpp>
+#include <boost/random/mersenne_twister.hpp>
+#include <boost/random/normal_distribution.hpp>
+#include <boost/random/variate_generator.hpp>
+#include <boost/ref.hpp>
+#include <boost/function.hpp>
+#include <boost/format.hpp>
+#include <boost/algorithm/linfinity_distance.hpp>
+#include <boost/algorithm/l2_distance_squared.hpp>
+#include <boost/iterator/vector2matrix_iterator.hpp>
+#include <boost/algorithm/l2_norm.hpp>
+#include <boost/math/ifgt/include.hpp>
+#include <boost/math/ifgt/detail/product_pdf.hpp>
+namespace boost{namespace math{namespace ifgt{
+template<
+typename RealType,
+std::size_t test_count,
+std::size_t SourceSize,
+std::size_t WeightSize = 2,
+typename Pdf_fun = product_pdf<math::normal_distribution<RealType> >,
+typename W_fun = l2_norm
+>
+class benchmark{
+ //For now the implementation of benchmark only allows
+ BOOST_STATIC_ASSERT(WeightSize == 2); //hint :W_fun = l2_norm
+
+ public:
+ typedef RealType value_type;
+ typedef std::vector<value_type> var_type;
+ typedef linfinity_distance::type linf_dist_type;
+ typedef l2_distance_squared::type l2_dist_sq_type;
+ typedef Pdf_fun pdf_fun_type;
+ typedef W_fun w_fun_type;
+ typedef timer timer_type;
+
+ typedef std::size_t size_type;
+ typedef std::vector<size_type> sizes_type;
+
+ static const size_type sdim = SourceSize;
+ static const size_type wdim = WeightSize;
+
+ template<typename G>
+ benchmark(const G& g,const pdf_fun_type& pdf_fun,const w_fun_type& w_fun)
+ :
+ cumulative_train_count_(1),
+ train_x_(0),//(train_count_ * sdim),
+ train_w_(0),//(train_count_ * (wdim-1)),
+ test_x_(test_count * sdim),
+ test_pdf_(test_count * 1),
+ test_w_(test_count * (wdim-1)),
+ est_pdf_(size(test_pdf_)),
+ est_w_(size(test_w_)),
+ pdf_fun_(pdf_fun),
+ w_fun_(w_fun),
+ acc_time_per_1e3_train_count_(0),
+ eval_rp_time_per_1e3_test_count_(0),
+ eval_nw_time_per_1e3_test_count_(0),
+ e0_rp_((value_type)(0)),
+ e0_nw_((value_type)(0)),
+ e1_rp_((value_type)(0)),
+ e1_nw_((value_type)(0)),
+ linf_dist_(),
+ l2_dist_sq_()
+ {
+ generate_data(g,size(test_x_),test_x_);
+ update_pdf(test_x_,test_pdf_);
+ update_w(test_x_,test_w_);
+ }
+
+ benchmark(const benchmark& that)
+ :
+ cumulative_train_count_(that.cumulative_train_count_),
+ train_x_(that.train_x_),
+ train_w_(that.train_w_),
+ test_x_(that.test_x_),
+ test_pdf_(that.test_pdf_),
+ test_w_(that.test_w),
+ est_pdf_(that.est_pdf_),
+ est_w_(that.est_w_),
+ pdf_fun_(that.pdf_fun_),
+ w_fun_(that.w_fun_),
+ acc_time_per_1e3_train_count_(that.acc_time_per_1e3_train_count_),
+ eval_rp_time_per_1e3_test_count_(
+ that.eval_rp_time_per_1e3_test_count_
+ ),
+ eval_nw_time_per_1e3_test_count_(
+ that.eval_nw_time_per_1e3_test_count_
+ ),
+ e0_rp_(that.e0_rp_),
+ e0_nw_(that.e0_nw_),
+ e1_rp_(that.e1_rp_),
+ e1_nw_(that.e1_nw),
+ linf_dist_(that.linf_dist_),
+ l2_dist_sq_(that.l2_dist_sq_)
+ {}
+
+ benchmark& operator=(const benchmark& that){
+ if(&that!=this){
+ cumulative_train_count_ = that.cumulative_train_count_;
+ train_x_ = that.train_x_;
+ train_w_ = that.train_w_;
+ test_x_ = that.test_x_;
+ test_pdf_ = that.test_pdf_;
+ test_w_ = that.test_w;
+ est_pdf_ = that.est_pdf_;
+ est_w_ = that.est_w_;
+ pdf_fun_ = that.pdf_fun_;
+ w_fun_ = that.w_fun_;
+ acc_time_per_1e3_train_count_
+ = that.acc_time_per_1e3_train_count_;
+ eval_rp_time_per_1e3_test_count_
+ = that.eval_rp_time_per_1e3_test_count_;
+ eval_nw_time_per_1e3_test_count_
+ = that.eval_nw_time_per_1e3_test_count_;
+ e0_rp_ = that.e0_rp_;
+ e0_nw_ = that.e0_nw_;
+ e1_rp_ = that.e1_rp_;
+ e1_nw_ = that.e1_nw_;
+ linf_dist_ = that.linf_dist_;
+ l2_dist_sq_ = that.l2_dist_sq_;
+ }
+ return *this;
+ }
+
+ template<typename Acc,typename G>
+ void accumulate(G& g,size_type n,Acc& acc){
+ train_x_.resize(sdim * n);
+ train_w_.resize((wdim-1)*n);
+ generate_data(g,n * sdim,train_x_);
+ update_w(train_x_,train_w_);
+ // accumulate
+ timer_.restart();
+ ifgt::for_each(
+ train_x_,
+ ifgt::make_rest_weights_wrapper(train_w_),//w[1] --> {w[0]=1,w[1]}
+ acc
+ );
+ acc_time_per_1e3_train_count_.push_back(
+ timer_.elapsed() * 1e3 / n
+ );
+ std::size_t cn = cumulative_train_count_.back();
+ cumulative_train_count_.push_back(n+cn);
+ }
+ template<typename Eval>
+ void estimate_pdf(Eval& eval){//by Rozenblatt--Parzen
+ timer_.restart();
+ ifgt::for_each(test_x_,est_pdf_,eval,
+ ifgt::call<ifgt::rozenblatt_parzen_estimate>());
+ eval_rp_time_per_1e3_test_count_.push_back(
+ timer_.elapsed() * 1e3/test_count
+ );
+
+ value_type e0 = linf_dist_(test_pdf_,est_pdf_);
+ value_type e1 = sqrt(l2_dist_sq_(test_pdf_,est_pdf_))
+ /((value_type)(test_count));
+
+ e0_rp_.push_back(e0);
+ e1_rp_.push_back(e1);
+ }
+
+ template<typename Eval>
+ void estimate_w(Eval& eval){//by Nadaraya--Watson
+ timer_.restart();
+ ifgt::for_each(test_x_,est_w_,eval,
+ ifgt::call<ifgt::nadaraya_watson_estimate>());
+ eval_nw_time_per_1e3_test_count_.push_back(
+ timer_.elapsed() * 1e3/test_count
+ );
+ value_type e0 = linf_dist_(test_w_,est_w_);
+ value_type e1 = sqrt(l2_dist_sq_(test_w_,est_w_))
+ /((value_type)(test_count));
+
+ e0_nw_.push_back(e0);
+ e1_nw_.push_back(e1);
+ }
+
+ void notation(std::ostream& os)const{
+ os << "M: # evaluated" << std::endl;
+ os << "N: # accumulated" << std::endl;
+ os << "err0(a,b) = max {|a - b|:m=0,...,M-1 }" << std::endl;
+ os << "err1(a,b) = sqrt sum {|a - b|^2:m=0,...,M-1 } / M"<< std::endl;
+ os << "e0_rp : err0(rp(y), pdf(y)) " << std::endl;
+ os << "e1_rp : err1(rp(y), pdf(y)) " << std::endl;
+ os << "e0_nw : err0(nw(y), w[1](y)) " << std::endl;
+ os << "e1_nw : err1(nw(y), w[1](y)) " << std::endl;
+ os << "Times are per N=M=1e3" << std::endl;
+ os << "ta : time accumulate {(x,w):i=0,...,N-1}" << std::endl;
+ os << "trp : time {rp(y):m = 0,..., M-1} " << std::endl;
+ os << "trp : time {nw(y):m = 0,..., M-1} " << std::endl;
+ }
+
+ void header(std::ostream& os)const{
+ std::string str = "%|1$-10| %|2$-10| %|3$-10| %|4$-10| %|5$-10|";
+ str+= "%|6$-10| %|7$-10| %|8$-10| ";
+ os << boost::format(str)
+ % "N"
+ % "e0_rp"
+ % "e1_rp"
+ % "e0_nw"
+ % "e1_nw"
+ % "ta"
+ % "trp"
+ % "tnw";
+ }
+ void statistics(std::ostream& os, unsigned i){
+ BOOST_ASSERT(size(e0_rp_)-i>0);
+ std::string str = "%|1$-10.3| %|2$-10.3e| %|3$-10.3e|";
+ str+= "%|4$-10.3e| %|5$-10.3e| %|6$-10.3e|";
+ str+= "%|7$-10.3e| %|8$-10.3e|";
+ os << boost::format(str)
+ % cumulative_train_count_[i+1] // 1
+ % e0_rp_[i] // 2
+ % e1_rp_[i] // 3
+ % e0_nw_[i] // 4
+ % e1_nw_[i] // 5
+ % acc_time_per_1e3_train_count_[i] // 6
+ % eval_rp_time_per_1e3_test_count_[i] // 7
+ % eval_nw_time_per_1e3_test_count_[i] ; // 8
+ }
+ private:
+ template<typename G,typename R>
+ void generate_data(G& g, size_type n,R& x){
+ std::generate_n(
+ begin(x),
+ n,
+ g
+ );
+ };
+ template<typename R0,typename R1>
+ void update_pdf(const R0& x,R1& pdf_out){
+ std::transform(
+ make_vector2matrix_iterator(begin(x),sdim),
+ make_end_vector2matrix_iterator(
+ begin(x),
+ end(x),
+ sdim
+ ),
+ begin(pdf_out),
+ bind<value_type>(
+ pdf_fun_,
+ _1
+ )
+ );
+ }
+ template<typename R0,typename R1>
+ void update_w(const R0& x, R1& w_out){
+ std::transform(
+ make_vector2matrix_iterator(begin(x),sdim),
+ make_end_vector2matrix_iterator(
+ begin(x),
+ end(x),
+ sdim
+ ),
+ begin(w_out),
+ bind<value_type>(
+ w_fun_,
+ _1
+ )
+ );
+ }
+
+ sizes_type cumulative_train_count_;
+
+ // train data set
+ //w[1] = f(x) + e. here e=0, f(x) = ||x|| = (x[0]^2+...+x[D-1]^2)^(1/2)
+ var_type train_x_;
+ var_type train_w_;
+
+ // test data
+ var_type test_x_;
+ var_type test_pdf_;
+ var_type test_w_;
+
+ // estimates
+ var_type est_pdf_; // estimate of pdf
+ var_type est_w_; // estimate of {E[w[j]|x]:j>0}
+
+ //
+ pdf_fun_type pdf_fun_;
+ w_fun_type w_fun_;
+
+ //
+ timer_type timer_;
+
+ //
+ var_type acc_time_per_1e3_train_count_ ;
+ var_type eval_rp_time_per_1e3_test_count_;
+ var_type eval_nw_time_per_1e3_test_count_;
+
+ //
+ var_type e0_rp_;
+ var_type e0_nw_;
+ var_type e1_rp_;
+ var_type e1_nw_;
+
+ linf_dist_type linf_dist_;
+ l2_dist_sq_type l2_dist_sq_;
+
+};
+
+}}}
+
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/call_gauss_transform.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/call_gauss_transform.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,38 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/call_gauss_transform.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_CALL_GAUSS_TRANSFORM_HPP_ER_2009
+#define BOOST_MATH_IFGT_CALL_GAUSS_TRANSFORM_HPP_ER_2009
+#include <stdexcept>
+#include <boost/math/ifgt/evaluator_base.hpp>
+namespace boost{namespace math{namespace ifgt{
+
+ template<typename Evaluator>
+ class gauss_transform{
+ public:
+ typedef Evaluator evaluator_type;
+ struct stride_traits{
+ typedef mpl::size_t<evaluator_type::source_size> stride0_type;
+ typedef mpl::size_t<evaluator_type::weight_size> stride1_type;
+ };
+ typedef typename evaluator_type::value_type value_type;
+ template<typename Common>
+ gauss_transform(const evaluator_base<Common,Evaluator>& e)
+ :eval_(static_cast<const Evaluator&>(e)){}
+ gauss_transform(const gauss_transform& that) : eval_(that.eval){}
+
+ template<typename R0,typename R1>
+ void operator()(const R0& targets,R1& ranges_out)const{
+ eval_.gauss_transform(targets,ranges_out);
+ };
+
+ private:
+ gauss_transform& operator=(const gauss_transform& that);
+ const Evaluator& eval_;
+ };
+
+}}}
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/call_nadaraya_watson_estimate.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/call_nadaraya_watson_estimate.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,49 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/call_nadaraya_watson_estimate.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_CALL_NADARAYA_WATSON_ESTIMATE_HPP_ER_2009
+#define BOOST_MATH_IFGT_CALL_NADARAYA_WATSON_ESTIMATE_HPP_ER_2009
+#include <stdexcept>
+#include <boost/mpl/size_t.hpp>
+#include <boost/math/ifgt/evaluator_base.hpp>
+namespace boost{namespace math{namespace ifgt{
+
+template<typename Evaluator>
+class nadaraya_watson_estimate{
+ public:
+ typedef Evaluator evaluator_type;
+ struct stride_traits{
+ typedef mpl::size_t<evaluator_type::source_size> stride0_type;
+ typedef mpl::size_t<evaluator_type::weight_size-1> stride1_type;
+ };
+ typedef typename evaluator_type::value_type value_type;
+
+ template<typename Common>
+ nadaraya_watson_estimate(const evaluator_base<Common,Evaluator>& e)
+ :eval_(static_cast<const Evaluator&>(e)){}
+ nadaraya_watson_estimate(const nadaraya_watson_estimate& that)
+ : eval_(that.eval){}
+
+ template<typename R0,typename R1>
+ void operator()(const R0& targets,R1 ranges_out)const{
+ eval_.nadaraya_watson_estimate(targets,ranges_out);
+ };
+
+ private:
+ nadaraya_watson_estimate& operator=(
+ const nadaraya_watson_estimate& that);
+ const Evaluator& eval_;
+};
+
+template<typename C,typename E>
+nadaraya_watson_estimate<E>
+make_nadaraya_watson_estimate(const evaluator_base<C,E>& e){
+ typedef nadaraya_watson_estimate<E> result_type;
+ return result_type(e);
+}
+
+}}}
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/call_rozenblatt_parzen_estimate.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/call_rozenblatt_parzen_estimate.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,44 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/call_rozenblatt_parzen_estimate.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_CALL_ROZENBLATT_PARZEN_HPP_ER_2009
+#define BOOST_MATH_IFGT_CALL_ROZENBLATT_PARZEN_HPP_ER_2009
+#include <stdexcept>
+#include <boost/math/ifgt/evaluator_base.hpp>
+namespace boost{namespace math{namespace ifgt{
+
+template<typename Evaluator>
+struct rozenblatt_parzen_estimate{
+ typedef typename Evaluator::value_type value_type;
+ public:
+ template<typename Common>
+ rozenblatt_parzen_estimate(const evaluator_base<Common,Evaluator>& e)
+ :eval_(static_cast<const Evaluator&>(e)){}
+
+ rozenblatt_parzen_estimate(
+ const rozenblatt_parzen_estimate& that)
+ : eval_(that.eval_){}
+
+ template<typename R0>
+ value_type operator()(const R0& target)const{
+ return eval_.rozenblatt_parzen_estimate(target);
+ };
+
+ private:
+ rozenblatt_parzen_estimate& operator=(
+ const rozenblatt_parzen_estimate& that);
+ const Evaluator& eval_;
+};
+
+template<typename C,typename E>
+rozenblatt_parzen_estimate<E>
+make_rozenblatt_parzen_estimate(const evaluator_base<C,E>& e){
+ typedef rozenblatt_parzen_estimate<E> result_type;
+ return result_type(e);
+}
+
+}}}
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/call_wrapper.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/call_wrapper.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,19 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/call_wrapper.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_CALL_WRAPPER_HPP_ER_2009
+#define BOOST_MATH_IFGT_CALL_WRAPPER_HPP_ER_2009
+#include <boost/mpl/assert.hpp>
+namespace boost{namespace math{namespace ifgt{
+
+ /// For use by for_each_evaluator
+ template<template<typename> class Method>
+ struct call{
+ call(){};
+ };
+
+}}}
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/cluster.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/cluster.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,306 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/cluster.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_CLUSTER_HPP_ER_2009
+#define BOOST_MATH_IFGT_CLUSTER_HPP_ER_2009
+#include <algorithm>
+#include <functional>
+#include <iterator>
+#include <vector>
+#include <cmath>
+#include <boost/range.hpp>
+#include <boost/assert.hpp>
+#include <boost/static_assert.hpp>
+#include <boost/bind.hpp>
+#include <boost/iterator/zip_iterator.hpp>
+#include <boost/iterator/iterator_traits.hpp>
+#include <boost/tuple/tuple.hpp>
+#include <boost/algorithm/l2_norm.hpp>
+#include <boost/math/monomials.hpp>
+#include <boost/math/multi_indexes_derived.hpp>
+#include <boost/math/ifgt/zscore.hpp>
+#include <boost/math/ifgt/tag.hpp>
+#include <boost/math/ifgt/coefficients.hpp>
+#include <boost/math/ifgt/truncation_degree_constant.hpp>
+#include <boost/math/ifgt/normal_kernel_properties.hpp>
+namespace boost{namespace math{namespace ifgt{
+
+ // A cluster is characterized by
+ // - A center
+ // - A bandwidth
+ // - A collection of coefficients (one for each weight sequence)
+
+ /// apply<TruncationDegreePolicy,T>::type
+ template<
+ unsigned SourceSize,
+ unsigned WeightSize,
+ typename TruncationDegreePolicy
+ = truncation_degree_constant<mpl::_1>,
+ typename Cont = std::vector<double> >
+ class cluster{
+ typedef monomials<Cont> monoms_type;
+ typedef typename monoms_type::result_type monoms_result_type;
+ typedef coefficients<SourceSize,Cont> coeffs_type;
+ typedef std::vector<coeffs_type> coll_coeffs_type;
+ typedef typename coll_coeffs_type::iterator coll_coeffs_iter_type;
+
+ public:
+ // C++ Standard (now) allows initialization of const static data
+ // members of an integral type inside their class.
+ static const std::size_t source_size = SourceSize;
+ static const std::size_t weight_size = WeightSize;
+ typedef Cont center_type;
+ typedef const center_type& result_of_center_type;
+ typedef typename monoms_type::value_type value_type;
+ typedef coeffs_type coefficients_type;
+ typedef coll_coeffs_type collection_coefficients_type;
+ typedef const collection_coefficients_type&
+ result_of_collection_coefficients_type;
+ typedef std::size_t sources_count_type;
+ typedef normal_kernel_properties<source_size,value_type> rpp_type;
+ typedef typename mpl::apply<TruncationDegreePolicy,value_type>::type
+ truncation_degree_policy_type;
+
+ template<typename ArgPack>
+ cluster(const ArgPack& arg)
+ :center_(source_size),
+ bandwidth_(arg[tag::bandwidth]),
+ max_radius_(arg[tag::max_cluster_radius]),
+ truncation_degree_policy_(arg),
+ inv_nc_((value_type)(1)/rpp_type::normalizing_constant(bandwidth())),
+ source_radius_((value_type)(0)),
+ radius_((value_type)(0)),
+ sources_count_(0),
+ coll_coeffs_(weight_size){
+ BOOST_ASSERT(size(arg[tag::center])-source_size==0);
+ copy(
+ begin(arg[tag::center]),
+ end(arg[tag::center]),
+ begin(center_)
+ );
+ }
+
+ template<typename R>
+ cluster(
+ const R& center,
+ value_type bandwidth_val,//cannot name it bandiwth coz mf
+ value_type max_radius,
+ const truncation_degree_policy_type& truncation_degree_policy)
+ :center_(source_size),
+ bandwidth_(bandwidth_val),
+ max_radius_(max_radius),
+ truncation_degree_policy_(truncation_degree_policy),
+ inv_nc_((value_type)(1)/rpp_type::normalizing_constant(bandwidth())),
+ source_radius_((value_type)(0)),
+ radius_((value_type)(0)),
+ sources_count_(0),
+ coll_coeffs_(weight_size){
+ BOOST_ASSERT(size(center)-source_size==0);
+ copy(
+ begin(center),
+ end(center),
+ begin(center_)
+ );
+ }
+
+ cluster(const cluster& that)
+ :center_(that.source_size),
+ bandwidth_(that.bandwidth_),
+ max_radius_(that.max_radius_),
+ truncation_degree_policy_(that.truncation_degree_policy_),
+ inv_nc_(that.inv_nc_),
+ source_radius_(that.source_radius_),
+ radius_(that.radius_),
+ sources_count_(that.sources_count_),
+ coll_coeffs_(that.coll_coeffs_){
+ copy(
+ begin(that.center_),
+ end(that.center_),
+ begin(center_)
+ );
+ }
+
+ cluster& operator=(const cluster& that){
+ if(&that!=this){
+ center_ = that.center_;
+ bandwidth_ = that.bandwidth_;
+ max_radius_ = that.max_radius_;
+ truncation_degree_policy_ = that.truncation_degree_policy_;
+ inv_nc_ = that.inv_nc_;
+ source_radius_ = that.source_radius_;
+ radius_ = that.radius_;
+ sources_count_ = that.sources_count_;
+ coll_coeffs_ = that.coll_coeffs_;
+ }
+ return *this;
+ }
+
+ result_of_center_type center()const{ return center_; }
+ value_type bandwidth()const{ return bandwidth_; }
+ /// Maximum allowed cluster radius
+ value_type max_radius()const{ return max_radius_; }
+
+ value_type inverse_normalizing_constant()const{ return inv_nc_; }
+
+ /// Maximum radius among collected sources
+ value_type radius()const{ return radius_; }
+ /// Number of sources collected
+ sources_count_type sources_count()const{ return sources_count_; }
+
+ /// Radius of the source in the last argument that was passed
+ value_type source_radius()const{ return source_radius_; }
+ bool source_radius_is_zero()const{
+ return source_radius() == (value_type)(0);
+ }
+ bool source_radius_exceeds_max_radius()const{
+ return (source_radius() > max_radius());
+ }
+
+ result_of_collection_coefficients_type
+ collection_coefficients()const{return coll_coeffs_; }
+
+ template<typename R0,typename R1>
+ bool operator()(const R0& source, const R1& weight){
+ typedef std::vector<unsigned> degrees_type;
+ typedef zscore<center_type> zscore_type;
+ typedef l2_norm_squared sqnorm_type;
+ static degrees_type degrees(weight_size,0);
+ //static zscore_type zscore(center());//not safe
+ static Cont zscore_val(source_size);
+ static sqnorm_type sqnorm;
+ zscore_type zscore(center());
+ zscore(source,bandwidth(),zscore_val);
+
+ value_type zscore_sqnorm = sqnorm(zscore_val);
+ source_radius_ = sqrt(zscore_sqnorm) * bandwidth();
+ bool do_collect = !(
+ source_radius_is_zero()
+ || source_radius_exceeds_max_radius() );
+
+ if(do_collect){
+ truncation_degree_policy_(
+ *this, source_radius(), weight, degrees);
+ (*this)(zscore_val, zscore_sqnorm, weight, degrees);
+ }
+ return do_collect;
+ }
+
+ // TODO check that
+ // range_iterator<R2>::type::value_type == unsigned
+ /// Returns true if source collected, false otherwise
+ template<typename R0,typename R1,typename R2>
+ void operator()(
+ const R0& zscore_val,
+ value_type zscore_sqnorm,
+ const R1& weight,
+ const R2& degrees
+ ){
+ typedef range_iterator<R0> ir0_type;
+ typedef range_iterator<R1> ir1_type;
+ typedef range_iterator<R2> ir2_type;
+ typedef typename ir0_type::type iter0_type;
+ typedef typename ir1_type::type iter1_type;
+ typedef typename ir2_type::type iter2_type;
+ typedef typename iter0_type::value_type value0_type;
+ typedef typename iter1_type::value_type value1_type;
+ typedef typename iter2_type::value_type value2_type;
+ BOOST_MPL_ASSERT((is_same<value0_type,value_type>));
+ BOOST_MPL_ASSERT((is_same<value1_type,value_type>));
+ BOOST_MPL_ASSERT((is_same<value2_type,unsigned>));
+ //static Cont zscore_val(source_size);
+ static Cont consts(source_size);
+ static monomials<Cont> monoms;
+ BOOST_ASSERT(size(zscore_val)-source_size==0);
+ BOOST_ASSERT(size(weight)-weight_size==0);
+ BOOST_ASSERT(size(degrees)-weight_size==0);
+ BOOST_ASSERT(zscore_sqnorm!=((value_type)(0)));
+
+ unsigned local_max_degree
+ = (*std::max_element(
+ begin(degrees),
+ end(degrees)));
+
+ monoms(zscore_val,local_max_degree);
+
+ for_each(
+ make_zip_iterator(
+ make_tuple(
+ begin(weight),
+ begin(coll_coeffs_),
+ begin(degrees)
+ )
+ ),
+ make_zip_iterator(
+ make_tuple(
+ end(weight),
+ end(coll_coeffs_),
+ end(degrees)
+ )
+ ),
+ op_factory::make(monoms(),zscore_sqnorm)
+ );
+
+ radius_
+ = (radius()<source_radius())
+ ? source_radius() : radius();
+
+ ++sources_count_;
+ }
+ //TODO remove. Necessary for temporary fix in copy constructor
+ // of derived classes
+ const truncation_degree_policy_type& truncation_degree_policy()const{
+ return truncation_degree_policy_;}
+
+ private:
+ template<typename Mons>
+ class op{
+ public:
+ typedef tuple<
+ const value_type&,
+ coeffs_type&,
+ const unsigned&> argument_type;
+ typedef void result_type;
+ typedef
+ typename range_iterator<const Mons>::type mons_ir_type;
+ op(const Mons& mons,value_type zscore_sqnorm)
+ :zscore_sqnorm_(zscore_sqnorm),mons_(mons){}
+ op(const op& that)
+ :zscore_sqnorm_(that.zscore_sqnorm_),mons_(that.mons_){}
+
+ void operator()(argument_type t){
+ value_type w = t.template get<0>();
+ coeffs_type& coeffs = t.template get<1>();
+ unsigned degree = t.template get<2>();
+ coeffs(zscore_sqnorm_,mons_,w,degree);
+ }
+ private:
+ op& operator=(const op& that);
+ value_type zscore_sqnorm_;
+ const Mons& mons_;
+ };//op
+ struct op_factory{
+ template<typename Mons>
+ static op<Mons> make(const Mons& mons,value_type zscore_sqnorm){
+ return op<Mons>(mons,zscore_sqnorm);
+ }
+ };
+ center_type center_;
+ value_type bandwidth_;
+ value_type max_radius_;
+ truncation_degree_policy_type truncation_degree_policy_;
+ value_type inv_nc_;
+
+ value_type source_radius_;
+ value_type radius_;
+ sources_count_type sources_count_;
+ coll_coeffs_type coll_coeffs_;
+ };
+
+
+}}}
+
+
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/cluster_call_center.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/cluster_call_center.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,21 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/cluster_call_center.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_CLUSTER_CALL_CENTER_HPP_ER_2009
+#define BOOST_MATH_IFGT_CLUSTER_CALL_CENTER_HPP_ER_2009
+namespace boost{namespace math{namespace ifgt{
+ template<typename Cc>
+ struct cluster_call_center{
+ typedef typename Cc::center_type center_type;
+ typedef const center_type& result_type;
+ typedef const Cc& argument_type;
+ cluster_call_center(){}
+ result_type operator()(argument_type arg)const{
+ return arg.center();
+ }
+ };
+}}}
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/cluster_evaluator.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/cluster_evaluator.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,141 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/cluster_evalutor.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_CLUSTER_EVALUATOR_HPP_ER_2009
+#define BOOST_MATH_IFGT_CLUSTER_EVALUATOR_HPP_ER_2009
+#include <stdexcept>
+#include <algorithm>
+#include <boost/range.hpp>
+#include <boost/bind.hpp>
+#include <boost/mpl/assert.hpp>
+#include <boost/mpl/apply.hpp>
+#include <boost/mpl/placeholders.hpp>
+#include <boost/utility.hpp>
+#include <boost/type_traits/remove_cv.hpp>
+#include <boost/type_traits/remove_reference.hpp>
+#include <boost/math/ifgt/tag.hpp>
+#include <boost/utility/dont_care.hpp>
+#include <boost/math/ifgt/coefficients_evaluator.hpp>
+#include <boost/iterator/range_same.hpp>
+
+namespace boost{ namespace math{namespace ifgt{
+
+ // cluster_evaluator takes a target and a range_out in its constructor
+ // remember that a cluster contains as many coefficients as
+ // weight sequences. The job of cluster_evaluator
+ // is to apply coefficients_evaluator to each
+ // coefficients and *add* the results to range_out
+
+ // \warning adds to range_out
+ // \warning The cutoff_radius_policy is passed by reference for efficiency
+ // R ForwarRange (that of the target)
+ // apply<CutoffRadiusPolicy,X>::type
+ template<
+ typename R,
+ typename RangeOut,
+ typename CutoffRadiusPolicy>
+ class cluster_evaluator{
+ public:
+ typedef typename range_value<R>::type value_type;
+ typedef RangeOut& ref_range_out_type;
+ private:
+ typedef typename mpl::apply<CutoffRadiusPolicy,value_type>::type
+ cutoff_radius_policy_type;
+ typedef const cutoff_radius_policy_type&
+ cref_cutoff_radius_policy_type;
+ public:
+ template<typename ArgPack>
+ cluster_evaluator(const ArgPack& args)
+ :range_out_(args[tag::range_out]),
+ coeffs_evaluator_(args[tag::target]),
+ cutoff_radius_policy_(args[tag::cutoff_radius_policy]){
+ // Do not replace by cutoff_radius_policy_(args)
+ }
+
+ //TODO remove
+ cluster_evaluator(
+ const R& target,
+ ref_range_out_type range_out,
+ cref_cutoff_radius_policy_type cutoff_radius_policy)
+ :range_out_(range_out),
+ coeffs_evaluator_(target),
+ cutoff_radius_policy_(cutoff_radius_policy){}
+
+ cluster_evaluator(const cluster_evaluator& that)
+ :range_out_(that.range_out_),
+ coeffs_evaluator_(that.coeffs_evaluator_),
+ cutoff_radius_policy_(that.cutoff_radius_policy_){}
+
+ template<typename Cluster>
+ void operator()(const Cluster& cluster){
+ static range_same policy;
+ return this->operator()(cluster, policy);
+ }
+
+ template<typename Cluster,typename SubsetRangePolicy>
+ void operator()(const Cluster& cluster,
+ const SubsetRangePolicy& subset_range_policy){
+ typedef typename Cluster::coefficients_type coeffs_type;
+ typedef SubsetRangePolicy subset_range_policy_type;
+ //compile warning:
+ //"comparison between signed and unsigned integer expressions|"
+ //why? source_size() is also a size_type
+ BOOST_ASSERT(
+ coeffs_evaluator_.target_size() - size(cluster.center()) == 0 );
+
+ value_type cutoff_radius = cutoff_radius_policy_(cluster);
+
+ typedef typename
+ Cluster::result_of_collection_coefficients_type arg_type;
+
+ typedef typename result_of<
+ subset_range_policy_type(arg_type)>::type subset_type;
+
+ subset_type subset
+ = subset_range_policy(cluster.collection_coefficients());
+
+ BOOST_ASSERT(size(subset) ==size(range_out_));
+
+ typename coeffs_type::result_type (coeffs_type::*f)() const =
+ &coeffs_type::operator();
+
+ std::transform(
+ begin(subset),
+ end(subset),
+ begin(range_out_),
+ begin(range_out_),
+ bind(
+ std::plus<value_type>(),
+ bind(
+ ref(coeffs_evaluator_),
+ cutoff_radius,
+ cref(cluster.center()),
+ cluster.bandwidth(),
+ cluster.inverse_normalizing_constant(),
+ bind(
+ f,
+ _1
+ ),
+ bind(
+ &coeffs_type::max_degree,
+ _1
+ )
+ ),
+ _2
+ )
+ );
+ }
+
+ private:
+ cluster_evaluator& operator=(const cluster_evaluator& that);
+ cluster_evaluator();
+ typedef coefficients_evaluator<R> coeffs_evaluator_type;
+ ref_range_out_type range_out_;
+ coeffs_evaluator_type coeffs_evaluator_;
+ cutoff_radius_policy_type cutoff_radius_policy_;
+ };
+}}}
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/clusters_evaluator.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/clusters_evaluator.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,89 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/clusters_evaluator.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_CLUSTERS_EVALUATOR_HPP_ER_2009
+#define BOOST_MATH_IFGT_CLUSTERS_EVALUATOR_HPP_ER_2009
+#include <algorithm>
+#include <boost/assert.hpp>
+#include <boost/static_assert.hpp>
+#include <boost/range.hpp>
+#include <boost/math/ifgt/cluster_evaluator.hpp>
+namespace boost{namespace math{namespace ifgt{
+
+ // Adds to a range the result from evaluating each cluster
+ template<
+ typename RealType,
+ typename CutoffRadiusPolicy
+ >
+ class clusters_evaluator{
+ public:
+ typedef RealType value_type;
+ typedef typename mpl::apply<CutoffRadiusPolicy,value_type>::type
+ cluster_radius_policy_type;
+ template<typename ArgPack>
+ clusters_evaluator(const ArgPack& args):cluster_radius_policy_(args){}
+
+ // make sure policy's default constructor is disabled
+ // if necessary
+ // so that compile error here
+ clusters_evaluator(){}
+
+ clusters_evaluator(const cluster_radius_policy_type&
+ cluster_radius_policy)
+ : cluster_radius_policy_(cluster_radius_policy){}
+ clusters_evaluator(const clusters_evaluator& that)
+ :cluster_radius_policy_(that.cluster_radius_policy_){}
+
+ clusters_evaluator& operator=(const clusters_evaluator& that){
+ if(&that!=this){
+ cluster_radius_policy_ = that.cluster_radius_policy_;
+ }
+ return *this;
+ }
+
+ template<typename R0,typename R1,typename R2>
+ void operator()(
+ const R0& target,const R1& clusters,R2& range_out)const{
+ static range_same policy;
+ return this->operator()(
+ target, clusters, range_out, policy
+ );
+ }
+
+ //\warning adds to range_out
+ template<typename R0,typename R1,typename R2,
+ typename SubsetRangePolicy>
+ void operator()(
+ const R0& target,const R1& clusters,R2& range_out,
+ const SubsetRangePolicy& subset_range_policy)const{
+ typedef typename range_value<R0>::type value0_type;
+ typedef typename range_value<R1>::type cluster_type;
+ typedef typename cluster_type::value_type value1_type;
+ typedef typename range_value<R2>::type value2_type;
+ BOOST_MPL_ASSERT((is_same<value0_type,value_type>));
+ BOOST_MPL_ASSERT((is_same<value1_type,value_type>));
+ BOOST_MPL_ASSERT((is_same<value2_type,value_type>));
+ typedef cluster_evaluator<R0,R2,CutoffRadiusPolicy> eval_type;
+
+ eval_type eval(target, range_out, cluster_radius_policy_);
+
+ for_each(
+ begin(clusters),
+ end(clusters),
+ bind<void>(
+ ref(eval),//not cref because eval modifies range_out
+ _1,
+ cref(subset_range_policy)
+ )
+ );
+ };
+
+ private:
+ cluster_radius_policy_type cluster_radius_policy_;
+ };
+
+}}}
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/coefficients.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/coefficients.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,150 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/coefficients.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_SOURCE_COEFFICIENTS_HPP_ER_2009
+#define BOOST_MATH_IFGT_SOURCE_COEFFICIENTS_HPP_ER_2009
+#include <algorithm>
+#include <functional>
+#include <iterator>
+#include <vector>
+#include <cmath>
+#include <boost/range.hpp>
+#include <boost/range/value_type.hpp>
+#include <boost/static_assert.hpp>
+#include <boost/bind.hpp>
+#include <boost/iterator/zip_iterator.hpp>
+#include <boost/iterator/iterator_traits.hpp>
+#include <boost/tuple/tuple.hpp>
+#include <boost/math/monomials.hpp>
+#include <boost/math/multi_indexes_derived.hpp>
+namespace boost{namespace math{namespace ifgt{
+
+ // coefficients keeps a set of coefficients and a max_degree
+ // Each time (zsqn, monoms, w,p) is passed to the object,
+ // the coefficients {C_\alpha:\alpha_0+...+\alpha_{SourceSz-1}<=p}
+ // are updated.
+ // Specifically,
+ // C_{\alpha} +=
+ // 2^{\alpha}/\alpha! * w exp(-zsqn ) * z^alpha
+ // The z^alpha's are the elements in monoms
+ // zsqn represents || z ||_2^2
+ //
+ // The reason for this set up (zsqn and monoms passed to the object)
+ // is that this operation may be repeated
+ // for different weight sequences (one coefficients for each).
+
+ template<unsigned SourceSz,typename Cont = std::vector<double> >
+ class coefficients{
+ typedef typename Cont::value_type value_type;
+ typedef multi_indexes_derived<SourceSz,multi_power2divfact_op>
+ mult_p2df;
+ typedef typename mult_p2df::iter_range_type p2df_ir_type;
+ typedef
+ typename range_iterator<p2df_ir_type>::type p2df_iter_type;
+ typedef
+ typename range_value<p2df_ir_type>::type p2df_value_type;
+ public:
+ typedef Cont coefficients_type;
+ typedef const coefficients_type& result_type;
+
+ coefficients():max_degree_(0){}
+ coefficients(const coefficients& that)
+ :max_degree_(that.max_degree_),coeffs_(that.coeffs_){}
+ coefficients& operator=(const coefficients& that){
+ if(&that!=this){
+ max_degree_ = that.max_degree_;
+ coeffs_ = that.coeffs_;
+ }
+ return *this;
+ }
+
+ unsigned max_degree()const{ return max_degree_; }
+
+ result_type operator()()const{ return coeffs_; }
+
+ template<typename R>
+ void operator()(
+ value_type zscore_sqnorm,
+ const R& monomials,
+ value_type weight,
+ unsigned degree
+ ){
+ typedef typename range_iterator<const R>::type monoms_iter_type;
+ typedef typename range_value<R>::type monoms_value_type;
+
+ value_type c = exp(-zscore_sqnorm) * weight;
+ p2df_ir_type p2df_ir = mult_p2df::get(degree);
+ std::size_t adv_by
+ = monomials_properties<>::number_degree_less_than(
+ degree,SourceSz);
+ if(adv_by>((std::size_t)(size(coeffs_)))){
+ coeffs_.resize(adv_by,(value_type)(0));
+ }
+
+ BOOST_ASSERT( adv_by <= ((std::size_t)(size(monomials))) );
+ BOOST_ASSERT( adv_by <= ((std::size_t)(size(p2df_ir))) );
+ monoms_iter_type m_e = begin(monomials);
+ p2df_iter_type p2df_e =
+ begin(p2df_ir);
+ std::advance(m_e,adv_by);
+ std::advance(p2df_e,adv_by);
+
+ inner_op<monoms_value_type> io(c);
+
+ transform(
+ make_zip_iterator(
+ make_tuple(
+ begin(monomials),
+ begin(p2df_ir)
+ )
+ ),
+ make_zip_iterator(
+ make_tuple(
+ m_e,
+ p2df_e
+ )
+ ),
+ begin(coeffs_),
+ begin(coeffs_),
+ bind(
+ std::plus<value_type>(),
+ bind(
+ io,//ref?
+ _1
+ ),
+ _2
+ )
+ );
+
+ if(max_degree()<degree){max_degree_ = degree;}
+ }
+
+ private:
+ template<typename V>
+ struct inner_op{
+ typedef tuple<
+ const V&,const p2df_value_type&> argument_type;
+
+ typedef V result_type;
+ explicit inner_op(result_type c_):c(c_){}
+ inner_op(const inner_op& that):c(that.c){}
+ inner_op& operator=(const inner_op& that){
+ if(&that!=this){ c=that.c;}
+ return *this;}
+ result_type operator()(argument_type t){
+ const V& mon = t.template get<0>();
+ const p2df_value_type& p2df = t.template get<1>();
+
+ return mon * ((result_type)(p2df)) * c;
+ }
+ result_type c;
+ };//inner_op
+
+ unsigned max_degree_;
+ coefficients_type coeffs_;
+ };
+}}}
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/coefficients_evaluator.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/coefficients_evaluator.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,154 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/coefficients_evaluator.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_COLLECTED_SOURCES_EVALUATOR_HPP_ER_2009
+#define BOOST_MATH_IFGT_COLLECTED_SOURCES_EVALUATOR_HPP_ER_2009
+#include <vector>
+#include <functional>
+#include <iterator>
+#include <cmath>
+#include <numeric>
+#include <stdexcept>
+#include <iostream> //TODO remove
+#include <boost/assert.hpp>
+#include <boost/static_assert.hpp>
+#include <boost/bind.hpp>
+#include <boost/range.hpp>
+#include <boost/math/ifgt/zscore.hpp>
+#include <boost/math/monomials_properties.hpp>
+namespace boost{namespace math{namespace ifgt{
+
+ // coefficients_evaluator keeps a ref to a target
+ // each time a (coefficients,center,bandwidth) is passed to it,
+ // it returns the gauss transform, if the target
+ // is within a cutoff radius, 0 otherwise
+ //
+ // This class is primarly intended for use within clusters_evaluator.hpp
+
+ /// R0 ForwardRange (that of the target)
+ template<typename R0>
+ class coefficients_evaluator{
+
+ typedef typename range_value<R0>::type value_type;
+ typedef std::vector<value_type> vec_type;
+ typedef typename range_size<R0>::type target_size_type;
+ public:
+ typedef value_type result_type;
+ coefficients_evaluator(const R0& target):target_(target){}
+
+ coefficients_evaluator(const coefficients_evaluator& that)
+ :target_(that.target_){}
+
+ const R0& target()const{return target_; }
+ target_size_type target_size()const{return size(target()); }
+
+ template<typename R1,typename Coeffs>
+ result_type operator()(
+ result_type cutoff_radius,
+ const R1& center,
+ result_type bandwidth,
+ result_type inv_nc,
+ const Coeffs& coeffs
+ ){
+ return (*this)(
+ cutoff_radius,
+ center,
+ bandwidth,
+ inv_nc,
+ coeffs(),
+ coeffs.max_degree()
+ );
+ }
+
+ template<typename R1,typename R2>
+ result_type operator()(
+ //ignore target if dist to center less than:
+ result_type cutoff_radius,
+ const R1& center,
+ result_type bandwidth,
+ result_type inv_nc,
+ const R2& coefficients,
+ unsigned degree
+ )const{
+ typedef typename range_value<R0>::type value0_type;
+ typedef typename range_value<R1>::type value1_type;
+ typedef typename range_value<R2>::type value2_type;
+ //TODO check that value0_type == ValueType, likewise 1,2
+ BOOST_MPL_ASSERT((is_same<value_type,value1_type>));
+ BOOST_MPL_ASSERT((is_same<value_type,value2_type>));
+ BOOST_ASSERT((target_size_type)(size(center))==target_size());
+ {
+ unsigned n = monomials_properties<>::number_degree_less_than(
+ degree,(unsigned)(target_size())
+ );
+ BOOST_ASSERT((size(coefficients)-n>=0)||(degree == 0));
+ }
+
+ result_type res = (result_type)(0);
+ if(degree>0){
+
+ static vec_type zscore_val(target_size());
+ static monomials<vec_type> monoms;
+ typedef zscore<vec_type> zscore_type;
+ {
+ zscore_type zscore(center);
+ zscore(target(),bandwidth,zscore_val);
+ }
+
+ result_type zscore_sqnorm = std::inner_product(
+ begin(zscore_val),
+ end(zscore_val),
+ begin(zscore_val),
+ (result_type)(0)
+ );
+
+ //TODO rename target_radius
+ result_type sq_d = zscore_sqnorm * (bandwidth * bandwidth);
+
+ // (sq_d<cutoff_radius * cutoff_radius) should suffice
+ // but (sq_d<cutoff_radius) necessary in case the client
+ // (usually cluster_evaluator) returns
+ // cutoff_radius = highest()
+ if( (sq_d<cutoff_radius)
+ || (sq_d<cutoff_radius * cutoff_radius) ){
+ result_type gauss_kernel
+ = inv_nc * exp(-zscore_sqnorm);
+ monoms(zscore_val,degree);
+ res = std::inner_product(
+ begin(coefficients),
+ end(coefficients),
+ begin(monoms()),
+ (result_type)(0),
+ std::plus<result_type>(),
+ bind(
+ std::multiplies<result_type>(),
+ gauss_kernel,
+ bind(
+ std::multiplies<result_type>(),
+ _1,
+ _2
+ )
+ )
+ );
+ }
+ }
+ return res;
+ }
+ private:
+ coefficients_evaluator& operator=(
+ const coefficients_evaluator& that);
+ const R0& target_;
+ };
+
+ template<typename R0>
+ coefficients_evaluator<R0>
+ make_coefficients_evaluator(const R0& target){
+ return coefficients_evaluator<R0>(target);
+ }
+
+}}}
+
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/cutoff_radius_none.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/cutoff_radius_none.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,32 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/cutoff_radius_non.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_CUTOFF_RADIUS_NONE_HPP_ER_2009
+#define BOOST_MATH_IFGT_CUTOFF_RADIUS_NONE_HPP_ER_2009
+#include <boost/numeric/conversion/bounds.hpp>
+#include <boost/limits.hpp>
+#include <boost/utility/dont_care.hpp>
+
+namespace boost{namespace math{namespace ifgt{
+ // Models CutoffRadiusPolicy
+ // This policy imposes no cutoff radius i.e. all clusters are used
+ // \warning the cutoff squared may be nan
+ template<typename RealType>
+ struct cutoff_radius_none{
+ typedef RealType result_type;
+ typedef RealType first_argument_type;
+ typedef RealType second_argument_type;
+ cutoff_radius_none(utility::dont_care){}
+ cutoff_radius_none(){}
+ template<typename Cluster>
+ result_type operator()(const Cluster& cluster)const{
+ static const result_type res
+ = numeric::bounds<result_type>::highest();
+ return res;
+ }
+ };
+}}}
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/cutoff_radius_rydg.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/cutoff_radius_rydg.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,82 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/cutoff_radius_rydg.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_CUTOFF_RADIUS_RYDG_HPP_ER_2009
+#define BOOST_MATH_IFGT_CUTOFF_RADIUS_RYDG_HPP_ER_2009
+#include <cmath>
+#include <boost/numeric/conversion/bounds.hpp>
+
+namespace boost{namespace math{namespace ifgt{
+
+ // Models CutoffRadiusPolicy
+ // Raykar2006a stands for Raykar et al. See Section 4.1:
+ // Automatically choosing the cut off radius for each cluster
+ template<typename RealType>
+ class cutoff_radius_rydg{
+ public:
+ typedef RealType result_type;
+ typedef RealType first_argument_type;
+ typedef RealType second_argument_type;
+ private:
+ static result_type default_max_distance_source_target(){
+ static const result_type r
+ = numeric::bounds<result_type>::highest();
+ return r;
+ }
+ public:
+
+ /// Suggestion: error_tolerance be in relation to the standard dev
+ /// of weight, sigma. Say sigma/1e3.
+ template<typename ArgPack>
+ cutoff_radius_rydg(const ArgPack& args)
+ :error_tol_(args[tag::error_tolerance]),
+ max_dist_source_target_(
+ args[tag::max_distance_source_target
+ |default_max_distance_source_target()]){}
+
+ cutoff_radius_rydg(
+ result_type error_tol,
+ result_type max_dist_source_target
+ ):error_tol_(error_tol),
+ max_dist_source_target_(max_dist_source_target){}
+
+ // default copy/assignement should work
+
+ /// epsilon (Raykar2006a)
+ result_type error_tolerance()const{ return error_tol_; }
+ /// R in Raykar2006a
+ result_type max_distance_source_target()const{
+ return max_dist_source_target_; }
+
+ /// Returns r_y^k (Raykar2006a)
+ template<typename Cluster>
+ result_type operator()(const Cluster& cluster)const{
+ return
+ this->operator()(cluster.bandwidth(),cluster.max_radius());
+ }
+
+ private:
+ /// Arguments: h and r_x^k (Raykar2006a)
+ result_type operator()(
+ result_type bandwidth,
+ result_type max_radius
+ )const{
+ static result_type c
+ = sqrt(log((result_type)(1)/error_tolerance()));
+
+ result_type a = max_distance_source_target();
+ result_type b = bandwidth * c;
+ result_type d = (a<b)? a : b;
+
+ return max_radius + d;
+ }
+ cutoff_radius_rydg();//dont implement
+ result_type error_tol_;
+ result_type max_dist_source_target_;
+ };
+
+}}}
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/detail/divider.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/detail/divider.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,31 @@
+//////////////////////////////////////////////////////////////////////////////
+// divider.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_DIVIDER_HPP_ER_2009
+#define BOOST_MATH_IFGT_DIVIDER_HPP_ER_2009
+namespace boost{namespace math{namespace ifgt{
+
+ template<typename RealType>
+ struct divider{
+ public:
+ typedef RealType value_type;
+ template<typename ArgPack>
+ divider(const ArgPack arg):x_(arg[tag::divide_by|(value_type)(2)]){}
+ divider(const divider& that):x_(that.x_){}
+ divider& operator=(const divider& that){
+ if(&that!=this){
+ x_ = that.x_;
+ }
+ return *this;
+ }
+
+ value_type operator()(value_type y){ return y/x_; }
+ private:
+ value_type x_;
+ };
+
+}}}
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/detail/product_pdf.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/detail/product_pdf.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,79 @@
+//////////////////////////////////////////////////////////////////////////////
+// product_pdf.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_PRODUCT_PDF_HPP_ER_2009
+#define BOOST_MATH_IFGT_PRODUCT_PDF_HPP_ER_2009
+#include <cmath>
+#include <numeric>
+#include <boost/function.hpp>
+#include <boost/bind.hpp>
+#include <boost/ref.hpp>
+#include <boost/math/distributions/chi_squared.hpp>
+namespace boost{namespace math{namespace ifgt{
+ //TODO not a part of ifgt per se, but used in some of the
+ // example/*.cpp files
+ template<typename Dist>
+ class product_pdf{
+ public:
+ typedef typename Dist::value_type result_type;
+ product_pdf(const Dist& dist):call_log_pdf_(dist){}
+ product_pdf(const product_pdf& that)
+ : call_log_pdf_(that.call_log_pdf_){}
+ product_pdf& operator=(const product_pdf& that){
+ if(&that!=this){
+ call_log_pdf_ = that.call_log_pdf_;
+ }
+ return *this;
+ }
+
+ template<typename R>
+ result_type operator()(const R& x)const{
+ result_type log_pdf = std::accumulate(
+ begin(x),
+ end(x),
+ (result_type)(0),
+ bind(
+ cref(call_log_pdf_),
+ _1,
+ _2
+ )
+ );
+ return exp(log_pdf);
+ };
+ private:
+ product_pdf();
+ class call_log_pdf{
+ public:
+ typedef typename product_pdf::result_type result_type;
+ typedef result_type first_argument_type;
+ typedef result_type second_argument_type;
+ call_log_pdf(const Dist& dist):dist_(dist){}
+ call_log_pdf(const call_log_pdf& that):dist_(that.dist_){}
+ call_log_pdf& operator()(const call_log_pdf& that){
+ if(&that!=this){
+ dist_ = that.dist_;
+ }
+ return dist_;
+ }
+
+ result_type operator()(result_type x,result_type y)const{
+ return x+log(pdf(dist_,y));
+ }
+ private:
+ call_log_pdf();
+ Dist dist_;
+ };
+ call_log_pdf call_log_pdf_;
+ };
+
+ template<typename Dist>
+ product_pdf<Dist> make_product_pdf(const Dist& dist){
+ return product_pdf<Dist>(dist);
+ };
+
+
+}}}
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/evaluator_base.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/evaluator_base.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,149 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/evaluator_base.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_EVALUATOR_BASE_HPP_ER_2009
+#define BOOST_MATH_IFGT_EVALUATOR_BASE_HPP_ER_2009
+#include <functional>
+#include <cmath>
+#include <vector>
+#include <algorithm>
+#include <boost/assert.hpp>
+#include <boost/bind.hpp>
+#include <boost/utility.hpp>
+#include <boost/iterator/range_left.hpp>
+#include <boost/iterator/range_right.hpp>
+#include <boost/iterator/range_same.hpp>
+#include <boost/math/ifgt/tag.hpp>
+namespace boost{namespace math{namespace ifgt{
+
+template<typename Common,typename Derived>
+class evaluator_base : public Common{
+ typedef Common common_type;
+
+ public:
+ typedef typename Common::value_type value_type;
+ static const std::size_t weight_size = Common::weight_size;
+
+ static const value_type default_nadaraya_watson_check_tolerance(){
+ static value_type val = 1e-15;
+ return val; };
+
+ template<typename ArgPack>
+ evaluator_base(const ArgPack& arg)
+ :common_type(arg),
+ nadaraya_watson_check_tolerance_(arg[
+ tag::nadaraya_watson_check_tolerance
+ |default_nadaraya_watson_check_tolerance()]){}
+
+ evaluator_base(const evaluator_base& that)
+ :common_type(that),nadaraya_watson_check_tolerance_(
+ that.nadaraya_watson_check_tolerance_){};
+
+ evaluator_base& operator=(const evaluator_base& that){
+ if(&that!=this){
+ common_type::operator=(that);
+ nadaraya_watson_check_tolerance_
+ = that.nadaraya_watson_check_tolerance_;
+ }
+ return *this;
+ }
+
+ /// \warning adds to range_out
+ template<typename R0,typename R1>
+ void gauss_transform(const R0& target,R1& range_out)const{
+ static range_same policy;
+ return this->gauss_transform(target,range_out,policy);
+ };
+
+ /// \warning adds to range_out
+ template<typename R0,typename R1,typename SubsetRangePolicy>
+ void gauss_transform(
+ const R0& target,
+ R1& range_out,
+ const SubsetRangePolicy& policy
+ )const{
+ const Derived& self = static_cast<const Derived&>(*this);
+ return self.call_impl(target,range_out,policy);
+ };
+
+ template<typename R0>
+ value_type rozenblatt_parzen_estimate(const R0& target)const{
+ typedef std::vector<value_type> range_type;
+ typedef range_left<1> policy_type;
+ static range_type range(1);
+ static policy_type policy;
+ range[0] = (value_type)(0);
+ this->gauss_transform(target,range,policy);
+ value_type count = (value_type)(this->sources_count());
+ range[0]/=count;
+ return range[0];
+ };
+
+ /// Requires weight[0]==1
+ /// Initializes range_out to zero and writes NW to it
+ template<typename R0,typename R1>
+ void nadaraya_watson_estimate(
+ const R0& target,R1& range_out)const{
+ typedef std::vector<value_type> range_type;
+ typedef range_right<weight_size-1> policy_type;
+ typedef typename result_of<policy_type(const range_type&)>::type
+ range_right_type;
+ static range_type range(weight_size);
+ static policy_type policy;
+ BOOST_ASSERT(size(range_out)-(weight_size-1)==0);
+ std::fill(
+ begin(range),
+ end(range),
+ (value_type)(0)
+ );
+ this->gauss_transform(target,range);
+ value_type denom = range[0];
+
+ range_right_type range_right = policy(range);
+ BOOST_ASSERT(size(range_right)-(weight_size-1)==0);
+ std::transform(
+ begin(range_right), end(range_right),
+ begin(range_out),
+ bind(std::divides<value_type>(), _1, denom)
+ );
+ };
+
+ /// Writes rozenblatt_parzen to range_out[0]
+ /// and NW to {range_out[j]:j=1,...,J-1}
+ /// Less practical but more efficient than each separately.
+ template<typename R0,typename R1>
+ void rozenblatt_parzen_and_nadaraya_watson_estimate(
+ const R0& target,R1& range_out)const{
+ typedef std::vector<value_type> range_type;
+ typedef range_right<weight_size-1> policy_type;
+ typedef typename range_iterator<R1>::type iter1_type;
+ typedef typename result_of<policy_type(const range_type&)>::type
+ range_right_type;
+ static range_type range(weight_size);
+ static policy_type policy;
+ std::fill(begin(range),end(range),(value_type)(0));
+ this->gauss_transform(target,range);
+ BOOST_ASSERT(size(range_out)-weight_size==0);
+ value_type denom = range[0];
+ value_type count = (value_type)(this->sources_count());;
+ iter1_type iter1 = begin(range_out);
+ (*iter1) = denom / count;
+ std::advance(iter1,1);
+ range_right_type range_right = policy(range);
+ std::transform(
+ begin(range_right), end(range_right), iter1,
+ bind(std::divides<value_type>(), _1, denom)
+ );
+ };
+
+ private:
+ evaluator_base();
+ value_type nadaraya_watson_check_tolerance_;
+
+};
+
+}}}
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/evaluator_common.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/evaluator_common.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,39 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/evaluator_common.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_EVALUATOR_COMMON_HPP_ER_2009
+#define BOOST_MATH_IFGT_EVALUATOR_COMMON_HPP_ER_2009
+#include <boost/math/ifgt/tag.hpp>
+namespace boost{namespace math{namespace ifgt{
+
+template<typename Accumulator>
+class evaluator_common{
+ public:
+ typedef Accumulator accumulator_type;
+ typedef typename accumulator_type::value_type value_type;
+ typedef typename accumulator_type::sources_count_type
+ sources_count_type;
+ typedef const accumulator_type& result_of_accumulator_type;
+
+ static const std::size_t source_size = accumulator_type::source_size;
+ static const std::size_t weight_size = accumulator_type::weight_size;
+
+ template<typename ArgPack>
+ evaluator_common(const ArgPack& arg):acc_(arg[tag::accumulator]){}
+
+ evaluator_common(const evaluator_common& that):acc_(that.acc_){}
+
+ sources_count_type sources_count()const{
+ return this->acc_.sources_count(); }
+ result_of_accumulator_type accumulator()const{ return this->acc_; }
+ private:
+ evaluator_common& operator=(const evaluator_common& that);
+ evaluator_common();
+ const accumulator_type& acc_;
+};
+
+}}}
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/exact_accumulator.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/exact_accumulator.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,141 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/exact_accumulator.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_EXACT_ACCUMULATOR_HPP_ER_2009
+#define BOOST_MATH_IFGT_EXACT_ACCUMULATOR_HPP_ER_2009
+#include <vector>
+#include <algorithm>
+#include <functional>
+#include <iterator>
+#include <vector>
+#include <cmath>
+#include <boost/range.hpp>
+#include <boost/static_assert.hpp>
+#include <boost/bind.hpp>
+#include <boost/iterator/zip_iterator.hpp>
+#include <boost/iterator/iterator_traits.hpp>
+#include <boost/tuple/tuple.hpp>
+#include <boost/math/ifgt/tag.hpp>
+#include <boost/math/ifgt/accumulator_base.hpp>
+#include <boost/math/ifgt/exact_contribution.hpp>
+namespace boost{namespace math{namespace ifgt{
+
+ template<
+ unsigned SourceSz,
+ unsigned WeightSize,
+ typename Cont = std::vector<double> >
+ class exact_accumulator
+ : public exact_accumulator_base<
+ exact_accumulator<SourceSz,WeightSize,Cont> >{
+ public:
+ typedef exact_contribution<SourceSz,WeightSize,Cont>
+ contribution_type;
+ typedef std::vector<contribution_type> contributions_type;
+ typedef typename contribution_type::value_type value_type;
+ typedef typename contributions_type::size_type size_type;
+ static const std::size_t source_size = SourceSz;
+ static const std::size_t weight_size = WeightSize;
+ typedef size_type sources_count_type;
+ typedef const contributions_type& cref_contributions_type;
+
+ exact_accumulator(value_type bandwidth);
+ template<typename ArgPack> exact_accumulator(const ArgPack& arg);
+ exact_accumulator(const exact_accumulator&);
+ exact_accumulator& operator=(const exact_accumulator&);
+
+ // Access
+ value_type active_bandwidth()const; //>models Accumulator
+ cref_contributions_type contributions()const;
+ sources_count_type sources_count()const;
+
+ // Update
+ void set(value_type bandwidth);
+
+ // Primarily intended for use internally
+ template<typename R0,typename R1>
+ void push_back_contribution(
+ const R0& source, const R1& weight,value_type bandwidth);
+
+ /// models Accumulator
+ template<typename R0,typename R1>
+ void operator()(const R0& source, const R1& weight);
+
+ private:
+ exact_accumulator();//dont implem
+ static size_type default_contributions_reserve_size(){
+ static size_type sz = 1e4; return sz; }
+ value_type bandwidth_;
+ contributions_type contributions_;
+ };
+ template<unsigned SourceSz, unsigned WeightSize, typename Cont>
+ exact_accumulator<SourceSz,WeightSize, Cont>::exact_accumulator(
+ value_type bandwidth):bandwidth_(bandwidth),contributions_(){
+ contributions_.reserve(default_contributions_reserve_size());
+ }
+
+ //TODO replace unsigned SourceSz by std::size_t SourceSz
+ template<unsigned SourceSz, unsigned WeightSize, typename Cont>
+ template<typename ArgPack>
+ exact_accumulator<SourceSz,WeightSize, Cont>::exact_accumulator(
+ const ArgPack& arg):bandwidth_(arg[tag::bandwidth]),
+ contributions_(){
+ contributions_.reserve(arg[tag::contributions_reserve_size|
+ default_contributions_reserve_size()]);
+ }
+
+ template<unsigned SourceSz, unsigned WeightSize, typename Cont>
+ exact_accumulator<SourceSz,WeightSize, Cont>::exact_accumulator(
+ const exact_accumulator& that):bandwidth_(that.bandwidth_),
+ contributions_(that.contributions_){}
+
+ template<unsigned SourceSz, unsigned WeightSize, typename Cont>
+ typename exact_accumulator<SourceSz,WeightSize, Cont>::exact_accumulator&
+ exact_accumulator<SourceSz,WeightSize, Cont>::operator=(
+ const exact_accumulator& that){
+ if(&that!=this){
+ bandwidth_ = that.bandwidth_;
+ contributions_ = that.contributions_; }
+ return *this; }
+
+ template<unsigned SourceSz, unsigned WeightSize, typename Cont>
+ typename exact_accumulator<SourceSz,WeightSize, Cont>::value_type
+ exact_accumulator<SourceSz,WeightSize, Cont>::active_bandwidth()const{
+ return bandwidth_; };
+
+ template<unsigned SourceSz, unsigned WeightSize, typename Cont>
+ typename exact_accumulator<SourceSz,WeightSize, Cont>::cref_contributions_type
+ exact_accumulator<SourceSz,WeightSize, Cont>::contributions()const{
+ return contributions_; };
+
+ template<unsigned SourceSz, unsigned WeightSize, typename Cont>
+ typename exact_accumulator<SourceSz,WeightSize, Cont>::sources_count_type
+ exact_accumulator<SourceSz,WeightSize, Cont>::sources_count()const{
+ return contributions_.size();
+ }
+
+ template<unsigned SourceSz, unsigned WeightSize, typename Cont>
+ void exact_accumulator<SourceSz,WeightSize, Cont>::set(
+ value_type bandwidth){ bandwidth_ = bandwidth; };
+
+ template<unsigned SourceSz, unsigned WeightSize, typename Cont>
+ template<typename R0,typename R1>
+ void exact_accumulator<SourceSz,WeightSize, Cont>::push_back_contribution(
+ const R0& source,
+ const R1& weight,
+ value_type bandwidth){
+ contributions_.push_back(
+ contribution_type(source,weight,bandwidth)
+ );
+ }
+
+ template<unsigned SourceSz,unsigned WeightSize, typename Cont>
+ template<typename R0,typename R1>
+ void exact_accumulator<SourceSz,WeightSize, Cont>::operator()(
+ const R0& source, const R1& weight){
+ push_back_contribution(source,weight,active_bandwidth()); }
+
+}}}
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/exact_contribution.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/exact_contribution.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,156 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/exact_contribution.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_EXACT_CONTRIBUTION_HPP_ER_2009
+#define BOOST_MATH_IFGT_EXACT_CONTRIBUTION_HPP_ER_2009
+#include <algorithm>
+#include <functional>
+#include <numeric>
+#include <cmath>
+#include <stdexcept>
+#include <boost/range.hpp>
+#include <boost/assert.hpp>
+#include <boost/utility.hpp>
+#include <boost/algorithm/l2_distance_squared.hpp>
+#include <boost/math/ifgt/normal_kernel_properties.hpp>
+#include <boost/iterator/range_same.hpp>
+
+namespace boost{namespace math{namespace ifgt{
+
+ // A source, weight, and bandwidth
+ // and associated operations.
+ template<
+ unsigned SourceSz,
+ unsigned WeightSz,
+ typename Cont = std::vector<double> > //TODO add cutoff radius policy
+ class exact_contribution{
+ public:
+ typedef typename Cont::value_type value_type;
+ typedef Cont source_type;
+ typedef Cont weight_type;
+ static const std::size_t source_size = SourceSz;
+ static const std::size_t weight_size = WeightSz;
+ template<typename R0,typename R1>
+ exact_contribution(
+ const R0& source,
+ const R1& weight,
+ value_type bandwidth
+ ):
+ source_(source_size),
+ weight_(weight_size),
+ bandwidth_(bandwidth)
+ {
+ BOOST_ASSERT(size(source_)-source_size==0);
+ BOOST_ASSERT(size(weight_)-weight_size==0);
+
+ copy(
+ begin(source),
+ end(source),
+ begin(source_)
+ );
+
+ copy(
+ begin(weight),
+ end(weight),
+ begin(weight_)
+ );
+ }
+
+ exact_contribution(const exact_contribution& that)
+ : source_(that.source_),
+ weight_(that.weight_),
+ bandwidth_(that.bandwidth_){}
+
+ exact_contribution& operator=(
+ const exact_contribution& that){
+ if(&that!=this){
+ source_ = that.source_;
+ weight_ = that.weight_;
+ bandwidth_ = that.bandwidth_;
+ }
+ return *this;
+ }
+
+ const source_type& source()const{ return source_; }
+ const source_type& weight()const{ return weight_; }
+ value_type bandwidth()const{ return bandwidth_; }
+ value_type normalizing_constant()const{
+ static value_type sqrt2 = sqrt((value_type)(2));
+ typedef normal_kernel_properties<
+ source_size,value_type> nkp_type;
+ static value_type result
+ = nkp_type::normalizing_constant(bandwidth()/sqrt2);
+ return result;
+ }
+ template<typename R>
+ value_type unnormalized_gauss_kernel(const R& target)const{
+ BOOST_ASSERT(size(target)-source_size==0);
+ typedef l2_distance_squared::type dist_type;
+ static dist_type dist;
+ static value_type squared_bandwidth
+ = (bandwidth() * bandwidth());
+ value_type zscore_sqnorm
+ = ( dist(source_,target) / squared_bandwidth );
+ return exp(-zscore_sqnorm);
+ };
+
+ template<typename R>
+ value_type gauss_kernel(const R& target)const{
+ return (unnormalized_gauss_kernel(target)/normalizing_constant());
+ };
+
+ template<typename R0,typename R1>
+ void evaluate(const R0& target,R1& range_out)const{
+ static range_same policy;
+ return evaluate(target,range_out,policy);
+ }
+ /// \warning *adds* to range_out
+ template<typename R0,typename R1,typename SubsetRangePolicy>
+ void evaluate(const R0& target,R1& range_out,
+ const SubsetRangePolicy& policy)const{
+ value_type gk = gauss_kernel(target);
+
+ typedef const weight_type& arg_type;
+ typedef typename result_of<SubsetRangePolicy(arg_type)>::type
+ subset_type;
+ subset_type subset = policy(weight_);
+
+ BOOST_ASSERT(
+ size(subset)==size(range_out)
+ );
+
+ transform(
+ begin(subset),
+ end(subset),
+ begin(range_out),
+ begin(range_out),
+ bind(
+ std::plus<value_type>(),
+ bind(
+ std::multiplies<value_type>(),
+ _1,
+ gk
+ ),
+ _2
+ )
+ );
+ }
+
+ template<typename R0>
+ value_type
+ nadaraya_watson_normalizing_constant(const R0& target)const{
+ return gauss_kernel(target);
+ }
+
+ private:
+ exact_contribution();//dont implement
+ source_type source_;
+ weight_type weight_;
+ value_type bandwidth_;
+ };
+
+}}}
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/exact_contributions_evaluator.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/exact_contributions_evaluator.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,76 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/exact_contributions_evaluator.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_EXACT_CONTRIBUTIONS_EVALUATOR_HPP_ER_2009
+#define BOOST_MATH_IFGT_EXACT_CONTRIBUTIONS_EVALUATOR_HPP_ER_2009
+#include <algorithm>
+#include <boost/range.hpp>
+#include <boost/ref.hpp>
+#include <boost/range/value_type.hpp>
+#include <boost/bind.hpp>
+#include <boost/ref.hpp>
+#include <boost/math/ifgt/exact_contribution.hpp>
+
+namespace boost{namespace math{namespace ifgt{
+
+ //\warning adds to range_out
+ template<typename R0>
+ class exact_contributions_evaluator{
+ typedef R0 contributions_type;
+ typedef typename range_value<R0>::type contribution_type;
+ public:
+ exact_contributions_evaluator(
+ const contributions_type& contributions
+ ):contributions_(contributions){}
+
+ exact_contributions_evaluator(
+ const exact_contributions_evaluator& that)
+ :contributions_(that.contributions_){}
+
+ template<typename R1,typename R2>
+ void operator()(const R1& target,R2& range_out)const{
+ typedef range_same policy_type;
+ static policy_type policy;
+ return this->operator()(target,range_out,policy);
+ }
+
+ // \warning *adds* to range_out
+ template<typename R1,typename R2,typename SubsetRangePolicy>
+ void operator()(const R1& target,R2& range_out,
+ const SubsetRangePolicy& policy)const{
+ void (contribution_type::*f) (
+ const R1&,R2&,const SubsetRangePolicy&) const
+ = &contribution_type::evaluate;
+ for_each(
+ begin(contributions_),
+ end(contributions_),
+ bind(
+ f,
+ _1,
+ cref(target),
+ ref(range_out),
+ cref(policy)
+ )
+ );
+ }
+
+ private:
+ exact_contributions_evaluator&
+ operator=(const exact_contributions_evaluator& that);
+ exact_contributions_evaluator();
+ const contributions_type& contributions_;
+ };
+
+ template<typename R0>
+ exact_contributions_evaluator<R0>
+ make_exact_contributions_evaluator(
+ const R0& contributions
+ ){
+ return exact_contributions_evaluator<R0>(contributions);
+ }
+
+}}}
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/exact_evaluator.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/exact_evaluator.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,47 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/exact_evaluator.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_EXACT_EVALUATOR_HPP_ER_2009
+#define BOOST_MATH_IFGT_EXACT_EVALUATOR_HPP_ER_2009
+#include <boost/math/ifgt/evaluator_common.hpp>
+#include <boost/math/ifgt/evaluator_base.hpp>
+#include <boost/math/ifgt/exact_contributions_evaluator.hpp>
+namespace boost{namespace math{namespace ifgt{
+
+template<typename ExactAccumulator>
+class exact_evaluator : public evaluator_base<
+ evaluator_common<ExactAccumulator>,
+ exact_evaluator<ExactAccumulator>
+ >{
+ typedef exact_evaluator<ExactAccumulator> self_type;
+ typedef evaluator_common<ExactAccumulator> c_type;
+ typedef evaluator_base<c_type,self_type> base_type;
+ //typedef typename base_type::value_type value_type;
+ public:
+
+ template<typename ArgPack>
+ exact_evaluator(const ArgPack& arg) :
+ base_type(arg){}
+
+ exact_evaluator(const exact_evaluator& that) :
+ base_type(that){}
+
+ // TODO privatize
+ template<typename R0,typename R1,typename SubsetRangePolicy>
+ void call_impl(const R0& target,R1& range_out,
+ const SubsetRangePolicy& policy)const{
+ make_exact_contributions_evaluator(
+ (this->accumulator()).contributions()
+ )(target, range_out,policy);
+ }
+
+ private:
+ exact_evaluator& operator=(const exact_evaluator& that);
+ exact_evaluator();
+};
+
+}}}
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/fast_accumulator.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/fast_accumulator.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,258 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/fast_accumulator.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_FAST_ACCUMULATOR_IMPL_HPP_ER_2009
+#define BOOST_MATH_IFGT_FAST_ACCUMULATOR_IMPL_HPP_ER_2009
+#include <algorithm>
+#include <vector>
+#include <functional>
+#include <numeric>
+#include <boost/range.hpp>
+#include <boost/bind.hpp>
+#include <boost/mpl/placeholders.hpp>
+#include <boost/mpl/apply.hpp>
+#include <iostream> //TODO remove
+#include <boost/math/ifgt/tag.hpp>
+#include <boost/math/ifgt/find_nearest_cluster.hpp>
+#include <boost/math/ifgt/truncation_degree_constant.hpp>
+#include <boost/math/ifgt/cutoff_radius_none.hpp>
+#include <boost/math/ifgt/accumulator_base.hpp>
+namespace boost{namespace math{namespace ifgt{
+
+
+ // fast_accumulator maintains a set of clusters, which is split
+ // into an inactive and an active range. each time a (source,weight)
+ // is passed to the object, only the active range is updated
+ // i.e. either the source is collected at the nearest cluster within
+ // that range,
+ // or the a new cluster is added to the active range,
+ // or the source is ignored if it equals the center of nearest cluster.
+ // each time (bandwith,max_radius) is reset, the beginning
+ // of the active range is to the end of the clusters range.
+ // This (optional) feature is useful in applications where the bandwidth
+ // is modified (usually decreased) in a piecewise-constant fashion.
+ // The client can also use push_back_cluster to add a center
+ // to the active range
+
+ /// apply<FindClusterPolicy,Cluster>::type
+ template<
+ typename Cluster,
+ typename FindClusterPolicy
+ = find_nearest_cluster<mpl::_1, l2_distance_squared>
+ >
+ class fast_accumulator:
+ public fast_accumulator_base<
+ fast_accumulator<Cluster,FindClusterPolicy> >{
+ public:
+ typedef Cluster cluster_type;
+ typedef typename cluster_type::center_type center_type;
+ typedef typename cluster_type::value_type value_type;
+ typedef typename cluster_type::coefficients_type
+ coefficients_type;
+ typedef typename cluster_type::collection_coefficients_type
+ collection_coefficients_type;
+ static const std::size_t source_size = cluster_type::source_size;
+ static const std::size_t weight_size = cluster_type::weight_size;
+ private:
+ typedef std::vector<cluster_type> clusters_type;
+ typedef typename clusters_type::size_type size_type;
+ static size_type default_clusters_reserve_size(){
+ static size_type sz = 1000; return sz; }
+ typedef typename clusters_type::difference_type diff_type;
+ typedef range_iterator<clusters_type> ri_type;
+ typedef typename ri_type::type iter_type;
+ typedef iterator_range<iter_type> active_range_type;
+ typedef typename cluster_type::truncation_degree_policy_type
+ truncation_degree_policy_type;
+ public:
+ typedef typename cluster_type::sources_count_type sources_count_type;
+
+ //Construct
+ fast_accumulator();
+ template<typename ArgPack>
+ fast_accumulator(const ArgPack& args);
+ fast_accumulator(const fast_accumulator& that);
+ //Assign
+ fast_accumulator& operator=(const fast_accumulator& that);
+
+ /// Triggers a new active range
+ void set(value_type bandwidth,value_type max_radius);
+ void set_truncation_degree_policy(
+ const truncation_degree_policy_type& arg);
+
+ //Access
+ value_type active_bandwidth()const;
+ value_type active_max_radius()const;
+ const clusters_type& clusters()const;
+ sources_count_type sources_count()const;
+
+ //Update
+ /// Used internally, but passing clusters directly also allowed.
+ template<typename R>
+ void push_back_cluster(const R& center);
+
+ template<typename R0,typename R1>
+ bool operator()(const R0& source, const R1& weight);
+
+ private:
+ active_range_type active_range();
+ diff_type active_dist_to_;
+ value_type active_bandwidth_;
+ value_type active_max_radius_;
+ truncation_degree_policy_type truncation_degree_policy_;
+ clusters_type clusters_;
+ };
+
+ template<typename Cluster, typename FindClusterPolicy>
+ fast_accumulator<Cluster,FindClusterPolicy>::fast_accumulator()
+ : active_dist_to_((diff_type)(0)),
+ active_bandwidth_((value_type)(0)),
+ active_max_radius_((value_type)(0)),
+ //make sure gen comp error if deflt constr inexistant
+ truncation_degree_policy_(),
+ clusters_()
+ { clusters_.reserve(default_clusters_reserve_size()); }
+
+ template<typename Cluster,typename FindClusterPolicy>
+ template<typename ArgPack>
+ fast_accumulator<Cluster,FindClusterPolicy>::fast_accumulator(
+ const ArgPack& args)
+ : active_dist_to_((diff_type)(0)),
+ active_bandwidth_(args[tag::bandwidth|(value_type)(0)]),
+ active_max_radius_(args[tag::max_cluster_radius|(value_type)(0)]),
+ truncation_degree_policy_(args){
+ clusters_.reserve(
+ args[tag::clusters_reserve_size
+ |default_clusters_reserve_size()]);
+ }
+
+ template<typename Cluster,typename FindClusterPolicy>
+ fast_accumulator<Cluster,FindClusterPolicy>::fast_accumulator(
+ const fast_accumulator& that)
+ : active_dist_to_(that.active_dist_to_),
+ active_bandwidth_(that.active_bandwidth_),
+ active_max_radius_(that.active_max_radius_),
+ clusters_(that.clusters_){}
+
+ template<typename Cluster,typename FindClusterPolicy>
+ typename fast_accumulator<Cluster,
+ FindClusterPolicy>::fast_accumulator&
+ fast_accumulator<Cluster,FindClusterPolicy>::operator=(
+ const fast_accumulator& that){
+ if(&that!=this){
+ active_dist_to_ = that.active_dist_to_;
+ active_bandwidth_ = that.active_bandwidth_;
+ active_max_radius_ = that.active_max_radius_;
+ clusters_ = that.clusters_;
+ }
+ return *this;
+ }
+
+ template<typename Cluster,typename FindClusterPolicy>
+ void fast_accumulator<Cluster,FindClusterPolicy>::set(
+ value_type bandwidth,value_type max_radius){
+ active_bandwidth_ = bandwidth;
+ active_max_radius_ = max_radius;
+ active_dist_to_
+ = std::distance(begin(clusters_), end(clusters_));
+ }
+
+ template<typename Cluster,typename FindClusterPolicy>
+ void
+ fast_accumulator<Cluster,
+ FindClusterPolicy>::set_truncation_degree_policy(
+ const truncation_degree_policy_type& truncation_degree_policy){
+ truncation_degree_policy_ = truncation_degree_policy;
+ }
+
+ template<typename Cluster,typename FindClusterPolicy>
+ typename fast_accumulator<Cluster,FindClusterPolicy>::value_type
+ fast_accumulator<Cluster,FindClusterPolicy>::active_bandwidth()const{
+ return active_bandwidth_; }
+
+ template<typename Cluster,typename FindClusterPolicy>
+ typename fast_accumulator<Cluster,FindClusterPolicy>::value_type
+ fast_accumulator<Cluster,FindClusterPolicy>::active_max_radius()const{
+ return active_max_radius_; }
+
+ template<typename Cluster,typename FindClusterPolicy>
+ typename fast_accumulator<Cluster,
+ FindClusterPolicy>::active_range_type
+ fast_accumulator<Cluster,FindClusterPolicy>::active_range(){
+ iter_type i = begin(clusters_);
+ std::advance(i,active_dist_to_);
+ return active_range_type(i,end(clusters_));
+ }
+
+ template<typename Cluster, typename FindClusterPolicy>
+ const typename fast_accumulator<Cluster,
+ FindClusterPolicy>::clusters_type&
+ fast_accumulator<Cluster,FindClusterPolicy>::clusters()const{
+ return clusters_; }
+
+ template<typename Cluster,typename FindClusterPolicy>
+ typename fast_accumulator<Cluster,
+ FindClusterPolicy>::sources_count_type
+ fast_accumulator<Cluster,FindClusterPolicy>::sources_count()const{
+ return std::accumulate(
+ begin(clusters_),
+ end(clusters_),
+ (sources_count_type)(0),
+ bind(
+ std::plus<value_type>(),
+ _1,
+ bind(
+ &cluster_type::sources_count,
+ _2
+ )
+ )
+ );
+ }
+
+ template<typename Cluster, typename FindClusterPolicy>
+ template<typename R>
+ void fast_accumulator<Cluster,FindClusterPolicy
+ >::push_back_cluster(const R& center){
+ clusters_.push_back(cluster_type(
+ center,
+ active_bandwidth(),
+ active_max_radius(),
+ truncation_degree_policy_
+ ));
+ }
+
+ template<typename Cluster, typename FindClusterPolicy>
+ template<typename R0,typename R1>
+ bool fast_accumulator<Cluster,FindClusterPolicy>::operator()(
+ const R0& source, const R1& weight){
+ static bool is_collected = false;
+ static active_range_type local_active_range = active_range();
+ if(active_bandwidth()>(value_type)(0)){
+ local_active_range = active_range();
+ if(empty(local_active_range)){
+ push_back_cluster(source);
+ }else{
+ typedef typename mpl::apply<FindClusterPolicy,
+ active_range_type>::type find_nc_type;
+
+ // TODO assigning local_active_range
+ // and creating a new find_nc on each call maybe costly
+ find_nc_type find_nc(local_active_range);
+ iter_type nn = find_nc(source);
+ is_collected = (*nn)(source,weight);
+ if(nn->source_radius_exceeds_max_radius()){
+ BOOST_ASSERT(!is_collected);
+ push_back_cluster(source);
+ }
+ }
+ }else{
+ throw std::runtime_error("fast_accumulator::collect");
+ }
+ return is_collected;
+ }
+
+}}}
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/fast_evaluator.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/fast_evaluator.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,58 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/fast_evaluator.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_FAST_EVALUATOR_HPP_ER_2009
+#define BOOST_MATH_IFGT_FAST_EVALUATOR_HPP_ER_2009
+#include <stdexcept>
+#include <boost/math/ifgt/cutoff_radius_none.hpp>
+#include <boost/math/ifgt/clusters_evaluator.hpp>
+#include <boost/math/ifgt/evaluator_common.hpp>
+#include <boost/math/ifgt/evaluator_base.hpp>
+namespace boost{namespace math{namespace ifgt{
+
+ template<
+ typename FastAccumulator,
+ typename CutoffRadiusPolicy = cutoff_radius_none<mpl::_1>
+ >
+ class fast_evaluator : public
+ evaluator_base<
+ evaluator_common<FastAccumulator>,
+ fast_evaluator<
+ FastAccumulator,
+ CutoffRadiusPolicy>
+ >{
+
+ typedef fast_evaluator<FastAccumulator,CutoffRadiusPolicy> self_type;
+ typedef evaluator_common<FastAccumulator> c_type;
+ typedef evaluator_base<c_type,self_type> base_type;
+// typedef typename base_type::value_type value_type;
+ typedef clusters_evaluator<typename fast_evaluator::value_type,
+ CutoffRadiusPolicy> eval_type;
+
+ public:
+
+ template<typename ArgPack>
+ fast_evaluator(const ArgPack& arg) : base_type(arg),eval_(arg){}
+
+ fast_evaluator(const fast_evaluator& that)
+ :base_type(that),eval_(that.eval_){}
+
+ //TODO find a way to privatize
+ template<typename R0,typename R1,typename SubsetRangePolicy>
+ void call_impl(const R0& target, R1& range_out,
+ const SubsetRangePolicy& policy)const{
+ return eval_(
+ target,this->accumulator().clusters(),range_out,policy);
+ }
+
+ private:
+ fast_evaluator& operator=(const fast_evaluator& that);
+ fast_evaluator();
+ eval_type eval_;
+ };
+
+}}}
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/find_nearest_cluster.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/find_nearest_cluster.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,53 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/find_nearest_cluster.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_FIND_NEAREST_CLUSTER_HPP_ER_2009
+#define BOOST_MATH_IFGT_FIND_NEAREST_CLUSTER_HPP_ER_2009
+#include <cmath>
+#include <stdexcept>
+#include <boost/function.hpp>
+#include <boost/range/value_type.hpp>
+#include <boost/mpl/assert.hpp>
+#include <boost/algorithm/find_nearest_neighbor.hpp>
+#include <boost/iterator/transform_iterator.hpp>
+#include <boost/math/ifgt/cluster.hpp>
+#include <boost/math/ifgt/cluster_call_center.hpp>
+
+namespace boost{namespace math{namespace ifgt{
+ template<typename R, typename Distance = l2_distance_squared>
+ class find_nearest_cluster{
+ typedef range_iterator<R> ri_type;
+ public:
+ typedef typename ri_type::type iter_type;
+ typedef typename range_value<R>::type cluster_type;
+
+ find_nearest_cluster(R& clusters):clusters_(clusters){}
+ find_nearest_cluster(const find_nearest_cluster& that)
+ :clusters_(that.clusters_){}
+ template<typename R0>
+ iter_type operator()(const R0& source){
+ typedef cluster_call_center<cluster_type> f_type;
+ typedef transform_iterator<f_type,iter_type>
+ iter_t_type;
+ typedef iterator_range<iter_t_type> ir_type;
+ f_type f;
+ ir_type ir(
+ iter_t_type(begin(clusters_),f),
+ iter_t_type(end(clusters_),f)
+ );
+
+ typedef find_nearest_neighbor<ir_type,Distance> find_type;
+
+ iter_type nn = find_type(ir)(source).base();
+ return nn;
+ }
+ private:
+ find_nearest_cluster& operator=(const find_nearest_cluster& that);
+ R& clusters_;
+ };
+
+}}}
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/for_each_accumulate.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/for_each_accumulate.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,28 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/for_each_accumulate.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_FOR_EACH_ACCUMULATE_HPP_ER_2009
+#define BOOST_MATH_IFGT_FOR_EACH_ACCUMULATE_HPP_ER_2009
+#include <boost/iterator/for_each_over_2_vector2matrix.hpp>
+#include <boost/math/ifgt/traits.hpp>
+namespace boost{namespace math{namespace ifgt{
+
+ template<typename R0, typename R1, typename Acc>
+ void for_each(
+ const R0& sources,
+ const R1& weights,
+ accumulator_base<Acc>& acc
+ ){
+ typedef Acc derived_type;
+ typedef traits<derived_type> the_traits;
+ typedef for_each_over_2_vector2matrix<the_traits> for_each_type;
+ derived_type& derived = static_cast<derived_type&>(acc);
+ for_each_type f(derived);
+ f(sources,weights);
+ };
+
+}}}
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/for_each_accumulate_rest_weights.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/for_each_accumulate_rest_weights.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,36 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/for_each_accumulate_rest_weights.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_FOR_EACH_ACCUMULATE_REST_WEIGHTS_HPP_ER_2009
+#define BOOST_MATH_IFGT_FOR_EACH_ACCUMULATE_REST_WEIGHTS_HPP_ER_2009
+#include <boost/iterator/for_each_over_2_vector2matrix.hpp>
+#include <boost/iterator/insert_element_every_n_step_iterator.hpp>
+#include <boost/math/ifgt/for_each_accumulate.hpp>
+#include <boost/math/ifgt/rest_weights_wrapper.hpp>
+
+namespace boost{namespace math{namespace ifgt{
+
+ template<typename R0, typename R1, typename Acc>
+ void for_each(
+ const R0& sources,
+ const rest_weights_wrapper<R1>& rw,
+ accumulator_base<Acc>& acc
+ ){
+
+ typedef typename result_of_make_range_insert_element_every_n_step<
+ const R1>::type weights_type;
+ std::size_t n = Acc::weight_size - 1;
+ BOOST_ASSERT(n>0);
+ weights_type weights = make_range_insert_element_every_n_step(
+ rw(),
+ 1.0,
+ n
+ );
+ for_each(sources,weights,acc);
+ }
+
+}}}
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/for_each_evaluate.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/for_each_evaluate.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,34 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/for_each_evaluate.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_IFGT_FOR_EACH_EVALUATE_HPP_ER_2009
+#define BOOST_IFGT_FOR_EACH_EVALUATE_HPP_ER_2009
+#include <algorithm>
+#include <boost/ref.hpp>
+#include <boost/iterator/for_each_over_2_vector2matrix.hpp>
+#include <boost/math/ifgt/traits.hpp>
+#include <boost/math/ifgt/evaluator_base.hpp>
+#include <boost/math/ifgt/call_wrapper.hpp>
+namespace boost{namespace math{namespace ifgt{
+
+template<
+ typename R0, typename R1,
+ typename Common, typename Derived,template<typename> class Method>
+void for_each(
+ const R0& targets,
+ R1& ranges_out,
+ const evaluator_base<Common,Derived>& e,
+ const call<Method>& method
+){
+ typedef Method<Derived> call_type;
+ typedef ifgt::call_traits<call_type> traits;
+ typedef for_each_over_2_vector2matrix<traits> for_each_type;
+ call_type call(e);
+ for_each_type f( call );
+ f(targets,ranges_out);
+}
+}}}
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/for_each_rozenblatt_parzen_estimate.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/for_each_rozenblatt_parzen_estimate.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,34 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/for_each_rozenblatt_parzen_estimate.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_FOR_EACH_ROZENBLATT_PARZEN_ESTIMATE_HPP_ER_2009
+#define BOOST_MATH_IFGT_FOR_EACH_ROZENBLATT_PARZEN_ESTIMATE_HPP_ER_2009
+#include <algorithm>
+#include <boost/bind.hpp>
+#include <boost/ref.hpp>
+#include <boost/iterator/vector2matrix_iterator.hpp>
+#include <boost/math/ifgt/call_rozenblatt_parzen_estimate.hpp>
+#include <boost/math/ifgt/call_wrapper.hpp>
+#include <boost/math/ifgt/evaluator_base.hpp>
+namespace boost{namespace math{namespace ifgt{
+template<typename R0, typename R1, typename Common,typename Derived>
+void for_each(
+ const R0& targets,
+ R1& ranges_out,
+ const evaluator_base<Common,Derived>& e,
+ const call<rozenblatt_parzen_estimate>& method
+){
+ const std::size_t stride = Derived::source_size;
+ transform(
+ make_vector2matrix_iterator(begin(targets),stride),
+ make_end_vector2matrix_iterator(begin(targets),end(targets),
+ stride),
+ begin(ranges_out),
+ make_rozenblatt_parzen_estimate(e)
+ );
+}
+}}}
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/include.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/include.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,28 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/include.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_INCLUDE_HPP_ER_2009
+#define BOOST_MATH_IFGT_INCLUDE_HPP_ER_2009
+#include <boost/math/ifgt/call_gauss_transform.hpp>
+#include <boost/math/ifgt/call_nadaraya_watson_estimate.hpp>
+#include <boost/math/ifgt/call_rozenblatt_parzen_estimate.hpp>
+#include <boost/math/ifgt/call_wrapper.hpp>
+#include <boost/math/ifgt/cutoff_radius_none.hpp>
+#include <boost/math/ifgt/cutoff_radius_rydg.hpp>
+#include <boost/math/ifgt/exact_accumulator.hpp>
+#include <boost/math/ifgt/exact_evaluator.hpp>
+#include <boost/math/ifgt/fast_accumulator.hpp>
+#include <boost/math/ifgt/fast_evaluator.hpp>
+#include <boost/math/ifgt/find_nearest_cluster.hpp>
+#include <boost/math/ifgt/for_each_accumulate.hpp>
+#include <boost/math/ifgt/for_each_accumulate_rest_weights.hpp>
+#include <boost/math/ifgt/for_each_evaluate.hpp>
+#include <boost/math/ifgt/for_each_rozenblatt_parzen_estimate.hpp>
+#include <boost/math/ifgt/optimal_bandwidth.hpp>
+#include <boost/math/ifgt/optimal_parameter_given_max_degree.hpp>
+#include <boost/math/ifgt/tag.hpp>
+#include <boost/math/ifgt/truncation_degree_constant.hpp>
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/normal_kernel_properties.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/normal_kernel_properties.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,39 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/normal_kernel_properties.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_NORMAL_KERNEL_PROPERTIES_HPP_ER_2009
+#define BOOST_MATH_IFGT_NORMAL_KERNEL_PROPERTIES_HPP_ER_2009
+#include <cmath>
+namespace boost{namespace math{namespace ifgt{
+ // Normal kernel : exp(-(1/2)||(x-y)/sigma||^2) (2*pi*sigma^2)^(-d/2)
+ // Note that the bandwidth is h = sqrt(2) * sigma
+ template<unsigned SourceSz,typename RealType = double>
+ struct normal_kernel_properties{
+ typedef RealType value_type;
+ static const std::size_t source_size = SourceSz;
+ static value_type
+ pi(){
+ static value_type
+ val = 3.14159265358979323846264338327950288419716939937510;
+ return val;
+ }
+ static value_type normalizing_constant(
+ value_type sigma
+ ){
+ static value_type d = (value_type)(SourceSz);
+ static value_type sqrt2pi = sqrt((value_type)(2)*pi());
+ return std::pow(sigma*sqrt2pi,d);
+ }
+ static value_type normalizing_constant(
+ value_type sigma,
+ std::size_t n
+ ){
+ return (value_type)(n) * normalizing_constant(sigma);
+ }
+ };
+
+}}}
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/optimal_bandwidth.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/optimal_bandwidth.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,48 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/optimal_bandwidth.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_OPTIMAL_BANDWIDTH_FOR_DENSITY_ESTIMATION_UNDER_NORMALITY_HPP_ER_2009
+#define BOOST_MATH_IFGT_OPTIMAL_BANDWIDTH_FOR_DENSITY_ESTIMATION_UNDER_NORMALITY_HPP_ER_2009
+#include <cmath>
+namespace boost{namespace math{namespace ifgt{
+// \warning the sart(2)*sigma adjustmnent is needed so that
+// the kernel is exp(-||(x-y)/h||^2) = exp(-(1/2)||(x-y)/sigma||^2)
+template<std::size_t Dimension, typename RealType = double>
+struct optimal_bandwidth{
+ typedef RealType value_type;
+ static value_type for_density_estimation_under_normality(
+ value_type sigma,
+ std::size_t count
+ );
+ ///Assumes sigma = 1;
+ static value_type for_density_estimation_under_normality(
+ std::size_t count);
+
+};
+
+template<std::size_t Dimension, typename RealType>
+typename optimal_bandwidth<Dimension,RealType>::value_type
+optimal_bandwidth<Dimension,RealType>::for_density_estimation_under_normality(
+ value_type sigma, std::size_t count){
+ static const value_type sqrt2 = sqrt(2.0);
+ static const value_type d = (value_type)(Dimension);
+ static const value_type e = -(value_type)(1)/((value_type)(d+4));
+ value_type n = (value_type)(count);
+ value_type a = (n * (value_type)(2+d))/ ((value_type)(4));
+ return sigma * std::pow( a, e ) * sqrt2;
+}
+
+template<std::size_t Dimension, typename RealType>
+typename optimal_bandwidth<Dimension,RealType>::value_type
+optimal_bandwidth<Dimension,RealType>::for_density_estimation_under_normality(
+std::size_t count){
+ return for_density_estimation_under_normality(
+ (value_type)(1),
+ count);
+}
+
+}}}
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/optimal_parameter_given_max_degree.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/optimal_parameter_given_max_degree.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,167 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/optimal_parameter_given_max_degree.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_OPTIMAL_PARAMETER_GIVEN_MAX_DEGREE_HPP_ER_2009
+#define BOOST_MATH_IFGT_OPTIMAL_PARAMETER_GIVEN_MAX_DEGREE_HPP_ER_2009
+#include <stdexcept>
+#include <boost/math/ifgt/detail/divider.hpp>
+#include <boost/math/ifgt/cutoff_radius_rydg.hpp>
+#include <boost/math/ifgt/truncation_properties.hpp>
+namespace boost{ namespace math{ namespace ifgt{
+
+/// Finds optimal max cluster radius recursively,
+/// given a cutoff radius policy and a maximum degree p.
+template<
+ typename RealType,
+ template<typename> class CutoffRadiusPolicy = cutoff_radius_rydg, //TODO lambda
+ template<typename> class UpdatePolicy = divider //TODO lambda
+>
+class optimal_parameter_given_max_degree{
+ typedef CutoffRadiusPolicy<RealType> cutoff_radius_policy_type;
+ typedef UpdatePolicy<RealType> update_policy_type;
+ public:
+ typedef RealType value_type;
+ typedef std::size_t size_type;
+ template<typename ArgPack> optimal_parameter_given_max_degree(
+ const ArgPack& arg);
+ optimal_parameter_given_max_degree(
+ const optimal_parameter_given_max_degree& that);
+ optimal_parameter_given_max_degree&
+ operator=(const optimal_parameter_given_max_degree& that);
+
+ value_type error_tolerance()const{return error_tol_; };
+ unsigned max_degree()const{ return max_degree_; };
+ value_type bandwidth()const{return bandwidth_; };//models ClusterPolicy
+ size_type max_recursions()const{ return max_recursions_; }
+
+ /// Set parameters
+ void operator()(
+ unsigned max_degree,
+ value_type bandwidth,
+ value_type start_cluster_radius_over_bandwidth = 2.0,
+ size_type max_recursions = 20);
+
+ /// Output
+ value_type cutoff_radius()const{ return cutoff_radius_; }
+ //models ClusterPolicy
+ value_type max_radius()const{ return max_cluster_radius_; }
+ value_type error()const { return error_; }
+ size_type recursions()const{ return recursions_; }
+
+ private:
+ void call_impl();
+ value_type error_tol_;
+ cutoff_radius_policy_type cutoff_radius_policy_;
+ update_policy_type update_policy_;
+ value_type bandwidth_;
+ unsigned max_degree_;
+ size_type max_recursions_;
+ value_type cutoff_radius_;
+ value_type max_cluster_radius_;
+ value_type error_;
+ size_type recursions_;
+};
+
+template<typename RealType,template<typename> class CutoffRadiusPolicy,
+template<typename> class UpdatePolicy>
+template<typename ArgPack>
+optimal_parameter_given_max_degree<RealType,CutoffRadiusPolicy,
+UpdatePolicy>::optimal_parameter_given_max_degree(
+ const ArgPack& arg)
+:error_tol_(arg[tag::error_tolerance]),
+cutoff_radius_policy_(arg),
+update_policy_(arg),
+bandwidth_(0.0),
+max_degree_(0.0),
+max_recursions_(0),
+cutoff_radius_(0.0),
+max_cluster_radius_(0.0),
+error_(0.0),
+recursions_(0){}
+
+template<typename RealType,template<typename> class CutoffRadiusPolicy,
+template<typename> class UpdatePolicy>
+optimal_parameter_given_max_degree<RealType,CutoffRadiusPolicy,
+UpdatePolicy>::optimal_parameter_given_max_degree(
+const optimal_parameter_given_max_degree& that)
+:error_tol_(that.error_tol_),
+cutoff_radius_policy_(that.cutoff_radius_policy_),
+update_policy_(that.update_policy_),
+bandwidth_(that.bandwidth_),
+max_degree_(that.max_degree_),
+max_recursions_(0),
+cutoff_radius_(that.cutoff_radius_),
+max_cluster_radius_(that.max_cluster_radius_),
+error_(that.error_),
+recursions_(0){}
+
+template<typename RealType, template<typename> class CutoffRadiusPolicy,
+template<typename> class UpdatePolicy>
+typename
+optimal_parameter_given_max_degree<RealType,CutoffRadiusPolicy,
+UpdatePolicy>::optimal_parameter_given_max_degree&
+optimal_parameter_given_max_degree<RealType,CutoffRadiusPolicy,
+UpdatePolicy>::operator=(
+ const optimal_parameter_given_max_degree& that){
+ if(&that!=this){
+ error_tol_ = that.error_tol_;
+ cutoff_radius_policy_ = that.cutoff_radius_policy_;
+ update_policy_ = that.update_policy_;
+ bandwidth_ = that.bandwidth_;
+ max_degree_ = that.max_degree_;
+ max_recursions_ = that.max_recursions_;
+ cutoff_radius_ = that.cutoff_radius_;
+ max_cluster_radius_ = that.max_cluster_radius_;
+ error_ = that.error_;
+ recursions_ = that.recursions_;
+ }
+ return *this;
+};
+
+template<typename RealType,template<typename> class CutoffRadiusPolicy,
+template<typename> class UpdatePolicy>
+void optimal_parameter_given_max_degree<RealType,CutoffRadiusPolicy,
+UpdatePolicy>::operator()(
+ unsigned max_degree,
+ value_type bandwidth,
+ value_type start_cluster_radius_over_bandwidth,
+ size_type max_recursions ){
+ bandwidth_ = bandwidth;
+ max_degree_ = max_degree;
+ max_recursions_ = max_recursions;
+ recursions_ = 0;
+ max_cluster_radius_ = start_cluster_radius_over_bandwidth
+ /this->bandwidth();
+ call_impl();
+};
+
+template<typename RealType,template<typename> class CutoffRadiusPolicy,
+template<typename> class UpdatePolicy>
+void optimal_parameter_given_max_degree<RealType,CutoffRadiusPolicy,
+UpdatePolicy>::call_impl(){
+ ++recursions_;
+ cutoff_radius_ = cutoff_radius_policy_(*this);
+ error_ = truncation_properties::max_error_bound(
+ max_degree(), bandwidth(),
+ max_radius(), cutoff_radius()
+ );
+
+ if(recursions()<max_recursions()){
+ if(error()>error_tolerance()){
+ max_cluster_radius_ = update_policy_(max_radius());
+ ++recursions_;
+ call_impl();
+ }
+ }else{
+ std::string str = "optimal_parameter_given_max_degree";
+ str += "max number of recursions exceeded";
+ throw std::runtime_error(str);
+ }
+}
+
+}}}
+
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/rest_weights_wrapper.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/rest_weights_wrapper.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,28 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/rest_weights_wrapper.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_REST_WEIGHTS_WRAPPER_HPP_ER_2009
+#define BOOST_MATH_IFGT_REST_WEIGHTS_WRAPPER_HPP_ER_2009
+namespace boost{namespace math{namespace ifgt{
+ /// This wrapper indicates that each weight in the sequence
+ /// has yet to be prepend by one in order to be NW compliant.
+ template<typename R1>
+ class rest_weights_wrapper{
+ public:
+ rest_weights_wrapper(const R1& r):r_(r){}
+ rest_weights_wrapper(const rest_weights_wrapper& that):r_(that.r_){}
+ const R1& operator()()const{ return r_; }
+ private:
+ rest_weights_wrapper();
+ const R1& r_;
+ };
+ template<typename R1>
+ rest_weights_wrapper<R1> make_rest_weights_wrapper(const R1& r){
+ typedef rest_weights_wrapper<R1> result_type;
+ return result_type(r);
+ }
+}}}
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/tag.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/tag.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,96 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/tag.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_TAG_HPP_ER_2009
+#define BOOST_MATH_IFGT_TAG_HPP_ER_2009
+#include <boost/parameter/parameters.hpp>
+#include <boost/parameter/keyword.hpp>
+namespace boost{namespace math{namespace ifgt{
+namespace{
+ struct bandwidth_{};
+ struct max_cluster_radius_{};
+ struct degree_{};
+ struct clusters_{};
+ struct center_{};
+ struct target_{};
+ struct range_out_{};
+ struct error_tolerance_{};
+ struct max_distance_source_target_{};
+ struct cutoff_radius_policy_{};
+ struct clusters_reserve_size_{};
+ struct contributions_reserve_size_{};
+ struct accumulator_{};
+ struct nadaraya_watson_check_tolerance_{};
+ struct divide_by_{};
+}
+namespace tag{
+
+ static
+ ::boost::parameter::keyword<bandwidth_>& bandwidth
+ = ::boost::parameter::keyword<bandwidth_>::get();
+ static
+ ::boost::parameter::keyword<max_cluster_radius_>&
+ max_cluster_radius
+ =
+ ::boost::parameter::keyword<max_cluster_radius_>::get();
+ static
+ ::boost::parameter::keyword<degree_>& degree
+ = ::boost::parameter::keyword<degree_>::get();
+ static
+ ::boost::parameter::keyword<clusters_>& clusters
+ = ::boost::parameter::keyword<clusters_>::get();
+ static
+ ::boost::parameter::keyword<center_>& center
+ = ::boost::parameter::keyword<center_>::get();
+ static
+ ::boost::parameter::keyword<target_>& target
+ = ::boost::parameter::keyword<target_>::get();
+ static
+ ::boost::parameter::keyword<range_out_>& range_out
+ = ::boost::parameter::keyword<range_out_>::get();
+ static
+ ::boost::parameter::keyword<error_tolerance_>& error_tolerance
+ = ::boost::parameter::keyword<error_tolerance_>::get();
+ static
+ ::boost::parameter::keyword<max_distance_source_target_>&
+ max_distance_source_target
+ =
+ ::boost::parameter::keyword<max_distance_source_target_>::get();
+ static
+ ::boost::parameter::keyword<cutoff_radius_policy_>&
+ cutoff_radius_policy
+ = ::boost::parameter::keyword<cutoff_radius_policy_>::get();
+ static
+ ::boost::parameter::keyword<clusters_reserve_size_>&
+ clusters_reserve_size
+ = ::boost::parameter::keyword<clusters_reserve_size_>::get();
+ static
+ ::boost::parameter::keyword<contributions_reserve_size_>&
+ contributions_reserve_size
+ =
+ ::boost::parameter::keyword<contributions_reserve_size_>::get();
+ static
+ ::boost::parameter::keyword<divide_by_>&
+ divide_by
+ = ::boost::parameter::keyword<
+ divide_by_>::get();
+ static
+ ::boost::parameter::keyword<accumulator_>&
+ accumulator
+ = ::boost::parameter::keyword<accumulator_>::get();
+
+ static
+ ::boost::parameter::keyword<nadaraya_watson_check_tolerance_>&
+ nadaraya_watson_check_tolerance
+ = ::boost::parameter::keyword<
+ nadaraya_watson_check_tolerance_>::get();
+// All this works as is, but
+// defining instead a keyword within say struct my_struct{}
+// has caused compile error "multiple definitions of tag::my_struct::keyword".
+// I have used defs within structs before w/o pbs (though templates) so why?
+}
+}}}
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/traits.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/traits.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,29 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/traits.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_TRAITS_HPP_ER_2009
+#define BOOST_MATH_IFGT_TRAITS_HPP_ER_2009
+#include <boost/mpl/size_t.hpp>
+namespace boost{namespace math{namespace ifgt{
+
+template<typename T>
+struct stride_traits{
+ typedef mpl::size_t<T::source_size> stride0_type;
+ typedef mpl::size_t<T::weight_size> stride1_type;
+};
+
+///These traits are for use by for_each_over_vector2matrix
+template<typename T>
+struct traits : public stride_traits<T>{
+ typedef T type;
+};
+
+template<typename T>
+struct call_traits : T::stride_traits{
+ typedef T type;
+};
+}}}
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/truncation_degree_constant.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/truncation_degree_constant.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,50 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/truncation_degree_constant.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_ASSIGN_DEGREES_CONSTANT_HPP_ER_2009
+#define BOOST_MATH_IFGT_ASSIGN_DEGREES_CONSTANT_HPP_ER_2009
+#include <algorithm>
+#include <boost/range.hpp>
+#include <boost/parameter/parameters.hpp>
+#include <boost/parameter/keyword.hpp>
+#include <boost/math/ifgt/tag.hpp>
+namespace boost{namespace math{namespace ifgt{
+
+ /// Models TruncationDegreesPolicy
+ template<typename RealType = double>
+ class truncation_degree_constant{
+ public:
+ typedef RealType value_type;
+
+ truncation_degree_constant():p_(20){}
+ explicit truncation_degree_constant(unsigned p):p_(p){}
+ //default copy/assign should be fine
+
+ template<typename ArgPack>
+ truncation_degree_constant(const ArgPack& args)
+ :p_(args[tag::degree|20]){}
+
+ unsigned degree()const{return p_;}
+
+ template<typename Cluster,typename R0,typename R1>
+ void operator()(
+ const Cluster& cluster,
+ value_type radius_source_to_center,
+ const R0& weight,
+ R1& degrees
+ ){
+ std::fill(
+ begin(degrees),
+ end(degrees),
+ degree()
+ );
+ }
+ private:
+ unsigned p_;
+ };
+
+}}}
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/truncation_properties.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/truncation_properties.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,158 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/truncation_properties.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_IFGT_TRUNCATION_PROPERTIES_HPP_ER_2009
+#define BOOST_MATH_IFGT_TRUNCATION_PROPERTIES_HPP_ER_2009
+#include <cmath>
+#include <boost/math/special_functions/factorials.hpp>
+namespace boost{namespace math{namespace ifgt{
+
+struct truncation_properties{
+ public:
+ //assumes weight = 1
+ template<typename RealType>
+ static RealType log_error_bound(
+ unsigned max_degree,
+ RealType bandwidth,
+ RealType source_to_center_distance,
+ RealType target_to_center_distance
+ );
+ template<typename RealType>
+ static RealType error_bound(
+ unsigned max_degree,
+ RealType bandwidth,
+ RealType source_to_center_distance,
+ RealType target_to_center_distance
+ );
+ template<typename RealType>
+ static RealType worst_target_to_center_distance(
+ unsigned max_degree,
+ RealType bandwidth,
+ RealType source_to_center_distance
+ );
+ template<typename RealType>
+ static RealType worst_target_to_center_distance(
+ unsigned max_degree,
+ RealType bandwidth,
+ RealType source_to_center_distance,
+ RealType cutoff_radius
+ );
+ template<typename RealType>
+ static RealType max_log_error_bound(
+ unsigned max_degree,
+ RealType bandwidth,
+ RealType source_to_center_distance,
+ RealType cutoff_radius
+ );
+ template<typename RealType>
+ static RealType max_error_bound(
+ unsigned max_degree,
+ RealType bandwidth,
+ RealType source_to_center_distance,
+ RealType cutoff_radius
+ );
+
+};
+
+template<typename RealType>
+RealType truncation_properties::log_error_bound(
+ unsigned max_degree,
+ RealType bandwidth,
+ RealType source_to_center_distance,
+ RealType target_to_center_distance
+){
+ typedef RealType value_type;
+ value_type p = max_degree;
+ value_type h = bandwidth;
+ value_type rx = source_to_center_distance;
+ value_type ry = target_to_center_distance;
+ value_type result =
+ max_degree * (
+ log(2)
+ + log(rx)
+ + log(ry)
+ - 2 * log(h )
+ ) - log ( factorial<value_type>(p) );
+ result -= std::pow( (rx - ry)/h, 2 );
+ return result;
+}
+
+template<typename RealType>
+RealType truncation_properties::error_bound(
+ unsigned max_degree,
+ RealType bandwidth,
+ RealType source_to_center_distance,
+ RealType target_to_center_distance
+){
+ return exp(log_error_bound(
+ max_degree,
+ bandwidth,
+ source_to_center_distance,
+ target_to_center_distance
+ ));
+}
+
+template<typename RealType>
+RealType
+truncation_properties::worst_target_to_center_distance(
+ unsigned max_degree,
+ RealType bandwidth,
+ RealType source_to_center_distance
+){
+ typedef RealType value_type;
+ value_type p = max_degree;
+ value_type h = bandwidth;
+ value_type rx = source_to_center_distance;
+ value_type result
+ = ( rx + sqrt(std::pow(rx,2) + 2 * p * std::pow(h,2)) )
+ /((value_type)(2));
+ return result;
+}
+
+template<typename RealType>
+RealType
+truncation_properties::worst_target_to_center_distance(
+ unsigned max_degree,
+ RealType bandwidth,
+ RealType source_to_center_distance,
+ RealType cutoff_radius
+){
+ typedef RealType value_type;
+ value_type result = worst_target_to_center_distance(
+ max_degree, bandwidth, source_to_center_distance);
+
+ return (result<cutoff_radius)? result : cutoff_radius;
+}
+
+template<typename RealType>
+RealType truncation_properties::max_log_error_bound(
+ unsigned max_degree,
+ RealType bandwidth,
+ RealType source_to_center_distance,
+ RealType cutoff_radius
+){
+ typedef RealType value_type;
+ value_type wry =
+ worst_target_to_center_distance(
+ max_degree, bandwidth, source_to_center_distance, cutoff_radius);
+ return log_error_bound(
+ max_degree, bandwidth, source_to_center_distance, wry);
+}
+
+template<typename RealType>
+RealType truncation_properties::max_error_bound(
+ unsigned max_degree,
+ RealType bandwidth,
+ RealType source_to_center_distance,
+ RealType cutoff_radius
+){
+ return exp(
+ max_log_error_bound(
+ max_degree, bandwidth, source_to_center_distance, cutoff_radius));
+}
+
+}}}
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/math/ifgt/zscore.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/math/ifgt/zscore.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,56 @@
+//////////////////////////////////////////////////////////////////////////////
+// ifgt/zscore.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_MATH_ZSCORE_HPP_ER_2009
+#define BOOST_MATH_ZSCORE_HPP_ER_2009
+#include <algorithm>
+#include <functional>
+#include <stdexcept>
+#include <boost/range.hpp>
+#include <boost/assert.hpp>
+#include <boost/bind.hpp>
+namespace boost{namespace math{
+
+ /// Maps x,y,h to (x-y)/h
+ template<typename R0>
+ class zscore{
+ typedef typename range_iterator<R0>::type iter_type;
+ public:
+ typedef typename iter_type::value_type value_type;
+ zscore(const R0& center_):center(center_){}
+ zscore(const zscore& that):center(that.zscore){}
+
+ template<typename R1,typename RealType>
+ void operator()(const R1& source,RealType bandwidth,R0& out){
+ BOOST_ASSERT(size(source)==size(center));
+
+ std::transform(
+ begin(source),
+ end(source),
+ begin(center),
+ begin(out),
+ bind(
+ std::divides<value_type>(),
+ bind(
+ std::minus<value_type>(),
+ _1,
+ _2
+ ),
+ (value_type)(bandwidth)
+ )
+ );
+ }
+ private:
+ zscore& operator=(const zscore& that);
+ zscore();
+ const R0& center;
+ };
+
+ template<typename R0>
+ zscore<R0> make_zscore(const R0& center){return zscore<R0>(center);}
+}}
+
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/utility/dont_care.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/utility/dont_care.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,20 @@
+//////////////////////////////////////////////////////////////////////////////
+// dont_care.cpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_UTILITY_DONT_CARE_HPP_ER_2009
+#define BOOST_UTILITY_DONT_CARE_HPP_ER_2009
+namespace boost{namespace utility{
+ //borrowed from Boost.Accumulator
+ struct dont_care
+ {
+ template<typename Args>
+ dont_care(Args const &)
+ {
+ }
+ };
+
+}}
+#endif

Added: sandbox/improved_fast_gauss_transform/boost/utility/nested_type.hpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/boost/utility/nested_type.hpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,24 @@
+///////////////////////////////////////////////////////////////////////////////
+// nested_type.hpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef BOOST_UTILITY_NESTED_TYPE_HPP_ER_2009
+#define BOOST_UTILITY_NESTED_TYPE_HPP_ER_2009
+#include <boost/mpl/has_xxx.hpp>
+#include <boost/mpl/eval_if.hpp>
+
+namespace boost{namespace utility{
+ BOOST_MPL_HAS_XXX_TRAIT_DEF(type);
+
+// See news://news.gmane.org:119/gma8hc$2a0$1@ger.gmane.org
+template<class T>
+struct nested_type :
+ boost::mpl::eval_if<has_type<T>,
+ T,
+ boost::mpl::identity<T> >
+{};
+
+}}
+#endif

Added: sandbox/improved_fast_gauss_transform/libs/math/ifgt/doc/readme.txt
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/libs/math/ifgt/doc/readme.txt 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,481 @@
+///////////////////////////////////////////////////////////////////////////////
+// improved_fast_gauss_transform
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+///////////
+/ Contact /
+///////////
+
+Please send questions or suggestions to erwann.rogard_at_[hidden]
+
+/////////////
+/ Overview /
+/////////////
+
+This collection of C++ classes implements the Improved Fast Gauss
+Transform (see Raykar2006a), and its exact counterpart. Applications
+include the Rozenblatt--Parzen (RP) density estimator and the Nadaraya--
+Watson (NW) conditional mean estimator.
+
+This package offers the specificity of allowing sequential accumulation of
+training data, and customization of data types and particulars of the
+algorithm through templates.
+
+Below is a description of the package, but individual *.hpp files contain
+additional information.
+
+Examples are in libs/ifgt/src/example.
+
+//////////////////
+/ Requirements /
+//////////////////
+
+Compiles under
+
+gcc version i686-apple-darwin9-gcc-4.0.1 (GCC) 4.0.1 (Apple Inc. build 5490)
+
+The compiler search path must include
+boost_1_37_0
+boost/sandbox/miscellanea_iterator_facade
+boost/sandbox/miscellanea_algorithm
+boost/sandbox/monomials_horner_rule
+boost/sandbox/improved_fast_gauss_transform
+
+///////////////
+/ Conventions /
+///////////////
+
+D source size
+a multivariate power
+|a| total degree of a |a| = sum{a(d):d=0,...,D-1}
+p a bound on |a|
+b coefficient
+
+x source dim(x)=D
+y target dim(y)=D
+w weight dim(w)=J
+(x,w) contribution
+h bandwith
+i,N contributions index i=0,...,N-1
+j,J weight index j=0,...,J-1
+m,M target index m=0,...,M-1
+
+rx max cluster radius
+ry cutoff radius
+c cluster center
+zx zscore of x (x-c)/h
+zy zscore of y (t-c)/h
+C a cluster (h,{b:a},cc)
+k,K index of clusters
+
+{x:i} is short for {x_i:i=0,1,2,...,N}
+
+Gu(y) Gauss kernel (unnormalized) exp(||x-y||/h)
+nc Normalizing constant (2 * pi * h^2 )^(d/2)
+G(y) Gauss kernel Gu(y)/nc
+gt(y) Gauss transform (u) at y sum{w G(y):i}
+rp(y) Rozenblatt Parzen estimate gt(y)/N
+rp(y)|w=1 RP estimate assuming w=1
+nw(y) Nadaray Watson estimate rp(y) / rp(y)|w=1
+
+
+
+The quantity nw(y) and rp(y) is an estimate of E[w|x=y].
+In the special case that w=1, rp(y) is an estimate of pdf(x=y).
+The quantity nw(y) and rp(y) is an estimate of E[w|x=y].
+By convention, if we are to use nw, we will take w[0]=1.
+
+Some literature defines the Gauss transform as sum{w Gu(y):i}.
+Our definition provides the appropriate scaling in case h varies with i.
+
+/////////////
+/ Algorithm /
+/////////////
+
+* The basic idea:
+
+First assume h_i = h, J=1 and fix {p_i:i}. The Fast Gauss transform
+is based on a multivariate Taylor expansion centered at c, such that the
+contribution of (x,y), can be separated into two components, one associated
+with (x,c), and the other with (c,t):
+ gta(y) = sum{b_a (1/h)^d exp|| (c-y)/h ||^2 (y-c/h)^a:|a|< p+1}
+converges to g(y) as min{p_i:i} -> infty, where
+ b_a = 2^a/a! sum{w_i exp(-||x_i-c||/h):i such that |a|<p_i+1}
+which has complexity O(N). In particular, b_a is independent of y,
+so to evaluate at {y:m=1,...,M}, total complexity is O(N+M).
+
+* Space division
+
+The greater ||x-c||, the greater p=min{p_i:i}, in oder
+to keep the approximation error within a given bound. To keep p manageable,
+the sources are split into K clusters, and the ifgt performed within each of
+them.
+
+Examples: Given a pre-set K, farthest point has complexity O(K N),
+two phase (Feder and Greene, 1988) has complexity O(N log(K)).
+
+* Degree truncation - the first source of error
+
+See, for example, Section 4.2 in Raykar2006a
+
+* Cutoff radius - the second source error
+
+To further lessen computations, fix ry, and ignore clusters such that
+||c_k-y||>ry. Note that their contribution is less than exp(-ry).
+
+As per Raykar2006a, Section 4.1, to keep the error within eps:
+ ry > rx + min(R,h sqrt(log(1/eps)),
+where R=max{||x-y||:all i,m}
+
+* Relaxing J=1
+
+We can allow multiple weight sequences i.e. each training data point is
+ (x,{w:j}). A separate approximation is needed for for each j, but common
+computations are mutualized.
+
+* Relaxing h_i=h
+
+If h_i be piecewise constant, the ifgt is computed within each set, and the
+sum taken over each of the set members. We call "active
+range" the range of clusters of the current h.
+
+* Accumulating data sequentially
+
+The training data may be available all at once, or it may be necessary to
+collect the sources' contributions sequentially (more general), whereby a new
+cluster is created, if the nearest to that of the current input is
+is farther than some pre--set value, rx.
+
+* Monomials:
+The IFGT coefficients and the approximation involve the monomials
+{z^a:|a|< p+1}, for z=(x-c)/h, and z=(y-c)/h, respectively. For efficiency,
+Hoerner's rule is used.
+
+* Number of clusters (or maximum cluster radius), K
+
+Given h, Raykar2006a, Section 4.3 provide an algorithm that find the optimal
+K, equivalently rx. This relies on certain assumptions and we have not
+implemented it yet (TODO). As of know, we fix max{p_i:i} and solve for rx.
+
+* Choosing h
+
+For w=1, under normality of x, the statistically optimal bandwidth is
+
+h_{N}=sigma(N(2+D)/4)^{-1/(D+4)} * sqrt(2)
+
+The sqrt(2) is needed because the kernel is
+exp(-||(x-y)/h||^2), not the usual normal kernel exp(-(1/2)||(x-y)/h||^2)
+
+A common approach is to transform x[i]<-x[i]/std(x[i]) allowing sigma = 1.
+
+See Silverman1986, Section 4.3.2
+
+* Subsetting
+We may want to evaluate the gt only for a subset of {j=0,...,J-1}. For
+that reason, evaluators take a subset range policy.
+
+//////////
+/ Design /
+//////////
+
+Our two main concepts are Accumulator and Evaluator. The first takes
+care of accumulating training data in a sequential fashion:
+ a(source,weight)
+The second, usually built upon the former (by const ref) takes care of
+evaluating gt(y):
+ e(target,range_out).
+Each of Accumulator and Evaluator inherit from a base class with corres-
+ponding names (CRTP), for overload resolution purposes (e.g. for_each) or to
+provide additional functionality (e.g. rp(y) and nw(y) in the case
+of an evaluator).
+
+The operations of accumulating and evaluating can be vectorized:
+ for_each(sources,weights,a)
+ for_each(targets,ranges_out,e,call<method>())
+where each of sources,weights, targets, and ranges_out are "flattened"
+matrices. To be able to use Nadaray Watson, we need to prepend
+each weight with 1:
+ for_each(sources,make_rest_weights_wrapper(w),a)
+
+Models of Accumulator include fast_accumulator and exact_accumulator.
+The former finds the nearest among a set of active clusters and forwards the
+job of collecting data to that cluster.
+
+The template class cluster is parameterized by a degree truncation policy,
+models of which include degree_truncation_rydg.
+
+Models of Evaluator include fast_evaluator and exact_evaluator. Naturally,
+a fast_evaluator is built around a fast_accumulator. It is parameterized by
+a cutoff radius policy, a model of which is cutoff_radius_rydg.
+
+End user classes may contain deeply nested components. Therefore,
+for convenience, their constructor takes an argument pack
+(see the Boost Parameter library) which is passed down into the class
+components. For example:
+ A a((
+ tag::bandwidth = h,
+ tag::max_cluster_radius = rx,
+ tag::degree = p
+ )
+ );
+
+Instead of template template parameters we use lambda
+expressions to make it easier to allow, for example, G<F<X,_1>,Y>
+
+//////////////
+/ Concepts /
+//////////////
+
+Cluster :
+Associated type:
+ result_of_center
+
+
+Accumulator (A)
+Associated types:
+ A has public base accumulator_base<A>
+ A::value_type
+ A::sources_count_type
+Valid expressions: let a an lvalue of A, and source and weight those
+of a ForwardRange,
+ a.source_size const std::size_t
+ a.weight_size const std::size_t
+ active_bandwidth() value_type
+ sources_count() sources_count_type
+ e(source, weight) void
+Models: fast_accumulator.hpp, exact_accumulator.hpp
+
+Evaluator (E), e an lvalue of E, C a type,
+ E has public base evaluator_base<C,E>
+Associated types:
+ E::sources_count_type
+Valid expressions: let target and range_out model ForwarRange.
+ e.sources_count() sources_count_type
+ e(target,range_out) void
+Models: fast_evaluator, exact_evaluator
+Comment: If A is the type of the accumulator that E builds upon,
+C = evaluator_common<A> takes care of sources_count.
+
+CuttoffRadiusPolicy (CR): f an an lvalue,
+Associated type:
+ value_type
+Valid expressions:
+h and rx instances of value_type
+ f(bandwidth,max_radius) returns a RealType
+
+TruncateDegreePolicy: f an lvalue,
+Associated type:
+ value_type
+Valid expression: rx a value_type,
+ f(C,rx,w,p) assigns a value to each member of p (one for each weight).
+
+FindClusterPolicy (F)
+Associate types:
+r and source lvalues of a forward range R
+F f(r) constructor
+f(x) returns an iterator pointing to an element of r.
+
+SubsetRangePolicy (S),
+Associated types: R0 a ForwardRange,
+ R= result_of<S(R0)>::type (another ForwardRange)
+
+Valid expressions:
+ f(r0) returns an lvalue of R
+
+////////////////////
+/ Files /
+////////////////////
+
+* accumulator_base.hpp
+* benchmark.hpp
+* call_gauss_transform.hpp
+* call_nadaraya_watson_estimate.hpp
+* call_rozenblatt_parzen_estimate.hpp
+* call_wrapper.hpp
+* cluster.hpp
+ vars: x, h, {{b:a}:j}
+ arg: x, {w:j},{p:j}
+ effect: updates {{b:a}:j}
+* cluster_call_center.hpp
+* cluster_evaluator.hpp
+ vars: t, {out:j}
+ arg: C
+ effect: adds {gt(t):j} to {out:j}
+* clusters_evaluator.hpp
+ arg: t, {C},{out:j}
+ effect: iterates over {C} and adds {gt(t):j} to {out:j}
+* coefficients.hpp
+ var: {b:a}
+ arg: ({zx^a},w)
+ effect: updates {b:a}
+* coefficients_evaluator.hpp
+ var: t
+ arg: r, cc,h, {b:a}
+ out: gt(t)
+* cutoff_radius_none.hpp
+* cutoff_radius_rydg.hpp
+* evaluator_base.hpp
+* evaluator_common.hpp
+* exact_accumulator.hpp
+* exact_contribution.hpp
+ var: x, {w:j}, h
+ member: evaluate
+ arg: t, out
+ effect: adds {gt(t):j} to out
+* exact_contributions_evaluator.hpp
+* exact_evaluator.hpp
+* fast_accumulator.hpp
+* fast_evaluator.hpp
+* find_nearest_cluster.hpp:
+ var: {cc}
+ arg: x
+ out: cc nearest to x
+* for_each_accumulate.hpp
+* for_each_evaluator.hpp
+* include.hpp
+* monomials.hpp:
+ arg: x, p
+ out: {x^a: |a|=p}
+* monomials_properties.hpp
+* multi_factorial.hpp
+ deprecated: redundant with multi_indexes_derived.
+* multi_indexes.hpp
+ member: get
+ arg: p
+ out: {a:|a|=p}
+* multi_indexes_derived.hpp
+* optimal_bandwidth.hpp
+* optimal_parameter_given_max_degree.hpp
+* normal_kernel_properties.hpp
+* tag.hpp
+* traits.hpp
+* truncation_degree_constant.hpp
+* truncation_properties.hpp
+
+////////
+/ TODO /
+////////
+This is an early version that begs for safeguards, optimization and testing.
+
+As far as design:
+- Customize the exact gauss transform by
+adding a policy for the cutoff radius as well.
+- Add helper classes to predict how parameters affect resource
+requirement (memory, in particular, easily becomes a problem) and precision.
+- For now, the bandwidth is optimized only in a statistical sense as a
+function of N. What is needed is a per-computation-time optimization.
+- Optimal p, K (see remark in Algorithm)
+
+/////////////
+/ Reference /
+/////////////
+
+@techreport{Raykar2006a,
+ title = {Fast computation of sums of gaussians in high dimensions},
+ author = {Vikas Chandrakant Raykar and Changjiang Yang and Ramani
+ Duraiswami and Nail Gumerov},
+ institution = {Department of Computer Science and Institute for Advanced
+ Computer
+ Studies, University of Maryland},
+ month = {April},
+ number = {CS-TR-4767/UMIACS-TR-2005-69},
+ url = {http://www.umiacs.umd.edu/~vikas/publications/CS-TR-4767.pdf},
+ year = {2006},
+ biburl =
+ {http://www.bibsonomy.org/bibtex/21702c6cf15de590feaed72075225555a/
+ marcoalvarez},
+ abstract = {Evaluating sums of multivariate Gaussian kernels is a key
+ computational
+ task in many problems in computational statistics and machine learning.
+ The computational cost of the direct evaluation of such sums scales
+ as the product of the number of kernel functions and the evaluation
+ points. The fast Gauss transform proposed by Greengard and Strain
+ (1991) is a e-exact approximation algorithm that reduces the computational
+ complexity of the evaluation of the sum of N Gaussians at M points
+ in d dimensions from O(MN) to O(M+N). However, the constant factor
+ in O(M+N) grows exponentially with increasing dimensionality d, which
+ makes the algorithm impractical for dimensions greater than three.
+ In this paper we present a new algorithm where the constant factor
+ is reduced to asymptotically polynomial order. The reduction is based
+ on a new multivariate Taylor's series expansion (which can act both
+ as a local as well as a far field expansion) scheme combined with
+ the efficient space subdivision using the k-center algorithm. The
+ proposed method differs from the original fast Gauss transform in
+ terms of a different factorization, effcient space subdivision, and
+ the use of point-wise error bounds. Algorithm details, error bounds,
+ procedure to choose the parameters and numerical experiments are
+ presented. As an example we shows how the proposed method can be
+ used for very fast e-exact multivariate kernel density estimation.},
+ timestamp = {2007.11.15}, owner = {Marco},
+ keywords = {Gaussian }
+}
+
+@book{Silverman1986,
+ title = "Density Estimation for Statistics and Data Analysis",
+ author = "B. W. Silverman and P. J. Green",
+ publisher = "Chapman and Hall",
+ address = "London",
+ year = 1986
+}
+
+/////////////////
+/ Sample output /
+/////////////////
+pdf(x) = prod{N(x[d]|0,1),d=0,...,D} esimated by Rozenblatt-Parzen
+f(x) = norm(x) = sqrt(x[0]^2+...+x[D]^2) estimated by Nadaraya-Watson
+
+M: # evaluated
+N: # accumulated
+err0(a,b) = max {|a - b|:m=0,...,M-1 }
+err1(a,b) = sqrt sum {|a - b|^2:m=0,...,M-1 } / M
+e0_rp : err0(rp(y), pdf(y))
+e1_rp : err1(rp(y), pdf(y))
+e0_nw : err0(nw(y), w[1](y))
+e1_nw : err1(nw(y), w[1](y))
+Quantities below are scaled to reflect N=M=1e3
+ta : time accumulate {(x,w):i=0,...,N-1}
+trp : time {rp(y):m = 0,..., M-1}
+trp : time {nw(y):m = 0,..., M-1}
+
+
+
+-> example_benchmark_exact
+dim(x) = 3
+dim(w) = 2
+bandwidth = 0.547937
+N e0_rp e1_rp e0_nw e1_nw
+ 2,500,000 1.004e-02 3.508e-04 3.901e-01 1.670e-02
+ 5,000,000 9.965e-03 3.487e-04 3.911e-01 1.672e-02
+ 7,500,000 9.996e-03 3.494e-04 3.914e-01 1.672e-02
+10,000,000 9.988e-03 3.490e-04 3.928e-01 1.673e-02
+
+N ta trp tnw
+ 2,500,000 1.033e-02 1.257e+04 1.235e+04
+ 5,000,000 1.036e-02 2.514e+04 2.472e+04
+ 7,500,000 1.310e-02 3.772e+04 3.702e+04
+10,000,000 7.707e-03 7.316e+03 6.385e+03
+<-
+
+ -> example_benchmark_fast
+cutoff radius = 0
+dim(x) = 3
+dim(w) = 2
+bandwidth =0.547937
+max_rx=0.456257
+
+
+N e0_rp e1_rp e0_nw e1_nw
+ 2,500,000 4.004e-02 1.793e-03 3.930e-01 1.672e-02
+ 5,000,000 4.002e-02 1.793e-03 3.925e-01 1.673e-02
+ 7,500,000 4.003e-02 1.793e-03 3.924e-01 1.673e-02
+10,000,000 4.003e-02 1.793e-03 3.936e-01 1.673e-02
+
+N ta trp tnw K
+2500000 2.133e-01 2.252e+01 4.439e+01 2469
+5000000 2.491e-01 2.540e+01 4.948e+01 2750
+7500000 2.617e-01 2.686e+01 5.291e+01 2901
+10000000 2.705e-01 2.978e+01 5.626e+01 3017
+
+<-

Added: sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/benchmark_exact.cpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/benchmark_exact.cpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,91 @@
+//////////////////////////////////////////////////////////////////////////////
+// example/benchmark_exact.cpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#include <boost/math/ifgt/benchmark.hpp>
+#include <libs/math/ifgt/src/example/benchmark_exact.h>
+void example_benchmark_exact(){
+ std::cout << " -> example_benchmark_exact" << std::endl;
+ std::string str =
+ "pdf(x) = prod{N(x[d]|0,1),d=0,...,D} estimated by Rozenblatt-Parzen";
+ str+= "and f(x) = norm(x) = sqrt(x[0]^2+...+x[D]^2) by Nadaraya-Watson";
+ std::cout << str << std::endl;
+
+ using namespace boost;
+ namespace math = boost::math;
+ namespace ifgt = math::ifgt;
+
+ typedef double value_type;
+ typedef std::vector<value_type> var_type;
+ typedef mt19937 urng_type;
+ typedef normal_distribution<value_type> nd_type;
+ typedef variate_generator<urng_type&,nd_type> gen_type;
+ typedef std::size_t size_type;
+ typedef std::vector<size_t> sizes_type;
+ typedef math::normal_distribution<value_type> univ_distr_type;
+ typedef ifgt::product_pdf<univ_distr_type> pdf_fun_type;
+ typedef l2_norm w_fun_type;
+ // fixed constants (do not modify)
+ const value_type sigma = 1.0;
+ const size_type wdim = 2;
+
+ // other constants (you may modify them)
+ // Warning: if you get a bad alloc error during exec,
+ // increase loop_count and try again
+ const size_type sdim = 3;
+ const size_type train_count = 2.5e3;
+ const size_type test_count = 1e2;
+ const size_type loop_count = 4;
+ const value_type bandwidth_adjust = 3;
+
+ // dependent constant
+ size_type all_train_count = train_count * loop_count;
+
+ typedef ifgt::benchmark<
+ value_type,
+ test_count,
+ sdim,
+ wdim,
+ pdf_fun_type,
+ w_fun_type> bench_type;
+ typedef ifgt::optimal_bandwidth<sdim,value_type> opt_bandwidth_type;
+ typedef ifgt::exact_accumulator<sdim,wdim,var_type> exact_acc_type;
+ typedef ifgt::exact_evaluator<exact_acc_type> exact_eval_type;
+
+ //pdf
+ // Data generator
+ urng_type urng(0);
+ gen_type gen(urng,nd_type());
+
+ bench_type bench(
+ gen,
+ pdf_fun_type(univ_distr_type(0,1)),
+ w_fun_type()
+ );
+
+ std::cout << "dim(x) = " << sdim << std::endl;
+ std::cout << "dim(w) = " << wdim << std::endl;
+
+ value_type opt_h
+ = opt_bandwidth_type::for_density_estimation_under_normality(
+ sigma,all_train_count);
+ opt_h *= bandwidth_adjust; //override
+
+ std::cout << "bandwidth ="
+ << opt_h << std::endl;
+
+ exact_acc_type exact_acc((ifgt::tag::bandwidth = opt_h ));
+ exact_eval_type exact_eval((ifgt::tag::accumulator = exact_acc));
+
+ bench.notation(std::cout);
+ bench.header(std::cout); std::cout << std::endl;
+ for(unsigned i = 0; i<loop_count; i++){
+ bench.accumulate(gen,train_count,exact_acc);
+ bench.estimate_pdf(exact_eval);
+ bench.estimate_w(exact_eval);
+ bench.statistics(std::cout,i); std::cout << std::endl;
+ }
+
+}

Added: sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/benchmark_exact.h
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/benchmark_exact.h 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,10 @@
+//////////////////////////////////////////////////////////////////////////////
+// example/benchmark_exact.h
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef LIBS_MATH_IFGT_EXAMPLE_BENCHMARK_EXACT_H_ER_2009
+#define LIBS_MATH_IFGT_EXAMPLE_BENCHMARK_EXACT_H_ER_2009
+void example_benchmark_exact();
+#endif

Added: sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/benchmark_fast.cpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/benchmark_fast.cpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,141 @@
+//////////////////////////////////////////////////////////////////////////////
+// example/benchmark_fast.cpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#include <boost/range.hpp>
+#include <boost/math/ifgt/benchmark.hpp>
+#include <libs/math/ifgt/src/example/benchmark_fast.h>
+void example_benchmark_fast(){
+ std::cout << " -> example_benchmark_fast" << std::endl;
+ std::string str =
+ "pdf(x) = prod{N(x[d]|0,1),d=0,...,D} estimated by Rozenblatt-Parzen";
+ str+= "and f(x) = norm(x) = sqrt(x[0]^2+...+x[D]^2) by Nadaraya-Watson";
+ std::cout << str << std::endl;
+
+ using namespace boost;
+ namespace math = boost::math;
+ namespace ifgt = math::ifgt;
+
+ typedef double value_type;
+ typedef std::vector<value_type> var_type;
+ typedef mt19937 urng_type;
+ typedef normal_distribution<value_type> nd_type;
+ typedef variate_generator<urng_type&,nd_type> gen_type;
+ typedef std::size_t size_type;
+ typedef std::vector<size_t> sizes_type;
+ typedef math::normal_distribution<value_type> univ_distr_type;
+ typedef ifgt::product_pdf<univ_distr_type> pdf_fun_type;
+ typedef l2_norm w_fun_type;
+ // fixed constants (do not modify)
+ const value_type sigma = 1.0;
+ const size_type wdim = 2;
+
+ // Other constants (you may modify them)
+ // Warning: if you get a bad alloc error during exec,
+ // increase bandwidth_adjust or train_count
+ const size_type sdim = 3;
+ const size_type train_count = 2.5e3;
+ const size_type test_count = 1e2;
+ const size_type loop_count = 4;
+ const value_type bandwidth_adjust = 3;
+
+ //
+ //const size_type sdim = 3;
+ //const size_type train_count = 2.5e5;
+ //const size_type test_count = 1e2;
+ //const size_type loop_count = 4;
+ //const value_type bandwidth_adjust = 4;
+
+ // fast addtional parameters
+ const bool do_cutoff_rydg = false;
+ const value_type error_tol = 1e-2;
+ const size_type max_p = 20;
+ const value_type start_rx_div_h = 2;
+ const size_type max_recursions = 20;
+
+
+ // dependent constant
+ size_type all_train_count = train_count * loop_count;
+ // additional statistics
+ sizes_type clusters_count(0);
+
+ typedef ifgt::benchmark<
+ value_type,
+ test_count,
+ sdim,
+ wdim,
+ pdf_fun_type,
+ w_fun_type> bench_type;
+ typedef ifgt::optimal_bandwidth<sdim,value_type> opt_bandwidth_type;
+ //fast types
+ typedef ifgt::optimal_parameter_given_max_degree<value_type>
+ opt_pars_type;
+ typedef ifgt::truncation_degree_constant<mpl::_1> trunc_degree;
+ typedef ifgt::cluster<sdim,wdim,trunc_degree,var_type> cluster_type;
+ typedef ifgt::find_nearest_cluster<mpl::_1, boost::l2_distance_squared>
+ find_type;
+ typedef ifgt::fast_accumulator<cluster_type,find_type> fast_acc_type;
+ typedef mpl::if_<
+ mpl::bool_<do_cutoff_rydg>,
+ mpl::identity<ifgt::cutoff_radius_rydg<mpl::_1> >,
+ mpl::identity<ifgt::cutoff_radius_none<mpl::_1> >
+ >::type cutoff_policy;
+
+ typedef ifgt::fast_evaluator<fast_acc_type,cutoff_policy> fast_eval_type;
+
+ //pdf
+ // Data generator
+ urng_type urng(0);
+ gen_type gen(urng,nd_type());
+
+ bench_type bench(
+ gen,
+ pdf_fun_type(univ_distr_type(0,1)),
+ w_fun_type()
+ );
+ std::cout << "do cutoff = " << do_cutoff_rydg << std::endl;
+ std::cout << "dim(x) = " << sdim << std::endl;
+ std::cout << "dim(w) = " << wdim << std::endl;
+
+ value_type opt_h
+ = opt_bandwidth_type::for_density_estimation_under_normality(
+ sigma,all_train_count);
+ opt_h *= bandwidth_adjust; //override
+
+ //fast specific parameters
+ opt_pars_type opt_pars((ifgt::tag::error_tolerance = error_tol));
+ opt_pars(max_p, opt_h, start_rx_div_h, max_recursions);
+ value_type max_rx = opt_pars.max_radius();
+
+ std::cout << "bandwidth =" << opt_h << std::endl;
+ std::cout << " max_rx=" << max_rx << std::endl;
+
+ fast_acc_type fast_acc(
+ (
+ ifgt::tag::bandwidth = opt_h,
+ ifgt::tag::max_cluster_radius = max_rx,
+ ifgt::tag::degree = max_p,
+ ifgt::tag::error_tolerance = error_tol
+ )
+ );
+ fast_eval_type fast_eval((
+ ifgt::tag::accumulator = fast_acc,
+ ifgt::tag::error_tolerance = error_tol));
+
+ bench.notation(std::cout);
+ std::cout << "K : # clusters" << std::endl;
+ bench.header(std::cout);
+ std::cout << boost::format("%|1$-10|")%"K" << std::endl;
+ for(unsigned i = 0; i<loop_count; i++){
+ bench.accumulate(gen,train_count,fast_acc);
+ clusters_count.push_back(boost::size(fast_acc.clusters()));
+ bench.estimate_pdf(fast_eval);
+ bench.estimate_w(fast_eval);
+ bench.statistics(std::cout,i);
+ std::cout
+ << boost::format("%|1$-10|")% clusters_count.back()<<std::endl;
+ }
+
+}

Added: sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/benchmark_fast.h
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/benchmark_fast.h 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,10 @@
+//////////////////////////////////////////////////////////////////////////////
+// example/benchmark_fast.h
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef LIBS_MATH_IFGT_EXAMPLE_BENCHMARK_FAST_H_ER_2009
+#define LIBS_MATH_IFGT_EXAMPLE_BENCHMARK_FAST_H_ER_2009
+void example_benchmark_fast();
+#endif

Added: sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/cluster.cpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/cluster.cpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,86 @@
+///////////////////////////////////////////////////////////////////////////////
+// example/cluster.cpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#include <algorithm>
+#include <vector>
+#include <boost/assign/std/vector.hpp>
+#include <boost/range.hpp>
+#include <iostream>
+#include <boost/mpl/placeholders.hpp>
+#include <boost/assign/std/vector.hpp>
+#include <boost/range.hpp>
+#include <boost/math/ifgt/cluster.hpp>
+#include <libs/math/ifgt/src/example/cluster.h>
+
+void example_cluster(){
+ std::cout << "-> cluster" << std::endl;
+
+ namespace mpl = boost::mpl;
+ namespace math = boost::math;
+ namespace ifgt = math::ifgt;
+ namespace tag = ifgt::tag;
+
+ const unsigned dim = 2;
+ const unsigned wdim = 2;
+ typedef double value_type;
+ typedef std::vector<value_type> var_type;
+ typedef ifgt::truncation_degree_constant<mpl::_1> trunc_degree;
+ typedef ifgt::cluster<dim,wdim,trunc_degree,var_type> cluster_type;
+
+ typedef std::vector<unsigned> ints_type;
+ typedef cluster_type::coefficients_type coeffs_type;
+ typedef cluster_type::collection_coefficients_type coll_coeffs_type;
+
+ const value_type bandwidth = 0.1;
+ const value_type max_cluster_radius = 1.0;
+ const unsigned degree = 20;
+ ints_type degrees;
+ var_type center;
+ var_type source;
+ var_type weight;
+ {
+ using namespace boost::assign;
+ center += (-0.1), 0.1;
+ source += (-0.2), 0.2;
+ weight += 1.0, 0.1;
+ }
+
+ cluster_type cluster(
+ (tag::center=center,
+ tag::bandwidth=bandwidth,
+ tag::max_cluster_radius = max_cluster_radius,
+ tag::degree = degree) );
+
+ cluster_type cluster_cpy(cluster);
+
+ cluster(source,weight);
+
+ {
+ const coeffs_type& sc
+ = cluster.collection_coefficients()[0];
+
+ std::cout << "weigths[0]: coeffs = ";
+ copy(
+ boost::begin(sc()),
+ boost::end(sc()),
+ std::ostream_iterator<value_type>(std::cout, " ")
+ ); std::cout << std::endl;
+ }
+ {
+ const coeffs_type& sc
+ = cluster.collection_coefficients()[1];
+
+ std::cout << "weigths[1]: coeffs = ";
+ copy(
+ boost::begin(sc()),
+ boost::end(sc()),
+ std::ostream_iterator<value_type>(std::cout, " ")
+ ); std::cout << std::endl;
+ }
+
+
+ std::cout << "<-" << std::endl;
+}

Added: sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/cluster.h
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/cluster.h 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,10 @@
+//////////////////////////////////////////////////////////////////////////////
+// example//cluster.h
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef LIBS_MATH_IFGT_EXAMPLE_CLUSTER_H_ER_2009
+#define LIBS_MATH_IFGT_EXAMPLE_CLUSTER_H_ER_2009
+void example_cluster();
+#endif

Added: sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/cluster_evaluator.cpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/cluster_evaluator.cpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,89 @@
+///////////////////////////////////////////////////////////////////////////////
+// example/cluster_evaluator.cpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#include <algorithm>
+#include <vector>
+#include <boost/assign/std/vector.hpp>
+#include <boost/range.hpp>
+#include <iostream>
+#include <boost/assign/std/vector.hpp>
+#include <boost/range.hpp>
+#include <boost/mpl/placeholders.hpp>
+#include <boost/math/ifgt/cluster.hpp>
+#include <boost/math/ifgt/truncation_degree_constant.hpp>
+#include <boost/math/ifgt/cluster_evaluator.hpp>
+#include <boost/math/ifgt/cutoff_radius_none.hpp>
+#include <boost/math/ifgt/tag.hpp>
+#include <libs/math/ifgt/src/example/cluster_evaluator.h>
+void example_cluster_evaluator(){
+ std::cout << "-> cluster_evaluator" << std::endl;
+ namespace mpl = boost::mpl;
+ namespace math = boost::math;
+ namespace ifgt = math::ifgt;
+ namespace tag = ifgt::tag;
+
+ const unsigned int dim = 2;
+ const unsigned int wdim = 2;
+ typedef double value_type;
+ typedef std::vector<value_type> var_type;
+
+ typedef ifgt::truncation_degree_constant<mpl::_1> trunc_degree;
+ typedef ifgt::cluster<dim,wdim,trunc_degree,var_type> cluster_type;
+
+ typedef cluster_type::coefficients_type coeffs_type;
+ typedef cluster_type::collection_coefficients_type coll_coeffs_type;
+ typedef ifgt::cutoff_radius_none<value_type> cutoff_policy_type;
+ typedef ifgt::cluster_evaluator<
+ var_type,var_type,ifgt::cutoff_radius_none<mpl::_1> > cluster_eval_type;
+
+ const value_type bandwidth = 0.1;
+ const value_type max_cluster_radius = 1.0;
+ const unsigned degree = 20;
+ //ints_type degrees;
+ var_type center;
+ var_type source;
+ var_type target;
+ var_type weight;
+ {
+ using namespace boost::assign;
+ center += (-0.1), 0.1;
+ source += (-0.2), 0.2;
+ target += -0.15, 0.15;
+ weight += 1.0, 0.1;
+ //degrees += 5, 5;
+ }
+
+ cluster_type cluster((
+ tag::center = center,
+ tag::bandwidth = bandwidth,
+ tag::max_cluster_radius = max_cluster_radius,
+ tag::degree = degree
+ ) );
+ std::cout << "before?" << std::endl;
+ cluster(source,weight);
+ std::cout << "after?" << std::endl;
+
+ var_type range_out(wdim,(value_type)(0));
+
+ cutoff_policy_type cutoff_policy;
+ cluster_eval_type cluster_eval(
+ (
+ tag::range_out = range_out,
+ tag::target = target,
+ tag::cutoff_radius_policy = cutoff_policy
+ )
+ );
+
+ cluster_eval(cluster);
+
+ copy(
+ boost::begin(range_out),
+ boost::end(range_out),
+ std::ostream_iterator<value_type>(std::cout, " ")
+ );
+
+ std::cout << "<-" << std::endl;
+}

Added: sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/cluster_evaluator.h
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/cluster_evaluator.h 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,10 @@
+//////////////////////////////////////////////////////////////////////////////
+// example/cluster_evaluator.h
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef LIBS_MATH_IFGT_EXAMPLE_CLUSTER_EVALUATOR_H_ER_2009
+#define LIBS_MATH_IFGT_EXAMPLE_CLUSTER_EVALUATOR_H_ER_2009
+void example_cluster_evaluator();
+#endif

Added: sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/clusters_evaluator.cpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/clusters_evaluator.cpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,96 @@
+///////////////////////////////////////////////////////////////////////////////
+// example/clusters_evaluator.cpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#include <algorithm>
+#include <boost/math/special_functions/binomial.hpp>
+#include <boost/assign/std/vector.hpp>
+#include <boost/range.hpp>
+#include <iostream>
+#include <boost/math/ifgt/fast_accumulator.hpp>
+#include <boost/math/ifgt/truncation_degree_constant.hpp>
+#include <boost/math/ifgt/clusters_evaluator.hpp>
+#include <boost/math/ifgt/cutoff_radius_none.hpp>
+#include <libs/math/ifgt/src/example/clusters_evaluator.h>
+
+void example_clusters_evaluator(){
+ std::cout << "-> example_clusters_evaluator" << std::endl;
+
+ namespace mpl = boost::mpl;
+ namespace math = boost::math;
+ namespace ifgt = math::ifgt;
+ namespace tag = ifgt::tag;
+
+ const unsigned int dim = 2;
+ const unsigned int wdim = 2;
+ typedef double value_type;
+ typedef std::vector<value_type> var_type;
+ typedef ifgt::truncation_degree_constant<mpl::_1> trunc_degree;
+ typedef ifgt::cluster<dim,wdim,trunc_degree,var_type> cluster_type;
+ typedef ifgt::find_nearest_cluster<mpl::_1, boost::l2_distance_squared>
+ find_type;
+
+ typedef ifgt::fast_accumulator<cluster_type,find_type> acc_type;
+
+ value_type bandwidth = 0.1;
+ value_type max_cluster_radius = 1.0;
+ unsigned degree = 20;
+ var_type center0;
+ var_type center1;
+ var_type center2;
+ var_type source1;
+ var_type source2;
+ var_type weight;
+ var_type target;
+ {
+ using namespace boost::assign;
+ center0 += -1.1, 0.1;
+ center1 += 0.2, -0.2;
+ center2 += 0.1, 0.1;
+ source1 += 0.15,-0.15;//nearest to center1
+ source2 += 0.11, 0.11;//nearest to center2
+ weight += 1.0, 0.1;
+ target += (-0.15),0.15;
+ }
+
+ acc_type acc((
+ tag::bandwidth = bandwidth,
+ tag::max_cluster_radius = max_cluster_radius,
+ tag::degree = degree
+ )
+ );
+
+ acc(center0, weight);
+ std::cout
+ << "sources_count = " << acc.sources_count()
+ << std::endl;
+ acc(center1, weight);
+ std::cout
+ << "sources_count = " << acc.sources_count()
+ << std::endl;
+ acc(center2, weight);
+ std::cout
+ << "sources_count = " << acc.sources_count()
+ << std::endl;
+ acc(source1, weight);
+ std::cout
+ << "sources_count = " << acc.sources_count()
+ << std::endl;
+ acc(source2, weight);
+ std::cout
+ << "sources_count = " << acc.sources_count()
+ << std::endl;
+
+
+ typedef ifgt::clusters_evaluator<
+ value_type,ifgt::cutoff_radius_none<mpl::_1> > clusters_evaluator_type;
+
+
+ clusters_evaluator_type clusters_evaluator;
+ //TODO to be continued or delete the whole file
+
+ std::cout << "<-" << std::endl;
+
+}

Added: sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/clusters_evaluator.h
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/clusters_evaluator.h 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,10 @@
+//////////////////////////////////////////////////////////////////////////////
+// example/clusters_evaluator.h
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef LIBS_MATH_IFGT_EXAMPLE_CLUSTERS_EVALUATOR_H_ER_2009
+#define LIBS_MATH_IFGT_EXAMPLE_CLUSTERS_EVALUATOR_H_ER_2009
+void example_clusters_evaluator();
+#endif

Added: sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/coefficients.cpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/coefficients.cpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,69 @@
+///////////////////////////////////////////////////////////////////////////////
+// example/coefficients.cpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#include <algorithm>
+#include <vector>
+#include <boost/assign/std/vector.hpp>
+#include <boost/range.hpp>
+#include <iostream>
+#include <boost/assign/std/vector.hpp>
+#include <boost/range.hpp>
+#include <boost/math/monomials.hpp>
+#include <boost/math/ifgt/zscore.hpp>
+#include <boost/math/ifgt/coefficients.hpp>
+#include <libs/math/ifgt/src/example/coefficients.h>
+
+void example_coefficients(){
+ std::cout << "-> coefficients" << std::endl;
+ using namespace boost;
+ const unsigned int dim = 2;
+ typedef double value_type;
+ typedef std::vector<value_type> var_type;
+ typedef std::vector<unsigned> ints_type;
+ typedef math::ifgt::coefficients<dim> coeffs_type;
+
+ double bandwidth = 0.1;
+ int degree = 5;
+ var_type center;
+ var_type source;
+ value_type weight = 0.1;
+ {
+ using namespace boost::assign;
+ center += (-0.1), 0.1;
+ source += (-0.2), 0.2;
+ }
+
+ var_type zscore_val(dim);
+ math::ifgt::monomials<var_type> monoms;
+ {
+ typedef math::zscore<var_type> zscore_type;
+ zscore_type zscore_op(center);
+ zscore_op(source,bandwidth,zscore_val);
+ }
+ value_type zscore_sqnorm = std::inner_product(
+ boost::begin(zscore_val),
+ boost::end(zscore_val),
+ boost::begin(zscore_val),
+ (value_type)(0)
+ );
+
+ monoms(zscore_val,degree);
+
+ coeffs_type coeffs;
+ coeffs(
+ zscore_sqnorm,
+ monoms(),
+ weight,
+ degree
+ );
+ copy(
+ boost::begin(coeffs()),
+ boost::end(coeffs()),
+ std::ostream_iterator<double>(std::cout, " ")
+ ); std::cout << std::endl;
+
+ std::cout << "<-" << std::endl;
+}

Added: sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/coefficients.h
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/coefficients.h 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,10 @@
+//////////////////////////////////////////////////////////////////////////////
+// example/coefficients.h
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef LIBS_MATH_IFGT_EXAMPLE_COLLECTED_SOURCES_H_ER_2009
+#define LIBS_MATH_IFGT_EXAMPLE_COLLECTED_SOURCES_H_ER_2009
+void example_coefficients();
+#endif

Added: sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/coefficients_evaluator.cpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/coefficients_evaluator.cpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,85 @@
+///////////////////////////////////////////////////////////////////////////////
+// example/coefficients_evaluator.cpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#include <algorithm>
+#include <vector>
+#include <boost/assign/std/vector.hpp>
+#include <boost/range.hpp>
+#include <iostream>
+#include <boost/assign/std/vector.hpp>
+#include <boost/range.hpp>
+#include <boost/math/monomials.hpp>
+#include <boost/math/ifgt/zscore.hpp>
+#include <boost/math/ifgt/normal_kernel_properties.hpp>
+#include <boost/math/ifgt/coefficients.hpp>
+#include <boost/math/ifgt/coefficients_evaluator.hpp>
+#include <libs/math/ifgt/src/example/coefficients_evaluator.h>
+
+void example_coefficients_evaluator(){
+ std::cout << "-> coefficients_evaluator" << std::endl;
+ using namespace boost::math;
+ using namespace boost::math::ifgt;
+ const unsigned int dim = 2;
+ typedef double value_type;
+ typedef std::vector<value_type> var_type;
+ typedef std::vector<unsigned> ints_type;
+ typedef coefficients<dim> cs_type;
+ typedef coefficients_evaluator<var_type> cs_evaluator_type;
+ typedef normal_kernel_properties<dim,value_type> rpp_type;
+
+ double bandwidth = 0.1;
+ int degree = 5;
+ double cutoff_radius = 0.1;
+ var_type center;
+ var_type source;
+ var_type target;
+ value_type weight = 0.1;
+ {
+ using namespace boost::assign;
+ center += (-0.1), 0.1;
+ source += (-0.2), 0.2;
+ target += (-0.15),0.15;
+ }
+
+ var_type zscore_val(dim);
+ monomials<var_type> monoms;
+ {
+ typedef zscore<var_type> zscore_type;
+ zscore_type zscore_op(center);
+ zscore_op(source,bandwidth,zscore_val);
+ }
+ value_type zscore_sqnorm = std::inner_product(
+ boost::begin(zscore_val),
+ boost::end(zscore_val),
+ boost::begin(zscore_val),
+ (value_type)(0)
+ );
+
+ monoms(zscore_val,degree);
+
+ cs_type cs;
+ cs(
+ zscore_sqnorm,
+ monoms(),
+ weight,
+ degree
+ );
+
+ cs_evaluator_type cse(target);
+ value_type inv_nc = rpp_type::normalizing_constant(bandwidth);
+ typedef cs_evaluator_type::result_type result_type;
+ result_type res = cse(
+ cutoff_radius,
+ center,
+ bandwidth,
+ inv_nc,
+ cs
+ );
+
+ std::cout << "res=" << res << std::endl;
+
+ std::cout << "<-" << std::endl;
+}

Added: sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/coefficients_evaluator.h
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/coefficients_evaluator.h 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,10 @@
+//////////////////////////////////////////////////////////////////////////////
+// example/coefficients_evaluator.h
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef LIBS_MATH_IFGT_EXAMPLE_COLLECTED_SOURCES_EVALUATOR_H_ER_2009
+#define LIBS_MATH_IFGT_EXAMPLE_COLLECTED_SOURCES_EVALUATOR_H_ER_2009
+void example_coefficients_evaluator();
+#endif

Added: sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/exact_gauss_transform.cpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/exact_gauss_transform.cpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,113 @@
+///////////////////////////////////////////////////////////////////////////////
+// example/exact_gauss_transform.cpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#include <algorithm>
+#include <boost/assign/std/vector.hpp>
+#include <boost/range.hpp>
+#include <iostream>
+#include <boost/math/ifgt/tag.hpp>
+#include <boost/math/ifgt/exact_contribution.hpp>
+#include <boost/math/ifgt/exact_accumulator.hpp>
+#include <boost/math/ifgt/exact_contributions_evaluator.hpp>
+#include <boost/math/ifgt/exact_evaluator.hpp>
+#include <libs/math/ifgt/src/example/exact_gauss_transform.h>
+
+void example_exact_gauss_transform(){
+
+ std::cout << "-> example_exact_gauss_transform" << std::endl;
+
+ namespace math = boost::math;
+ namespace ifgt = math::ifgt;
+ namespace tag = ifgt::tag;
+
+ const unsigned int dim = 2;
+ const unsigned int wdim = 2;
+
+ typedef double value_type;
+ typedef std::vector<double> var_type;
+ typedef std::vector<unsigned> ints_type;
+ typedef ifgt::exact_contribution<dim,wdim,var_type> exact_contrib_type;
+ typedef std::vector<exact_contrib_type> contributions_type;
+ typedef ifgt::exact_accumulator<dim,wdim,var_type> exact_acc_type;
+ typedef ifgt::exact_evaluator<exact_acc_type> exact_eval_type;
+
+ const double bandwidth = 0.5;
+ var_type center0;
+ var_type center1;
+ var_type center2;
+ var_type source1;
+ var_type source2;
+ var_type weight;
+ var_type target;
+ {
+ using namespace boost::assign;
+ center0 += -1.1, 0.1;
+ center1 += 0.2, -0.2;
+ center2 += 0.1, 0.1;
+ source1 += 0.15,-0.15;//nearest to center1
+ source2 += 0.11, 0.11;//nearest to center2
+ weight += 1.0, 0.1;
+ target += (-0.15),0.15;
+ }
+
+ var_type range_out0(dim,(value_type)(0));
+ var_type range_out1(dim,(value_type)(0));
+
+ exact_acc_type exact_acc(
+ (tag::bandwidth = bandwidth, tag::contributions_reserve_size = 2));
+ exact_acc(center0,weight);
+ exact_acc(center1,weight);
+ exact_acc(center2,weight);
+ exact_acc(source1,weight);
+ exact_acc(source2,weight);
+
+ make_exact_contributions_evaluator(exact_acc.contributions())(
+ target,range_out0);
+
+ std::cout << "gauss transform" << std::endl;
+ std::cout << "(contributions evaluator)" << std::endl;
+ copy(
+ boost::begin(range_out0),
+ boost::end(range_out0),
+ std::ostream_iterator<value_type>(std::cout," ")
+ ); std::cout << std::endl;
+
+ exact_eval_type exact_eval((tag::accumulator = exact_acc));
+ exact_eval.gauss_transform(target,range_out1);
+
+ std::cout << "(exact evaluator)" << std::endl;
+ copy(
+ boost::begin(range_out1),
+ boost::end(range_out1),
+ std::ostream_iterator<value_type>(std::cout," ")
+ ); std::cout << std::endl;
+
+ BOOST_ASSERT(
+ range_out0 == range_out1
+ );
+
+ //-> TODO remove
+ std::cout << "contributions()[0].source()" << std::endl;
+ std::cout << "(exact_acc.)" << std::endl;
+ copy(
+ boost::begin(exact_acc.contributions()[0].source()),
+ boost::end(exact_acc.contributions()[0].source()),
+ std::ostream_iterator<value_type>(std::cout, " ")
+ ); std::cout << std::endl;
+ std::cout << "(exact_eval.accumulator().)" << std::endl;
+ copy(
+ boost::begin(exact_eval.accumulator().contributions()[0].source()),
+ boost::end(exact_eval.accumulator().contributions()[0].source()),
+ std::ostream_iterator<value_type>(std::cout, " ")
+ ); std::cout << std::endl;
+ //<-
+
+ value_type hat_pdf = exact_eval.rozenblatt_parzen_estimate(target);
+ std::cout << "rp estimate" << hat_pdf << std::endl;
+
+ std::cout << "<-" << std::endl;
+
+};

Added: sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/exact_gauss_transform.h
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/exact_gauss_transform.h 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,10 @@
+//////////////////////////////////////////////////////////////////////////////
+// example/exact_gauss_transform.h
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef LIBS_MATH_IFGT_EXAMPLE_EXACT_GAUSS_TRANSFORM_H_ER_2009
+#define LIBS_MATH_IFGT_EXAMPLE_EXACT_GAUSS_TRANSFORM_H_ER_20090
+void example_exact_gauss_transform();
+#endif

Added: sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/fast_accumulator.cpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/fast_accumulator.cpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,120 @@
+///////////////////////////////////////////////////////////////////////////////
+// example/fast_accumulator.cpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#include <algorithm>
+#include <boost/math/special_functions/binomial.hpp>
+#include <boost/assign/std/vector.hpp>
+#include <boost/range.hpp>
+#include <iostream>
+#include <boost/math/ifgt/fast_accumulator.hpp>
+#include <boost/math/ifgt/truncation_degree_constant.hpp>
+#include <libs/math/ifgt/src/example/fast_accumulator.h>
+void example_fast_accumulator(){
+ std::cout << "-> example_fast_accumulator" << std::endl;
+ namespace math = boost::math;
+ namespace ifgt = math::ifgt;
+ namespace tag = ifgt::tag;
+ namespace mpl = boost::mpl;
+
+ const unsigned int dim = 2;
+ const unsigned int wdim = 2;
+ typedef std::vector<double> var_type;
+ typedef ifgt::truncation_degree_constant<mpl::_1> trunc_degree;
+ typedef ifgt::cluster<dim,wdim,trunc_degree,var_type> cluster_type;
+ typedef ifgt::find_nearest_cluster<mpl::_1, boost::l2_distance_squared>
+ find_type;
+
+ typedef ifgt::fast_accumulator<cluster_type,find_type> acc_type;
+
+ double bandwidth = 0.1;
+ double max_cluster_radius = 1.0;
+ unsigned degree = 20;
+ var_type center0;
+ var_type center1;
+ var_type center2;
+ var_type source1;
+ var_type source2;
+ var_type weight;
+ {
+ using namespace boost::assign;
+ center0 += -1.1, 0.1;
+ center1 += 0.2, -0.2;
+ center2 += 0.1, 0.1;
+ source1 += 0.15,-0.15;//nearest to center1
+ source2 += 0.11, 0.11;//nearest to center2
+ weight += 1.0, 0.1;
+ }
+
+ acc_type acc((
+ tag::bandwidth = bandwidth,
+ tag::max_cluster_radius = max_cluster_radius,
+ tag::degree = degree
+ )
+ );
+
+ acc(center0, weight);
+ std::cout
+ << "sources_count = " << acc.sources_count()
+ << std::endl;
+ acc(center1, weight);
+ std::cout
+ << "sources_count = " << acc.sources_count()
+ << std::endl;
+ acc(center2, weight);
+ std::cout
+ << "sources_count = " << acc.sources_count()
+ << std::endl;
+ acc(source1, weight);
+ std::cout
+ << "sources_count = " << acc.sources_count()
+ << std::endl;
+ acc(source2, weight);
+ std::cout
+ << "sources_count = " << acc.sources_count()
+ << std::endl;
+
+// typedef clusters_type::iterator iter_type;
+
+// for(
+// iter_type i = boost::begin(clusters);
+// i<boost::end(clusters);
+// i++
+// ){
+// std::cout
+// << "cluster # "
+// << std::distance(boost::begin(clusters),i)
+// << std::endl;
+
+// {
+// const coeffs_type& sc
+// = i->collection_coefficients()[0];
+
+// std::cout << "1st weight: coeffs=";
+// copy(
+// boost::begin(sc()),
+// boost::end(sc()),
+// std::ostream_iterator<double>(std::cout," ")
+// ); std::cout << std::endl;
+
+// }
+// {
+// const coeffs_type& sc
+// = i->collection_coefficients()[1];
+
+// std::cout << "2nd weight: coeffs=";
+// copy(
+// boost::begin(sc()),
+// boost::end(sc()),
+// std::ostream_iterator<double>(std::cout," ")
+// ); std::cout << std::endl;
+
+// }
+// }//only cluster #1 should have coefficients
+
+
+ std::cout << "<-" << std::endl;
+
+}

Added: sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/fast_accumulator.h
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/fast_accumulator.h 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,10 @@
+//////////////////////////////////////////////////////////////////////////////
+// example/fast_accumulator.h
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef LIBS_MATH_IFGT_EXAMPLE_FAST_ACCUMULATOR_H_ER_200901
+#define LIBS_MATH_IFGT_EXAMPLE_FAST_ACCUMULATOR_H_ER_200901
+void example_fast_accumulator();
+#endif

Added: sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/fast_evaluator.cpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/fast_evaluator.cpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,102 @@
+///////////////////////////////////////////////////////////////////////////////
+// example/fast_evaluator.cpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#include <boost/assign/std/vector.hpp>
+#include <boost/range.hpp>
+#include <iostream>
+#include <boost/math/ifgt/fast_accumulator.hpp>
+#include <boost/math/ifgt/truncation_degree_constant.hpp>
+#include <boost/math/ifgt/cutoff_radius_none.hpp>
+#include <boost/math/ifgt/fast_evaluator.hpp>
+#include <libs/math/ifgt/src/example/fast_evaluator.h>
+
+
+void example_fast_evaluator(){
+
+ std::cout << "-> example_fast_evaluator" << std::endl;
+ namespace math = boost::math;
+ namespace ifgt = math::ifgt;
+ namespace tag = ifgt::tag;
+ namespace mpl = boost::mpl;
+
+ const unsigned int dim = 2;
+ const unsigned int wdim = 2;
+ typedef double value_type;
+ typedef std::vector<value_type> var_type;
+ typedef ifgt::truncation_degree_constant<mpl::_1> trunc_degree;
+ typedef ifgt::cluster<dim,wdim,trunc_degree,var_type> cluster_type;
+ typedef ifgt::find_nearest_cluster<mpl::_1, boost::l2_distance_squared>
+ find_type;
+ typedef ifgt::fast_accumulator<cluster_type,find_type> acc_type;
+
+ value_type bandwidth = 0.1;
+ value_type max_cluster_radius = 1.0;
+ unsigned degree = 20;
+ var_type center0;
+ var_type center1;
+ var_type center2;
+ var_type source1;
+ var_type source2;
+ var_type weight;
+ var_type target;
+ {
+ using namespace boost::assign;
+ center0 += -1.1, 0.1;
+ center1 += 0.2, -0.2;
+ center2 += 0.1, 0.1;
+ source1 += 0.15,-0.15;//nearest to center1
+ source2 += 0.11, 0.11;//nearest to center2
+ weight += 1.0, 0.1;
+ target += (-0.15),0.15;
+ }
+
+ acc_type acc((
+ tag::bandwidth = bandwidth,
+ tag::max_cluster_radius = max_cluster_radius,
+ tag::degree = degree
+ )
+ );
+
+ acc(center0, weight);
+ acc(center1, weight);
+ acc(center2, weight);
+ acc(source1, weight);
+ acc(source2, weight);
+ std::cout
+ << "sources_count = " << acc.sources_count()
+ << std::endl;
+
+ typedef ifgt::fast_evaluator<acc_type,
+ ifgt::cutoff_radius_none<mpl::_1> > eval_type;
+ eval_type eval((tag::accumulator = acc));
+
+ var_type range_out(wdim);
+
+ eval.gauss_transform(target,range_out);
+
+ std::cout << "gauss_transform =" << std::endl;
+ std::copy(
+ boost::begin(range_out),
+ boost::end(range_out),
+ std::ostream_iterator<value_type>(std::cout, " ")
+ );
+
+ value_type rpe = eval.rozenblatt_parzen_estimate(target);
+ std::cout << "rozenblatt_parzen_estimate =" << rpe << std::endl;
+
+ var_type range_out_nw(wdim-1);
+ eval.nadaraya_watson_estimate(target,range_out_nw); //wrong size
+
+ std::cout << "nadaya_watson_estimate =" << std::endl;
+ std::copy(
+ boost::begin(range_out_nw),
+ boost::end(range_out_nw),
+ std::ostream_iterator<value_type>(std::cout, " ")
+ );
+
+ std::cout << "<-" << std::endl;
+
+};

Added: sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/fast_evaluator.h
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/fast_evaluator.h 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,10 @@
+//////////////////////////////////////////////////////////////////////////////
+// example/fast_evaluator.h
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef LIBS_MATH_IFGT_FAST_EVALUATOR_H_ER2008
+#define LIBS_MATH_IFGT_FAST_EVALUATOR_H_ER2008
+void example_fast_evaluator();
+#endif

Added: sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/find_nearest_cluster.cpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/find_nearest_cluster.cpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,86 @@
+///////////////////////////////////////////////////////////////////////////////
+// example/find_nearest_cluster.cpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#include <algorithm>
+#include <vector>
+#include <boost/assign/std/vector.hpp>
+#include <boost/range.hpp>
+#include <iostream>
+#include <boost/assign/std/vector.hpp>
+#include <boost/range.hpp>
+#include <boost/math/ifgt/cluster.hpp>
+#include <boost/math/ifgt/find_nearest_cluster.hpp>
+#include <libs/math/ifgt/src/example/find_nearest_cluster.h>
+
+void example_find_nearest_cluster(){
+ std::cout << "-> find_nearest_cluster" << std::endl;
+
+ namespace mpl = boost::mpl;
+ namespace math = boost::math;
+ namespace ifgt = math::ifgt;
+ namespace tag = ifgt::tag;
+
+ const unsigned dim = 2;
+ const unsigned wdim = 2;
+ typedef double value_type;
+ typedef std::vector<value_type> var_type;
+ typedef ifgt::truncation_degree_constant<mpl::_1> trunc_degree;
+ typedef ifgt::cluster<dim,wdim,trunc_degree,var_type> cluster_type;
+ typedef std::vector<cluster_type> clusters_type;
+ typedef boost::l2_distance_squared dist_type;
+ typedef ifgt::find_nearest_cluster<clusters_type,dist_type> find_nc_type;
+
+ typedef cluster_type::coefficients_type coeffs_type;
+ typedef cluster_type::collection_coefficients_type coll_coeffs_type;
+
+ const value_type bandwidth = 0.1;
+ const value_type max_cluster_radius = 0.2;
+ const unsigned degree = 20;
+ var_type center0;
+ var_type center1;
+ var_type center2;
+ var_type source;
+ {
+ using namespace boost::assign;
+ center0 += -1.1, 0.1;
+ center1 += 0.2, -0.2;
+ center2 += 0.1, 0.1;
+ source += 0.15,-0.15;
+ }
+
+ typedef mpl::apply<trunc_degree,value_type>::type trunc_degree_type;
+
+ clusters_type clusters;
+ clusters.push_back(cluster_type(
+ center0,
+ bandwidth,
+ max_cluster_radius,
+ trunc_degree_type(degree)));
+ clusters.push_back(
+ cluster_type((
+ tag::center = center1,
+ tag::bandwidth = bandwidth,
+ tag::max_cluster_radius = max_cluster_radius,
+ tag::degree = degree
+ )));
+ clusters.push_back(
+ cluster_type((
+ tag::center = center2,
+ tag::bandwidth = bandwidth,
+ tag::max_cluster_radius = max_cluster_radius,
+ tag::degree = degree
+ )));
+
+ find_nc_type find_nc(clusters);
+ typedef find_nc_type::iter_type iter_type;
+ iter_type nc = find_nc(source);
+ std::cout
+ << "dist to nc = "
+ << std::distance(boost::begin(clusters),nc)
+ << std::endl;
+
+ std::cout << "<-" << std::endl;
+}

Added: sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/find_nearest_cluster.h
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/find_nearest_cluster.h 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,10 @@
+//////////////////////////////////////////////////////////////////////////////
+// example/find_nearest_cluster.h
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef LIBS_MATH_IFGT_EXAMPLE_FIND_NEAREST_CLUSTER_H_ER_2009
+#define LIBS_MATH_IFGT_EXAMPLE_FIND_NEAREST_CLUSTER_H_ER_2009
+void example_find_nearest_cluster();
+#endif

Added: sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/for_each_accumulate.cpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/for_each_accumulate.cpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,109 @@
+///////////////////////////////////////////////////////////////////////////////
+// example/for_each_accumulate.cpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#include <algorithm>
+#include <iostream>
+#include <boost/assign/std/vector.hpp>
+#include <boost/range.hpp>
+#include <boost/math/ifgt/exact_accumulator.hpp>
+#include <boost/math/ifgt/fast_accumulator.hpp>
+#include <boost/math/ifgt/truncation_degree_constant.hpp>
+#include <boost/math/ifgt/for_each_accumulate.hpp>
+#include <boost/math/ifgt/for_each_accumulate_rest_weights.hpp>
+#include <libs/math/ifgt/src/example/for_each_accumulate.h>
+
+void example_for_each_accumulate(){
+ std::cout << "-> example_for_each_fast_accumulator" << std::endl;
+ namespace math = boost::math;
+ namespace ifgt = math::ifgt;
+ namespace tag = ifgt::tag;
+ namespace mpl = boost::mpl;
+
+ const unsigned int dim = 2;
+ const unsigned int wdim = 2;
+ typedef std::vector<double> var_type;
+ typedef ifgt::truncation_degree_constant<mpl::_1> trunc_degree;
+ typedef ifgt::cluster<dim,wdim,trunc_degree,var_type> cluster_type;
+ typedef ifgt::find_nearest_cluster<mpl::_1, boost::l2_distance_squared>
+ find_type;
+
+ typedef ifgt::fast_accumulator<cluster_type,find_type> fast_acc_type;
+ typedef ifgt::exact_accumulator<dim,wdim,var_type> exact_acc_type;
+
+
+ double bandwidth = 0.1;
+ double max_cluster_radius = 1.0;
+ unsigned degree = 20;
+ var_type sources;
+ var_type weights;
+ var_type r_weights;
+ {
+ using namespace boost::assign;
+ sources += -1.1, 0.1;
+ weights += 1.0, 0.1;
+ r_weights += 0.1;
+
+ sources += 0.2, -0.2;
+ weights += 1.0, 0.1;
+ r_weights += 0.1;
+
+ sources += 0.1, 0.1;
+ weights += 1.0, 0.1;
+ r_weights += 0.1;
+
+ sources += 0.15,-0.15;
+ weights += 1.0, 0.1;
+ r_weights += 0.1;
+
+ sources += 0.11, 0.11;
+ weights += 1.0, 0.1;
+ r_weights += 0.1;
+ }
+
+ std::cout << "sources count" << std::endl;
+
+ fast_acc_type fast_acc((
+ tag::bandwidth = bandwidth,
+ tag::max_cluster_radius = max_cluster_radius,
+ tag::degree = degree
+ )
+ );
+ ifgt::for_each(sources,weights,fast_acc);
+ std::cout
+ << "fast (weights) " << fast_acc.sources_count()
+ << std::endl;
+
+ fast_acc_type fast_acc2((
+ tag::bandwidth = bandwidth,
+ tag::max_cluster_radius = max_cluster_radius,
+ tag::degree = degree
+ )
+ );
+ ifgt::for_each(sources,
+ ifgt::make_rest_weights_wrapper(r_weights),fast_acc2);
+ std::cout
+ << "fast (rest weights) " << fast_acc2.sources_count()
+ << std::endl;
+
+ exact_acc_type exact_acc((tag::bandwidth = bandwidth));
+ ifgt::for_each(sources,weights,exact_acc);
+
+ std::cout
+ << "exact (weights) " << exact_acc.sources_count()
+ << std::endl;
+
+ exact_acc_type exact_racc((tag::bandwidth = bandwidth));
+ ifgt::for_each(sources,
+ ifgt::make_rest_weights_wrapper(r_weights),exact_racc);
+
+ std::cout
+ << "exact (rest weights) " << exact_racc.sources_count()
+ << std::endl;
+
+
+ std::cout << "<-" << std::endl;
+
+};

Added: sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/for_each_accumulate.h
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/for_each_accumulate.h 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,10 @@
+//////////////////////////////////////////////////////////////////////////////
+// example/for_each_accumulate.h
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef LIBS_MATH_IFGT_EXAMPLE_FOR_EACH_ACCUMULATE_H_ER_2009
+#define LIBS_MATH_IFGT_EXAMPLE_FOR_EACH_ACCUMULATE_H_ER_2009
+void example_for_each_accumulate();
+#endif

Added: sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/for_each_evaluate.cpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/for_each_evaluate.cpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,141 @@
+///////////////////////////////////////////////////////////////////////////////
+// example/for_each_evaluate.cpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#include <algorithm>
+#include <iostream>
+#include <boost/assign/std/vector.hpp>
+#include <boost/range.hpp>
+#include <boost/math/ifgt/exact_accumulator.hpp>
+#include <boost/math/ifgt/fast_accumulator.hpp>
+#include <boost/math/ifgt/fast_evaluator.hpp>
+#include <boost/math/ifgt/truncation_degree_constant.hpp>
+#include <boost/math/ifgt/for_each_accumulate_rest_weights.hpp>
+#include <boost/math/ifgt/for_each_evaluate.hpp>
+#include <boost/math/ifgt/call_wrapper.hpp>
+#include <boost/math/ifgt/call_gauss_transform.hpp>
+#include <boost/math/ifgt/call_rozenblatt_parzen_estimate.hpp>
+#include <boost/math/ifgt/call_nadaraya_watson_estimate.hpp>
+#include <boost/math/ifgt/for_each_rozenblatt_parzen_estimate.hpp>
+#include <boost/math/ifgt/for_each_evaluate.hpp>
+#include <libs/math/ifgt/src/example/for_each_evaluate.h>
+
+void example_for_each_evaluate(){
+ std::cout << "-> example_for_each_evaluate" << std::endl;
+
+ namespace math = boost::math;
+ namespace ifgt = math::ifgt;
+ namespace tag = ifgt::tag;
+ namespace mpl = boost::mpl;
+
+ const unsigned int dim = 2;
+ const unsigned int wdim = 2;
+ typedef double value_type;
+ typedef std::vector<value_type> var_type;
+ typedef ifgt::truncation_degree_constant<mpl::_1> trunc_degree;
+ typedef ifgt::cluster<dim,wdim,trunc_degree,var_type> cluster_type;
+ typedef ifgt::find_nearest_cluster<mpl::_1, boost::l2_distance_squared>
+ find_type;
+
+ typedef ifgt::fast_accumulator<cluster_type,find_type> fast_acc_type;
+ typedef ifgt::exact_accumulator<dim,wdim,var_type> exact_acc_type;
+
+
+ value_type bandwidth = 0.1;
+ value_type max_cluster_radius = 1.0;
+ unsigned degree = 20;
+ var_type sources;
+ var_type weights;
+ {
+ using namespace boost::assign;
+ sources += -1.1, 0.1;
+ weights += 1.0, 0.1;
+ sources += 0.2, -0.2;
+ weights += 1.0, 0.1;
+ sources += 0.1, 0.1;
+ weights += 1.0, 0.1;
+ sources += 0.15,-0.15;
+ weights += 1.0, 0.1;
+ sources += 0.11, 0.11;
+ weights += 1.0, 0.1;
+ }
+
+ fast_acc_type fast_acc((
+ tag::bandwidth = bandwidth,
+ tag::max_cluster_radius = max_cluster_radius,
+ tag::degree = degree
+ )
+ );
+
+ ifgt::for_each(sources,weights,fast_acc);
+
+ exact_acc_type exact_acc((tag::bandwidth = bandwidth));
+ ifgt::for_each(sources,weights,exact_acc);
+
+ var_type targets;
+ {
+ using namespace boost::assign;
+ targets += (-0.15),0.15;//1
+ targets += (-0.15),0.15;//2
+ targets += (-0.15),0.15;//3
+ targets += (-0.15),0.15;//4
+ targets += (-0.15),0.15;//5
+ }
+ var_type::size_type targets_count = boost::size(targets)/dim;
+
+
+ typedef ifgt::cutoff_radius_none<mpl::_1> cutoff_policy_type;
+ typedef ifgt::fast_evaluator<fast_acc_type, cutoff_policy_type>
+ fast_eval_type;
+ fast_eval_type eval((tag::accumulator = fast_acc));
+
+ var_type ranges_out( targets_count * wdim);
+ std::cout << "gauss_transform->" << std::endl;
+ ifgt::for_each(
+ targets,
+ ranges_out,
+ eval,
+ ifgt::call<ifgt::gauss_transform>()
+ );
+ copy(
+ boost::begin(ranges_out),
+ boost::end(ranges_out),
+ std::ostream_iterator<value_type>(std::cout," ")
+ );
+ std::cout << "<-" << std::endl;
+ std::cout << "rozenblatt_parzen_estimate->" << std::endl;
+
+ var_type rp_out( targets_count * 1);
+ ifgt::for_each(
+ targets,
+ rp_out,
+ eval,
+ ifgt::call<ifgt::rozenblatt_parzen_estimate>()
+ );
+ copy(
+ boost::begin(rp_out),
+ boost::end(rp_out),
+ std::ostream_iterator<value_type>(std::cout," ")
+ );
+ std::cout << "<-" << std::endl;
+
+ var_type nw_out( targets_count * (wdim-1) );
+ std::cout << "nadaraya_watston_estimate->" << std::endl;
+ ifgt::for_each(
+ targets,
+ nw_out,
+ eval,
+ ifgt::call<ifgt::nadaraya_watson_estimate>()
+ );
+ copy(
+ boost::begin(nw_out),
+ boost::end(nw_out),
+ std::ostream_iterator<value_type>(std::cout," ")
+ );
+
+ std::cout << "<-" << std::endl;
+
+ std::cout << "<-" << std::endl;
+};

Added: sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/for_each_evaluate.h
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/for_each_evaluate.h 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,10 @@
+//////////////////////////////////////////////////////////////////////////////
+// example/for_each_evaluate.h
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef LIBS_MATH_EXAMPLE_FOR_EACH_EVALUATE_H_ER_2009
+#define LIBS_MATH_EXAMPLE_FOR_EACH_EVALUATE_H_ER_2009
+void example_for_each_evaluate();
+#endif

Added: sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/main.cpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/main.cpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,42 @@
+//////////////////////////////////////////////////////////////////////////////
+// example/main.cpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#include <libs/math/ifgt/src/example/truncation_properties.h>
+#include <libs/math/ifgt/src/example/cluster.h>
+#include <libs/math/ifgt/src/example/find_nearest_cluster.h>
+#include <libs/math/ifgt/src/example/truncation_degree_constant.h>
+#include <libs/math/ifgt/src/example/coefficients.h>
+#include <libs/math/ifgt/src/example/coefficients_evaluator.h>
+#include <libs/math/ifgt/src/example/cluster_evaluator.h>
+#include <libs/math/ifgt/src/example/clusters_evaluator.h>
+#include <libs/math/ifgt/src/example/fast_accumulator.h>
+#include <libs/math/ifgt/src/example/fast_evaluator.h>
+#include <libs/math/ifgt/src/example/exact_gauss_transform.h>
+#include <libs/math/ifgt/src/example/for_each_accumulate.h>
+#include <libs/math/ifgt/src/example/for_each_evaluate.h>
+#include <libs/math/ifgt/src/example/benchmark_exact.h>
+#include <libs/math/ifgt/src/example/benchmark_fast.h>
+
+int main()
+{
+
+ example_truncation_properties();
+ example_coefficients();
+ example_coefficients_evaluator();
+ example_cluster();
+ example_find_nearest_cluster();
+ example_truncation_degree_constant();
+ example_cluster_evaluator();
+ example_clusters_evaluator();
+ example_exact_gauss_transform();
+ example_fast_accumulator();
+ example_fast_evaluator();
+ example_for_each_accumulate();
+ example_for_each_evaluate();
+ example_benchmark_exact();
+ example_benchmark_fast();
+ return 0;
+}

Added: sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/truncation_degree_constant.cpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/truncation_degree_constant.cpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,54 @@
+///////////////////////////////////////////////////////////////////////////////
+// example/truncation_degree_constant.cpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#include <algorithm>
+#include <vector>
+#include <iostream>
+#include <boost/assign/std/vector.hpp>
+#include <boost/range.hpp>
+#include <boost/math/ifgt/truncation_degree_constant.hpp>
+#include <boost/math/ifgt/tag.hpp>
+#include <libs/math/ifgt/src/example/truncation_degree_constant.h>
+void example_truncation_degree_constant(){
+ std::cout << "-> example_truncation_degree_constant" << std::endl;
+
+ using namespace boost::math::ifgt;
+
+ const unsigned degree = 5;
+ const unsigned dim = 2;
+
+ typedef truncation_degree_constant<> truncation_degree_type;
+ typedef std::vector<unsigned> degrees_type;
+ typedef double value_type;
+ typedef std::vector<value_type> var_type;
+
+ value_type bandwidth = 0.1;
+ value_type sqrt_dist_source_to_center = 0.1;
+
+ truncation_degree_type a((tag::degree = degree));
+
+ var_type weight;
+ {
+ using namespace boost::assign;
+ weight += 1.0, 0.1;
+ }
+ degrees_type degrees(dim);
+
+ a(
+ bandwidth,
+ sqrt_dist_source_to_center,
+ weight,
+ degrees
+ );
+
+ copy(
+ boost::begin(degrees),
+ boost::end(degrees),
+ std::ostream_iterator<unsigned>(std::cout, " ")
+ );//should return degree, degree
+
+ std::cout << "<-" << std::endl;
+}

Added: sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/truncation_degree_constant.h
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/truncation_degree_constant.h 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,10 @@
+//////////////////////////////////////////////////////////////////////////////
+// example/truncation_degree_constant.h
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef LIBS_MATH_IFGT_EXAMPLE_ASSIGN_DEGREES_CONSTANT_H_ER_2009
+#define LIBS_MATH_IFGT_EXAMPLE_ASSIGN_DEGREES_CONSTANT_H_ER_2009
+void example_truncation_degree_constant();
+#endif

Added: sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/truncation_properties.cpp
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/truncation_properties.cpp 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,77 @@
+///////////////////////////////////////////////////////////////////////////////
+// example/truncation_properties.cpp
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#include <iostream>
+#include <iterator>
+#include <boost/range.hpp>
+#include <boost/math/ifgt/include.hpp>
+#include <libs/math/ifgt/src/example/truncation_properties.h>
+
+void example_truncation_properties(){
+ std::cout << "->example_truncation_properties" << std::endl;
+ using namespace boost;
+ namespace math = boost::math;
+ namespace ifgt = math::ifgt;
+ typedef double value_type;
+ typedef std::vector<value_type> vec_type;
+ typedef ifgt::truncation_properties prop_type;
+ typedef ifgt::optimal_parameter_given_max_degree<value_type>
+ opt_pars_type;
+
+ value_type h = 0.5; //Raykar2006 probably meant h=0.5 in Figure 6
+ value_type rx = 0.5;
+ value_type ry = 1e15;
+ value_type eps = 1e-3;
+// unsigned max_p = 40;
+ unsigned n = 100;
+ opt_pars_type opt_pars((ifgt::tag::error_tolerance = eps));
+
+ value_type log10 = log((value_type)(10));
+ vec_type wry;
+ vec_type log10_err;
+ vec_type log10_err_bound;
+// vec_type opt_max_radius;
+ value_type value;
+ for(unsigned i = 0; i<n; i++){
+ value_type ryi = 5 * (value_type)(i+1)/((value_type)(n));
+ value = prop_type::log_error_bound(15, h, rx, ryi)/log10;
+ log10_err.push_back(value);
+ value_type rxi = (value_type)(i+1)/((value_type)(n));
+ value = prop_type::worst_target_to_center_distance(40, h, rxi);
+ wry.push_back(value);
+ value = prop_type::max_log_error_bound(40, h, rxi, ry)/log10;
+ log10_err_bound.push_back(value);
+
+ };
+
+ std::cout << " ->log10_errs - Figure 5a in Raykar2006a" << std::endl;
+ std::copy(
+ begin(log10_err),
+ end(log10_err),
+ std::ostream_iterator<value_type>(std::cout," ")
+ );
+ std::cout << " <-" << std::endl;
+
+ std::cout << " -> worst ry" << std::endl;
+ std::copy(
+ begin(wry),
+ end(wry),
+ std::ostream_iterator<value_type>(std::cout," ")
+ );
+ std::cout << " <-" << std::endl;
+
+ std::cout
+ << " ->log10_err_bounds - Figure 6a in Raykar2006a"
+ << std::endl;
+ std::copy(
+ begin(log10_err_bound),
+ end(log10_err_bound),
+ std::ostream_iterator<value_type>(std::cout," ")
+ );
+ std::cout << " <-" << std::endl;
+
+ std::cout << "<-" << std::endl;
+};

Added: sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/truncation_properties.h
==============================================================================
--- (empty file)
+++ sandbox/improved_fast_gauss_transform/libs/math/ifgt/src/example/truncation_properties.h 2009-02-10 04:15:13 EST (Tue, 10 Feb 2009)
@@ -0,0 +1,10 @@
+//////////////////////////////////////////////////////////////////////////////
+// example/truncation_properties.h
+// (C) Copyright 2009 Erwann Rogard
+// Use, modification and distribution are 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)
+#ifndef LIBS_MATH_IFGT_TRUNCATION_PROPERTIES_H_ER_2009
+#define LIBS_MATH_IFGT_TRUNCATION_PROPERTIES_H_ER_2009
+void example_truncation_properties();
+#endif

Next message: pbristow_at_[hidden]: "[Boost-commit] svn:boost r51175 - sandbox/SOC/2007/visualization/libs/svg_plot/test"
Previous message: erwann.rogard_at_[hidden]: "[Boost-commit] svn:boost r51173 - sandbox/miscellanea_algorithm/boost/algorithm"

Date view	Thread view	Subject view	Author view

Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk