# -*- coding: utf-8 -*-
"""
TA_gen2_basic5.py:
Simulation tool of timescale generation

- Input datafiles are read from the list './clock_file_list.txt'.
- Input clock data are prepared in th edirectry '../Data_Clock".
- Output of menu 1-3 are saved in the './TA_RESULT'.

Created on 31 May, 2025
@author: Y. Hanado
"""

import numpy as np
import matplotlib.pyplot as plt
import datetime
import allantools
import os
import shutil

# ==============================================================
# DEF : calculation ADEV by "allantools"
def adev_calc(data_list, menu_n, x_limit):                                          

    line = []
    name = []
    tt = []                                                                       
    xx = []                                                                       
    taus = []
    adevs = []
    data_adev = []
    data_adev_list = []
    MJDspan = []                                                                  

  # search the MJD span of calclated TA   
    for nclk, line in enumerate(data_list):
        if (menu_n =='3'):
            name = line[0]
            if ('TA' in name):
                tt = line[1][0]
                xx = line[1][1]  
                t = []
                x = []
                for i in range(len(tt)):
                    #if ((xx[i]!=0.0) & (np.isnan(xx[i])!=True)):                      
                    if ((xx[i] < x_limit) & (np.isnan(xx[i])!=True)):                   
                        t.append(tt[i])
                        x.append(xx[i])
  
    name = []
    tt = []                                
    xx = []  
  # calculate ADEV for each data of data_list
    for nclk, line in enumerate(data_list):
        name = line[0]
        tt = line[1][0]                                                           
        xx = line[1][1]                                                           
     
      # fix MJD span for ADEV   
        t = []
        x = []
        p_data = []
        E_flag = False
        for i in range(len(tt)):
          # Separate data to portion by before/after missing data
            if ((xx[i] < x_limit) & (np.isnan(xx[i])!=True)):                         
                t.append(tt[i])
                x.append(xx[i])
                E_flag = True            
            else:
                if (E_flag == True):
                    MJD0 = int(t[0])              
                    MJD1 = int(t[-1])             
                    MJDspan = [MJD0, MJD1]        
                    p_data.append([t, x, MJDspan])
                    E_flag = False
                    t = []
                    x = []
          
        if (len(t) != 0):             
            MJD0 = int(t[0])              
            MJD1 = int(t[-1])     
            MJDspan = [MJD0, MJD1]
        # p_data is the portion data list deparated by the data missing 
            p_data.append([t, x, MJDspan])

      # Adev calclation for each portion data  
        x = []                                                                    
        t = []                                                                    
        for i in range(len(p_data)):                                                    
            if (i > 0):                                                           
                name = name + "*" + str(i)                                        
            t = p_data[i][0]                                                      
            x = p_data[i][1]                                                       
            MJDspan = p_data[i][2]                                                
            t_int = float(t[1] - t[0])*86400        
          # np.nan?                                                                 
            npx=np.isnan(x)                                                         
            if ((any(xx >= x_limit for xx in x) != True) 
                & (npx.any() != True)):                                              
              # Adev by allantools
                (taus, adevs, errors, ns) = allantools.oadev(x)
                taus = taus*t_int
                adevs = adevs/t_int
                data_adev = np.array([taus, adevs])
              # data_adev_list[] is [[CS**, [tau, ADEV]],,,]
                data_adev_list.append([name, data_adev, MJDspan])                 
            else:
                print("!!!!! " + name + "  is not available for allantools.")

    return data_adev_list

# ==============================================================
# DEF : plot ADEV results 
def adev_plot(outfile,detdata,wgtfile,all_wgt_list,MJDall,log,comment, 
              menu_n, x_limit):                                                        

  # data_list[] is the data from outfile <== "all_*.dat"
    data_list = []
    data_list = read_jointdata(outfile, MJDall, x_limit)
  
  # pickup only the weighted clock if menu_n=3                                    
    data_list_wgt = []
    if (menu_n == '3'):
        for nclk, line in enumerate(data_list):
          # data_list[-3]=wgt, data_list[-2]=flg                                  
            if (line[-3] != 0 or line[-2] == 3):                                     
                data_list_wgt.append(line)
        data_list = data_list_wgt        
  
  # calculate ADEV by allantools
  # data_adev_list[] is [["UTCk-CSxx", [[tau,,,],[ADEV,,,]]],,,]
    data_adev_list = []
    data_adev_list = adev_calc(data_list, menu_n, x_limit)                           

  # format Adev data-set for plotting   ..........................                
    tau_list = data_adev_list[0][1][0]
    for nclk in range(1,len(data_adev_list)):
        tau_nclk = data_adev_list[nclk][1][0]
        if (max(tau_nclk) > max(tau_list)):
            tau_list.append(max(tau_nclk), max(tau_nclk))                         
    for n in range(1,len(data_adev_list)):
        tau_clk= data_adev_list[n][1][0]
        adv_clk= data_adev_list[n][1][1]
        if (max(tau_clk) < max(tau_list)):
            tau_clk = np.append(tau_clk, max(tau_list))
            adv_clk = np.append(adv_clk, np.nan)
        data_adev_list[n][1] = [tau_clk,adv_clk]
        
  # plot graph
    fig = plot_data(data_list, data_adev_list, menu_n)
    plt.show(block=False)

  # plot `all_weight.dat'                                                            
    w_fig = ""                                                                       
    if (menu_n == '3'):                                                              
        w_fig = plot_weight(all_wgt_list)                                             
        plt.show(block=False)                                                        
 
  # save? or not?
    #print ('\n' + comment)                                                        
    while (save_yn := input ("????? Save figure? (y/n) = ")) not in ('y','n'):                  
        print("!!!! Input is not correct. Please try again!" + '\n')                            
    if (save_yn == "y"):
        log, fig = save_logfigdat(outfile, data_adev_list, detdata, wgtfile,                    
                                  log, fig, w_fig)                                   

    return

