matrix.h

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

			 * referenced by the input iterator \a it.  Pointers are a form
			 * of input iterator, so one can use this constructor to
			 * initialize a matrix object from a c-style array.  The caller
			 * is responsible for supplying an iterator that won't be
			 * exhausted too soon.
       *
       * \param rows The number of rows in the Matrix.
       * \param cols The number of columns in the Matrix.
       * \param it The input iterator to read from.
       *
       * \see Matrix()
       * \see Matrix(T_type)
       * \see Matrix(uint, uint, bool, T_type)
       * \see Matrix(const std::string&)
       * \see Matrix(const Matrix&)
			 * \see Matrix(const Matrix<T_type, O, S> &)
       * \see Matrix(const Matrix<S_type, O, S> &)
       * \see Matrix(const Matrix<T_type, O, S>&, uint, uint, uint, uint)
       *
       * \throw scythe_alloc_error (Level 1)
       *
       * \b Example:
       * \include example.matrix.constructor.iterator.cc
       */
			template <typename T_iterator>
			Matrix (uint rows, uint cols, T_iterator it)
				:	Base (rows, cols),
					DBRef (rows * cols)
			{
        // TODO again, should probably use iterator
				for (uint i = 0; i < Base::size(); ++i) {
					data_[Base::index(i)] = *it; // we know data_ is contig
					++it;
				}
			}

      /*! \brief File constructor.
       *
       * Constructs a matrix from the contents of a file.  The
       * standard input file format is a simple rectangular text file
       * with one matrix row per line and spaces delimiting values in
       * a row.  Optionally, one can also use Scythe's old file format
       * which is a space-delimited, row-major ordered, list of values
       * with row and column lengths in the first two slots.
       *
       * \param path The path of the input file.
       * \param oldstyle Whether or not to use Scythe's old file format.
       *
       * \see Matrix()
       * \see Matrix(T_type)
       * \see Matrix(uint, uint, bool, T_type)
       * \see Matrix(uint, uint, T_iterator)
       * \see Matrix(const Matrix&)
			 * \see Matrix(const Matrix<T_type, O, S> &)
       * \see Matrix(const Matrix<S_type, O, S> &)
       * \see Matrix(const Matrix<T_type, O, S>&, uint, uint, uint, uint)
       * \see save(const std::string&)
       *
       * \throw scythe_alloc_error (Level 1)
       * \throw scythe_file_error (Level 1)
       * \throw scythe_bounds_error (Level 3)
       *
       * \b Example:
       * \include example.matrix.constructor.file.cc
       */
      Matrix (const std::string& path, bool oldstyle=false)
        : Base (),
          DBRef ()
      {
        std::ifstream file(path.c_str());
        SCYTHE_CHECK_10(! file, scythe_file_error,
            "Could not open file " << path);

        if (oldstyle) {
          uint rows, cols;
          file >> rows >> cols;
          resize(rows, cols);
          std::copy(std::istream_iterator<T_type> (file), 
                    std::istream_iterator<T_type>(), begin_f<Row>());
        } else {
          std::string row;

          unsigned int cols = -1;
          std::vector<std::vector<T_type> > vals;
          unsigned int rows = 0;
          while (std::getline(file, row)) {
            std::vector<T_type> column;
            std::istringstream rowstream(row);
            std::copy(std::istream_iterator<T_type> (rowstream),
                 std::istream_iterator<T_type>(),
                 std::back_inserter(column));

            if (cols == -1)
              cols = (unsigned int) column.size();

            SCYTHE_CHECK_10(cols != column.size(), scythe_file_error,
                "Row " << (rows + 1) << " of input file has "
                << column.size() << " elements, but should have " << cols);

            vals.push_back(column);
            rows++;
          }

          resize(rows, cols);
          for (unsigned int i = 0; i < rows; ++i)
            operator()(i, _) = Matrix<T_type>(1, cols, vals[i].begin());
        }
      }

      /* Copy constructors. Uses template args to set up correct
       * behavior for both concrete and view matrices.  The branches
       * are no-ops and get optimized away at compile time.
       *
       * We have to define this twice because we must explicitly
       * override the default copy constructor; otherwise it is the
       * most specific template in a lot of cases and causes ugliness.
       */

      /*! \brief Default copy constructor.
       *
       * Copy constructing a concrete Matrix makes an exact copy of M
       * in a new data block.  On the other hand, copy constructing a
       * view Matrix generates a new Matrix object that references (or
       * views) M's existing data block.
       *
       * \param M The Matrix to copy or make a view of.
       *
       * \see Matrix()
       * \see Matrix(T_type)
       * \see Matrix(uint, uint, bool, T_type)
       * \see Matrix(uint, uint, T_iterator)
       * \see Matrix(const std::string&)
			 * \see Matrix(const Matrix<T_type, O, S> &)
       * \see Matrix(const Matrix<S_type, O, S> &)
       * \see Matrix(const Matrix<T_type, O, S>&, uint, uint, uint, uint)
       * \see copy()
       * \see copy(const Matrix<T_type, O, S> &)
       * \see reference(const Matrix<T_type, O, S> &)
       *
       * \throw scythe_alloc_error (Level 1)
       *
       * \b Example:
       * \include example.matrix.constructor.copy.cc
       */

      Matrix (const Matrix& M)
				:	Base (M), // this call deals with concrete-view conversions
					DBRef ()
			{
        if (STYLE == Concrete) {
          referenceNew(M.size());
          scythe::copy<ORDER,ORDER>(M, *this);
        } else // STYLE == View
          referenceOther(M);
			}

      /*! \brief Cross order and/or style copy constructor.
       *
       * Copy constructing a concrete Matrix makes an exact copy of M
       * in a new data block.  On the other hand, copy constructing a
       * view Matrix generates a new Matrix object that references (or
       * views) M's existing data block.
       *
       * This version of the copy constructor extends Matrix(const
       * Matrix &) by allowing the user to make concrete copies and
       * views of matrices that have matrix_order or matrix_style that
       * does not match that of the constructed Matrix.  That is, this
       * constructor makes it possible to create views of concrete
       * matrices and concrete copies of views, row-major copies of
       * col-major matrices, and so on.
       *
       * \param M The Matrix to copy or make a view of.
       *
       * \see Matrix()
       * \see Matrix(T_type)
       * \see Matrix(uint, uint, bool, T_type)
       * \see Matrix(uint, uint, T_iterator)
       * \see Matrix(const std::string&)
       * \see Matrix(const Matrix&)
       * \see Matrix(const Matrix<S_type, O, S> &)
       * \see Matrix(const Matrix<T_type, O, S>&, uint, uint, uint, uint)
       * \see copy()
       * \see copy(const Matrix<T_type, O, S> &)
       * \see reference(const Matrix<T_type, O, S> &)
       *
       * \throw scythe_alloc_error (Level 1)
       *
       * \b Example:
       * \include example.matrix.constructor.crosscopy.cc
       */

      template <matrix_order O, matrix_style S>
			Matrix (const Matrix<T_type, O, S> &M)
				:	Base (M), // this call deals with concrete-view conversions
					DBRef ()
			{
        if (STYLE == Concrete) {
          referenceNew(M.size());
          scythe::copy<ORDER,ORDER> (M, *this);
        } else // STYLE == View
          referenceOther(M);
			}

      /*! \brief Cross type copy constructor
       *
       * The type conversion copy constructor takes a reference to an
       * existing matrix containing elements of a different type than
       * the constructed matrix and creates a copy. This constructor
       * will only work if it is possible to cast elements from the
       * copied matrix to the type of elements in the constructed
       * matrix.
       *
       * This constructor always creates a deep copy of the existing
       * matrix, even if the constructed matrix is a view. It is
       * impossible for a matrix view with one element type to
       * reference the data block of a matrix containing elements of a
       * different type. 
       * 
       * \param M The Matrix to copy.
       *
       * \see Matrix()
       * \see Matrix(T_type)
       * \see Matrix(uint, uint, bool, T_type)
       * \see Matrix(uint, uint, T_iterator)
       * \see Matrix(const std::string&)
       * \see Matrix(const Matrix&)
			 * \see Matrix(const Matrix<T_type, O, S> &)
       * \see Matrix(const Matrix<T_type, O, S>&, uint, uint, uint, uint)
       *
       * \throw scythe_alloc_error (Level 1)
       *
       * \b Example:
       * \include example.matrix.constructor.convcopy.cc
       */
			template <typename S_type, matrix_order O, matrix_style S>
			Matrix (const Matrix<S_type, O, S> &M)
				:	Base(M), // this call deal with concrete-view conversions
					DBRef (M.size())
			{
        scythe::copy<ORDER,ORDER> (M, *this);
			}

      /*! \brief Submatrix constructor
       *
       * The submatrix constructor takes a reference to an existing
       * matrix and a set of indices, and generates a new Matrix
       * object referencing the submatrix described by the indices.
       * One can only construct a submatrix with a view template and
       * this constructor will throw an error if one tries to use it
       * to construct a concrete matrix.
       *
       * \note
       * The submatrix-returning operators provide the same
       * functionality as this constructor with less obtuse syntax.
       * Users should generally employ these methods instead of this
       * constructor.
       *
       * \param M  The Matrix to view.
       * \param x1 The first row coordinate, \a x1 <= \a x2.
       * \param y1 The first column coordinate, \a y1 <= \a y2.
       * \param x2 The second row coordinate, \a x2 > \a x1.
       * \param y2 The second column coordinate, \a y2 > \a y1.
       *
       * \see Matrix()
       * \see Matrix(T_type)
       * \see Matrix(uint, uint, bool, T_type)
       * \see Matrix(uint, uint, T_iterator)
       * \see Matrix(const std::string&)
       * \see Matrix(const Matrix&)
			 * \see Matrix(const Matrix<T_type, O, S> &)
       * \see Matrix(const Matrix<S_type, O, S> &)
       * \see operator()(uint, uint, uint, uint)
       * \see operator()(uint, uint, uint, uint) const
       * \see operator()(all_elements, uint)
       * \see operator()(all_elements, uint) const
       * \see operator()(uint, all_elements)
       * \see operator()(uint, all_elements) const
       * \see reference(const Matrix<T_type, O, S> &)
       *
       * \throw scythe_style_error (Level 0)
       * \throw scythe_alloc_error (Level 1)
       */
      template <matrix_order O, matrix_style S>
      Matrix (const Matrix<T_type, O, S> &M,
          uint x1, uint y1, uint x2, uint y2)
        : Base(M, x1, y1, x2, y2),
          DBRef(M, Base::index(x1, y1))
      {
        /* Submatrices always have to be views, but the whole
         * concrete-view thing is just a policy maintained by the
         * software.  Therefore, we need to ensure this constructor
         * only returns views.  There's no neat way to do it but this
         * is still a compile time check, even if it only reports at
         * run-time.
         */
        if (STYLE == Concrete) {
          SCYTHE_THROW(scythe_style_error, 
              "Tried to construct a concrete submatrix (Matrix)!");
        }
      }

    public:
      /**** DESTRUCTOR ****/

      /*!\brief Destructor. 
       */
      ~Matrix() {}

      /**** COPY/REFERENCE METHODS ****/

			/* Make this matrix a view of another's data. If this matrix's
			 * previous datablock is not viewed by any other object it is
			 * deallocated.  Concrete matrices cannot be turned into views
       * at run-time!  Therefore, we generate an error here if *this
       * is concrete.
			 */

      /*!\brief View another Matrix's data.
       *
       * This modifier makes this matrix a view of another's data.
       * The action detaches the Matrix from its current view; if no
       * other Matrix views the detached DataBlock, it will be
       * deallocated.  
       *
       * Concrete matrices cannot convert into views at
       * run time.  Therefore, it is an error to invoke this method on
       * a concrete Matrix.
       *
       * \param M The Matrix to view.
       *
       * \see Matrix(const Matrix&)
			 * \see Matrix(const Matrix<T_type, O, S> &)
       * \see Matrix(const Matrix<S_type, O, S> &)
       * \see copy()
       * \see copy(const Matrix<T_type, O, S> &)
       *
       * \throw scythe_style_error (Level 0)
       *
       * \b Example:
       * \include example.matrix.reference.cc
       */
      template <matrix_order O, matrix_style S>
			inline void reference (const Matrix<T_type, O, S> &M)
			{
        if (STYLE == Concrete) {
          SCYTHE_THROW(scythe_style_error, 
              "Concrete matrices cannot reference other matrices");
        } else {
          referenceOther(M);
          mimic(M);
        }
			}

      /*!\brief Create a copy of this matrix.
       *
       * Creates a deep copy of this Matrix.  The returned concrete
       * matrix references a newly created DataBlock that contains
       * values that are identical to, but distinct from, the values
       * contained in the original Matrix.
       *
       * \see Matrix(const Matrix&)
			 * \see Matrix(const Matrix<T_type, O, S> &)
       * \see Matrix(const Matrix<S_type, O, S> &)
       * \see copy(const Matrix<T_type, O, S> &)
       * \see reference(const Matrix<T_type, O, S> &)
       *
       * \throw scythe_alloc_error (Level 1)
       *
       * \b Example:
       * \include example.matrix.copy.cc
       */
			inline Matrix<T_type, ORDER, Concrete> copy () const
			{
				Matrix<T_type, ORDER> res (Base::rows(), Base::cols(), false);
				std::copy(begin_f(), end_f(), res.begin_f());