# -*- coding: utf-8 -*-
"""
Created on Wed Aug 17 15:45:26 2022

todo:
    extract the timescales
    plot the maximum swelling level in both directions at different time scales

there is some reference data for CA in the dvs evaluation of agrana CMS coating of the Cerdia samples

@author: csag4207
"""
from os import environ
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from numpy.polynomial import Polynomial #https://numpy.org/doc/stable/reference/routines.polynomials.html

from scipy import optimize
#from matplotlib.ticker import NullFormatter #useful for logit scale
from lmfit import Model, Parameters # conf_interval #, Parameter, report_fit
from PIL import Image

environ['matplotlib']= 'qt' #to automatically plot the figures in a separate window

#%% import the data

# with open # would be better to include to actually close the file afterwards

df = pd.ExcelFile(r'C:\Users\csag4207\Documents\forschung\arbeit\2021\cerdia\data\211119_viscose_ca\211119_viscose_ca_rerun_twice.xlsx')

buffer={}


for i in df.sheet_names:
    if 'Setup' != i and 'equ' not in i:
        buffer[i]=pd.read_excel(df, i)
    if i == 'Setup':
        meta=pd.read_excel(df, i, header=22, na_values=0)

# the 'equ' cannot be imported since it is only pictures 

#%% cut the important part from the meta files

samples = meta.loc[0:23].fillna(0)
samples = samples.drop(samples[samples['name'] == 0].index)
# samples = samples.assign(lowest=0)

measurement = meta.loc[27:, 'sample pos.':'mol mass'].rename(columns=meta.loc[26, 'sample pos.':'mol mass']).fillna(0)
measurement = measurement.drop(measurement[measurement['RH start'] == 0].index)

data={}

for i in buffer:
    if i == 'Ref':
        data[i] = buffer[i].loc[48:].rename(columns=buffer[i].loc[47, :])
    if i != 'Ref':
        try:
            data[samples.loc[int(i[1:])]['name']] = buffer[i].loc[48:].rename(columns=buffer[i].loc[47])
        except KeyError:
            pass

del buffer

# see for: https://stackoverflow.com/questions/47060980/renaming-the-column-names-of-pandas-dataframe-is-not-working-as-expected-pytho
for i in data:
    data[i].rename(columns={'net weight (corrected)       in mg':'net weight'}, inplace=True)
    # data[i].drop(columns=[nan, nan]) # of no importance at the moment

#%% check to see if it worked

spot = 'Ref'

fig, ax1 = plt.subplots(constrained_layout=True)
ax1.plot(data[spot]['time in h'], data[spot]['weight in mg'])

ax1.set_ylabel('reference pan /mg')
ax1.set_xlabel('time /h')


ax2 = ax1.twinx()
ax2.plot(data[spot]['time in h'], data[spot]['RH read in %'], color='k')
ax2.plot(data[spot]['time in h'], data[spot]['RH set in %'], color='b', linestyle=':', alpha=0.7)
ax2.plot(data[spot]['time in h'], data[spot]['temp read in °C'], color='r')

ax2.set_ylabel('RH /% and temp /°C')
plt.show()

#%% so check the excel first: the lowest values should coincide with the min in the data df

for i in data:
    if 'Ref' not in i:
        print(i, data[i]['net weight'].min())
# this seems to coincide - the minima values could be extracted from here

# what else to check ??

#%% plot the data as needed

plt.figure()
for i in data:
    if 'Ref' != i:
        plt.plot(data[i]['time in hours'], data[i]['dm based on lowest net weight in  %'], label=i)
plt.legend()
plt.show()

# in the first 90% RH data, sth must have happened either a window opened or sth else. There was quite some water flux into the samples, therefore the second run of 90% RH was done

#%% plot the relative humidity equilibrium values

linestyles=[':', '-.', '--']
cmap = plt.get_cmap('tab10')
colors = [cmap(i) for i in np.linspace(0,1,4)]

spot='CA' # choose between C and KF

plt.figure()
k=0
for i in data:
    if spot in i:
        if 'citric' not in i:
            plt.plot(data[i]['RH read at equili-brium in %'], data[i]['dm at equilibrium based on lowest net weight in %'], label=i, color= colors[k])
        if 'citric' in i:
            plt.plot(data[i]['RH read at equili-brium in %'], data[i]['dm at equilibrium based on lowest net weight in %'], linestyle=':', label=i, color=colors[k], alpha=0.8)
            k=k+1  
