C++ mat::set_size方法代码示例

///
/// \brief Vespucci::Math::Quantification::FindBandwidthMat
/// \param X
/// \param abscissa
/// \param min
/// \param max
/// \param midlines
/// \param baselines
/// \return
/// Finds the bandwidth of every column of a arma::matrix.
arma::vec Vespucci::Math::Quantification::FindBandwidthMat(const arma::mat &X, arma::vec abscissa, double &min, double &max, arma::mat &midlines, arma::mat &baselines, arma::uvec &boundaries)
{
    double delta = std::abs(abscissa(1) - abscissa(0));
    arma::uvec left_bound = find(((min-delta) <= abscissa) && (abscissa <= (min+delta)));
    arma::uvec right_bound = find(((max-delta) <= abscissa) && (abscissa <= (max+delta)));

    arma::uword min_index = left_bound(0);
    arma::uword max_index = right_bound(0);
    boundaries << min_index << arma::endr << max_index;

    min = abscissa(min_index);
    max = abscissa(max_index);

    arma::uword size = abscissa.subvec(min_index, max_index).n_elem;
    arma::vec results(X.n_cols);
    midlines.set_size(X.n_cols, size);
    baselines.set_size(X.n_cols, size);
    arma::vec midline;
    arma::vec baseline;
    for (arma::uword i = 0; i < X.n_cols; ++i){
        results(i) = FindBandwidth(X.col(i), min_index, max_index, midline, baseline, delta);
        midlines.row(i) = midline;
        baselines.row(i) = baseline;
    }

    return results;
}

void FastMKS<KernelType, TreeType>::Search(TreeType* queryTree,
                                           const size_t k,
                                           arma::Mat<size_t>& indices,
                                           arma::mat& kernels)
{
  // If either naive mode or single mode is specified, this must fail.
  if (naive || singleMode)
  {
    throw std::invalid_argument("can't call Search() with a query tree when "
        "single mode or naive search is enabled");
  }

  // No remapping will be necessary because we are using the cover tree.
  indices.set_size(k, queryTree->Dataset().n_cols);
  kernels.set_size(k, queryTree->Dataset().n_cols);
  kernels.fill(-DBL_MAX);

  Timer::Start("computing_products");
  typedef FastMKSRules<KernelType, TreeType> RuleType;
  RuleType rules(referenceSet, queryTree->Dataset(), indices, kernels,
      metric.Kernel());

  typename TreeType::template DualTreeTraverser<RuleType> traverser(rules);

  traverser.Traverse(*queryTree, *referenceTree);

  Log::Info << rules.BaseCases() << " base cases." << std::endl;
  Log::Info << rules.Scores() << " scores." << std::endl;

  Timer::Stop("computing_products");
}

/**
 * Create a Lovasz-Theta initial point.
 */
void CreateLovaszThetaInitialPoint(const arma::mat& edges,
                                   arma::mat& coordinates)
{
  // Get the number of vertices in the problem.
  const size_t vertices = max(max(edges)) + 1;

  const size_t m = edges.n_cols + 1;
  float r = 0.5 + sqrt(0.25 + 2 * m);
  if (ceil(r) > vertices)
    r = vertices; // An upper bound on the dimension.

  coordinates.set_size(vertices, ceil(r));

  // Now we set the entries of the initial matrix according to the formula given
  // in Section 4 of Monteiro and Burer.
  for (size_t i = 0; i < vertices; ++i)
  {
    for (size_t j = 0; j < ceil(r); ++j)
    {
      if (i == j)
        coordinates(i, j) = sqrt(1.0 / r) + sqrt(1.0 / (vertices * m));
      else
        coordinates(i, j) = sqrt(1.0 / (vertices * m));
    }
  }
}

arma::mat Vespucci::Math::Transform::cwtPeakAnalysis(const arma::mat &X,
                              std::string wavelet, arma::uword qscale,
                              double threshold, std::string threshold_method,
                              arma::mat &transform)
{
    transform.set_size(X.n_rows, X.n_cols);
    arma::uvec scales(1);
    scales(0) = qscale;
    arma::uword i;
    arma::umat peak_positions;
    arma::mat peak_extrema(X.n_rows, X.n_cols);
    arma::vec spectrum, current_transform, dtransform, peak_magnitudes;
    try{
        for (i = 0; i < X.n_cols; ++i){
            spectrum = X.col(i);
            current_transform = Vespucci::Math::Transform::cwt(spectrum, wavelet, scales);
            transform.col(i) = current_transform;
            dtransform = Vespucci::Math::diff(transform, 1);
            dtransform.insert_rows(0, 1, true); //dtransform(i) is derivative at transform(i)
            peak_positions = Vespucci::Math::PeakFinding::FindPeakPositions(transform, dtransform,
                                                         threshold, threshold_method,
                                                         peak_magnitudes);
            peak_extrema.col(i) = Vespucci::Math::PeakFinding::PeakExtrema(X.n_elem, peak_positions);
        }
    }
    catch(std::exception e){
        std::cerr << "CWTPeakAnalysis" << std::endl;
        std::cerr << "i = " << i << std::endl;
        throw(e);
    }
    return peak_extrema;
}

