""" this script takes extracted stress relaxation data and, 1) finds zero force reference location, calculates various moduli and saves them to the output xml file. If multiple files provided, mean and SD for moduli and thickness values are also calculated and stored. data is plotted and saved

# input needed =  stress_relaxation_repeatability.txt :  contains name of text file with extracted stress relaxation data, name of output csv file generated from data extraction file, sample thickness/length, area of cross section depending on the test (thickness * 1 for tension, 19.625 for compression (pi * r * r = 3.14 * 2.5 * 2.5)) : information for atleast one test is expected """


# -----------------------------------------

import sys
import matplotlib.pyplot as plt
import csv
import numpy as np
import statistics
import ast
import operator

def USAGE(argv):
 print
 print ('USAGE: ') + argv[0] + ' <at least one data file required>'
 print

def main(argv):

 if len(argv)<2:
        USAGE(argv)
        sys.exit(1)
        

 input_filename = argv[1]
 infile = open(input_filename, 'r')
 
 def find_nearest(array, value):
        array = np.asarray(array)
        idx = (np.abs(array - value)).argmin()
        return  idx
 
 sample_names = []
 output_file_names = []
 dimensions = []
 area_cross_section = []
 
 for line in infile: # read each line to extract individual test information
     line = line.rstrip('\n')
     line=line.split()
     sample_names.append(line[0])
     output_file_names.append(line[1])
     dimensions.append(line[2])
     area_cross_section.append(line[3])
 

 num = len(sample_names)
 
 output_filename = sample_names[0][0:18]
 input_tags = ""
 for i in range (0,num):
     input_tags = input_tags + sample_names[i][20:21] # grabbing test numbers from input file names 

 time = [[] for x in range (num)]
 disp = [[] for x in range (num)]
 load_raw = [[] for x in range (num)]
 load_filtered = [[] for x in range (num)]
 norm_pnt_sr = [] 
 strain = [[] for x in range (num)]
 stress = [[] for x in range (num)]
 inst_disp = [[] for x in range (num)]
 inst_load = [[] for x in range (num)]
 eqbm_disp = [[] for x in range (num)]
 eqbm_load = [[] for x in range (num)]
 inst_stress = [[] for x in range (num)]
 inst_strain = [[] for x in range (num)]
 eqbm_stress = [[] for x in range (num)]
 eqbm_strain = [[] for x in range (num)]

 
 average_inst_mod = []
 average_eqbm_mod = []
 high_strain_inst_mod = []
 high_strain_eqbm_mod = []
 low_strain_inst_mod = []
 low_strain_eqbm_mod = []

 per_rlx = [[] for x in range (num)]

 
 for i in range(1,num+1):
    # extract data from each file 
    f_name = sample_names[i-1]
    infile_1  = open(f_name,'r')    
    for line in infile_1:
      line = line.rstrip('\n')
      line = line.split()
      line = map(float,line)
      time[i-1].append(line[0]) 
      disp[i-1].append(line[1])
      load_raw[i-1].append(line[2]) 
      load_filtered[i-1].append(line[3])      
    infile.close() 
    
    plt.figure(1)
    plt.plot(disp[i-1], load_filtered[i-1], label = ('test' + str(i)))
    plt.legend(loc=2)
    plt.legend(loc=2)
    plt.xlabel('disp,mm')
    plt.ylabel('load,gf')

    # extract information from output csv files to ditionary
    f_output = output_file_names[i-1]
    with open(f_output, mode='r') as f_output:
       reader = csv.reader(f_output)
       mydict = dict(reader)
    norm_pnt_sr = mydict["Experimental reference location sr, 10g"]

    
    # normalizing disp
    disp[i-1] = [ x - float(norm_pnt_sr) for x in disp[i-1]]
    
    plt.figure(2)
    plt.plot(disp[i-1], load_filtered[i-1], label = ('test' + str(i)))
    plt.legend(loc=2)
    plt.legend(loc=2)
    plt.title('Experimental reference location')
    plt.xlabel('disp,mm')
    plt.ylabel('load,gf')
    
    # pick reference location at actual 10g load
    
    idx = find_nearest(load_filtered[i-1], 10)
    ref_disp = disp[i-1][idx]
    ref_load = load_filtered[i-1][idx]
    disp[i-1] = [x - disp[i-1][idx] for x in disp[i-1]] 
    
    
    plt.figure(3)
    plt.plot(disp[i-1], load_filtered[i-1], label = ('test' + str(i)))
    plt.legend(loc=2)
    plt.legend(loc=2)
    plt.title('Actual 10g location')
    plt.xlabel('disp')
    plt.ylabel('load, g')
      
    # remove all data before the adjusted reference position 
    idx = find_nearest(disp[i-1], 0.0)    
    disp[i-1] = disp[i-1][idx:]
    load_filtered[i-1] = load_filtered[i-1][idx:]
    
    # adjust loads to 0 
    load_filtered[i-1] = [x - load_filtered[i-1][0] for x in load_filtered[i-1]]
    
    plt.figure(4)
    plt.plot(disp[i-1], load_filtered[i-1], label = ('test' + str(i)))
    plt.legend(loc=2)
    plt.legend(loc=2)
    plt.xlabel('disp')
    plt.ylabel('load, g')
    

    
    # calculate stresses and strains
    strain[i-1] = [c/ float(dimensions[i-1]) for c in disp[i-1]]
    stress[i-1] = [c/float(area_cross_section[i-1]) for c in load_filtered[i-1]]
    stress[i-1] = [c*0.0098 for c in stress[i-1]]      # gf to N
  
    plt.figure(5)
    plt.plot(strain[i-1], stress[i-1],'--', label = ('test' + str(i)))
    plt.legend(loc=2)
    plt.legend(loc=2)
    plt.xlabel('strain')
    plt.ylabel('stress, MPa')


    # ----------------------

    # working on instantaneous stress and strain values to get various moduli, % changes and COV over multiple files

    inst_disp_tmp = []
    eqbm_disp_tmp = []
    inst_load_tmp = []
    inst_load_tmp = []
    inst_disp_tmp =  ast.literal_eval(mydict['inst_disp'])  # extract list as string and convert to list from csv
    inst_load_tmp =  ast.literal_eval(mydict['inst_load'])
    eqbm_disp_tmp =  ast.literal_eval(mydict['eqbm_disp'])  # extract list as string and convert to list from csv
    eqbm_load_tmp =  ast.literal_eval(mydict['eqbm_load'])
    

    inst_disp[i-1] = [x - float(norm_pnt_sr) for x in inst_disp_tmp]
    inst_disp[i-1] = [x - ref_disp for x in inst_disp[i-1]]
    inst_load[i-1] = [x - ref_load for x in inst_load_tmp]

    eqbm_disp[i-1] = [x - float(norm_pnt_sr) for x in eqbm_disp_tmp]
    eqbm_disp[i-1] = [x - ref_disp for x in eqbm_disp[i-1]]
    eqbm_load[i-1] = [x - ref_load for x in eqbm_load_tmp]
 
 
    inst_strain[i-1] = [c / float(dimensions[i-1]) for c in inst_disp[i-1]]
    inst_stress[i-1] = [c / float(area_cross_section[i-1]) for c in inst_load[i-1]] 
    inst_stress[i-1] = [c*0.0098 for c in inst_stress[i-1]]      # gf to N   
    inst_strain[i-1].insert(0,0)
    inst_stress[i-1].insert(0,0) 


    eqbm_strain[i-1] = [c/ float(dimensions[i-1]) for c in eqbm_disp[i-1]]
    eqbm_stress[i-1] = [c/float(area_cross_section[i-1]) for c in eqbm_load[i-1]] 
    eqbm_stress[i-1] = [c*0.0098 for c in eqbm_stress[i-1]]      # gf to N   
    eqbm_strain[i-1].insert(0,0)
    eqbm_stress[i-1].insert(0,0)

    # plot inst and eqbm stress-strains and save fig
    plt.figure(5)
    plt.plot(inst_strain[i-1],inst_stress[i-1], 'r*',  ms=8, label = 'instantaneous' if i==0 else "")
    plt.plot(eqbm_strain[i-1],eqbm_stress[i-1], 'bo',  ms=8, label = 'equilibrium' if i==0 else "")
    plt.legend(loc=2)
    plt.legend(loc=2)
    plt.xlabel('strain')
    plt.ylabel('stress, MPa')
    load_disp_Plot = output_filename + '_' + input_tags  + '_stress_strain_repeat'+ '.svg'
    plt.savefig(load_disp_Plot)
    
    
    # calculate moduli for each file and save in a list for repeatability calculations    
    (avg_inst_mod_fit_all,inst_inter) = np.polyfit(inst_strain[i-1],inst_stress[i-1],1)
    (avg_eqbm_mod_fit_all,eqbm_inter) = np.polyfit(eqbm_strain[i-1],eqbm_stress[i-1],1)
    
    (hs_inst_mod,inst_inter) = np.polyfit(inst_strain[i-1][2:],inst_stress[i-1][2:],1)
    (hs_eqbm_mod,eqbm_inter) = np.polyfit(eqbm_strain[i-1][2:],eqbm_stress[i-1][2:],1)
    
    ls_inst_mod = (inst_stress[i-1][1] - inst_stress[i-1][0]) / (inst_strain[i-1][1]-inst_strain[i-1][0])       
    ls_eqbm_mod = (eqbm_stress[i-1][1] - eqbm_stress[i-1][0]) / (eqbm_strain[i-1][1]-eqbm_strain[i-1][0]) 
    
    
    average_inst_mod.append(avg_inst_mod_fit_all)   # SAVE to csv 
    average_eqbm_mod.append(avg_eqbm_mod_fit_all)   # SAVE to csv
    high_strain_inst_mod.append(hs_inst_mod)        # SAVE to csv
    high_strain_eqbm_mod.append(hs_eqbm_mod)        # SAVE to csv
    low_strain_inst_mod.append(ls_inst_mod)        # SAVE to csv
    low_strain_eqbm_mod.append(ls_eqbm_mod)        # SAVE to csv
   
    # percent relaxation at each strain level  
    per_rlx[i-1] = map(operator.sub, inst_stress[i-1], eqbm_stress[i-1])
    per_rlx[i-1][1:] = map(operator.div, per_rlx[i-1][1:], inst_stress[i-1][1:])
    per_rlx[i-1][1:] = [x * 100 for x in per_rlx[i-1][1:]]                      # SAVE to csv
 print 'per relax =', per_rlx

 
 
