mad.hpp

Go to the documentation of this file.

#ifndef QUICKSTATS_MAD_HPP
#define QUICKSTATS_MAD_HPP
 
#include <cstddef>
#include <cmath>
#include <limits>
 
#include "median.hpp"
 
namespace quickstats {
 
template<typename Output_ = double>
struct MadOptions {
    Output_ placeholder = nan_if_available_else_zero<Output_>();
 
    bool difference_between_infinities_is_zero = false;
};
 
template<typename Output_ = double, typename Input_>
Output_ mad(const std::size_t num_total, Input_* const ptr, const Input_ median, const MadOptions<Output_>& options) {
    // Check for infinities to avoid shenanigans from trying to order NaN deviations from Inf - Inf.
    // This entire conditional should be optimized out if Input_ has no infinities.
    if (std::isinf(median)) {
        if (!options.difference_between_infinities_is_zero) {
            return options.placeholder;
        } 
 
        // We use inf_if_available_else_max() just to get it to compile if Input_ doesn't have infs.
        // At this point, Input_ must support infinities otherwise we wouldn't have gotten here.
        for (std::size_t i = 0; i < num_total; ++i) {
            ptr[i] = (median == ptr[i] ? 0 : inf_if_available_else_max<Input_>());
        }
 
    } else {
        for (std::size_t i = 0; i < num_total; ++i) {
            ptr[i] = std::abs(ptr[i] - median);
        }
    }
 
    MedianOptions<Output_> medopt;
    medopt.placeholder = options.placeholder;
    return ::quickstats::median<Output_>(num_total, ptr, medopt);
}
 
template<typename Output_ = double, typename Input_>
Output_ mad(const std::size_t num_total, const std::size_t num_non_zero, Input_* const values, const Input_ median, const MadOptions<Output_>& options) {
    // Check for infinities to avoid shenanigans from trying to order NaN deviations from Inf - Inf.
    // This entire conditional should be optimized out if Input_ has no infinities.
    if (std::isinf(median)) {
        if (!options.difference_between_infinities_is_zero) {
            return options.placeholder;
        }
 
        // We use inf_if_available_else_max() just to get it to compile if Input_ doesn't have infs.
        // At this point, Input_ must support infinities otherwise we wouldn't have gotten here.
        for (std::size_t i = 0; i < num_non_zero; ++i) {
            values[i] = (median == values[i] ? 0 : inf_if_available_else_max<Input_>());
        }
 
    } else {
        for (std::size_t i = 0; i < num_non_zero; ++i) {
            values[i] = std::abs(values[i] - median);
        }
    }
 
    // It is also possible to implement the sparse MAD by subtracting 'abs(median)' from the absolute deviations,
    // computing the sparse median of the difference, and then adding 'abs(median)' back to the result.
    // We don't do this as the subtraction and addition introduces some numerical error,
    // which isn't that consequential in practice but interferes with exact comparisons to the dense results in our tests.
 
    MedianOptions<Output_> medopt;
    medopt.placeholder = options.placeholder;
    return median_internal<Output_>(num_total, num_non_zero, values, std::abs(median), medopt);
}
 
template<typename Float_>
Float_ scale_mad_to_sd(const Float_ x) {
    return x * 1.4826;
}
 
// Backwards compatibility.
template<typename Output_ = double, typename Input_>
Output_ mad(const std::size_t num_total, Input_* const ptr, const Input_ median) {
    return mad(num_total, ptr, median, MadOptions<Output_>());
}
 
template<typename Output_ = double, typename Input_>
Output_ mad(const std::size_t num_total, const std::size_t num_non_zero, Input_* const values, const Input_ median) {
    return mad(num_total, num_non_zero, values, median, MadOptions<Output_>());
}
 
template<typename Output_ = double, typename Input_>
Output_ mad_with_infinities(const std::size_t num_total, Input_* const ptr, const Input_ median) {
    MadOptions<Output_> opt;
    opt.difference_between_infinities_is_zero = true;
    return mad(num_total, ptr, median, opt);
}
 
template<typename Output_ = double, typename Input_>
Output_ mad_with_infinities(const std::size_t num_total, const std::size_t num_non_zero, Input_* const values, const Input_ median) {
    MadOptions<Output_> opt;
    opt.difference_between_infinities_is_zero = true;
    return mad(num_total, num_non_zero, values, median, opt);
}
}
 
#endif