from orangescape.model.engine.exception.Error import FormulaArgumentError
from orangescape.model.engine.exception.ErrorMessage import *
from orangescape.model.engine.util.NumericUtil import is_a_number

from MathandTrig import SUM,COMBIN
import math
import operator
import re

def AVERAGE(argList):
    return  reduce(operator.add, argList,0)/len(argList)

def COUNT(argList):
    return len([i for i in argList if is_a_number(str(i))])

def MAX(argList):
    if(len(argList)>0):
        return max(argList)
    return 0

def MIN(argList):
    return min(argList)

def MEDIAN(argList):
    """Returns the median of a sequence of numbers."""
    argList.sort()
    length = len(argList)
    if (length % 2) != 0:
        return argList[length/2]
    else:
        return (argList[(length/2)]+argList[(length/2)-1])/2

def MODE(argList):
    # Returns the most common value in a Array Value.
    # If there are several values with the same frequency, it returns the smallest value.
    if len(argList)<=1:
        #return -2146826246.0  # #N/A Error in Excel
        raise FormulaArgumentError, Engine2005 %(len(argList),'MODE')
    setArgList=set(argList)    # Get Unique Number
    MaxCountNum = 1
    ModeNum=argList[0]
    for num in setArgList:
        cnt = argList.count(num)
        if MaxCountNum<cnt:
            MaxCountNum=cnt
            ModeNum=num
        elif MaxCountNum==cnt and num<ModeNum:
            ModeNum=num
    if MaxCountNum==1:        # If there is no Duplicate value exists, ERROR
        #return -2146826246.0  # #N/A Error in Excel
        raise FormulaArgumentError, Engine2005 %(len(argList),'MODE')
    return ModeNum

def STDEV(argList):
    if len(argList) < 2:
        raise FormulaArgumentError,Engine1003 % ('STDEV')
    #STDEV  - Standard deviation for sample
    #STEP 1 - Find the arithmetic mean (ex: ((3 + 7 + 7 + 19) / 4 = 9  )
    #STEP 2 - Subtract each value from mean,  this is deviation from mean (ex : (3 - 9) = (-9) )
    #STEP 3 - square each of the deviation,( ex : pow(-6,2) = 36 )
    #STEP 4 - Find the mean of this values [ ex :(36 + 4 + 4 + 100) /  " (4 - 1) " = 36)]
    #STEP 5 - Square root of the mean ( ex: sqrt (36) = 6)
    final_list = filter(lambda arg: is_a_number(str(arg)) ,argList)
    temp = 0.0
    avg = reduce(operator.add, argList,0)/len(argList)
    for arg in final_list:
        temp = temp + pow((arg-avg),2)
    return math.sqrt(temp / (len(final_list)-1))

def STDEVP(argList):
    #STDEV  - Standard deviation
    #STEP 1 - Find the arithmetic mean (ex: ((3 + 7 + 7 + 19) / 4 = 9  )
    #STEP 2 - Subtract each value from mean,  this is deviation from mean (ex : (3 - 9) = (-9) )
    #STEP 3 - square each of the deviation,( ex : pow(-6,2) = 36 )
    #STEP 4 - Find the mean of this values [ ex :(36 + 4 + 4 + 100) /  " (4 ) " = 36)]  In previous case we devided with ( n -1 ) This is the difference
    #STEP 5 - Square root of the mean ( ex: sqrt (36) = 6)
    avg = AVERAGE(argList)
    return math.sqrt(reduce(operator.add,[pow(arg-avg,2) for arg in argList]) / len(argList))

def VAR(argList):
    mean = SUM(argList)/len(argList)
    result = 0.0
    for arg in argList:
        result = result + (pow((arg-mean),2))
    return result/(len(argList)-1)

def VARP(argList):
    mean = SUM(argList)/len(argList)
    result = 0.0
    for arg in argList:
        result = result + (pow((arg-mean),2))
    return result/(len(argList))

def AVEDEV(argList):
    sum = 0
    temp = 0
    for arg in argList:
        sum += arg
    mean = sum/len(argList)
    for arg in argList:
        temp += abs(arg - mean)

    return temp/len(argList)

def AVERAGEA(argList):
    sum = 0.0
    final_list = filter(lambda arg: arg is not None ,argList)
    sum = reduce(lambda sum,arg: sum+arg,(filter(lambda arg: is_a_number(arg) ,final_list)),0)
    return sum/len(final_list)

def GetBeta(x, alpha, beta):
    if beta == 1.0:
        return math.pow (x, alpha)
    elif alpha == 1.0:
        return math.pow((1.0 - x), beta) - 1.0

    fEps = 1.0E-8
    if (x < (alpha+1.0)/(alpha+beta+1.0)) :
        bReflect = False
        fA = alpha
        fB = beta
    else :
        bReflect = True
        fA = beta;
        fB = alpha;
        x = 1.0 - x;
    if (x < fEps):
        cf = 0.0
    else :
        a1 = 1.0
        b1 = 1.0
        b2 = 1.0 - (fA + fB )*x/(fA+1.0)
        if (b2 == 0.0) :
            a2 = b2
            fnorm = 1.0;
            cf = 1.0;
        else :
            a2 = 1.0
            fnorm = 1.0/b2
            cf = a2*fnorm
        cfnew = 1.0
        for j in range(1,100):
            rm =  j
            apl2m = fA + 2.0*rm;
            d2m = rm*(float(fB)-float(rm))*x/((apl2m-1.0)*apl2m);
            d2m1 = -(fA+rm)*(fA+fB+rm)*x/(apl2m*(apl2m+1.0));
            a1 = (a2+d2m*a1)*fnorm;
            b1 = (b2+d2m*b1)*fnorm;
            a2 = a1 + d2m1*a2*fnorm;
            b2 = b1 + d2m1*b2*fnorm;
            if (b2 != 0.0) :
                fnorm = 1.0/b2;
                cfnew = a2*fnorm;
                if (abs(float(cf)-float(cfnew))/cf < fEps):
                   j = 101;
                else:
                   cf = cfnew;
        if (fB < fEps):
           b1 = 1.0E30
        else:
            b1 =  math.exp(float(GetLogGamma(fA))+float(GetLogGamma(fB))-float(GetLogGamma(fA+fB)))
            if b1 == 0.0:
                raise FormulaArgumentError,Engine2008 % ('BETADIST')
        cf *= math.pow(x, fA)*math.pow(1.0-x,float(fB))/(fA*b1)
    if (bReflect):
        return 1.0 - cf
    else:
        return cf

