conf.py

import numpy as np


def test(flipped=False):
    bt = np.arange(265, 275)
    if flipped is False:
        thr = np.array([267, 270, 273, 1, 1])
    else:
        thr = np.array([273, 270, 267, 1, 1])
    c = conf_test(bt, thr)
    print(c)


def test_dble(flipped=False):
    bt = np.arange(260, 282)
    if flipped is False:
        thr = np.array([264, 267, 270, 273, 276, 279, 1, 1])
    else:
        thr = np.array([279, 276, 273, 270, 267, 264, 1, 1])
    c = conf_test_dble(bt, thr)
    print(c)


def conf_test(rad, thr):
    '''
    Assuming a linear function between min and max confidence level, the plot below shows
    how the confidence (y axis) is computed as function of radiance (x axis).
    This case illustrates alpha < gamma, obviously in case alpha > gamma, the plot would be
    flipped.
                       gamma
    c  1                 ________
    o  |                /
    n  |               /
    f  |              /
    i  |     beta    /
    d 1/2    |....../
    e  |           /
    n  |          /
    c  |         /
    e  0________/
       |      alpha
    --------- radiance ---------->
    '''

    radshape = rad.shape
    rad = rad.reshape(np.prod(radshape))

    if thr.ndim == 1:
        thr = np.full((rad.shape[0], 4), thr[:4]).T

    coeff = np.power(2, (thr[3] - 1))
    hicut = thr[2, :]
    beta = thr[1, :]
    locut = thr[0, :]
    power = thr[3, :]
    confidence = np.zeros(rad.shape)

    alpha, gamma = np.empty(rad.shape), np.empty(rad.shape)
    flipped = np.zeros(rad.shape)

    gamma[hicut > locut] = thr[2, hicut > locut]
    alpha[hicut > locut] = thr[0, hicut > locut]
    flipped[hicut > locut] = 0
    gamma[hicut < locut] = thr[0, hicut < locut]
    alpha[hicut < locut] = thr[2, hicut < locut]
    flipped[hicut < locut] = 1

    # Rad between alpha and beta
    range_ = 2. * (beta - alpha)
    s1 = (rad - alpha) / range_
    idx = np.nonzero((rad <= beta) & (flipped == 0))
    confidence[idx] = coeff[idx] * np.power(s1[idx], power[idx])
    idx = np.nonzero((rad <= beta) & (flipped == 1))
    confidence[idx] = 1.0 - coeff[idx] * np.power(s1[idx], power[idx])

    # Rad between beta and gamma
    range_ = 2. * (beta - gamma)
    s1 = (rad - gamma) / range_
    idx = np.nonzero((rad > beta) & (flipped == 0))
    confidence[idx] = 1.0 - coeff[idx] * np.power(s1[idx], power[idx])
    idx = np.nonzero((rad > beta) & (flipped == 1))
    confidence[idx] = coeff[idx] * np.power(s1[idx], power[idx])

    # Rad outside alpha-gamma interval
    confidence[(rad > gamma) & (flipped is False)] = 1
    confidence[(rad < alpha) & (flipped is False)] = 0
    confidence[(rad > gamma) & (flipped is True)] = 0
    confidence[(rad < alpha) & (flipped is True)] = 1

    confidence[confidence > 1] = 1
    confidence[confidence < 0] = 0

    return confidence


