matrix.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.
 * --------------------------------------------------------------------
 *  scythe's/matrix.h
 *
 */

/*!
 * \file matrix.h
 * \brief Definitions of Matrix and related classes and functions.
 *
 * This file contains the definitions of the Matrix, Matrix_base, and
 * associated classes.  It also contains a number of external
 * functions that operate on Matrix objects, such as mathematical
 * operators.
 *
 * Many of the arithmetic and logical operators in this file are
 * implemented in terms of overloaded template definitions to provide
 * both generic and default versions of each operation.  Generic
 * templates allow (and force) the user to fully specify the 
 * template type of the returned matrix object (row or column order,
 * concrete or view) while default templates automatically return
 * concrete matrices with the ordering of the first or only Matrix
 * argument to the function.  Furthermore, we overload binary
 * functions to provide scalar by Matrix operations, in addition to
 * basic Matrix by Matrix arithmetic and logic.  Therefore,
 * definitions for multiple versions appear in the functions list
 * below.  We adopt the convention of providing explicit documentation
 * for only the most generic Matrix by Matrix version of each of these
 * operators and describe the behavior of the various overloaded
 * versions in these documents.
 */


#ifndef SCYTHE_MATRIX_H
#define SCYTHE_MATRIX_H

#include <climits>
#include <iostream>
#include <iomanip>
#include <string>
#include <sstream>
#include <fstream>
#include <iterator>
#include <algorithm>
//#include <numeric>
#include <functional>
#include <list>

#ifdef SCYTHE_COMPILE_DIRECT
#include "defs.h"
#include "algorithm.h"
#include "error.h"
#include "datablock.h"
#include "matrix_random_access_iterator.h"
#include "matrix_forward_iterator.h"
#include "matrix_bidirectional_iterator.h"
#ifdef SCYTHE_LAPACK
#include "lapack.h"
#endif
#else
#include "scythestat/defs.h"
#include "scythestat/algorithm.h"
#include "scythestat/error.h"
#include "scythestat/datablock.h"
#include "scythestat/matrix_random_access_iterator.h"
#include "scythestat/matrix_forward_iterator.h"
#include "scythestat/matrix_bidirectional_iterator.h"
#ifdef SCYTHE_LAPACK
#include "scythestat/lapack.h"
#endif
#endif

namespace scythe {

  namespace { // make the uint typedef local to this file
	  /* Convenience typedefs */
	  typedef unsigned int uint;
  }

  /* Forward declare the matrix multiplication functions for *= to use
   * within Matrix proper .
   */

  template <typename T_type, matrix_order ORDER, matrix_style STYLE,
            matrix_order L_ORDER, matrix_style L_STYLE,
            matrix_order R_ORDER, matrix_style R_STYLE>
  Matrix<T_type, ORDER, STYLE>
  operator* (const Matrix<T_type, L_ORDER, L_STYLE>& lhs,
             const Matrix<T_type, R_ORDER, R_STYLE>& rhs);


  template <matrix_order L_ORDER, matrix_style L_STYLE,
            matrix_order R_ORDER, matrix_style R_STYLE, typename T_type>
  Matrix<T_type, L_ORDER, Concrete>
  operator* (const Matrix<T_type, L_ORDER, L_STYLE>& lhs,
             const Matrix<T_type, R_ORDER, R_STYLE>& rhs);

	/* forward declaration of the matrix class */
	template <typename T_type, matrix_order ORDER, matrix_style STYLE>
	class Matrix;

  /*!  \brief A helper class for list-wise initialization.  
   *
   * This class gets used behind the scenes to provide listwise
   * initialization for Matrix objects.  This documentation is mostly
   * intended for developers.
   *
   * The Matrix class's assignment operator returns a ListInitializer
   * object when passed a scalar.  The assignment operator binds before
   * the comma operator, so this happens first, no matter if there is
   * one scalar, or a list of scalars on the right hand side of the
   * assignment sign.  The ListInitializer constructor keeps an iterator
   * to the Matrix that created it and places the initial item at the
   * head of a list.  Then the ListInitializer comma operator gets
   * called 0 or more times, appending items to the list. At this
   * point the ListInitializer object gets destructed because the
   * expression is done and it is just a temporary.  All the action is
   * in the destructor where the list is copied into the Matrix with
   * R-style recycling.
   *
   * To handle chained assignments, such as A = B = C = 1.2 where A,
   * B, and C are matrices, correctly, we encapsulate the Matrix
   * population sequence that is typically called by the destructor in
   * the private function populate, and make Matrix a friend of this
   * class.  The Matrix class contains an assignment operator for
   * ListInitializer objects that calls this function.  When a call
   * like "A = B = C = 1.2" occurs the compiler first evaluates C =
   * 1.2 and returns a ListInitializer object.  Then, the
   * ListInitializer assignment operator in the Matrix class (being
   * called on B = (C = 1.2)) forces the ListInitializer object to
   * populate C early (it would otherwise not occur until destruction
   * at the end of th entire call) by calling the populate method and
   * then does a simple Matrix assignment of B = C and the populated C
   * and the chain of assignment proceeds from there in the usual
   * fashion.
   *
   * Based on code in Blitz++ (http://www.oonumerics.org/blitz/) by
   * Todd Veldhuizen.  Blitz++ is distributed under the terms of the
   * GNU GPL.
   */

  template<typename T_elem, typename T_iter, 
           matrix_order O, matrix_style S>
  class ListInitializer {
    // An unbound friend
    template <class T, matrix_order OO, matrix_style SS>
    friend class Matrix;
    
