import copy, math, shutil, os, json, itertools, decimal, sklearn, scipy, warnings, re, multiprocessing
import time

from . import utils, ihs, molecule, geometry, ML, inverseMethod, feff, plotting, funcModel, uiControls, sampling, exafs
from .funcModel import FuncModel
import numpy as np
import pandas as pd
import ipywidgets as widgets
from ipywidgets import HBox, VBox
from IPython.display import display

decimal.getcontext().rounding = decimal.ROUND_CEILING
feff_exe = ''


def tryFloat(text):
    try:
        return float(text)
    except:
        return text


def unique(list1):
    unique_list=[]
    for x in list1:
        if x not in unique_list:
            unique_list.append(x)
    return unique_list


def truncate(x, ndigits):
    return float(round(decimal.Decimal(str(x)), ndigits))


def run_FEFF(folder, file_inp):
    feff.runLocal(folder, feffVersion=6)


def create_input(folder, file_xyz, key):
    """
    This function creates the FEFF input file. The key parameter are the FEFF key terms (e.g. RMAX, AFOLP, HOLE ,... etc. see http://leonardo.phys.washington.edu/feff/Docs/feff7/feff7-4.html#pri)
    """
    m = molecule.Molecule(file_xyz)
    feff.generateInputUniqueAbsorber(m, folder, additional='', absorber=0, feffVersion=6, separatePotentials=False, **key)


def find_relevant_paths(lines, constraints):
    """
    read the file.dat generated by feff
    """
    lines1=[]# stands for matrix lines1
    lines2=[]# stands for matrix lines2
    file_names=[]
    r_path=[]# relevant per relevant paths
    for row in lines: #first selection
        num=lines.index(row)
        line=lines[num]
        line_selected= list(map(tryFloat,line.strip().split())) #selected lines
        lines1.append(line_selected)
    
    for i in range(len(lines1)):
        if 'file' in lines1[i]:
            lines2.append(lines1[i+1:len(lines1)])
    lines2=lines2[0]
    
    for j in range(len(lines2)):
        file_names.append(lines2[j][0])
        if lines2[j][2]>=constraints['Amplitude']: #for the definition of the amplitude see eq. A3 in https://doi.org/10.1103/PhysRevB.52.2995
            r_path.append(lines2[j]) #append all the relevant paths
    
    file_names=[]
    amplitude=[]
    degeneracy=[]
    n_legs=[]
    r_eff=[]
    for k in range(len(r_path)):
        file_names.append(r_path[k][0])
        amplitude.append(r_path[k][2])
        degeneracy.append(r_path[k][3])
        n_legs.append(r_path[k][4])
        r_eff.append(r_path[k][5])
    relevant_paths = pd.DataFrame({'file_names':file_names,'amplitude':amplitude,
                             'degeneracy':degeneracy,'n_legs':n_legs,'r_eff':r_eff})
    return relevant_paths


def read_feff_path(lines):
    """
    returns the x,y,z coordinates of a given path in the related FEFFNNNN.dat file (e.g. for the Cu->O1->O2 it will return the components of the Cu->O1 and O1->O2 vectors)
    """
    lines1=[]# stands for matrix lines1 structure
    lines2=[]# stands for matrix lines2 structure
    lines3=[]
    for row in lines:
        num=lines.index(row)
        line=lines[num]
        line_selected= list(map(tryFloat,line.strip().split()))
        lines1.append(line_selected)

    for i in range(len(lines1)):
        if 'x' in lines1[i]:
            lines2.append(lines1[i+1:len(lines1)])
    lines2=lines2[0]

    for j in range(len(lines2)):
        if 'k' in lines2[j]:
            break
        else:
            lines3.append(lines2[j])

    x=[]
    y=[]
    z=[]
    for k in range(len(lines3)):
        x.append(truncate(lines3[k][0],4)) #truncation at the 4th digit (usually FEFF does it)
        y.append(truncate(lines3[k][1],4))
        z.append(truncate(lines3[k][2],4))
    feff_file_atoms=pd.DataFrame({'x':x,'y':y,'z':z})
    return feff_file_atoms


def get_path_name(path_atoms, all_atoms, all_atom_names):
    path = []
    for ip in range(len(path_atoms)):
        atom_ind = np.argmin(np.linalg.norm(all_atoms-path_atoms[ip], axis=1))
        path.append(all_atom_names[atom_ind])
    return '_'.join(path)


def recoverDLR(folder,file):
    """
    DLR stands for degeneracies, legs and reff. This function read a folder path (e.g. 0000, 0001, ... etc.) and the file.dat contained in it and extract the DLR for the user specified FEFFNNNN.dat file (e.g. feff0001.dat)
    """
    f=open(os.path.join(folder,'files.dat'),"r")
    lines=f.readlines()
    f.close()
    lines1=[]
    for row in lines: #first selection
        num=lines.index(row)
        line=lines[num]
        line_selected= list(map(tryFloat,line.strip().split())) #selected lines
        lines1.append(line_selected)
    for i in range(len(lines1)):
        if file in lines1[i]:
            deg=lines1[i][3]
            nlegs=lines1[i][4]
            reff=lines1[i][5]
    dlr=pd.DataFrame({'file':[file],'deg':[deg],'nlegs':[nlegs],'reff':[reff]})
    return dlr


def extract_FEFF_attributes(lines):
    """
    Given the row lines of a FEFFNNNN.dat file, this function returns the FEFF calculations (e.g. phases, backscattering amplitude, ...etc.) For more ingo look http://leonardo.phys.washington.edu/feff/Docs/feff7/feff7-6.html#ss6.8
    """
    lines1=[]
    attributes=[]
    for row in lines:
        num=lines.index(row)
        line=lines[num]
        line_selected= list(map(tryFloat,line.strip().split()))
        lines1.append(line_selected)
    for i in range(len(lines1)):
        if 'k' in lines1[i]:
            attributes.append(lines1[i+1:len(lines1)])
            attributes=attributes[0]
            c0=[]
            c1=[]
            c2=[]
            c3=[]
            c4=[]
            c5=[]
            c6=[]
            for j in range(len(attributes)):
                c0.append(attributes[j][0])
                c1.append(attributes[j][1])
                c2.append(attributes[j][2])
                c3.append(attributes[j][3])
                c4.append(attributes[j][4])
                c5.append(attributes[j][5])
                c6.append(attributes[j][6])
    dattr=pd.DataFrame({'k':c0,'real[2*phc]':c1,'mag[feff]':c2,'phase[feff]':c3,'red_factor':c4,'lambda':c5,'real[p]':c6})
    return dattr


def return_path_names(path):
    """
    It returns the names of the .csv files attributed to the paths (e.g. Ru.0_C.1_O.2.csv->Ru.0_C.1_O.2, indicating that you are considering the multiple scattering Ru0->C1->C2)
    """
    path_names = [os.path.splitext(filename)[0] for filename in sorted(os.listdir(path))]
    return path_names


# functions used in the module find_items

def names_equal(string1, string2):
    """
    Given two string names (e.g. Cu.0_C.1_C.2 and Cu.0_C.2_C.3), it simply compares two path names. Thi command makes the inversion of the name too (e.g. mirror simmetry: if we have Cu.0_C.1_C.2 and Cu.0_C.2_C.1 these two paths are the same for simmetry)
    """
    split_string1=string1.split('_')
    split_string2=string2.split('_')
   
    str_1_revert=copy.copy(split_string1)
    str_1_revert=str_1_revert[1:]
    str_1_revert=str_1_revert[::-1]
    str_1_revert.insert(0,split_string1[0])
    
    return split_string1==split_string2 or str_1_revert==split_string2


def len_equal(string1, string2):
    """
    Given two string names (e.g. Cu.0_C.1_C.2 and Cu.0_C.2_C.3) It compare path lengths
    """
    return len(string1.split('_')) == len(string2.split('_'))

    
def paths_equal(string1, string2):
    """
    It is similar to the function "compare names" in row  362. In addition here the path names are totally disasambled (e.g. Cu.0_C.1_C.2 becomes [[Cu, 0], [C,1],[C,2]])
    """
    split_string1=string1.split('_')
    split_string2=string2.split('_')
        
    str_1=[]
    str_2=[]
    
    for j1 in range(len(split_string1)):
        str_1.append(split_string1[j1].split('.'))
    for j2 in range(len(split_string2)):
        str_2.append(split_string2[j2].split('.'))

    names1=[]
    names2=[]
    
    for k1 in range(len(str_1)):
        names1.append(str_1[k1][0])
    for k2 in range(len(str_2)):
        names2.append(str_2[k2][0])
        
    names1_revert=copy.copy(names1)
    names1_revert=names1_revert[1:]
    names1_revert=names1_revert[::-1]
    names1_revert.insert(0,names1[0])
    
    return names1==names2 or names1_revert==names2
    
    
# functions used for angles calculations

def angle_between(v1, v2):
    """
    Given two arrays (e.g. two single path coordiantes) this function returns their angle
    """
    v1_u = geometry.normalize(v1)
    v2_u = geometry.normalize(v2)
    if np.linalg.norm(v1_u)==0 or np.linalg.norm(v2_u)==0:
        return float('nan')
    else:
        return np.arccos(np.clip(np.dot(v1_u, v2_u), -1.0, 1.0))


def path_angles(path):
    path=copy.copy(path)
    v=path.values
    diff=[x - v[i - 1] for i, x in enumerate(v)][1:]
    diff.append(-1*v[len(v)-1])
    
    for i1 in range(np.shape(diff)[0]):
        for i2 in range(np.shape(diff)[1]):
            if abs(diff[i1][i2])<=10**-9:
                diff[i1][i2]=0
            else:
                continue          
    angles_rad=[]
    for i in range(len(diff)-1):
        angles_rad.append(angle_between(diff[i+1],diff[i]))
    
    angles=np.zeros(len(angles_rad))
    for j in range(len(angles)):
        if math.isnan(angles_rad[j]):
            angles[j]=float('nan')
        else:
            angles[j]=angles_rad[j]*180/np.pi
    return angles


