smath.h

使用R语言的马尔科夫链蒙特卡洛模拟（MCMC）源代码程序。

/* 
 * Scythe Statistical Library Copyright (C) 2000-2002 Andrew D. Martin
 * and Kevin M. Quinn; 2002-present Andrew D. Martin, Kevin M. Quinn,
 * and Daniel Pemstein.  All Rights Reserved.
 *
 * This program is free software; you can redistribute it and/or
 * modify under the terms of the GNU General Public License as
 * published by Free Software Foundation; either version 2 of the
 * License, or (at your option) any later version.  See the text files
 * COPYING and LICENSE, distributed with this source code, for further
 * information.
 * --------------------------------------------------------------------
 *  scythestat/smath.h
 *
 */

/*!
 * \file smath.h
 * \brief Definitions for functions that perform common mathematical
 * operations on every element of a Matrix.
 * 
 * \note As is the case throughout the library, we provide both
 * general and default template definitions of the Matrix-returning
 * functions in this file, explicitly providing documentation for only
 * the general template versions. As is also often the case, Doxygen
 * does not always correctly add the default template definition to
 * the function list below; there is always a default template
 * definition available for every function.
 *
 */

#ifndef SCYTHE_MATH_H
#define SCYTHE_MATH_H

#ifdef SCYTHE_COMPILE_DIRECT
#include "matrix.h"
#include "algorithm.h"
#include "error.h"
#else
#include "scythestat/matrix.h"
#include "scythestat/algorithm.h"
#include "scythestat/error.h"
#endif

#include <cmath>
#include <numeric>
#include <set>

namespace scythe {

  namespace {
    typedef unsigned int uint;
  }

/* Almost every function in this file follows one of the two patterns
 * described by these macros.  The first macro handles single-argument
 * functions.  The second handles two-matrix-argument functions (or
 * scalar-matrix, matrix-scalar.  The second macro also permits
 * cross-type operations (these are limited only by the capabilities
 * of the underlying functions).
 */
#define SCYTHE_MATH_OP(NAME, OP)                                      \
  template <matrix_order RO, matrix_style RS, typename T,             \
            matrix_order PO, matrix_style PS>                         \
  Matrix<T,RO,RS>                                                     \
  NAME (const Matrix<T,PO,PS>& A)                                     \
  {                                                                   \
    Matrix<T,RO,RS> res(A.rows(), A.cols(), false);                   \
    std::transform(A.begin_f(), A.end_f(), res.begin_f(), OP);        \
    return res;                                                       \
  }                                                                   \
                                                                      \
  template <typename T, matrix_order O, matrix_style S>               \
  Matrix<T,O,Concrete>                                                \
  NAME (const Matrix<T,O,S>& A)                                       \
  {                                                                   \
    return NAME<O,Concrete>(A);                                       \
  }

#define SCYTHE_MATH_OP_2ARG(NAME, OP)                                 \
  template <matrix_order RO, matrix_style RS, typename T,             \
            matrix_order PO1, matrix_style PS1,                       \
            matrix_order PO2, matrix_style PS2, typename S>           \
  Matrix<T,RO,RS>                                                     \
  NAME (const Matrix<T,PO1,PS1>& A, const Matrix<S,PO2,PS2>& B)       \
  {                                                                   \
    SCYTHE_CHECK_10 (A.size() != 1 && B.size() != 1 &&                \
        A.size() != B.size(), scythe_conformation_error,              \
        "Matrices with dimensions (" << A.rows()                      \
        << ", " << A.cols()                                           \
        << ") and (" << B.rows() << ", " << B.cols()                  \
        << ") are not conformable");                                  \
                                                                      \
    Matrix<T,RO,RS> res;                                              \
                                                                      \
    if (A.size() == 1) {                                              \
      res.resize2Match(B);                                            \
      std::transform(B.template begin_f<RO>(), B.template end_f<RO>(),\
          res.begin_f(), std::bind1st(OP, A(0)));                     \
    } else if (B.size() == 1) {                                       \
      res.resize2Match(A);                                            \
      std::transform(A.template begin_f<RO>(), A.template end_f<RO>(),\
                     res.begin_f(), std::bind2nd(OP, B(0)));          \
    } else {                                                          \
      res.resize2Match(A);                                            \
      std::transform(A.template begin_f<RO>(), A.template end_f<RO>(),\
                     B.template begin_f<RO>(), res.begin_f(), OP);    \
    }                                                                 \
                                                                      \
    return res;                                                       \
  }                                                                   \
                                                                      \
  template <typename T, matrix_order PO1, matrix_style PS1,           \
                        matrix_order PO2, matrix_style PS2,           \
                        typename S>                                   \
  Matrix<T,PO1,Concrete>                                              \
  NAME (const Matrix<T,PO1,PS1>& A, const Matrix<S,PO2,PS2>& B)       \
  {                                                                   \
    return NAME<PO1,Concrete>(A, B);                                  \
  }                                                                   \
                                                                      \
  template<matrix_order RO, matrix_style RS, typename T,              \
           matrix_order PO, matrix_style PS, typename S>              \
  Matrix<T,RO,RS>                                                     \
  NAME (const Matrix<T,PO,PS>& A, S b)                                \
  {                                                                   \
    return NAME<RO,RS>(A, Matrix<S,RO,Concrete>(b));                  \
  }                                                                   \
                                                                      \
  template <typename T, typename S, matrix_order PO, matrix_style PS> \
  Matrix<T,PO,Concrete>                                               \
  NAME (const Matrix<T,PO,PS>& A, S b)                                \
  {                                                                   \
    return NAME<PO,Concrete>(A, Matrix<S,PO,Concrete>(b));            \
  }                                                                   \
                                                                      \
  template<matrix_order RO, matrix_style RS, typename T,              \
           matrix_order PO, matrix_style PS, typename S>              \
  Matrix<T,RO,RS>                                                     \
  NAME (T a, const Matrix<S,PO,PS>& B)                                \
  {                                                                   \
    return NAME<RO,RS>(Matrix<S, RO,Concrete>(a), B);                 \
  }                                                                   \
                                                                      \
  template <typename T, typename S, matrix_order PO, matrix_style PS> \
  Matrix<T,PO,Concrete>                                               \
  NAME (T a, const Matrix<S,PO,PS>& B)                                \
  {                                                                   \
    return NAME<PO,Concrete>(Matrix<S,PO,Concrete>(a), B);            \
  }


