#! /usr/local/bin/python3
from __future__ import division
from future.builtins import input
from datetime import datetime
from lxml import etree
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
from matplotlib import cm
import pandas as pd
import numpy as np
import seaborn as sns
from scipy import stats
# from skimage import exposure, feature, io, transform, filters
import glob, os, json, sys, math, warnings
import ebovchan as ebc
pd.set_option('display.width', 1000)
pd.options.display.max_rows = 999
sns.set()
#*********************************************************************************************#
#
#           FUNCTIONS
#
#*********************************************************************************************#

#*********************************************************************************************#
def nano_csv_reader(chip_name, spot_data, csv_list):
    min_corr = input("\nWhat is the correlation cutoff for particle count?"+
                     " (choose value between 0.5 and 1)\t")
    if min_corr == "": min_corr = 0.75
    min_corr = float(min_corr)
    contrast_window = input("\nEnter the minimum and maximum percent contrast values," +
                            " separated by a comma (for VSV, 0-6% works well)\t")
    assert isinstance(contrast_window, str)
    contrast_window = contrast_window.split(",")
    cont_0 = (float(contrast_window[0])/100)+1
    cont_1 = (float(contrast_window[1])/100)+1
    min_corr_str = str("%.2F" % min_corr)
    particles_list = []
    particle_dict = {}
    nano_csv_list = [csvfile for csvfile in csv_list if csvfile.split(".")[-2].isdigit()]
    for csvfile in nano_csv_list: ##This pulls particle data from the CSVs generated by nanoViewer
        csv_data = pd.read_table(
                             csvfile, sep = ',',
                             error_bad_lines = False, usecols = [1,2,3,4,5],
                             names = ("contrast", "correlation", "x", "y", "slice")
                             )
        filtered = csv_data[(csv_data['contrast'] <= cont_1)
                    & (csv_data['contrast'] > cont_0)
                    & (csv_data['correlation'] >= min_corr)][['contrast','correlation']]
        particles = len(filtered)
        csv_id = csvfile.split(".")[1] + "." + csvfile.split(".")[2]
        particle_dict[csv_id] = list(round((filtered.contrast - 1) * 100, 4))
        particles_list.append(particles)
        print('File scanned: '+ csvfile + '; Particles counted: ' + str(particles))
        particle_count_col = str('particle_count_'+ min_corr_str
                           + '_' + contrast_window[0]
                           + '_' + contrast_window[1]+ '_')
    spot_data[particle_count_col] = particles_list
    #for row in spot_data.iterrows():
    filtered_density = spot_data[particle_count_col] / spot_data.area * 0.001
    spot_data = pd.concat([spot_data, filtered_density.rename('kparticle_density')], axis = 1)
    dict_file = pd.io.json.dumps(particle_dict)
    with open('../virago_output/' + chip_name + '/' + chip_name
              + '_particle_dict_' + min_corr_str + 'corr.txt', 'w') as f:
              f.write(dict_file)
    print("Particle dictionary file generated")

    return min_corr, spot_data, particle_dict, contrast_window
#*********************************************************************************************#
def csv_fixer(csv_list, txtcheck):
    csvcheck = set([(".".join(file[:-1])+'.csv') for file in txtcheck if (len(file) > 2)
                    and (file[2].isdigit())])
    miss_csv = list(csvcheck.difference(csv_list))
    blank_csv = np.array([[1,0,0,0,0,1]])
    for csv in miss_csv:
        # if vir_toggle is True:
        #     vcount_str = '.0.vcount'
        #     csv = csv.split(".")
        #     csv = ".".join(csv[:3]) + vcount_str + '.csv'
        print("Missing particle data... Generating blank: ", csv)
        np.savetxt(csv, blank_csv, delimiter=",", fmt=('%i','%i','%i','%i','%i','%i'))
    csv_list = sorted(glob.glob('*.csv'))
    return csv_list
