# -*- coding: utf-8 -*-
"""
Created on Sun Feb 27 13:21:40 2022

@author: Rene
"""

import pandas as pd
import numpy as np

def GT_Probs(a,b,c,d,n1,n3,alpha):
    denominator = 2*alpha*(b-a)*n3+(alpha+1)*n1*n3*(c-b)+2*n1*(d-c)
    numerator1 = 2*alpha*(b-a)*n3
    numerator2 = (alpha+1)*n1*n3*(c-b)
    numerator3 = 2*n1*(d-c)
    pi1 = numerator1/denominator
    pi2 = numerator2/denominator
    pi3 = numerator3/denominator
    df = pd.DataFrame({'pi_prob':[pi1,pi2,pi3]})
    return df

def first_branch_pdf(x,a,b,n1):
    value = 0
    if a < b:
        if (a<x) & (x<=b):
            y = (x-a)/(b-a)
            value = n1*y**(n1-1)/(b-a)
    return value 

def first_branch_cdf(x,a,b,n1):
    value = 0
    if a < b: 
        if (a<x) & (x<=b):
            y = (x-a)/(b-a)
            value = y**n1
        elif (x>=b):
            value = 1
        else:
            value = 0
    return value

def first_branch_cdf_inv(y,a,b,n1):
    if (y==0):
        value = a
    elif (y==1):
        value = b
    else:    
        value = y**(1/n1)
        value = (b-a)*value+a
    return value

def second_branch_pdf(x,b,c,alpha):
    value = 0
    if b < c:
        if (b<x) & (x<=c):
            y = (x-b)/(c-b)
            factor = 2/((alpha+1)*(c-b))
            y = (1-alpha)*y+alpha
            value = factor*y
    return value 

def second_branch_cdf(x,b,c,alpha):
    value = 0
    if b < c: 
        if (b<=x) & (x<=c):
            y = (x-b)/(c-b)
            factor1 = 2*alpha/(alpha+1)
            factor2 = (1-alpha)/(1+alpha)
            value = factor2*(y**2) + factor1*y
        elif (x>=c):
            value = 1
        else:
            value = 0
    return value

def second_branch_cdf_inv(y,b,c,alpha):
    if (alpha==1):
       x = (c-b)*y+b
    elif ((alpha>=0) & (alpha<1)):
        factor = (1+alpha)/(1-alpha)
        term2 = y+(alpha**2/(1-alpha**2))
        term1 = (factor*term2)**0.5
        y_value = (alpha/(alpha-1))+term1
        x = y_value*(c-b)+b
    elif ((alpha>1) & (alpha<=2)):
        factor = (1+alpha)/(1-alpha)
        term2 = y+(alpha**2/(1-alpha**2))
        term1 = (factor*term2)**0.5
        y_value = (alpha/(alpha-1))-term1 
        x = y_value*(c-b)+b
    return x
    

def third_branch_pdf(x,c,d,n3):
    value = 0
    if c < d:
        if (c<x) & (x<=d):
            y = (d-x)/(d-c)
            value = n3*y**(n3-1)/(d-c)
    return value 

def third_branch_cdf(x,c,d,n3):
    value = 0
    if c < d: 
        if (c<x) & (x<=d):
            y = (d-x)/(d-c)
            value = 1-y**n3
        elif (x>=d):
            value = 1
        else:
            value = 0
    return value


def third_branch_cdf_inv(y,c,d,n3):
    if (y==0):
        value = d
    elif (y==1):
        value = c
    else:    
        value = y**(1/n3)
        value = d-(d-c)*value
    return value