# ==============================================================
# DEF : average clock data
def do_average(infile_all, MJDall, log, x_limit):
   
  # initialize
    wgt = []
    x = []
    MJD = []
    avedata = []
    data_ave_list = []
    all_wgt_list = []                                                                 

  # data_ave_list[] is the data from infile_all <== "all_*_clk.dat"
    data_ave_list = read_jointdata(infile_all, MJDall, x_limit)
 
  # read data from infile_all <== "all_*_clk.dat"
    with open(infile_all, "r") as infile:
        for i, strline in enumerate(infile):
          # line = [column1, column2,,,,]
            line = strline.split()
            if (line[0] =='fname'):                                                    
                clk = [line[j].split('.')[0].split('-')[1] for j in range(1, len(line))]                   
            if (line[0] =='weight'):                                                   
                wgt = [float(line[j]) for j in range(1, len(line))]                    
            if (line[0] =='cid'):                                                      
                cid = [int(line[j]) for j in range(1, len(line))]                                     
            if (line[0].isdigit()):                                                    
                t = float(line[0])
                if (t >= MJDall[0]) & (t <= MJDall[-1]):
                    x = [float(line[j]) for j in range(1, len(line))] 
                  # caluclate average
                    ave = 0
                    sum_wgt = 0
                    wgt0 = [w for w in wgt]
                    for k in range(len(x)):
                        wgt0[k] = wgt[k] 
                        if ( x[k]<x_limit and np.isnan(x[k])!=True ):                  
                            ave = ave + wgt0[k] * x[k]
                            sum_wgt = sum_wgt + wgt0[k]
                        else:                                                          
                            wgt0[k] = 0.0                                               
                  # if all weights are 0                                             
                    if (sum_wgt == 0):
                        print(t, " !!!!! No available clock in averaging !!!!")
                    else:
                        ave = ave/sum_wgt                      
                        wgt1 = [w/sum_wgt for w in wgt0]                                
                  # MJD[] is [MJD0,,,] and avedata[] is []
                    MJD.append(t)
                    avedata.append(ave)
                    all_wgt_list.append([t,clk,wgt1,cid])                              

  # data_ave_list[] is [["UTCk-TA", [MJD, avedata], "0", "3"],[],,]               
    ave_data = np.array([MJD, avedata])
    data_ave_list.append(["UTCk-TA", ave_data, "0.0", "3", "0"])                  

  # joint the clock data & their average ==> outfile (= "all_*_ave.dat")
    outfile = infile_all.split(".")[0] + "_ave" + ".dat"
    outfile_p, log = joint_rawdata(outfile, MJDall, data_ave_list, log)
    
  # save the clock weights data
    wgtfile = save_weight(all_wgt_list)                                                

  # comment to log[]n
    comment = "**3** Weighted average *****" + '\n'
    log_0 = "::::: Input data file  : " + infile_all + '\n'
    log_1 = "::::: Saved joint data : " + outfile_p + '\n'            
    log_2 = "::::: Weight of clocks : " + wgtfile + '\n'
    log = log + [comment + log_0 + log_1 + log_2]

    return outfile_p, log, comment, wgtfile, all_wgt_list                        