def GammaHelp(x) :
    c = [ 76.18009173, -86.50532033    , 24.01409822,-1.231739516,   0.120858003E-2, -0.536382E-5]
    if (x >= 1.0) :
        reflect = False
        x -= 1.0;
    else:
        reflect = True
        x = 1.0 - x;
    res = 1.0;
    anum = x;
    for i in xrange(0,6):
        anum += 1.0
        res += c[i]/anum
    res *= 2.506628275
    return res,x,reflect

def GetLogGamma(x):
    x = float (x)
    G,x,bReflect = GammaHelp(x);
    G = (x+0.5)*math.log(x+5.5)+math.log(G)-(x+5.5);
    if (bReflect):
        G = math.log(math.pi*x) - G - math.log(math.sin(math.pi*x));
    return G

def BETADIST(argList):
    x = argList[0]
    alpha = argList[1]
    beta = argList[2]
    if alpha <= 0 or beta <= 0:
        raise FormulaArgumentError,Engine2007 % ('BETADIST' , 'lesser',0)
    if len(argList) == 5:
        if argList[3] == None :
            fA = 0
        else:
            fA = argList[3]
        fB = argList[4]
    else:
        fA = 0
        fB = 1
    if x < fA or x > fB or fA == fB:
        raise FormulaArgumentError,Engine2009 % ('argument values', 'BETADIST' )
    kum = True
    scale = (x - fA) /  (fB - fA)
    return abs(GetBeta (scale, alpha, beta))

def BINOMDIST(argList):
    x = int(argList[0])
    n = int(argList[1])
    p = argList[2]
    kum = argList[3]

    if  n < 0.0 or x < 0.0 or x > n or p < 0.0 or p > 1.0 :
        raise FormulaArgumentError,Engine2009 % ('argument values', 'BETADIST' )

    if  kum :
        if (n == x):
            res = 1.0
        else:
            q = 1.0 - p
            factor = math.pow(q, n)
            if (factor == 0.0) :
                factor = pow(p, n)
                if (factor == 0.0):
                    raise Exception
                else:
                    res = 1.0 - factor
                    max = (n - x) - 1
                    i = 0
                    while i < max and i > factor:
                        factor *= (n-i)/(i+1)*q/p
                        res -= factor
                        i +=1
                    if (res < 0.0):
                        res = 0.0
            else:
                res = factor
                max =  x
                i = 0
                while i < max and factor > 0.0 :
                    factor *= float(n-i)/(i+1)*p/q
                    res += factor
                    i +=1
    else:
        q = 1.0 - p
        factor = math.pow(q, n)
        if (factor == 0.0):
            factor = math.pow(p, n)
            if (factor == 0.0):
                raise FormulaArgumentError,Engine2009 % ('argument value', 'BETADIST')
            else:
                max = (n - x)
                while (i < max and factor > 0.0):
                    factor *= float(n-i)/(i+1)*q/p
                    res = factor
                    i += 1
        else:
            max =  x
            i =0
            while i < max and factor > 0.0:
                factor *= float(n-i)/(i+1)*float(p)/q
                res = factor
                i += 1
    return res

def gaussinv (xx):
    x = float (xx)
    q = x-0.5

    if (abs(q)<=0.425):
        t=0.180625-q*q;
        z=q*(((((((t*2509.0809287301226727+33430.575583588128105)*t+67265.770927008700853)*t+45921.953931549871457)*t+13731.693765509461125)*t+1971.5909503065514427)*t+133.14166789178437745)*t+3.387132872796366608)/(((((((t*5226.495278852854561+28729.085735721942674)*t+39307.89580009271061)*t+21213.794301586595867)*t+5394.1960214247511077)*t+687.1870074920579083)*t+42.313330701600911252)*t+1.0)
    else:
        if (q>0):  t=1-x
        else :   t=x
        t = math.sqrt(-math.log(t))
        if(t<=5.0):
              t+=-1.6
              z=(((((((t*7.7454501427834140764e-4+0.0227238449892691845833)*t+0.24178072517745061177)*t+1.27045825245236838258)*t+3.64784832476320460504)*t+5.7694972214606914055)*t+4.6303378461565452959)*t+1.42343711074968357734)/(((((((t*1.05075007164441684324e-9+5.475938084995344946e-4)*t+0.0151986665636164571966)*t+0.14810397642748007459)*t+0.68976733498510000455)*t+1.6763848301838038494)*t+2.05319162663775882187)*t+1.0)
        else:
              t+=-5.0
              z=(((((((t*2.01033439929228813265e-7+2.71155556874348757815e-5)*t+0.0012426609473880784386)*t+0.026532189526576123093)*t+0.29656057182850489123)*t+1.7848265399172913358)*t+5.4637849111641143699)*t+6.6579046435011037772)/(((((((t*2.04426310338993978564e-15+1.4215117583164458887e-7)*t+1.8463183175100546818e-5)*t+7.868691311456132591e-4)*t+0.0148753612908506148525)*t+0.13692988092273580531)*t+0.59983220655588793769)*t+1.0)
        if(q<0.0): z=-z
    return z

