import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
from scipy import stats


# a)
m44 = np.loadtxt("m44.dat")
N = len(m44)
n = N**2
# meshgrid vytvori matice x a y s rovnakymi rozmermi ako ma m44
x, y = np.meshgrid(range(N), range(N))

def surface(x, y, z, d):
    """
    Plot 3D surface graph.
    
    Args:
        x: X-ove hodnoty vo forme numpy.ndarray.
        y: X-ove hodnoty vo forme numpy.ndarray.
        z: X-ove hodnoty vo forme numpy.ndarray.
        d: Velkost kroku pri vykresleni 3D grafu. 
    """
    fig = plt.figure(figsize=(15, 15))
    ax = fig.add_subplot(111, projection='3d')
    # rstride a cstride specifikuje jemnost delenia, cmap farebnu mapu
    ax.plot_surface(x, y, z, rstride=d, cstride=d, cmap=cm.coolwarm, vmax=3000)
    ax.view_init(30, -80)

# Realnu hviezdokopu mozme vidiet na obrazku m44.png.
surface(x, y, m44, 2)


# b)
def calculate_statistics(x):
    """
    Vypocita plno roznych statistickych parametrov a vypis ich.
    
    Args:
        x: Vstupne data vo forme listu.
    
    Return: Robustna odchylka.
    """
    n = len(x)
    sigma = np.sqrt(np.sum(np.power(x - np.average(x), 2))/(n - 1))
    gamma1 = np.sum(np.power(x - np.average(x), 3))/n/sigma**3
    gamma2 = ((n + 1)/(n - 1)*np.sum(np.power(x - np.average(x), 4))/n/sigma**4
              - 3)
    print("priemer: ", np.average(x))
    print("median: ", np.median(x))
    print("modus: ", stats.mode(x))
    print("standardna odchylka: ", sigma)
    print("std: ", np.std(x))
    print("sikmost: ", gamma1)
    print("skewness: ", stats.skew(x))
    print("spicatost: ", gamma2)
    print("kurtosis: ", stats.kurtosis(x))
    # robustna odchylka
    sigma_r = 1.482*np.median(np.abs(x - np.median(x)))
    print("robustna odchylka: ", sigma_r, "\n")
    return sigma_r

# vektorizuje array
m44_arr = m44.flatten()
sigma_r = calculate_statistics(m44_arr)


# c)
m44_corr = m44[np.where(abs(m44 - np.median(m44)) < 3.5*sigma_r)]

sigma_r = calculate_statistics(m44_corr)


# d)
# Pozadie (background) je subor vsetkych pixelov bez odlahlych bodov.
# Intenzita pozadia je najvhodnejsie definovana ako median pozadia - 1722.
# Hodnotu sumu najlepsie popisuje robustna odchylka pozadia - 190.
# Je mozne pouzit aj menej vhodnu standardnu odchylku pozadia s hodnotou 197.


# e)
def my_median(x):
    """
    Vypocita median.
    
    Args:
        x: List hodnot.
    
    Return: Hodnota medianu.
    """
    x = np.sort(x)
    n = len(x)
    if n%2 == 1:
        return x[int((n - 1)/2)]
    else:
        return 0.5*(x[int(n/2)] + x[int(n/2) - 1])

print('Hodnota mojho medianu je {}.'.format(my_median(m44_corr)))

# f)
def my_modus(x):
    """
    Najde najcastejsie sa vyskytujucu hodnotu zo suboru.
    
    Args:
        x: List hodnot.
    
    Returns: Hodnota modusu a jeho vyskyt.
    """
    x = np.sort(x)
    N = 1   # najvyssi pocet vyskytov 
    i = 0
    result = x[0]   # aktualny modus
    while i <= len(x):
        n = 1
        while i + n < len(x) and x[i] == x[i + n]:
            n += 1
        # Ak najdeme novy modus s castejsim vyskytom, menime N a result.
        if n > N:
            N = n
            result = x[i]
        i += n
    return result, N

my_mode, N = my_modus(m44_arr)
print('Hodnota mojho modusu je {} a jeho vyskyt {}.'.format(my_mode, N))

# g)
m44_smooth = m44.copy()
background = np.median(m44_arr)
N = len(m44_smooth[0])
for i in range(N):
    for j in range(N):
        if np.abs(m44_smooth[i][j] - background) < 3.5*sigma_r:
            m44_smooth[i][j] = background
surface(x, y, m44_smooth, 1)

# h)
def inst_mag(matrix, background, x, y, r, sigma):
    """
    Vypocita instrumentalnu hviezdnu velkost.
    
    Args:
        matrix: Datova sada vo formate numpy.ndarray.
        background: Hodnota priemerneho sumu.
        x: X-ova suradnica pixelu stredu hviezdy.
        y: Y-ova suradnica pixelu stredu hviezdy.
        r: Stvorcovy 'polomer' okruhu hviezdy, z ktoreho urcujeme vyslednu
           intenzitu hviezdy.
        sigma: Hodnota deviacie sumu.
    
    Returns: Hodnota instrumentalnej hviezdnej velkosti s dolnou a hornou
             chybou.
    """
    n = (2*r + 1)**2   # pocet pixelov vyrezu z okolia hviezdy
    # Od celkovej intenzity odcitame intenzitu sumu.
    I = np.sum(matrix[x - r : x + r + 1, y - r : y + r + 1]) - n*background
    dI = np.sqrt(n)*sigma
    i_mag = -2.5*np.log10(I)
    di_mag1 = -2.5*np.log10(I + dI) + 2.5*np.log10(I)   # vypocet hornej chyby
    di_mag2 = -2.5*np.log10(I) + 2.5*np.log10(I - dI)   # vypocet dolnej chyby
    return i_mag, di_mag1, -di_mag2

# i)
def find_stars(matrix, n):
    """
    Najde n najjasnejsich pixelov spolu s ich poziciou.
    
    Args:
        matrix: Datova sada vo formate numpy.ndarray.
        n: Pocet najjasnejsich pixelov, ktore potrebujeme.
    
    Returns: List n trojic, kde prva hodnota je intenzita, druha x-ova
             a tretia y-ova suradnica.
    """
    sequence = np.sort(matrix.flatten())   # zoradi prvky matice vzostupne
    sequence = sequence[::-1]   # zoradi vektor zostupne
    result = []
    for i in range(n):
        x, y = np.where(matrix == sequence[i])
        result.append([sequence[i], x, y])
    return np.array(result)

# Najdeme pozicie najjasnejsich pixelov a identifikujeme ich so stredmi hviezd.
print(find_stars(m44, 20))
# Pozicia gamma Cnc je priblizne [x, y] = [193, 19].
# Urcime instrumentalne hviezdne velkosti gamma Cnc a m44.
imag_gCnc = inst_mag(m44, 1722, 193, 19, 3, sigma_r)
imag_M44 = inst_mag(m44, 1722, 125, 65, 30, sigma_r)
# Urcime limitnu instrumentalnu hviezdnu velkost pre pomer signal:sum = 5.
imag_lim = -2.5*np.log10(5*sigma_r)
# Okalibrujeme jasnosti porovnanim instrumentalnej a zdanlivej hviezdnej
# velkosti pre gamma Cnc.
mag_corr = 4.65 - imag_gCnc[0]
print('instrumentalna hviezdna velkost gamma Cnc:', imag_gCnc)
print('instrumentalna hviezdna velkost M44:', imag_M44)
print('zdanliva hviezdna velkost M44:', imag_M44[0] + mag_corr)
print('limitna hviezdna velkost:', imag_lim + mag_corr)