#*********************************************************************************************#
# def density_normalizer(spot_data, pass_counter, spot_list):
#     """Particle count normalizer so pass 1 = 0 particle density"""
#     normalized_density = []
#     for spot in spot_list:
#         normspot = [val[0] for val in spot_data.spot_number.iteritems() if int(val[1]) == spot]
#         x = 0
#         while x < pass_counter - 1:
#             if all(np.isnan(spot_data.kparticle_density[normspot])):
#                 print("Missing Value")
#                 normalized_density = normalized_density + ([np.nan] * pass_counter)
#                 break
#             elif np.isnan(spot_data.kparticle_density[normspot[x]]):
#                     print("Missing value for Pass " + str(x + 1))
#                     normalized_density = normalized_density + [np.nan]
#                     x += 1
#             else:
#                 norm_d = [
#                           (spot_data.kparticle_density[normspot[scan]]
#                            - spot_data.kparticle_density[normspot[x]])
#                           for scan in np.arange(x,pass_counter,1)
#                          ]
#                 #print(norm_d)
#                 normalized_density = normalized_density + norm_d
#                 break
#     #print(normalized_density)
#     return normalized_density
# #*********************************************************************************************#
#*********************************************************************************************#
# def chip_file_reader(xml_file): ##XML file reader, reads the chip file for IRIS experiment
#     xml_raw = etree.iterparse(xml_file)
#     chip_dict = {}
#     chip_file = []
#     for action, elem in xml_raw:
#         if not elem.text:
#             text = "None"
#         else:
#             text = elem.text
#         #print(elem.tag + " => " + text)
#         chip_dict[elem.tag] = text
#         if elem.tag == "spot":
#             chip_file.append(chip_dict)
#             chip_dict = {}
#     return chip_file
#*********************************************************************************************#
#
#    CODE BEGINS HERE
#
#*********************************************************************************************#
##Point to the correct directory
retval = os.getcwd()
print("\nCurrent working directory is:\n %s" % retval)
iris_path = input("\nPlease type in the path to the folder that contains the IRIS data:\n")
os.chdir(iris_path.strip('"'))

txt_list = sorted(glob.glob('*.txt'))
pgm_list = sorted(glob.glob('*.pgm'))
csv_list = sorted(glob.glob('*.csv'))
xml_list = sorted(glob.glob('*/*.xml'))
if not xml_list: xml_list = sorted(glob.glob('../*/*.xml'))
chip_name = pgm_list[0].split(".")[0]

mirror_file = str(glob.glob('*000.pgm')).strip("'[]'")
if mirror_file:
    pgm_list.remove(mirror_file)
    # mirror = io.imread(mirror_file)
    print("Mirror file detected")
    mirror_toggle = True
else: print("Mirror file absent"); mirror_toggle = False

# zslice_count = max([int(pgmfile.split(".")[3]) for pgmfile in pgm_list])
txtcheck = [file.split(".") for file in txt_list]
iris_txt = [".".join(file) for file in txtcheck if (len(file) >= 3) and (file[2].isalpha())]
nv_txt = [".".join(file) for file in txtcheck if (len(file) > 3) and (file[2].isdigit())]

### Important Value
if nv_txt: pass_counter = max([int(file[2]) for file in txtcheck if (len(file) > 3)])
###

xml_file = [file for file in xml_list if chip_name in file]
chip_file = ebc.chip_file_reader(xml_file[0])
intro = chip_file[0]
#*********************************************************************************************#

mAb_dict, mAb_dict_rev = ebc.dejargonifier(chip_file)

#*********************************************************************************************#
spot_counter = len([key for key in mAb_dict])##Important
if sys.platform == 'win32': folder_name = iris_path.split("\\")[-1]
elif sys.platform == 'darwin': folder_name = iris_path.split("/")[-1]
else: folder_name = ''
if len(folder_name.split("_")) == 2:
    sample_name = folder_name.split("_")[-1]
else:
    sample_name = input("\nPlease enter a sample descriptor (e.g. VSV-MARV@1E6 PFU/mL)\n")

if not os.path.exists('../virago_output/'+ chip_name): os.makedirs('../virago_output/' + chip_name)
##Varibles
averaged_data = []
normalized_density = ([])
spot_labels = []

#*********************************************************************************************#
# Text file Parser
#*********************************************************************************************#

