C++ MatrixX::row方法代码示例

GMMExpectationMaximization::uint GMMExpectationMaximization::execute(const MatrixX & dataset)
{  
  const uint data_count = dataset.rows();
  const uint num_gaussians = m_means.size();
  const uint dim = dataset.cols();

  MatrixX pxi(data_count,num_gaussians);
  MatrixX pix(data_count,num_gaussians);
  VectorX pxidatatot(data_count);
  VectorX weights(num_gaussians);
  VectorX ex(data_count);
  MatrixX ts(dim,dim);
  VectorX dif(dim);

  Real prev_log_likelyhood = 1.0;
  
  uint it_num;
  for (it_num = 0; it_num < m_max_iterations; it_num++)
  {
    for (uint g = 0; g < num_gaussians; g++)
      weights[g] = m_weights[g];

    for (uint d = 0; d < data_count; d++)
      for (uint g = 0; g < num_gaussians; g++)
        pxi(d,g) = gauss(m_means[g],m_covs[g],dataset.row(d).transpose());

    pxidatatot = pxi * weights;
    Real log_likelyhood = pxidatatot.array().log().sum() / Real(data_count);

    if (it_num != 0 && (std::abs(log_likelyhood / prev_log_likelyhood - 1.0) < m_termination_threshold))
      break;
    prev_log_likelyhood = log_likelyhood;

    for (uint d = 0; d < data_count; d++)
      pix.row(d) = (pxi.row(d).transpose().array() * weights.array()).transpose() / pxidatatot[d];
    
    ex = pix.colwise().sum();

    for(uint g = 0; g < num_gaussians; g++)
    {
      m_weights[g] = ex[g] / Real(data_count);

      m_means[g] = (dataset.transpose() * pix.col(g)) / ex[g];

      ts = MatrixX::Zero(dim,dim);
      for (uint d = 0; d < data_count; d++)
      {
        dif = dataset.row(d).transpose() - m_means[g];
        ts.noalias() += (dif * dif.transpose()) * pix(d,g);
      }
      m_covs[g] = (ts / ex[g]) + MatrixX::Identity(dim,dim) * m_epsilon;
    }

    // interruption point here
    if (m_termination_handler && m_termination_handler->isTerminated())
      return it_num;
  }

  return it_num;
}

void run_test(int dim, int num_elements)
{
  using std::abs;
  typedef typename internal::traits<MatrixType>::Scalar Scalar;
  typedef Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> MatrixX;
  typedef Matrix<Scalar, Eigen::Dynamic, 1> VectorX;

  // MUST be positive because in any other case det(cR_t) may become negative for
  // odd dimensions!
  const Scalar c = abs(internal::random<Scalar>());

  MatrixX R = randMatrixSpecialUnitary<Scalar>(dim);
  VectorX t = Scalar(50)*VectorX::Random(dim,1);

  MatrixX cR_t = MatrixX::Identity(dim+1,dim+1);
  cR_t.block(0,0,dim,dim) = c*R;
  cR_t.block(0,dim,dim,1) = t;

  MatrixX src = MatrixX::Random(dim+1, num_elements);
  src.row(dim) = Matrix<Scalar, 1, Dynamic>::Constant(num_elements, Scalar(1));

  MatrixX dst = cR_t*src;

  MatrixX cR_t_umeyama = umeyama(src.block(0,0,dim,num_elements), dst.block(0,0,dim,num_elements));

  const Scalar error = ( cR_t_umeyama*src - dst ).norm() / dst.norm();
  VERIFY(error < Scalar(40)*std::numeric_limits<Scalar>::epsilon());
}

