ide.h

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

   * matrix to invert, ahead of time.
   *
   * \param A The symmetric positive definite matrix to invert.
   * \param M The lower triangular matrix from the Cholesky decomposition of \a A.
   *
   * \see invpd(const Matrix<T, PO, PS>&)
   * \see inv(const Matrix<T,PO1,PS1>&, const Matrix<T,PO2,PS2>&, const Matrix<T,PO3,PS3>&, const Matrix<unsigned int,PO4,PS4>&)
   * \see inv(const Matrix<T, PO, PS>&)
   * \see cholesky(const Matrix<T, PO, PS>&)
   *
   * \throw scythe_alloc_error (Level 1)
   * \throw scythe_null_error (Level 1)
   * \throw scythe_conformation_error (Level 1)
   * \throw scythe_dimension_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>
  invpd (const Matrix<T,PO1,PS1>& A, const Matrix<T,PO2,PS2>& M)
  {
    SCYTHE_CHECK_10(A.isNull(), scythe_null_error,
        "A is NULL")
    SCYTHE_CHECK_10(! A.isSquare(), scythe_dimension_error,
        "A is not square")
    SCYTHE_CHECK_10(A.rows() != M.cols() || A.cols() != M.rows(), 
        scythe_conformation_error, "A and M do not conform");
      
    // for chol_solve block
    T *y = new T[A.rows()];
    T *x = new T[A.rows()];
    Matrix<T, RO, Concrete> b(A.rows(), 1); // full of zeros
    Matrix<T, RO, Concrete> null;
    
    // For final answer
    Matrix<T, RO, Concrete> Ainv(A.rows(), A.cols(), false);

    for (uint k = 0; k < A.rows(); ++k) {
      b[k] = (T) 1;

      solve(M, null, b, x, y);

      b[k] = (T) 0;
      for (uint l = 0; l < A.rows(); ++l)
        Ainv(l,k) = x[l];
    }

    delete[] y;
    delete[] x;

    SCYTHE_VIEW_RETURN(T, RO, RS, Ainv)
  }

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

  /*!\brief Calculate the inverse of a symmetric positive definite
   * matrix.
  *
  * This function returns the inverse of a symmetric positive definite
  * matrix, using cholesky() to do the necessary decomposition. This
  * method is significantly faster than the generalized inverse
  * function.
  *
  * \param A The symmetric positive definite matrix to invert.
  *
  * \see invpd(const Matrix<T, PO1, PS1>&, const Matrix<T, PO2, PS2>&)
  * \see inv (const Matrix<T,PO1,PS1>&, const Matrix<T,PO2,PS2>&, const Matrix<T,PO3,PS3>&, const Matrix<unsigned int,PO4,PS4>&)
  * \see inv (const Matrix<T, PO, PS>&)
  *
  * \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)
  */
  template <matrix_order RO, matrix_style RS, typename T,
            matrix_order PO, matrix_style PS>
  Matrix<T, RO, RS>
  invpd (const Matrix<T, PO, PS>& A)
  { 
    // Cholesky checks to see if A is square and symmetric
  
    return invpd<RO,RS>(A, cholesky<RO,Concrete>(A));
  }

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

  /* This code is based on  Algorithm 3.4.1 of Golub and Van Loan 3rd
   * edition, 1996. Major difference is in how the output is
   * structured.  Returns the sign of the row permutation (used by
   * det).  Internal function, doesn't need doxygen.
   */
  namespace {
    template <matrix_order PO1, matrix_style PS1, typename T,
              matrix_order PO2, matrix_order PO3, matrix_order PO4>
    inline T
    lu_decomp_alg(Matrix<T,PO1,PS1>& A, Matrix<T,PO2,Concrete>& L, 
                  Matrix<T,PO3,Concrete>& U, 
                  Matrix<unsigned int, PO4, Concrete>& perm_vec)
    {
      if (A.isRowVector()) {
        L = Matrix<T,PO2,Concrete> (1, 1, true, 1); // all 1s
        U = A;
        perm_vec = Matrix<uint, PO4, Concrete>(1, 1);  // all 0s
        return (T) 0;
      }
      
      L = U = Matrix<T, PO2, Concrete>(A.rows(), A.cols(), false);
      perm_vec = Matrix<uint, PO3, Concrete> (A.rows() - 1, 1, false);

      uint pivot;
      T temp;
      T sign = (T) 1;

      for (uint k = 0; k < A.rows() - 1; ++k) {
        pivot = k;
        // find pivot
        for (uint i = k; i < A.rows(); ++i) {
          if (std::fabs(A(pivot,k)) < std::fabs(A(i,k)))
            pivot = i;
        }
        
        SCYTHE_CHECK_20(A(pivot,k) == (T) 0, scythe_type_error,
            "Matrix is singular");

        // permute
        if (k != pivot) {
          sign *= -1;
          for (uint i = 0; i < A.rows(); ++i) {
            temp = A(pivot,i);
            A(pivot,i) = A(k,i);
            A(k,i) = temp;
          }
        }
        perm_vec[k] = pivot;

        for (uint i = k + 1; i < A.rows(); ++i) {
          A(i,k) = A(i,k) / A(k,k);
          for (uint j = k + 1; j < A.rows(); ++j)
            A(i,j) = A(i,j) - A(i,k) * A(k,j);
        }
      }

      L = A;

      for (uint i = 0; i < A.rows(); ++i) {
        for (uint j = i; j < A.rows(); ++j) {
          U(i,j) = A(i,j);
          L(i,j) = (T) 0;
          L(i,i) = (T) 1;
        }
      }
      return sign;
    }
  }

  /* Calculates the LU Decomposition of a square Matrix */

  /* Note that the L, U, and perm_vec must be concrete. A is passed by
   * value, because it is changed during the decomposition.  If A is a
   * view, it will get mangled, but the decomposition will work fine.
   * Not sure what the copy/view access trade-off is, but passing a
   * view might speed things up if you don't care about messing up
   * your matrix.
   */
    /*! \brief LU decomposition of a square matrix.
     *
     * This function performs LU decomposition. That is, given a
     * non-singular square matrix \a A and three matrix references, \a
     * L, \a U, and \a perm_vec, lu_decomp fills the latter three
     * matrices such that \f$LU = A\f$. This method does not actually
     * calculate the LU decomposition of \a A, but of a row-wise
     * permutation of \a A. This permutation is recorded in perm_vec.
     *
     * \note Note that \a L, \a U, and \a perm_vec must be concrete.
     * \a A is passed by value because the function modifies it during
     * the decomposition.  Users should generally avoid passing Matrix
     * views as the first argument to this function because this
     * results in modification to the Matrix being viewed.
     *
     * \param A Non-singular square matrix to decompose.
     * \param L Lower triangular portion of LU decomposition of A.
     * \param U Upper triangular portion of LU decomposition of A.
     * \param perm_vec Permutation vector recording the row-wise permutation of A actually decomposed by the algorithm.
     *
     * \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_null_error (Level 1)
     * \throw scythe_dimension_error (Level 1)
     * \throw scythe_type_error (Level 2)
     */
  template <matrix_order PO1, matrix_style PS1, typename T,
            matrix_order PO2, matrix_order PO3, matrix_order PO4>
  void
  lu_decomp(Matrix<T,PO1,PS1> A, Matrix<T,PO2,Concrete>& L, 
            Matrix<T,PO3,Concrete>& U, 
            Matrix<unsigned int, PO4, Concrete>& perm_vec)
  {
    SCYTHE_CHECK_10(A.isNull(), scythe_null_error,
        "A is NULL")
    SCYTHE_CHECK_10(! A.isSquare(), scythe_dimension_error,
        "Matrix A not square");

    lu_decomp_alg(A, L, U, perm_vec);
  }

  /* lu_solve overloaded: you need A, b + L, U, perm_vec from
   * lu_decomp.
   *
   */
    /*! \brief Solve \f$Ax=b\f$ for x via forward and backward
     * substitution, given the results of a LU decomposition.
     *
     * This function solves the system of equations \f$Ax = b\f$ via
     * forward and backward substitution and LU decomposition. \a A
     * must be a non-singular square matrix for this method to work.
     * This function requires the actual LU decomposition to be
     * performed ahead of time; by lu_decomp() for example.
     *
     * 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 Non-singular square Matrix to decompose, passed by reference.
     * \param b Column vector with as many rows as \a A.
     * \param L Lower triangular portion of LU decomposition of \a A.
     * \param U Upper triangular portion of LU decomposition of \a A.
     * \param perm_vec Permutation vector recording the row-wise permutation of \a A actually decomposed by the algorithm.
     *
     * \see lu_solve (Matrix<T,PO1,PS1>, const Matrix<T,PO2,PS2>&)
     * \see lu_decomp(Matrix<T,PO1,PS1>, Matrix<T,PO2,Concrete>&, Matrix<T,PO3,Concrete>&, Matrix<unsigned int, PO4, Concrete>&)
   * \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> &)
     *
     * \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_order PO4, matrix_style PS4,
            matrix_order PO5, matrix_style PS5>
  Matrix<T, RO, RS>
  lu_solve (const Matrix<T,PO1,PS1>& A, const Matrix<T,PO2,PS2>& b,
            const Matrix<T,PO3,PS3>& L, const Matrix<T,PO4,PS4>& U,
            const Matrix<unsigned int, PO5, PS5> &perm_vec) 
  {
    SCYTHE_CHECK_10(A.isNull(), scythe_null_error,
        "A is NULL")
    SCYTHE_CHECK_10(! b.isColVector(), scythe_dimension_error,
        "b is not a column vector");
    SCYTHE_CHECK_10(! A.isSquare(), scythe_dimension_error,
        "A is not square");
    SCYTHE_CHECK_10(A.rows() != b.rows(), scythe_conformation_error,
        "A and b have different row sizes");
    SCYTHE_CHECK_10(A.rows() != L.rows() || A.rows() != U.rows() ||
                    A.cols() != L.cols() || A.cols() != U.cols(),
                    scythe_conformation_error,
                    "A, L, and U do not conform");
    SCYTHE_CHECK_10(perm_vec.rows() + 1 != A.rows(),
        scythe_conformation_error,
        "perm_vec does not have exactly one less row than A");

    T *y = new T[A.rows()];
    T *x = new T[A.rows()];
    
    Matrix<T,RO,Concrete> bb = row_interchange(b, perm_vec);
    solve(L, U, bb, 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_order PO4, matrix_style PS4,
            matrix_order PO5, matrix_style PS5>
  Matrix<T, PO1, Concrete>
  lu_solve (const Matrix<T,PO1,PS1>& A, const Matrix<T,PO2,PS2>& b,
            const Matrix<T,PO3,PS3>& L, const Matrix<T,PO4,PS4>& U,
            const Matrix<unsigned int, PO5, PS5> &perm_vec) 
  {
    return lu_solve<PO1,Concrete>(A, b, L, U, perm_vec);
  }

    /*! \brief Solve \f$Ax=b\f$ for x via forward and backward
     * substitution, using LU decomposition 
     *
     * This function solves the system of equations \f$Ax = b\f$ via
     * forward and backward substitution and LU decomposition. \a A
     * must be a non-singular square matrix for this method to work.
     *
     * \param A A non-singular square Matrix to decompose.
     * \param b A column vector with as many rows as \a A.
     *
     * \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_decomp(Matrix<T,PO1,PS1>, Matrix<T,PO2,Concrete>&, Matrix<T,PO3,Concrete>&, Matrix<unsigned int, PO4, Concrete>&)
   * \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> &)
   *
     * \throw scythe_null_error (Level 1)
     * \throw scythe_dimension_error (Level 1)
     * \throw scythe_conformation_error (Level 1)
     * \throw scythe_type_error (Level 2)
     */
  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>
  lu_solve (Matrix<T,PO1,PS1> A, const Matrix<T,PO2,PS2>& b)
  {
    // step 1 compute the LU factorization 
    Matrix<T, RO, Concrete> L, U;
    Matrix<uint, RO, Concrete> perm_vec;
    lu_decomp_alg(A, L, U, perm_vec);

    return lu_solve<RO,RS>(A, b, L, U, perm_vec);
  }

  template <typename T, matrix_order PO1, matrix_style PS1,
            matrix_order PO2, matrix_style PS2>
  Matrix<T,PO1,Concrete>
  lu_solve (Matrix<T,PO1,PS1> A, const Matrix<T,PO2,PS2>& b)
  {
    // Slight code rep here, but very few lines
   
    // step 1 compute the LU factorization 
    Matrix<T, PO1, Concrete> L, U;
    Matrix<uint, PO1, Concrete> perm_vec;
    lu_decomp_alg(A, L, U, perm_vec);
    
    return lu_solve<PO1,Concrete>(A, b, L, U, perm_vec);