Boost C++ Libraries Home Libraries People FAQ More
PrevUpHomeNext

為特殊函數作圖，繪製輪廓以及生成測試數據(Graphing, Profiling, and Generating Test Data for Special Functions)

設計類test_data 以及相關的輔助函數使得你可以只使用幾行代碼就可以做到：

概要 
namespace boost{ namespace math{ namespace tools{

enum parameter_type
{
   random_in_range = 0,
   periodic_in_range = 1,
   power_series = 2,
   dummy_param = 0x80,
};

template <class T>
struct parameter_info;

template <class T>
parameter_info<T> make_random_param(T start_range, T end_range, int n_points);

template <class T>
parameter_info<T> make_periodic_param(T start_range, T end_range, int n_points);

template <class T>
parameter_info<T> make_power_param(T basis, int start_exponent, int end_exponent);

template <class T>
bool get_user_parameter_info(parameter_info<T>& info, const char* param_name);

template <class T>
class test_data
{
public:
   typedef std::vector<T> row_type;
   typedef row_type value_type;
private:
   typedef std::set<row_type> container_type;
public:
   typedef typename container_type::reference reference;
   typedef typename container_type::const_reference const_reference;
   typedef typename container_type::iterator iterator;
   typedef typename container_type::const_iterator const_iterator;
   typedef typename container_type::difference_type difference_type;
   typedef typename container_type::size_type size_type;

   // 構造:
   test_data(){}
   template <class F>
   test_data(F func, const parameter_info<T>& arg1);

   // 插入:
   template <class F>
   test_data& insert(F func, const parameter_info<T>& arg1);

   template <class F>
   test_data& insert(F func, const parameter_info<T>& arg1, 
                     const parameter_info<T>& arg2);

   template <class F>
   test_data& insert(F func, const parameter_info<T>& arg1, 
                     const parameter_info<T>& arg2, 
                     const parameter_info<T>& arg3);

   void clear();

   // 訪問:
   iterator begin();
   iterator end();
   const_iterator begin()const;
   const_iterator end()const;
   bool operator==(const test_data& d)const;
   bool operator!=(const test_data& d)const;
   void swap(test_data& other);
   size_type size()const;
   size_type max_size()const;
   bool empty()const;

   bool operator < (const test_data& dat)const;
   bool operator <= (const test_data& dat)const;
   bool operator > (const test_data& dat)const;
   bool operator >= (const test_data& dat)const;
};

template <class charT, class traits, class T>
std::basic_ostream<charT, traits>& write_csv(
            std::basic_ostream<charT, traits>& os,
            const test_data<T>& data);

template <class charT, class traits, class T>
std::basic_ostream<charT, traits>& write_csv(
            std::basic_ostream<charT, traits>& os,
            const test_data<T>& data,
            const charT* separator);

template <class T>
std::ostream& write_code(std::ostream& os,
                         const test_data<T>& data, 
                         const char* name);
                         
}}} // namespaces
說明

這些工具由下面的級數例子進行了很好的說明說明：

test_data 的功能分為以下幾個部分：

例一：輸出數據用於繪圖

假如，假設我們想要在 -3 到 100之間為函數lgamma繪圖，可以通過下面的代碼做到：

#include <boost/math/tools/test_data.hpp>
#include <boost/math/special_functions/gamma.hpp>

int main()
{
   using namespace boost::math::tools;
   
   // 生成一個對像來保存這些數據:
   test_data<double> data;
   
   // 在 -3 到 100之間均勻插入500個點:
   double (*pf)(double) = boost::math::lgamma;
   data.insert(pf, make_periodic_param(-3.0 + 0.00001, 100.0, 500));
   
   // 以 csv 格式打印:
   write_csv(std::cout, data, ", ");
   return 0;
}

當將結果進行繪圖時，圖像為：

例2：生成測試數據

作為第二個例子，假設我們想要為一個特殊函數生成高精度的測試數據。因為許多測試函數有兩個或更多的相互獨立的參數，在沒有耗費巨大的計算量的情況下生成千兆字節的測試數據很難有效地完全覆蓋所有可能的參數空間。一個好的方法是為用戶提供一個可以在一個指定區間上快速生成測試數據的工具。

在這個例子中，我們將使用精度為1000bit的NTL::RR 庫為β函數生成測試數據。我們使用lgamma函數有意地實現了一個相當簡單的β函數而沒有使用庫中實現的通用版本的β函數，並依賴於高精度的數據類型來獲取精度至少為128bit的結果。通過這種方式，我們的測試數據就會完全獨立於我們在β函數內部所使用的所有的技巧。