def CONFIDENCE(argList):
    alpha = argList[0]
    sigma = argList[1]
    n = argList[2]

    if  sigma <= 0.0 or alpha <= 0.0 or alpha >= 1.0 or n < 1:
        raise FormulaArgumentError,Engine2009 % ('arguments','CONFIDENCE')

    g = gaussinv (1.0 - alpha / 2.0)
    return  g * sigma / math.sqrt (n)

def COUNTA(argList):
    return len(filter(lambda arg: arg != None , argList ))

def COUNTBLANK(argList):
    return len(filter(lambda arg: arg != None and arg!= '',argList))

#def COUNTIF(argList):
#    criteria =  argList.pop()
#    if re.search('\?',criteria) and not re.search('[a-bA-B0-9*]',criteria):
#        n = len(criteria)
#        criteria = '[a-z]'+'{'+`n`+'}'
#        return len(filter(lambda arg:(re.match(criteria,str(arg))),argList))
#    elif re.search('[*?]',criteria) :
#        criteria = criteria.replace('*', '.*')
#        criteria = criteria.replace('?', '.')
#        criteria = criteria.lower()
#        text_list = filter(lambda arg:type(arg) != float,argList)
#        return len(filter(lambda arg:(re.search(criteria,arg.lower())),text_list))
#
#    elif re.search('[<>,>,<,<=,>=]\Z', criteria):
#        number_list = filter(lambda arg:type(arg) == float,argList)
#        return len(filter(lambda arg: eval(criteria+`arg`), number_list))
#
#    elif re.search('[<>,>,<,<=,>=]^', criteria):
#        number_list = filter(lambda arg:type(arg) == float,argList)
#        return len(filter(lambda arg: eval(`arg`+criteria), number_list))
#
#    else:
#        return len(filter(lambda arg: eval(`criteria`.lower() +'=='+ `arg`.lower()), argList))

def DEVSQ(argList):
    avg =  (reduce(operator.add, argList,0.0))/len(argList)
    return reduce(operator.add,map(lambda arg: abs(arg - avg)**2,argList) ,0.0)

def EXPONDIST (argList):
    x = argList[0]
    lambdas = argList[1]
    kum = argList[2]
    result = 0.0
    if not(lambdas > 0.0) :
        raise FormulaArgumentError,Engine2017 % (2,'EXPONDIST',lambdas)
    if not(x >= 0.0) :
        raise FormulaArgumentError,Engine2017 % (2,'EXPONDIST',x)

    ex = math.exp ((lambdas * -1)* x)
    if (kum == 0)  :
        if not (x < 0.0) :
          result = lambdas * ex
    else :
        if (x > 0.0) :
            result = 1.0 - ex
    return result

def FISHER(argList):
    x = argList[0]
    if x >= 1 or x<= -1:
        raise FormulaArgumentError,Engine2009 % ('argument value', 'FISHER')
    return 0.5 * math.log((1+x)/(1-x))

def FISHERINV(argList):
    if argList[0] > 25:
        return 1.0
    x = math.exp(2*argList[0])
    return (x-1)/(x+1)

def GEOMEAN(argList):
    if len(filter(lambda arg: arg <= 0,argList)) > 0: raise FormulaArgumentError,Engine2009 %('argument value','GEOMEAN')
    prod = reduce(operator.mul, argList,1)
    return math.pow(prod, (1.0 / len(argList)))

def HARMEAN(argList):
    if len(filter(lambda arg: arg <= 0,argList)) > 0: raise FormulaArgumentError,Engine2009 %('argument value','GEOMEAN')
    return 1 /(reduce(operator.add, map(lambda arg : 1.0/arg,argList)) * 1.0/len(argList))

def HYPGEOMDIST(argList):
    sample_s = int(argList[0])
    number_sample = int(argList[1])
    population_s = int(argList[2])
    number_population = int(argList[3])
    if  sample_s < 0 or number_sample <= 0 or population_s <= 0 or number_population <= 0 :
        raise FormulaArgumentError,Engine2007 % ('HYPGEOMDIST','lesser',0)
    if ( sample_s > population_s ) :
        raise FormulaArgumentError,Engine2017 % (1 ,'HYPGEOMDIST',sample_s)
    if( number_sample > number_population ):
        raise FormulaArgumentError,Engine2017 % (2,'HYPGEOMDIST' ,number_sample)
    if not sample_s < number_sample or not sample_s < population_s:
        raise FormulaArgumentError,Engine2017 % (1,'HYPGEOMDIST',sample_s)

    d1 = COMBIN([population_s,sample_s])
    temp1 = number_population-population_s
    temp2 = number_sample - sample_s
    if  temp1 < 1: temp1 = 1
    if temp2 < 1: temp2 = 1
    d2 = COMBIN([temp1,temp2])
    d3 = COMBIN([ number_population,number_sample])
    return float(d1*d2) / d3

def LARGE(argList):
    k = int(math.ceil(argList[len(argList)-1]))
    array = argList[0:len(argList)-1]
    array.sort(reverse=True)
    if k > len(array) or (k <= 0):
        raise FormulaArgumentError,Engine2017 % (2,'LARGE',k )
    return array[k-1]

def MAXA(argList):
    number_list = filter(lambda arg:is_a_number(str(arg)), argList)
    if len(number_list) < len(argList):
        number_list = filter(lambda number: number >= 0, number_list)
    if len(number_list) == 0 :number_list.append(0)
    number_list = map(lambda number: float(number), number_list)
    return max(number_list)