# ==============================================================
# DEF : de-trend for eack clock file 
def do_detrend(infile_all, MJDall, log, x_limit):

  # print the comment
    comment = "**2** Full-span De-trend of each clock data *****" + '\n'               
    print(comment) 
    
  # data_list[] is the data from infile_all <== "all_raw_clk.dat"
    data_list = []
    data_list = read_jointdata(infile_all, MJDall, x_limit)

  # de-trend for each clock file
    det_clk = []
    data_det_list = []
    log_rate = []                                                                 
    for nclk in range(len(data_list)):
        name = data_list[nclk][0]
        name = name.split('.')[0] + '_det'
        t0 = data_list[nclk][1][0]
        x0 = data_list[nclk][1][1]
      # adjust the data range to match with MJDall[]        
        t = np.array([t0[n] for n in range(len(t0)) 
                      if (t0[n] >= MJDall[0]) & (t0[n] <= MJDall[-1])])
        x = np.array([x0[n] for n in range(len(t0)) 
                      if (t0[n] >= MJDall[0]) & (t0[n] <= MJDall[-1])])
        wgt = data_list[nclk][2]
        flg = data_list[nclk][3]
        cid = data_list[nclk][4]                                                  
        nx = []
        xint = []
        rate = []
      # data in the full data range for linear fitting 
        tt = np.array([t[n] for n in range(len(t)) if (x[n] < x_limit)])
        xx = np.array([x[n] for n in range(len(t)) if (x[n] < x_limit)])
      # linear fit
        slope, xint = np.polyfit(tt, xx, 1)
        
      # np.nan?                      
        if (np.isnan(slope) == False) and (np.isnan(xint) == False):                
          # nx[] is the remains from the fitting line
            nx = x - (xint + slope * t)
            log_rate.append([name, slope, xint, tt[0], tt[-1]])                         
            rate = np.array([slope for i in range(len(t))])  # all elements of rate[] are the same in full-range
          # ndata[] is the remains for data_det_list[]
            ndata = [t, nx]
          # ndata2[] is for det_clk[]
            ndata2 = [name, t, x, rate, nx]
          # det_clk[] is [[name, [MJD], [rawdata], [rate], [remain]],,,]
            det_clk.append(ndata2)
        else:
            print("!!!!! " + name + "  is not available for de-trend") 
    
        data_det_list.append([name, ndata, wgt, flg, cid])
    
  # joint clock data and output ==> "all_det_clk.dat"
    outfile = "all_det_clk.dat"
    all_det_data, log = joint_rawdata(outfile, MJDall, data_det_list, log)
        
  # comment to log[]n
    log_0 = "::::: input data file  : " + infile_all + '\n'
    log_1 = "::::: saved joint data : " + outfile + '\n'
    log_2 = "::::: clock rate       :    int         slope     period for slope" +'\n'                  
    for i in range(len(log_rate)):                                                
        log_slope = '{:+.6e}'.format(log_rate[i][1])                              
        log_xint  = '{:+.6e}'.format(log_rate[i][2])                              
        log_t0    = '{:.1f}'.format(log_rate[i][3])                              
        log_tz    = '{:.1f}'.format(log_rate[i][4])                                  
        log_2 = log_2 + "      " + log_rate[i][0] + " : " + str(log_slope)               
        log_2 = log_2 + " " + str(log_xint)                                       
        log_2 = log_2 + " " + str(log_t0) + " - " + str(log_tz) +'\n'                  
    log = log + [comment + log_0 + log_1 + log_2]                                  

    return all_det_data, data_det_list, det_clk, log, comment

# ==============================================================
# DEF: joint single column of clock data 
def joint_rawdata(outfile, MJDall, data_list, log):
 
  # set MJDall[] as the first element of all_data[]
    all_data = []
    all_data.append(MJDall)     # all_data = [[56xx0,56xx1,,,]]
    output_line0 = ["fname "]                                                      
    output_line1 = ["weight"]
    output_line2 = ["flag " ]
    output_line3 = ["cid  " ]                                                      
       
  # read clock data from data_list[]
    for nclk, line in enumerate(data_list):
        name = line[0]          # name is "UTCk-CSxx"
        t = line[1][0]          # t[] is list of MJD
        x = line[1][1]          # x[] is list of UTCk_CS data
        wgt = line[2]           # wgt is a value of weight (not list)
        flg = line[3]           # flg is a value (not list)
        cid = line[4]           # cid is a value of clock id (not list)                 
                                 
        output_line0.append(" " + name + " ")
        output_line1.append("      " + str(wgt) + "      ")
        output_line2.append("       " + str(flg) + "       ")
        output_line3.append("       " + str(cid) + "       ")                      
        
      # Missing data => np.nan                                                     
        xx=[]
        for i in range(len(MJDall)):
            aruyo = (MJDall[i] in t)
            if (aruyo == True): 
                for j in range(len(t)):
                    if (t[j] == MJDall[i]): 
                        xx.append(x[j])
                        break
            else: 
                xx.append(np.nan)
            
      # all_data[] is a joint of each column of clock data
      # all_data=[[MJD,,,],[c#1data,,],[C#2data,,,],,,]
        col_clock = []
        for i in range(len(MJDall)):
            col_clock = np.hstack((col_clock,xx[i])) # col_clock is a column of UTCk-CS data
        all_data.append(col_clock)                   # all_data is multi-columns of UTCk-CS data
       
  # save the joint data file ==> "temp.dat"
    n_lin = len(col_clock)
    n_col = len(all_data)
    output = np.transpose(all_data)
    np.savetxt("temp.dat", output.reshape((n_lin,n_col)),
               fmt='  '.join(['%i'] + ['%+.6e']*(n_col-1)))

  # save the joint data file with title & weight & flag ==> outfile
    with open(outfile, "w") as outf:
        outf.writelines(output_line0 + ['\n'])
        outf.writelines(output_line1 + ['\n'])
        outf.writelines(output_line2 + ['\n'])
        outf.writelines(output_line3 + ['\n'])                                    
        with open("temp.dat", "r") as inf:
            for line in inf:
                outf.writelines(line)

    return outfile, log

