enpy_initest.cc 3.33 KB
Newer Older
davidkep's avatar
davidkep committed
1 2 3 4 5 6 7 8 9 10 11
//
//  enpy_initest.cc
//  pense
//
//  Created by David Kepplinger on 2019-01-30.
//  Copyright © 2019 David Kepplinger. All rights reserved.
//

#include <algorithm>

#include <nsoptim.hpp>
12
#include "constants.hpp"
davidkep's avatar
davidkep committed
13 14 15 16 17 18 19 20 21 22 23 24 25
#include "rcpp_utils.hpp"
#include "enpy_initest.hpp"

namespace {
constexpr int kDefaultCfgMaxIt = 1;  //!< Maximum number of iterations.
constexpr double kKeepPscProportion = 0.5;  //!< Proportion of observations to keep based on PSCs.
constexpr bool kUseResidualThreshold = false;  //!< Use a fixed threshold instead of a proportion
                                               //!< to screen observations based on their
                                               //!< residuals.
constexpr double kKeepResidualsProportion = 0.5;  //!< Proportion of observations to keep based
                                                  //!< on the residuals.
constexpr double kKeepResidualsThreshold = 2;  //!< Fixed threshold to keep observations based
                                               //!< on the residuals.
26 27
constexpr double kRetainBestFactor = 1.1;  //!< Retain not only the candidates from the last iteration,
                                           //!< but also those that are within this factor of the best candidate.
davidkep's avatar
davidkep committed
28 29 30 31 32 33 34 35 36

inline arma::uword HashUpdate(const arma::uword hash, const arma::uword value) noexcept;
}  // namespace

namespace pense {
namespace enpy_initest_internal {
PyConfiguration ParseConfiguration(const Rcpp::List& config) noexcept {
  return PyConfiguration{
    GetFallback(config, "max_it", kDefaultCfgMaxIt),
37
    GetFallback(config, "eps", kDefaultConvergenceTolerance),
davidkep's avatar
davidkep committed
38 39 40
    GetFallback(config, "keep_psc_proportion", kKeepPscProportion),
    GetFallback(config, "use_residual_threshold", kUseResidualThreshold),
    GetFallback(config, "keep_residuals_proportion", kKeepResidualsProportion),
41 42
    GetFallback(config, "keep_residuals_threshold", kKeepResidualsThreshold),
    GetFallback(config, "retain_best_factor", kRetainBestFactor)
davidkep's avatar
davidkep committed
43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61
  };
}

arma::uword HashIndexVector(const arma::uvec& vector) noexcept {
  arma::uword hash = vector.n_elem;
  for (const auto &val : vector) {
    hash ^= HashUpdate(hash, val);
  }
  return hash;
}

arma::uword HashSequence(const arma::uword to) noexcept {
  arma::uword hash = to + 1;
  for (arma::uword val = 0; val <= to; ++val) {
    hash ^= HashUpdate(hash, val);
  }
  return hash;
}

davidkep's avatar
davidkep committed
62
arma::uvec GetResidualKeepIndices(const arma::vec& residuals, const double mscale_est,
davidkep's avatar
davidkep committed
63 64
                                  const PyConfiguration& config, arma::uvec* all_indices) {
  if (config.use_residual_threshold) {
davidkep's avatar
davidkep committed
65
    const double resid_threshold = config.keep_residuals_threshold * mscale_est;
davidkep's avatar
davidkep committed
66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82
    return arma::find(residuals <= resid_threshold);
  } else {
    const arma::uword n_resid_keep = std::max<arma::uword>(residuals.n_elem * config.keep_residuals_proportion,
                                                           kMinObs);
    std::partial_sort(all_indices->begin(), all_indices->begin() + n_resid_keep,
                      all_indices->end(), IndexCompAbsoluteAscending<arma::vec>(residuals));
    return arma::sort(all_indices->head_rows(n_resid_keep));
  }
}
}  // namespace enpy_initest_internal
}  // namespace pense

namespace {
inline arma::uword HashUpdate(const arma::uword hash, const arma::uword value) noexcept {
  return value + 0x9e3779b9 + (hash << 6) + (hash >> 2);
}
}  // namespace