def compare_angles(path1, path2, sAngle=2, secondCheckA=False, limitsA=None):
    """
    sAngle is the number od digit allowed.  The second check is used to determine if the truncation is still fine
    """
    if limitsA is None: limitsA = [0,2]
    angles1=path_angles(path1)
    angles2=path_angles(path2)
    check = 0
    for i in range(len(angles1)):
        if math.isnan(angles1[i]) and math.isnan(angles2[i]):
            check = check+1
        if (not math.isnan(angles1[i])) and not math.isnan(angles2[i]):
            da = abs(truncate(angles1[i]-angles2[i], sAngle))
            if da == 0:  # before it was 4
                check = check+1
            else:
                if secondCheckA:
                    if limitsA[0] <= da <= limitsA[1]:
                        check = check+1
                    else:
                        continue
                else:
                    continue
        if math.isnan(angles1[i]) and not math.isnan(angles2[i]):
            continue
        elif not math.isnan(angles1[i]) and math.isnan(angles2[i]):
            continue
    return check == len(angles1)


def compare_path(folder1,name1,folder2,name2):
    path1=pd.read_csv(os.path.join(folder1,name1)+'.csv')
    path2=pd.read_csv(os.path.join(folder2,name2)+'.csv')
    if names_equal(name1, name2):
        return True
    else:
        if not len_equal(name1, name2):
            return False
        elif len_equal(name1, name2) and not paths_equal(name1, name2):
            return False
        elif len_equal(name1, name2) and paths_equal(name1, name2):
            return compare_angles(path1,path2)


def checkReff(name1, name2, molecule, sReff=2, secondCheckR=False, limitsR=None, debug=False):
    """
    sReff is the truncation over the half path length. The second check is used to determine if the truncation is still fine
    """
    if limitsR is None: limitsR = [0,2]
    l1=name1.split('_')
    l2=name2.split('_')
    path1=molecule.loc[l1]
    path2=molecule.loc[l2]
    v1=path1.values
    v2=path2.values
    diff1=[x - v1[i - 1] for i, x in enumerate(v1)][1:]
    diff2=[x - v2[i - 1] for i, x in enumerate(v2)][1:]
    diff1.append(-1*v1[len(v1)-1])
    diff2.append(-1*v2[len(v2)-1])
    p1_len=np.zeros(len(diff1))
    p2_len=np.zeros(len(diff2))
    for i in range(len(diff1)):
        p1_len[i]=0.5*np.linalg.norm(diff1[i])
    reff1 = sum(p1_len)
    for j in range(len(diff2)):
        p2_len[j]=0.5*np.linalg.norm(diff2[j])
    reff2 = sum(p2_len)
    delta = truncate(abs(reff1-reff2),sReff)  # before it was fixed to 3
    if debug:
        print(f'delta={delta} limitsR={limitsR} reff1={reff1} reff2={reff2} sReff={sReff} p1_len={p1_len} p2_len={p2_len}')
    if delta == 0:
        return True
    else:
        if secondCheckR:
            return limitsR[0] <= delta <= limitsR[1]
        else:
            return False


def revert_name(string1):     
    split_string1=string1.split('_')    
    names1=split_string1       
    names1_revert=names1[1:]
    names1_revert=names1_revert[::-1]
    names1_revert.insert(0,names1[0])
    return ('_').join(names1_revert)
      

def compare_path0(name1, name2, molecule, sAngle=2, sReff=2, secondCheckR=False, limitsR=None, secondCheckA=False, limitsA=None):
    """
    molecule is ALWAYS the primitive (i.e. the unperturbed one)
    This is a very important function and it is used for the comparison
    """
    if limitsA is None: limitsA = [0, 2]
    if limitsR is None: limitsR = [0, 2]
    l1=name1.split('_')
    l2=name2.split('_')
    path1=molecule.loc[l1]
    path1.index=np.arange(len(l1))
    path2=molecule.loc[l2]
    path2.index=np.arange(len(l2))
    
    l2Rev=revert_name(name2).split('_')
    path2Rev=molecule.loc[l2Rev]
    path2Rev.index=np.arange(len(l2Rev))
    
    if names_equal(name1, name2):
        return True
    else:
        if not len_equal(name1, name2):
            return False
        else:
            if not paths_equal(name1, name2):
                return False
            else:
                angles_eq = compare_angles(path1,path2,sAngle,secondCheckA,limitsA)
                # debug = name1 == 'Ti.0_Cl.12' and name2 == 'Ti.0_Cl.3'
                debug = False
                reff_eq = checkReff(name1,name2,molecule,sReff,secondCheckR,limitsR, debug)
                # if debug:
                #     print(angles_eq, reff_eq)
                if (not angles_eq) and not reff_eq:
                    return False
                elif (not angles_eq) and reff_eq:
                    return compare_angles(path1,path2Rev,sAngle,secondCheckA,limitsA)
                elif angles_eq and not reff_eq:
                    return False
                else:  # angles_eq and reff_eq:
                    return True


def calculate_reff(name, molecule):
    """
    internal directories. The molecule is NOT the primitive (unperturbed) but the deformed one
    """
    molecule.index=molecule['name']
    molecule = molecule.drop(['name','number'], axis=1, errors='ignore')
    l=name.split('_')
    path=molecule.loc[l]
    v=path.values
    diff=[x - v[i - 1] for i, x in enumerate(v)][1:]
    diff.append(-1*v[len(v)-1])
    p_len=np.zeros(len(diff))
    for i in range(len(diff)):
        p_len[i]=0.5*np.linalg.norm(diff[i])
    reff=truncate(sum(p_len),4)
    return reff


def getAttributes(folder,step):
    subfolders=sorted([ f.path for f in os.scandir(folder) if f.is_dir() ])
    aamP=[]
    pphS=[]
    llaM=[]
    cchI=[]
    numberFiles=[] #e.g. 
    nameFolders=[] #e.g. 000, 001, ..., etc
    for m in range(len(subfolders)):
        nameFolder=os.path.basename(subfolders[m])
        name_paths=[]
        for file in sorted(os.listdir(subfolders[m])):
            if file=='molecule.csv':
                continue
            else:
                files_sep=file.split('_')
                if 'attr' in files_sep:
                    name_paths.append('_'.join(files_sep[:len(files_sep)-2]))
        ampInPath=[]
        phInPath=[]
        lambInPath=[]
        chiInPath=[]
        for name in name_paths:
            attr0=pd.read_csv(subfolders[m]+os.sep+name+'.csv')
            attr1=pd.read_csv(subfolders[m]+os.sep+name+'_attr_'+'.csv')
            assert len(attr1) == 1
            reff=attr1['reff'][0]
            deg=attr1['deg'][0]
            
            k0=attr0['k'].values
            amplitude0=attr0['mag[feff]'].values*attr0['red_factor'].values
            phase0=attr0['real[2*phc]'].values+attr0['phase[feff]'].values
            lamb0=attr0['lambda'].values
            
            k=np.arange(np.min(k0),np.max(k0)+step,step)
            amplitude=np.interp(k,k0,amplitude0)
            phase=np.interp(k,k0,phase0)
            lamb=np.interp(k,k0,lamb0)
            
            cchi1=np.zeros(len(k))
            cchi2=np.zeros(len(k))
            chi=np.zeros(len(k))
            
            for i in range(len(k)):
                if k[i]==0:
                    continue
                else:
                    cchi1[i]=(deg*amplitude[i])/(k[i]*(reff)**2)
                    cchi2[i]=np.exp(-2*reff/lamb[i])*np.sin(2*k[i]*reff+phase[i])
                    chi[i]=cchi1[i]*cchi2[i]
            #chi[0] = 2*chi[1] - chi[2]
            
            ampInPath.append(amplitude)
            phInPath.append(phase)
            lambInPath.append(lamb)
            chiInPath.append(chi)
    
            
        aamP.append(ampInPath)
        pphS.append(phInPath)
        llaM.append(lambInPath)
        cchI.append(chiInPath) #not summed
        
        numberFiles.append(len(name_paths))
        nameFolders.append(str(nameFolder))
        
    diag=pd.DataFrame({'name':nameFolders,'path contained':numberFiles})
        
    amP0=[]
    phS0=[]
    laM0=[]
    for i1 in range(len(aamP)):
        for i2 in range(len(aamP[i1])):
            amP0.append(aamP[i1][i2])
            phS0.append(pphS[i1][i2])
            laM0.append(llaM[i1][i2])

    lk=len(np.arange(0,20+step,step))
    chI=np.zeros((lk,len(subfolders)))
        
    for j in range(np.shape(chI)[1]):
        chI[:,j]=np.sum(cchI[j],axis=0)
        
    return diag,amP0,phS0,laM0,chI
    

def getParams(fileName):
    with open(fileName, 'r') as f: params0 = json.load(f)
    return [p[0] for p in params0], [p[1] for p in params0]


def isolateReff(folder):
    name_path=os.path.basename(folder)
    folder_paths=sorted([ f.path for f in os.scandir(folder) if f.is_dir() ])
    numberFolders=len(folder_paths)
    v=0
    for folder2 in folder_paths:
        if len(os.listdir(folder2))>3:
            continue
        elif len(os.listdir(folder2))<=3:
            v=v+1
    if v!=numberFolders:
        print('In folder:',name_path,'there are multiple paths')
    elif v==numberFolders:
        Reff=[]
        for m in range(len(folder_paths)):
            for files in sorted(os.listdir(folder_paths[m])):
                files_sep=files.split('_')
                if 'attr' in files_sep:
                    attr1=pd.read_csv(folder_paths[m]+os.sep+files)
                    Reff.append(attr1['reff'])
                else:
                    continue
        return Reff



################ start: Used to separe the paths from the pathGroups ##########################

def searchName(nome, lista):
    ll=len(lista)
    v=0
    for l in lista:
        if names_equal(nome, l)!= True:v+=1
    if v==ll: 
        return False
    else: 
        return True


def serchNameReff(nome,lista,lista_reff,folder):
    molecule=pd.read_csv(os.path.join(folder,'molecule.csv'))
    reff_nome=calculate_reff(nome,molecule) #reff della lista dei paths presenti nella cartella
    diff=np.zeros(len(lista_reff))
    for i in range(len(lista_reff)):diff[i]=abs(reff_nome-lista_reff[i])
    min_diff=np.min(diff)
    position=list(diff).index(min_diff)
    return lista[position]
################ end: Used to separe the paths from the pathGroups ##########################