# ==============================================================
# DEF: assign the MJD range within available one  
def MJD_limit(MJDall, log):
    
  # input the param of MJD range as "MJD1 MJD2" from the console 
    chk1 = False
    while chk1 == False:
        try:
            comment = ("\n**1** Assign the MJD range for calculation in [" 
                       + str(MJDall[0]) + "," + str(MJDall[-1]) + '].\n')
            MJD1, MJD2 = map(int, input(comment + "\n??1?? Input MJD1 "
                                        "and MJD2 with spacing = ").split()) 
        except ValueError:
            print('!!! input ValueError, try again !!!')
        else:
            print('\n**1** Selected range is : [', MJD1, ' - ', MJD2, ']')
            if MJD1 < MJDall[0]:
                print("!!!!! MJD1 should be equal or larger than ",           
                      MJDall[0], ". !!!!!")                                        
            elif MJD2 > MJDall[-1]:
                print("!!!!! MJD2 should be equal or smaller than ",          
                      MJDall[-1], ". !!!!!")                                        
            elif MJD1 > MJD2:
                print("!!!!! MJD1 should be smaller than MJD2. !!!!!")               
            else:
                chk1 = True

  # MJDall[] is common MJD range for all clockdata 
    MJDall = [MJDall[i] for i in range(len(MJDall)) 
        if (MJDall[i] >= MJD1) & (MJDall[i] <= MJD2)]
    MJDall = np.array(MJDall)

  # comment to log[]n
    comment = "**1** Assign the MJD range for calculation *****" + '\n'
    log_0 = ("::::: Selected MJD     :  [" + str(MJD1) + " - " + str(MJD2) 
             + ']\n')
    log = log + [comment + log_0]

    return MJDall, log

# ==============================================================
# DEF : maximum MJD for all selected clock file
def MJD_max(data_list):
    
  # search the Maximum MJD range over the data_list
    for nclk, line in enumerate(data_list):
        MJD_min_0 = int(np.amin(line[1][0]))
        MJD_max_0 = int(np.amax(line[1][0]))
        if (nclk == 0):
            MJD_min = MJD_min_0
            MJD_max = MJD_max_0
        if (MJD_min_0 <= MJD_min):
            MJD_min = MJD_min_0
        if (MJD_max_0 >= MJD_max):
            MJD_max = MJD_max_0

  # make MJDall[] in range of (MJD_min, MJD_max+1)
    MJDall = np.array([i for i in range(MJD_min, MJD_max+1, 1)])

    return MJDall