    public:
      ListInitializer (T_elem val, T_iter begin, T_iter end, 
                       Matrix<T_elem,O,S>* matrix)
        : vals_ (),
          iter_ (begin),
          end_ (end),
          matrix_ (matrix),
          populated_ (false)
      {
        vals_.push_back(val);
      }

      ~ListInitializer ()
      {
        if (! populated_)
          populate();
      }

      ListInitializer &operator, (T_elem x)
      {
        vals_.push_back(x);
        return *this;
      }

    private:
      void populate ()
      {
        typename std::list<T_elem>::iterator vi = vals_.begin();

        while (iter_ < end_) {
          if (vi == vals_.end())
            vi = vals_.begin();
          *iter_ = *vi;
          ++iter_; ++vi;
        }

        populated_ = true;
      }

      std::list<T_elem> vals_;
      T_iter iter_;
      T_iter end_;
      Matrix<T_elem, O, S>* matrix_;
      bool populated_;
  };
	
  /*! \brief Matrix superclass.
   *
   * The Matrix_base class handles Matrix functionality that doesn't
   * need to be templatized with respect to data type.  This helps
   * reduce code bloat by reducing replication of code for member
   * functions that don't rely on templating.  Furthermore, it
   * hides all of the implementation details of indexing.  The
   * constructors of this class are protected and end-users should
   * always work with full-fledged Matrix objects.
   *
   * The public functions in this class generally provide Matrix
   * metadata like information on dimensionality and size.
	 */

  template <matrix_order ORDER = Col, matrix_style STYLE = Concrete>
  class Matrix_base
  {
    protected:
      /**** CONSTRUCTORS ****/

      /* Default constructor */
      Matrix_base ()
        : rows_ (0),
          cols_ (0),
          rowstride_ (0),
          colstride_ (0),
          storeorder_ (ORDER)
      {}

      /* Standard constructor */
      Matrix_base (uint rows, uint cols)
        : rows_ (rows),
          cols_ (cols),
          storeorder_ (ORDER)
      {
        if (ORDER == Col) {
          rowstride_ = 1;
          colstride_ = rows;
        } else {
          rowstride_ = cols;
          colstride_ = 1;
        }
      }

      /* Copy constructors 
       *
       * The first version handles matrices of the same order and
       * style.  The second handles matrices with different
       * orders/styles.  The same- templates are more specific,
       * so they will always catch same- cases.
       *
       */

      Matrix_base (const Matrix_base &m)
        : rows_ (m.rows()),
          cols_ (m.cols()),
          rowstride_ (m.rowstride()),
          colstride_ (m.colstride())
      {
        if (STYLE == View)
          storeorder_ = m.storeorder();
        else
          storeorder_ = ORDER;
      }

      template <matrix_order O, matrix_style S>
      Matrix_base (const Matrix_base<O, S> &m)
        : rows_ (m.rows()),
          cols_ (m.cols())
      {
        if (STYLE == View) {
          storeorder_ = m.storeorder();
          rowstride_ = m.rowstride();
          colstride_ = m.colstride();
         } else {
          storeorder_ = ORDER;
          if (ORDER == Col) {
            rowstride_ = 1;
            colstride_ = rows_;
          } else {
            rowstride_ = cols_;
            colstride_ = 1;
          }
         }
      }

      /* Submatrix constructor */
      template <matrix_order O, matrix_style S>
      Matrix_base (const Matrix_base<O, S> &m,
          uint x1, uint y1, uint x2, uint y2)
        : rows_ (x2 - x1 + 1),
          cols_ (y2 - y1 + 1),
          rowstride_ (m.rowstride()),
          colstride_ (m.colstride()),
          storeorder_ (m.storeorder())
      {
        /* 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.  Of course, we should never get this far.
         */
        if (STYLE == Concrete) {
          SCYTHE_THROW(scythe_style_error, 
              "Tried to construct a concrete submatrix (Matrix_base)!");
        }
      }


      /**** DESTRUCTOR ****/

      ~Matrix_base ()
      {}

      /**** OPERRATORS ****/

      // I'm defining one just to make sure we don't get any subtle
      // bugs from defaults being called.
      Matrix_base& operator=(const Matrix_base& m)
      {
        SCYTHE_THROW(scythe_unexpected_default_error,
          "Unexpected call to Matrix_base default assignment operator");
      }

      /**** MODIFIERS ****/

			/* Make this Matrix_base an exact copy of the matrix passed in.
			 * Just like an assignment operator but can be called from a derived
			 * class w/out making the = optor public and doing something
			 * like
			 * *(dynamic_cast<Matrix_base *> (this)) = M;
			 * in the derived class.
       *
       * Works across styles, but should be used with care
			 */
      template <matrix_order O, matrix_style S>
			inline void mimic (const Matrix_base<O, S> &m)
			{
				rows_ = m.rows();
				cols_ = m.cols();
				rowstride_ = m.rowstride();
				colstride_ = m.colstride();
        storeorder_ = m.storeorder();
			}

      /* Reset the dimensions of this Matrix_base.
       *
       * TODO This function is a bit of an interface weakness.  It
       * assumes a resize always means a fresh matrix (concrete or
       * view) with no slack between dims and strides. This happens to
       * always be the case when it is called, but it tightly couples
       * Matrix_base and extending classes.  Not a big issue (since
       * Matrix is probably the only class that will ever extend this)
       * but something to maybe fix down the road.
       */
			inline void resize (uint rows, uint cols)
			{
				rows_ = rows;
				cols_ = cols;

        if (ORDER == Col) {
          rowstride_ = 1;
          colstride_ = rows;
        } else {
          rowstride_ = cols;
          colstride_ = 1;
        }

        storeorder_ = ORDER;
			}