#!/usr/bin/python2
'''
The following script explains how to compute inequality indices.

First some formulas from Amartya Sen's and James E. Foster's book
"On Economic Inequality" (1973/1997) are implemented:
(1) Entropy of expected information content (2.11 on pg.35, A. Sen)
(2) Theil's measure (redundancy of the expected information content)
    (2.12 on pg.35, A. Sen)
(3) Atkinson family (Chapter A.2.2, pg.128, J. E. Foster)
(4) Generalized entropy class (Chapter A.4.1, pg.140, J. E. Foster)

Then further inequality computation routines are shown:
(5) Redundancy R (R[1] = Theil[1])
(6) Gini-Hoover index for societies devided into two a-fractiles:
    Gini index and Hoover index computed from Theil index
(7) Gini Index,Hoover Index, three versions of the Theil[1] redundancy
    (without any inequality aversion parameters) and the Lorenz curve.
    This routine is wrapped into the class GiniHooverTheil(), where I
    may add some functionality later.

The last chapter allows you to run experiments with the previous function and
the GiniHooverTheil() class. You should delete that section (or at least any
code which runs experiments) in case you want to use this script in your own
application script. Such a script simply would start with two lines:

#!/usr/bin/python
from onOEI import *

The first line is recommended for UNIX like OS environments and won't do any
harm in other environments. "onOEI" may have to be replaced by the name of the
version you are using. At the end of the listing in the "run demo" section,
calls which still are active, may have to be turned into comments.
(8) For your experiments
    
G.Kluge, Munich, 2008-11-17

'''
###############################################################################
##                  Preparation of some tools and constants                  ##
###############################################################################
#
# NON-PROGRAMMERS DON'T NEED TO READ THIS SECTION
#
#------------------------------------------------------------------------------
# required functions
from math import e, exp, fabs, modf, log as ln
from pprint import pprint

#--- constants ----------------------------------------------------------------
# constants
infinity = 9e999
ln2 = ln(2)
binary = 1/ln2
natural = 1 # 1/log(2.71828182846)
decimal = 1/ln(10)

#--- functions ----------------------------------------------------------------
# logarithm with configurable base
def log(x,*logtype):
    try:
        basefactor = logtype[0]
    except:
        basefactor = 1
    return ln(x) * basefactor
'''application example:
print log (16,(binary)) # yields 4.0
print log(16) # yields 2.77258872224
print log (16,(decimal)) # yields 1.20411998266
#'''

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
def error_if_not_in_range01(value):
    if (value <= 0) or (value > 1): # check range
        raise Exception, \
            str(value) + ' is not in [0,1)!'
def error_if_not_1(value):
    if (value != 1):
        # not 100%
        raise Exception, 'Sum ' + str(value) + ' is not 1!'

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
def sum_of_table(x):
    sum = 0.0
    for x_i in x:
        error_if_not_in_range01(x_i)
        sum += x_i
    error_if_not_1(sum)
    return sum



#--- remarks ------------------------------------------------------------------
# Differences to OIE:
# - My n starts from 0: n = n_OIE - 1

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Other remarks:
# - In programming, a more elegant syntax and use of structures offered by
#   the Python language has been sacrified to understandability, as I want to
#   make this script also understandable to readerw who do not know python.


###############################################################################
##           Some formulas from ON ECOOMIC INEQUALITY in Python              ##
###############################################################################

# ===== 2.11 on pg.35 (A. Sen) ================================================
# (1) Entropy of expected information content
# =============================================================================
# (Base of logarithm is configurable. Natural logarithm is default.)
# (x is a single column table.)
def H(x): # natural logarithm is default
    n = len(x)
    entropy = 0.0 # initiate entropy
    sum = 0.0
    for x_i in x: # work on all x[i]
        error_if_not_in_range01(x_i)
        sum += x_i # add up probabilities for test
        group_negentropy = x_i*log(x_i)
        entropy += group_negentropy # add group negentropies to entropy
        # remember: ln(1/x) = -ln(x)
    error_if_not_1(sum)
    return -entropy # turn negentropy into entropy