class Diagnosys:
    def __init__(self, folder):
        self.cwd=os.path.join(os.getcwd(), folder)
        parameters=pd.read_csv(self.cwd+os.sep+'parameters.csv') #parameters used in the IHS grid saved in the "folder" directory
        self.param=parameters.values
        self.paramNames=list(parameters.columns)
    
    def readAttributes(self):
        """
        Widgets for viewing different spectra from sample. Plot amplitudes, phases and path signals
        """
        #folder_paths=sorted([ f.path for f in os.scandir(os.path.join(self.workingFolder)) if f.is_dir() ])
        names_path=sorted([ f.name for f in os.scandir(self.cwd) if f.is_dir() ]) #names of the paths in the folder directory
        
        wd1=widgets.Dropdown(options=names_path,value=names_path[0],description='Paths:',disabled=False)
        cb=widgets.Checkbox(value=False,description='Separe Deformations',disabled=False,indent=False)
        def selectPath(**args): 
            if cb.value==False:
                pathName=wd1.value
                attributes=[]
                for file in sorted(os.listdir(os.path.join(self.cwd,pathName))):
                    attributes.append(file)
                attributes.remove('report.csv')
                wd2=widgets.Dropdown(options=attributes,value=attributes[0],description='Attributes:',disabled=False)
                sl1=widgets.IntSlider(value=0,min=0,max=5,step=1,description='k_power:',disabled=False,continuous_update=False,orientation='horizontal',readout=True,readout_format='d')
                def selectAttr1(**args):
                    fileAttr=wd2.value
                    kpow=sl1.value
                    attr=np.loadtxt(os.path.join(self.cwd+'//'+pathName,fileAttr))
                    length=np.shape(attr)[0]
                    step=20/length
                    k=np.arange(0,20,step)
                    fig,ax=plt.subplots(figsize=(13,8))
                    for i in range(np.shape(attr)[1]):
                        ax.plot(k,k**kpow*attr[:,i])
                    ax.set_title(fileAttr[:-4])
                    ax.set_xlabel('k')
                    ax.set_ylabel(fileAttr[:-4])
                
                
                tab1=VBox(children=[wd2,sl1])
                out1=widgets.interactive_output(selectAttr1,{wd2.description:wd2,sl1.description:sl1})
                display(tab1,out1)
            
            elif cb.value==True:
                pathName=wd1.value
                report=pd.read_csv(os.path.join(self.cwd+'//'+pathName,'report.csv'))
                namesReport=list(report['name'].values)
                filesReport=list(report['path contained'].values)
                attributes=[]
                for file in sorted(os.listdir(os.path.join(self.cwd,pathName))):
                    attributes.append(file)
                attributes.remove('report.csv')
                attributes.remove('chi.txt')
                wd2=widgets.Dropdown(options=attributes,value=attributes[0],description='Attributes:',disabled=False)
                sl1=widgets.IntSlider(value=0,min=0,max=5,step=1,description='k_power:',disabled=False,continuous_update=False,orientation='horizontal',readout=True,readout_format='d')
                wd3=widgets.Dropdown(options=namesReport,value=namesReport[0],description='Folder:',disabled=False)
                def selectAttr2(**args):
                    fileAttr=wd2.value
                    kpow=sl1.value
                    folderName=wd3.value
                    attr=np.loadtxt(os.path.join(self.cwd+'//'+pathName,fileAttr))
                    length=np.shape(attr)[0]
                    step=20/length
                    k=np.arange(0,20,step)
                    attr_sep=[]
                    for i in filesReport:
                        attr_sep.append(attr[:,0:i])
                        attr=np.delete(attr, np.s_[0:i], axis=1)
                    
                     #_=plt.plot(k,attr_sep[folderName])
                    print(np.shape(attr_sep[folderName]))
                    fig,ax=plt.subplots(figsize=(15,8))
                    for j in range(np.shape(attr_sep[folderName])[1]):
                        ax.plot(k,k**kpow*attr_sep[folderName].flatten())
                    
                    
                tab2=VBox(children=[wd2,wd3,sl1])
                out2=widgets.interactive_output(selectAttr2,{wd2.description:wd2,sl1.description:sl1,wd3.description:wd3})
                display(tab2,out2)
        
                        
        tab0=HBox(children=[wd1,cb])
        out0=widgets.interactive_output(selectPath,{wd1.description:wd1,cb.description:cb})
        display(tab0,out0)
        
    def CrossValidation(self,group,regressor):
        rb=widgets.RadioButtons(options=['Run for each path','Run for all paths'],description='Cross Validation',disabled=False)
        def run_CV(**args):
            chose=rb.value
            if chose=='Run for all paths':
                buttonRunCVForAll=widgets.Button(description='Run CV')
                @buttonRunCVForAll.on_click
                def actioning_CV(selfbutton1):
                    scores=[]
                    for name in group:
                        chi=np.loadtxt(self.cwd+os.sep+name+os.sep+'chi.txt')
                        kf = sklearn.model_selection.KFold(n_splits=10, shuffle=True, random_state=0)
                        res = np.zeros(np.transpose(chi).shape)
                        for train_index, test_index in kf.split(self.param):
                            X_train, X_test = self.param[train_index], self.param[test_index]
                            y_train, y_test = np.transpose(chi)[train_index], np.transpose(chi)[test_index]
                            regressor.fit(X_train, y_train)
                            res[test_index] = regressor.predict(X_test)
                        score=sklearn.metrics.r2_score(np.transpose(chi),res)*100
                        scores.append(score)
                    fig,ax=plt.subplots(figsize=(15,8))
                    ax.bar(np.arange(len(group)),scores,label=group)
                    ax.set_xticks(np.arange(len(group)))
                    print(scores)
                    buttonRunCVForAll.disabled=True
                display(buttonRunCVForAll)
            elif chose=='Run for each path':
                print('ciao')
        out=widgets.interactive_output(run_CV,{rb.description:rb})
        display(rb,out)


def rectW(k, kmin, kmax):
    W = np.zeros(len(k))
    W[(kmin<=k) & (k<=kmax)] = 1
    return W


def convert2RSpaceByLarch(k, chi, kmin=3, kmax=11, dk=1, dk2=None, window='kaiser'):
    from larch import Group
    from larch.xafs import xftf
    rdf = Group(k=k, chi=chi)
    xftf(rdf, kmin=kmin, kmax=kmax, dk=dk, dk2=dk2, window=window)
    return rdf.r, rdf.chir_mag