# ==============================================================
# DEF : plot the original & Allan deviation
def plot_data(data_list1, data_list2, menu_n):

  # set the plotting area
    fig = plt.figure(figsize=(12, 7))   

  # left box position is set as [left-x, ground-y, x-width, ywidth]
    rect = [0.05, 0.1, 0.35, 0.5] #0307
    xlabel = 'measurement time (MJD)'
    ylabel = 'time difference (sec)'
    ax1 = plt.axes(rect)
    ax1 = plot_set(ax1, xlabel, ylabel)

  # right box position is set as [left-x, ground-y, x-width, ywidth]
    rect = [0.48, 0.1, 0.35, 0.5] #0307
    xlabel = 'averaging time (sec)'
    ylabel = 'allan deviation'
    ax2 = plt.axes(rect)
    ax2 = plot_set(ax2, xlabel, ylabel)

  # set the colortable                                                            
    colors = plt.cm.tab20(np.linspace(0, 1, 20))                                  

  # set the fig-param in left box
    clk_param = []                                                                   
    for i in range(len(data_list1)):
        name1 = data_list1[i][0]
        data1 = data_list1[i][1]
        cwgt1 = data_list1[i][2]                                                     
        cflg1 = data_list1[i][3]                                                     
        cid1  = data_list1[i][4]                                                     
        c_color1 = colors[cid1]                                                      
        if ("TA" in name1): c_color1 = "black"                                    
        if (cflg1 == 1):                                                          
            refclk = name1.split('-')[1].split('.')[0]                               
        clk_param.append([name1,c_color1,cwgt1, cflg1, cid1])                        
        figx = data1[0]
        figy = data1[1]
        figy = [None if k == 1000 else k for k in figy]
        ax1.set_xlim(xmin=data1[0][0], xmax=data1[0][-1])
        ax1.plot(figx, figy, color=c_color1, marker="o", label=name1) 
        ax1.legend(bbox_to_anchor=(0.28, 1.05, 0.38, 0.2), loc='center', ncol=3) 
                  
    if (menu_n == 0):                                                                
        ax１.text(0.01, 0.01, f'${refclk}$ = reference clock', fontsize=10,                  
                 color='dimgray',transform=ax1.transAxes) 

  # set the fig-param in right box
    for i in range(len(data_list2)):
        name2 = data_list2[i][0]
        data2 = data_list2[i][1]
        for j in range(len(clk_param)):                                              
          # if multi ADEV in one clock                                            
            if ('*' in name2):                                                    
                name2a = name2.split('*') [0]                                     
                name2b = name2.split('*') [1]                                     
            else:                                                                 
                name2a = name2                                                    
                name2b = 0                                                        
            if (name2a == clk_param[j][0]):                                       
                c_color2 = clk_param[j][1]                                           
            if ("TA" in name2): c_color2 = "black"                                    
        figx = data2[0]
        figy = data2[1]
        ax2.set_xlim(xmin=data2[0][0], xmax=data2[0][-1])
        ax2.loglog(figx, figy, color=c_color2, marker=f'${name2b}$', label=name2)  
        ax2.legend(bbox_to_anchor=(0.28, 1.05, 0.38, 0.2), loc='center', ncol=3)  

    if (menu_n == '3'):                                                        
        ax2.text(0.01, 0.01, 'All ADEVs are calculated within the MJDs' 
                 +' of wgtd-TA.', fontsize=9, color='red',transform=ax2.transAxes) 

    return fig

# ==============================================================
# DEF : set the basic plotting parameter
def plot_set(ax, xlabel, ylabel):

  # setting for the figure parameters
    ax.grid(which='both', color='gray', linestyle='dotted')
    ax.set_xlabel('x: ' + xlabel)
    ax.set_ylabel('y: ' + ylabel)
    ax.ticklabel_format(style='sci', axis='y', scilimits=(0, 0))

    return ax

# ==============================================================                     
# DEF : plot the weights 
def plot_weight(all_wgt_list):                                                       

  # set the plotting area                                                            
    fig = plt.figure(figsize=(12, 7))                                                 

  # left box position is set as [left-x, ground-y, x-width, ywidth]                  
    rect = [0.05, 0.1, 0.35, 0.5]                                                     
    xlabel = 'time (MJD)'                                                            
    ylabel = 'weight'                                                                
    ax1 = plt.axes(rect)                                                               
    ax1 = plot_set(ax1, xlabel, ylabel)                                              

  # set the colortable                                                               
    colors = plt.cm.tab20(np.linspace(0, 1, 20))                                     
  
  # reconstract the lists from all_wgt_list
    clk = all_wgt_list[0][1]                                                         
    cid = all_wgt_list[0][3]                                                         
    mjd = [t[0] for t in all_wgt_list]                                               
    wgt = [t[2] for t in all_wgt_list]                                               
    wgts = []                                                                        
    for n in range(len(clk)):                                                        
        xx = [x[n] for x in wgt]                                                     
        wgts.append(xx)                                                              
    data = [mjd, wgts]                                                               
  
  # set the fig-param in left box                                                    
    for i in range(len(clk)):                                                        
        c_color1 = colors[cid[i]]                                                                 
        figx = data[0]                                                               
        figy = data[1][i]                                                            
        ax1.set_xlim(xmin=figx[0], xmax=figx[-1])                                     
        ax1.plot(figx, figy, color=c_color1, marker=".", label=clk[i])                  
        ax1.legend(bbox_to_anchor=(0.28, 1.01, 0.38, 0.2), loc='center', ncol=4)                  

    return fig                                                                       