arma::mat Vespucci::Math::Transform::cwt_spdbc_mat(const arma::mat &X, std::string wavelet, arma::uword qscale,
                            double threshold, std::string threshold_method,
                            arma::uword window_size, arma::field<arma::umat> &peak_positions,
                            arma::mat &baselines)
{
    baselines.set_size(X.n_rows, X.n_cols);
    arma::vec baseline;
    arma::vec spectrum;
    arma::vec current_corrected;
    arma::mat corrected(X.n_rows, X.n_cols);
    arma::umat current_peakpos;
    peak_positions.set_size(X.n_cols);
    arma::uword i;
    try{
        for (i = 0; i < X.n_cols; ++i){
            spectrum = X.col(i);
            current_corrected = Vespucci::Math::Transform::cwt_spdbc(spectrum, wavelet,
                                                   qscale, threshold,
                                                   threshold_method, window_size,
                                                   current_peakpos, baseline);

            peak_positions(i) = current_peakpos;
            baselines.col(i) = baseline;
            corrected.col(i) = current_corrected;
        }
    }catch(std::exception e){
        std::cerr << std::endl << "exception! cwt_spdbc_mat" << std::endl;
        std::cerr << "i = " << i << std::endl;
        std::cerr << e.what();
        throw(e);
    }

    return corrected;
}

/**
 * Calculates and stores the gradient values given a set of parameters.
 */
void SoftmaxRegressionFunction::Gradient(const arma::mat& parameters,
                                         arma::mat& gradient) const
{
  // Calculate the class probabilities for each training example. The
  // probabilities for each of the classes are given by:
  // p_j = exp(theta_j' * x_i) / sum(exp(theta_k' * x_i))
  // The sum is calculated over all the classes.
  // x_i is the input vector for a particular training example.
  // theta_j is the parameter vector associated with a particular class.
  arma::mat probabilities;
  GetProbabilitiesMatrix(parameters, probabilities);

  // Calculate the parameter gradients.
  gradient.set_size(parameters.n_rows, parameters.n_cols);
  if (fitIntercept)
  {
    // Treating the intercept term parameters.col(0) seperately to avoid
    // the cost of building matrix [1; data].
    arma::mat inner = probabilities - groundTruth;
    gradient.col(0) =
      inner * arma::ones<arma::mat>(data.n_cols, 1) / data.n_cols +
      lambda * parameters.col(0);
    gradient.cols(1, parameters.n_cols - 1) =
      inner * data.t() / data.n_cols +
      lambda * parameters.cols(1, parameters.n_cols - 1);
  }
  else
  {
    gradient = (probabilities - groundTruth) * data.t() / data.n_cols +
               lambda * parameters;
  }
}

void AugLagrangianTestFunction::Gradient(const arma::mat& coordinates,
                                         arma::mat& gradient)
{
  // f'_x1(x) = 12 x_1 + 4 x_2
  // f'_x2(x) = 4 x_1 + 6 x_2
  gradient.set_size(2, 1);

  gradient[0] = 12 * coordinates[0] + 4 * coordinates[1];
  gradient[1] = 4 * coordinates[0] + 6 * coordinates[1];
}

void LogisticRegression<MatType>::Classify(const MatType& dataset,
                                           arma::mat& probabilities) const
{
  // Set correct size of output matrix.
  probabilities.set_size(2, dataset.n_cols);

  probabilities.row(1) = 1.0 / (1.0 + arma::exp(-parameters(0) - dataset.t() *
      parameters.subvec(1, parameters.n_elem - 1))).t();
  probabilities.row(0) = 1.0 - probabilities.row(1);
}

 inline static void Cluster(const MatType& data,
                            const size_t clusters,
                            arma::mat& centroids)
 {
   centroids.set_size(data.n_rows, clusters);
   for (size_t i = 0; i < clusters; ++i)
   {
     // Randomly sample a point.
     const size_t index = math::RandInt(0, data.n_cols);
     centroids.col(i) = data.col(index);
   }
 }

/**
 * Remove a certain set of rows in a matrix while copying to a second matrix.
 *
 * @param input Input matrix to copy.
 * @param rowsToRemove Vector containing indices of rows to be removed.
 * @param output Matrix to copy non-removed rows into.
 */