#===== 2.12 on pg.35 (A. Sen) =================================================
# (2) Theil's measure (redundancy of the expected information content)
#==============================================================================
# (x is a single column table.)
def T(x): # natural logarithm is default
    n = len(x)
    maximum_entropy = log(n)
    actual_entropy = H(x)
    redundancy = maximum_entropy - actual_entropy
    inequality = 1 - exp(-redundancy)
    return redundancy,inequality,'Theil\'s measure T[1]'

#===== Chapter A.2.2, pg.128 (J. E. Foster) ===================================
# (3) Atkinson family
#==============================================================================
def atkinson_family(x,epsilon):
    should_be_1 = sum_of_table(x) # checks table
    n_x = len(x)
    u_x = 1.0/n_x # mean group income (here Foster used u, I use u_x)
    if (epsilon == 0):
        message = 'epsilon = 0'
        product = 1.0
        for x_i in x: # go through list
            product *= x_i/u_x
        equality = product ** (1.0/n_x)
    elif (epsilon <= 1):
        message = 'epsilon = ' + str(epsilon) +  ' <=1 and not 0'
        sum = 0.0
        for x_i in x: # go through list
            sum += (x_i/u_x) ** epsilon
        equality = (sum/n_x) ** (1.0/epsilon)
    else:
        # error: epsilon is >1
        raise Exception, \
            'epsilon = ' \
            + str(epsilon) + ' is out of range!'
    inequality = 1 - equality # from Atkinson family
    redundancy = -ln(equality) # Kluge
    return redundancy,inequality,'Atkinson family with ' + message

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
def atkinson_family_demo(group_income_pairs,epsilons):
    printout = []
    for epsilon in epsilons:
        printout.append(atkinson_family(group_income_pairs,epsilon))
    return printout
#print atkinson_family_demo((0.01,0.09,0.2,0.3,0.4),(0.25,1))

# You will find such "demo" functions also in the calculation functions below.
# You can delete them or keep them. I used them for first tests and just left
# them in the code in order to give examples how to call the calculation
# functions.

#===== Chapter A.4.1, pg.140 (J. E. Foster) ===================================
# (4) Generalized entropy class
#==============================================================================
def generalized_entropy(x,alpha):
    should_be_1 = sum_of_table(x) # checks table
    n_x = len(x) # (here Foster used n, I use n_x)
    u_x = 1.0/n_x # mean group income (here Foster used u, I use u_x)
    if (alpha == 0):
        message = 'Theil I[0](x) is the \"second\" Theil measure'
        sum = 0.0
        for x_i in x: # go through list
            sum += ln(u_x/x_i)
        redundancy = sum/n_x
    elif (alpha == 1):
        message = 'Theil I[1](x) is the  \"first\" Theil measure'
        sum = 0.0
        for x_i in x: # go through list
            sum += (x_i/u_x) * ln(x_i/u_x)
        redundancy = sum/n_x
    else:
        message = 'Theil I[' + str(alpha) \
            + '](x) from generalized entropy class'
        sum = 0.0
        for x_i in x: # go through list
            sum += 1 - ((x_i/u_x) ** alpha)
        redundancy = sum/alpha/(1-alpha)/n_x
    inequality = 1 - exp(-redundancy)
    return redundancy,inequality,message

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
def generalized_entropy_demo(group_incomes,alphas):
    printout = []
    for alpha in alphas:
        printout.append(generalized_entropy(group_incomes,alpha))
    return printout
#print generalized_entropy_demo((0.01,0.09,0.2,0.3,0.4),(0.25,1,2))

###############################################################################
##                    Other inequality computation tools                     ##
###############################################################################

#===== (Kluge) ================================================================
# (5) Redundancy R
#==============================================================================
def R(x,beta):
    should_be_1 = sum_of_table(x) # checks table
    n_x = len(x) # (here Foster used n, I use n_x)
    u_x = 1.0/n_x # mean group income (here Foster used u, I use u_x)
    sum = 0.0
    for x_i in x: # go through list
        sum += (x_i * ln(x_i/u_x))
    if (beta == 1):
        redundancy = sum
        message = 'R[1](x) = T[1](x)'
        inequality = 1 - exp(-redundancy)
    else:
        if (beta > 0):
            inequality = (1 - exp(-sum)) ** beta
            redundancy = -ln(1-inequality)
            message = 'R[' + str(beta) + '](x)'
        else:
            # beta is 0
            inequality = 1
            redundancy = infinity
            message = 'R[0](x)'
    return redundancy,inequality,message

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
def R_demo(group_incomes,betas):
    printout = []
    for beta in betas:
        printout.append(R(group_incomes,beta))
    return printout