plt.legend()
plt.show()

# the more the fiber itself takes, the lower the effect, possibly also as there is less CMS on the single fiber (sample weight ~ 30 mg) compared to ~100 mg in the CA case
# 

#%%implement the parallel exponential model
# use Kohler 2003 and Thybring 2019 etc. 

def peak(x, a, b, c, d, e, f, g, h):
    return (1-a*np.exp(-(x-h)/b))*c + (1-d*np.exp(-(x-h)/e))*f + g

# def peak(x, a, b, c, d, e, f):
#     return a*np.exp(-x/b)+c + d*np.exp(-x/e)+ f

pek = Model(peak)
parapek = pek.make_params()

eval_pek = np.vectorize(peak)

parapek = {}
parapek['positive'] = Parameters()
parapek['positive'].add_many(('a', 2, True, 0, 75), ('b', 0.1, True, 0.001, 20), ('c', 3, True , 0, 1000), ('d', 2, True, 0, 75), ('e', 4, True, 0, 40), ('f', 0.3, True , 0, 1000), ('g', 25, True , 10, 400)) 
parapek['negative'] = Parameters()
parapek['negative'].add_many(('a', -2, True, -75, 0), ('b', 0.1, True, 0.001, 20), ('c', 3, True , 0, 1000), ('d', -2, True, -75, 0), ('e', 4, True, 0, 40), ('f', 0.3, True , 0, 1000), ('g', 25, True , 10, 400)) 


#%% one first has to divide the data according to the temp / RH set values
# however the number of rows an the index should be similar for all the different data sets

index = {}

for i in data:
    index[i]=[]
    if 'Ref' not in i:
        for j in data[i]['RH set in %'].index[:-1]:
            if data[i]['RH set in %'][j]- data[i]['RH set in %'][j+1] !=0:
                index[i].append(j)
    
result = {}

for i in data:
    if'Ref' not in i:
        result[i]={}
        print('working on ' + i)
        for j in np.arange(len(index[i])):
            parapek['positive'].add('h', data[i]['time in hours'].loc[index[i][j]], False)
            parapek['negative'].add('h', data[i]['time in hours'].loc[index[i][j]], False)
            if j != np.arange(len(index[i]))[-1]:
                if data[i]['net weight'][index[i][j]] <= data[i]['net weight'][index[i][j+1]]:
                    result[i][j]=pek.fit(data[i]['net weight'].loc[index[i][j]:index[i][j+1]], parapek['positive'], x=data[i]['time in hours'].loc[index[i][j]:index[i][j+1]])
                if data[i]['net weight'][index[i][j]] > data[i]['net weight'][index[i][j+1]]:
                    result[i][j]=pek.fit(data[i]['net weight'].loc[index[i][j]:index[i][j+1]], parapek['negative'], x=data[i]['time in hours'].loc[index[i][j]:index[i][j+1]])
            if j == np.arange(len(index[i]))[-1]:
                if data[i]['net weight'][index[i][j]] <= data[i]['net weight'].iloc[-1]:
                    result[i][j]=pek.fit(data[i]['net weight'].loc[index[i][j]:], parapek['positive'], x=data[i]['time in hours'].loc[index[i][j]:])
                if data[i]['net weight'][index[i][j]] > data[i]['net weight'].iloc[-1]:
                    result[i][j]=pek.fit(data[i]['net weight'].loc[index[i][j]:], parapek['negative'], x=data[i]['time in hours'].loc[index[i][j]:])
                    
                # result[i][j]=pek.fit(data[i]['net weight'].loc[index[i][j]:], parapek['positive'], x=data[i]['time in hours'].loc[index[i][j]:])

#%% see the results of the fitting routine

spot = list(data.keys())[2]

