from helita.sim import rh15d
import matplotlib.pyplot as plt
import numpy as np
import os
import warnings     # ignore tedious warnings
warnings.filterwarnings("ignore")


##############################################################


def load_rh_data(folder, print_attributes=False):
    # reset IPython kernel
    try:
        # works only if in IPython kernel
        from IPython import get_ipython
        get_ipython().magic("reset -sf")
    except:
        pass


    # defines variables that should be global
    global DATA, WAVELENGTH, WAVELENGTH_INDEX, WAVELENGTH_SELECTED
    #global
    DATA = rh15d.Rh15dout(folder)
    WAVELENGTH = DATA.files[4].wavelength
    #wl_idx = data.files[4].wavelength_indices
    WAVELENGTH_INDEX = np.where((WAVELENGTH.values>279.401921894) &
                                (WAVELENGTH.values<280.501526399))[0]
    WAVELENGTH_SELECTED = WAVELENGTH[WAVELENGTH_INDEX]


    # print attributes from RH code
    if print_attributes==True:
        for i, file in enumerate(DATA.files):
            print("\n\n--------------------------\n\n")
            print("data.files[index]:      ", i)
            print(file)


##############################################################


# get all output folders within current folder
def get_output_folders():
    folders = os.walk(".")
    output_folders = []
    for path in folders:
        if not "wrong" in path[0] and all(x in path[0] for x in ["x","y","dx","dy"]):
            output_folders.append(path[0].replace("./", ""))
    return output_folders


##############################################################


# plots intensity from RH code
def plot_intensity():
    nx = DATA.files[4].nx
    ny = DATA.files[4].ny
    colx,coly = np.random.randint(nx), np.random.randint(ny)

    """for i in np.linspace(0,nx-1,nx).astype("int"):
        for j in np.linspace(0,ny-1,nx).astype("int"):
            I = DATA.files[4].intensity[i, j, WAVELENGTH_INDEX]
            plt.plot(WAVELENGTH_SELECTED, I/I[0], color="black", alpha=0.01)
    """

    mean_I = np.mean(DATA.files[4].intensity[:,:,WAVELENGTH_INDEX], axis=(0,1))

    plt.plot(WAVELENGTH_SELECTED, mean_I, ".-")
    plt.xlabel("Wavelength (nm)")
    plt.ylabel(r"Intensity $(I/I_0)$")
    plt.grid()
    plt.show()


##############################################################


def kmean_model(data=False, k_create_model=False, test_model=False,
                test_profiles=False, wavelength=False):
    #if not data:
    #    assert False, "Did not input data!"
    if I_data.ndim != 2:
        assert False, "Wrong dimension on input data! Must be two-dimensianal."


    # if 'test_model' is not defined, train a model instead
    if not test_model:
        from rh_kmean import create_kmean_from_data
        kmean = create_kmean_from_data(data, k_value=k_create_model)

    # test a trained model
    else:
        from rh_kmean import test_kmean_model
        test_kmean_model(test_model,
                         test_profiles=test_profiles,
                         wavelength=wavelength,
                         k_nearest=False)


##############################################################


