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.
    
    x: X-axis values in a form of numpy.ndarray.
    y: Y-axis values in a form of numpy.ndarray.
    z: Z-axis values in a form of numpy.ndarray.
    d: Step size when plotting rows and columns of array. 
    """
    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)
    
surface(x, y, m44, 2)


# b)
def calculate_statistics(x):
    """
    Do a whole bunch of statistics.
    
    array: Input dataset in a form of list.
    
    return: Floating value of robust deviation.
    """
    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 = m44.flatten()
sigma_r = calculate_statistics(m44)


# 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 - 1721
# Hodnotu sumu najlepsie popisuje robustna odchylka pozadia, ktorej hodnota je 190.
# Je mozne pouzit aj menej vhodnu standardnu odchylku pozadia s hodnotou 197.


# e)
def my_median(x):
    """
    Calculate the median, which is the middle value of sorted dataset.
    
    x: List or array.
    
    return: Floating value of median.
    """
    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(my_median(m44_corr))