# ==============================================================
# DEF : read jointed all clock data
def read_jointdata(infile_all, MJDall, x_limit):

  # data_org[] is the data of infile_all <== "all_*_clk.dat"
    data = []
    data_list = []
    data_org = np.genfromtxt(infile_all, dtype=str, invalid_raise=False)

  # data_org[] is the data of infile_all
    for ncol in range(len(data_org[0])):  # data_org[0] is ['MJD' 'UTCk-CS**.txt' ,,,]
        name = data_org[0][ncol]
        if (ncol == 0):                   # in the case of 'MJD'                  
            wgt = 0
            flg = 0
            cid = 0                                                               
            t = [float(data_org[n][ncol]) for n in range(4, len(data_org))]
        else:                             # in the case of clock data
            wgt = float(data_org[1][ncol])
            flg = int(data_org[2][ncol])
            cid = int(data_org[3][ncol])                                          
            x = [float(data_org[n][ncol]) for n in range(4, len(data_org))]
            data = np.array([t, x])
          # data_list[] is [[name, [MJD, UTCk-CS**], weight, flag, cid],[],,,]
            data_list.append([name, data, wgt, flg, cid])                         
          # message of phase jump
            #if (x_limit in x):                                                      
            if (any(xx >= x_limit for xx in x) == True):                             
                print("\n" + "!!!!! " + name + "  includes outliers")

    return data_list

# ==============================================================                     
# DEF : read parameter file                                                          
def read_param_file():                                                               
    
    with open('param_list.txt', "r") as infile:                                      
        for line in infile:                                                          
            param = line.split('#')[0]                                               
            param_nam = param.split('=')[0]                                          
            param_val = param.split('=')[1]                                          
            param_nam = param_nam.strip()                                            
            param_val = param_val.strip()                                             
            if (param_nam == 't_intvl'): t_intvl = int(param_val)                    
            if (param_nam == 'x_limit'): x_limit = float(param_val)                  
    
    return t_intvl, x_limit                                                          

# ==============================================================
# DEF : read raw clock data 
#def read_rawdata():                                                                 
def read_rawdata(t_intvl, x_limit):                                                  

  # read the clockdata file name from 'clock_list.txt'
    in_clocklist = 'clock_file_list.txt'   # list of clock data
    clock_list = []
    with open(in_clocklist, "r") as infile:
        for line in infile:
            if (line[0] != "#"):
                cnum = int(line.split()[0])                                       
                name = line.split()[1]                                             
                wgt = float(line.split()[3])                                      
                flg = int(line.split()[5])                                        
                clock_list.append([name, wgt, flg, cnum])                         

  # If too-many or no reference clock => error message                               
    flg1_clk_count = 0                                                               
    for i in range(len(clock_list)):                                                 
        if (clock_list[i][-2] == 1):                                                 
            flg1_clk_count = flg1_clk_count + 1                                      
    if (flg1_clk_count == 0):                                                        
        print("!!!!! No reference clock was assigned in clock_file_list.txt => ")
        print("      Set the 1st clock as the reference clock.")                     
        clock_list[0][-2] = 1                                                        
        print("      ",clock_list[0],'\n')                                           
    if (flg1_clk_count > 1):                                                         
        print("!!!!! Too many reference clock is assigned in clock_file_list.txt!" + '\n')
        print("      Set the 1st clock as the reference clock.")                     
        clock_list[0][-2] = 1                                                        
        print(clock_list[0],'\n')                                                    
  
  # initialize
    data = []
    t_raw = []
    x_raw = []
    namelist = []
    data_raw = []
    data_raw_list = []

  # read clock data file from ../Data_Clock/****.txt (<== clock_list.txt)
    for line in clock_list:
        clk_name = "../Data_Clock/" + line[0]
        namelist.append(line[0])
        clk_wgt = line[1]
        clk_flg = line[2]
        clk_id  = line[3]                                                         
        data = np.genfromtxt(clk_name, usecols=(0, 1), invalid_raise=False)
        t_raw = [data[n][0] for n in range(len(data))]
        x_raw = [data[n][1] for n in range(len(data))]
      # data_raw[] is [[MJD,,,],[phase,,,]]
        data_raw = np.array([t_raw, x_raw])
      # data_raw_list[] is [["UTCk-CSxx", [[MJD,,,],[phase,,,]], wgt, flg, id], [],,]
        data_raw_list.append([line[0], data_raw, clk_wgt, clk_flg, clk_id])               
        
  # comment to log[]n
    comment = ('**0** Read and plot initial clock data and Allan deviation' 
               '*****' + '\n')
    log_0 = "::::: Input clocklist  : " + in_clocklist + '\n'
    log_1 = "::::: Input clock file : " + str(namelist) + '\n'
    log_2 = "::::: data time interval (t_intval): " + str(t_intvl) + '(sec)\n'                  
    log_3 = "::::: Outlier limit (x_limit))     : " + str(x_limit) + '\n'                       
    log = [comment + log_0 + log_1 + log_2 + log_3]                                  
    return data_raw_list, log, comment