def get_intensity_from_all_folders(plot_converge_box=False):
    folders = get_output_folders()

    colx = []               # x coordinate of converged columns
    coly = []               # y coordinate of converged columns
    colx_nc = []            # x columns that did not converge
    coly_nc = []            # y columns that did not converge
    I = []                  # intensities


    #folders.pop(3)      # as the folder is currently empty
    #folders.pop(1)
    #folders.append(".")
    for f in folders:
        print(f)

    # goes through all folders that is
    # output from rh code
    for f_i, folder in enumerate(folders):
        load_rh_data(folder=folder, print_attributes=False)
        nx = DATA.files[4].nx       # number of columns in x
        ny = DATA.files[4].ny       # number of columns in y
        xnum = DATA.files[3].xnum.values    # column index x from model
        ynum = DATA.files[3].ynum.values    # column index y from model


        # making sure that same columns are
        # not added multiple times
        if f_i != 0:
            for i in range(nx):
                for j in range(ny):
                    add_col = True

                    # if column is already added, do not add again
                    for k in range(len(colx)):
                        if xnum[i] == colx[k] and ynum[j] == coly[k]:
                            add_col = False
                            break

                    if add_col and ~np.isnan(DATA.files[4].intensity[i,j,:].values).any():
                        colx.append(xnum[i])
                        coly.append(ynum[j])
                        I.append(DATA.files[4].intensity[i, j, WAVELENGTH_INDEX].values)
                    elif add_col:
                        colx_nc.append(xnum[i])
                        coly_nc.append(ynum[j])


        # if first output folder, add all columns
        else:
            for i in range(nx):
                for j in range(ny):
                    if ~np.isnan(DATA.files[4].intensity[i,j,:].values).any():
                        colx.append(xnum[i])
                        coly.append(ynum[j])
                        I.append(DATA.files[4].intensity[i, j, WAVELENGTH_INDEX].values)
                    else:
                        colx_nc.append(xnum[i])
                        coly_nc.append(ynum[j])


    # remove elemnts from col(x/y)_nc that later
    # has converged (from other folders/runs)
    pop = []
    for i in range(len(colx_nc)):
        for j in range(len(colx)):
            if colx_nc[i] == colx[j] and coly_nc[i] == coly[j]:
                pop.append(i)
                #colx_nc.pop(i)
                #coly_nc.pop(i)
    # must pop aftewards, else remove wrong index
    pop.reverse()
    for i in pop:
        colx_nc.pop(i)
        coly_nc.pop(i)



    # plot columns that have and have not
    # converged in model
    if plot_converge_box:
        plt.close()
        fig,ax = plt.subplots()
        ax.set_title("Converged synthetic spectra of model: %.2f %%"
                    % (100*len(colx)/1024**2))
        ax.plot([0,1024,1024,0,0], [0,0,1024,1024,0], "--", color="black", label="Model box", lw=0.5)
        ax.scatter(colx, coly, marker="s", label="Converged columns", s=20, color="green")
        ax.scatter(colx_nc, coly_nc, marker="s", label="Non-converged columns", s=20, color="red")
        ax.legend()
        ax.grid()
        fig.show()


    # return intensities, converged columns
    # and non-converged columns
    return np.array(I), np.array([colx,coly]), np.array([colx_nc,coly_nc])


##############################################################








get_intensity_from_all_folders(plot_converge_box=True)
#get_output_folders()


#output_folder = "x0-1024_y0-1024_dxdy25-25"
#load_rh_data(folder=output_folder, print_attributes=False)
#plot_intensity()

#I_data = DATA.files[4].intensity[:,:,WAVELENGTH_INDEX].values
#I_data = I_data[~np.isnan(I_data).any(axis=2)]


"""k = 5
kmean_model(data=I_data, k_create_model=k)
kmean_model(test_model="KMEAN_MODELS/1641_%i"%k,
            test_profiles=I_data,
            wavelength=WAVELENGTH_SELECTED)"""







#

from helita.sim import rh15d
import matplotlib.pyplot as plt
import numpy as np
import os
import warnings     # ignore tedious warnings
warnings.filterwarnings("ignore")


##############################################################


def load_rh_data(folder, print_attributes=False):
    # reset IPython kernel
    try:
        # works only if in IPython kernel
        from IPython import get_ipython
        get_ipython().magic("reset -sf")
    except:
        pass


    # defines variables that should be global
    global DATA, WAVELENGTH, WAVELENGTH_INDEX, WAVELENGTH_SELECTED
    #global
    DATA = rh15d.Rh15dout(folder)
    WAVELENGTH = DATA.files[4].wavelength
    #wl_idx = data.files[4].wavelength_indices
    WAVELENGTH_INDEX = np.where((WAVELENGTH.values>279.401921894) &
                                (WAVELENGTH.values<280.501526399))[0]
    WAVELENGTH_SELECTED = WAVELENGTH[WAVELENGTH_INDEX]


    # print attributes from RH code
    if print_attributes==True:
        for i, file in enumerate(DATA.files):
            print("\n\n--------------------------\n\n")
            print("data.files[index]:      ", i)
            print(file)


