la.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/la.h
 *
 */

/*!
 * \file la.h
 * \brief Definitions and implementations for functions that perform
 * common linear algebra manipulations on Scythe Matrix objects.
 *
 * This file provides a number of common linear algebraic functions
 * for use with the Matrix class.  These functions include common
 * operations such as transposition, a number of utility functions for
 * creating useful matrices like the identity matrix, and efficient
 * implementations for common operations like the cross-product.
 *
 * \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.
 */

#ifndef SCYTHE_LA_H
#define SCYTHE_LA_H

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

#include <numeric>
#include <algorithm>
#include <set>

namespace scythe {

  namespace {
    typedef unsigned int uint;
  }

  /* Matrix transposition */

  /*!\brief Transpose a Matrix.
   *
   * This function transposes \a M, returning a Matrix \a R where each
   * element of \a M, \f$M_ij\f$ is placed in position \f$R_ji\f$.
   * Naturally, the returned Matrix has M.cols() rows and M.rows()
   * columns.
   *
   * \param M The Matrix to transpose.
   *
   * \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>
  t (const Matrix<T,PO,PS>& M)
  {
    uint rows = M.rows();
    uint cols = M.cols();
    Matrix<T,RO,Concrete> ret(cols, rows, false);
    if (PO == Col)
      copy<Col,Row>(M, ret);
    else
      copy<Row,Col>(M, ret);

    SCYTHE_VIEW_RETURN(T, RO, RS, ret)
  }

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

  /* Ones matrix generation */
  
  /*! 
	 * \brief Create a matrix of ones.
	 *
	 * This function creates a matrix of ones, with the given dimensions
	 * \a rows and \a cols.
	 *
	 * \param rows The number of rows in the resulting Matrix.
	 * \param cols The number of columns in the resulting Matrix.
	 *
	 * \see eye (unsigned int k)
   *
   * \throw scythe_alloc_error (Level 1)
	 */
  template <typename T, matrix_order O, matrix_style S>
  Matrix<T,O,S>
  ones (unsigned int rows, unsigned int cols)
  {
    return Matrix<T,O,S> (rows, cols, true, (T) 1);
  }

  template <typename T, matrix_order O>
  Matrix<T, O, Concrete>
  ones (unsigned int rows, unsigned int cols)
  {
    return ones<T, O, Concrete>(rows, cols);
  }

  template <typename T>
  Matrix<T, Col, Concrete>
  ones (unsigned int rows, unsigned int cols)
  {
    return ones<T, Col, Concrete>(rows, cols);
  }

  inline Matrix<double, Col, Concrete>
  ones (unsigned int rows, unsigned int cols)
  {
    return ones<double, Col, Concrete>(rows, cols);
  }

  /* Identity  Matrix generation */

  // This functor contains the working parts of the eye algorithm.
  namespace {
    template <class T> struct eye_alg {
        T operator() (uint i, uint j) {
          if (i == j)
            return (T) 1.0;
          return (T) 0.0;
        }
    };
  }
  
  /*!\brief Create a \a k by \a k identity Matrix.
   *
   * This function creates a \a k by \a k Matrix with 1s along the
   * diagonal and 0s on the off-diagonal.  This template is overloaded
   * multiple times to provide default type, matrix_order, and
   * matrix_style.  The default call to eye returns a Concrete Matrix
   * containing double precision floating point numbers, in
   * column-major order.  The user can write explicit template calls
   * to generate matrices with other orders and/or styles.
   *
   * \param k The dimension of the identity Matrix.
   *
   * \see diag(const Matrix<T,O,S>& M)
   * \see ones(unsigned int rows, unsigned int cols)
   *
   * \throw scythe_alloc_error (Level 1)
   * 
   */
  template <typename T, matrix_order O, matrix_style S>
  Matrix<T,O,S>
  eye (unsigned int k)
  {
    Matrix<T,O,Concrete> ret(k, k, false);
    for_each_ij_set(ret, eye_alg<T>());
    SCYTHE_VIEW_RETURN(T, O, S, ret)
  }

  template <typename T, matrix_order O>
  Matrix<T, O, Concrete>
  eye (uint k)
  {
    return eye<T, O, Concrete>(k);
  }

  template <typename T>
  Matrix<T, Col, Concrete>
  eye (uint k)
  {
    return eye<T, Col, Concrete>(k);
  }

  inline Matrix<double, Col, Concrete>
  eye (uint k)
  {
    return eye<double, Col, Concrete>(k);
  }

  /* Create a k x 1 vector-additive sequence matrix */