#print R_demo((0.01,0.09,0.2,0.3,0.4),(0.25,1,2))

#===== (Kluge) ================================================================
# (6) Gini-Hoover index: Gini index and Hoover index computed from Theil index
#==============================================================================
# For societies which are devided into two a-fractiles, the Gini index and the
# Hoover index are equal. They can be computed from the Theil index.
#------------------------------------------------------------------------------

# ----- required libraries -----
# from math import exp, e, fabs, log as ln   # has been done previously already
from math import sqrt, tanh

# ----- constants required by this function -----
ineqTheilAt1 = 0.647918229 # ghFromTheil(1) = 0.647918229
oneMinusiTA1 = 1-ineqTheilAt1
approxfactor = -4.63 # in ghApproximation are two terms:
    # - tanh(theil/2.0)*ineqTheilAt1 is "squeezed" 0.5 times on the x scale   
    # - (1-exp(approxfactor*theil))*oneMinusiTA1 is "squeezed"-4.63
    #   times on the x scale.

# ----- general functions ----
asinh = lambda x: ln(sqrt(x*x+1.0)+x)
atanh = lambda x: ln((1.0+x)/(1.0-x))/2.0
atanh2 = lambda x: ln((1.0+x)/(1.0-x))

# ----- special functions -----
# Theil index computed from GiniHoover index
theilFromGh = lambda gh: atanh2(gh)*gh # atanh2(gh)*gh
# GiniHoover index approximated from Theil index
ghApproximation = lambda theil: tanh(theil/2.0)*ineqTheilAt1 \
                  +(1-exp(approxfactor*theil))*oneMinusiTA1
# new GiniHoover approximated from Theil and old GiniHoover
ghRecursive = lambda theil,gh: tanh(theil/gh/2)

def ghFromTheil(theil):
    # ----- computation -----
    theil = float(theil)
    if (theil <= 0.0): return 0.0,theil # theil = 0, error = theil
    theil = min(theil,16) # all Theil indices > 16 are almost 1
    gh_0 = ghApproximation(theil) # start value for GiniHoover index
    t_0 = theilFromGh(gh_0) # start value for Theil index
    gh_1 = ghRecursive(theil,gh_0) # first try for GiniHoover index
    for loopCount in xrange(1000): # limit loop count (for safety)
        t_1 = theilFromGh(gh_1) # first and next tries for Theil index
        if fabs(theil-t_1) < 0.000001: break # loop exit
        try: slope = (gh_1-gh_0)/(t_1-t_0) # for interpolation
        except ZeroDivisionError: break # annother loop exit
        delta = ((t_1+t_0)*slope -gh_1 -gh_0)/2.0 # for interpolation
        gh_0,t_0 = gh_1,t_1 # memorize new values
        gh_1 = max(0,theil*slope - delta) # 1st & next tries f. GiniHoover
        #print '2: gh0=%f, t0=%f, gh1=%f, t1=%f' % (gh_0,t_0,gh_1,t_1) # test
    return gh_1,theil-t_1 # return with GiniHoover and Theil deviation

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
def ghFromTheil_demo():
    for theil in xrange(0,50):
        gh,errTheil = ghFromTheil(theil/10.0)
        print 'gh=%f, t=%f, err=%f' % (gh,theil,errTheil) ### test
#ghFromTheil_demo()


