Boost-Commit :

Date view	Thread view	Subject view	Author view

Subject: [Boost-commit] svn:boost r52002 - in sandbox/conditionally_specified_distribution: conditionally_specified_distribution/boost/conditionally_specified_distribution/function conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter/dependent conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter/dependent/detail conditionally_specified_distribution/boost/conditionally_specified_distribution/sampler/adaptive_rejection conditionally_specified_distribution/boost/conditionally_specified_distribution/support conditionally_specified_distribution/boost/conditionally_specified_distribution/support/detail conditionally_specified_distribution/libs/conditionally_specified_distribution/src conditionally_specified_distribution/libs/conditionally_specified_distribution/src/doc conditionally_specified_distribution/libs/conditionally_specified_distribution/src/example/sampler math_limit/boost/math_limit
From: erwann.rogard_at_[hidden]
Date: 2009-03-26 18:24:28

Author: e_r
Date: 2009-03-26 18:24:23 EDT (Thu, 26 Mar 2009)
New Revision: 52002
URL: http://svn.boost.org/trac/boost/changeset/52002

Log:
Fixed bug support(), added as_string() and misc in libs
Added:
   sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/support/pair_independent.hpp (contents, props changed)
Text files modified:
   sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/function/normal.hpp | 1
   sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter/dependent/detail/regression_coefficient.hpp | 28 ++++++-
   sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter/dependent/normal_given_normal.hpp | 20 +++++
   sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter/normal.hpp | 18 ++++
   sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter/pair_independent.hpp | 16 ++++
   sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/sampler/adaptive_rejection/sampler.hpp | 57 ++++++++++++----
   sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/support/detail/negative_max_to_max.hpp | 8 -
   sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/support/include.hpp | 5 +
   sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/support/normal.hpp | 9 ++
   sandbox/conditionally_specified_distribution/conditionally_specified_distribution/libs/conditionally_specified_distribution/src/doc/readme.txt | 31 +++++++-
   sandbox/conditionally_specified_distribution/conditionally_specified_distribution/libs/conditionally_specified_distribution/src/example/sampler/gibbs_ars_bivariate_normal.cpp | 64 +++++++++++++----
   sandbox/conditionally_specified_distribution/conditionally_specified_distribution/libs/conditionally_specified_distribution/src/main.cpp | 2
   sandbox/conditionally_specified_distribution/math_limit/boost/math_limit/miscellanea.hpp | 139 +++++++++++++++++++++++++++++----------
   13 files changed, 307 insertions(+), 91 deletions(-)

Modified: sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/function/normal.hpp
==============================================================================
--- sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/function/normal.hpp (original)
+++ sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/function/normal.hpp 2009-03-26 18:24:23 EDT (Thu, 26 Mar 2009)
@@ -9,6 +9,7 @@
#define BOOST_CONDITIONALLY_SPECIFIED_DISTRIBUTION_FUNCTION_NORMAL_HPP_ER_2009
#include <cmath>
#include <boost/math/special_functions/fpclassify.hpp>
+#include <boost/math/constants/constants.hpp>
namespace boost{
namespace conditionally_specified_distribution{
namespace function{

Modified: sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter/dependent/detail/regression_coefficient.hpp
==============================================================================
--- sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter/dependent/detail/regression_coefficient.hpp (original)
+++ sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter/dependent/detail/regression_coefficient.hpp 2009-03-26 18:24:23 EDT (Thu, 26 Mar 2009)
@@ -8,8 +8,10 @@
#ifndef BOOST_CONDITIONALLY_SPECIFIED_DISTRIBUTION_PARAMETER_DEPENDENT_REGRESSION_COEFFICIENT_HPP_ER_2009
#define BOOST_CONDITIONALLY_SPECIFIED_DISTRIBUTION_PARAMETER_DEPENDENT_REGRESSION_COEFFICIENT_HPP_ER_2009
#include <vector>
+#include <string>
#include <boost/range.hpp>
#include <boost/bind.hpp>
+#include <boost/format.hpp>
#include <boost/type_traits.hpp>
#include <boost/mpl/assert.hpp>
#include <boost/mpl/identity.hpp>
@@ -18,6 +20,7 @@
#include <boost/utility/dont_care.hpp>
#include <boost/utility/remove_qualifier.hpp>
#include <boost/utility/assert_is_base_of.hpp>
+#include <boost/utility/container_to_string.hpp>
#include <boost/shared_features/feature/keyword.hpp>
#include <boost/shared_features/parameter.hpp>
#include <boost/shared_features/depends_on.hpp>
@@ -86,17 +89,14 @@
>::type
>::type value_type;
             typedef std::vector<value_type> container_type;
-
             public:
             struct result_of_range_skip_position_covariate_dot_coefficient
              : mpl::identity<const container_type&>
             {};
-
             // outside interface
             struct function_value : mpl::identity<value_type>{};
             struct function_argument : mpl::identity<value_type>{};
         };
-
     template<typename Id, typename Dataset,typename Coefficients>
     class regression_coefficient
     : public regression_coefficient_base<
@@ -115,8 +115,6 @@
         :position_(pos){
                 //Warning : state_ default initializes
         }
-
-
         template<typename Args>
         regression_coefficient(const Args& args)
         :position_(args[shared_features::kwd<Id>::position|0]){
@@ -188,6 +186,8 @@
             return this->matching_covariate(u);
         }

+ //TODO range_matching_covariate
+
         // {<x[-i],b[j]> : i =1,...,n}
         typename conditionally_specified_distribution::result_of
             ::range_skip_position_covariate_dot_coefficient<
@@ -222,6 +222,24 @@
                 the_op
             );
         }