spot_data_nv = pd.DataFrame([])
spot_list = [int(file[1]) for file in txtcheck if (len(file) > 2) and (file[2].isalpha())]
scanned_spots = set(np.arange(1,spot_counter+1,1))
missing_spots = scanned_spots.difference(spot_list)
miss_txt = 1
for txtfile in iris_txt:
    if miss_txt in missing_spots:
        print('Missing text file:  ' + str(miss_txt))
        miss_list = pd.Series(list(str(miss_txt))*pass_counter)
        blanks = pd.DataFrame(np.zeros((pass_counter,3)))
        blanks.insert(0,'spot_number', miss_list)
        miss_txt += 1

    txtdata = pd.read_table(txtfile, sep = ':', error_bad_lines = False,
                            header = None, index_col = 0, usecols = [0, 1])
    pass_labels = [
                    row for row in txtdata.index
                    if row.startswith('pass_time')
                    ]
    if not nv_txt: pass_counter = int(len(pass_labels)) ##If nanoViewer hasn't run on data

    spot_idxs = pd.Series(list(txtdata.loc['spot_index']) * pass_counter)
    pass_list = pd.Series(np.arange(1,pass_counter + 1))
    spot_types = pd.Series(list([mAb_dict[int(txtfile.split(".")[1])]]) * pass_counter)

    times_s = pd.Series(txtdata.loc[pass_labels].values.flatten().astype(np.float))
    times_min = round(times_s / 60,2)
    pass_diff = pass_counter - len(pass_labels)
    if pass_diff > 0:
        times_min = times_min.append(pd.Series(np.zeros(pass_diff)), ignore_index = True)
    print('File scanned:  ' + txtfile)
    miss_txt += 1
    spot_data_solo = pd.concat([spot_idxs.rename('spot_number').astype(int),
                                pass_list.rename('scan_number').astype(int),
                                times_min.rename('scan_time'),
                                spot_types.rename('spot_type')], axis = 1)
    spot_data_nv = spot_data_nv.append(spot_data_solo, ignore_index = True)

spot_data_vir = spot_data_nv.copy()

area_col = []
for txtfile in nv_txt:
    if int(txtfile.split(".")[1]) in missing_spots:
        print("Did not scan " + txtfile + "; data missing")
    else:
        txtdata = pd.read_table(txtfile, sep = ':', error_bad_lines = False,
                                header = None, index_col = 0, usecols = [0, 1])
        area = float(txtdata.loc['area'])
        area_col.append(area)
        print('File scanned:  ' + txtfile)
area_col = pd.Series(area_col, name = 'area')

spot_data_nv['area'] = area_col
spot_data_nv.scan_time.replace(0, np.nan, inplace = True)

spot_labels = [[val]*(pass_counter) for val in mAb_dict.values()]

# [spot_set.add(val) for val in mAb_dict.values()]
# spot_set = list(spot_set)
spot_set = []
for val in mAb_dict.values():
    if val not in spot_set: spot_set.append(val)

#*********************************************************************************************#
# CSV Reader
#*********************************************************************************************#


csv_list = csv_fixer(csv_list, txtcheck)

olddata_toggle = 'no'
if os.path.exists('../virago_output/' + chip_name + '/' + chip_name + '_spot_data_nv.csv'):
    olddata_toggle = input("\nPre-existing data exists. Do you want use this data? (y/n)\n"
        "NOTE: Selecting 'no' will overwrite old data)\t")
    assert isinstance(olddata_toggle, str)
    if olddata_toggle.lower() in ('yes', 'y'):
        ##This reads in the pre-existing data files and stores the needed variables
        print("\nUsing pre-existing data...\n")
        preexist_csv = '../virago_output/' + chip_name + '/' + chip_name + '_spot_data_nv.csv'
        spot_data_nv = pd.read_table(preexist_csv, sep = ',', error_bad_lines = False)
        count_info = [col for col in spot_data_nv.columns if col.startswith('particle_count')]
        count_info = str(count_info).split("_")
        min_corr = float(count_info[2])
        contrast_window = count_info[3:5]
        print(str(contrast_window[0])+"-"+contrast_window[1])
        min_corr_str = str("%.2F" % min_corr)
        print("Correlation value is: "+ min_corr_str +" or greater\n")
        with open('../virago_output/' + chip_name + '/'
            + chip_name + '_particle_dict_'
            + min_corr_str + 'corr.txt', 'r') as infile:
            particle_dict = json.load(infile)
    elif olddata_toggle.lower() in ('no', 'n'):
        print("")
        min_corr, spot_data_nv, particle_dict, contrast_window = ebc.nano_csv_reader(chip_name, spot_data_nv, csv_list)
