##########################################################
#                                                        #
#    Utility methods for Environmental Contours          #
#                                                        #
##########################################################

from math import *
from numba import jit


##########################################################
#    Gaussian distribution methods                       #
##########################################################
@jit(nopython=True)
def ran_gaussian_pdf(x):
    return exp(-x * x / 2.0) / 2.50662827463100029


@jit(nopython=True)
def ran_gaussian_cdf(x):
    if (x >= 0.0):
        t = 1.0/(1.0 + 0.33267 * x)
        return 1.0 - ran_gaussian_pdf(x) * (0.4361836*t - 0.1201676*t*t + 0.9372980*t*t*t)
    else:
        t = 1.0/(1.0 - 0.33267 * x)
        return ran_gaussian_pdf(x) * (0.4361836*t - 0.1201676*t*t + 0.9372980*t*t*t)


@jit(nopython=True)
def ran_gaussian_sdf(x):
    if (x >= 0.0):
        t = 1.0/(1.0 + 0.33267 * x)
        return ran_gaussian_pdf(x) * (0.4361836*t - 0.1201676*t*t + 0.9372980*t*t*t)
    else:
        t = 1.0/(1.0 - 0.33267 * x)
        return 1.0 - ran_gaussian_pdf(x) * (0.4361836*t - 0.1201676*t*t + 0.9372980*t*t*t)


@jit(nopython=True)
def ran_gaussian_invcdf(u):
        if (u < 0.5):
            t = sqrt(log(1.0/(u*u)))
            return -t + (2.515517 + 0.802853*t + 0.010328*t*t) / (1.0 + 1.432788*t + 0.189269*t*t + 0.001308*t*t*t)
        elif (u == 0.5):
            return 0.0;
        else:
            t = sqrt(log(1.0/((1.0 - u)*(1.0 - u))))
            return t - (2.515517 + 0.802853*t + 0.010328*t*t) / (1.0 + 1.432788*t + 0.189269*t*t + 0.001308*t*t*t)

      
@jit(nopython=True)
def ran_gaussian_invsdf(u):
        if (u < 0.5):
            t = sqrt(log(1.0/(u*u)))
            return t - (2.515517 + 0.802853*t + 0.010328*t*t) / (1.0 + 1.432788*t + 0.189269*t*t + 0.001308*t*t*t)
        elif (u == 0.5):
            return 0.0;
        else:
            t = sqrt(log(1.0/((1.0 - u)*(1.0 - u))))
            return -t + (2.515517 + 0.802853*t + 0.010328*t*t) / (1.0 + 1.432788*t + 0.189269*t*t + 0.001308*t*t*t)

       
##########################################################
#    Exceedence probability methods                      #
##########################################################

# Calculate the exceedence probability for a given return period (years)
# The sampling rate is the number of seastate samples pr. day (e.g., 8)
@jit(nopython=True)
def getExceedenceProbability(samplingRate, returnPeriod):
    return 1.0 / (returnPeriod * samplingRate * 365.25)

# Calculate the inceedence probability for a given return period (years)
# The sampling rate is the number of seastate samples pr. day (e.g., 8)
@jit(nopython=True)
def getNonExceedenceProbability(samplingRate, returnPeriod):
    return 1.0 - 1.0 / (returnPeriod * samplingRate * 365.25)

# The inner radius factor denotes the ratio between the
# radius of the circle which we sample outside when we
# estimate the form curve, and the radius of the circle 
# used to produce the Rosenblatt contour.
# 
# The inner radius factor should always be less than 1.0.
# Note, however, that if the radius factor is too close to 1.0,
# the boundary of the set of excluded sampling points, i.e, 
# the boundary of the transformed inner circle, may interfer 
# with the form curve, and thus yielding incorrect estimates. 
# 
# A factor of 0.95 appears to produce good results. 
# A lower value will also work, but then we need more 
# simulations in order to obtain stable results.
INNER_RADIUS_FACTOR = 0.95;
    
def setInnerRadiusFactor(irf):
    global INNER_RADIUS_FACTOR
    INNER_RADIUS_FACTOR = irf


# The arma inner radius factor denotes the ratio between the
# radius of the circle which we sample outside when we
# estimate the arma form curve, and the radius of the circle 
# used to produce the Rosenblatt contour.
# 
# The arma inner radius factor should always be less than 1.0.
# Note, however, that if the radius factor is too close to 1.0,
# the boundary of the set of excluded sampling points, i.e, 
# the boundary of the transformed inner circle, may interfer 
# with the arma form curve, and thus yielding incorrect estimates. 
# 
# A factor of 0.80 appears to produce good results. 
# A lower value will also work, but then we need more 
# simulations in order to obtain stable results.
ARMA_INNER_RADIUS_FACTOR = 0.80;
    
def setArmaInnerRadiusFactor(irf):
    global ARMA_INNER_RADIUS_FACTOR
    ARMA_INNER_RADIUS_FACTOR = irf


# Get the radius of the circle corresponding to a given returnPeriod.
@jit(nopython=True)
def getRadius(samplingRate, returnPeriod):
    qe = getNonExceedenceProbability(samplingRate, returnPeriod)
    return ran_gaussian_invcdf(qe)


# Use this method to calculate the radius of circle inside the circle
# corresponding to the given returnPeriod. This inner circle can be
# used when sampling points in order to calculate a form plot. By
# sampling points outside this inner circle only, we can reduce the
# number of samples needed to obtain a stable form estimate.
@jit(nopython=True)
def getInnerRadius(samplingRate, returnPeriod):
    qe = 1.0 - getExceedenceProbability(samplingRate, returnPeriod)
    return INNER_RADIUS_FACTOR * ran_gaussian_invcdf(qe)
    