def MINA(argList):
    number_list = filter(lambda arg: is_a_number(str(arg)), argList)
    # To achieve excel behaviour which is if you give string as a argument value then it will give the mina value of negative numbers,and vice versa in maxa
    if len(number_list) < len(argList):
        number_list = filter(lambda number: number < 0, number_list)
    if len(number_list) == 0: number_list.append(0)
    return min(map(lambda arg: float(arg) ,number_list))

def NEGBINOMDIST (argList):
    x = int(argList[0])
    r = int(argList[1])
    p = argList[2]

    if ( r < 1.0 or x < 0.0 or p < 0.0 or p >= 1.0 ):
        raise FormulaArgumentError,Engine2009 % ('argument value', 'NEGBINOMDIST')

    q = 1.0 - p
    res = math.pow(p,r)
    i = 0
    while (i < x) :
        res *= (i+r)/(i+1.0)*q
        i += 1

    return res

def PERMUT(argList):
    number = int(argList[0])
    number_chosen = int(argList[1])
    if number == 0 and number_chosen == 0:
        return 1
    if number <= 0:
        raise FormulaArgumentError,Engine2017 % (1,'PERMUT',number)
    if number_chosen < 0:
        raise FormulaArgumentError,Engine2017 % (1,'PERMUT',number_chosen)
    if number < number_chosen:
        raise FormulaArgumentError,Engine2007 % ('PERMUT','greater',number)
    if number_chosen == 0:
        return 1
    return reduce(operator.mul, xrange(number ,number - number_chosen , -1))

def POISSON(argList):
    x = int(argList[0])
    lambdaa = argList[1]
    kum = argList[2]
    if x == 0 and lambdaa == 0 and kum == 0 :
        return 1
    if  lambdaa < 0.0 or x < 0.0 :
        raise FormulaArgumentError,Engine2007 % ('POISSION','lesser',0)

    ex = math.exp (lambdaa * -1)

    if kum == 0:
        if lambdaa == 0:
          result = 0
        else:
            result = ex * math.pow (lambdaa, x) / math.factorial (x)
    else :
        if lambdaa == 0:
          result =  1
        else:
            result = 1.0
            fFak = 1.0
            nEnd = int(x)
            for i in xrange(1,nEnd+1):
                fFak = fFak *i
                result += math.pow(lambdaa, i) / fFak
            result = result * ex
    return result

def PERCENTILE(argList):
    #Ref : http://www.haiweb.org/medicineprices/manual/quartiles_iTSS.pdf
    p = argList[len(argList) - 1]
    array = argList[0:-1]
    if p < 0 or p > 1:
        raise FormulaArgumentError,Engine2017 % (2,'PERCENTILE',p)
    if p == 1:
        return max(array)
    if len(array) == 1:
        return argList[0]
    n = len(array)
    array.sort()
    j = int((n - 1) * p)
    g =  float((n - 1) * p) - j
    return (1-g) *array[j] + g* array[j+1]

def QUARTILE(argList):
    #0 Minimum value
    #1 First quartile (25th percentile)
    #2 Median value (50th percentile)
    #3 Third quartile (75th percentile)
    #4 Maximum value
    quart = int(argList[len(argList) - 1])
    array = argList[:-1]
    if quart not in (0,1,2,3,4):
        raise FormulaArgumentError,Engine2017%(2,'QUARTILE',quart)
    if quart == 0:
        return min(array)
    if len(array) == 1:
        return array[0]
    if quart == 1:
        p = 0.25
        n = len(array)
        array.sort()
        j = int((n - 1) * p)
        g =  float((n - 1) * p) - j
        return (1-g) *array[j] + g* array[j+1]
    if quart == 2:
        p = 0.50
        n = len(array)
        array.sort()
        j = int((n - 1) * p)
        g =  float((n - 1) * p) - j
        return (1-g) *array[j] + g* array[j+1]
    if quart == 3:
        p = 0.75
        n = len(array)
        array.sort()
        j = int((n - 1) * p)
        g =  float((n - 1) * p) - j
        return (1-g) *array[j] + g* array[j+1]
    if quart == 4:
        return max(array)

def SKEW(argList):
    n = len(argList)
    stdDev = STDEV(argList)
    mean = AVERAGE(argList)
    if len(filter(lambda arg: arg == 0,argList)) == len(argList):
        raise FormulaArgumentError,Engine2019 % (len(argList),'SKEW',0)
    #Formula for SKEW is n / (n-1) / (n-2) * SUM(((x-mean)/s)3)
    #So I devided the formula into 2 parts t1 will calculate 'n / (n-1) / (n-2)' ie first part and t2 will calculate 2nd part

    t1 = n/((n-1.0)*(n-2.0))

    def sum(arg,mean,sd):
        return math.pow(((arg - mean)/stdDev) , 3)

    t2 = float(reduce(operator.add, [sum(arg,mean,stdDev) for arg in argList]))
    return t1 * t2
def SMALL(argList):
    k = int(math.floor(argList[len(argList)-1]))
    array = argList[0:len(argList)-1]
    if k <= 0 or k > len(array):
        raise FormulaArgumentError,Engine2017 % (2,'SMALL',k)
    array.sort()
    return array[k-1]

def STANDARDIZE(argList):
    x = argList[0]
    mean = argList[1]
    stdDev = argList[2]
    if stdDev <= 0:
        raise FormulaArgumentError,Engine2017 % (3,'STANDARDIZE ',stdDev)
    return (x-mean)/stdDev

def STDEVA(argList):
    if len(argList) == 1:
        return 0
    def changeArgumentValue(arg):
        index = argList.index(arg)
        argList[index] = 0
    strList = filter(lambda arg: not is_a_number(str(arg)),argList)
    [changeArgumentValue(arg) for arg in strList]

    temp = 0.0
    avg = reduce(operator.add, argList)/len(argList)
    for arg in argList:
        if not is_a_number(str(arg)): arg = 0.0
        temp = temp + pow((arg-avg),2)
    return round(math.sqrt(temp / (len(argList)-1)),12)