else:
    min_corr, spot_data_nv, particle_dict, contrast_window = ebc.nano_csv_reader(chip_name,
                                                                                 spot_data_nv,
                                                                                 csv_list)

min_corr_str = str("%.2F" % min_corr)
if not os.path.exists('../virago_output/'+ chip_name + '/vcounts'):
    os.makedirs('../virago_output/' + chip_name + '/vcounts')
os.chdir('../virago_output/'+ chip_name + '/vcounts')
# vir_csv_list = sorted(glob.glob(chip_name +'*.vcount.csv'))
# if pgm_toggle is True:
#     vir_csv_list = [".".join(csv.split(".")[:3])+'.csv' for csv in vir_csv_list]
#     vir_csv_list = csv_fixer(vir_csv_list,txtcheck, vir_toggle = True)
#os.chdir(iris_path.strip('"'))

# if len(vir_csv_list) == (len(iris_txt) * pass_counter):
#     particle_count_vir, contrast_window, particle_dict = virago_csv_reader(chip_name, vir_csv_list, vir_toggle = True)
#
#     area_list = np.array([(float(csvfile.split(".")[-3]) / 1e6) for csvfile in vir_csv_list])
#     spot_data_vir['area_sqmm'] = area_list
#
#     particle_count_col = str('particle_count_0'
#                              + '_' + contrast_window[0]
#                              + '_' + contrast_window[1] + '_')
#     spot_data_vir[particle_count_col] = particle_count_vir
#
#     kparticle_density = np.round(np.array(particle_count_vir) / area_list * 0.001,3)
#     spot_data_vir['kparticle_density'] = kparticle_density
# # if (olddata_toggle.lower() not in ('yes', 'y')):
#     nv_vir_toggle = input("Would you like to use nanoViewer or VIRAGO data? (type N or V)\n")
#     assert isinstance(nv_vir_toggle, str)
#     if nv_vir_toggle.lower() in ('n','nanoViewer'):
#         print("Using nanoViewer data...")
#         spot_data = spot_data_nv
#     else:
#         print("Using VIRAGO data...")
#         spot_data = spot_data_vir
#         min_corr_str = ""
# else:
spot_data = spot_data_nv
os.chdir(iris_path.strip('"'))

#####################################################################
# Histogram generator
#####################################################################
#--------------------------------------------------------------------
# spots_to_hist = input("Which spots would you like to generate histograms for?\t")
# hist_norm = False
# # hist_norm_toggle = input("Do you want to normalize the counts to a percentage? (y/[n])")
# # if hist_norm_toggle.lower() in ('y','yes'): hist_norm = True
# spots_to_hist = spots_to_hist.split(",")
# print(spots_to_hist)
# spots_to_hist2 = [int(spot) for spot in spots_to_hist]
# #cont_0 = float(contrast_window[0])
# cont_1 = float(contrast_window[1])
# for spot in spots_to_hist:
#     hist_dict = {}
#     for key in sorted(particle_dict.keys()):
#         hist_spot = int(key.split(".")[0])
#
#         if hist_spot == int(spot): hist_dict[key] = particle_dict[key]
#     nrows = 2
#     ncols = math.ceil(pass_counter / 2)
#     fig = plt.figure()
#     plt.axis('off')
#
#     if hist_norm_toggle == False:
#         fig.text(0.06,0.6,"Particle Counts " + min_corr_str, fontsize = 10, rotation = 'vertical')
#     elif hist_norm_toggle == True:
#         fig.text(0.06,0.6,"Particle Frequency" + min_corr_str, fontsize = 10, rotation = 'vertical')
#     fig.text(.4,0.04,"Percent Contrast", fontsize = 10)
#
#     for key in sorted(hist_dict.keys()):
#
#         hist_pass = int(key.split(".")[1])
#         sbplt = fig.add_subplot(nrows,ncols,hist_pass)
#         (hist_vals, bins, patches) = sbplt.hist([hist_dict[key]],
#                                      100, range = [0,cont_1], color = ['#0088FF'],
#                                      rwidth = 1, alpha = 0.75, normed = hist_norm)
#         plt.xticks(np.arange(0, cont_1+1,2), size = 5)
#         if max(hist_vals) <= 50: grads = 5
#         elif max(hist_vals) <= 50: grads = 10
#         else: grads = 25
#         if hist_norm_toggle == False: plt.yticks(np.arange(0, (max(hist_vals)) + 10 , grads), size = 5)
#
#         plt.title("Pass "+ str(hist_pass), size = 5)
#         plt.grid(True, alpha = 0.5)
#
#         if hist_norm_toggle == False: sbplt.set_ylim(0, (max(hist_vals)) + 10)
#     plot_title = str("Particle Contrast Distribution of " + chip_name + " "
#                     + spot_labels[int(spot)-1][0]
#                     + " Spot " + spot)
#     plt.suptitle(plot_title, size = 12)
#     plt.subplots_adjust(wspace = 0.25)
#
#     plt.savefig('../virago_output/' + chip_name + '/' + chip_name + '_spot-' + spot
#                 +  '_histo.png', bbox_inches = 'tight', dpi = 300)
#     print('File generated: ' +  chip_name + '_spot-' + spot + '_histo.png')
#     plt.close()

