C++ Matrix::rows方法代码示例

    inline Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>
    wishart_rng(const double nu,
                const Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic>& S,
                RNG& rng) {

      static const char* function = "stan::prob::wishart_rng(%1%)";

      using stan::math::check_size_match;
      using stan::math::check_positive;

      check_positive(function,nu,"degrees of freedom",(double*)0);
      check_size_match(function, 
                       S.rows(), "Rows of scale parameter",
                       S.cols(), "columns of scale parameter",
                       (double*)0);

      Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> B(S.rows(), S.cols());
      B.setZero();

      for(int i = 0; i < S.cols(); i++) {
        B(i,i) = std::sqrt(chi_square_rng(nu - i, rng));
        for(int j = 0; j < i; j++)
          B(j,i) = normal_rng(0,1,rng);
      }

      return stan::math::multiply_lower_tri_self_transpose(S.llt().matrixL() * B);
    }

    inline Eigen::Matrix<typename return_type<T1, T2>::type, 
                         Eigen::Dynamic, Eigen::Dynamic>
    append_col(const Eigen::Matrix<T1, R1, C1>& A,
               const Eigen::Matrix<T2, R2, C2>& B) {

      using Eigen::Dynamic;
      using Eigen::Matrix; 
      using stan::math::check_size_match;

      int Arows = A.rows();
      int Brows = B.rows();
      int Acols = A.cols();
      int Bcols = B.cols();
      check_size_match("append_col",
                       "rows of A", Arows, 
                       "rows of B", Brows);
      
      Matrix<typename return_type<T1, T2>::type, Dynamic, Dynamic>
        result(Arows, Acols+Bcols);
      for (int j = 0; j < Acols; j++)
        for (int i = 0; i < Arows; i++)
          result(i, j) = A(i, j);
          
      for (int j = Acols, k = 0; k < Bcols; j++, k++)
        for (int i = 0; i < Arows; i++)
          result(i, j) = B(i, k);
      return result;
    }

    inline 
    Eigen::Matrix<fvar<T>,R1,C2>
    mdivide_right(const Eigen::Matrix<fvar<T>,R1,C1> &A,
                  const Eigen::Matrix<double,R2,C2> &b) {
      
      using stan::math::multiply;      
      using stan::math::mdivide_right;
      stan::math::check_square("mdivide_right(%1%)",b,"b",(double*)0);
      stan::math::check_multiplicable("mdivide_right(%1%)",A,"A",
                                      b,"b",(double*)0);

      Eigen::Matrix<T,R2,C2> deriv_b_mult_inv_b(b.rows(),b.cols());
      Eigen::Matrix<T,R1,C1> val_A(A.rows(),A.cols()); 
      Eigen::Matrix<T,R1,C1> deriv_A(A.rows(),A.cols());

      for (int j = 0; j < A.cols(); j++) {
        for(int i = 0; i < A.rows(); i++) {
          val_A(i,j) = A(i,j).val_;
          deriv_A(i,j) = A(i,j).d_;
        }
      }

      return stan::agrad::to_fvar(mdivide_right(val_A, b), 
                                  mdivide_right(deriv_A, b));
    }

 inline bool check_corr_matrix(const char* function,
               const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
               const char* name,
               T_result* result,
               const Policy&) {
   if (!check_size_match(function, 
       y.rows(), "Rows of correlation matrix",
       y.cols(), "columns of correlation matrix",
       result, Policy())) 
     return false;
   if (!check_positive(function, y.rows(), "rows", result, Policy()))
     return false;
   if (!check_symmetric(function, y, "y", result, Policy()))
     return false;
   for (typename Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>::size_type
          k = 0; k < y.rows(); ++k) {
     if (fabs(y(k,k) - 1.0) > CONSTRAINT_TOLERANCE) {
       std::ostringstream message;
       message << name << " is not a valid correlation matrix. " 
               << name << "(" << k << "," << k 
               << ") is %1%, but should be near 1.0";
       T_result tmp 
         = policies::raise_domain_error<T_y>(function,
                                             message.str().c_str(),
                                             y(k,k), Policy());
       if (result != 0)
         *result = tmp;
       return false;
     }
   }
   if (!check_pos_definite(function, y, "y", result, Policy()))
     return false;
   return true;
 }

    inline 
    fvar<T>
    determinant(const Eigen::Matrix<fvar<T>, R, C>& m) {
      using stan::math::multiply;
      using stan::math::inverse;

      stan::math::check_square("determinant", "m", m);
      Eigen::Matrix<T,R,C> m_deriv(m.rows(), m.cols());
      Eigen::Matrix<T,R,C> m_val(m.rows(), m.cols());

      for(size_type i = 0; i < m.rows(); i++) {
        for(size_type j = 0; j < m.cols(); j++) {
          m_deriv(i,j) = m(i,j).d_;
          m_val(i,j) = m(i,j).val_;
        }
      }

      Eigen::Matrix<T,R,C> m_inv = inverse<T>(m_val);
      m_deriv = multiply(m_inv, m_deriv);

      fvar<T> result;
      result.val_ = m_val.determinant();
      result.d_ = result.val_ * m_deriv.trace();

      // FIXME:  I think this will overcopy compared to retur fvar<T>(...);
      return result;
    }

    inline int
    categorical_rng(const Eigen::Matrix<double, Eigen::Dynamic, 1>& theta,
                    RNG& rng) {
      using boost::variate_generator;
      using boost::uniform_01;

      static const char* function("categorical_rng");

      check_simplex(function, "Probabilities parameter", theta);

      variate_generator<RNG&, uniform_01<> >
        uniform01_rng(rng, uniform_01<>());

      Eigen::VectorXd index(theta.rows());
      index.setZero();

      for (int i = 0; i < theta.rows(); i++) {
        for (int j = i; j < theta.rows(); j++)
          index(j) += theta(i, 0);
      }

      double c = uniform01_rng();
      int b = 0;
      while (c > index(b, 0))
        b++;
      return b + 1;
    }

    typename boost::math::tools::promote_args<T_prob,T_prior_sample_size>::type
    dirichlet_log(const Eigen::Matrix<T_prob,Eigen::Dynamic,1>& theta,
		  const Eigen::Matrix<T_prior_sample_size,Eigen::Dynamic,1>& alpha) {
      static const char* function = "stan::prob::dirichlet_log(%1%)";
      using boost::math::lgamma;
      using boost::math::tools::promote_args;
      using stan::math::check_consistent_sizes;
      using stan::math::check_positive;
      using stan::math::check_simplex;
      using stan::math::multiply_log;
      
      typename promote_args<T_prob,T_prior_sample_size>::type lp(0.0);      
      check_consistent_sizes(function, theta, alpha,
                             "probabilities", "prior sample sizes",
                             &lp);
      check_positive(function, alpha, "prior sample sizes", &lp);
      check_simplex(function, theta, "probabilities", &lp);

      if (include_summand<propto,T_prior_sample_size>::value) {
        lp += lgamma(alpha.sum());
        for (int k = 0; k < alpha.rows(); ++k)
          lp -= lgamma(alpha[k]);
      }
      if (include_summand<propto,T_prob,T_prior_sample_size>::value)
        for (int k = 0; k < theta.rows(); ++k) 
          lp += multiply_log(alpha[k]-1, theta[k]);
      return lp;
    }

    inline bool check_corr_matrix(const std::string& function,
                                  const std::string& name,
                                  const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y) {
      using Eigen::Dynamic;
      using Eigen::Matrix;
      using stan::math::index_type;

      typedef typename index_type<Matrix<T_y,Dynamic,Dynamic> >::type size_t;

      check_size_match(function, 
                       "Rows of correlation matrix", y.rows(), 
                       "columns of correlation matrix", y.cols());
      
      check_positive(function, "rows", y.rows());

      check_symmetric(function, "y", y);
      
      for (size_t k = 0; k < y.rows(); ++k) {
        if (!(fabs(y(k,k) - 1.0) <= CONSTRAINT_TOLERANCE)) {
          std::ostringstream msg;
          msg << "is not a valid correlation matrix. " 
              << name << "(" << stan::error_index::value + k 
              << "," << stan::error_index::value + k 
              << ") is "; ;
          dom_err(function, name, y(k,k),
                  msg.str(),
                  ", but should be near 1.0");
          return false;
        }
      }
      stan::error_handling::check_pos_definite(function, "y", y);
      return true;
    }

