C++ Mat::col方法代码示例

void Split(const arma::Mat<T>& input,
           const arma::Row<U>& inputLabel,
           arma::Mat<T>& trainData,
           arma::Mat<T>& testData,
           arma::Row<U>& trainLabel,
           arma::Row<U>& testLabel,
           const double testRatio)
{
  const size_t testSize = static_cast<size_t>(input.n_cols * testRatio);
  const size_t trainSize = input.n_cols - testSize;
  trainData.set_size(input.n_rows, trainSize);
  testData.set_size(input.n_rows, testSize);
  trainLabel.set_size(trainSize);
  testLabel.set_size(testSize);

  const arma::Col<size_t> order =
      arma::shuffle(arma::linspace<arma::Col<size_t>>(0, input.n_cols - 1,
                                                      input.n_cols));

  for (size_t i = 0; i != trainSize; ++i)
  {
    trainData.col(i) = input.col(order[i]);
    trainLabel(i) = inputLabel(order[i]);
  }

  for (size_t i = 0; i != testSize; ++i)
  {
    testData.col(i) = input.col(order[i + trainSize]);
    testLabel(i) = inputLabel(order[i + trainSize]);
  }
}

arma::Col<double> compute_column_rms(const arma::Mat<double>& data) {
	const long n_cols = data.n_cols;
	arma::Col<double> rms(n_cols);
    for (long i=0; i<n_cols; ++i) {
        const double dot = arma::dot(data.col(i), data.col(i));
        rms(i) = std::sqrt(dot / (data.col(i).n_rows-1));
    }
	return std::move(rms);
}

void enforce_positive_sign_by_column(arma::Mat<double>& data) {
	for (long i=0; i<long(data.n_cols); ++i) {
		const double max = arma::max(data.col(i));
		const double min = arma::min(data.col(i));
		bool change_sign = false;
		if (std::abs(max)>=std::abs(min)) {
			if (max<0) change_sign = true;
		} else {
			if (min<0) change_sign = true;
		}
		if (change_sign) data.col(i) *= -1;
	}
}