# hist_dict2 = {}
# for key in particle_dict.keys():
#     if int(key.split(".")[0]) in spots_to_hist2:
#         hist_dict2[key] = particle_dict[key]

vhf_colormap = ('#e41a1c','#377eb8','#4daf4a',
            '#984ea3','#ff7f00','#ffff33',
            '#a65628','#f781bf','gray','black')

histogram_df = ebc.histogrammer(particle_dict, spot_counter, baselined = True)

mean_histogram_df = ebc.histogram_averager(histogram_df, mAb_dict_rev, pass_counter)

ebc.combo_histogram_fig(mean_histogram_fg, colormap = vhf_colormap, corr = min_corr_str)

#     # for scan in range(1,pass_counter+1):
#     # int(key.split(".")[1]) == scan:
#     sns.distplot(particle_dict[key], bins = int(cont_1 * 10),
#                                  color = 'blue', kde = False, norm_hist = False,
#                                  hist_kws = {
#                                               "histtype":"step",
#                                               "linewidth":1,
#                                               "alpha":1})
#     c += 1
#     plot_title2 = str("Particle Contrast Distribution of " + chip_name + str(key))
#     plt.title(plot_title2)
#
#     plt.ylabel('Particle Count\n'+ 'Correlation Value >=' + min_corr_str, size = 8)
#     plt.yticks(range(0,101,10), size = 8)
#     plt.xlabel("Percent Contrast", size = 8)
#     plt.xticks(np.arange(0,10.5,.5), size = 8)
#     plt.axvline(x = 1, linestyle = '--', color = 'k', linewidth = 1)
#     plt.axvline(x = 6, linestyle = '--', color = 'k', linewidth = 1)

#
#
#     print('File generated: ' +  chip_name + '_histo_sns_' + str(key) + '.png')
#     plt.clf(); plt.close('all')

# for scan in hist_scans:
#     for n, key in enumerate(hist_dict2):
#         sns.distplot(hist_dict2[key], bins = int(cont_1 * 10), axlabel = "Percent Contrast",
#                                      color = vhf_colormap[n], kde = False, norm_hist = True,
#                                      hist_kws = {
#                                                   "histtype":"step",
#                                                   "linewidth":1,
#                                                   "alpha":1})
#         plt.title(plot_title)
#

#
#         print('File generated: ' +  chip_name + '_spot-' + spot + '_histo_sns.png')
#         plt.close()

#*********************************************************************************************#
# Particle count normalizer so pass 1 = 0 particle density
#*********************************************************************************************#

normalized_density = ebc.density_normalizer(spot_data, pass_counter, spot_list)
len_diff = len(spot_data) - len(normalized_density)
if len_diff != 0:
    normalized_density = np.append(np.asarray(normalized_density),np.full(len_diff, np.nan))
spot_data['normalized_density'] = normalized_density
print(spot_data)
#*********************************************************************************************#
def spot_remover(spot_data):
    excise_toggle = input("Would you like to remove any spots from the dataset? (y/(n))\t")
    assert isinstance(excise_toggle, str)
    if exise_toggle.lower() in ('y','yes'):
        spots_to_excise = input("Which spots? (Separate all spot numbers by a comma)\t")
        spots_to_excise = spots_to_excise.split(",")