def TDIST(argList):
    X  = argList[0]
    degrees_freedom = int(argList[1])
    tails = int(argList[2])

    if degrees_freedom <  1 :
        raise FormulaArgumentError,Engine2017 % (2,'TDIST',degrees_freedom)
    if tails not in (1,2):
        raise FormulaArgumentError,Engine2017 % (3,'TDIST',tails)
    if X < 0:
        raise FormulaArgumentError,Engine2017 % (1,'TDIST',X)

    arg = float(degrees_freedom)/(degrees_freedom + X * X)

    R = GetBeta (arg,(degrees_freedom/2.0) ,0.5) * 0.5

    if tails == 1 :
        return R
    return R * 2

def WEIBULL(argList):
    # Note: This function gives 6 digit accuracy with excel
    X =argList[0]
    alpha = argList[1]
    beta = argList[2]
    cumulative = argList[3]

    if X < 0:
        raise FormulaArgumentError,Engine2017 % (1,'WEIBULL',X)
    if alpha <= 0 :
        raise FormulaArgumentError,Engine2017 % (2,'WEIBULL',alpha)
    if beta <= 0 :
        raise FormulaArgumentError,Engine2017 % (3,'WEIBULL',beta)

    try:
        if cumulative == 0 :
            return (alpha/pow(beta,alpha))*pow(X,(alpha - 1)) * math.exp((-1*pow((X/beta),alpha)))
        return 1- math.exp(-1 * (pow((X/beta),alpha)))
    except OverflowError:
        raise FormulaArgumentError,Engine2008 % ('WEIBULL')

def VARA(argList):
    def changeArgumentValue(arg):
        index = argList.index(arg)
        argList[index] = 0
    strList = filter(lambda arg: not is_a_number(str(arg)),argList)
    [changeArgumentValue(arg) for arg in strList]

    if len(argList) < 2:
        raise FormulaArgumentError,Engine2010 % ('Not enough','VARA')
    mean = SUM(argList)/ len(argList)
    result = 0.0
    for arg in argList:
        result = result + (pow((arg-mean),2))
    return result/(len(argList)-1)

def VARPA(argList):
    if len(argList) == 1:
        return 0
    def changeArgumentValue(arg):
        index = argList.index(arg)
        argList[index] = 0
    strList = filter(lambda arg: not is_a_number(str(arg)),argList)
    [changeArgumentValue(arg) for arg in strList]

    mean = SUM(argList)/ len(argList)
    result = 0.0
    for arg in argList:
        result = result + (pow((float(arg)-mean),2))
    return result/(len(argList))

def STDEVPA(argList):
    if len(argList) == 1:
        return 0
    def changeArgumentValue(arg):
        index = argList.index(arg)
        argList[index] = 0
    strList = filter(lambda arg: not is_a_number(str(arg)),argList)
    [changeArgumentValue(arg) for arg in strList]

    avg = AVERAGE(argList)
    return math.sqrt(reduce(operator.add,[pow(arg-avg,2) for arg in argList]) / len(argList))

def TRIMMEAN(argList):
    percent = argList[len(argList)-1]
    array   = argList[0:len(argList)-1]

    Count = len(array)
    if percent < 0 or percent >= 1 :
        raise FormulaArgumentError,Engine2017 % (2,'TRIMMEAN',percent)
    if Count == 0:
        raise FormulaArgumentError,Engine2009 % ('argument value' ,'TRIMMEAN')
    trim = math.floor((percent * Count)/2)
    array.sort()
    res = 0.0
    item = 0
    for i in xrange(int(0+trim),int(Count-trim)):
        res += array[i]
        item += 1

    return res/item


#------------------- TINV ----------------------------------------


# The largest argument for which <code>logGamma(x)</code> is representable in the machine.

# Relative machine precision.
EPS = 2.22e-16
# The smallest positive floating-point number such that 1/xminin is machine representable.
XMININ = 2.23e-308
# Square root of 2 * pi
SQRT2PI = 2.5066282746310005024157652848110452530069867406099
LOGSQRT2PI = math.log(SQRT2PI);
# Rough estimate of the fourth root of logGamma_xBig
lg_frtbig = 2.25e76
pnt68 = 0.6796875
# lower value = higher precision
PRECISION = 4.0*EPS
ACCURACY = 10**-7
MAX_ITERATIONS = 10000

def findRoot(value, x_low, x_high, function):
    """Use the bisection method to find root such that function(root) == value."""

    guess = (x_high + x_low) / 2.0
    v = function(guess)
    difference = v - value
    i = 0
    while abs(difference) > ACCURACY and i < MAX_ITERATIONS:
        i += 1
        if difference > 0:
            x_high = guess
        else:
            x_low = guess

        guess = (x_high + x_low) / 2.0
        v = function(guess)
        difference = v - value

    return guess

def betaFraction(x, p, q):
    """Evaluates of continued fraction part of incomplete beta function.

    Based on an idea from Numerical Recipes (W.H. Press et al, 1992)."""

    sum_pq  = p + q
    p_plus  = p + 1.0
    p_minus = p - 1.0
    h = 1.0-sum_pq*x/p_plus;
    if abs(h) < XMININ:
        h = XMININ
    h = 1.0/h
    frac = h
    m = 1
    delta = 0.0
    c = 1.0

    while m <= MAX_ITERATIONS and abs(delta-1.0) > PRECISION:
        m2 = 2*m

        # even index for d
        d=m*(q-m)*x/((p_minus+m2)*(p+m2))
        h=1.0+d*h
        if abs(h) < XMININ: h=XMININ
        h=1.0/h;
        c=1.0+d/c;
        if abs(c) < XMININ: c=XMININ
        frac *= h*c;

        # odd index for d
        d = -(p+m)*(sum_pq+m)*x/((p+m2)*(p_plus+m2))
        h=1.0+d*h
        if abs(h) < XMININ: h=XMININ;
        h=1.0/h
        c=1.0+d/c
        if abs(c) < XMININ: c = XMININ
        delta=h*c
        frac *= delta
        m += 1

    return frac