void calculateJacobian(const arma::Mat<std::complex<double> >& myOffsets,
		       arma::Mat<double>& myJacobian, 
		       arma::Col<std::complex<double> >& myTargetsCalculated, 
		       arma::Col<std::complex<double> >& myCurrentGuess, 
		       void myCalculateDependentVariables(const arma::Mat<std::complex<double> >&, const arma::Col<std::complex<double> >&, arma::Col<std::complex<double> >&))
{
	//Calculate a temporary, unperturbed target evaluation, such as is needed for solving for the updated guess 
	//formula
	arma::Col<std::complex<double> > unperturbedTargetsCalculated(NUMDIMENSIONS);
	unperturbedTargetsCalculated.fill(0.0);
	myCalculateDependentVariables(myOffsets, myCurrentGuess, unperturbedTargetsCalculated);
	std::complex<double> oldGuessValue(0.0, 0.0);

	//Each iteration fills a column in the Jacobian
	//The Jacobian takes this form:
	//
	//	dF0/dx0 dF0/dx1 
	//	dF1/dx0 dF1/dx1
	//
	for(int j = 0; j< NUMDIMENSIONS; j++)
	{
		//Store old element value, perturb the current value
		oldGuessValue = myCurrentGuess[j];
		myCurrentGuess[j] += std::complex<double>(0.0, PROBEDISTANCE);

		//Evaluate functions for perturbed guess
		myCalculateDependentVariables(myOffsets, myCurrentGuess, myTargetsCalculated);

		//The column of the Jacobian that goes with the independent variable we perturbed
		//can be determined using the finite-difference formula
		//The arma::Col allows this to be expressed as a single vector operation
		//note slice works as: std::slice(start_index, number_of_elements_to_access, index_interval_between_selections)
		//std::cout << "Jacobian column " << j << " with:" << std::endl;
		//std::cout << "myTargetsCalculated" << std::endl;
		//std::cout << myTargetsCalculated << std::endl;
		//std::cout << "unperturbedTargetsCalculated" << std::endl;
		//std::cout << unperturbedTargetsCalculated << std::endl;
		myJacobian.col(j) = arma::imag(myTargetsCalculated);
	       	myJacobian.col(j) *= pow(PROBEDISTANCE, -1.0);
		//std::cout << "The jacobian: " << std::endl;
		//std::cout << myJacobian << std::endl;

		myCurrentGuess[j] = oldGuessValue;
	}

	//Reset to unperturbed, so we dont waste a function evaluation
	myTargetsCalculated = unperturbedTargetsCalculated;
}

  void Backward(const arma::Cube<eT>& input,
                const arma::Mat<eT>& gy,
                arma::Cube<eT>& g)
  {
    // Generate a cube using the backpropagated error matrix.
    arma::Cube<eT> mappedError = arma::zeros<arma::cube>(input.n_rows,
        input.n_cols, input.n_slices);

    for (size_t s = 0, j = 0; s < mappedError.n_slices; s+= gy.n_cols, j++)
    {
      for (size_t i = 0; i < gy.n_cols; i++)
      {
        arma::Col<eT> temp = gy.col(i).subvec(
            j * input.n_rows * input.n_cols,
            (j + 1) * input.n_rows * input.n_cols - 1);

        mappedError.slice(s + i) = arma::Mat<eT>(temp.memptr(),
            input.n_rows, input.n_cols);
      }
    }

    arma::Cube<eT> derivative;
    Deriv(input, derivative);
    g = mappedError % derivative;
  }

    inline ParallelKinematicMachine3PUPS::ParallelKinematicMachine3PUPS() noexcept {
      setMinimalActiveJointActuations({0.39, 0.39, 0.39});
      setMaximalActiveJointActuations({0.95, 0.95, 0.95});

      setEndEffectorJointPositions({
        -0.025561381023353, 0.086293776138137, 0.12,
        0.087513292835791, -0.021010082747031, 0.12,
        -0.061951911812438, -0.065283693391106, 0.12});

      setRedundantJointStartPositions({
        -0.463708870031622, 0.417029254828353, -0.346410161513775,
        0.593012363818459, 0.193069033993384, -0.346410161513775,
        -0.129303493786837, -0.610098288821738, -0.346410161513775});

      setRedundantJointEndPositions({
        -0.247202519085512, 0.292029254828353, 0.086602540378444,
        0.376506012872349, 0.068069033993384, 0.086602540378444,
        -0.129303493786837, -0.360098288821738, 0.086602540378444});

      redundantJointStartToEndPositions_ = redundantJointEndPositions_ - redundantJointStartPositions_;
      redundantJointIndicies_ = arma::find(arma::any(redundantJointStartToEndPositions_));

      redundantJointAngles_.set_size(3, redundantJointIndicies_.n_elem);

      for (std::size_t n = 0; n < redundantJointIndicies_.n_elem; ++n) {
        const double& redundantJointXAngle = std::atan2(redundantJointStartToEndPositions_(1, n), redundantJointStartToEndPositions_(0, n));
        const double& redundantJointYAngle = std::atan2(redundantJointStartToEndPositions_(2, n), redundantJointStartToEndPositions_(1, n));
        redundantJointAngles_.col(n) = arma::Col<double>::fixed<3>({std::cos(redundantJointXAngle) * std::cos(redundantJointYAngle), std::sin(redundantJointXAngle) * std::cos(redundantJointYAngle), std::sin(redundantJointYAngle)});
      }
    }

  void FeedBackward(const arma::Cube<eT>& inputActivation,
                    const arma::Mat<eT>& error,
                    arma::Cube<eT>& delta)
  {
    delta = delta % mask * scale;

    // Generate a cube from the error matrix.
    arma::Cube<eT> mappedError = arma::zeros<arma::cube>(inputActivation.n_rows,
        inputActivation.n_cols, inputActivation.n_slices);

    for (size_t s = 0, j = 0; s < mappedError.n_slices; s+= error.n_cols, j++)
    {
      for (size_t i = 0; i < error.n_cols; i++)
      {
        arma::Col<eT> temp = error.col(i).subvec(
            j * inputActivation.n_rows * inputActivation.n_cols,
            (j + 1) * inputActivation.n_rows * inputActivation.n_cols - 1);

        mappedError.slice(s + i) = arma::Mat<eT>(temp.memptr(),
            inputActivation.n_rows, inputActivation.n_cols);
      }
    }

    delta = mappedError;
  }