# calculate and save individual moduli, mean, SD and coeff of var, individual thickness, mean, SD and coeff of var to a  csv file    
 
# mean and SD 
 avg_inst_mod_repeat_mean = statistics.mean(average_inst_mod)   # SAVE to csv
 avg_inst_mod_repeat_SD = statistics.stdev(average_inst_mod)    # SAVE to csv
 
 avg_eqbm_mod_repeat_mean = statistics.mean(average_eqbm_mod)   # SAVE to csv
 avg_eqbm_mod_repeat_SD = statistics.stdev(average_eqbm_mod)    # SAVE to csv
 
 avg_hs_inst_mod_repeat = statistics.mean(high_strain_inst_mod) # SAVE to csv
 SD_hs_inst_mod_repeat = statistics.stdev(high_strain_inst_mod) # SAVE to csv
 
 avg_hs_eqbm_mod_repeat = statistics.mean(high_strain_eqbm_mod) # SAVE to csv
 SD_hs_eqbm_mod_repeat = statistics.stdev(high_strain_eqbm_mod) # SAVE to csv
# 
# avg_ls_inst_mod_repeat = statistics.mean(low_strain_inst_mod) # SAVE to csv
# SD_ls_inst_mod_repeat = statistics.stdev(low_strain_inst_mod) # SAVE to csv
# 
# avg_ls_eqbm_mod_repeat = statistics.mean(low_strain_eqbm_mod) # SAVE to csv
# SD_ls_eqbm_mod_repeat = statistics.stdev(low_strain_eqbm_mod) # SAVE to csv
 
 
 
 dimensions = [float(x) for x in dimensions]
 mean_dimensions = statistics.mean(dimensions)                  # SAVE to csv
 SD_dimensions = statistics.stdev(dimensions)                   # SAVE to csv