fig, axs = plt.subplots()
# axis_numbers = [[0,0], [0,1], [1,0], [1,1]]
# k = -1
axs.scatter(data[spot]['time in hours'], data[spot]['net weight'])
for j in np.arange(len(index[i])):
    if j != np.arange(len(index[i]))[-1]:
        axs.plot(data[spot]['time in hours'].loc[index[spot][j]:index[spot][j+1]], result[spot][j].best_fit, 'b-', label='best fit')
        axs.plot(data[spot]['time in hours'].loc[index[spot][j]:index[spot][j+1]], result[spot][j].init_fit, 'r:', label='init fit')
    if j == np.arange(len(index[i]))[-1]:
        axs.plot(data[spot]['time in hours'].loc[index[spot][j]:], result[spot][j].best_fit, 'b-', label='best fit')
        axs.plot(data[spot]['time in hours'].loc[index[spot][j]:], result[spot][j].init_fit, 'r:', label='init fit')
# axs.legend(loc='best')
axs.set_title(spot)
axs.set_xlabel('time /h')
axs.set_ylabel('mass /mg')
plt.show()

# seems to work so the time scales could be compared
# it seems on the second run there is two opposing e-functions: the fast positive (overshooting) and the slow negative (downregulating) - adapt the parameter accordingly to have both the same sign, though fit quality might decrease

# it is difficult to set the sign of one parameter in the fitting routine according to the other. Is it important? it might be wrong in terms of magnitude, but do the timescales fit? Or would they be shifted in terms of the equivalent sign?

#corrected as it checks first whether up or down and then uses respective restricted parameter

#%% bring the similar thing for the cerdia data

cmap = plt.get_cmap('PuRd') #get a gradual colormap
colors = [cmap(i) for i in np.array([0.2, 0.29375, 1])]

selection_1 = ['viscose KF 28 dtex', '2.5pc NaOH', '10pc NaOH']


plt.figure()
# axs[1]=axs[0].sharey(True)
# axs[2]=axs[0].sharey(True)
k=0
for i in selection_1:
    plt.plot(data[i]['RH read at equili-brium in %'].iloc[2:11], data[i]['dm at equilibrium based on lowest net weight in %'].iloc[2:11], color=colors[k])
    k=k+1
plt.xlabel('RH /%')
plt.ylabel(r'$\Delta$ m /%')
plt.ylim([0, 50])
plt.xlim([0, 95])
plt.legend()
plt.show()

# plt.savefig('C:/Users/csag4207/Documents/forschung/arbeit/2021/cerdia/pic/ca-visc/dvs_fiber_diff-naoh-pre-treatment.png', dpi=500)
# to_save= Image.open('C:/Users/csag4207/Documents/forschung/arbeit/2021/cerdia/pic/ca-visc/dvs_fiber_diff-naoh-pre-treatment.png').convert('RGB')
# to_save.save('C:/Users/csag4207/Documents/forschung/arbeit/2021/cerdia/pic/ca-visc/dvs_fiber_diff-naoh-pre-treatment.jpg')

#%% from arbeit/2022/agrana/data/220802_cms_viscose
# use the 30°C data and fit a polynomial
# calculate from 75 to 89% RH 
# do a respective error

ca_rh = np.array([[4, 0.3552], [26, 2.1355], [50, 4.3486], [74, 8.4706]]) #equilibrium RH dm values

order_2nd = Polynomial.fit(ca_rh[:,0], ca_rh[:,1], deg=2, full=True) #, full=True
plot_2nd = Polynomial(order_2nd)
x = np.arange(2,91)

plt.figure()
plt.scatter(ca_rh[:,0], ca_rh[:,1])
plt.plot(x, order_2nd[0](x), linestyle=':', color='r')
plt.xlim([0, 90])
plt.show()



# the value at 89%RH can be easily calculated
# since the error determination is lengthy, one can use a fixed value for the DVS
value_seeked = order_2nd[0](89)

#%% calculate one error for all samples as this is assume statistically from vibration

selection_2 = ['2.5pc NaOH', '2pc CA viscose', '4pc CA viscose ', '6pc CA viscose', '6pc CA twice', '8pc CA viscose', '8pc CA twice', '10pc CA viscose']

error_cal = []

for i in selection_2:
    error_cal.append([data[i]['dm at equilibrium based on lowest net weight in %'].loc[50], data[i]['dm at equilibrium based on lowest net weight in %'].loc[58]])

error_est = np.mean(error_cal) #use the flattened array

#%% plot the maximum swelling at 90%

