ide.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/ide.h
 *
 *
 */
 
/*!  \file ide.h
 *
 * \brief Definitions for inversion and decomposition functions that
 * operate on Scythe's Matrix objects.
 *
 * This file provides a number of common inversion and decomposition
 * routines that operate on Matrix objects.  It also provides related
 * functions for solving linear systems of equations and calculating
 * the determinant of a Matrix.
 *
 * Scythe will use LAPACK/BLAS routines to perform these operations on
 * concrete column-major matrices of double-precision floating point
 * numbers if LAPACK/BLAS is available and you compile your program
 * with the SCYTHE_LAPACK flag enabled.
 *
 * \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.
 */

/* TODO: This interface exposes the user to too much implementation.
 * We need a solve function and a solver object.  By default, solve
 * would run lu_solve and the solver factory would return lu_solvers
 * (or perhaps a solver object encapsulating an lu_solver).  Users
 * could choose cholesky when appropriate.  Down the road, qr or svd
 * would become the default and we'd be able to handle non-square
 * matrices.  Instead of doing an lu_decomp or a cholesky and keeping
 * track of the results to repeatedly solve for different b's with A
 * fixed in Ax=b, you'd just call the operator() on your solver object
 * over and over, passing the new b each time.  No decomposition
 * specific solvers (except as toggles to the solver object and
 * solve function).  We'd still provide cholesky and lu_decomp.  We
 * could also think about a similar approach to inversion (one
 * inversion function with an option for method).
 *
 * If virtual dispatch in C++ wasn't such a performance killer (no
 * compiler optimization across virtual calls!!!) there would be an
 * obvious implementation of this interface using simple polymorphism.
 * Unfortunately, we need compile-time typing to maintain performance
 * and makes developing a clean interface that doesn't force users to
 * be template wizards much harder.  Initial experiments with the
 * Barton and Nackman trick were ugly.  The engine approach might work
 * a bit better but has its problems too.  This is not going to get
 * done for the 1.0 release, but it is something we should come back
 * to.
 *
 */

#ifndef SCYTHE_IDE_H
#define SCYTHE_IDE_H

#ifdef SCYTHE_COMPILE_DIRECT
#include "matrix.h"
#include "error.h"
#include "defs.h"
#ifdef SCYTHE_LAPACK
#include "lapack.h"
#include "stat.h"
#endif
#else
#include "scythestat/matrix.h"
#include "scythestat/error.h"
#include "scythestat/defs.h"
#ifdef SCYTHE_LAPACK
#include "scythestat/lapack.h"
#include "scythestat/stat.h"
#endif
#endif

#include <cmath>
#include <algorithm>

namespace scythe {

  namespace {
    typedef unsigned int uint;
  }
  
  /*! 
   * \brief Cholesky decomposition of a symmetric positive-definite
   * matrix.
   *
   * This function performs Cholesky decomposition.  That is, given a
   * symmetric positive definite Matrix, \f$A\f$, cholesky() returns a
   * lower triangular Matrix \f$L\f$ such that \f$A = LL^T\f$.  This
   * function is faster than lu_decomp() and, therefore, preferable in
   * cases where one's Matrix is symmetric positive definite.
   *
   * \param A The symmetric positive definite Matrix to decompose.
   *
   * \see chol_solve(const Matrix<T,PO1,PS1> &, const Matrix<T,PO2,PS2> &)
   * \see chol_solve(const Matrix<T,PO1,PS1> &, const Matrix<T,PO2,PS2> &, const Matrix<T,PO3,PS3> &)
   * \see lu_decomp(Matrix<T,PO1,PS1>, Matrix<T,PO2,Concrete>&, Matrix<T,PO3,Concrete>&, Matrix<unsigned int, PO4, Concrete>&)
   *
   * \throw scythe_alloc_error (Level 1)
   * \throw scythe_dimension_error (Level 1)
   * \throw scythe_null_error (Level 1)
   * \throw scythe_type_error (Level 2)
   * \throw scythe_alloc_error (Level 1)
   *
   */
  template <matrix_order RO, matrix_style RS, typename T,
            matrix_order PO, matrix_style PS>
  Matrix<T, RO, RS>
  cholesky (const Matrix<T, PO, PS>& A)
  {
    SCYTHE_CHECK_10(! A.isSquare(), scythe_dimension_error,
        "Matrix not square");
    SCYTHE_CHECK_10(A.isNull(), scythe_null_error,
        "Matrix is NULL");
    // Rounding errors can make this problematic.  Leaving out for now
    //SCYTHE_CHECK_20(! A.isSymmetric(), scythe_type_error,
    //    "Matrix not symmetric");
    
    Matrix<T,RO,Concrete> temp (A.rows(), A.cols(), false);
    T h;
    
    if (PO == Row) { // row-major optimized
      for (uint i = 0; i < A.rows(); ++i) {
        for (uint j = i; j < A.cols(); ++j) {
          h = A(i,j);
          for (uint k = 0; k < i; ++k)
            h -= temp(i, k) * temp(j, k);
          if (i == j) {
            SCYTHE_CHECK_20(h <= (T) 0, scythe_type_error,
                "Matrix not positive definite");

            temp(i,i) = std::sqrt(h);
          } else {
            temp(j,i) = (((T) 1) / temp(i,i)) * h;
            temp(i,j) = (T) 0;
          }
        }
      }
    } else { // col-major optimized
      for (uint j = 0; j < A.cols(); ++j) {
        for (uint i = j; i < A.rows(); ++i) {
          h = A(i, j);
          for (uint k = 0; k < j; ++k)
            h -= temp(j, k) * temp(i, k);
          if (i == j) {
            SCYTHE_CHECK_20(h <= (T) 0, scythe_type_error,
                "Matrix not positive definite");
            temp(j,j) = std::sqrt(h);
          } else {
            temp(i,j) = (((T) 1) / temp(j,j)) * h;
            temp(j,i) = (T) 0;
          }
        }
      }
    }
  
    SCYTHE_VIEW_RETURN(T, RO, RS, temp)
  }