arma::Row<T> column_apply(const arma::Mat<T>& matrix, Functor functor) {
  arma::Row<T> result(matrix.n_cols);
  std::size_t index = 0;
  std::generate(result.begin(), result.end(), [&matrix, &index, &functor] {
    return functor(matrix.col(index++));
  });
  return result;
}

void normalize_by_column(arma::Mat<double>& data, const arma::Col<double>& rms) {
	if (data.n_cols != rms.n_elem)
		throw std::range_error("Number of elements of rms is not equal to the number of columns of data");
	for (long i=0; i<long(data.n_cols); ++i) {
        if (rms(i)==0)
        	throw std::runtime_error("At least one of the entries of rms equals to zero");
        data.col(i) *= 1./rms(i);
    }
}

// Predict the rating for a group of user/item combinations.
void CF::Predict(const arma::Mat<size_t>& combinations,
                 arma::vec& predictions) const
{
  // First, for nearest neighbor search, stretch the H matrix.
  arma::mat l = arma::chol(w.t() * w);
  arma::mat stretchedH = l * h; // Due to the Armadillo API, l is L^T.

  // Now, we must determine those query indices we need to find the nearest
  // neighbors for.  This is easiest if we just sort the combinations matrix.
  arma::Mat<size_t> sortedCombinations(combinations.n_rows,
                                       combinations.n_cols);
  arma::uvec ordering = arma::sort_index(combinations.row(0).t());
  for (size_t i = 0; i < ordering.n_elem; ++i)
    sortedCombinations.col(i) = combinations.col(ordering[i]);

  // Now, we have to get the list of unique users we will be searching for.
  arma::Col<size_t> users = arma::unique(combinations.row(0).t());

  // Assemble our query matrix from the stretchedH matrix.
  arma::mat queries(stretchedH.n_rows, users.n_elem);
  for (size_t i = 0; i < queries.n_cols; ++i)
    queries.col(i) = stretchedH.col(users[i]);

  // Now calculate the neighborhood of these users.
  neighbor::KNN a(stretchedH);
  arma::mat distances;
  arma::Mat<size_t> neighborhood;

  a.Search(queries, numUsersForSimilarity, neighborhood, distances);

  // Now that we have the neighborhoods we need, calculate the predictions.
  predictions.set_size(combinations.n_cols);

  size_t user = 0; // Cumulative user count, because we are doing it in order.
  for (size_t i = 0; i < sortedCombinations.n_cols; ++i)
  {
    // Could this be made faster by calculating dot products for multiple items
    // at once?
    double rating = 0.0;

    // Map the combination's user to the user ID used for kNN.
    while (users[user] < sortedCombinations(0, i))
      ++user;

    for (size_t j = 0; j < neighborhood.n_rows; ++j)
      rating += arma::as_scalar(w.row(sortedCombinations(1, i)) *
          h.col(neighborhood(j, user)));
    rating /= neighborhood.n_rows;

    predictions(ordering[i]) = rating;
  }
}