##############################################################


# get all output folders within current folder
def get_output_folders():
    folders = os.walk(".")
    output_folders = []
    for path in folders:
        if not "wrong" in path[0] and all(x in path[0] for x in ["x","y","dx","dy"]):
            output_folders.append(path[0].replace("./", ""))
    return output_folders


##############################################################


# plots intensity from RH code
def plot_intensity():
    nx = DATA.files[4].nx
    ny = DATA.files[4].ny
    colx,coly = np.random.randint(nx), np.random.randint(ny)

    """for i in np.linspace(0,nx-1,nx).astype("int"):
        for j in np.linspace(0,ny-1,nx).astype("int"):
            I = DATA.files[4].intensity[i, j, WAVELENGTH_INDEX]
            plt.plot(WAVELENGTH_SELECTED, I/I[0], color="black", alpha=0.01)
    """

    mean_I = np.mean(DATA.files[4].intensity[:,:,WAVELENGTH_INDEX], axis=(0,1))

    plt.plot(WAVELENGTH_SELECTED, mean_I, ".-")
    plt.xlabel("Wavelength (nm)")
    plt.ylabel(r"Intensity $(I/I_0)$")
    plt.grid()
    plt.show()


##############################################################


def kmean_model(data=False, k_create_model=False, test_model=False,
                test_profiles=False, wavelength=False):
    #if not data:
    #    assert False, "Did not input data!"
    if I_data.ndim != 2:
        assert False, "Wrong dimension on input data! Must be two-dimensianal."


    # if 'test_model' is not defined, train a model instead
    if not test_model:
        from rh_kmean import create_kmean_from_data
        kmean = create_kmean_from_data(data, k_value=k_create_model)

    # test a trained model
    else:
        from rh_kmean import test_kmean_model
        test_kmean_model(test_model,
                         test_profiles=test_profiles,
                         wavelength=wavelength,
                         k_nearest=False)


##############################################################


def get_intensity_from_all_folders(plot_converge_box=False):
    folders = get_output_folders()

    colx = []               # x coordinate of converged columns
    coly = []               # y coordinate of converged columns
    colx_nc = []            # x columns that did not converge
    coly_nc = []            # y columns that did not converge
    I = []                  # intensities


    #folders.pop(3)      # as the folder is currently empty
    #folders.pop(1)
    #folders.append(".")
    for f in folders:
        print(f)

    # goes through all folders that is
    # output from rh code
    for f_i, folder in enumerate(folders):
        load_rh_data(folder=folder, print_attributes=False)
        nx = DATA.files[4].nx       # number of columns in x
        ny = DATA.files[4].ny       # number of columns in y
        xnum = DATA.files[3].xnum.values    # column index x from model
        ynum = DATA.files[3].ynum.values    # column index y from model


        # making sure that same columns are
        # not added multiple times
        if f_i != 0:
            for i in range(nx):
                for j in range(ny):
                    add_col = True

                    # if column is already added, do not add again
                    for k in range(len(colx)):
                        if xnum[i] == colx[k] and ynum[j] == coly[k]:
                            add_col = False
                            break

                    if add_col and ~np.isnan(DATA.files[4].intensity[i,j,:].values).any():
                        colx.append(xnum[i])
                        coly.append(ynum[j])
                        I.append(DATA.files[4].intensity[i, j, WAVELENGTH_INDEX].values)
                    elif add_col:
                        colx_nc.append(xnum[i])
                        coly_nc.append(ynum[j])


        # if first output folder, add all columns
        else:
            for i in range(nx):
                for j in range(ny):
                    if ~np.isnan(DATA.files[4].intensity[i,j,:].values).any():
                        colx.append(xnum[i])
                        coly.append(ynum[j])
                        I.append(DATA.files[4].intensity[i, j, WAVELENGTH_INDEX].values)
                    else:
                        colx_nc.append(xnum[i])
                        coly_nc.append(ynum[j])


    # remove elemnts from col(x/y)_nc that later
    # has converged (from other folders/runs)
    pop = []
    for i in range(len(colx_nc)):
        for j in range(len(colx)):
            if colx_nc[i] == colx[j] and coly_nc[i] == coly[j]:
                pop.append(i)
                #colx_nc.pop(i)
                #coly_nc.pop(i)
    # must pop aftewards, else remove wrong index
    pop.reverse()
    for i in pop:
        colx_nc.pop(i)
        coly_nc.pop(i)



    # plot columns that have and have not
    # converged in model
    if plot_converge_box:
        plt.close()
        fig,ax = plt.subplots()
        ax.set_title("Converged synthetic spectra of model: %.2f %%"
                    % (100*len(colx)/1024**2))
        ax.plot([0,1024,1024,0,0], [0,0,1024,1024,0], "--", color="black", label="Model box", lw=0.5)
        ax.scatter(colx, coly, marker="s", label="Converged columns", s=20, color="green")
        ax.scatter(colx_nc, coly_nc, marker="s", label="Non-converged columns", s=20, color="red")
        ax.legend()
        ax.grid()
        fig.show()


    # return intensities, converged columns
    # and non-converged columns
    return np.array(I), np.array([colx,coly]), np.array([colx_nc,coly_nc])