+
+ template<typename Args>
+ std::string
+ as_string(const Args& args)const{
+ typedef utility::container_to_string to_str_t;
+ std::string str = "[position = %1%\n";
+ format f(str);
+ f%(this->position(args));
+ str = f.str();
+ str += "range_coefficient =";
+ str += to_str_t::get_indexed(this->range_coefficient(args));
+ str += "\nrange_skip_position_covariate_dot_coefficient = ";
+ str += to_str_t::get_indexed(
+ this->range_skip_position_covariate_dot_coefficient(args)
+ );
+ str +="]";
+ return str;
+ }
         private:
         void alert_if(mpl::bool_<true>,int i)const{
             BOOST_ASSERT(i>-1);

Modified: sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter/normal.hpp
==============================================================================
--- sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter/normal.hpp (original)
+++ sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter/normal.hpp 2009-03-26 18:24:23 EDT (Thu, 26 Mar 2009)
@@ -7,9 +7,12 @@
////////////////////////////////////////////////////////////////////////////
#ifndef BOOST_CONDITIONALLY_SPECIFIED_DISTRIBUTION_PARAMETER_NORMAL_HPP_ER_2009
#define BOOST_CONDITIONALLY_SPECIFIED_DISTRIBUTION_PARAMETER_NORMAL_HPP_ER_2009
+#include <string>
#include <boost/mpl/apply.hpp>
#include <boost/parameter/keyword.hpp>
+#include <boost/format.hpp>
#include <boost/utility/assert_is_base_of.hpp>
+#include <boost/utility/dont_care.hpp>
#include <boost/conditionally_specified_distribution/keyword/parameter.hpp>
#include <boost/conditionally_specified_distribution/result_of/include.hpp>
#include <boost/conditionally_specified_distribution/crtp/include.hpp>
@@ -76,8 +79,7 @@
             return *this;
         }

- template<typename Args>
- void update(const Args& args){}
+ void update(utility::dont_care){}

         typename result_of::mu<super_t>::type
         mu() const {
@@ -97,6 +99,18 @@
         typename result_of::sigma<super_t>::type
         sigma(const Args& args) const { return this->sigma(); }

+ std::string as_string()const{
+ std::string str = "parameter::impl::normal : ";
+ str += " mu = %1%";
+ str += " sigma = %2%";
+ format f(str); f%mu_%sigma_;
+ return f.str();
+ }
+
+ std::string as_string(utility::dont_care)const{
+ return this->as_string();
+ }
+
         private:
         typename super_t::value_type mu_;
         typename super_t::value_type sigma_;

Modified: sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter/pair_independent.hpp
==============================================================================
--- sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter/pair_independent.hpp (original)
+++ sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/parameter/pair_independent.hpp 2009-03-26 18:24:23 EDT (Thu, 26 Mar 2009)
@@ -7,6 +7,7 @@
////////////////////////////////////////////////////////////////////////////
#ifndef BOOST_CONDITIONALLY_SPECIFIED_DISTRIBUTION_PAIR_INDEPENDENT_HPP_ER_2009
#define BOOST_CONDITIONALLY_SPECIFIED_DISTRIBUTION_PAIR_INDEPENDENT_HPP_ER_2009
+#include <string>
#include <boost/mpl/apply.hpp>
#include <boost/mpl/vector.hpp>
#include <boost/parameter/keyword.hpp>
@@ -85,7 +86,7 @@
             return *this;
         }