int
DynamicsTrappenberg::update_y(
    const arma::Mat< double > x,
    const arma::Mat< double > u,
    const mxArray *theta,
    const mxArray *ptheta,
    arma::Mat< double > &y)
{

    y = x.col(INDEX_NODES);

    return 0;

}

  /**
   * Given the sufficient statistics of a proposed split, calculate the
   * information gain if that split was to be used.  The 'counts' matrix should
   * contain the number of points in each class in each column, so the size of
   * 'counts' is children x classes, where 'children' is the number of child
   * nodes in the proposed split.
   *
   * @param counts Matrix of sufficient statistics.
   */
  static double Evaluate(const arma::Mat<size_t>& counts)
  {
    // Calculate the number of elements in the unsplit node and also in each
    // proposed child.
    size_t numElem = 0;
    arma::vec splitCounts(counts.n_elem);
    for (size_t i = 0; i < counts.n_cols; ++i)
    {
      splitCounts[i] = arma::accu(counts.col(i));
      numElem += splitCounts[i];
    }

    // Corner case: if there are no elements, the gain is zero.
    if (numElem == 0)
      return 0.0;

    arma::Col<size_t> classCounts = arma::sum(counts, 1);

    // Calculate the gain of the unsplit node.
    double gain = 0.0;
    for (size_t i = 0; i < classCounts.n_elem; ++i)
    {
      const double f = ((double) classCounts[i] / (double) numElem);
      if (f > 0.0)
        gain -= f * std::log2(f);
    }

    // Now calculate the impurity of the split nodes and subtract them from the
    // overall gain.
    for (size_t i = 0; i < counts.n_cols; ++i)
    {
      if (splitCounts[i] > 0)
      {
        double splitGain = 0.0;
        for (size_t j = 0; j < counts.n_rows; ++j)
        {
          const double f = ((double) counts(j, i) / (double) splitCounts[i]);
          if (f > 0.0)
            splitGain += f * std::log2(f);
        }

        gain += ((double) splitCounts[i] / (double) numElem) * splitGain;
      }
    }

    return gain;
  }

  static double Evaluate(const arma::Mat<size_t>& counts)
  {
    // We need to sum over the difference between the un-split node and the
    // split nodes.  First we'll calculate the number of elements in each split
    // and total.
    size_t numElem = 0;
    arma::vec splitCounts(counts.n_cols);
    for (size_t i = 0; i < counts.n_cols; ++i)
    {
      splitCounts[i] = arma::accu(counts.col(i));
      numElem += splitCounts[i];
    }

    // Corner case: if there are no elements, the impurity is zero.
    if (numElem == 0)
      return 0.0;

    arma::Col<size_t> classCounts = arma::sum(counts, 1);

    // Calculate the Gini impurity of the un-split node.
    double impurity = 0.0;
    for (size_t i = 0; i < classCounts.n_elem; ++i)
    {
      const double f = ((double) classCounts[i] / (double) numElem);
      impurity += f * (1.0 - f);
    }

    // Now calculate the impurity of the split nodes and subtract them from the
    // overall impurity.
    for (size_t i = 0; i < counts.n_cols; ++i)
    {
      if (splitCounts[i] > 0)
      {
        double splitImpurity = 0.0;
        for (size_t j = 0; j < counts.n_rows; ++j)
        {
          const double f = ((double) counts(j, i) / (double) splitCounts[i]);
          splitImpurity += f * (1.0 - f);
        }

        impurity -= ((double) splitCounts[i] / (double) numElem) *
            splitImpurity;
      }
    }

    return impurity;
  }

arma::Mat< double >
DynamicsTrappenberg::sigma(
    const arma::Mat< double > x,
    const arma::Mat< double > u,
    const mxArray *theta,
    const mxArray *ptheta)
{
    arma::Mat< double > beta(mxGetPr(mxGetField(theta, 0, "beta")), 
        x.n_rows, 1, 1);
    arma::Mat< double > baseline(mxGetPr(mxGetField(theta, 0, "theta")), 
        x.n_rows, 1, 1);
   
    arma::Mat< double > dx = 1/(1 + exp(-x.col(INDEX_NODES) % beta + 
                baseline));

    return dx;

}

    /**
     * Convert the features of images into histogram
     * @param features features of the image
     * @param code_book code book of the data sets
     * @return histogram of bovw
     */
    Hist describe(arma::Mat<T> const &features,
                  arma::Mat<T> const &code_book) const
    {
        arma::Mat<T> dist(features.n_cols, code_book.n_cols);
        for(arma::uword i = 0; i != features.n_cols; ++i){
            dist.row(i) = euclidean_dist(features.col(i),
                                         code_book);
        }
        //dist.print("dist");

        Hist hist = create_hist(dist.n_cols,
                                armd::is_two_dim<Hist>::type());
        for(arma::uword i = 0; i != dist.n_rows; ++i){
            arma::uword min_idx;
            dist.row(i).min(min_idx);
            ++hist(min_idx);
        }
        //hist.print("\nhist");

        return hist;
    }