lg_d1 = -0.5772156649015328605195174
lg_d2 = 0.4227843350984671393993777
lg_d4 = 1.791759469228055000094023
lg_p1 = [ 4.945235359296727046734888, 201.8112620856775083915565, 2290.838373831346393026739, 11319.67205903380828685045,
          28557.24635671635335736389, 38484.96228443793359990269, 26377.48787624195437963534, 7225.813979700288197698961 ]

lg_q1 = [ 67.48212550303777196073036, 1113.332393857199323513008, 7738.757056935398733233834, 27639.87074403340708898585,
          54993.10206226157329794414, 61611.22180066002127833352, 36351.27591501940507276287, 8785.536302431013170870835 ]

lg_p2 = [ 4.974607845568932035012064, 542.4138599891070494101986, 15506.93864978364947665077, 184793.2904445632425417223,
          1088204.76946882876749847,  3338152.967987029735917223, 5106661.678927352456275255, 3074109.054850539556250927 ]

lg_q2 = [ 183.0328399370592604055942, 7765.049321445005871323047, 133190.3827966074194402448, 1136705.821321969608938755,
          5267964.117437946917577538, 13467014.54311101692290052, 17827365.30353274213975932, 9533095.591844353613395747 ]

lg_p4 = [ 14745.02166059939948905062, 2426813.369486704502836312, 121475557.4045093227939592, 2663432449.630976949898078,
          29403789566.34553899906876, 170266573776.5398868392998, 492612579337.743088758812,  560625185622.3951465078242 ]

lg_q4 = [ 2690.530175870899333379843, 639388.5654300092398984238, 41355999.30241388052042842, 1120872109.61614794137657,
          14886137286.78813811542398, 101680358627.2438228077304, 341747634550.7377132798597, 446315818741.9713286462081 ]

lg_c = [ -0.001910444077728, 8.4171387781295e-4, -5.952379913043012e-4, 7.93650793500350248e-4, -0.002777777777777681622553,
         0.08333333333333333331554247, 0.0057083835261 ]

LOG_GAMMA_X_MAX_VALUE = 2.55e305

def logGamma(x):
    y = float(x)
    if y < 0.0 or y > LOG_GAMMA_X_MAX_VALUE:
        # Bad arguments
        return float("inf")

    if y <= EPS:
        return -math.log(y)

    if y <= 1.5:
        if (y < pnt68):
            corr = -math.log(y)
            xm1 = y
        else:
            corr = 0.0;
            xm1 = y - 1.0;

        if y <= 0.5 or y >= pnt68:
            xden = 1.0;
            xnum = 0.0;
            for i in xrange(8):
                xnum = xnum * xm1 + lg_p1[i];
                xden = xden * xm1 + lg_q1[i];
            return corr + xm1 * (lg_d1 + xm1 * (xnum/xden));
        else:
            xm2 = y - 1.0;
            xden = 1.0;
            xnum = 0.0;
            for i in xrange(8):
                xnum = xnum * xm2 + lg_p2[i];
                xden = xden * xm2 + lg_q2[i];
            return corr + xm2 * (lg_d2 + xm2 * (xnum / xden));

    if (y <= 4.0):
        xm2 = y - 2.0;
        xden = 1.0;
        xnum = 0.0;
        for i in xrange(8):
            xnum = xnum * xm2 + lg_p2[i];
            xden = xden * xm2 + lg_q2[i];
        return xm2 * (lg_d2 + xm2 * (xnum / xden));

    if y <= 12.0:
        xm4 = y - 4.0;
        xden = -1.0;
        xnum = 0.0;
        for i in xrange(8):
            xnum = xnum * xm4 + lg_p4[i];
            xden = xden * xm4 + lg_q4[i];
        return lg_d4 + xm4 * (xnum / xden);

    assert y <= lg_frtbig
    res = lg_c[6];
    ysq = y * y;
    for i in xrange(6):
        res = res / ysq + lg_c[i];
    res /= y;
    corr = math.log(y);
    res = res + LOGSQRT2PI - 0.5 * corr;
    res += y * (corr - 1.0);
    return res

def logBeta(p, q):
    assert p > 0
    assert q > 0
    if p <= 0 or q <= 0 or p + q > LOG_GAMMA_X_MAX_VALUE:
        return 0

    return logGamma(p)+logGamma(q)-logGamma(p+q)

def incompleteBeta(x, p, q):
    assert 0 <= x <= 1
    assert p > 0
    assert q > 0

    # Range checks to avoid numerical stability issues?
    if x <= 0.0:
        return 0.0
    if x >= 1.0:
        return 1.0
    if p <= 0.0 or q <= 0.0 or (p+q) > LOG_GAMMA_X_MAX_VALUE:
        return 0.0

    beta_gam = math.exp(-logBeta(p,q) + p*math.log(x) + q*math.log(1.0-x))
    if x < (p+1.0)/(p+q+2.0):
        return beta_gam*betaFraction(x, p, q)/p
    else:
        return 1.0-(beta_gam*betaFraction(1.0-x,q,p)/q)

def StudentTCDF(degree_of_freedom, X):

    A = 0.5 * incompleteBeta(float(degree_of_freedom)/(degree_of_freedom+X*X), 0.5*degree_of_freedom, 0.5)
    if X > 0:
        return 1 - A
    return A