- template<typename Set,typename Args>
+ template<typename Args>
         void update(const Args& args){
             this->update_impl(args[shared_features::kwd_set],args);
         }
@@ -101,6 +102,19 @@
         second(const Args& args) const {
             return second_impl(args[shared_features::kwd_set],args);
         }
+
+ template<typename Args>
+ std::string as_string(const Args& args)const{
+ std::string str = "parameter::pair_independent :";
+ str+= "\n first : ";
+ str+= first_impl(
+ args[shared_features::kwd_set],args).as_string(args);
+ str+= "\n second : ";
+ str+= second_impl(
+ args[shared_features::kwd_set],args).as_string(args);
+ return str;
+ }
+
         private:
         template<typename Set,typename Args>
         typename result_of::first<super_t(const Args&)>::type

Modified: sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/sampler/adaptive_rejection/sampler.hpp
==============================================================================
--- sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/sampler/adaptive_rejection/sampler.hpp (original)
+++ sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/sampler/adaptive_rejection/sampler.hpp 2009-03-26 18:24:23 EDT (Thu, 26 Mar 2009)
@@ -35,6 +35,19 @@

namespace impl{

+ template<template<typename,typename> class Cont>
+ struct apply_reserve{
+ template<typename Impl,typename T>
+ void operator()(Impl& impl, T n){}
+ };
+ template<>
+ struct apply_reserve<std::vector>{
+ template<typename Impl,typename T>
+ void operator()(Impl& impl, T n){
+ impl.reserve(n);
+ }
+ };
+
     template<typename Par>
     class adaptive_rejection_sampler_base{
         public:
@@ -67,10 +80,13 @@
         typedef kwd<> kwd;

         std::size_t static default_max_init_recursion(){
- std::size_t static result = 100;
+ static std::size_t result = 100;
+ return result;
+ }
+ std::size_t static default_reserve(){
+ static std::size_t result = 10;
             return result;
         }
-
         public:

         template<typename Args>
@@ -78,11 +94,13 @@
         : init_(args),
         max_init_recursion_(
          args[kwd::ars_max_init_recursion|default_max_init_recursion()]
- ){}
+ ),
+ reserve_(default_reserve()){}

         adaptive_rejection_sampler(const adaptive_rejection_sampler& that)
         : init_(that.init_),
- max_init_recursion_(that.max_init_recursion_){}
+ max_init_recursion_(that.max_init_recursion_),
+ reserve_(that.reserve_){}

         //TODO check
         adaptive_rejection_sampler&
@@ -90,6 +108,7 @@
             if(&that!=this){
                 init_ = that.init_;
                 max_init_recursion_ = that.max_init_recursion_;
+ reserve_ = that.reserve_;
             }
             return *this;
         }
@@ -127,11 +146,6 @@
             typedef typename
                 Set::template result_of_extract<State>::type extr_state_t;
             typedef function::detail::adapter<extr_par_t,Args> adapter_t;
-// typedef ::boost::random::adaptive_rejection_sampler::simulator<
-// adapter_t,
-// Container,
-// Allocator
-// > impl_t;
             typedef ::boost::adaptive_rejection_sampling::sampler<
                 adapter_t,
                 Container,
@@ -139,12 +153,24 @@
> impl_t;
             const extr_par_t& par = set.template extract<Par>();
             adapter_t adapter (par,args);
+ //TODO (*)
+ //static impl_t impl(
+ // adapter,
+ // this->max_init_recursion_
+ //);
+
             impl_t impl(
                 adapter,
                 this->max_init_recursion_
             );
+
             interval_t interval = (this->init_)();
- // Have had px!=0 fails here with gibbs rte::regression_coeff
+ // TODO (*)
+ // This causes problem due I think to dependence adapter
+ // on Args (wrap around an abstract base class?)
+ // impl.reset_distribution_function(adapter);
+ typedef apply_reserve<Container> apply_res_t;
+ apply_res_t()(impl,reserve_);
             impl.initialize(interval.first,interval.second);

             gen_t gen(urng,unif01_t());