讓我們定義生成測試數據的函數對像：

#include <boost/math/tools/ntl.hpp>
#include <boost/math/special_functions/gamma.hpp>
#include <boost/math/tools/test_data.hpp>
#include <fstream>

#include <boost/math/tools/test_data.hpp>

using namespace boost::math::tools;

struct beta_data_generator
{
   NTL::RR operator()(NTL::RR a, NTL::RR b)
   {
      //
      // 如果我們拋出一個定義域錯誤的異常，那麼 test_data 
      // 將會忽略這個輸出點，我們將會使用這種方式來過濾
      //  a < b 的情況，因為β函數
      // 在a和b上是對象的:
      //
      if(a < b)
         throw std::domain_error("");
         
      // very naively calculate spots with lgamma:
      NTL::RR g1, g2, g3;
      int s1, s2, s3;
      g1 = boost::math::lgamma(a, &s1);
      g2 = boost::math::lgamma(b, &s2);
      g3 = boost::math::lgamma(a+b, &s3);
      g1 += g2 - g3;
      g1 = exp(g1);
      g1 *= s1 * s2 * s3;
      return g1;
   }
};

為了生成數據，我們需要輸出將進行測試的 a 和 b的定義域：函數get_user_parameter_info 函數是出於這個目的而設計的。main函數的開始看起來是這樣的：

// 設置 RR 庫的精度:
NTL::RR::SetPrecision(1000); // bits.
NTL::RR::SetOutputPrecision(40); // 十進制數字個數.

parameter_info<NTL::RR> arg1, arg2;
test_data<NTL::RR> data;

std::cout << "Welcome.\n"
   "This program will generate spot tests for the beta function:\n"
   "  beta(a, b)\n\n";

bool cont;
std::string line;

do{
   // 提示用戶 a 和 b 的定義域:
   get_user_parameter_info(arg1, "a");
   get_user_parameter_info(arg2, "b");
   
   // 生成測試數據:
   data.insert(beta_data_generator(), arg1, arg2);
   
   // 判斷用戶是否希望進行更多的測試:
   std::cout << "Any more data [y/n]?";
   std::getline(std::cin, line);
   boost::algorithm::trim(line);
   cont = (line == "y");
}while(cont);
[Caution] 注意

在這裡，一個潛在的障礙需要提及：當有兩個或更多參數的時候，函數test_data<>::insert 交付生成測試數據的一個矩陣，所以，如果我們有兩個參數，並且對於每個參數我們需要1000個測試數據，那麼會有總共會有一百萬個測試點 。不要說沒有提醒過你!

這裡是程序的最後一步，將測試數據寫入到文件中：

std::cout << "Enter name of test data file [default=beta_data.ipp]";
std::getline(std::cin, line);
boost::algorithm::trim(line);
if(line == "")
   line = "beta_data.ipp";
std::ofstream ofs(line.c_str());
write_code(ofs, data, "beta_data");

測試數據的格式類似於下面：

#define SC_(x) static_cast<T>(BOOST_JOIN(x, L))
   static const boost::array<boost::array<T, 3>, 1830>
   beta_med_data = {
      SC_(0.4883005917072296142578125),
      SC_(0.4883005917072296142578125),
      SC_(3.245912809500479157065104747353807392371), 
      SC_(3.5808107852935791015625),
      SC_(0.4883005917072296142578125),
      SC_(1.007653173802923954909901438393379243537), 
      /* ... lots of rows skipped */
};

第五行開始的兩個數據是我們作遞給函數對象的參數，而最後一個數據是從函數對像中返回的結果。如果我們的函數對像返回一個元組(tuple)，那麼除了輸入參數之外 ，對於元組中的每個元素我們還可以有一個入口 (entry)。

開始的 #define 滿足兩個目的：

#define SC_(x) lexical_cast<T>(BOOST_STRINGIZE(x))

來確保在轉換為類型T之前沒有數值截斷。注意，這並不是缺省的設置，因為當表太大時，對於編譯器將是一項艱難的工作。

例 3:為連分數的收斂和精度繪圖

另一方面，假設我們想要為一個連分數的收斂和誤差繪圖。作為一個例子，我們將使用上部不完全γ函數(upper incomplete gamma function)的連分數，在每次調用這個函數對像時，它返回連分數的下一個 aN 和 bN ： 