  template <typename T, matrix_order O, matrix_style S>
  Matrix<T, O, Concrete>
  cholesky (const Matrix<T,O,S>& A)
  {
    return cholesky<O,Concrete>(A);
  }

  namespace {
    /* This internal routine encapsulates the  
     * algorithm used within chol_solve and lu_solve.  
     */
    template <typename T,
              matrix_order PO1, matrix_style PS1,
              matrix_order PO2, matrix_style PS2,
              matrix_order PO3, matrix_style PS3>
    inline void
    solve(const Matrix<T,PO1,PS1>& L, const Matrix<T,PO2,PS2>& U,
          Matrix<T,PO3,PS3> b, T* x, T* y)
    {
      T sum;

      /* TODO: Consider optimizing for ordering.  Experimentation
       * shows performance gains are probably minor (compared col-major
       * with and without lapack solve routines).
       */
      // solve M*y = b
      for (uint i = 0; i < b.size(); ++i) {
        sum = T (0);
        for (uint j = 0; j < i; ++j) {
          sum += L(i,j) * y[j];
        }
        y[i] = (b[i] - sum) / L(i, i);
      }
          
      // solve M'*x = y
      if (U.isNull()) { // A= LL^T
        for (int i = b.size() - 1; i >= 0; --i) {
          sum = T(0);
          for (uint j = i + 1; j < b.size(); ++j) {
            sum += L(j,i) * x[j];
          }
          x[i] = (y[i] - sum) / L(i, i);
        }
      } else { // A = LU
        for (int i = b.size() - 1; i >= 0; --i) {
          sum = T(0);
          for (uint j = i + 1; j < b.size(); ++j) {
            sum += U(i,j) * x[j];
          }
          x[i] = (y[i] - sum) / U(i, i);
        }
      }
    }
  }

