utils.hpp

// Copyright © 2020 Thomas Nagler
//
// This file is part of the wdm library and licensed under the terms of
// the MIT license. For a copy, see the LICENSE file in the root directory
// or https://github.com/tnagler/wdm/blob/master/LICENSE.
 
#pragma once
 
#include <algorithm>
#include <cmath>
#include <numeric>
#include <stdexcept>
#include <string>
#include <vector>
 
namespace wdm {
 
namespace utils {
 
inline double
normalCDF(double x)
{
  return std::erfc(-x / std::sqrt(2)) / 2;
}
 
inline double
linear_interp(const double& x,
              const std::vector<double>& grid,
              const std::vector<double>& values)
{
  // find upper end point of interval
  size_t i = 1;
  while (x > grid[i])
    i++;
 
  // linear interpolation
  double w = (x - grid[i - 1]) / (grid[i] - grid[i - 1]);
  return w * values[i - 1] + (1 - w) * values[i];
}
 
inline void
check_sizes(const std::vector<double>& x,
            const std::vector<double>& y,
            const std::vector<double>& weights)
{
  if (y.size() != x.size())
    throw std::runtime_error("x and y must have the same size.");
  if ((weights.size() > 0) && (weights.size() != y.size()))
    throw std::runtime_error("x, y, and weights must have the same size.");
}
 
inline std::vector<double>
pow(const std::vector<double>& x, size_t n)
{
  std::vector<double> res(x.size(), 1.0);
  if (n > 0) {
    for (size_t i = 0; i < x.size(); i++) {
      for (size_t j = 0; j < n; j++) {
        res[i] *= x[i];
      }
    }
  }
 
  return res;
}
 
inline double
sum(const std::vector<double>& x)
{
  double res = 0.0;
  for (size_t i = 0; i < x.size(); i++)
    res += x[i];
  return res;
}
 
inline double
perm_sum(const std::vector<double>& x, size_t k)
{
  if (k == 0)
    return 1.0;
  double s = 0;
  for (size_t i = 1; i <= k; i++)
    s += std::pow(-1.0, i - 1) * perm_sum(x, k - i) * sum(pow(x, i));
  return s / k;
}
 
inline double
effective_sample_size(size_t n, const std::vector<double>& weights)
{
  double n_eff;
  if (weights.size() == 0) {
    n_eff = static_cast<double>(n);
  } else {
    n_eff = std::pow(sum(weights), 2);
    n_eff /= sum(pow(weights, 2));
  }
 
  return n_eff;
}
 
inline std::vector<size_t>
invert_permutation(const std::vector<size_t>& perm)
{
  std::vector<size_t> inv_perm(perm.size());
  for (size_t i = 0; i < perm.size(); i++)
    inv_perm[perm[i]] = i;
  return inv_perm;
}
 
inline std::vector<size_t>
get_order(const std::vector<double>& x, bool ascending = true)
{
  size_t n = x.size();
  std::vector<size_t> perm(n);
  for (size_t i = 0; i < n; i++)
    perm[i] = i;
  auto sorter = [&](size_t i, size_t j) {
    if (ascending)
      return (x[i] < x[j]);
    else
      return (x[i] > x[j]);
  };
  std::sort(perm.begin(), perm.end(), sorter);
 
  return perm;
}
 
inline void
sort_all(std::vector<double>& x,
         std::vector<double>& y,
         std::vector<double>& weights)
{
  size_t n = x.size();
  std::vector<size_t> order(n);
  for (size_t i = 0; i < n; i++)
    order[i] = i;
  auto sorter_with_tie_break = [&](size_t i, size_t j) {
    return (x[i] < x[j]) || ((x[i] == x[j]) && (y[i] < y[j]));
  };
  std::sort(order.begin(), order.end(), sorter_with_tie_break);
 
  std::vector<double> xx(n), yy(n);
  for (size_t i = 0; i < n; i++) {
    xx[i] = x[order[i]];
    yy[i] = y[order[i]];
  }
 
  // sort weights accordingly
  std::vector<double> w = weights;
  if (weights.size() > 0) {
    for (size_t i = 0; i < n; i++) {
      w[i] = weights[order[i]];
    }
  }
 
  x = xx;
  y = yy;
  weights = w;
}
 
inline double
count_ties_v(const std::vector<double>& x, const std::vector<double>& weights)
{
  bool weighted = (weights.size() > 0);
  double count = 0.0, w1 = 0.0, w2 = 0.0;
  size_t reps = 1;
  for (size_t i = 1; i < x.size(); i++) {
    if ((x[i] == x[i - 1])) {
      if (weighted) {
        if (reps == 1) {
          w1 = weights[i - 1];
          w2 = w1 * w1;
        }
        w1 += weights[i];
        w2 += std::pow(weights[i], 2);
      }
      reps++;
    } else if (reps > 1) {
      if (weighted) {
        count += (w1 * w1 - w2) * (2 * w1 + 5);
      } else {
        count += reps * (reps - 1) * (2 * reps + 5);
      }
      reps = 1;
    }
  }
 
  if (reps > 1) {
    if (weighted) {
      count += (w1 * w1 - w2) * (2 * w1 + 5);
    } else {
      count += reps * (reps - 1) * (2 * reps + 5);
    }
  }
 
  return count;
}
 
inline double
count_tied_pairs(const std::vector<double>& x,
                 const std::vector<double>& weights)
{
  bool weighted = (weights.size() > 0);
  double count = 0.0, w1 = 0.0, w2 = 0.0;
  size_t reps = 1;
  for (size_t i = 1; i < x.size(); i++) {
    if ((x[i] == x[i - 1])) {
      if (weighted) {
        if (reps == 1) {
          w1 = weights[i - 1];
          w2 = w1 * w1;
        }
        w1 += weights[i];
        w2 += std::pow(weights[i], 2);
      }
      reps++;
    } else if (reps > 1) {
      if (weighted) {
        count += (w1 * w1 - w2) / 2.0;
      } else {
        count += reps * (reps - 1) / 2.0;
      }
      reps = 1;
    }
  }
 
  if (reps > 1) {
    if (weighted) {
      count += (w1 * w1 - w2) / 2.0;
    } else {
      count += reps * (reps - 1) / 2.0;
    }
  }
 
  return count;
}
 
inline double
count_tied_triplets(const std::vector<double>& x,
                    const std::vector<double>& weights)
{
  bool weighted = (weights.size() > 0);
  double count = 0.0, w1 = 0.0, w2 = 0.0, w3 = 0.0;
  size_t reps = 2;
  for (size_t i = 2; i < x.size(); i++) {
    if ((x[i] == x[i - 1]) && (x[i] == x[i - 2])) {
      if (weighted) {
        if (reps == 1) {
          w1 = weights[i - 1];
          w2 = std::pow(weights[i - 1], 2);
          w3 = std::pow(weights[i - 1], 3);
        }
        w1 += weights[i];
        w2 += std::pow(weights[i], 2);
        w3 += std::pow(weights[i], 3);
      }
      reps++;
    } else if (reps > 2) {
      if (weighted) {
        count += (std::pow(w1, 3) - 3 * w2 * w1 + 2 * w3) / 6.0;
      } else {
        count += reps * (reps - 1) * (reps - 2) / 6.0;
      }
      reps = 1;
    }
  }
 
  if (reps > 2) {
    if (weighted) {
      count += (std::pow(w1, 3) - 3 * w2 * w1 + 2 * w3) / 6.0;
    } else {
      count += reps * (reps - 1) * (reps - 2) / 6.0;
    }
  }
 
  return count;
}
 
inline double
count_joint_ties(const std::vector<double>& x,
                 const std::vector<double>& y,
                 const std::vector<double>& weights)
{
  bool weighted = (weights.size() > 0);
  double count = 0.0, w1 = 0.0, w2 = 0.0;
  size_t reps = 1;
  for (size_t i = 1; i < x.size(); i++) {
    if ((x[i] == x[i - 1]) && (y[i] == y[i - 1])) {
      if (weighted) {
        if (reps == 1) {
          w1 = weights[i - 1];
          w2 = weights[i - 1] * weights[i - 1];
        }
        w1 += weights[i];
        w2 += weights[i] * weights[i];
      }
      reps++;
    } else if (reps > 1) {
      if (weighted) {
        count += (w1 * w1 - w2) / 2.0;
      } else {
        count += reps * (reps - 1) / 2.0;
      }
      reps = 1;
    }
  }
 
  if (reps > 1) {
    if (weighted) {
      count += (w1 * w1 - w2) / 2.0;
    } else {
      count += reps * (reps - 1) / 2.0;
    }
  }
 
  return count;
}
 
inline void
merge(std::vector<double>& vec,
      const std::vector<double>& vec1,
      const std::vector<double>& vec2,
      std::vector<double>& weights,
      const std::vector<double>& weights1,
      const std::vector<double>& weights2,
      double& count)
{
  double w_acc = 0.0, w1_sum = 0.0;
  bool weighted = (weights.size() > 0);
  if (weighted) {
    for (size_t i = 0; i < weights1.size(); i++)
      w1_sum += weights1[i];
  }
  size_t i, j, k;
  for (i = 0, j = 0, k = 0; i < vec1.size() && j < vec2.size(); k++) {
    if (vec1[i] <= vec2[j]) {
      vec[k] = vec1[i];
      if (weighted) {
        weights[k] = weights1[i];
        w_acc += weights1[i];
      }
      i++;
    } else {
      vec[k] = vec2[j];
      if (weighted) {
        weights[k] = weights2[j];
        count += weights2[j] * (w1_sum - w_acc);
      } else {
        count += vec1.size() - i;
      }
      j++;
    }
  }
 
  while (i < vec1.size()) {
    vec[k] = vec1[i];
    if (weighted)
      weights[k] = weights1[i];
    k++;
    i++;
  }
 
  while (j < vec2.size()) {
    vec[k] = vec2[j];
    if (weighted)
      weights[k] = weights2[j];
    k++;
    j++;
  }
}
 
inline void
merge_sort(std::vector<double>& vec,
           std::vector<double>& weights,
           double& count)
{
  if (vec.size() > 1) {
    size_t n = vec.size();
    std::vector<double> vec1(vec.begin(), vec.begin() + n / 2);
    std::vector<double> vec2(vec.begin() + n / 2, vec.end());
 
    n = weights.size();
    std::vector<double> weights1(weights.begin(), weights.begin() + n / 2);
    std::vector<double> weights2(weights.begin() + n / 2, weights.end());
 
    merge_sort(vec1, weights1, count);
    merge_sort(vec2, weights2, count);
    merge(vec, vec1, vec2, weights, weights1, weights2, count);
  }
}
 
inline void
merge_count_per_element(std::vector<double>& vec,
                        const std::vector<double>& vec1,
                        const std::vector<double>& vec2,
                        std::vector<double>& weights,
                        const std::vector<double>& weights1,
                        const std::vector<double>& weights2,
                        std::vector<double>& counts,
                        const std::vector<double>& counts1,
                        const std::vector<double>& counts2)
{
  double w_acc = 0.0;
  bool weighted = (weights.size() > 0);
  double w1_sum = 0.0;
  if (weighted) {
    for (size_t i = 0; i < weights1.size(); i++)
      w1_sum += weights1[i];
  }
  size_t i, j, k;
  for (i = 0, j = 0, k = 0; i < vec1.size() && j < vec2.size(); k++) {
    if (vec1[i] > vec2[j]) {
      vec[k] = vec1[i];
      counts[k] = counts1[i];
      if (weighted) {
        weights[k] = weights1[i];
        w_acc += weights1[i];
      }
      i++;
    } else {
      vec[k] = vec2[j];
      if (weighted) {
        counts[k] = counts2[j] + w1_sum - w_acc;
        weights[k] = weights2[j];
      } else {
        counts[k] = counts2[j] + vec1.size() - i;
      }
      j++;
    }
  }
 
  while (i < vec1.size()) {
    vec[k] = vec1[i];
    if (weighted)
      weights[k] = weights1[i];
    counts[k] = counts1[i];
    k++;
    i++;
  }
 
  while (j < vec2.size()) {
    vec[k] = vec2[j];
    if (weighted)
      weights[k] = weights2[j];
    counts[k] = counts2[j];
    k++;
    j++;
  }
}
 
inline void
merge_sort_count_per_element(std::vector<double>& vec,
                             std::vector<double>& weights,
                             std::vector<double>& counts)
{
  if (vec.size() > 1) {
    size_t n = vec.size();
    std::vector<double> vec1(vec.begin(), vec.begin() + n / 2);
    std::vector<double> vec2(vec.begin() + n / 2, vec.end());
 
    n = weights.size();
    std::vector<double> weights1(weights.begin(), weights.begin() + n / 2);
    std::vector<double> weights2(weights.begin() + n / 2, weights.end());
 
    n = counts.size();
    std::vector<double> counts1(counts.begin(), counts.begin() + n / 2);
    std::vector<double> counts2(counts.begin() + n / 2, counts.end());
 
    merge_sort_count_per_element(vec1, weights1, counts1);
    merge_sort_count_per_element(vec2, weights2, counts2);
    merge_count_per_element(
      vec, vec1, vec2, weights, weights1, weights2, counts, counts1, counts2);
  }
}
 
} 
 
} // end wdm