def InverseStudentT(degree_of_freedom, probability):
    assert 0 <= probability <= 1

    if probability == 1:
        return float("inf")
    if probability == 0:
        return float("-inf")
    if probability == 0.5:
        return 0.0

    def f(x):
        return StudentTCDF(degree_of_freedom, x)

    return findRoot(probability, -10**4, 10**4, f)

def TINV(argList):
    # Note: This function gives 6 digit accuracy with excel
    #Ref : http://evanjones.ca/statistics.py
    """p: 1-confidence (eg. 0.05 = 95% confidence)"""
    p = argList[0]
    degree_of_freedom = int(argList[1])
    if p <= 0 or  p > 1:
        raise FormulaArgumentError,Engine2017 % (1,'TINV',p)
    if degree_of_freedom < 1:
        raise FormulaArgumentError,Engine2017 % (2,'TINV',degree_of_freedom)

    confidence = 1 - p
    return InverseStudentT(degree_of_freedom, (1.0+confidence)/2.0)

def NORMDIST(argList):
    # Ref :http://unk1911.blogspot.com/2005/10/php-normal-cumulative-distribution.html
    # Note: This function gives 6 digit accuracy with excel
    X = argList[0]
    mean = argList[1]
    stdDev = argList[2]
    cumulative = argList[3]
    if stdDev <= 0:
        raise FormulaArgumentError,Engine2017 % (3,'NORMDIST',stdDev)

    if cumulative == 0:
        p1 = 1.0 / (math.sqrt(2* math.pi) * stdDev)
        p2 = math.exp(-1*(pow((X - mean),2.0) / (2*pow(stdDev,2))))
        return p1*p2
    else:
        sigma = stdDev
        res = 0
        x = (X - mean)/sigma

        if (x == 0) :
            res = 0.5;
        else :
            oor2pi = 1.0 / (math.sqrt(2 * math.pi))
            t = 1.0 / (1 + 0.2316419 * abs(x))
            t *= oor2pi* math.exp(-0.5 * x * x)* (0.31938153 + t* (-0.356563782 + t* (1.781477937 + t* (-1.821255978 + t * 1.330274429))))

            if (x >= 0) :
              res = 1 - t
            else:
              res = t
        return res

def norminv_sub(q):
    if(q == 0.50):
        return 0.0

    q = 1.0 - q

    if q > 0.0 and q < 0.5 :
        p =  q
    else:
        p = 1.0 - q

    t = math.sqrt(abs(math.log (1.0/math.pow(p, 2.0))))

    c0 = 2.515517
    c1 = 0.802853
    c2 = 0.010328

    d1 = 1.432788
    d2 = 0.189269
    d3 = 0.001308

    X = t - (c0 + c1*t + c2*pow(t, 2.0)) /(1.0 + d1*t + d2*pow(t, 2.0) + d3*pow(t,3.0))

    if q > 0.5 :
        X *= -1.0
    return X

def NORMINV(argList):
    # Ref :http://unk1911.blogspot.com/search?q=norminv
    # Note : It gives only 3 decimal point accuracy when compared with excel

    q = argList[0]
    mu = argList[1]
    sigma = argList[2]

    if q <= 0 or q >= 1:
        raise FormulaArgumentError,Engine2017 % (1,'NORMINV',q)
    if sigma <= 0 :
        raise FormulaArgumentError,Engine2017 % (2,'NORMINV',sigma)

    return mu + sigma * norminv_sub(q)

def NORMSDIST(argList):
    z = argList[0]
    normDistList = [z,0,1,1]
    return NORMDIST(normDistList)

def NORMSINV(argList):
    probability = argList[0]
    if probability <= 0 or probability >= 1:
        raise FormulaArgumentError,Engine2017 % (1,'NORMSINV',probability)
    return NORMINV([probability,0,1])

def CORREL(argList):
    if len(argList[0]) != len(argList[1]) :
        raise FormulaArgumentError,Engine2022 % ('CORREL')
    avg1 = (reduce(lambda a,b: a + b, argList[0])) / float(len(argList[0]))
    avg2 = (reduce(lambda a,b: a + b, argList[1])) / float(len(argList[1]))
    top = 0.0
    bottom1 = 0
    bottom2 = 0
    for i in xrange(len(argList[0])):
        top += ((argList[0][i] - avg1) * ( argList[1][i] - avg2))
        bottom1 += (argList[0][i] - avg1)**2
        bottom2 += (argList[1][i] - avg2)**2
    bottom  = math.sqrt( bottom1 * bottom2 )
    if bottom == 0 or top == 0.0:
        raise FormulaArgumentError,Engine2009 % ('arrays','CORREL')
    return top / bottom

def gamma(x):

    #Ref : http://www.crbond.com/download/gamma.cpp

    g = [1.0, 0.5772156649015329, -0.6558780715202538, -0.420026350340952e-1, 0.1665386113822915,-0.421977345555443e-1, -0.9621971527877e-2, 0.7218943246663e-2, -0.11651675918591e-2, -0.2152416741149e-3, 0.1280502823882e-3, -0.201348547807e-4,    -0.12504934821e-5, 0.1133027232e-5,    -0.2056338417e-6, 0.6116095e-8,    0.50020075e-8, -0.11812746e-8, 0.1043427e-9, 0.77823e-11, -0.36968e-11,    0.51e-12, -0.206e-13, -0.54e-14, 0.14e-14]

    if (x > 171.0) :return 1e308
    if  x == int(x) :
        if (x > 0.0) :
            ga = 1.0
            for i in xrange(2,int(x)): #(i=2 ; i<x; i++):
                ga *= i
        else:
            ga = 1e308
    else:
        if ( abs(x) > 1.0):
            z = abs(x)
            m = int(z)
            r = 1.0
            for k in xrange(1,m+1) : #(k=1;k<=m;k++)
                r *= (z-k)
            z -= m
        else:
            z = x
        gr = g[24]
        for k in xrange(23,-1,-1): #(k=23;k>=0;k--)
            gr = gr*z+g[k]
        ga = 1.0 / (gr*z)
        if (abs(x) > 1.0):
            ga *= r
            if (x < 0.0):
                ga = (-1 * math.pi) /(x*ga*sin(math.pi*x))
    return ga

