normalize.h

Go to the documentation of this file.

 
#pragma once
 
#include "jit_normalization.h"
#include "math.h"
#include <Eigen/SparseCore>
 
namespace slope {
 
template<typename T>
void
computeCenters(Eigen::VectorXd& x_centers, const T& x, const std::string& type)
{
  int p = x.cols();
 
  if (type == "manual") {
    if (x_centers.size() != p) {
      throw std::invalid_argument("Invalid dimensions in centers");
    }
 
    if (!x_centers.allFinite()) {
      throw std::invalid_argument("Centers must be finite");
    }
 
  } else if (type == "mean") {
    x_centers = means(x);
  } else if (type == "min") {
    x_centers = mins(x);
  } else if (type != "none") {
    throw std::invalid_argument("Invalid centering type");
  }
}
 
template<typename T>
void
computeScales(Eigen::VectorXd& x_scales, const T& x, const std::string& type)
{
  int p = x.cols();
 
  if (type == "manual") {
    if (x_scales.size() != p) {
      throw std::invalid_argument("Invalid dimensions in scales");
    }
    if (!x_scales.allFinite()) {
      throw std::invalid_argument("Scales must be finite");
    }
  } else if (type == "sd") {
    x_scales = stdDevs(x);
  } else if (type == "l1") {
    x_scales = l1Norms(x);
  } else if (type == "l2") {
    x_scales = l2Norms(x);
  } else if (type == "max_abs") {
    x_scales = maxAbs(x);
  } else if (type == "range") {
    x_scales = ranges(x);
  } else if (type != "none") {
    throw std::invalid_argument("Invalid scaling type");
  }
}
 
template<typename T>
JitNormalization
normalize(Eigen::MatrixBase<T>& x,
          Eigen::VectorXd& x_centers,
          Eigen::VectorXd& x_scales,
          const std::string& centering_type,
          const std::string& scaling_type,
          const bool modify_x)
{
  const int p = x.cols();
 
  computeCenters(x_centers, x, centering_type);
  computeScales(x_scales, x, scaling_type);
 
  // Handle zero-variance columns by setting their scale to 1.0
  // to avoid division by zero during normalization
  if (scaling_type != "none" && scaling_type != "manual") {
    for (int j = 0; j < p; ++j) {
      if (std::abs(x_scales(j)) == 0.0) {
        x_scales(j) = 1.0;
      }
    }
  }
 
  bool center = centering_type != "none";
  bool scale = scaling_type != "none";
  bool center_jit = center && !modify_x;
  bool scale_jit = scale && !modify_x;
 
  JitNormalization jit_normalization;
 
  if (center_jit && scale_jit) {
    jit_normalization = JitNormalization::Both;
  } else if (center_jit) {
    jit_normalization = JitNormalization::Center;
  } else if (scale_jit) {
    jit_normalization = JitNormalization::Scale;
  } else {
    jit_normalization = JitNormalization::None;
  }
 
  if (modify_x && (center || scale)) {
    for (int j = 0; j < p; ++j) {
      if (center) {
        x.col(j).array() -= x_centers(j);
      }
      if (scale) {
        x.col(j).array() /= x_scales(j);
      }
    }
  }
 
  return jit_normalization;
}
 
template<typename T>
JitNormalization
normalize(Eigen::SparseMatrixBase<T>& x,
          Eigen::VectorXd& x_centers,
          Eigen::VectorXd& x_scales,
          const std::string& centering_type,
          const std::string& scaling_type,
          const bool)
{
  computeCenters(x_centers, x, centering_type);
  computeScales(x_scales, x, scaling_type);
 
  // Handle zero-variance columns by setting their scale to 1.0
  // to avoid division by zero during normalization
  const int p = x.cols();
  if (scaling_type != "none" && scaling_type != "manual") {
    for (int j = 0; j < p; ++j) {
      if (std::abs(x_scales(j)) == 0.0) {
        x_scales(j) = 1.0;
      }
    }
  }
 
  bool center = centering_type != "none";
  bool scale = scaling_type != "none";
  bool center_jit = center;
  bool scale_jit = scale;
 
  JitNormalization jit_normalization;
 
  if (center_jit && scale_jit) {
    jit_normalization = JitNormalization::Both;
  } else if (center_jit) {
    jit_normalization = JitNormalization::Center;
  } else if (scale_jit) {
    jit_normalization = JitNormalization::Scale;
  } else {
    jit_normalization = JitNormalization::None;
  }
 
  // TODO: Implement in-place scaling for sparse matrices.
  // if (modify_x && scaling_type != "none") {
  //   for (int j = 0; j < x.cols(); ++j) {
  //     for (Eigen::SparseMatrix<double>::InnerIterator it(x, j); it; ++it) {
  //       it.valueRef() = it.value() / x_scales(j);
  //     }
  //   }
  // }
 
  return jit_normalization;
}
 
std::tuple<Eigen::VectorXd, Eigen::MatrixXd>
rescaleCoefficients(const Eigen::VectorXd& beta0,
                    const Eigen::SparseMatrix<double>& beta,
                    const Eigen::VectorXd& x_centers,
                    const Eigen::VectorXd& x_scales,
                    const bool intercept);
 
} // namespace slope