selection_3 = ['viscose KF 28 dtex'] + selection_2

labels=['viscose', '0%', '2%', '4%', '6%', '6%', '8%', '8%', '10%', 'CA']

cmap = plt.get_cmap('PuRd')
colors = [cmap(i) for i in np.linspace(0.2, 0.95, selection_3.__len__())]

plt.figure()
k=0
# plt.text(1, 16, s='darker shade: 2.5% citric acid cross-linked', fontsize=25)
for i in selection_3:
    if 'twice' not in i and '6pc' not in i and '8pc' not in i:
        plt.bar(i.split(' ')[0], data[i]['dm at equilibrium based on lowest net weight in %'].loc[54], yerr=error_est, color=colors[k])
    if '6pc' in i or '8pc' in i:
        if 'twice' not in i:
            plt.bar(i.split(' ')[0], data[i]['dm at equilibrium based on lowest net weight in %'].loc[54], yerr=error_est, color=colors[k], align='edge', width=-0.4)
        if 'twice' in i:
            plt.bar(i.split(' ')[0], data[i]['dm at equilibrium based on lowest net weight in %'].loc[54], yerr=error_est, color='k', align='edge', width=0.4, alpha=0.5)
    k=k+1
plt.bar('CA', value_seeked, yerr=1, color=cmap(0.99))
plt.title('Equ. moisture gain at 89% RH')
plt.xticks(rotation=60)
plt.ylabel('mass gain /%')
plt.subplots_adjust(bottom=0.25)
plt.ylim([10, 31])
plt.tight_layout()
plt.xticks([i.split(' ')[0] for i in selection_3]+['CA'], labels)
plt.show()

# plt.savefig('C:/Users/csag4207/Documents/forschung/arbeit/2021/cerdia/pic/ca-visc/dvs_max-water-regain.png', dpi=500)
# to_save= Image.open('C:/Users/csag4207/Documents/forschung/arbeit/2021/cerdia/pic/ca-visc/dvs_max-water-regain.png').convert('RGB')
# to_save.save('C:/Users/csag4207/Documents/forschung/arbeit/2021/cerdia/pic/ca-visc/dvs_max-water-regain.jpg')


#%% plot the values of the initial 66%RH moisture gain

colors = [cmap(i) for i in np.linspace(0.29375, 0.95, selection_2.__len__())]

plt.figure()
k=0
# plt.text(1, 16, s='darker shade: 2.5% citric acid cross-linked', fontsize=25)
for i in selection_2:
    if 'twice' not in i:
        plt.bar(i.split(' ')[0], data[i]['dm at equilibrium based on lowest net weight in %'].loc[53], color=colors[k])
    if 'twice' in i:
        plt.bar(i.split(' ')[0], data[i]['dm at equilibrium based on lowest net weight in %'].loc[53], color='k', alpha=0.5)
    k=k+1
plt.bar('viscose', data['viscose KF 28 dtex']['dm at equilibrium based on lowest net weight in %'].loc[53])
plt.title('Equ. moisture gain at 67% RH')
plt.xticks(rotation=60)
plt.ylabel('mass gain /%')
plt.subplots_adjust(bottom=0.25)
plt.ylim([10, 15])
plt.show()

#%%

plt.figure()
k=0
# plt.text(1, 16, s='darker shade: 2.5% citric acid cross-linked', fontsize=25)
for i in selection_2:
    plt.plot(data[i]['time in hours'].loc[740:1094], data[i]['dm based on lowest net weight in  %'].loc[740:1094], color=colors[k])
        # plt.bar(i, data[i + ' citric acid']['dm at equilibrium based on lowest net weight in %'].loc[55], color='k', alpha=0.4)
    k=k+1
# plt.title('dm based on lowest')

plt.ylabel('mass gain /%')
plt.xlabel('time /h')
plt.ylim([0, 31])
plt.xlim([115, 175])
# plt.tight_layout()
plt.show()

# plt.savefig('C:/Users/csag4207/Documents/forschung/arbeit/2021/cerdia/pic/ca-visc/dvs_dm_based_on_lowest.png', dpi=500)
# to_save= Image.open('C:/Users/csag4207/Documents/forschung/arbeit/2021/cerdia/pic/ca-visc/dvs_dm_based_on_lowest.png').convert('RGB')
# to_save.save('C:/Users/csag4207/Documents/forschung/arbeit/2021/cerdia/pic/ca-visc/dvs_dm_based_on_lowest.jpg')