void run_fixed_size_test(int num_elements)
{
  using std::abs;
  typedef Matrix<Scalar, Dimension+1, Dynamic> MatrixX;
  typedef Matrix<Scalar, Dimension+1, Dimension+1> HomMatrix;
  typedef Matrix<Scalar, Dimension, Dimension> FixedMatrix;
  typedef Matrix<Scalar, Dimension, 1> FixedVector;

  const int dim = Dimension;

  // MUST be positive because in any other case det(cR_t) may become negative for
  // odd dimensions!
  const Scalar c = abs(internal::random<Scalar>());

  FixedMatrix R = randMatrixSpecialUnitary<Scalar>(dim);
  FixedVector t = Scalar(50)*FixedVector::Random(dim,1);

  HomMatrix cR_t = HomMatrix::Identity(dim+1,dim+1);
  cR_t.block(0,0,dim,dim) = c*R;
  cR_t.block(0,dim,dim,1) = t;

  MatrixX src = MatrixX::Random(dim+1, num_elements);
  src.row(dim) = Matrix<Scalar, 1, Dynamic>::Constant(num_elements, Scalar(1));

  MatrixX dst = cR_t*src;

  Block<MatrixX, Dimension, Dynamic> src_block(src,0,0,dim,num_elements);
  Block<MatrixX, Dimension, Dynamic> dst_block(dst,0,0,dim,num_elements);

  HomMatrix cR_t_umeyama = umeyama(src_block, dst_block);

  const Scalar error = ( cR_t_umeyama*src - dst ).array().square().sum();

  VERIFY(error < Scalar(10)*std::numeric_limits<Scalar>::epsilon());
}

GMMExpectationMaximization::Real GMMExpectationMaximization::getBIC(const MatrixX & dataset) const
{
  const uint dim = dataset.cols();
  const uint num_gaussians = m_means.size();

  Real number_of_parameters = (num_gaussians * dim * (dim + 1) / 2) + num_gaussians * dim + num_gaussians - 1;

  uint data_count = dataset.rows();
  Real sum = 0.0;
  
  for(uint i = 0; i < data_count; i++)
    sum += log(expectation(dataset.row(i).transpose()));

  return -sum + (number_of_parameters / 2.0) * log(Real(data_count));
}

void GMMExpectationMaximization::autoInitializeByEqualIntervals(uint num_gaussians,uint col,const MatrixX & dataset)
{
  uint data_count = dataset.rows();
  uint dim = dataset.cols();

  std::vector<std::vector<uint> > index(num_gaussians);
  for(uint g = 0; g < num_gaussians; g++)
    index[g].reserve(data_count / num_gaussians);

  m_weights.clear();
  m_weights.resize(num_gaussians);
  m_means.clear();
  m_means.resize(num_gaussians,VectorX::Zero(dim));
  m_covs.clear();
  m_covs.resize(num_gaussians,MatrixX::Zero(dim,dim));

  // find max and min value for column col
  Real cmax = dataset(0,col);
  Real cmin = dataset(0,col);
  for(uint n = 1; n < data_count; n++)
  {
    if (dataset(n,col) > cmax) cmax = dataset(n,col);
    if (dataset(n,col) < cmin) cmin = dataset(n,col);
  }
  Real cspan = cmax - cmin;

  for(uint n = 0; n < data_count; n++) 
  {
    // compute gaussian index to which this point belongs
    uint gi = uint((dataset(n,col) - cmin) / (cspan + 1.0) * Real(num_gaussians));

    // sum the points to obtain means
    m_means[gi] += dataset.row(n);

    index[gi].push_back(n);
  }

  for (uint g = 0; g < num_gaussians; g++)
  {
    uint popsize = index[g].size();
    // avoid division by zero: if no samples are available, initialize to something from somewhere
    if (popsize == 0)
    {
      m_means[g] = dataset.row(g % data_count);
      m_covs[g] = MatrixX::Identity(dim,dim);
      m_weights[g] = 1.0f / Real(num_gaussians);
      continue;
    }

    // average by popsize
    m_means[g] /= Real(popsize);
    // same weight for all gaussians
    m_weights[g] = 1.0f / Real(num_gaussians);
     
    // compute covariance matrix
    for (uint p = 0; p < popsize; p++)
    {
      const Eigen::VectorXf & r = dataset.row(index[g][p]);
      const Eigen::VectorXf & m = m_means[g];
      m_covs[g] += (r - m) * (r - m).transpose();
    }

    m_covs[g] /= Real(popsize);
    m_covs[g] += MatrixX::Identity(dim,dim) * m_epsilon;
  }
}