def GAMMALN(argList):
    if argList[0] <= 0:
        raise FormulaArgumentError,Engine2017 % (1,'GAMMALN',argList[0])
    return math.log(gamma(argList[0]))

def KURT(argList):
    """ Formula for this function, present in Excel documentation. Here, I have just implemented the formula"""
    n = float(len(argList))
    if n < 4:
        raise FormulaArgumentError, Engine2010 % ('You have passed less than four arguments','KURT')

    p1 = (n*(n+1.0))/((n-1.0)*(n-2.0)*(n-3.0))

    stddev = STDEV(argList)
    avg = reduce(operator.add, argList,0)/len(argList)

    p2 = reduce(operator.add ,  map(lambda arg: pow( ( (arg - avg) / stddev  ) ,4) , argList) )

    p3 = (3.0*pow((n-1),2.0))/((n-2.0)*(n-3.0))
    return (p1*p2)-p3

def LOGNORMDIST(argList):
    X = argList[0]
    mean = argList[1]
    stdDev = argList[2]
    cumulative = 1
    if X <= 0:
        raise FormulaArgumentError,Engine2017 % (1,'LOGNORMDIST',X)
    if stdDev <= 0:
        raise FormulaArgumentError,Engine2017 % (3,'LOGNORMDIST',stdDev)
    x = math.log(X)
    return NORMDIST([x,mean,stdDev,cumulative])

def PEARSON(argList):
    #Ref: Excel
    array1 = argList[0]
    array2 = argList[1]
    n = len(array1)
    if n == 0:
        raise FormulaArgumentError, Engine2010 % ('First array has length zero which you have passed', 'PEARSON')
    elif n != len(array2):
        raise FormulaArgumentError, Engine2022 % ('PEARSON')
    x_mean = reduce(operator.add, array1) / n
    y_mean = reduce(operator.add, array2) / n

    top = reduce(operator.add, [(array1[i] - x_mean)*(array2[i] - y_mean) for i in xrange(0,n)])
    p1 = reduce(operator.add, [(array1[i]-x_mean)**2 for i in xrange(0,n)])
    p2 = reduce(operator.add, [(array2[i]-y_mean)**2 for i in xrange(0,n)])
    bottom = math.sqrt(p1*p2)
    return top/bottom

def RSQ(argList):
    #Ref: Excel
    array1 = argList[0]
    array2 = argList[1]
    n = len(array1)
    if n == 0:
        raise FormulaArgumentError, Engine2010 % ('First array has length zero which you have passed', 'RSQ')
    elif n != len(array2):
        raise FormulaArgumentError, Engine2022 % ('RSQ')
    return (PEARSON(argList))**2.0

def RANK(argList):
    #Ref: Excel
    number = argList[0]
    ref = argList[1]
    if len(argList)>2: order = argList[2]
    else: order = 0
    if len(filter (lambda i : number == i,ref)) == 0:
        raise FormulaArgumentError,Engine2017 % (1,'RANK',number)
    if order == 0: ref.sort(reverse=True)
    else: ref.sort(reverse=False)

    return ref.index(number)+1

def SLOPE(argList):
    #Ref: Excel
    knownys = argList[0]
    knownxs = argList[1]
    length = len(knownys)

    if len(knownxs)<2 or length<2:
        raise FormulaArgumentError,Engine2010 % ('Invalid Arguments', 'SLOPE')
    elif len(knownxs) != length:
        raise FormulaArgumentError,Engine2022 % ('SLOPE')

    y_mean = reduce(operator.add,knownys) / length
    x_mean = reduce(operator.add,knownxs) / length
    top = reduce(operator.add,[(knownxs[i] - x_mean)*(knownys[i] - y_mean) for i in xrange(0,length) ])
    bottom = reduce(operator.add, [(knownxs[i] - x_mean)**2 for i in xrange(0,length)])

    return top/bottom

def STEXY(argList):
    #Ref: Excel
    knownys = argList[0]
    knownxs = argList[1]
    length = len(knownys)

    if len(knownxs)<3 or length<3:
        raise FormulaArgumentError,Engine2010 % ('Invalid Arguments', 'STEXY')
    elif len(knownxs) != length:
        raise FormulaArgumentError,Engine2022 % ('STEXY')

    y_mean = reduce(operator.add,knownys) / length
    x_mean = reduce(operator.add,knownxs) / length
    top = reduce(operator.add,[(knownxs[i] - x_mean)*(knownys[i] - y_mean) for i in xrange(0,length) ])
    bottom = reduce(operator.add, [(knownxs[i] - x_mean)**2 for i in xrange(0,length)])
    p1 =  reduce(operator.add, [(knownys[i] - y_mean)**2 for i in xrange(0,length)])
    return math.sqrt( ( 1.0 / (length-2) ) * ( p1-(top**2 / bottom) ) )

def ZTEST(argList):
    #Ref: Excel
    #Note: It gives only have six digit accuracy with excel result
    array = argList[0]
    mean = argList[1]
    n = len(array)
    if n == 0:
        raise FormulaArgumentError, Engine2009 % ('Invalid array value', 'ZTEST')
    x_mean = reduce(operator.add, array) / n
    if len(argList)>2:
        sigma = argList[2]
        return 1 - NORMSDIST([(x_mean-mean) / ( sigma / math.sqrt(n) )])
    return 1 - NORMSDIST([(x_mean-mean) / ( STDEV(array)/ math.sqrt(n) )])