# ==============================================================
# DEF: joint each column of clock data 
def save_jointdata(outfile, data_list, t_title, log):

  # title of the time as "Tau(s)"
    output_line0 = ["   " + str(t_title) + "  "]
    output_line1 = ["  MJDspan   "]                                                

  # ...............................................................                
  # First culumn of all_data[] is values of Tau
  # get the maximum column from the ADEVs of all clocks
    all_data = []
    maxt = []
    for nclk, line in enumerate(data_list):
        if (len(line[1][0]) > len(maxt)):
            maxt = line[1][0]    
    all_data = [maxt]
  # ...............................................................                
  
  # read Adev & clock data from data_list[]
    for nclk, line in enumerate(data_list):
        name = line[0]                # name = "Tau(s)" or "UTCk-CSxx"
        #t = line[1][0]                # t is list of t_title(MJD or Tau)
        x = line[1][1]                # x is list of UTCk_CS data or ADEV
        MJD0 = line[2][0]                               
        MJD1 = line[2][1]                              
        output_line0.append(name + " ")
        output_line1.append(str(MJD0) + "-" + str(MJD1) + "    ") 
      # np.transpose needs the same length of all columns
        lenx = len(maxt) - len(x)
        if (lenx > 0):
          for i in range(lenx):                   
            x = np.append(x,np.nan)
      # all_data[] is [[tau], [Adev of clock1], [Adev of clock2],,]
        col_clock = []
        for i in range(len(maxt)):
            col_clock = np.hstack((col_clock,x[i])) # col_clock is a column of each adev
        all_data.append(col_clock)                  # all_data is bunch of multi-columns

  # save the joint data file ==> "temp2.dat"
    n_lin = len(col_clock)
    n_col = len(all_data)
    output = np.transpose(all_data)
    np.savetxt("temp2.dat", output.reshape((n_lin,n_col)),
               fmt='  '.join(['%9d'] + ['%+.6e']*(n_col-1)))

  # save the joint data file with title ==> "all_adev.dat"
    with open(outfile, "w") as outf:
        outf.writelines(output_line0 + ['\n'] + output_line1 + ['\n'])                      
        with open("temp2.dat", "r") as inf:
            for line in inf:
                outf.writelines(line)

    return outfile, log

# ==============================================================
# DEF : save data & log
def save_logfigdat(alldata, data_adev_list, detdata, wgtfile, log, fig, w_fig):                   

  # make the directory
    dt_str, dt_com, path = set_dirpath("./TA_RESULT/")                                
    try:
        os.makedirs(path)
    except FileExistsError:
        print("Directry already exists.")
    
  # set the output filename
    output_log = path + "log_" + dt_str + ".txt"
    output_fig = path + "fig_" + dt_str + ".pdf"
    output_adv = path + "all_adev_" + dt_str + ".dat"
    output_dat = path + "" + alldata.split('.')[0] + "_" + dt_str + ".dat" 
    output_wgt = path + wgtfile.split('.')[0] + "_" + dt_str + ".dat"                    
    output_wfg = path + "wgt_fig_" + dt_str + ".pdf"                                     

  # message of the output path and filename
    print("")
    print("::::: save logfile = ", output_log)
    print("::::: save figfile = ", output_fig)
    print("::::: save alldata = ", output_dat)
    print("::::: save ADEVdat = ", output_adv)

  # save logfile to './TA_RESULT/log_mmdd_hhmm.txt'
    np.savetxt(output_log, np.transpose(log), fmt="%s", 
               header = "Data log at " + dt_com + '\n')

  # save figure to './TA_RESULT/fig_mmdd_hhmm.pdf'
    fig.savefig(output_fig, bbox_inches='tight')
    if (w_fig != ""):                                                                 
        w_fig.savefig(output_wfg, bbox_inches='tight')                               

  # copy the alldata(='all_*.dat') to './TA_RESULT/all_*_mmdd_hhmm.dat' 
    shutil.copyfile(alldata, output_dat)

  # save the adev result in 'all_adev.dat'
    outfile, log = save_jointdata("all_adev.dat", data_adev_list, "Tau(s)", log)
  # copy the outfile(='all_adev.dat') to './TA_RESULT/all_adev_mmdd_hhmm.dat' 
    shutil.copyfile(outfile, output_adv)

  # save each de-trend clockdata
    data = []
    if (len(detdata) != 0):
        for nclk, line in enumerate(detdata):
            name = line[0]
            name = path + "/" + name + "_" + dt_str
            data = [line[1], line[2], line[3], line[4]]
            title = "MJD     raw data        rate         detrended"
          # save each de-trended data to './TA_RESULT/UTCk-CS**_det_mmdd_hhmm.dat'
            np.savetxt(name + ".dat", np.transpose(data), header = title,
                       fmt='  '.join(['%5d'] + ['%+.6e']*3))

  # copy the wgtfile(='all_weight.dat') to './TA_RESULT/all_weight.dat' 
    if (wgtfile != ""):                                              # yuk0250521
        print("::::: save weights = ", output_wgt)                                    
        shutil.copyfile(wgtfile, output_wgt)                                          
        
    return log, fig