  /* calc the inverse cosine of each element of a Matrix */
  
 /*! 
	* \brief Calculate the inverse cosine of each element of a Matrix
	*
	* This function calculates the inverse cosine of each element in a Matrix
	*
	* \param A The matrix whose inverse cosines are of interest.
	*
	* \see tan()
	* \see tanh()
	* \see sin()
	* \see sinh()
	* \see cos()
	* \see cosh()
	* \see acosh()
	* \see asin()
	* \see asinh()
	* \see atan()
	* \see atanh()
	* \see atan2()
	*/

  SCYTHE_MATH_OP(acos, ::acos)
  
  /* calc the inverse hyperbolic cosine of each element of a Matrix */
   /*! 
	* \brief Calculate the inverse hyperbolic cosine of each element of a Matrix
	*
	* This function calculates the inverse hyperbolic cosine of each element
	* in a Matrix
	*
	* \param A The matrix whose inverse hyperbolic cosines are of interest.
	*
	* \see tan()
	* \see tanh()
	* \see sin()
	* \see sinh()
	* \see cos()
	* \see cosh()
	* \see acos()
	* \see asin()
	* \see asinh()
	* \see atan()
	* \see atanh()
	* \see atan2()
	*/

  SCYTHE_MATH_OP(acosh, ::acosh)

  /* calc the inverse sine of each element of a Matrix */
  
   /*! 
	* \brief Calculate the inverse sine of each element of a Matrix
	*
	* This function calculates the inverse sine of each element
	* in a Matrix
	*
	* \param A The matrix whose inverse sines are of interest.
	*
	* \see tan()
	* \see tanh()
	* \see sin()
	* \see sinh()
	* \see cos()
	* \see cosh()
	* \see acos()
	* \see acosh()
	* \see asinh()
	* \see atan()
	* \see atanh()
	* \see atan2()
	*/

  SCYTHE_MATH_OP(asin, ::asin)
  
  /* calc the inverse hyperbolic sine of each element of a Matrix */
  
  /*! 
	* \brief Calculate the inverse hyperbolic sine of each element of a Matrix
	*
	* This function calculates the inverse hyperbolic sine of each element
	* in a Matrix
	*
	* \param A The matrix whose inverse hyperbolic sines are of interest.
	*
	* \see tan()
	* \see tanh()
	* \see sin()
	* \see sinh()
	* \see cos()
	* \see cosh()
	* \see acos()
	* \see acosh()
	* \see asin()
	* \see atan()
	* \see atanh()
	* \see atan2()
	*/
	
  SCYTHE_MATH_OP(asinh, ::asinh)
  
  /* calc the inverse tangent of each element of a Matrix */
  
   /*! 
	* \brief Calculate the inverse tangent of each element of a Matrix
	*
	* This function calculates the inverse tangent of each element
	* in a Matrix
	*
	* \param A The matrix whose inverse tangents are of interest.
	*
	* \see tan()
	* \see tanh()
	* \see sin()
	* \see sinh()
	* \see cos()
	* \see cosh()
	* \see acos()
	* \see acosh()
	* \see asin()
	* \see asin()
	* \see atanh()
	* \see atan2()
	*/
	
  SCYTHE_MATH_OP(atan, ::atan)
  
  /* calc the inverse hyperbolic tangent of each element of a Matrix */
   /*! 
	* \brief Calculate the inverse hyperbolic tangent of each element of a Matrix
	*
	* This function calculates the inverse hyperbolic tangent of each element
	* in a Matrix
	*
	* \param A The matrix whose inverse hyperbolic tangents are of interest.
	*
	* \see tan()
	* \see tanh()
	* \see sin()
	* \see sinh()
	* \see cos()
	* \see cosh()
	* \see acos()
	* \see acosh()
	* \see asin()
	* \see asinh()
	* \see atan()
	* \see atan2()
	*/
	