@@ -153,17 +179,20 @@
             extr_state_t& extr_state = set.template extract<State>();
             extr_state.set(args,x);
             (this->init_).update(impl);
+ if(reserve_<impl.size_range_point()){
+ reserve_ = impl.size_range_point();
+ }
         }
         mutable init_t init_;
         typename super_t::size_type max_init_recursion_;
+ mutable typename super_t::size_type reserve_;
     };
}//impl

namespace feature{
- // Warning : pair_quantile_averaged_over_iterations is non_markov
- // Consequential theretically? Dunno, but from a practical
- // standpoint, less likely to lead to stability issues
- // such as exp(-inf).
+ // Warning : pair_quantile_averaged_over_iterations is likely
+ // to require fewer iterations than pair_quantile.
+ // Non-markovian : consequential for convergence?

     template<
         typename Par,

Modified: sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/support/detail/negative_max_to_max.hpp
==============================================================================
--- sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/support/detail/negative_max_to_max.hpp (original)
+++ sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/support/detail/negative_max_to_max.hpp 2009-03-26 18:24:23 EDT (Thu, 26 Mar 2009)
@@ -21,13 +21,11 @@
     support_negative_max_to_max(
         const D& d
     ){
- typedef typename function_support<D>::type fs;
+ typedef typename function_support<D>::type sup_t;
+ typedef typename function_argument<D>::type arg_t;

         using boost::math::tools::max_value;
- return fs(
- -max_value<typename fs::first_type>(),
- max_value<typename fs::second_type>()
- );
+ return sup_t(-max_value<arg_t>(),max_value<arg_t>());
     };

}//detail

Modified: sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/support/include.hpp
==============================================================================
--- sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/support/include.hpp (original)
+++ sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/support/include.hpp 2009-03-26 18:24:23 EDT (Thu, 26 Mar 2009)
@@ -8,7 +8,10 @@
#ifndef BOOST_CONDITIONALLY_SPECIFIED_DISTRIBUTION_SUPPORT_INCLUDE_HPP_ER_2009
#define BOOST_CONDITIONALLY_SPECIFIED_DISTRIBUTION_SUPPORT_INCLUDE_HPP_ER_2009

+#include <boost/conditionally_specified_distribution/support/pair_independent.hpp>
#include <boost/conditionally_specified_distribution/support/normal.hpp>
-#include <boost/conditionally_specified_distribution/support/fwd.hpp>
+
+// Bad idea:
+//#include <boost/conditionally_specified_distribution/support/fwd.hpp>

#endif

Modified: sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/support/normal.hpp
==============================================================================
--- sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/support/normal.hpp (original)
+++ sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/support/normal.hpp 2009-03-26 18:24:23 EDT (Thu, 26 Mar 2009)
@@ -7,6 +7,7 @@
////////////////////////////////////////////////////////////////////////////
#ifndef BOOST_CONDITIONALLY_SPECIFIED_DISTRIBUTION_SUPPORT_NORMAL_HPP_ER_2009
#define BOOST_CONDITIONALLY_SPECIFIED_DISTRIBUTION_SUPPORT_NORMAL_HPP_ER_2009
+#include <boost/utility/dont_care.hpp>
#include <boost/conditionally_specified_distribution/result_of/function_support.hpp>
#include <boost/conditionally_specified_distribution/crtp/normal.hpp>
#include <boost/conditionally_specified_distribution/support/detail/negative_max_to_max.hpp>
@@ -24,6 +25,14 @@
return detail::support_negative_max_to_max(d_);
}

+template <typename D>
+inline
+typename function_support<D>::type
+support(const crtp::normal<D>& d,utility::dont_care)
+{
+ return support(d);
+}
+
}
}
}