##IN PROGRESS



# -------------------------------------------------------------------
#####################################################################
# This gets the average values and standard deviations for each spot type
#####################################################################
#--------------------------------------------------------------------

scan_series = spot_data.scan_number
for k, val in enumerate(spot_set):
    x = 1
    for val in pass_labels:
        data_slice = spot_data[['spot_type', 'scan_time', 'kparticle_density',
                                'normalized_density']][(scan_series == x)
                                & (spot_data['spot_type'] == spot_set[k])]
        scan_time_mean = round(data_slice['scan_time'].mean(),2)
        filt_density_mean = round(data_slice['kparticle_density'].mean(),2)
        filt_density_std = round(np.std(data_slice['kparticle_density']),2)
        norm_density_mean = round(data_slice['normalized_density'].mean(),2)
        norm_density_std = round(np.std(data_slice['normalized_density']),4)
        avg_data = (spot_set[k],
                    spot_data.loc[x - 1,'scan_number'],
                    scan_time_mean,
                    filt_density_mean,
                    filt_density_std,
                    norm_density_mean,
                    norm_density_std)
        averaged_data.append(avg_data)
        x += 1
averaged_data = pd.DataFrame(averaged_data,
                             columns =  ['spot_type', 'scan_number', 'avg_scan_time',
                                         'avg_kparticle_density', 'kparticle_density_std',
                                         'avg_normalized_density', 'normalized_density_std']
                             )
# -------------------------------------------------------------------
#####################################################################
# Asks whether the time series should be set such that Time 0 == 0 particle density
#####################################################################
#--------------------------------------------------------------------
baseline_toggle = input("Do you want the time series chart normalized to baseline? ([y]/n)\t")
assert isinstance(baseline_toggle, str)
if baseline_toggle.lower() in ('no', 'n'):
    filt_toggle = 'kparticle_density'
    avg_filt_toggle = 'avg_kparticle_density'
    stdev_filt_toggle = 'kparticle_density_std'
else:
    filt_toggle = 'normalized_density'
    avg_filt_toggle = 'avg_normalized_density'
    stdev_filt_toggle = 'normalized_density_std'
    print("Normalizing...")

# -------------------------------------------------------------------
#####################################################################
# Time Series Generator
#####################################################################
#--------------------------------------------------------------------

fig = plt.figure(figsize = (8,6))
ax1 = fig.add_subplot(111)
n,c = 1,0
for key in mAb_dict.keys():
    time_x = spot_data[spot_data['spot_number'] == key]['scan_time'].reset_index(drop = True)
    density_y = spot_data[spot_data['spot_number'] == key][filt_toggle].reset_index(drop = True)
    while n > 1:
        if mAb_dict[n-1] != mAb_dict[n]:
            c += 1
            break
        else:
            break
    ax1.plot(time_x, density_y, marker = '+', linewidth = 1,
                 color = vhf_colormap[c], alpha = 0.4, label = '_nolegend_')
    n += 1
ax2 = fig.add_subplot(111)

for n, spot in enumerate(spot_set):
    avg_data = averaged_data[averaged_data['spot_type'].str.contains(spot)]
    avg_time_x = avg_data['avg_scan_time']
    avg_density_y = avg_data[avg_filt_toggle]
    errorbar_y = avg_data[stdev_filt_toggle]
    ax2.errorbar(avg_time_x, avg_density_y,
                    yerr = errorbar_y, marker = 'o', label = spot_set[n],
                    linewidth = 2, elinewidth = 1, capsize = 3,
                    color = vhf_colormap[n], alpha = 0.9, aa = True)

#min_corr_str = str(0)
ax2.legend(loc = 'upper left', fontsize = 8, ncol = 1)
plt.xlabel("Time (min)", color = 'gray')
plt.ylabel('Particle Density (kparticles/sq. mm)\n'+ contrast_window[0]+'-'+contrast_window[1]
            + '% Contrast, Correlation Value >=' + min_corr_str, color = 'gray')
plt.xticks(np.arange(0, max(spot_data.scan_time) + 1, 5), color = 'gray')
plt.yticks(color = 'gray')
plt.title(chip_name + ' Time Series of ' + sample_name)