def fourier_transform_my(k, chi, kmin, kmax, A):
    from scipy.interpolate import Rbf
    w = np.ones(k.shape)
    i = k < kmin + A
    w[i] = 0.5 * (1 - np.cos(np.pi * (k[i] - kmin) / A))
    w[k < kmin] = 0
    i = k > kmax - A
    w[i] = 0.5 * (1 + np.cos(np.pi * (k[i] - kmax + A) / A))
    w[k > kmax] = 0
    M = k.size
    delta = (k[-1] - k[0]) / M
    m = np.arange(0, M // 2)
    wm = 2 * np.pi * m / M / delta
    # print('chi: ' + str(chi.size) + '\n')
    ft = delta * np.exp(complex(0, 1) * wm * k[0]) * np.fft.fft(chi * w)[:M // 2]
    rbfi = Rbf(wm, np.abs(ft))  # radial basis function interpolator instance
    freqw = np.linspace(wm[0], wm[-1], 500)
    rbfft = rbfi(freqw)  # interpolated values
    freqw /= 2  # почему-то нужно - чтобы с Андреа и Larch совпадало
    rbfft /= 2  # почему-то нужно - чтобы с Андреа и Larch совпадало
    return freqw, rbfft, wm, ft


def convert2RSpaceMy(k, chi, kmin=3, kmax=11, A=1):
    R, ft, _, _ = fourier_transform_my(k, chi, kmin, kmax, A)
    return R, ft


class ExafsPredictor:
    def __init__(self, folder, subfolder='dataset'):
        self.workingFolder = folder + os.sep + subfolder
        self.pathGroups = [f for f in sorted(os.listdir(self.workingFolder)) if os.path.isdir(self.workingFolder+os.sep+f)]
        params = pd.read_csv(self.workingFolder+os.sep+'parameters.csv')
        self.sample = {}
        self.k = np.arange(0, 20 + 0.05, 0.05)
        for g in self.pathGroups:
            spectra = np.loadtxt(self.workingFolder + os.sep + g + os.sep + 'chi.txt').T
            if len(spectra.shape) == 1: spectra = spectra.reshape(1,-1)
            self.sample[g] = ML.Sample(params=params, spectra=spectra, energy=self.k, spType='chi')
        self.regressor = None
        self.SMALL = 1e-6
        self.ETOK = 0.2625
        self.noParams = params.shape[1] == 1 and params.columns[0] == 'dummy'
        if self.noParams: self.spectrum = {g:self.sample[g].getSpectrum(0, spType='chi', returnIntensityOnly=True) for g in self.pathGroups}

    def fit(self, method='RBF'):
        """
        paramNames are the selected paths present in the directory folder
        """
        if self.noParams: return
        self.regressor = {}
        for g in self.pathGroups:
            regressor = inverseMethod.getMethod(method)
            # ft = np.fft.fft(self.sample[g].spectra.to_numpy(), axis=1)
            # ft = np.hstack((np.real(ft), np.imag(ft)))
            # regressor.fit(self.sample[g].params, ft)
            regressor.fit(self.sample[g].params, self.sample[g].spectra)
            self.regressor[g] = regressor

    def predict(self, geomParams, pathGroupName):
        if isinstance(geomParams, dict):
            p = np.array([geomParams[pn] for pn in self.sample[pathGroupName].paramNames])
        else:
            if isinstance(geomParams, list): geomParams = np.array(geomParams)
            assert len(geomParams) == len(self.sample[pathGroupName].paramNames)
            p = geomParams
        if len(p.shape) == 1:  p = p.reshape(1, -1)
        chi = self.regressor[pathGroupName].predict(p).reshape(-1)
        # chi = chi[:len(chi)//2] + 1j*chi[len(chi)//2:]
        # res = np.real(np.fft.ifft(chi))
        # assert len(res) == len(self.k)
        # return res
        return chi

    def calcCV(self, method='RBF', CVcount=2):
        res = {}
        for g in self.pathGroups:
            regressor = inverseMethod.getMethod(method)
            sample = samplePreprocessor(self.sample[g])
            res[g],_,_ = ML.crossValidation(regressor, sample.params, sample.spectra, CVcount, nonUniformSample=True)
        return res

    def calcSingleSpectrum(self, geomParams, exafsParams, pathGroupName):
        """
        :param geomParams: (dict or array) are those, listed in dataset/parameters.csv file
        :param exafsParams: dict with keys: enot, S02, sigma2
        """
        if self.noParams:
            chi0 = self.spectrum[pathGroupName]
        else:
            assert self.regressor is not None, 'Call fit method first'
            chi0 = self.predict(geomParams, pathGroupName)
        # generate the shift in the k-space (enot parameter effect)
        en = self.k * self.k - exafsParams['enot'] * self.ETOK
        if min(abs(en)) < self.SMALL:
            en[abs(en) < 2 * self.SMALL] = self.SMALL
        q = np.sign(en) * np.sqrt(abs(en))  # new shifted k-values

        chi1 = exafsParams['S02'] * chi0 * np.exp(-2 * exafsParams['sigma2'] * self.k ** 2)  # chi signal corresponding to the selected path
        chi2 = np.interp(q, self.k, chi1)  # chi signal interpolated over the q grid corresponding to the selected path
        chi2[0] = 2 * chi2[1] - chi2[2]
        return chi2

    def calcTotalSpectrum(self, geomParams, exafsParams):
        """
        geomParams are common for all path groups, but exafsParams is a dict of dicts {pathGroupName:exafsParamsOfGroup}
        """
        matrix_chi = np.zeros((len(self.k), len(self.pathGroups)))
        for i in range(len(self.pathGroups)):
            pathGroupName = self.pathGroups[i]
            matrix_chi[:, i] = self.calcSingleSpectrum(geomParams, exafsParams[pathGroupName], pathGroupName)
        total = matrix_chi.sum(axis=1)
        return total

    def buildExafsParams(self, enot, S02, sigma2):
        """
        Build dict of dicts for the given scalars (common params) or dicts. enot and S02 must be common
        """
        p = {'enot': enot}
        if not isinstance(S02, dict): S02 = {g: S02 for g in self.pathGroups}
        if not isinstance(sigma2, dict): sigma2 = {g:sigma2 for g in self.pathGroups}
        exafsParams = {g:copy.deepcopy(p) for g in self.pathGroups}
        for g in self.pathGroups:
            exafsParams[g]['S02'] = S02[g]
            exafsParams[g]['sigma2'] = sigma2[g]
        return exafsParams

    def getInverseMethodEstimator(self, name, expSpectrum, geomParamBounds, kPower=2, exafsParamsFunc=None, exafsParamsFuncArgBounds=None, multipleS02=False, dependentParams=None, RSpaceParams=None, fitSpectrumInterval=None, plotSpectrumInterval=None):
        """
        Returns estimator class to use in inverse method
        :param expSpectrum: experiment spectrum
        :param exafsParamsFunc: func(arg1=val1, arg2=val2, ...) returns input args (dict) for buildExafsParams
        :param exafsParamsFuncArgBounds: dict{arg1:[min,max], arg2:[min,max],...}
        :param multipleS02: True means make S02 different for different path groups
        :param dependentParams: dict{name: func(exafsParamDict)} functions for calculating dependent params
        :param RSpaceParams: None - means fit in K-space, dict of convert2RSpace args (even empty) means fit in R-space
        """
        if dependentParams is None: dependentParams = {}
        if exafsParamsFuncArgBounds is None:
            assert exafsParamsFunc is None
            exafsParamsFuncArgBounds = {'enot': [-8,8]}
            if multipleS02:
                for g in self.pathGroups:
                    exafsParamsFuncArgBounds[f'S02_{g}'] = [0.5, 1.2]
            else: exafsParamsFuncArgBounds['S02'] = [0.5, 1.2]
            for g in self.pathGroups:
                exafsParamsFuncArgBounds[f'sigma2_{g}'] = [5e-4, 2e-2]
        else: exafsParamsFuncArgBounds = copy.deepcopy(exafsParamsFuncArgBounds)
        assert set(dependentParams.keys()) <= set(exafsParamsFuncArgBounds.keys())
        for dep in dependentParams:
            if dep in exafsParamsFuncArgBounds: del exafsParamsFuncArgBounds[dep]

        def exafsParamsDefaultFunc(**args):
            sigma2 = {}
            for g in self.pathGroups:
                assert f'sigma2_{g}' in args, f'Add sigma2_{g} to the exafsParamsFuncArgBounds of getInverseMethodEstimator'
                sigma2[g] = args[f'sigma2_{g}']
            assert 'enot' in args, 'Add enot to the exafsParamsFuncArgBounds of getInverseMethodEstimator'
            if multipleS02:
                S02 = {}
                for g in self.pathGroups:
                    assert f'S02_{g}' in args, f'Add S02_{g} to the exafsParamsFuncArgBounds of getInverseMethodEstimator'
                    S02[g] = args[f'S02_{g}']
            else:
                assert 'S02' in args, 'Add S02 to the exafsParamsFuncArgBounds of getInverseMethodEstimator'
                S02 = args['S02']
            return {'enot':args['enot'], 'S02':S02, 'sigma2':sigma2}
        if exafsParamsFunc is None:
            exafsParamsFunc = exafsParamsDefaultFunc

        res = copy.deepcopy(self)
        res.name = name
        expSpectrum = copy.deepcopy(expSpectrum)
        expSpectrum.intensity *= expSpectrum.energy**kPower
        if RSpaceParams is not None:
            R, ft = exafs.convert2RSpace(expSpectrum.energy, expSpectrum.intensity, **RSpaceParams)
            expSpectrum = utils.Spectrum(R, ft)
        res.diffFrom = None
        exafsParamsFuncArgNames = sorted(list(exafsParamsFuncArgBounds.keys()))
        g = self.pathGroups[0]
        res.paramRanges = copy.deepcopy(exafsParamsFuncArgBounds)
        if self.noParams:
            res.paramNames = exafsParamsFuncArgNames
        else:
            res.paramNames = list(self.sample[g].paramNames) + exafsParamsFuncArgNames
            for p in geomParamBounds: res.paramRanges[p] = geomParamBounds[p]
        ab = [expSpectrum.energy[0], expSpectrum.energy[-1]]
        res.fitSpectrumInterval = fitSpectrumInterval if fitSpectrumInterval is not None else ab
        res.plotSpectrumInterval = plotSpectrumInterval if plotSpectrumInterval is not None else ab

        def predictSpectrumFunc(arg):
            assert len(arg.shape) == 2
            assert arg.shape[0] == 1
            arg = arg[0]
            ng = 0 if self.noParams else len(self.sample[g].paramNames)
            ne = len(exafsParamsFuncArgBounds)
            if not self.noParams:
                assert len(arg) == ng+ne, f"{len(arg)} != {ng}+{ne} = " + str(self.sample[g].paramNames) + ' + ' + str(list(exafsParamsFuncArgBounds.keys()))
                geomParams = arg[:ng]
            else: geomParams = None
            exafsParamsFuncArg = {exafsParamsFuncArgNames[i]:arg[ng+i] for i in range(ne)}
            for name in dependentParams: exafsParamsFuncArg[name] = dependentParams[name](exafsParamsFuncArg)
            exafsParams = self.buildExafsParams(**exafsParamsFunc(**exafsParamsFuncArg))
            totSp = self.calcTotalSpectrum(geomParams, exafsParams) * self.k**kPower
            if RSpaceParams is not None:
                R, ft = exafs.convert2RSpace(self.k, totSp, **RSpaceParams)
                return np.interp(expSpectrum.energy, R, ft)
            else:
                return np.interp(expSpectrum.energy, self.k, totSp)
        res.predict = predictSpectrumFunc
        res.expSpectrum = expSpectrum
        return res

    def getFuncModel(self, expSpectrum, geomParamBounds=None, kPower=2, exafsParamsFunc=None, exafsParamsFuncArgBounds=None, multipleS02=False, dependentParams=None, RSpaceParams=None, additionalFunc=None, distToExperiment=None, moleculeConstructor=None, **kwargs):
        """
        Returns FuncModel for EXAFS fitting, paramProperties, fixedParams(sample, expSpectrum), defaultParams(smooth and ML method) to use in FuncModelSliders
        :param expSpectrum: experiment spectrum
        :param exafsParamsFunc: func(arg1=val1, arg2=val2, ...) returns input args (dict) for buildExafsParams
        :param exafsParamsFuncArgBounds: dict{arg1:[min,max], arg2:[min,max],...}
        :param multipleS02: True means make S02 different for different path groups
        :param dependentParams: dict{name: func(exafsParamDict)} functions for calculating dependent params
        :param RSpaceParams: None - means fit in K-space, dict of convert2RSpace args (even empty) means fit in R-space
        :param distToExperiment: func(theorySp, expSp, params) -> distance to experiment
        :returns: funcModel
        """
        if not self.noParams:
            assert geomParamBounds is not None
        estimator = self.getInverseMethodEstimator(name='', expSpectrum=expSpectrum, geomParamBounds=geomParamBounds, kPower=kPower, exafsParamsFunc=exafsParamsFunc, exafsParamsFuncArgBounds=exafsParamsFuncArgBounds, multipleS02=multipleS02, dependentParams=dependentParams, RSpaceParams=RSpaceParams)

        def getComplexPart(data, part):
            if part == 'abs': return np.abs(data)
            elif part == 'real': return np.real(data)
            elif part == 'imag': return np.imag(data)
            elif part == 'arg': return np.angle(data)
            else: assert False, f'Unknown complex part {part}'

        def exafsPredictorFunc(thisFuncModel, params):
            expEnergy = estimator.expSpectrum.energy
            arg = np.array([[params[pn] for pn in estimator.paramNames]])
            theory = estimator.predict(arg)
            assert len(theory) == len(expEnergy)
            interval = params['interval']
            complexPart = params['complex part']
            thisFuncModel.data['exp'] = FuncModel.createDataItem('plot', x=expEnergy, y=getComplexPart(estimator.expSpectrum.intensity, complexPart), order=-1, color='black', lw=2, save=False)
            thisFuncModel.data['theory'] = FuncModel.createDataItem('plot', x=expEnergy, y=getComplexPart(theory, complexPart), order=1, save=False)
            thisFuncModel.data['exp complex'] = FuncModel.createDataItem('plot', x=expEnergy, y=estimator.expSpectrum.intensity, plot=False)
            thisFuncModel.data['theory complex'] = FuncModel.createDataItem('plot', x=expEnergy, y=theory, plot=False)
            def xlim(ax):
                ax.set_xlim(interval)
                plotting.updateYLim(ax)
            thisFuncModel.data['xlim'] = FuncModel.createDataItem('custom', plotter=xlim, save=False, order=1000)
            if distToExperiment is None:
                error = utils.rFactorSp(utils.Spectrum(expEnergy,theory), estimator.expSpectrum, p=1, interval=interval)
            else:
                error = distToExperiment(utils.Spectrum(expEnergy,theory), estimator.expSpectrum, params)
            thisFuncModel.data['error'] = FuncModel.createDataItem('text', str='err = %.3g' % error, order=1001)
            if moleculeConstructor is not None:
                geomParamsDict = {pn: params[pn] for pn in geomParamBounds}
                thisFuncModel.data['molecule'] = FuncModel.createDataItem('lazy', generator=lambda: FuncModel.createDataItem('text', str=moleculeConstructor(geomParamsDict).export_xyz_string(), filePostfix='[label].xyz', plot=False))
            if additionalFunc is not None:
                val = additionalFunc(thisFuncModel, params)
                if val is not None: return val
            return error

        paramProperties = funcModel.ParamProperties()
        for pn in estimator.paramNames: paramProperties[pn] = {'type': 'float', 'domain': estimator.paramRanges[pn]}
        paramProperties['interval'] = {'type': 'range', 'domain': estimator.plotSpectrumInterval, 'default':[0.5,4]}
        paramProperties['complex part'] = {'type': 'list', 'domain': ['abs','real','imag','arg'], 'default': 'abs'}
        model = funcModel.FuncModel(function=exafsPredictorFunc, paramProperties=paramProperties, **kwargs)
        return model


def read_xyz_file_with_name_index(file):
    """
    it reads the .xyz file name and it returns the related xyz dataframe as .csv
    """
    m = molecule.Molecule(file)
    df = pd.DataFrame({'x': m.atom[:,0], 'y': m.atom[:,1], 'z': m.atom[:,2]})
    names = [f'{m.atomName[i]}.{i}' for i in range(m.atom.shape[0])]
    df.index = names
    return df


def samplePreprocessor(sample):
    def spectrumPreprocessor(spectrum):
        spectrum.y *= spectrum.x ** 2
        R, spR = exafs.convert2RSpace(spectrum.x, spectrum.y)
        return utils.Spectrum(R, np.abs(spR)).limit([0.5, 4])

    if isinstance(sample, ML.Sample):
        spectra = [spectrumPreprocessor(sample.getSpectrum(i)) for i in range(sample.getLength())]
        return ML.Sample(params=sample.params, spectra=spectra)
    else:
        return spectrumPreprocessor(sample)


class MSPathGroupManager:
    def __init__(self, workingFolder, constraints=None, feffCards=None, paramRanges=None, moleculeConstructor=None, molecule=None, sampleSize=None, maxAdaptiveSize=500, debug=False):
        """
        :param workingFolder:
        :param constraints: dict {'Amplitude': ..., MaxPathLenForShells:..., AtomNamesForShells:[...]} - constrains for path filtering
        :param feffCards: default - {'RMAX': '5.0', 'NLEG': '4', 'CRITERIA': '0 0'} if there are no constraints
        :param paramRanges: dict{name:[a,b],...}
        :param moleculeConstructor: function(paramDict) -> molecule
        :param molecule: Molecule instance or xyz-file name for managing paths for one molecule (when paramRanges and moleculeConstructor is None)
        :param sampleSize: if None run adaptive sampling else - IHS
        :param debug: print debug info
        """
        assert molecule is None or (paramRanges is None and moleculeConstructor is None)
        assert molecule is not None or (paramRanges is not None and moleculeConstructor is not None)
        self.managmentHistory = []
        if constraints is None: constraints = {}
        self.constraints = constraints
        self.debug = debug
        if feffCards is None: feffCards = {'RMAX': '5.0', 'NLEG': '4', 'CRITERIA': '0 0'}
        if 'MaxPathLenForShells' in constraints and constraints['MaxPathLenForShells'] == 1:
            feffCards['NLEG'] = 3 # for Fe->C->Fe (back) (NLEG - atomCount)
        if 'Amplitude' in constraints:
            feffCards['CRITERIA'] = f'{constraints["Amplitude"]} {constraints["Amplitude"]}'
        self.workingFolder = workingFolder
        if os.path.exists(self.workingFolder): shutil.rmtree(self.workingFolder)
        os.makedirs(self.workingFolder, exist_ok=True)
        self.feffCards = feffCards
        self.paramRanges = paramRanges
        self.paramNames = sorted(list(paramRanges.keys())) if paramRanges is not None else None
        self.moleculeConstructor = moleculeConstructor
        self.ir_names = None
        self.primitivePath = self.workingFolder + os.sep + 'primitive.xyz'
        if molecule is None:
            primitive = moleculeConstructor({p:np.mean(paramRanges[p]) for p in self.paramNames})
            primitive.export_xyz(self.primitivePath)
        else:
            if isinstance(molecule, str):
                shutil.copyfile(molecule, self.primitivePath)
            else:
                molecule.export_xyz(self.primitivePath)
        self.primitive = read_xyz_file_with_name_index(self.primitivePath)
        self.datasetFolder = None
        self.create_extreme_cases()
        self.run_FEFF_for_extreme_deformations(self.feffCards)
        self.filterExtremePaths(constraints)
        # secondCheckR and secondCheckA enable check equality of distances (R) and angles (A) by intervals limitsR and limitsA
        # so the number of pathGroups becomes less.
        self.defaultPathGroupParams = {'secondCheckR':True, 'secondCheckA':True, 'limitsR':[0, 0.1], 'limitsA':[0, 0.1], 'sAngle':2, 'sReff':2}
        self.gatherIrreduciblePaths(**self.defaultPathGroupParams)
        if debug: print('Irreducible paths:', self.ir_names)
        self.shells = None
        self.constructShells()
        if debug: print('Auto constructed shells:', self.shells)
        if self.paramRanges is None:
            self.generateDefForOneMolecule()
        else:
            sampleFolder = self.workingFolder+os.sep+'sample_calc'
            if sampleSize is None:
                self.generateDeformationsAdaptive(folder=sampleFolder, maxSampleSize=maxAdaptiveSize)  #, maxError=0.01)
            else:
                sampling.generateInputFiles(self.paramRanges, self.moleculeConstructor, sampleCount=sampleSize, method='IHS', spectralProgram='feff6', spectrCalcParams=self.feffCards, folder=sampleFolder)
                sampling.calcSpectra(spectralProgram='feff6', runType='local', calcSampleInParallel=multiprocessing.cpu_count(), folder=sampleFolder)

    def create_extreme_cases(self, gridDimSize=2):
        """
        starting from the geometryParamRanges is generates a simple grid involving the maxima dn minima deformations plus the null deformation of the initial molecule
        """
        set_folder = self.workingFolder + os.sep + 'extreme'
        if self.paramNames is None:
            # only one molecule is given
            os.makedirs(set_folder + os.sep + 'extreme_0', exist_ok=True)
            shutil.copyfile(self.primitivePath, set_folder + os.sep + 'extreme_0' + os.sep + 'molecule.xyz')
            return
        list_values=[]
        for p in self.paramNames:
            list_values.append(np.linspace(*self.paramRanges[p], gridDimSize))
        matrix_combinations = np.vstack(np.meshgrid(*list_values)).reshape(len(list_values),-1).T
        if gridDimSize == 2:
            center = np.array([[np.mean(self.paramRanges[p]) for p in self.paramNames]])
            matrix_combinations = np.append(matrix_combinations, center, axis=0)

        # Creates folders and put inside of each of them the molecule .xyz file
        dictionary=copy.copy(self.paramRanges)
        # Creates folders and put inside of each of them the molecule .xyz file
        for h in range(np.shape(matrix_combinations)[0]):
            os.makedirs(set_folder+os.sep+'extreme_'+str(h),exist_ok=True)
            for m in range(len(self.paramNames)):
                    dictionary[self.paramNames[m]]=matrix_combinations[h][m]
            molecule = self.moleculeConstructor(dictionary)
            molecule.export_xyz(set_folder+os.sep+'extreme_'+str(h)+os.sep+'molecule.xyz')
            
    def run_FEFF_for_extreme_deformations(self, key):
        """
        Run FEFF for the extreme deformations generated through the function "create_extreme_cases".
        """
        set_folder= self.workingFolder + os.sep + 'extreme'
        folder_paths=sorted([ f.path for f in os.scandir(set_folder) if f.is_dir() ]) #path of subfolders in folder
        for i in range(len(folder_paths)):
            ii = folder_paths[i].split('_')[-1]
            if ii == 'paths': continue
            file = utils.findFile(folder_paths[i], postfix='.xyz')
            create_input(folder_paths[i], file, key)
        for i in range(len(folder_paths)):
            ii = folder_paths[i].split('_')[-1]
            if ii == 'paths': continue
            file = ii+'.inp'
            run_FEFF(folder_paths[i], file)
                    
    def filterExtremePaths(self, constraints):
        """
        Among all the generated paths using the extreme deformation grid, if finds those paths satisfying some user defined constraints, e.g. AMPLITUDE higher than 15%
        Inside the folder extreme deformations, a new folder called selected paths will be generated and it will contain the csv files with the most intense paths for every extreme deformation.
        """
        set_folder= self.workingFolder + os.sep + 'extreme'
        folder_paths = sorted([ f.path for f in os.scandir(set_folder) if f.name != 'extreme_paths' ])
        all_extreme_paths = {}
        for m in range(len(folder_paths)):
            folder_ind = int(folder_paths[m].split('_')[-1])
            f=open(folder_paths[m]+os.sep+'files.dat',"r")
            feff_file_lines=f.readlines()
            f.close()
            data_paths = find_relevant_paths(feff_file_lines, constraints) #file_names (e.g. feff0001.dat) - amplitude - deg.- n_legs -r_eff
            all_atoms_info = read_xyz_file_with_name_index(folder_paths[m]+os.sep+'molecule.xyz') #Returns the xyz input file: x - y - z - name (e.g. Ru0) number (e.g. 0)
            all_atoms = all_atoms_info.loc[:,['x','y','z']].to_numpy()
            all_atom_names = all_atoms_info.index
            extreme_paths = []
            for i in range(len(data_paths['file_names'])):
                with open(os.path.join(folder_paths[m],data_paths['file_names'][i]),"r") as feff_file:
                    lines_feff = feff_file.readlines()
                path_atoms = read_feff_path(lines_feff).to_numpy()
                extreme_paths.append(get_path_name(path_atoms, all_atoms, all_atom_names))
            all_extreme_paths[folder_ind] = extreme_paths
        utils.saveData(all_extreme_paths, set_folder+os.sep+'extreme_paths.json')
    
    def gatherIrreduciblePaths(self, sAngle=2, sReff=2, secondCheckR=False, limitsR=None, secondCheckA=False, limitsA=None):
        """
        IMPORTANT: This function compares the most intense paths generated in the extreme deformations grid and it will provide a list of paths that will be the most intense and not equivalent overall the deformations.
        Two paths are equivalent if evaluated on a reference molecule (e.g. the initial, unperturbed) they have the SAME Reff (half path length), the SAME paths angles and the SAME elements in order or in reverse order.
        """
        if limitsR is None: limitsR = [0, 2]
        if limitsA is None: limitsA = [0, 2]
        extreme_paths = utils.loadData(self.workingFolder + os.sep + 'extreme' + os.sep + 'extreme_paths.json')
        all_paths = set()
        for i in extreme_paths: all_paths = all_paths.union(set(extreme_paths[i]))
        ir_names = []
        for name1 in all_paths:
            found_equal = False
            for name0 in ir_names:
                if compare_path0(name0, name1, self.primitive, sAngle=sAngle, sReff=sReff, secondCheckR=secondCheckR, limitsR=limitsR, secondCheckA=secondCheckA, limitsA=limitsA):
                    found_equal = True
                    break
            if not found_equal: ir_names.append(name1)
            # else: print('For path', name1, 'we found equal', name0)
        ir_names = sorted(ir_names, key=lambda s: len(s.split('_')))
        self.ir_names = ir_names
        return ir_names

    def constructShells(self):
        shells = {}
        for ir in self.ir_names:
            if 'MaxPathLenForShells' in self.constraints and len(ir.split('_'))-1 > self.constraints['MaxPathLenForShells']:
                continue
            name = copy.deepcopy(ir)
            name = name.replace('.','')
            name = re.sub(r"\d", "", name)
            if 'AtomNamesForShells' in self.constraints:
                if len(set(name.split('_')) - set(self.constraints['AtomNamesForShells'])) > 0: continue
            if name in shells: shells[name]['paths'].append(ir)
            else: shells[name] = {'paths':[ir]}
        shells = [{'name':s, **shells[s]} for s in shells]
        self.shells = sorted(shells, key=lambda s: len(s['name'].split('_')))
        return self.shells
    
    def plot_items(self): #Warning, old don't use it
        fig, axs = plt.subplots(len(self.intense_paths), 1,figsize=(10,30))
        for k in range(len(self.intense_paths)):
            amplitudes=np.zeros(len(self.savedItems))
            for i in range(len(self.savedItems)):
                item_split=self.savedItems[i].split('_')
                for j in range(len(self.intense_paths[k]['paths'].values)):
                    names_path=self.intense_paths[k]['paths'].values[j]
                    names_path_split=names_path.split('_')
                    if item_split==names_path_split:
                        amplitudes[i]=self.intense_paths[k]['amplitude'].values[j]
            axs[k].bar(np.arange(len(self.savedItems)),amplitudes,label=self.extremeNames[k])
            axs[k].set_xticks(np.arange(len(self.savedItems)))
            axs[k].legend()
            if k==len(self.intense_paths)-1:
                axs[k].set_xticklabels(self.savedItems,rotation=90)

    def generateDefForOneMolecule(self, key=None, folder='sample_calc'):
        assert self.paramNames is None
        folder = self.workingFolder + os.sep + folder + os.sep + '0'
        if key is None or key == self.feffCards:
            shutil.copytree(self.workingFolder + os.sep + 'extreme' + os.sep + 'extreme_0', folder)
        else:
            if os.path.exists(folder): shutil.rmtree(folder)
            os.makedirs(folder, exist_ok=True)
            shutil.copyfile(self.primitivePath, folder + os.sep + 'molecule.xyz')
            create_input(folder, folder + os.sep + 'molecule.xyz', key)
        with open(os.path.join(folder, 'params.txt'), 'w') as f: json.dump([['dummy',0]], f)

    def generateDeformationsAdaptive(self, folder, maxSampleSize=500, maxError=0.01, calcSampleInParallel=4, seed=0):
        if self.debug: print('Run adaptive sampling')
        sampling.sampleAdaptively(self.paramRanges, self.moleculeConstructor, self.feffCards, maxError=maxError, spectralProgram='feff6', samplePreprocessor=samplePreprocessor, workingFolder=folder, seed=seed, outputFolder=folder+'_result', debugFolder=folder+'_debug', runConfiguration={'calcSampleInParallel': calcSampleInParallel}, adaptiveSamplerParams={'initialIHSDatasetSize':100}, maxSampleSize=maxSampleSize, debugLevel=1 if self.debug else 0)

    def extractSelectedPathsForOneFolder(self, shells, inputFolder, outputFolder):
        all_atoms_info = read_xyz_file_with_name_index(inputFolder + os.sep + 'molecule.xyz')
        all_atoms = all_atoms_info.loc[:, ['x', 'y', 'z']].to_numpy()
        all_atom_names = all_atoms_info.index
        for file1 in sorted(os.listdir(inputFolder)):
            if not file1.endswith(".dat"): continue
            if file1 in ['files.dat', 'nstar.dat', 'paths.dat', 'chi.dat', 'crit.dat']: continue
            if file1[:4] != 'feff':
                if feff.isSuccessful(inputFolder):
                    assert file1[:4] == 'feff', f'{inputFolder}/{file1} must be path file'
                else:
                    warnings.warn(f'Feff error in folder {inputFolder}')
                    continue
            with open(os.path.join(inputFolder, file1), "r") as feff_file:
                lines_feff = feff_file.readlines()
            path_atoms = read_feff_path(lines_feff).to_numpy()
            path_name = get_path_name(path_atoms, all_atoms, all_atom_names)
            for i in range(len(shells)):
                sh = shells[i]
                for p in ['sAngle', 'sReff', 'secondCheckR', 'secondCheckA', 'limitsA', 'limitsR']:
                    if p not in sh: sh[p] = self.defaultPathGroupParams[p]
                isEqual = [compare_path0(name0, path_name, self.primitive, sAngle=sh['sAngle'], sReff=sh['sReff'], secondCheckR=sh['secondCheckR'], limitsR=sh['limitsR'], secondCheckA=sh['secondCheckA'], limitsA=sh['limitsA']) for name0 in sh['paths']]
                if np.any(isEqual):
                    folder_deform_path = outputFolder +os.sep+ sh['name'] + os.sep + os.path.split(inputFolder)[-1]
                    attributes0 = extract_FEFF_attributes(lines_feff)
                    attributes1 = recoverDLR(inputFolder, file1)
                    if not os.path.exists(folder_deform_path):
                        os.makedirs(folder_deform_path)
                    attributes0.to_csv(folder_deform_path + os.sep + path_name + '.csv', index=False, header=True)
                    attributes1.to_csv(folder_deform_path + os.sep + path_name + '_attr_' + '.csv', index=False, header=True)
                    if not os.path.isfile(folder_deform_path + os.sep + 'molecule.csv'):
                        all_atoms_info.to_csv(folder_deform_path + os.sep + 'molecule.csv', index=True, index_label='name', header=True)
    
    def extractSelectedPaths(self, shells=None, folder='sample_calc'):
        """
        For each folder inside the "sample folder", this function reads each path and extract only those paths specified in the "shell" argument. The paths, the paths attributes and the deformed molecule, will be save in the "Training" dataset under every subfolder (e.g. 0000, 0001, ...)
        """
        if shells is None: shells = self.shells
        self.managmentHistory.append({'function': 'extractSelectedPaths', 'shells':copy.deepcopy(shells)})
        folder = self.workingFolder + os.sep + folder
        folder_deform = sorted([f.path for f in os.scandir(folder) if f.is_dir()]) #folders in "folder"
        folder_ML=os.path.join(self.workingFolder, 'training')
        if os.path.exists(folder_ML): shutil.rmtree(folder_ML)
        os.makedirs(folder_ML, exist_ok=True)
        for shell in shells:
            os.makedirs(os.path.join(folder_ML,shell['name']), exist_ok=True)
        for m in range(len(folder_deform)):
            self.extractSelectedPathsForOneFolder(shells, inputFolder=folder_deform[m], outputFolder=folder_ML)

    def separateAllPaths(self, mdic):
        """
        mdic is a dictionary containing the grouped numbered elements e.g mdic={'Ru':['Ru.0'],'C':['C.1','C.2', ...]}
        """
        primitive=self.primitive
        for name1 in self.ir_names:
            nameSplit=name1.split('_')
            elements=[]
            for name2 in nameSplit:
                elements.append(name2.split('.')[0]) # e.g. Ru, C, O without position numbers
            selected=[]
            for name2 in elements:
                selected.append(mdic[name2]) # e.g selected [['Ru.0'],['C.1','C.2','C.3'], ...]
            comb=list(itertools.product(*selected))
            r,c=np.shape(comb)
            pathsComb=[]
            for i in range(r):
                singlePathComb=[]
                for j in range(c):
                    singlePathComb.append(comb[i][j])
                pathsComb.append(singlePathComb) # e.g. pathsComb:  [['Ru.0', 'C.1'], ['Ru.0', 'C.2'], ['Ru.0', 'C.3']]
            listPaths=[]
            for name3 in pathsComb:
                first=name3[0]
                if first.split('.')[1]=='0':
                    listPaths.append('_'.join(name3)) #e.g. listPath=['Ru.0_C.1', 'Ru.0_C.2', 'Ru.0_C.3']
                else:
                    continue
            listPathsRev=[]
            for name4 in listPaths:
                listPathsRev.append(revert_name(name4))
            listDiff=[]
            for name5 in listPathsRev:
                if name5 not in listPaths:
                    listDiff.append(name5)
            if listDiff==[]:
                listAll=listPaths
            else:
                listAll=listPaths+listDiff

            selected_path=[]
            for name6 in listAll:
                if compare_path0(name1,name6,primitive)==True:
                    selected_path.append(name6)
                else:
                    continue
            print(selected_path)
            self.separateSelectedPaths(selected_path)
        self.deleteRedundantFolder()

    def renamePathGroup(self, oldName, newName):
        self.managmentHistory.append({'function': 'renamePathGroup', 'oldName':oldName, 'newName':newName})
        set_folder = os.path.join(self.workingFolder, 'training')
        assert os.path.exists(set_folder)
        assert os.path.exists(set_folder+os.sep+oldName), f'Path group {oldName} doesn\'t exist'
        assert not os.path.exists(set_folder+os.sep+newName), f'Path group {newName} already exists'
        os.rename(set_folder+os.sep+oldName, set_folder+os.sep+newName)

    def deletePathGroup(self, name):
        self.managmentHistory.append({'function': 'deletePathGroup', 'name':name})
        set_folder = os.path.join(self.workingFolder, 'training')
        assert os.path.exists(set_folder)
        assert os.path.exists(set_folder+os.sep+name), f'Path group {name} doesn\'t exist'
        shutil.rmtree(set_folder+os.sep+name)
        
    def getPathGroupNames(self):
        set_folder = os.path.join(self.workingFolder, 'training')
        assert os.path.exists(set_folder)
        return os.listdir(set_folder)

    def separateSelectedPaths(self, group, selectedPaths):
        """
        The group is the group of paths where you want to make the separation. The function will search inside every subfolder of the Training folder those paths indicated in the selected_paths
        If one path is not present because of the symmetry, the function will search, among the present paths, the one with the Reff as close to the one corresponding to the searched path.

        My question to Andrea: Do we really need this step? Why we can't set shells beforehand in extractSelectedPaths to separate these paths?
        Answer: That is a good point. I thought about that too time ago, however I see just a possible issue. Suppose that you set three paths, e.g. Cu.0->C.1, Cu.0->C.2 and Cu.0->C.3. There could be a structure where, for two paths their distances would be similar and in this case we would have: Cu.0->C.1 with degeneration 2, Cu.0->C.2 with degeneration 0 (no paths) and Cu.0->C.3 with degeneration 1. Since, for each geometry we need to find at least one path, we should define already at this stage the criterion for path splitting present in c.separeSelectedPaths('Ti_Cl_g', paths). I think that it is possible to do it but it should make the code slower since the function should run on every defined group of paths.
        """
        self.managmentHistory.append({'function': 'separateSelectedPaths', 'group':group, 'selectedPaths':copy.deepcopy(selectedPaths)})
        set_folder=os.path.join(self.workingFolder, 'training')
        folder_inside=sorted([f.path for f in os.scandir(os.path.join(set_folder,group)) if f.is_dir() ])
        
        path_directory=os.path.join(set_folder, selectedPaths['group'])
        if os.path.exists(path_directory): shutil.rmtree(path_directory)
        os.makedirs(path_directory, exist_ok=True)
        
        for m in range(len(folder_inside)):
            
            name_inside=os.path.basename(folder_inside[m])
            if os.path.exists(os.path.join(path_directory,name_inside)):
                shutil.rmtree(os.path.join(path_directory,name_inside))
            os.makedirs(os.path.join(path_directory,name_inside), exist_ok=True)
            
            name_paths=[]
            reff_name_paths=[]
            for files in sorted(os.listdir(folder_inside[m])):
                files_sep=files.split('_')
                if 'attr' in files_sep:
                    name_paths.append('_'.join(files_sep[:len(files_sep)-2]))
                    reff_name_paths.append(pd.read_csv(folder_inside[m]+os.sep+files)['reff'].values[0])
            
            for name in selectedPaths['paths']:
                if searchName(name, name_paths):
                    attr0=pd.read_csv(os.path.join(folder_inside[m],name+'.csv'))
                    attr1=pd.read_csv(os.path.join(folder_inside[m],name+'_attr_'+'.csv'))
                    attr1['deg']=attr1['deg']-1
                    attr1.to_csv(os.path.join(folder_inside[m],name+'_attr_'+'.csv'),index = False, header=True)
                    attr2=copy.copy(attr1)
                    attr2['deg']=1
                    attr2.to_csv(os.path.join(os.path.join(path_directory,name_inside),name+'_attr_'+'.csv'),index = False, header=True)
                    attr0.to_csv(os.path.join(os.path.join(path_directory,name_inside),name+'.csv'),index = False, header=True)
                    if attr1['deg'].values[0]==0:
                        os.remove(os.path.join(folder_inside[m],name+'_attr_'+'.csv'))
                        os.remove(os.path.join(folder_inside[m],name+'.csv'))
                else:
                    name2=serchNameReff(name,name_paths,reff_name_paths,folder_inside[m])
                    attr0=pd.read_csv(os.path.join(folder_inside[m],name2+'.csv'))
                    attr1=pd.read_csv(os.path.join(folder_inside[m],name2+'_attr_'+'.csv'))
                    attr1['deg']=attr1['deg']-1
                    attr1.to_csv(os.path.join(folder_inside[m],name2+'_attr_'+'.csv'),index = False, header=True)
                    attr2=copy.copy(attr1)
                    attr2['deg']=1
                    attr2.to_csv(os.path.join(os.path.join(path_directory,name_inside),name+'_attr_'+'.csv'),index = False, header=True)
                    attr0.to_csv(os.path.join(os.path.join(path_directory,name_inside),name+'.csv'),index = False, header=True)
                    if attr1['deg'].values[0]==0:
                        os.remove(os.path.join(folder_inside[m],name2+'_attr_'+'.csv'))
                        os.remove(os.path.join(folder_inside[m],name2+'.csv'))

    def applyManagementHistory(self, history):
        for action in history:
            f = action['function']
            if f == 'separateSelectedPaths':
                self.separateSelectedPaths(group=action['group'], selectedPaths=action['selectedPaths'])
            elif f == 'renamePathGroup':
                self.renamePathGroup(oldName=action['oldName'], newName=action['newName'])
            elif f == 'deletePathGroup':
                self.deletePathGroup(name=action['name'])
            elif f == 'extractSelectedPaths':
                self.extractSelectedPaths(shells=action['shells'])
            else: assert False, 'Unknown function name'

    def deleteRedundantFolder(self):
        """
        If some folders are empty because all the paths have been separed using the separeSelectedPaths, this function will remove them
        """
        set_folder= self.workingFolder + os.sep + 'training'
        folder_inside0=sorted([f.path for f in os.scandir(os.path.join(self.workingFolder, set_folder)) if f.is_dir()])
        for m in range(len(folder_inside0)):
            folder_inside1=sorted([f.path for f in os.scandir(os.path.join(self.workingFolder, folder_inside0[m])) if f.is_dir()])
            v=0
            for n in range(len(folder_inside1)):
                if len(os.listdir(folder_inside1[n]))==1 and sorted(os.listdir(folder_inside1[n]))[0]=='molecule.csv': v=v+1
            if v==len(folder_inside1):
                print('Folder:',folder_inside0[m],'was deleted')
                shutil.rmtree(folder_inside0[m])
    
    def createMLDataset(self, folder='dataset', folder0='sample_calc', step=0.05, Reff=False):
        """
        folder Dataset, folder0=sample. This function creates a dataset folder containing all the training set for the path indicated in the shell dictionary
        """
        self.datasetFolder = os.path.join(self.workingFolder, folder)
        if os.path.exists(self.datasetFolder):
            shutil.rmtree(self.datasetFolder)
        os.makedirs(self.datasetFolder, exist_ok=True)
        training= self.workingFolder + os.sep + 'training'
        folder_paths = sorted([ f.path for f in os.scandir(training) if f.is_dir() ])
        for m in range(len(folder_paths)):
            name_path=os.path.basename(folder_paths[m])
            if os.path.exists(os.path.join(self.workingFolder, folder + os.sep + name_path)):
                shutil.rmtree(os.path.join(self.workingFolder, folder + os.sep + name_path))
            os.makedirs(os.path.join(self.workingFolder, folder + os.sep + name_path), exist_ok=True)
            diag, amp, phs, lam, chi = getAttributes(folder_paths[m], step)
            np.savetxt(self.datasetFolder + os.sep + name_path + os.sep + 'amplitude.txt', np.transpose(amp))
            np.savetxt(self.datasetFolder + os.sep + name_path + os.sep + 'phase.txt', np.transpose(phs))
            np.savetxt(self.datasetFolder + os.sep + name_path + os.sep + 'lambda.txt', np.transpose(lam))
            np.savetxt(self.datasetFolder + os.sep + name_path + os.sep + 'chi.txt', chi)
            diag.to_csv(self.datasetFolder + os.sep + name_path + os.sep + 'report.csv', index=False, header=True)

        goodFolders = sorted(os.listdir(folder_paths[0]))
        valueS=[]
        for fold in goodFolders:
            name,value = getParams(self.workingFolder+os.sep+folder0+os.sep+fold+os.sep+'params.txt')
            #nameS.append(name)
            valueS.append(value)
        assert len(valueS) > 0, 'No good folders in '+self.workingFolder+os.sep+folder0
        self.names=name
        self.values=valueS
        
        parameters=pd.DataFrame(valueS)
        parameters.columns=name
        parameters.to_csv(self.workingFolder + os.sep + folder + os.sep + 'parameters.csv', index=False, header=True)
        
        if Reff:
            for folder1 in folder_paths:
                namefolder1=os.path.basename(folder1)
                reff=isolateReff(folder1)
                np.savetxt(folder+os.sep+namefolder1+os.sep+'Reff.txt', reff)
        utils.saveData(self.paramRanges, self.workingFolder + os.sep + folder + os.sep + 'paramRanges.json')
        return parameters

    def check(self, outputFolder, method='RBF', CVcount=2, feffVersion=6, showInNotebook=True):
        feffCards = self.feffCards
        exafsPredictor = ExafsPredictor(self.datasetFolder, subfolder='')
        k = exafsPredictor.k
        exParam = exafsPredictor.buildExafsParams(enot=0, S02=1, sigma2=0)
        if self.paramNames is None:
            feffSp = feff.parseOneFolder(self.workingFolder + os.sep + 'extreme' + os.sep + 'extreme_0')
            feffSp.chi *= feffSp.k ** 2
            R, spR = exafs.convert2RSpace(feffSp.k, feffSp.chi)
            chi = exafsPredictor.calcTotalSpectrum(None, exafsParams=exParam)
            chi = chi * k ** 2
            R2, spR2 = exafs.convert2RSpace(k, chi)
            plotting.plotToFile(axisMatrix=[[(feffSp.k, feffSp.chi, 'feff', k, chi, 'prediction')],
                                            [(R, np.abs(spR), 'feff', R2, np.abs(spR2), 'prediction')]],
                                axisMatrixKw=[[{'xlim':[0, 12], 'xlabel':'k'}],
                                              [{'xlim':[0.5, 4], 'xlabel':'R(A)'}]],
                                fileName=outputFolder + os.sep +'check.png', showInNotebook=showInNotebook)
            return
        print('Relative to constant prediction error:', exafsPredictor.calcCV(method=method, CVcount=CVcount))
        exafsPredictor.fit(method=method)
        ab = self.paramRanges
        rng = np.random.default_rng(0)
        toCheck = []
        for i in range(20):
            toCheck.append({pn:rng.uniform(*ab[pn]) for pn in ab})
        checkResult = [None]*len(toCheck)
        for i in range(len(toCheck)):
            p = toCheck[i]
            mol = self.moleculeConstructor(p)
            folder = feff.generateInput(mol, feffVersion=feffVersion, **feffCards)
            feff.runLocal(folder, feffVersion=feffVersion)
            feffSp = feff.parseOneFolder(folder)
            feffSp.chi *= feffSp.k ** 2
            R, spR = exafs.convert2RSpace(feffSp.k, feffSp.chi)
            chi = exafsPredictor.calcTotalSpectrum(geomParams=p, exafsParams=exParam)
            chi = chi * k ** 2
            R2, spR2 = exafs.convert2RSpace(k, chi)
            spR2_ = np.interp(R, R2, spR2)
            ind = (0.5<=R) & (R<=4)
            error = np.linalg.norm(spR[ind]-spR2_[ind]) / np.linalg.norm(spR[ind])
            checkResult[i] = {'param':p, 'molecule':mol, 'feff-k':(feffSp.k, feffSp.chi), 'feff-r':(R, np.abs(spR)), 'prediction-k':(k,chi), 'prediction-r':(R2, np.abs(spR2)), 'error':error}
        checkResult = sorted(checkResult, key=lambda r: r['error'], reverse=True)
        wf = outputFolder+os.sep+'worst'
        feff.generateInput(checkResult[0]['molecule'], folder=wf, feffVersion=feffVersion, **feffCards)
        checkResult[0]['molecule'].export_xyz(wf+os.sep+'molecule.xyz')
        feff.runLocal(wf, feffVersion=feffVersion)
        for i, r in enumerate(checkResult):
            if i == 0:
                oneMolManager = MSPathGroupManager(workingFolder=outputFolder+os.sep+'worst_path_manager', constraints=self.constraints, feffCards=self.feffCards, molecule=wf+os.sep+'molecule.xyz', debug=False)
                oneMolManager.applyManagementHistory(self.managmentHistory)
                oneMolManager.createMLDataset()
                oneMolManager.check(outputFolder+os.sep+'worst_path_manager'+os.sep+'check', showInNotebook=False)
                oneMolData = plotting.readPlottingFile(outputFolder+os.sep+'worst_path_manager'+os.sep+'check'+os.sep+'check.txt')
                def plotMoreFunction0(ax): ax.plot(oneMolData['prediction x'], oneMolData['prediction y'], label='chosen paths')
                def plotMoreFunction1(ax): ax.plot(oneMolData['prediction x_1'], oneMolData['prediction y_1'], label='chosen paths')
            else: plotMoreFunction0,plotMoreFunction1 = None,None
            title = 'The worst approximation. Error=%.3g' % r['error'] + ' ' + utils.dict2str(r['param'])
            plotting.plotToFile(axisMatrix=[[(*r['feff-k'], 'feff', *r['prediction-k'], 'prediction')],
                                            [(*r['feff-r'], 'feff', *r['prediction-r'], 'prediction')]],
                                axisMatrixKw=[[{'xlim':[0, 12], 'xlabel':'k', 'title':title, 'plotMoreFunction':plotMoreFunction0}],
                                              [{'xlim':[0.5, 4], 'xlabel':'R(A)', 'plotMoreFunction':plotMoreFunction1}]],
                                fileName=outputFolder + os.sep + utils.zfill(i, len(checkResult)) + '_check.png', showInNotebook=showInNotebook)
            if utils.isJupyterNotebook(): break


def fitExafsParams(xyzFile, expSpectrum, constraints, debug, workingFolder, defaultParams={}, multipleS02=False, separatePathGroups=[], renamePathGroups={}, deletePathGroups=[], kPower=2, RSpaceParams={'kmin':2, 'kmax':10}, exafsParamsFuncArgBounds={'enot':[-10,10], 'S02':[0, 2], 'sigma2':[0.001,0.04]}):
    """
    Set at first debug=True to select multiple scattering paths and separate path groups if necessary.
    Then set debug=False to get sliders

    :param separatePathGroups: break apart existed path group [('existedPathGroupName', {'group':'newGroupName', 'paths':[path names]}), ...]
    :param renamePathGroups: dict {'existedPathGroupName': 'newPathGroupName'} to rename path groups
    :param deletePathGroups: list of path group names to delete
    """
    pathManager = MSPathGroupManager(workingFolder=workingFolder, constraints=constraints, molecule=xyzFile, debug=debug)
    pathManager.extractSelectedPaths()
    for existedPathGoupName, newPathInfo in separatePathGroups:
        pathManager.separateSelectedPaths(existedPathGoupName, newPathInfo)
    for oldPathGoupName, newPathGroupName in renamePathGroups.items():
        pathManager.renamePathGroup(oldPathGoupName, newPathGroupName)
    for name in deletePathGroups: pathManager.deletePathGroup(name)
    if debug:
        print('Final path groups:', pathManager.getPathGroupNames())
    pathManager.createMLDataset()
    if debug:
        pathManager.check(workingFolder+os.sep+'check')
        return
    exafsPredictor = ExafsPredictor(workingFolder)
    exafsParamsFuncArgBounds = copy.deepcopy(exafsParamsFuncArgBounds)
    if multipleS02:
        for g in exafsPredictor.pathGroups: exafsParamsFuncArgBounds[f'S02_{g}'] = exafsParamsFuncArgBounds['S02']
        del exafsParamsFuncArgBounds['S02']
    for g in exafsPredictor.pathGroups:
        exafsParamsFuncArgBounds[f'sigma2_{g}'] = exafsParamsFuncArgBounds['sigma2']
    del exafsParamsFuncArgBounds['sigma2']
    exafsModel = exafsPredictor.getFuncModel(expSpectrum=expSpectrum, kPower=kPower, exafsParamsFuncArgBounds=exafsParamsFuncArgBounds, multipleS02=multipleS02, RSpaceParams=RSpaceParams, name=os.path.splitext(os.path.split(xyzFile)[-1])[0])
    return uiControls.FuncModelSliders(funcModel=exafsModel, debug=False, defaultParams=defaultParams)


def sampleExafs(moleculeConstructor, paramRanges, constraints, debug, workingFolder, separatePathGroups=[], renamePathGroups={}, deletePathGroups=[], sampleSize=None, maxAdaptiveSize=500):
    """
    Set at first debug=True to select multiple scattering paths and separate path groups if necessary.
    Then set debug=False to get sliders

    :param separatePathGroups: break apart existed path group [('existedPathGroupName', {'group':'newGroupName', 'paths':[path names]}), ...]
    :param renamePathGroups: dict {'existedPathGroupName': 'newPathGroupName'} to rename path groups
    :param deletePathGroups: list of path group names to delete
    """
    pathManager = MSPathGroupManager(workingFolder=workingFolder, constraints=constraints, paramRanges=paramRanges, moleculeConstructor=moleculeConstructor, debug=debug, sampleSize=sampleSize, maxAdaptiveSize=maxAdaptiveSize)
    pathManager.extractSelectedPaths()
    for existedPathGoupName, newPathInfo in separatePathGroups:
        pathManager.separateSelectedPaths(existedPathGoupName, newPathInfo)
    for oldPathGoupName, newPathGroupName in renamePathGroups.items():
        pathManager.renamePathGroup(oldPathGoupName, newPathGroupName)
    for name in deletePathGroups: pathManager.deletePathGroup(name)
    if debug:
        print('Final path groups:', pathManager.getPathGroupNames())
    pathManager.createMLDataset()
    pathManager.check(workingFolder+os.sep+'check', method='RBF')


def fitExafsByStructureSliders(sampleFolder, expSpectrum, moleculeConstructor=None, defaultParams=None, multipleS02=False, kPower=2, RSpaceParams={'kmin':2, 'kmax':10}, exafsParamsFuncArgBounds={'enot':[-10,10], 'S02':[0, 2], 'sigma2':[0.001,0.04]}):
    exafsModel = FuncModel.createExafsFittingModel(sampleFolder, expSpectrum, moleculeConstructor=moleculeConstructor, multipleS02=multipleS02, kPower=kPower, RSpaceParams=RSpaceParams, exafsParamsFuncArgBounds=exafsParamsFuncArgBounds)
    if defaultParams is None: defaultParams = {}
    for p in ['enot', 'S02', 'sigma2']:
        uiControls.addNoFit(p, exafsModel.paramProperties, defaultParams)
    return uiControls.FuncModelSliders(funcModel=exafsModel, debug=False, defaultParams=defaultParams)


def fitXanesAndExafsSimultaneously(project, xanesSampleFolder, exafsSampleFolder, exafsSpectrum, defaultParams=None, multipleS02=False, kPower=2, RSpaceParams={'kmin':2, 'kmax':10}, exafsParamsFuncArgBounds={'enot':[-10,10], 'S02':[0, 2], 'sigma2':[0.001,0.04]}):
    exafsModel = FuncModel.createExafsFittingModel(exafsSampleFolder, exafsSpectrum, moleculeConstructor=project.moleculeConstructor, multipleS02=multipleS02, kPower=kPower, RSpaceParams=RSpaceParams, exafsParamsFuncArgBounds=exafsParamsFuncArgBounds)
    xanesModel = FuncModel.createXanesFittingModel(project, ML.readSample(xanesSampleFolder), name='xanes')
    combinedModel = exafsModel.combineWith(xanesModel, newName=project.name, valueMergeRule=[1, 1], dataMergeRule='add axes', createfigParams={'nrows': 2})
    if defaultParams is None: defaultParams = {}
    for p in ['Gamma_hole', 'Gamma_max', 'Ecent', 'Efermi', 'Elarg', 'shift', 'norm', 'enot', 'S02', 'sigma2']:
        uiControls.addNoFit(p, exafsModel.paramProperties, defaultParams)
    return uiControls.FuncModelSliders(funcModel=combinedModel, debug=False, defaultParams=defaultParams)