  // The seqa algorithm
  namespace {
    template <typename T> struct seqa_alg {
      T cur_; T inc_;
      seqa_alg(T start, T inc) : cur_ (start), inc_ (inc)  {}
      T operator() () { T ret = cur_; cur_ += inc_; return ret; }
    };
  }
  
 /*! 
	* \brief Create a \a rows x 1 vector-additive sequence Matrix.
	*
	* This function creates a \a rows x 1 Matrix \f$v\f$, where 
  * \f$v_i = \mbox{start} + i \cdot \mbox{incr}\f$.
  *
  * This function is defined by a series of templates.  This template
  * is the most general, requiring the user to explicitly instantiate
  * the template in terms of element type, matrix_order and
  * matrix_style.  Further versions allow for explicit instantiation
  * based just on type and matrix_order (with matrix_style defaulting
  * to Concrete) and just on type (with matrix_style defaulting to
  * Col).  Finally, the default version of th function generates
  * column-major concrete Matrix of doubles.
	*
	* \param start Desired start value.
	* \param incr Amount to add in each step of the sequence.
	* \param rows Total number of rows in the Matrix.
  *
   * \throw scythe_alloc_error (Level 1)
	*/

  template <typename T, matrix_order O, matrix_style S>
  Matrix<T,O,S>
  seqa (T start, T incr, uint rows)
  {
    Matrix<T,O,Concrete> ret(rows, 1, false);
    generate(ret.begin_f(), ret.end_f(), seqa_alg<T>(start, incr));
    SCYTHE_VIEW_RETURN(T, O, S, ret)
  }

  template <typename T, matrix_order O>
  Matrix<T, O, Concrete>
  seqa (T start, T incr, uint rows)
  {
    return seqa<T, O, Concrete>(start, incr, rows);
  }

  template <typename T>
  Matrix<T, Col, Concrete>
  seqa (T start, T incr, uint rows)
  {
    return seqa<T, Col, Concrete>(start, incr, rows);
  }

  inline Matrix<double, Col, Concrete>
  seqa (double start, double incr, uint rows)
  {
    return seqa<double, Col, Concrete>(start, incr, rows);
  }

  /* Uses the STL sort to sort a Matrix in ascending row-major order */
  
  /*! 
	 * \brief Sort a Matrix.
	 *
   * This function returns a copy of \a M, sorted in ascending order.
   * The sorting order is determined by the template parameter
   * SORT_ORDER or, by default, to matrix_order of \a M.
   *
	 * \param M The Matrix to sort.
   *
   * \see sortc
   * 
   * \throw scythe_alloc_error (Level 1)
	 */
  template <matrix_order SORT_ORDER,
            matrix_order RO, matrix_style RS, typename T,
            matrix_order PO, matrix_style PS>
  Matrix<T, RO, RS>
  sort (const Matrix<T, PO, PS>& M)
  {
    Matrix<T,RO,Concrete> ret = M;

    std::sort(ret.template begin<SORT_ORDER>(), 
              ret.template end<SORT_ORDER>());

    SCYTHE_VIEW_RETURN(T, RO, RS, ret)
  }

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

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

  /*!\brief Sort the columns of a Matrix.
   *
   * This function returns a copy of \a M, with each column sorted in
   * ascending order.
	 *
	 * \param M The Matrix to sort.
   * 
   * \see sort
   *
   * \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>
  sortc (const Matrix<T, PO, PS>& M)
  {
    Matrix<T,RO,Concrete> ret = M;

    // TODO need to figure out a way to do fully optimized
    // vector iteration
    for (uint col = 0; col < ret.cols(); ++col) {
      Matrix<T,PO,View> column = ret(_, col);
      std::sort(column.begin(), column.end());
    }

    SCYTHE_VIEW_RETURN(T, RO, RS, ret)
  }

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

  /* Column bind two matrices */
  
  /*! 
	 * \brief Column bind two matrices.
	 *
	 * This function column binds two matrices, \a A and \a B.
	 *
	 * \param A The left-hand Matrix.
	 * \param B The right-hand Matrix.
	 *
	 * \see rbind(const Matrix<T,PO1,PS1>& A, 
   *            const Matrix<T,PO2,PS2>& B) 
   *
   * \throw scythe_conformation_error (Level 1)
   * \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>
  cbind (const Matrix<T,PO1,PS1>& A, const Matrix<T,PO2,PS2>& B)
  {
    SCYTHE_CHECK_10(A.rows() != B.rows(), scythe_conformation_error,
        "Matrices have different numbers of rows");

    Matrix<T,RO,Concrete> ret(A.rows(), A.cols() + B.cols(), false);
    std::copy(B.template begin_f<Col>(), B.template end_f<Col>(),
              std::copy(A.template begin_f<Col>(), 
                        A.template end_f<Col>(), 
                        ret.template begin_f<Col>()));
    SCYTHE_VIEW_RETURN(T, RO, RS, ret)
  }


  template <typename T, matrix_order PO1, matrix_style PS1,
            matrix_order PO2, matrix_style PS2>
  Matrix<T,PO1,Concrete>
  cbind (const Matrix<T,PO1,PS1>& A, const Matrix<T,PO2,PS2>& B)
  {
    return cbind<PO1,Concrete>(A, B);
  }

  /* Row bind two matrices */
  
