import numpy as np
import matplotlib.pyplot as plt

def my_hist(A, n):
    """
    Create histogram plot.
    
    Args:
        A: Data array.
        n: Number of columns in histogram.
    """
    dx = (max(A) - min(A))/(n - 1)
    x = np.linspace(min(A), max(A), n)
    counts = []
    for i in x:
        counts.append(sum(map(lambda a: abs(a - i) < dx/2, A)))
    plt.bar(x, counts, edgecolor='k', width=dx)

A = np.loadtxt('hd188041.dat', usecols=7)
my_hist(A, 25)
plt.xlabel("m [mag]")
plt.ylabel("counts")
plt.title('hd188041')
nA = len(A)
sA = np.sqrt(sum((A - np.mean(A))**2)/(nA - 1.46))
lambdaA = 0.22*np.sqrt(nA)*((nA + 1)/(nA - 1)*np.mean(
        (A - np.mean(A))**4)/nA/sA**4 - 3)
lambdaA2 = np.sqrt(nA)*(1.084 - 1.608*np.median(abs(A - np.median(A)))/sA)
print('Hodnota Lambda_kurt pre hviezdu HD188041 je {}.'.format(lambdaA))
print('Hodnota Lambda_mad pre hviezdu HD188041 je {}.'.format(lambdaA2))

B = np.loadtxt('hr7575.dat', usecols=7)
plt.figure()
my_hist(B, 25)
plt.xlabel("m [mag]")
plt.ylabel("counts")
plt.title('hr7575')
nB = len(B)
sB = np.sqrt(sum((B - np.mean(B))**2)/(nB - 1.46))
lambdaB = 0.22*np.sqrt(nB)*((nB + 1)/(nB - 1)*np.mean(
        (B - np.mean(B))**4)/sB**4 - 3)
lambdaB2 = np.sqrt(nB)*(1.084 - 1.608*np.median(abs(B - np.median(B)))/sB)
print('Hodnota Lambda_kurt pre hviezdu HR7575 je {}.'.format(lambdaB))
print('Hodnota Lambda_mad pre hviezdu HR7575 je {}.'.format(lambdaB2))

# Podla indexov normality by mala byt hviezda HD188041 nepremenna a HR7575
# premenna, ale podla vzhladu histogramu by to malo byt naopak. Vysoke hodnoty
# indexov pre HR7575 su sposobene odlahlymi bodmi, ktore silne skresluju
# statistiku. V skutocnosti sa jedna o chemicky pekuliarnu hviezdu, ktora je
# premenna vo filtri v (data HD188041), ale nepremenna vo filtri u (data
# HR7575), vid obrazok hd188041.png.