# coeff of var 
 cov_avg_inst_mod = (avg_inst_mod_repeat_SD / avg_inst_mod_repeat_mean) * 100   # SAVE to csv
 cov_avg_eqbm_mod = (avg_eqbm_mod_repeat_SD / avg_eqbm_mod_repeat_mean) * 100   # SAVE to csv
 
 cov_avg_hs_inst_mod = (SD_hs_inst_mod_repeat/ avg_hs_inst_mod_repeat) * 100    # SAVE to csv
 cov_avg_hs_eqbm_mod = (SD_hs_eqbm_mod_repeat / avg_hs_eqbm_mod_repeat) * 100   # SAVE to csv
 
# cov_avg_ls_inst_mod = (SD_ls_inst_mod_repeat/ avg_ls_inst_mod_repeat) * 100    # SAVE to csv
# cov_avg_ls_eqbm_mod = (SD_ls_eqbm_mod_repeat / avg_ls_eqbm_mod_repeat) * 100   # SAVE to csv
 
 cov_dimensions = (SD_dimensions / mean_dimensions) * 100       # SAVE to csv


# error in  percent relaxation at each strain level  (redo to loop instead of calulating several times)
 a = []
 b = []
 c = []

 for i in range(0,num):
     a.append(per_rlx[i][1])
     b.append(per_rlx[i][2])
     c.append(per_rlx[i][3])
  
 mean_PR_strain1 = statistics.mean(a)
 SD_PR_strain1 = statistics.stdev(a)
 cov_PR_strain1 = (SD_PR_strain1 / mean_PR_strain1) * 100 # SAVE to csv
 mean_PR_strain2 = statistics.mean(b)
 SD_PR_strain2 = statistics.stdev(b)
 cov_PR_strain2 = (SD_PR_strain2 / mean_PR_strain2) * 100 # SAVE to csv
 mean_PR_strain3 = statistics.mean(c)
 SD_PR_strain3 = statistics.stdev(c)
 cov_PR_strain3 = (SD_PR_strain3 / mean_PR_strain3) * 100 # SAVE to csv
 