#===== (Kluge) ================================================================
# (7) Gini index, Hoover index, three versions of Theil[1] index, Lorenz curve
#==============================================================================
def threeInequalities(x_pairs):

    table                = []
    theiltable           = []
    rownumber            = -1
    societysize          = 0.0
    societyincome        = 0.0
    
    theiltableerror      = 0
    grouptableItoGerror  = 0
    grouptableGtoIerror  = 0
    grouptableSerror     = 0

    for groupdata_i in x_pairs:
        groupsize_i,groupincome_i = groupdata_i[0],groupdata_i[1]
        # get optional group theil indices (if available)                    #T
        try: theiltable.append(groupdata_i[2])                               #T
        except IndexError: theiltableerror = 1                               #T
        # get basic group data
        societysize      += groupsize_i
        societyincome    += groupincome_i
        percapitaincome_i = groupincome_i/float(groupsize_i)
        row_i = (percapitaincome_i, \
                 groupsize_i, \
                 groupincome_i)
        table.append(row_i)
    table.sort() # required for Lorenz curve and Gini only
    theiltable = list(theiltable)

    lorenzcurve            = [(0.0,0.0)]
    lorenzsize_previous    = 0.0
    lorenzincome_previous  = 0.0

    median                 = 0.0
    groupsize_previous     = 0.0
    groupincome_previous   = 0.0
    tablelength            = len(table)
    oddlength              = tablelength%2

    sum_for_gini           = 0.0
    sum_for_hoover         = 0.0
    sum_for_theilItoG      = 0.0
    sum_for_theilGtoI      = 0.0
    sum_for_theilS         = 0.0

    for percapitaincome_i,groupsize_i,groupincome_i in table:
        rownumber += 1

        lorenzsize_i          = lorenzsize_previous + groupsize_i
        lorenzsize_previous   = lorenzsize_i
        lorenzincome_i        = lorenzincome_previous + groupincome_i
        lorenzincome_previous = lorenzincome_i
        lorenzcurve.append((lorenzsize_i,lorenzincome_i))

        if (median == 0):
            if (tablelength > 2):
                if (lorenzsize_i > societysize/2.0):
                    medianinfo = str(rownumber+1) + '/' + str(tablelength)
                    if oddlength: # odd table length
                        median = percapitaincome_i
                    else:         # even table length
                        median = (groupincome_previous + groupincome_i) \
                                 /float(groupsize_previous + groupsize_i)
                        medianinfo = str(rownumber) + '/' + str(tablelength)\
                                     + '+' + medianinfo
                groupsize_previous = groupsize_i
                groupincome_previous = groupincome_i
            else:
                median = societyincome/societysize
                medianinfo = 'irrelevant'

        deviation_i        = groupincome_i/float(societyincome) \
                             - groupsize_i/float(societysize)
        gini_i             = (lorenzincome_i*2.0 - groupincome_i) * groupsize_i
        hoover_i           = fabs(deviation_i)
        theilItoG_i        = ln(percapitaincome_i) * groupsize_i
        theilGtoI_i        = -ln(percapitaincome_i) * groupincome_i
        theilS_i           = ln(percapitaincome_i) * deviation_i

        sum_for_gini      += gini_i
        sum_for_hoover    += hoover_i
        sum_for_theilItoG += theilItoG_i
        sum_for_theilGtoI += theilGtoI_i
        sum_for_theilS    += theilS_i

        try: sum_for_theilItoG  -= theiltable[rownumber][0]                  #T
        except IndexError: grouptableItoGerror = 1                           #T
        try: sum_for_theilGtoI  -= theiltable[rownumber][1]                  #T
        except IndexError: grouptableGtoIerror = 1                           #T
        try: sum_for_theilS     += theiltable[rownumber][2]                  #T
        except IndexError: grouptableSerror = 1                              #T
    
    percapitaincome = societyincome/float(societysize)
    gini = 1 - float(sum_for_gini)/societysize/societyincome
    hoover = sum_for_hoover/2.0

#   Redundancy = maximum entropy    minus    effective entropy
    theilItiG  = ln(percapitaincome)  - sum_for_theilItoG/float(societysize)
    theilGtiI  = -ln(percapitaincome) - sum_for_theilGtoI/float(societyincome)

    if grouptableSerror: theilS = (theilItiG + theilGtiI)/2.0                #T
    else:                theilS = sum_for_theilS/2.0                         #T
    
    results = (societysize, \
               societyincome, \
               percapitaincome, \
               medianinfo, \
               median, \
               gini*100, \
               hoover*100, \
               theilItiG, \
               theilGtiI, \
               theilS)
    format = 'Size of society:  %.3f\n' + \
             'Income of society: %.3f\n' + \
             'Per capita income: %.3f\n' + \
             'Per capita median (row %s): %.3f\n' + \
             'Gini: %.1f%%' + '\n' + \
             'Hoover: %.1f%%' +'\n' + \
             'Theil[1] income2groups: %.3f\n' + \
             'Theil[1] groups2income: %.3f\n' + \
             'Theil[1] symmetrized: %.3f'
    message = format % results
    
    if (theiltableerror==0):                                                 #T
        message = message + '\n(Theil table status: %d%d%d)' \
          % (1-grouptableItoGerror,1-grouptableGtoIerror,1-grouptableSerror) #T

    return results, \
           message, \
           table, \
           lorenzcurve