#%%print the respective timescales 

for i in data:
    if 'Ref' not in i:
        plt.figure()
        for j in result[i]:
            plt.bar(j, result[i][j].best_values['b'], color='r', align='edge', width=-0.4)
            plt.bar(j, result[i][j].best_values['e'], color='k', align='edge', width=0.4, alpha=0.4)
        plt.yscale('log')
        plt.title(i)
        plt.show()

#%%
spotting=5

for i in data:
    if 'Ref' not in i:
        print(i, result[i][spotting].best_values['b'], result[i][spotting].best_values['e'])
    

#%% plot the fitting results
# check if the faster time scale is the more important 

rh_settings_cum=[26, 50, 74, 98, 122, 146]

axises = [[0,0], [0,1], [1,0], [1,1], [2,0], [2,1], [3,0], [3,1]]

fig, axs=plt.subplots(nrows=4, ncols=2)

k = 0
for i in selection_2:
    for j in np.arange(4,10):
        axs[axises[k][0],axises[k][1]].bar(rh_settings_cum[j-4], result[i][j].best_values['c'], color='r')
        axs[axises[k][0],axises[k][1]].bar(rh_settings_cum[j-4], np.multiply(-1, result[i][j].best_values['f']), color='k', alpha=0.4)
    axs[axises[k][0],axises[k][1]].title.set_text(i)
    axs[axises[k][0],axises[k][1]].axhline(0)
    axs[axises[k][0],axises[k][1]].set_ylim([-15, 25])
    k=k+1


plt.suptitle('magnitude of fast(red) and slow(-black) funtion in PEK (30°C)', fontsize=14)
plt.show()

#so the fast and the slow timescales are 


#%% similar for Cerdia

axises = [[0,0], [0,1], [1,0], [1,1]]

fig, axs=plt.subplots(nrows=2, ncols=2)

k = 0
for i in selection_2:
    for j in np.arange(4,10):
        axs[axises[k][0],axises[k][1]].bar(rh_settings_cum[j-4], result[i][j].best_values['c'], color='r')
        axs[axises[k][0],axises[k][1]].bar(rh_settings_cum[j-4], np.multiply(-1, result[i][j].best_values['f']), color='k', alpha=0.4)
    axs[axises[k][0],axises[k][1]].title.set_text(i)
    axs[axises[k][0],axises[k][1]].axhline(0)
    k=k+1
plt.suptitle('magnitude of fast(red) and slow(-black) funtion in PEK (30°C)', fontsize=14)
plt.show()

#so the fast and the slow timescales are 
#%%
rh_settings=[26, 50, 74, 50, 26, 2]

fig, ax1 = plt.subplots()
ax2=ax1.twinx()
k=0
for i in selection_2:
    buff=[]
    for j in np.arange(4, 10):
        buff.append([rh_settings[j-4], result[i][j].best_values['b']])#np.sqrt(np.diag(result[i][j].covar)[1])
    ax1.scatter(np.array(buff)[:3,0], np.array(buff)[:3,1], marker='x', label=i, color=colors[k])
    ax1.plot(np.array(buff)[:,0], np.array(buff)[:,1], label=i, color=colors[k])
    ax2.scatter(np.array(buff)[3:,0], np.array(buff)[3:,1], marker='o', label=i, color=colors[k], alpha=0.5)
    k=k+1
ax1.set_ylabel('sorption timescales /h')
ax2.set_ylabel('desorption timescales /h')
# ax1.set_yscale('log')
# ax1.set_ylim([0.1, 0.5])
# ax2.set_yscale('log')
# ax2.set_ylim([0.1, 0.5])
plt.title('sorption time scales at 30°C (adsorp. x, desorp. o)')
plt.legend()
plt.show()

#%% check to see if a tiny plot brings the differences in the speed constants

selection_4=['2.5pc NaOH', '6pc CA viscose']

plt.figure()
for i in selection_4:
    plt.plot(data[i]['time in hours'], data[i]['dm based on lowest net weight in  %'], label=i)
plt.legend()
plt.show()