# errors in strains at each strain cluster  (redo to loop instead of calulating several times)
 
 a = []
 b = []
 c = []

 for i in range(0,num):
     a.append(inst_strain[i][1])
     b.append(inst_strain[i][2])
     c.append(inst_strain[i][3])
 print 'a=', a
 print 'b=', b
 print 'c=', c
  
  
 mean_strain1 = statistics.mean(a)
 SD_strain1 = statistics.stdev(a)
 cov_strain1 = (SD_strain1 / mean_strain1) * 100    # SAVE to csv
 mean_strain2 = statistics.mean(b)
 SD_strain2 = statistics.stdev(b)
 cov_strain2 = (SD_strain2 / mean_strain2) * 100    # SAVE to csv
 mean_strain3 = statistics.mean(c)
 SD_strain3 = statistics.stdev(c)
 cov_strain3 = (SD_strain3 / mean_strain3) * 100    # SAVE to csv
 
 
# errors in stresses (both inst and eqbm) at each strain cluster 



# write all results to csv 

 output_file = output_filename + '_' + input_tags + '_output_stress_relaxation' # repeatability output file
 out_datafile_id = open(output_file + '.csv', 'w')
 out_writer = csv.writer(out_datafile_id)

 out_writer.writerow(['average_inst_mod_for_each_test', average_inst_mod ])
 out_writer.writerow(['average_eqbm_mod_for_each_test', average_eqbm_mod ])
 out_writer.writerow(['high_strain_inst_mod_for_each_test', high_strain_inst_mod ])
 out_writer.writerow(['high_strain_eqbm_mod_for_each_test', high_strain_eqbm_mod ])
 
 out_writer.writerow(['mean_of_average_inst_mod_for_test_set', avg_inst_mod_repeat_mean ])
 out_writer.writerow(['SD_of_average_inst_mod_for_test_set', avg_inst_mod_repeat_SD ])
 
 out_writer.writerow(['mean_of_average_eqbm_mod_for_test_set', avg_eqbm_mod_repeat_mean ])
 out_writer.writerow(['SD_of_average_eqbm_mod_for_test_set', avg_eqbm_mod_repeat_SD ])
 
 out_writer.writerow(['mean_of_high_strain_inst_mod_for_test_set', avg_hs_inst_mod_repeat ])
 out_writer.writerow(['SD_of_high_strain_inst_mod_for_test_set', SD_hs_inst_mod_repeat ])
 
 out_writer.writerow(['mean_of_high_strain_eqbm_mod_for_test_set', avg_hs_eqbm_mod_repeat ])
 out_writer.writerow(['SD_of_high_strain_eqbm_mod_for_test_set', SD_hs_eqbm_mod_repeat ])
 