# Remark: The Theil index (more precise: the Theil redundancy) can be        #T
# computed based on a society dataset. That set is a table containing        #T
# (groupsize,groupincome) for each group (each fractile) in the society.     #T
# That is the standard situation.                                            #T
# However, also datasets with per group data like                            #T
# (groupsize,groupincome,(groupTheilItoG,groupTheilGtoI,groupTheilS)) may    #T
# be available as an option. This code contains elements required to deal    #T
# which such extended per-group data. This can be done, because the Theil    #T
# redundancy is "decomposable". (Similar support for the Hoover index has    #T
# not been implemented yet. That may happen later. The Gini index is not     #T
# decomposable.                                                              #T

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
def threeInequalities_demo(x_pairs):
    results,message,table,lorenzcurve = (threeInequalities(x_pairs))
    print
    print 'Table:'
    pprint(table)
    print
    print 'Lorenz curve:'
    pprint(lorenzcurve)
    print
    print message
#threeInequalities_demo(((40,47.5),(9,27),(1,23),(50,2.5)))   

def groupTheil_demo():    # ---- not well tested yet ---
    # Experimental: Not yet for use.
    somedata = ( \
     ((40,47.5),(9,27),(1,23),(50,2.5)), \
    \
     ((40,47.5,(0.00,0.01,0.02)), \
      (9,27,(0.10,0.11,0.12)), \
      (1,23,(0.20,0.21,0.22)), \
      (50,2.5,(0.30,0.31,0.32))),
    \
     ((40,47.5,(0.00,0.01,0.02)), \
      (9,27,(0.10,0.11,0.12)), \
      (1,23,(0.20,0.21)), \
      (50,2.5,(0.30,0.31,0.32))) )

    for society in somedata:
        results,message,table,lorenzcurve = \
          (threeInequalities(society))
        print message
        print


#===== (Kluge) ================================================================
# (7a) Gini index, Hoover index, three versions of Theil[1] index, Lorenz curve
#==============================================================================
# The function threeInequalities() has been written in a way which does not
# burden non-programmers with programming techniques. The following class
# implements an inequality object. If you read this code just for understanding
# the formulas, you can skip reading the code of class GiniHooverTheil.

class GiniHooverTheil:
    def __init__(self,x_pairs):

        # Instead of gropdata with the structure (groupsize,groupincome) also
        # a-fractile formats can be processed.
        try:
            a,b = x_pairs
            dummy = a*b
            x_pairs = ((a,b),(b,a))
        except:
            pass
        try:
            a = fabs(float(x_pairs))
            s = str(a).lstrip('0')
            for digit in '123456789':
                s = s.replace(digit,'0')
            b = eval('1'+s) - a # complement
            x_pairs = ((a,b),(b,a))
        except:
            pass
        try:
            if (len(x_pairs) == 1):
                a,b = x_pairs[0]
                x_pairs = ((a,b),(b,a))
        except:
            pass

        # compute distribution statistics
        self.table = []
        try:
            results,self.message,table,lorenzcurve = \
                    threeInequalities(x_pairs)
        except:
            results = (-1,-1,-1,'error',-1,-1,-1,-1,-1,-1)
            self.message = 'error'
            table =[]
            lorenzcurve =[0]
        #-------------------------------------------------

        # prepare output data
        self.societysize, \
        self.societyincome, \
        self.percapitaincome, \
        self.medianinfo, \
        self.median, \
        self.gini, \
        self.hoover, \
        self.theilItoG, \
        self.theilGtoI, \
        self.theilS = results

        # table
        lorenzcurve.reverse()
        table.reverse()
        dummy = lorenzcurve.pop() #remove (0.0,0.0)
        while len(table):
            r0,r1,r2 = table.pop()
            l0,l1 = lorenzcurve.pop()
            self.table.append((r0,r1,r2,l0,l1))

        # compute Plato index (that's my name for the GiniTheil index)
        self.plato,errTheil = ghFromTheil(self.theilS)

        # display in format like the one for the "Pareto principle"
        self.aFractile = self.plato/2.0+0.5

        # output
        self.message = \
            self.message.replace('Hoover', \
                'Plato: %.1f%% (%.0f:%.0f)\nHoover' \
                % (self.plato*100,self.aFractile*100,100-self.aFractile*100))

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
def giniHooverTheil_demo(x_pairs):
    # generate inequality object
    incomedistribution = GiniHooverTheil(x_pairs)
    # print some data provided by the object
    print incomedistribution.message