def GT_pdf(x,a,b,c,d,n1,n3,alpha):
    GTProbs = GT_Probs(a,b,c,d,n1,n3,alpha)
    if isinstance(x, float) or isinstance(x, int):
        n_points = 1    
        pdf_values = np.zeros(n_points)
        for i in np.arange(0,n_points):
            pdf_1 = first_branch_pdf(x,a,b,n1)
            pdf_2 = second_branch_pdf(x,b,c,alpha)
            pdf_3 = third_branch_pdf(x,c,d,n3)
            pdf = GTProbs['pi_prob'][0]*pdf_1+GTProbs['pi_prob'][1]*pdf_2+GTProbs['pi_prob'][2]*pdf_3
        pdf_values[0] = pdf
    else:   
        n_points = len(x)
        pdf_values = np.zeros(n_points)
        for i in np.arange(0,n_points):
            pdf_1 = first_branch_pdf(x[i],a,b,n1)
            pdf_2 = second_branch_pdf(x[i],b,c,alpha)
            pdf_3 = third_branch_pdf(x[i],c,d,n3)
            pdf = GTProbs['pi_prob'][0]*pdf_1+GTProbs['pi_prob'][1]*pdf_2+GTProbs['pi_prob'][2]*pdf_3
            pdf_values[i] = pdf
    return pdf_values

def GT_cdf(x,a,b,c,d,n1,n3,alpha):
    GTProbs = GT_Probs(a,b,c,d,n1,n3,alpha)
    if isinstance(x, float) or isinstance(x, int):
        n_points = 1    
        cdf_values = np.zeros(n_points)
        for i in np.arange(0,n_points):
            cdf_1 = first_branch_cdf(x,a,b,n1)
            cdf_2 = second_branch_cdf(x,b,c,alpha)
            cdf_3 = third_branch_cdf(x,c,d,n3)
            pdf = GTProbs['pi_prob'][0]*cdf_1+GTProbs['pi_prob'][1]*cdf_2+GTProbs['pi_prob'][2]*cdf_3
        cdf_values[0] = pdf
    else:   
        n_points = len(x)
        cdf_values = np.zeros(n_points)
        for i in np.arange(0,n_points):
            cdf_1 = first_branch_cdf(x[i],a,b,n1)
            cdf_2 = second_branch_cdf(x[i],b,c,alpha)
            cdf_3 = third_branch_cdf(x[i],c,d,n3)
            cdf = GTProbs['pi_prob'][0]*cdf_1+GTProbs['pi_prob'][1]*cdf_2+GTProbs['pi_prob'][2]*cdf_3
            cdf_values[i] = cdf
    return cdf_values


def GT_cdf_inv(y,a,b,c,d,n1,n3,alpha):
    GTProbs = GT_Probs(a,b,c,d,n1,n3,alpha)
    pi_1 = GTProbs['pi_prob'][0]
    pi_2 = GTProbs['pi_prob'][1]
    pi_3 = GTProbs['pi_prob'][2]
    if (isinstance(y, float) or isinstance(y, int)):
        n_points = 1    
        quantile_values = np.zeros(n_points)
        for i in np.arange(0,n_points):
            if ((0<=y) & (y<=pi_1)):
                y_value = y/pi_1
                quantile = first_branch_cdf_inv(y_value,a,b,n1)
            elif ((pi_1<=y) & (y<=(pi_1+pi_2))):
                y_value = (y-pi_1)/pi_2
                quantile = second_branch_cdf_inv(y_value,b,c,alpha)
            elif (((pi_1+pi_2)<=y) & (y<=1)):
                y_value = (1-y)/pi_3
                quantile = third_branch_cdf_inv(y_value,c,d,n3)    
            quantile_values[0] = quantile
    else:   
        n_points = len(y)
        quantile_values = np.zeros(n_points)
        for i in np.arange(0,n_points):
            if ((0<=y[i]) & (y[i]<=pi_1)):
                y_value = y[i]/pi_1
                quantile = first_branch_cdf_inv(y_value,a,b,n1)
            elif ((pi_1<=y[i]) & (y[i]<=(pi_1+pi_2))):
                y_value = (y[i]-pi_1)/pi_2
                quantile = second_branch_cdf_inv(y_value,b,c,alpha)
            elif (((pi_1+pi_2)<=y[i]) & (y[i]<=1)):
                y_value = (1-y[i])/pi_3
                quantile = third_branch_cdf_inv(y_value,c,d,n3)    
            quantile_values[i] = quantile
    return quantile_values