# Use this method to calculate the radius of circle inside the circle
# corresponding to the given returnPeriod. This arma inner circle can be
# used when sampling points in order to calculate an arma form plot. By
# sampling points outside this inner circle only, we can reduce the
# number of samples needed to obtain a stable form estimate.
@jit(nopython=True)
def getArmaInnerRadius(samplingRate, returnPeriod):
    qe = 1.0 - getExceedenceProbability(samplingRate, returnPeriod)
    return ARMA_INNER_RADIUS_FACTOR * ran_gaussian_invcdf(qe)


# Use this method to calculate the probability that a random point is outside
# the circle corresponding to a given returnPeriod.
@jit(nopython=True)
def getCircleExceedenceProbability(samplingRate, returnPeriod):
    r = getRadius(samplingRate, returnPeriod)
    return exp(-0.5 * r * r)


# Use this method to calculate the probability that a random point is inside
# the circle corresponding to a given returnPeriod.    
@jit(nopython=True)
def getCircleNonExceedenceProbability(samplingRate, returnPeriod):
    r = getRadius(samplingRate, returnPeriod)
    return 1.0 - exp(-0.5 * r * r)
    

# Use this method to calculate the probability that a random point is outside
# the inner circle corresponding to a given returnPeriod.
@jit(nopython=True)
def getInnerCircleExceedenceProbability(samplingRate, returnPeriod):
    ir = getInnerRadius(samplingRate, returnPeriod)
    return exp(-0.5 * ir * ir)
    

# Use this method to calculate the probability that a random point is inside
# the inner circle corresponding to a given returnPeriod.    
@jit(nopython=True)
def getInnerCircleNonExceedenceProbability(samplingRate, returnPeriod):
    ir = getInnerRadius(samplingRate, returnPeriod)
    return 1.0 - exp(-0.5 * ir * ir)


# Use this method to calculate the probability that a random point is outside
# the arma inner circle corresponding to a given returnPeriod.    
@jit(nopython=True)
def getArmaInnerCircleExceedenceProbability(samplingRate, returnPeriod):
    ir = getArmaInnerRadius(samplingRate, returnPeriod)
    return exp(-0.5 * ir * ir)


# Use this method to calculate the probability that a random point is inside
# the arma inner circle corresponding to a given returnPeriod.    
@jit(nopython=True)
def getArmaInnerCircleNonExceedenceProbability(samplingRate, returnPeriod):
    ir = getArmaInnerRadius(samplingRate, returnPeriod)
    return 1.0 - exp(-0.5 * ir * ir)
    

# Use this method to calculate the "adjusted exceedence probability"
# P_e'. This probability is the probability that the projection of 
# a random point is exceeding the critical value C(theta) given that the 
# point is outside of the transformed inner circle. This adjusted
# probability can be used to estimate C(theta) given that all the
# sampled points are outside the inner circle.
# 
# Note that P_e = exp(-r^2/2) P_e', where the factor exp(-r^2/2)
# is the probability that a random point is outside the transformed
# inner circle.
@jit(nopython=True)
def getAdjExceedenceProbability(samplingRate, returnPeriod):
    ir = getInnerRadius(samplingRate, returnPeriod)
    pe = getExceedenceProbability(samplingRate, returnPeriod)
    return exp(0.5 * ir * ir) * pe / returnPeriod
    

# Use this method to calculate the "adjusted arma exceedence probability"
# P_e'. This probability is the probability that the projection of 
# a random point is exceeding the critical value C(theta) given that the 
# point is outside of the transformed arma inner circle. This adjusted
# probability can be used to estimate C(theta) given that all the
# sampled points are outside the arma inner circle.
# 
# Note that P_e = exp(-r^2/2) P_e', where the factor exp(-r^2/2)
# is the probability that a random point is outside the transformed
# arma inner circle.
@jit(nopython=True)
def getAdjArmaExceedenceProbability(samplingRate, returnPeriod):
    ir = getArmaInnerRadius(samplingRate, returnPeriod)
    pe = getExceedenceProbability(samplingRate, returnPeriod)
    return exp(0.5 * ir * ir) * pe / returnPeriod
    

# Use this method to calculate the "adjusted nonexceedence probability"
# Q_e'. This probability is the probability that the projection of 
# a random point is not exceeding the critical value C(theta) given that the 
# point is outside of the transformed inner circle. This adjusted
# probability can be used to estimate C(theta) given that all the
# sampled points are outside the inner circle.
@jit(nopython=True)
def getAdjNonExceedenceProbability(samplingRate, returnPeriod):
    ir = getInnerRadius(samplingRate, returnPeriod)
    pe = getExceedenceProbability(samplingRate, returnPeriod)
    return 1.0 - exp(0.5 * ir * ir) * pe
    

# Use this method to calculate the "adjusted arma nonexceedence probability"
# Q_e'. This probability is the probability that the projection of 
# a random point is not exceeding the critical value C(theta) given that the 
# point is outside of the transformed arma inner circle. This adjusted
# probability can be used to estimate C(theta) given that all the
# sampled points are outside the arma inner circle.
@jit(nopython=True)
def getAdjArmaNonExceedenceProbability(samplingRate, returnPeriod):
    ir = getArmaInnerRadius(samplingRate, returnPeriod)
    pe = getExceedenceProbability(samplingRate, returnPeriod)
    return 1.0 - exp(0.5 * ir * ir) * pe