  /*! 
	 * \brief Row bind two matrices.
	 *
	 * This function row binds two matrices, \a A and \a B.
	 *
	 * \param A The upper Matrix.
	 * \param B The lower Matrix.
	 *
	 * \see cbind(const Matrix<T,PO1,PS1>& A, 
   *            const Matrix<T,PO2,PS2>& B) 
   *
   * \throw scythe_alloc_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<T,RO,RS>
  rbind (const Matrix<T,PO1,PS1>& A, const Matrix<T,PO2,PS2>& B)
  {
    SCYTHE_CHECK_10(A.cols() != B.cols(), scythe_conformation_error,
        "Matrices have different numbers of columns");

    Matrix<T,RO,Concrete> ret(A.rows() + B.rows(), A.cols(), false);
    std::copy(B.template begin_f<Row>(), B.template end_f<Row>(),
              std::copy(A.template begin_f<Row>(), 
                        A.template end_f<Row>(), 
                        ret.template begin_f<Row>()));
    SCYTHE_VIEW_RETURN(T, RO, RS, ret)
  }

  template <typename T, matrix_order PO1, matrix_style PS1,
            matrix_order PO2, matrix_style PS2>
  Matrix<T,PO1,Concrete>
  rbind (const Matrix<T,PO1,PS1>& A, const Matrix<T,PO2,PS2>& B)
  {
    return rbind<PO1,Concrete>(A, B);
  }

  /* Calculates the order of each element in a Matrix */

  // Functor encapsulating the meat of the algorithm
  namespace {
    template <class T,matrix_order O,matrix_style S> struct order_alg {
      Matrix<T,O> M_;
      order_alg (const Matrix<T,O,S>& M) : M_ (M) {}
      uint operator() (T x) {
        Matrix<bool,O> diff = (M_ < x);
        return std::accumulate(diff.begin_f(), diff.end_f(), (uint) 0);
      }
    };
  }
    
  /*! 
	 * \brief Calculate the rank-order of each element in a Matrix.
	 *
	 * This function calculates the rank-order of each element in a
   * Matrix, returning a Matrix in which the \e i'th element
   * indicates the order position of the \e i'th element of \a M.
   * The returned Matrix contains unsigned integers.
	 * 
	 * \param M A column vector.
   *
   * \throw scythe_alloc_error (Level 1)
	 */

  /* NOTE This function used to only work on column vectors.  I see no
   * reason to maintain this restriction.
   */
  template <matrix_order RO, matrix_style RS, typename T,
            matrix_order PO, matrix_style PS>
  Matrix<unsigned int, RO, RS>
  order (const Matrix<T, PO, PS>& M)
  {
    Matrix<uint, RO, Concrete> ranks(M.rows(), M.cols(), false);
    std::transform(M.begin_f(), M.end_f(), ranks.template begin_f<PO>(),
                   order_alg<T, PO, PS>(M));
    SCYTHE_VIEW_RETURN(uint, RO, RS, ranks)
  }

  template <typename T, matrix_order O, matrix_style S>
  Matrix<unsigned int, O, Concrete>
  order (const Matrix<T,O,S>& M)
  {
    return order<O,Concrete>(M);
  }
  
  /* Selects all the rows of Matrix A for which binary column vector e
   * has an element equal to 1
   */
   
  /*! 
	 * \brief Locate rows for which a binary column vector equals 1
	  
   * This function identifies all the rows of a Matrix \a M for which
   * the binary column vector \a e has an element equal to 1,
   * returning a Matrix 
	 
	 * \param M The Matrix of interest.
	 * \param e A boolean column vector.
	 *
	 * \see unique(const Matrix<T>& M)
   *
   * \throw scythe_conformation_error (Level 1)
   * \throw scythe_dimension_error (Level 1)
   * \throw scythe_alloc_error (Level 1)
	 */
  template <matrix_order RO, matrix_style RS, typename T,