 /*!\brief Solve \f$Ax=b\f$ for x via backward substitution, given a
  * lower triangular matrix resulting from Cholesky decomposition
  *
  * This function solves the system of equations \f$Ax = b\f$ via
  * backward substitution. \a L is the lower triangular matrix generated
  * by Cholesky decomposition such that \f$A = LL'\f$.
  *
  * This function is intended for repeatedly solving systems of
  * equations based on \a A.  That is \a A stays constant while \a
  * b varies.
  *
  * \param A A symmetric positive definite Matrix.
  * \param b A column vector with as many rows as \a A.
  * \param M The lower triangular matrix from the Cholesky decomposition of \a A.
  *
  * \see chol_solve(const Matrix<T,PO1,PS1>&, const Matrix<T,PO2,PS2>&)
  * \see cholesky(const Matrix<T, PO, PS>&)
  * \see lu_solve (const Matrix<T,PO1,PS1>&, const Matrix<T,PO2,PS2>&, const Matrix<T,PO3,PS3>&, const Matrix<T,PO4,PS4>&, const Matrix<unsigned int, PO5, PS5>&)
  * \see lu_solve (Matrix<T,PO1,PS1>, const Matrix<T,PO2,PS2>&)
  *
  * \throw scythe_alloc_error (Level 1)
  * \throw scythe_null_error (Level 1)
  * \throw scythe_dimension_error (Level 1)
  * \throw scythe_conformation_error (Level 1)
  *
  */
  template <matrix_order RO, matrix_style RS, typename T,
            matrix_order PO1, matrix_style PS1,
            matrix_order PO2, matrix_style PS2,
            matrix_order PO3, matrix_style PS3>
  Matrix<T,RO,RS>
  chol_solve (const Matrix<T,PO1,PS1>& A, const Matrix<T,PO2,PS2>& b,
              const Matrix<T,PO3,PS3>& M)
  {
    SCYTHE_CHECK_10(A.isNull(), scythe_null_error,
        "A is NULL")
    SCYTHE_CHECK_10(! b.isColVector(), scythe_dimension_error,
        "b must be a column vector");
    SCYTHE_CHECK_10(A.rows() != b.rows(), scythe_conformation_error,
        "A and b do not conform");
    SCYTHE_CHECK_10(A.rows() != M.rows(), scythe_conformation_error,
        "A and M do not conform");
    SCYTHE_CHECK_10(! M.isSquare(), scythe_dimension_error,
        "M must be square");

    T *y = new T[A.rows()];
    T *x = new T[A.rows()];
    
    solve(M, Matrix<>(), b, x, y);

    Matrix<T,RO,RS> result(A.rows(), 1, x);
     
    delete[]x;
    delete[]y;
   
    return result;
  }

  template <typename T, matrix_order PO1, matrix_style PS1,
            matrix_order PO2, matrix_style PS2,
            matrix_order PO3, matrix_style PS3>
  Matrix<T,PO1,Concrete>
  chol_solve (const Matrix<T,PO1,PS1>& A, const Matrix<T,PO2,PS2>& b,
              const Matrix<T,PO3,PS3>& M)
  {
    return chol_solve<PO1,Concrete>(A,b,M);
  }

 /*!\brief Solve \f$Ax=b\f$ for x via backward substitution, 
  * using Cholesky decomposition
  *
  * This function solves the system of equations \f$Ax = b\f$ via
  * backward substitution and Cholesky decomposition. \a A must be a
  * symmetric positive definite matrix for this method to work. This
  * function calls cholesky() to perform the decomposition.
  *
  * \param A A symmetric positive definite matrix.
  * \param b A column vector with as many rows as \a A.
  *
  * \see chol_solve(const Matrix<T,PO1,PS1>&, const Matrix<T,PO2,PS2>&, const Matrix<T,PO3,PS3>&)
  * \see cholesky(const Matrix<T, PO, PS>&)
  * \see lu_solve (const Matrix<T,PO1,PS1>&, const Matrix<T,PO2,PS2>&, const Matrix<T,PO3,PS3>&, const Matrix<T,PO4,PS4>&, const Matrix<unsigned int, PO5, PS5>&)
  * \see lu_solve (Matrix<T,PO1,PS1>, const Matrix<T,PO2,PS2>&)
  *
  * \throw scythe_alloc_error (Level 1)
  * \throw scythe_null_error (Level 1)
  * \throw scythe_conformation_error (Level 1)
  * \throw scythe_dimension_error (Level 1)
  * \throw scythe_type_error (Level 2)
  * \throw scythe_alloc_error (Level 1)
  *
  */
  template <matrix_order RO, matrix_style RS, typename T,
            matrix_order PO1, matrix_style PS1,
            matrix_order PO2, matrix_style PS2>
  Matrix<T,RO,RS>
  chol_solve (const Matrix<T,PO1,PS1>& A, const Matrix<T,PO2,PS2>& b)
  {
    /* NOTE: cholesky() call does check for square/posdef of A,
     * and the overloaded chol_solve call handles dimensions
     */
  
    return chol_solve<RO,RS>(A, b, cholesky<RO,Concrete>(A));
  }
 
  template <typename T, matrix_order PO1, matrix_style PS1,
            matrix_order PO2, matrix_style PS2>
  Matrix<T,PO1,Concrete>
  chol_solve (const Matrix<T,PO1,PS1>& A, const Matrix<T,PO2,PS2>& b)
  {
    return chol_solve<PO1,Concrete>(A, b);
  }

  
  /*!\brief Calculates the inverse of a symmetric positive definite
   * matrix, given a lower triangular matrix resulting from Cholesky
   * decomposition.
   *
   * This function returns the inverse of a symmetric positive
   * definite matrix. Unlike the one-parameter version, this function
   * requires the caller to perform Cholesky decomposition on the