# ==============================================================                      
# DEF : merge weight                                                                  
def save_weight(all_wgt_list):                                                        

    clk = all_wgt_list[0][1]                                                          
    all_w_data = []                                                                   
    for t in range(len(all_wgt_list)):                                                
        mjd = all_wgt_list[t][0]                                                      
        wgt = all_wgt_list[t][2]                                                      
        all_w_data.append([mjd,wgt])                                                  
        
    wfile = 'all_weight.dat'                                                          
    with open(wfile, "w") as f:                                                       
        line0 = "  MJD   " + " ".join(map(str, clk)) + "\n"                                                                                   
        f.write(line0)                                                                
        for xx in all_w_data:                                                         
            xx_wgt = [f"{x:.2f}" for x in xx[1]]                                      
            line = f"{xx[0]} " + " ".join(map(str, xx_wgt)) + "\n"                    
            f.write(line)        
                                                      
    return wfile

# ==============================================================
# DEF : make the data directry                                                            
def set_dirpath(dirname):
    
  # get now-time
    dt_now = str(datetime.datetime.now())
    dt_mmdd = dt_now[5:7] + dt_now[8:10]
    dt_hhmmss = dt_now[11:13] + dt_now[14:16] + dt_now[17:19]
    dt_str = str(dt_mmdd)+"_"+str(dt_hhmmss)
    dt_com = (dt_now[0:4] + "." + dt_now[5:7] + "." + dt_now[8:10] 
              + "_" + dt_now[11:13] + "h" + dt_now[14:16] + "m"
              + "_" + dt_now[17:19] + "s")
    path = dirname + dt_str + "/"                                                    
    
    return dt_str, dt_com, path

# ==============================================================
# main program
# ==============================================================
# select the process menu

menu_n = 0                                                                           
data_list = []                                                                       
clk_initparam = []                                                                   

# read parameter file                                                                
t_intvl, x_limit = read_param_file()                                                 

# read initial clock data
print("\n" + "..... Reading initial data of 'UTC(k) - Clock' ")                   
data_list, log, comment = read_rawdata(t_intvl, x_limit)                             

# Search maximum MJD span (instead of common span)
MJDall = MJD_max(data_list)

# joint clock data and output ==> "all_raw_clk.dat"
print("..... Formatting the initial data for processing ")
outfile_raw = "all_raw_clk.dat"
outfile_raw, log = joint_rawdata(outfile_raw, MJDall, data_list, log)                  

# Show the initial status of all clocks
print("\n************ Initial status of all clocks ************")
menu_n = 0
outfile = outfile_raw                                                                  
adev_plot(outfile, "","",[], MJDall, log, comment, menu_n, x_limit)                     

# select data processing while int(menu_n) < 10:
chk0 = False                                                                       
while (chk0 == False):                                                            

    print('')
    print('=======================================')
    print('      Menu of the TA processing        ')
    print('---------------------------------------')
    print('  0.  Plot the "UTC(k)-Clock" data'     )                
    print('  1.  Assign the MJD range'            )
    print('  2.  De-trend "UTC(k)-Clock" data'     )                              
    print('  3.  Make TA with simple averaging'    )
    print(' 10.  Exit'                             )
    print('=======================================')
    
    menu_n = input('\n????? Input menu number ? >> ')
    if (menu_n not in ['0','1','2','3','10']):                                     
        print ("\n!!!!! Input is wrong, Try again! ")                             

    if menu_n == '0':
        data_list = []
        data_list, log, comment = read_rawdata(t_intvl, x_limit)                     
        MJDall = MJD_max(data_list)
        outfile, log = joint_rawdata(outfile_raw, MJDall, data_list, log)
        adev_plot(outfile, "","",[], MJDall, log, comment, menu_n, x_limit)                  

    if menu_n == '1':
        MJDall, log = MJD_limit(MJDall, log)                                      
        outfile, log = joint_rawdata(outfile, MJDall, data_list, log)
        adev_plot(outfile, "","",[], MJDall, log, comment, menu_n, x_limit)                  

    if menu_n == '2':
        outfile, data_list, det_clk, log, comment = do_detrend(outfile, 
                                                    MJDall, log, x_limit)
        adev_plot(outfile, det_clk,"",[], MJDall, log, comment, menu_n, x_limit)                  

    if menu_n == '3':
        outfile_p, log, comment, wgtfile, all_wgt_list = do_average(outfile, MJDall, log, x_limit) 
        adev_plot(outfile_p, "",wgtfile, all_wgt_list, MJDall, log, comment, menu_n, x_limit) 

    if menu_n == '10':
        chk0 == 0                                                                 
        break