# out_writer.writerow(['mean_of_low_strain_inst_mod_for_test_set', avg_ls_inst_mod_repeat ])
# out_writer.writerow(['SD_of_low_strain_inst_mod_for_test_set', SD_ls_inst_mod_repeat ])
# 
# out_writer.writerow(['mean_of_low_strain_eqbm_mod_for_test_set', avg_ls_eqbm_mod_repeat ])
# out_writer.writerow(['SD_of_low_strain_eqbm_mod_for_test_set', SD_ls_eqbm_mod_repeat ])
 
 
 out_writer.writerow(['coeff_of_var_of_average_inst_mod_for_test_set', cov_avg_inst_mod ])
 out_writer.writerow(['coeff_of_var_of_average_eqbm_mod_for_test_set', cov_avg_eqbm_mod ])
 
 out_writer.writerow(['coeff_of_var_of_high_strain_inst_mod_for_test_set', cov_avg_hs_inst_mod ])
 out_writer.writerow(['coeff_of_var_of_high_strain_eqbm_mod_for_test_set', cov_avg_hs_eqbm_mod ])
 
# out_writer.writerow(['coeff_of_var_of_low_strain_inst_mod_for_test_set', cov_avg_ls_inst_mod ])
# out_writer.writerow(['coeff_of_var_of_low_strain_eqbm_mod_for_test_set', cov_avg_ls_eqbm_mod ])
 
 
 out_writer.writerow(['mean_percent_relaxation_at_first_strain_level_for_test_set',  mean_PR_strain1 ])
 out_writer.writerow(['SD_percent_relaxation_at_first_strain_level_for_test_set',  SD_PR_strain1 ])
 out_writer.writerow(['mean_percent_relaxation_at_second_strain_level_for_test_set',  mean_PR_strain2 ])
 out_writer.writerow(['SD_percent_relaxation_at_second_strain_level_for_test_set',  SD_PR_strain2 ])
 out_writer.writerow(['mean_percent_relaxation_at_third_strain_level_for_test_set',  mean_PR_strain3 ]) 
 out_writer.writerow(['SD_percent_relaxation_at_third_strain_level_for_test_set',  SD_PR_strain3 ]) 

 out_writer.writerow(['coeff_of_var_percent_relaxation_at_first_strain_level_for_test_set',  cov_PR_strain1 ])
 out_writer.writerow(['coeff_of_var_percent_relaxation_at_second_strain_level_for_test_set',  cov_PR_strain2 ])
 out_writer.writerow(['coeff_of_var_percent_relaxation_at_third_strain_level_for_test_set',  cov_PR_strain3 ])
 
 out_writer.writerow(['mean_for_first_strain_level_for_test_set',  mean_strain1 ])
 out_writer.writerow(['SD_for_first_strain_level_for_test_set',  SD_strain1 ])
 out_writer.writerow(['mean_for_second_strain_level_for_test_set',  mean_strain2 ])
 out_writer.writerow(['SD_for_second_strain_level_for_test_set',  SD_strain2 ])
 out_writer.writerow(['mean_for_third_strain_level_for_test_set',  mean_strain3 ])
 out_writer.writerow(['SD_for_third_strain_level_for_test_set',  SD_strain3 ])
 
 out_writer.writerow(['coeff_of_var_for_first_strain_level_for_test_set',  cov_strain1 ])
 out_writer.writerow(['coeff_of_var_for_second_strain_level_for_test_set',  cov_strain2 ])
 out_writer.writerow(['coeff_of_var_for_second_strain_level_for_test_set',  cov_strain3 ])
 

 out_writer.writerow(['mean_dimensions, mm', mean_dimensions ])
 out_writer.writerow(['SD_dimensions, mm', SD_dimensions ])
 out_writer.writerow(['cov_dimensions, mm', cov_dimensions ])

 plt.show()
 
if __name__=="__main__":
 main(sys.argv)