/*
 * Copyright 2018 Google LLC.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice,
 *    this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * 3. Neither the name of the copyright holder nor the names of its
 *    contributors may be used to endorse or promote products derived from this
 *    software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 *
 */

// Core mathematical routines.
//
// A quick note on terminology here.
//
// There are a bunch kinds of probabilities used commonly in genomics:
//
// -- pError: the probability of being wrong.
// -- pTrue: the probability of being correct.
//
// Normalized probabilities vs. unnormalized likelihoods:
//
// -- Normalized probabilities: p_1, ..., p_n such that sum(p_i) == 1 are said
//    said to be normalized because they represent a valid probability
//    distribution over the states 1 ... n.
// -- Unnormalized likelihoods: p_1, ..., p_n where sum(p_i) != 1. These are not
//    normalized and so aren't a valid probabilities distribution.
//
// To add even more complexity, probabilities are often represented in three
// semi-equivalent spaces:
//
// -- Real-space: the classic space, with values ranging from [0.0, 1.0]
//    inclusive.
// -- log10-space: If p is the real-space value, in log10-space this would be
//    represented as log10(p). How the p == 0 case is handled is often function
//    dependent, which may accept/return -Inf or not handle the case entirely.
// -- Phred-scaled: See https://en.wikipedia.org/wiki/Phred_quality_score for
//    more information. But briefly, the Phred-scale maintains resolution in the
//    lower parts of the probability space using integer quality scores (though
//    using ints is optional, really). The phred-scale is defined as
//
//      `phred(p) = -10 * log10(p)`
//
//    where p is a real-space probability.
//
// The code in math.h dealing with probabilities is very explicit about what
// kind probability and representation is expects and returns, as unfortunately
// these are all represented commonly as doubles in C++. Though tempting to
// address this issue with classic software engineering practices like creating
// a Probability class, in practice this is extremely difficult to do as this
// code is often performance critical and the low-level mathematical operations
// used in this code (e.g., log10) don't distiguish themselves among the types
// of probabilities.
#ifndef THIRD_PARTY_NUCLEUS_UTIL_MATH_H_
#define THIRD_PARTY_NUCLEUS_UTIL_MATH_H_

#include <vector>

namespace nucleus {

// Converts Phred scale to probability scale. Phred value must be >= 0.
double PhredToPError(int phred);

// Converts Phred scale to log10 scale. Phred value must be >= 0.
double PhredToLog10PError(int phred);

// Converts Phred scale to pError probability scale.
// Note: There is no Phred Scale equivalent for PError = 0 (would be
// infinity), so this function does not accept PError == 0.
double PErrorToPhred(double perror);
int PErrorToRoundedPhred(double perror);

// Converts probability space to Log10 space.
// Note: There is no Phred Scale equivalent for probability = 0 (would be
// infinity), so this function does not accept probability == 0.
double PErrorToLog10PError(double perror);

// Converts Log10 scale to Phred scale.
double Log10PErrorToPhred(double log10_perror);
int Log10PErrorToRoundedPhred(double log10_perror);

// Converts a Log10(ptrue) value into a phred-scaled value of 1 - 10^log10p.
//
// This operation is common when you've got a probability of an event occurring,
// p, and you want to emit the Phred-equivalent of it being wrong, which is
// -10 * log10(1 - p). The operation 1 - p can easily underflow, causing the us
// to evaluate log10(0), leading to an infinite value. In that case, the
// function returns value_if_not_finite.
double Log10PTrueToPhred(double log10_ptrue, double value_if_not_finite);

// Converts Log10 scale to real scale.
double Log10ToReal(double log10_probability);

// Takes the maximum value (remember, likelihoods are in log10 space and are all
// negative values) and subtract it from all genotype likelihoods so that the
// most likely likelihood is 0. This gives a bit more resolution in the
// conversion.
std::vector<double> ZeroShiftLikelihoods(
    const std::vector<double>& likelihoods);

}  // namespace nucleus

#endif  // THIRD_PARTY_NUCLEUS_UTIL_MATH_H_