#giniHooverTheil_demo(((40,47.5),(9,27),(1,23),(50,2.5)))

def germany2004_demo():
    # table for net incomes in Germany (2004)
    # column(1): group size
    # column(2): per capita income (converted into group income)
    incomedistribution2004de = \
        ([(groupsize,groupsize*percapitaincome) \
        for groupsize,percapitaincome in (\
        (5043009,653),\
        (1789731,3664),\
        (1658405,6233),\
        (1609017,8726),\
        (1436832,11219),\
        (1371636,13752),\
        (2833254,17510),\
        (3059569,22539),\
        (3099288,27471),\
        (3743779,33517),\
        (3866269,43156),\
        (4971433,70066),\
        (419001,163964),\
        (87869,331804),\
        (21729,670095),\
        (9688,2740910))])
    giniHooverTheil_demo(incomedistribution2004de)
#germany2004_demo()

'''
Size of society:  35020509.000
Income of society: 1052552946587.000
Per capita income: 30055.330
Per capita median (row 7/16+8/16): 20121.070
Gini: 51.2%
Plato: 53.6% (77:23)
Hoover: 36.6%
Theil[1] income2groups: 0.737
Theil[1] groups2income: 0.547
Theil[1] symmetrized: 0.642
'''
def aFractiles_demo():
    # generate inequality object
    paretoPrinciple = GiniHooverTheil(0.8)
    print paretoPrinciple.message
    print
    paretoPrinciple = GiniHooverTheil(0.05)
    print paretoPrinciple.message
    print
    paretoPrinciple = GiniHooverTheil(600)
    print paretoPrinciple.message
    print
    paretoPrinciple = GiniHooverTheil((800,200))
    print paretoPrinciple.message
    print
    paretoPrinciple = GiniHooverTheil([(800,200)])
    print paretoPrinciple.message
    print
#aFractiles_demo()

def aFractiles2_demo():
    afractiles = (
      0.620,  # Theil - Hoover at minimum
      0.731,  # Theil = Hoover
      0.740,  # Theil = 0.5
      0.800,  # "Pareto principle"
      0.824,  # Theil = 1.0
      0.917,  # Theil = 2.0
      0.984)  # Theil = 4.0
    for dataset in afractiles:
      inequalityobject = GiniHooverTheil(dataset)
      print inequalityobject.message
      print
#aFractiles2_demo()

#===== (Kluge) ================================================================
# (8) For your experiments
#==============================================================================
###############################################################################
##                          Some data to play with                           ##
###############################################################################

# Below the playfield starts. You can removr it in case you do not want to
# run experiments in that area. The code in the playfield is not reqyired
# by the calculation functions above this area.

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# example for some propabilities or shares
x = (0.01,0.09,0.2,0.3,0.4)
sum_x = 0.0
for x_i in x:
    sum_x += x_i

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# example for some group incomes
group_incomes = (20000,40000,60000,80000)
sum_group_incomes = 0.0
for x_i in group_incomes:
    sum_group_incomes += x_i

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# example for some pairs of group size shares and group income shares
x_pairs = []
sum_x = 0
for x_i in x:
    x_pairs.append((1/float(len(x)),x_i))
    sum_x += x_i