Added: sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/support/pair_independent.hpp
==============================================================================
--- (empty file)
+++ sandbox/conditionally_specified_distribution/conditionally_specified_distribution/boost/conditionally_specified_distribution/support/pair_independent.hpp 2009-03-26 18:24:23 EDT (Thu, 26 Mar 2009)
@@ -0,0 +1,34 @@
+/////////////////////////////////////////////////////////////////////////////
+// pair_independent.hpp
+//
+// Copyright 2009 Erwann Rogard. Distributed under the Boost
+// Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+////////////////////////////////////////////////////////////////////////////
+#ifndef BOOST_CONDITIONALLY_SPECIFIED_DISTRIBUTION_SUPPORT_PAIR_INDEPENDENT_HPP_ER_2009
+#define BOOST_CONDITIONALLY_SPECIFIED_DISTRIBUTION_SUPPORT_PAIR_INDEPENDENT_HPP_ER_2009
+#include <boost/conditionally_specified_distribution/result_of/function_support.hpp>
+#include <boost/conditionally_specified_distribution/crtp/pair_independent.hpp>
+namespace boost{
+namespace conditionally_specified_distribution{
+namespace function{
+
+template <typename D,typename Args>
+inline
+typename function_support<D>::type
+support(const crtp::pair_independent<D>& d,const Args& args)
+{
+ const D& d_ = static_cast<const D&>(d);
+
+ // TODO :
+ // pdf(x) = pdf_1(x)pdf_2(x) = 0 if either pdf_1(x) = 0 or pdf_2(y) =0
+ // so support (x,y) is intersection(suppor_1,support_2)
+
+ // For now, assume same support:
+ return support(d_.first(args));
+}
+
+}
+}
+}
+#endif // BOOST_CONDITIONALLY_SPECIFIED_DISTRIBUTION_SUPPORT_PAIR_INDEPENDENT_HPP_ER_2009

Modified: sandbox/conditionally_specified_distribution/conditionally_specified_distribution/libs/conditionally_specified_distribution/src/doc/readme.txt
==============================================================================
--- sandbox/conditionally_specified_distribution/conditionally_specified_distribution/libs/conditionally_specified_distribution/src/doc/readme.txt (original)
+++ sandbox/conditionally_specified_distribution/conditionally_specified_distribution/libs/conditionally_specified_distribution/src/doc/readme.txt 2009-03-26 18:24:23 EDT (Thu, 26 Mar 2009)
@@ -229,7 +229,6 @@
////////////////////////
/ Output from main.cpp /
////////////////////////
-
->example_parameter_normal_given_normal
mu_y|x=0.1
beta_y|x=0.75
@@ -276,12 +275,32 @@
1000000 samples from (x,y) jointly normal, mu_x = -0.1 sigma_x = 1
rho_xy = 0.9 mu_y = 0.1 sigma_y = 1.5 using the decomposition (y|x, x|y)
and an ars sampler for each component.
+i=0
+i=100000
+i=200000
+i=300000
+i=400000
+i=500000
+i=600000
+i=700000
+i=800000
+i=900000
sample statistics :
- mu_x : -0.0976351
- sigma_x : 0.972748
- mu_y : 0.103237
- sigma_y : 1.45947
- rho_xy :0.896274
+ mu_x : -0.0984778
+ sigma_x : 0.974992
+ mu_y : 0.102567
+ sigma_y : 1.46175
+ rho_xy :0.896784
<-
->example_function_pair_independent
OK <-
+
+
+///////////
+/ History /
+///////////
+
+On 2009/03/18, added:
+
+parameter::dependent::detail::regression_coefficient::as_string(args)
+

Modified: sandbox/conditionally_specified_distribution/conditionally_specified_distribution/libs/conditionally_specified_distribution/src/example/sampler/gibbs_ars_bivariate_normal.cpp
==============================================================================
--- sandbox/conditionally_specified_distribution/conditionally_specified_distribution/libs/conditionally_specified_distribution/src/example/sampler/gibbs_ars_bivariate_normal.cpp (original)
+++ sandbox/conditionally_specified_distribution/conditionally_specified_distribution/libs/conditionally_specified_distribution/src/example/sampler/gibbs_ars_bivariate_normal.cpp 2009-03-26 18:24:23 EDT (Thu, 26 Mar 2009)
@@ -8,8 +8,10 @@
#include <numeric>
#include <vector>
#include <cmath>
+#include <string>
#include <boost/assert.hpp>
#include <boost/bind.hpp>
+#include <boost/format.hpp>
#include <boost/random/mersenne_twister.hpp>
#include <boost/type_traits.hpp>
#include <boost/mpl/size_t.hpp>
@@ -24,6 +26,7 @@