def conf_test_dble(rad, coeffs):
    # '''
    #            gamma1                         gamma2
    #    c  1_______                               ________
    #    o  |       \                             /
    #    n  |        \                           /
    #    f  |         \                         /
    #    i  |          \   beta1       beta2   /
    #    d 1/2          \....|          |...../
    #    e  |            \                   /
    #    n  |             \                 /
    #    c  |              \               /
    #    e  0               \_____________/
    #       |             alpha1       alpha2
    #    --------------------- radiance ------------------------->
    # '''

    coeffs = np.array(coeffs)
    radshape = rad.shape
    rad = rad.reshape(np.prod(radshape))
    confidence = np.zeros(rad.shape)

    alpha1, gamma1 = np.empty(rad.shape), np.empty(rad.shape)
    alpha2, gamma2 = np.empty(rad.shape), np.empty(rad.shape)

    if coeffs.ndim == 1:
        coeffs = np.full((rad.shape[0], 7), coeffs[:7]).T

    gamma1 = coeffs[0, :]
    beta1 = coeffs[1, :]
    alpha1 = coeffs[2, :]
    alpha2 = coeffs[3, :]
    beta2 = coeffs[4, :]
    gamma2 = coeffs[5, :]
    power = coeffs[6, :]

    coeff = np.power(2, (power - 1))
    # radshape = rad.shape
    # rad = rad.reshape((rad.shape[0]*rad.shape[1]))

    # ## Find if interval between inner cutoffs passes or fails test

    # Inner region fails test

    # Value is within range of lower set of limits
    range_ = 2. * (beta1 - alpha1)
    s1 = (rad - alpha1) / range_
    idx = np.nonzero((rad <= alpha1) & (rad >= beta1) & (alpha1 - gamma1 > 0))
    confidence[idx] = coeff[idx] * np.power(s1[idx], power[idx])

    range_ = 2. * (beta1 - gamma1)
    s1 = (rad - gamma1) / range_
    idx = np.nonzero((rad >= gamma1) & (rad < beta1) & (alpha1 - gamma1 > 0))
    confidence[idx] = 1.0 - coeff[idx] * np.power(s1[idx], power[idx])

    # Value is within range of upper set of limits
    range_ = 2. * (beta2 - alpha2)
    s1 = (rad - alpha2) / range_
    idx = np.nonzero((rad > alpha1) & (rad <= beta2) & (alpha1 - gamma1 > 0))
    confidence[idx] = coeff[idx] * np.power(s1[idx], power[idx])

    range_ = 2. * (beta2 - gamma2)
    s1 = (rad - gamma2) / range_
    idx = np.nonzero((rad > alpha1) & (rad > beta2) & (alpha1 - gamma1 > 0))
    confidence[idx] = 1.0 - coeff[idx] * np.power(s1[idx], power[idx])

    # Check for value beyond function range
    confidence[(alpha1 - gamma1 > 0) & (rad > alpha1) & (rad < alpha2)] = 0
    confidence[(alpha1 - gamma1 > 0) & ((rad < gamma1) | (rad > gamma2))] = 1

    ###

    # Inner region passes test
    print("I NEED TO REVIEW THIS TO WRITE IT MORE CLEARLY")
    # FOR NOW ALPHA AND GAMMA ARE SWITCHED BECAUSE OF HOW THE ARRAYS ARE DEFINED.
    # THINK ON HOW THIS COULD BE WRITTEN SO THAT IT'S EASIER TO UNDERSTAND (AND DEBUG)
    # Value is within range of lower set of limits
    range_ = 2 * (beta1 - alpha1)
    s1 = (rad - alpha1) / range_
    idx = np.nonzero((rad > alpha1) & (rad <= gamma1) & (rad <= beta1) & (alpha1 - gamma1 <= 0))
    confidence[idx] = 1.0 - coeff[idx] * np.power(s1[idx], power[idx])

    range_ = 2 * (beta1 - gamma1)
    s1 = (rad - gamma1) / range_
    idx = np.nonzero((rad > alpha1) & (rad <= gamma1) & (rad > beta1) & (alpha1 - gamma1 <= 0))
    confidence[idx] = coeff[idx] * np.power(s1[idx], power[idx])

    # Values is within range of upper set of limits
    range_ = 2 * (beta2 - alpha2)
    s1 = (rad - alpha2) / range_
    idx = np.nonzero((rad > gamma2) & (rad < alpha2) & (rad >= beta2) & (alpha1 - gamma1 <= 0))
    confidence[idx] = 1.0 - coeff[idx] * np.power(s1[idx], power[idx])

    range_ = 2 * (beta2 - gamma2)
    s1 = (rad - gamma2) / range_
    idx = np.nonzero((rad > gamma2) & (rad < alpha2) & (rad < beta2) & (alpha1 - gamma1 <= 0))
    confidence[idx] = coeff[idx] * np.power(s1[idx], power[idx])

    confidence[(alpha1 - gamma1 <= 0) & ((rad > gamma1) | (rad < gamma2))] = 0
    confidence[(alpha1 - gamma1 <= 0) & (rad <= alpha1) & (rad >= alpha2)] = 1

    confidence[confidence > 1] = 1
    confidence[confidence < 0] = 0

    return confidence


if __name__ == "__main__":
    test_dble()