from IPython.display import SVG
SVG(filename='Gini_Coefficient.svg')

def gini_coefficient(data,sort=True):
    "Compute the gini coefficient for some distribution"
    slope = float(sum(data)) / len(data) 
    running_sum = 0
    i = 0
    area = 0
    if sort:
        data = sorted(data)
    for d in data:
        i += 1
        running_sum += d
        area += slope*i - running_sum
    return (2. * area) / (running_sum * len(data))

import numpy as np
gini_coefficient(np.ones(100))

gini_coefficient(xrange(10**5))

ethersale_data = np.load('ethersale.npz')['data']
ethersale_data.sort()

%matplotlib inline
%config InlineBackend.figure_format = 'svg'

import matplotlib.pyplot as plt

def Lorenz_Curve(data):
    "Calculates the Lorenz Curve from a wealth distribution"
    total = sum(data)
    out = [0]
    for d in sorted(data):
        out.append(out[-1] + d / total)
    return out

X = [i / float(len(ethersale_data)) for i in range(len(ethersale_data)+1)]
Y = Lorenz_Curve(ethersale_data)

f = plt.figure()
ax = f.add_subplot(111, aspect='equal')
ax.yaxis.tick_right()
ax.yaxis.set_label_position("right")
ax.yaxis.set_ticks_position('both')
plt.title('Ethereum Genesis Lorenz Curve')
plt.plot(X, X, label='Line of Equality')
plt.plot(X, Y, label='Lorenz Curve')
plt.xlabel('Percentage of Stakeholders')
plt.ylabel('Percentage of Stake')
plt.legend(loc='upper left')
plt.show()

print "Gini Coefficient:", gini_coefficient(ethersale_data)

ginis = [gini_coefficient(ethersale_data)]
data = np.copy(ethersale_data)
sample_frequency = 200000
for t in xrange(1,20000001):
    idx = np.random.randint(len(data))
    data[idx] += 7
    if idx != len(data)-1 and data[idx] > data[idx+1]:
        data[idx], data[idx+1] = (data[idx+1], data[idx])
    if t % sample_frequency == 0:
        ginis.append(gini_coefficient(data,sort=False))

X = sample_frequency * np.arange(len(ginis))
X_in_Years = X * 12. / (60*60*24*365) # Assuming 12 second timesteps

from scipy.optimize import curve_fit 
def model(x, alpha):
    return ginis[0] / (x / alpha + 1)
(alpha,),_ = curve_fit(model,X,ginis)
Y = model(X, alpha)

import pylab
pylab.ylim([0,1])
plt.plot(X_in_Years, ginis, color='blue', label='Simulation')
plt.plot(X_in_Years, Y, color='red', linestyle='--', linewidth=2, label='Model Fit')
plt.title('Egalitarian Mining')
plt.ylabel('Gini Coefficient')
plt.xlabel('Years')
plt.legend(loc='upper right')
plt.show()
print "⍺ ≈", alpha

ginis = [gini_coefficient(ethersale_data)]
data = np.copy(ethersale_data)
csum = np.cumsum(ethersale_data)
sample_frequency = 200000
for t in xrange(1,20000001):
    top_stakes = csum > np.random.randint(csum[-1])
    csum += top_stakes * 7
    idx = np.argwhere(top_stakes)[0,0]
    data[idx] += 7
    if idx != len(data)-1 and data[idx] > data[idx+1]:
        csum[idx], csum[idx+1] = (csum[idx+1] - data[idx], csum[idx] + data[idx+1]) 
        data[idx], data[idx+1] = (data[idx+1], data[idx])
    if t % sample_frequency == 0:
        ginis.append(gini_coefficient(data,sort=False))

X = sample_frequency * np.arange(len(ginis))
X_in_Years = X * 12. / (60*60*24*365)

import pylab
plt.plot(X_in_Years, ginis)
plt.title('Weighted Mining')
plt.ylabel('Gini Coefficient')
plt.xlabel('Years')
plt.show()

def slasher_sim(reward=0.03, sample_frequency=1000, steps=500000):
    ginis = [gini_coefficient(ethersale_data)]
    data = np.copy(ethersale_data)
    csum = np.cumsum(ethersale_data)
    reward = 0.03
    sample_frequency = 1000

    for t in xrange(1,steps+1):
        new_csum = np.copy(csum)
        for _ in xrange(15):
            top_stakes = csum > np.random.randint(csum[-1])
            new_csum += top_stakes * reward
            idx = np.argwhere(top_stakes)[0,0]
            data[idx] += reward
            if idx != len(data)-1 and data[idx] > data[idx+1]:
                new_csum[idx], new_csum[idx+1] = (new_csum[idx+1] - data[idx], new_csum[idx] + data[idx+1]) 
                data[idx], data[idx+1] = (data[idx+1], data[idx])
        csum = new_csum
        if t % sample_frequency == 0:
            ginis.append(gini_coefficient(data,sort=False))

    return ginis

small_ginis = slasher_sim(reward=0.03)
large_ginis = slasher_sim(reward=0.47)

X = 1000 * np.arange(len(small_ginis))
X_in_Days = X * 12 / (60.*60*24)

import pylab
plt.plot(X_in_Days, small_ginis, label="0.02x Genesis Supply")
plt.plot(X_in_Days, large_ginis, label="Regular Issuance")
plt.title('Slasher Ghost')
plt.ylabel('Gini Coefficient')
plt.xlabel('Days')
plt.legend(loc='upper left')
plt.show()