#include <boost/shared_features/feature/keyword.hpp>
#include <boost/shared_features/set.hpp>
+#include <boost/shared_features/parameter.hpp>
#include <boost/shared_features/mpl_features.hpp>
#include <boost/shared_features/feature/scalar.hpp>
#include <boost/shared_features/functor/sample.hpp>
@@ -40,7 +43,7 @@
     std::cout << "->example_sampler_gibbs_ars_bivariate_normal" << std::endl;

     using namespace boost;
- namespace shared = shared_features;
+ namespace shared = boost::shared_features;
     namespace csd = conditionally_specified_distribution;
     namespace par = csd::parameter;
     namespace sampl = csd::sampler;
@@ -186,24 +189,51 @@
     std::cout << the_format.str() << std::endl;

- for(unsigned i=0; i<n; ++i){
- //if(i%1e4==0){
- // std::cout << "i=" << i << std::endl;
- //}
- set.visit_if<is_sampler_x>(
- shared::functor::sample(),
- (sampl::kwd<>::random_number_generator = urng)
+ try{
+ for(unsigned i=0; i<n; ++i){
+ unsigned mod = 1e5;
+ if(i%mod==0){
+ std::cout << "i=" << i << std::endl;
+ }
+ try{
+ set.visit_if<is_sampler_x>(
+ shared::functor::sample(),
+ (sampl::kwd<>::random_number_generator = urng)
+ );
+ }catch(std::exception& e){
+ std::string str = "x|y : ";
+ str+= e.what();
+ throw std::runtime_error(str);
+ }
+ try{
+ draws_x.push_back(
+ set.extract<state_x>()()
+ );
+ set.visit_if<is_sampler_y>(
+ shared::functor::sample(),
+ (sampl::kwd<>::random_number_generator = urng)
+ );
+ }catch(std::exception& e){
+ std::string str = "y|x : ";
+ str+= e.what();
+ throw std::runtime_error(str);
+ }
+ draws_y.push_back(
+ set.extract<state_y>()()
+ );
+ }
+ }catch(std::exception& e){
+ std::string str = e.what();
+ str += "y|x=";
+ str += set.extract<par_y_given_x>().as_string(
+ (shared::kwd_set = set)
         );
- draws_x.push_back(
- set.extract<state_x>()()
- );
- set.visit_if<is_sampler_y>(
- shared::functor::sample(),
- (sampl::kwd<>::random_number_generator = urng)
- );
- draws_y.push_back(
- set.extract<state_y>()()
+ str += "x|y=";
+ str += set.extract<par_x_given_y>().as_string(
+ (shared::kwd_set = set)
         );
+
+ std::cout << str << std::endl;
     }

     std::cout << "sample statistics : " << std::endl;

Modified: sandbox/conditionally_specified_distribution/conditionally_specified_distribution/libs/conditionally_specified_distribution/src/main.cpp
==============================================================================
--- sandbox/conditionally_specified_distribution/conditionally_specified_distribution/libs/conditionally_specified_distribution/src/main.cpp (original)
+++ sandbox/conditionally_specified_distribution/conditionally_specified_distribution/libs/conditionally_specified_distribution/src/main.cpp 2009-03-26 18:24:23 EDT (Thu, 26 Mar 2009)
@@ -19,7 +19,6 @@
{
     example_parameter_normal_given_normal();

-
     //normal (from simple to complex):
     example_sampler_exact_normal_given_normal();
     example_sampler_exact_bivariate_normal();
@@ -29,6 +28,5 @@

     example_function_pair_independent();

-
     return 0;
}

Modified: sandbox/conditionally_specified_distribution/math_limit/boost/math_limit/miscellanea.hpp
==============================================================================
--- sandbox/conditionally_specified_distribution/math_limit/boost/math_limit/miscellanea.hpp (original)
+++ sandbox/conditionally_specified_distribution/math_limit/boost/math_limit/miscellanea.hpp 2009-03-26 18:24:23 EDT (Thu, 26 Mar 2009)
@@ -6,56 +6,119 @@
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#ifndef BOOST_MATH_LIMITS_MISCELLANEA_HPP_ER_2009
#define BOOST_MATH_LIMITS_MISCELLANEA_HPP_ER_2009
+#include <iostream>
+#include <cmath>
+#include <stdexcept>
+#include <boost/assert.hpp>
+#include <boost/format.hpp>
+#include <boost/mpl/identity.hpp>
+#include <boost/static_assert.hpp>
+#include <boost/math/special_functions/next.hpp>
#include <boost/math/tools/precision.hpp>
+#include <boost/numeric/conversion/converter.hpp>
+#include <boost/math/special_functions/fpclassify.hpp>
namespace boost{
namespace math_limit{

+ //Discussed in
+ //comp.lang.c++/browse_thread/thread/976701405e6ad0cf?hl=en#
         template<typename RealType>
- struct limits{
- static RealType log_max_value(){
- return math::tools::log_max_value<RealType>();
+ struct exp_support{
+ static RealType highest(){
+ struct check{
+ // Warning: fails with RealType = long double with
+ // i686-apple-darwin9-gcc-4.0.1 (GCC) 4.0.1
+ check(RealType x){
+ BOOST_ASSERT(!math::isinf(exp(x)));
+ }
+ };
+ static RealType result
+ = math::tools::log_max_value<RealType>();
+ static check const the_check( result );
+ return result;
             }
- static RealType log_min_value(){
- return math::tools::log_min_value<RealType>();
- }
- static RealType negative_infinity(){
- return -infinity();
- }
- static RealType min_value(){
- //beware math::tools::min_value<RealType>() is the
- // min positive value (for floats);
- return math::tools::min_value<RealType>();
- }
- static RealType infinity(){
- return math::tools::max_value<RealType>();
- }
- static RealType is_negative_infinity(RealType x){
- return !(x>negative_infinity());
- }
- static RealType is_infinity(RealType x){
- return !(x<infinity());
+
+ static RealType lowest(){
+ struct check{
+ //TODO is this what needs to be checked?
+ check(RealType x){
+ BOOST_ASSERT(!math::isnan(exp(x)));
+ }
+ };
+ static RealType result
+ = math::tools::log_min_value<RealType>();
+ static check const the_check( result );
+ return result;
             }
         };

- template<typename RealType>
- RealType
- safer_exp(RealType a){
- typedef math_limit::limits<RealType> limits_t;
- RealType e = 0.0;
- RealType lmx = limits_t::log_max_value();
- RealType lmi = limits_t::log_min_value();
- if(a < lmx){
- if(a > lmi){
- e = exp(a); // -inf < a < inf
+
+ // Warning : if possible avoid truncated expression
+ // Better to dispatch on inf
+ template<typename RealType>
+ RealType
+ truncated_exp(RealType a){
+ typedef math_limit::exp_support<RealType> sup_t;
+ RealType e = 0.0;
+ RealType h = sup_t::highest();
+ RealType l = sup_t::lowest();
+ if(a < h){
+ if(a > l){
+ e = exp(a); // -inf < a < inf
+ }else{
+ e = 0.0;
+ }
             }else{
- e = 0.0;
+ e = math::tools::max_value<RealType>(); // a = inf
             }
- }else{
- e = limits_t::infinity(); // a = inf
+ BOOST_ASSERT(!std::isnan(e));
+ return e;
         }
- BOOST_ASSERT(!std::isnan(e));
- return e;
- }
+
+
+ struct scale_to_max{
+ //news://news.gmane.org:119/gq3755$vb9$1@ger.gmane.org
+ template<class T>// (x/c)+(y/c) <= max
+ static T plus(
+ const T& x,
+ const T& y
+ ) {
+ double max = (std::numeric_limits<T>::max)();
+ if(x > max - y) {
+ double c = x/max + y/max;
+ while(!boost::math::isinf(x/c + y/c)) {
+ c = boost::math::float_prior(c);
+ }
+ while(boost::math::isinf(x/c + y/c)) {
+ c = boost::math::float_next(c);
+ }
+ return(c);
+ } else {
+ return(static_cast<T>(1));
+ }
+ }
+ template<class T>// (x/c)+(y/c) <= max
+ static T divides(
+ const T& x,
+ const T& y
+ ) {
+ double max = (std::numeric_limits<T>::max)();
+ if(x / max > y) {
+ double c = (x/max)/y;
+ while(!boost::math::isinf((x/c)/y)) {
+ c = boost::math::float_prior(c);
+ }
+ while(boost::math::isinf((x/c)/y)) {
+ c = boost::math::float_next(c);
+ }
+ return(c);
+ } else {
+ return(static_cast<T>(1));
+ }
+ }
+
+ };
+

}
}

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