  SCYTHE_MATH_OP(atanh, ::atanh)
  
  /* calc the angle whose tangent is y/x  */
  
   /*! 
	* \brief Calculate the angle whose tangent is y/x
	*
	* This function calculates the angle whose tangent is y/x, given two 
	* matrices A and B (where y is the ith element of A, and x is the jth element
	* of matrix B).
	*
	* \param A The matrix of y values 
	* \param B The matrix of x values
	*
	* \see tan()
	* \see tanh()
	* \see sin()
	* \see sinh()
	* \see cos()
	* \see cosh()
	* \see acos()
	* \see acosh()
	* \see asin()
	* \see asinh()
	* \see atan()
	* \see atanh()
	*/

  SCYTHE_MATH_OP_2ARG(atan2, std::ptr_fun(::atan2))

  /* calc the cube root of each element of a Matrix */
   /*! 
	* \brief Calculate the cube root of each element of a Matrix
	*
	* This function calculates the cube root of each element
	* in a Matrix
	*
	* \param A The matrix whose cube roots are of interest.
	*
	* \see sqrt()
	*/

  SCYTHE_MATH_OP(cbrt, ::cbrt)
  
  /* calc the ceil of each element of a Matrix */
  /*! 
	* \brief Calculate the ceiling of each element of a Matrix
	*
	* This function calculates the ceiling of each element
	* in a Matrix
	*
	* \param A The matrix whose ceilings are of interest.
	*
	* \see floor()
	*/

  SCYTHE_MATH_OP(ceil, ::ceil)
  
  /* create a matrix containing the absval of the first input and the
   * sign of the second
   */
    /*! 
	* \brief Create a matrix containing the absolute value of the first input
	* and the sign of the second input
	*
	* This function creates a matrix containing the absolute value of the first
	* input, a matrix called A, and the sign of the second input, matrix B.
	*
	* \param A The matrix whose absolute values will comprise the resultant matrix.
	* \param B The matrix whose signs will comprise the resultant matrix
	*/

  SCYTHE_MATH_OP_2ARG(copysign, std::ptr_fun(::copysign))
  
  /* calc the cosine of each element of a Matrix */
    
 /*! 
	* \brief Calculate the cosine of each element of a Matrix
	*
	* This function calculates the cosine of each element in a Matrix
	*
	* \param A The matrix whose cosines are of interest.
	*
	* \see tan()
	* \see tanh()
	* \see sin()
	* \see sinh()
	* \see cosh()
	* \see acos()
	* \see acosh()
	* \see asin()
	* \see asinh()
	* \see atan()
	* \see atanh()
	* \see atan2()
	*/

  SCYTHE_MATH_OP(cos, ::cos)
  
  /* calc the hyperbolic cosine of each element of a Matrix */
   /*! 
	* \brief Calculate the hyperbolic cosine of each element of a Matrix
	*
	* This function calculates the hyperbolic cosine of each element in a Matrix
	*
	* \param A The matrix whose hyperbolic cosines are of interest.
	*
	* \see tan()
	* \see tanh()
	* \see sin()
	* \see sinh()
	* \see cos()
	* \see acos()
	* \see acosh()
	* \see asin()
	* \see asinh()
	* \see atan()
	* \see atanh()
	* \see atan2()
	*/

  SCYTHE_MATH_OP(cosh, ::cosh)
  
  /* calc the error function of each element of a Matrix */
   /*! 
	* \brief Calculate the error function of each element of a Matrix
	*
	* This function calculates the error function of each element in a Matrix
	*
	* \param A The matrix whose error functions are of interest.
	*
	* \see erfc()
	*/

  SCYTHE_MATH_OP(erf, ::erf)
  
  /* calc the complementary error function of each element of a Matrix */
   /*! 
	* \brief Calculate the complementary error function of each element of a Matrix
	*
	* This function calculates the complemenatry error function of each 
	* element in a Matrix
	*
	* \param A The matrix whose complementary error functions are of interest.
	*
	* \see erf()
	*/

  SCYTHE_MATH_OP(erfc, ::erfc)
  
  /* calc the vaue e^x of each element of a Matrix */
   /*! 
	* \brief Calculate the value e^x for each element of a Matrix
	*
	* This function calculates the value e^x for each element of a matrix, where
	* x is the ith element of the matrix A
	*
	* \param A The matrix whose elements are to be exponentiated.
	*
	* \see expm1()
	*/

  SCYTHE_MATH_OP(exp, ::exp)
  
  /* calc the exponent - 1 of each element of a Matrix */
  /*! 
	* \brief Calculate the value e^(x-1) for each element of a Matrix
	*
	* This function calculates the value e^(x-1) for each element of a matrix, where
	* x is the ith element of the matrix A
	*
	* \param A The matrix whose elements are to be exponentiated.
	*
	* \see exp()
	*/