#pprint(x_pairs)

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
epsilons = (1,0.999,0.5,0.0001,0,-0.0001,-0.5,-0.999,-1,-5,-100)

alphas = epsilons

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
a_fractiles = (
    (0.620,0.380),  # Theil - Hoover at minimum
    (0.731,0.269),  # Theil = Hoover
    (0.740,0.260),  # Theil = 0.5
    (0.800,0.200),  # "Pareto principle"
    (0.824,0.176),  # Theil = 1.0
    (0.917,0.083),  # Theil = 2.0
    (0.984,0.016))  # Theil = 4.0


###############################################################################
##                                  Play                                     ##
###############################################################################
#------------------------------------------------------------------------------
def demo0(x,alpha,epsilon):
    if alpha == 1:
        print T(x)
    if (epsilon <= 1):
        print atkinson_family(x,epsilon)
    print generalized_entropy(x,alpha)
    print R(x,alpha)
#demo0((0.01,0.29,0.3,0.4),1,1)

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
def demo1():
    print '\n=== H(x) demo ==='
    print H(x)

    print '\n=== T(x) demo ==='
    print T(x)

    print '\n=== atkinson_family demo with different epsilons ==='
    for list in atkinson_family_demo(x,epsilons):
        print list

    print '\n=== generalized entropy demo with different alphas ==='
    for list in generalized_entropy_demo(x,alphas):
        print list

    print '\n=== straight redundancy demo with different alphas ==='
    for list in R_demo(x,alphas):
        print list
#demo1()

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
def demo2():
    x = (0.01,0.29,0.3,0.4)
    epsilons = (2,1.5,1,0.5,0,-0.5,-1,-1.5,-2,-4)
    for epsilon in epsilons:
        print
        alpha = epsilon
        demo0(x,alpha,epsilon)
#demo2()

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
#http://www.poorcity.richcity.org/calculator/?quantiles=1,1|1,2|1,4|1,8|1,16
def make_csl0(): # for analysis in spreadsheet
    #exponent = 1 # Gini 0.465
    exponent = 0.5 # Gini 0.264
    x0 = (1,2,4,8,16)
    sum = 0.0
    xtest = []
    for x_i in x0:
        sum += x_i**exponent #
    for x_i in x0:
        xtest.append(x_i**exponent/sum)
    #print T(xtest)
    #print xtest
    print 'beta;Atkinson;ineq(Theil);ineq(R)'
    for beta in xrange(-20,31):
        line = ''
        alpha = beta/10.0
        epsilon = alpha
        if epsilon <= 1:
            atkinson = atkinson_family(xtest,epsilon)[1]
        line = '%f;%f;%f;%f' % \
            (alpha, \
            atkinson, \
            generalized_entropy(xtest,alpha)[1], \
            R(xtest,alpha)[1])
        print line

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
def demoThreeInequalities():
    threeInequalities_demo(((40,47.5),(9,27),(1,23),(50,2.5)))
#demoThreeInequalities()
'''
http://www.poorcity.richcity.org/calculator/?quantiles=50,2.5|40,47.5|9,27|1,23
100 quantile elements, 4 quantiles

Mean:                           1.000
Median:         44.4%           0.556 (#1/4, #2/4)
 
           Inequality         Welfare
1-e^-TheilT:    64.1%           2.786   (1/Welfare)
1-e^-TheilL:    72.7%           0.273
1-e^-TheilS:    68.7%           0.313
Gini:           64.5%           0.355
Plato:          68.8%           0.312   Pareto: 844/156
100%-SOEP:      78.5%           0.215
Hoover:         47.5%           0.525

Theil-T Redundancy:         1.025
Theil-L Redundancy:         1.299
Symmetric Redundancy:       1.162
Inequality Issuization:    +0.687
'''

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
def demoParetoPrinciple(): #(Old. Now there is direct support fot a-fractiles.)
    for a,b in a_fractiles:
        print '----------'
        threeInequalities_demo(((a,b),(b,a)))
    print '----------'
#demoParetoPrinciple()

#----- run demo ---------------------------------------------------------------

##make_csl0()

##demoThreeInequalities()

##demoParetoPrinciple()

germany2004_demo()

##aFractiles_demo()

##groupTheil_demo()

 
# [23475]