void mlpack::math::RemoveRows(const arma::mat& input,
                              const std::vector<size_t>& rowsToRemove,
                              arma::mat& output)
{
  const size_t nRemove = rowsToRemove.size();
  const size_t nKeep = input.n_rows - nRemove;

  if (nRemove == 0)
  {
    output = input; // Copy everything.
  }
  else
  {
    output.set_size(nKeep, input.n_cols);

    size_t curRow = 0;
    size_t removeInd = 0;
    // First, check 0 to first row to remove.
    if (rowsToRemove[0] > 0)
    {
      // Note that this implies that n_rows > 1.
      output.rows(0, rowsToRemove[0] - 1) = input.rows(0, rowsToRemove[0] - 1);
      curRow += rowsToRemove[0];
    }

    // Now, check i'th row to remove to (i + 1)'th row to remove, until i is the
    // penultimate row.
    while (removeInd < nRemove - 1)
    {
      const size_t height = rowsToRemove[removeInd + 1] -
          rowsToRemove[removeInd] - 1;

      if (height > 0)
      {
        output.rows(curRow, curRow + height - 1) =
            input.rows(rowsToRemove[removeInd] + 1,
                       rowsToRemove[removeInd + 1] - 1);
        curRow += height;
      }

      removeInd++;
    }

    // Now that i is the last row to remove, check last row to remove to last
    // row.
    if (rowsToRemove[removeInd] < input.n_rows - 1)
    {
      output.rows(curRow, nKeep - 1) = input.rows(rowsToRemove[removeInd] + 1,
          input.n_rows - 1);
    }
  }
}

/**
 * Calculate the gradient.
 */
void RosenbrockFunction::Gradient(const arma::mat& coordinates,
                                  arma::mat& gradient)
{
  // f'_{x1}(x) = -2 (1 - x1) + 400 (x1^3 - (x2 x1))
  // f'_{x2}(x) = 200 (x2 - x1^2)

  double x1 = coordinates[0];
  double x2 = coordinates[1];

  gradient.set_size(2, 1);
  gradient[0] = -2 * (1 - x1) + 400 * (std::pow(x1, 3) - x2 * x1);
  gradient[1] = 200 * (x2 - std::pow(x1, 2));
}

void LogisticRegressionFunction<MatType>::Gradient(
    const arma::mat& parameters,
    arma::mat& gradient) const
{
  // Regularization term.
  arma::mat regularization;
  regularization = lambda * parameters.tail_cols(parameters.n_elem - 1);

  const arma::rowvec sigmoids = (1 / (1 + arma::exp(-parameters(0, 0)
      - parameters.tail_cols(parameters.n_elem - 1) * predictors)));

  gradient.set_size(arma::size(parameters));
  gradient[0] = -arma::accu(responses - sigmoids);
  gradient.tail_cols(parameters.n_elem - 1) = (sigmoids - responses) *
      predictors.t() + regularization;
}

void NormalizeColByMax(const arma::mat &input,
                       arma::mat &output)
{
  output.set_size(input.n_rows, input.n_cols);
  for (arma::uword i = 0; i != input.n_cols; ++i)
  {
    const double max = arma::max(arma::abs(input.col(i)));
    if (max != 0.0)
    {
      output.col(i) = input.col(i) / max;
    }
    else
    {
      output.col(i) = input.col(i);
    }
  }
}

void DualTreeBoruvka<MetricType, TreeType>::EmitResults(arma::mat& results)
{
  // Sort the edges.
  std::sort(edges.begin(), edges.end(), SortFun);

  Log::Assert(edges.size() == data.n_cols - 1);
  results.set_size(3, edges.size());

  // Need to unpermute the point labels.
  if (!naive && ownTree && tree::TreeTraits<TreeType>::RearrangesDataset)
  {
    for (size_t i = 0; i < (data.n_cols - 1); i++)
    {
      // Make sure the edge list stores the smaller index first to
      // make checking correctness easier
      size_t ind1 = oldFromNew[edges[i].Lesser()];
      size_t ind2 = oldFromNew[edges[i].Greater()];

      if (ind1 < ind2)
      {
        edges[i].Lesser() = ind1;
        edges[i].Greater() = ind2;
      }
      else
      {
        edges[i].Lesser() = ind2;
        edges[i].Greater() = ind1;
      }

      results(0, i) = edges[i].Lesser();
      results(1, i) = edges[i].Greater();
      results(2, i) = edges[i].Distance();
    }
  }
  else
  {
    for (size_t i = 0; i < edges.size(); i++)
    {
      results(0, i) = edges[i].Lesser();
      results(1, i) = edges[i].Greater();
      results(2, i) = edges[i].Distance();
    }
  }
} // EmitResults

  /**
   * Initialize the dictionary by adding together three random observations from
   * the data, and then normalizing the atom.  This implementation is simple
   * enough to be included with the definition.
   *
   * @param data Dataset to initialize the dictionary with.
   * @param atoms Number of atoms in dictionary.
   * @param dictionary Dictionary to initialize.
   */
  static void Initialize(const arma::mat& data,
                         const size_t atoms,
                         arma::mat& dictionary)
  {
    // Set the size of the dictionary.
    dictionary.set_size(data.n_rows, atoms);

    // Create each atom.
    for (size_t i = 0; i < atoms; ++i)
    {
      // Add three atoms together.
      dictionary.col(i) = (data.col(math::RandInt(data.n_cols)) +
          data.col(math::RandInt(data.n_cols)) +
          data.col(math::RandInt(data.n_cols)));

      // Now normalize the atom.
      dictionary.col(i) /= norm(dictionary.col(i), 2);
    }
  }