void expect_matrix_eq(const Eigen::Matrix<T1,Eigen::Dynamic,Eigen::Dynamic>& a,
                      const Eigen::Matrix<T2,Eigen::Dynamic,Eigen::Dynamic>& b) {
  EXPECT_EQ(a.rows(), b.rows());
  EXPECT_EQ(a.cols(), b.cols());
  for (int i = 0; i < a.rows(); ++i)
    for (int j = 0; j < a.cols(); ++j)
      EXPECT_FLOAT_EQ(value_of(a(i,j)), value_of(b(i,j)));
}

void vector_test(const Eigen::Matrix<Real, NbRows, NbCols>& A, const Eigen::Matrix<Real, NbRows, NbCols>& B, Accumulator& Result)
{
    for(int i = 0, k = 0; i != A.rows() && k != B.rows(); ++i, ++k)
        for(int j = 0, l = 0; j != A.cols() && l != B.cols(); ++j, ++l)
            test(A(i, j), B(k, l), Result);

    Result.exact(A.rows() == B.rows() && A.cols() == B.cols());
};

 inline Eigen::Matrix<fvar<T>,R,1> 
 rows_dot_self(const Eigen::Matrix<fvar<T>,R,C>& x) {
   Eigen::Matrix<fvar<T>,R,1> ret(x.rows(),1);
   for (size_type i = 0; i < x.rows(); i++) {
     Eigen::Matrix<fvar<T>,1,C> crow = x.row(i);
     ret(i,0) = dot_self(crow);
   }
   return ret;
 }

 inline Eigen::Matrix<fvar<T>,R,C>
 divide(const Eigen::Matrix<fvar<T>, R,C>& v, const fvar<T>& c) {
   Eigen::Matrix<fvar<T>,R,C> res(v.rows(),v.cols());
   for(int i = 0; i < v.rows(); i++) {
     for(int j = 0; j < v.cols(); j++)
       res(i,j) = v(i,j) / c;
   }
   return res;
 }

 inline Eigen::Matrix<fvar<typename stan::return_type<T1,T2>::type>,R,C>
 divide(const Eigen::Matrix<T1, R,C>& v, const fvar<T2>& c) {
   Eigen::Matrix<fvar<typename stan::return_type<T1,T2>::type>,R,C> res(v.rows(),v.cols());
   for(int i = 0; i < v.rows(); i++) {
     for(int j = 0; j < v.cols(); j++)
       res(i,j) = to_fvar(v(i,j)) / c;
   }
   return res;
 }

 inline 
 Eigen::Matrix<fvar<T>,R1,C1> 
 multiply(const Eigen::Matrix<fvar<T>, R1, C1>& m, const fvar<T>& c) {
   Eigen::Matrix<fvar<T>,R1,C1> res(m.rows(),m.cols());
   for(int i = 0; i < m.rows(); i++) {
     for(int j = 0; j < m.cols(); j++)
       res(i,j) = c * m(i,j);
   }
   return res;
 }

 inline Eigen::Matrix<double, R1, 1>
 rows_dot_product(const Eigen::Matrix<double, R1, C1>& v1, 
                     const Eigen::Matrix<double, R2, C2>& v2) {
   validate_matching_sizes(v1,v2,"rows_dot_product");
   Eigen::Matrix<double, R1, 1> ret(v1.rows(),1);
   for (size_type j = 0; j < v1.rows(); ++j) {
     ret(j) = v1.row(j).dot(v2.row(j));
   }
   return ret;
 }