##############################################################








get_intensity_from_all_folders(plot_converge_box=True)
#get_output_folders()


#output_folder = "x0-1024_y0-1024_dxdy25-25"
#load_rh_data(folder=output_folder, print_attributes=False)
#plot_intensity()

#I_data = DATA.files[4].intensity[:,:,WAVELENGTH_INDEX].values
#I_data = I_data[~np.isnan(I_data).any(axis=2)]


"""k = 5
kmean_model(data=I_data, k_create_model=k)
kmean_model(test_model="KMEAN_MODELS/1641_%i"%k,
            test_profiles=I_data,
            wavelength=WAVELENGTH_SELECTED)"""







#

# ignore tedious warnings ############################################################## # reset IPython kernel # works only if in IPython kernel # defines variables that should be global #global #wl_idx = data.files[4].wavelength_indices # print attributes from RH code ############################################################## # get all output folders within current folder ############################################################## # plots intensity from RH code for i in np.linspace(0,nx-1,nx).astype("int"):         for j in np.linspace(0,ny-1,nx).astype("int"):             I = DATA.files[4].intensity[i, j, WAVELENGTH_INDEX]             plt.plot(WAVELENGTH_SELECTED, I/I[0], color="black", alpha=0.01) ############################################################## #if not data: #    assert False, "Did not input data!" # if 'test_model' is not defined, train a model instead # test a trained model ############################################################## # x coordinate of converged columns # y coordinate of converged columns # x columns that did not converge # y columns that did not converge # intensities #folders.pop(3)      # as the folder is currently empty #folders.pop(1) #folders.append(".") # goes through all folders that is # output from rh code # number of columns in x # number of columns in y # column index x from model # column index y from model # making sure that same columns are # not added multiple times # if column is already added, do not add again # if first output folder, add all columns # remove elemnts from col(x/y)_nc that later # has converged (from other folders/runs) #colx_nc.pop(i) #coly_nc.pop(i) # must pop aftewards, else remove wrong index # plot columns that have and have not # converged in model # return intensities, converged columns # and non-converged columns ############################################################## #get_output_folders() #output_folder = "x0-1024_y0-1024_dxdy25-25" #load_rh_data(folder=output_folder, print_attributes=False) #plot_intensity() #I_data = DATA.files[4].intensity[:,:,WAVELENGTH_INDEX].values #I_data = I_data[~np.isnan(I_data).any(axis=2)] k = 5 kmean_model(data=I_data, k_create_model=k) kmean_model(test_model="KMEAN_MODELS/1641_%i"%k,             test_profiles=I_data,             wavelength=WAVELENGTH_SELECTED) #