plt.axhline(linestyle = '--', color = 'gray')
plot_name = chip_name + '_timeseries_' + min_corr_str + 'corr.png'

plt.savefig('../virago_output/' + chip_name + '/' +  plot_name,
            bbox_inches = 'tight', pad_inches = 0.1, dpi = 300)
print('File generated: ' + plot_name)
csv_spot_data = str('../virago_output/' + chip_name + '/' + chip_name + '_spot_data.csv')
spot_data.to_csv(csv_spot_data, sep = ',')
#plt.show()
plt.clf(); plt.close('all')
print('File generated: '+ csv_spot_data)
# -------------------------------------------------------------------
#####################################################################
# Bar Plot Generator
#####################################################################
#--------------------------------------------------------------------
first_scan = min(scan_series)
last_scan = max(scan_series)
baseline = (spot_data[scan_series == first_scan][['spot_type', 'kparticle_density']]).reset_index(drop = True)
post_scan = pd.Series(spot_data[scan_series == last_scan]['kparticle_density'],
                      name = 'post_scan').reset_index(drop = True)
difference = pd.Series(spot_data[scan_series == last_scan]['normalized_density'],
                    name = 'difference').reset_index(drop = True)
barplot_data = pd.concat([baseline, post_scan, difference], axis = 1)
csv_bardata = str('../virago_output/' + chip_name + '/' + chip_name + '_barplot_data.csv')
barplot_data.to_csv(csv_bardata, sep =',')
#barplot_data.kparticle_density = barplot_data.kparticle_density * -1
baseline_avg, post_scan_avg, baseline_std, post_scan_std, diff_avg, diff_std = [],[],[],[],[],[]
for spot in spot_set:
    avg_data = barplot_data[barplot_data['spot_type'].str.contains(spot)]
    baseline_avg.append(np.mean(avg_data.kparticle_density))
    baseline_std.append(np.std(avg_data.kparticle_density))

    post_scan_avg.append(np.mean(avg_data.post_scan))
    post_scan_std.append(np.std(avg_data.post_scan))

    diff_avg.append(np.mean(avg_data.difference))
    diff_std.append(np.std(avg_data.difference))
fig,axes = plt.subplots(nrows = 1, ncols = 1, figsize = (5,4), sharey = True)
fig.subplots_adjust(left=0.08, right=0.98, wspace=0)
plt.suptitle("Experiment "+ chip_name + "- Final Scan difference versus Inital Scan\n"
             + "Sample Conditions: " + sample_name, size = 12)
#
# bar1 = axes[0].bar(np.arange(len(spot_set)), baseline_avg, width = 0.45, color = vhf_colormap[0],
#                    tick_label = spot_set, yerr = baseline_std,capsize = 4, alpha = 0.75)
#
# bar2 = axes[0].bar(np.arange(len(spot_set)) + 0.45, post_scan_avg, width = 0.45,
#                    color = vhf_colormap[1], tick_label = spot_set, yerr = post_scan_std,
#                    capsize = 4, alpha = 0.75)

axes.set_ylabel('Particle Density (kparticles/sq. mm)\n' + contrast_window[0]
            +'-'+contrast_window[1] + '% Contrast, Correlation Value >='
            + min_corr_str, color = 'k', size = 8)
bar3 = axes.bar(np.arange(len(spot_set)) + (0.45/2), diff_avg, width = 0.5,
                   color = vhf_colormap[3],tick_label = spot_set, yerr = diff_std, capsize = 4)
#plt.xticks(np.arange(len(spot_set)) + (0.45/2), spot_set, rotation = 45, size = 6)
# for ax in axes:
axes.yaxis.grid(True)
axes.set_xticklabels(spot_set, rotation = 45, size = 6)
axes.set_xlabel("Antibody", color = 'k', size = 8)

#plt.figlegend((bar3),'Final Scan Difference',loc = 'upper right', fontsize = 10)
barplot_name = (chip_name + '_' + contrast_window[0] + '-' + contrast_window[1]
                + '_'+  min_corr_str + '_barplot.png')
plt.savefig('../virago_output/' + chip_name + '/' + barplot_name,
            bbox_inches = 'tight', pad_inches = 0.1, dpi = 300)
print('File generated: '+ barplot_name)
#plt.show()
plt.clf(); plt.close('all')