template <class T>
struct upper_incomplete_gamma_fract
{
private:
   T z, a;
   int k;
public:
   typedef std::pair<T,T> result_type;

   upper_incomplete_gamma_fract(T a1, T z1)
      : z(z1-a1+1), a(a1), k(0)
   {
   }

   result_type operator()()
   {
      ++k;
      z += 2;
      return result_type(k * (a - k), z);
   }
};

我們想要測試連分數的相對誤差和收斂速率，因此我們寫出一個將兩個值做為一個元組返回的函數： test_data 類將為我們取出這個元組，並為元組中的每個數據生成一列數據 (除了輸出參數之外)：

#include <boost/math/tools/test_data.hpp>
#include <boost/math/tools/test.hpp>
#include <boost/math/special_functions/gamma.hpp>
#include <boost/math/tools/ntl.hpp>
#include <boost/tr1/tuple.hpp>

template <class T>
struct profile_gamma_fraction
{
   typedef std::tr1::tuple<T, T> result_type;

   result_type operator()(T val)
   {
      using namespace boost::math::tools;
      // 估測真值，使用任意精度的算術
      // 以及 NTL::RR庫:
      NTL::RR rval(val);
      upper_incomplete_gamma_fract<NTL::RR> f1(rval, rval);
      NTL::RR true_val = continued_fraction_a(f1, 1000);
      //
      // 在double精確度獲得近似值以及要求的迭代次數:
      boost::uintmax_t iters = std::numeric_limits<boost::uintmax_t>::max();
      upper_incomplete_gamma_fract<T> f2(val, val);
      T found_val = continued_fraction_a(f2, std::numeric_limits<T>::digits, iters);
      //
      // 算出相對誤差， 以10的-5次方為單位:
      T err = real_cast<T>(relative_error(true_val, NTL::RR(found_val)) / std::numeric_limits<T>::epsilon());
      //
      // 將結果作為一個元組(tuple)返回:
      return std::tr1::make_tuple(err, iters);
   }
};

將這個函數對像傳入到類 test_data 中可以快速生成 csv 數據，對於任何我們感興趣的類型T：

int main()
{
   using namespace boost::math::tools;
   // 生成一個對像來保存這些數據:
   test_data<double> data;
   // 在 0之後 到 100之間均勻插入500個點:
   data.insert(profile_gamma_fraction<double>(), make_periodic_param(0.01, 100.0, 100));
   // 以 csv 格式打印:
   write_csv(std::cout, data, ", ");
   return 0;
}

這一次沒有必要生成圖像，開始的一些輸出結果為：

a and z,  Error/epsilon,  Iterations required

0.01,     9723.14,        4726
1.0099,   9.54818,        87
2.0098,   3.84777,        40
3.0097,   0.728358,       25
4.0096,   2.39712,        21
5.0095,   0.233263,       16

所以，非常明顯，對於小的a和z不應當使用這個分數。

參考

這些工具中的大部分都在上面的例子中進行了介紹，我們將只對非成員函數增加下面的說明：

template <class T> parameter_info<T> make_random_param(T start_range, T end_range, int n_points);

指示類 test_data 在區間 [start_range,end_range]中只測試 n_points 個隨機值。

template <class T>
parameter_info<T> make_periodic_param(T start_range, T end_range, int n_points);

指示類 test_data 在區間 [start_range,end_range]中只測試 n_points 個均勻分佈的值。

template <class T>
parameter_info<T> make_power_param(T basis, int start_exponent, int end_exponent);

指示類 test_data 測試形式為 + R * 2expon 的數據，其中每個expon 都在區間 [start_exponent, end_exponent]中，而每個R 是一個在 [0.5, 1]中的隨機值。

template <class T>
bool get_user_parameter_info(parameter_info<T>& info, const char* param_name);

提示用戶參數的範圍以及使用形式。

最後，如果我們不想參數包含在輸出結果中，我們可以通過設置"dummy parameter"來指示 test_data 類：

parameter_info<double> p = make_random_param(2.0, 5.0, 10);
p.type |= dummy_param;

在函數對像在將這個參數傳遞給正在測試的函數之前通過某種形式將參數進行了變化時，這就很有用，通常，函數對像會以一個元組的形式(tuple)將轉換參數和轉換結果輸出，所以，這就沒有必要使初始的偽參數(pseudo-parameter)出現在程序的輸出結果中。

Copyright ?2006 , 2007, 2008 John Maddock, Paul A. Bristow, Hubert Holin, Xiaogang Zhang and Bruno Lalande

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)

PrevUpHomeNext