# import a few packages
%matplotlib notebook
from thermocouples_reference import thermocouples
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import sympy as sp
from scipy import optimize, interpolate, signal


typeC=thermocouples['C']

# make sure you are in the same dir as the file
# read in the file and drop Na cols
df = pd.read_excel('Type C Table 4-2-18.xlsx')
df.dropna(axis=1, inplace=True)
df.head()

# NIST has values calibrated for T > 273.15 K, lets find the Tref based on these points
# I am using Kelvin for all T. The CJC is quoted in deg C.
tempdf = df.query('T>273.15')
tempdf.head()

# Let's find the T_ref by using this function to take the TypeC mV and the T to find the Tref
def find_Tref(mV, T):
    x = np.arange(290, 301, 0.01)
    x = x[::-1] # lets reverse x
    i = 1
    while typeC.inverse_KmV(mV, Tref=x[i]) - T >= 0:
        i += 1
#     print(x[i])
    return x[i]

# This isn't the fastest way to do things, but since its just a short amount of rows, lets iterate over the mV and T
# to find Tref
Treflist=[]
for idx in tempdf.index:
    # print(idx)
    Treflist.append(find_Tref(mV=tempdf['TypeCmV'][idx], T=tempdf['T'][idx]))


print( ['%0.2f'% x for x in Treflist])

# now average the Trefs:
avg_Tref = np.mean(Treflist)
print(avg_Tref)
# I will use this Tref for further calcs
Tref_emf = typeC.emf_mVK(avg_Tref)
print(Tref_emf)
# The Tref_emf value is very close to the value in the table at 273.15 K, so we'll use this value to correct the new values
# The value taken at 273.15 K was during the cooling process and is likely to be less accurate than the room temperature value
# across these multiple observations

# The emf correction for 273.15 K is then: calibrated_emf = raw_emf + Tref_emf

# Let's add this to the df we initially imported...

df['TypeC_calib_mV'] = df['TypeCmV'] + Tref_emf

df.head()
# Compared to the NIST table, we appear to be off at most a little less than 1 deg K


# Had we used the CJC temperature as a proxy for room temp, we would've been even more off.
# compare the TypeCmV using Tref = CJC vs using Tref = 294.67:
print(typeC.emf_mVK(291.22, Tref =(25.26+273.15)))
print(typeC.emf_mVK(291.22, Tref =avg_Tref))

#  Let's visualize these results
plt.plot(df['T'], df['TypeC_calib_mV'], 'o', ms=0.5 )
plt.xlabel('Temperature (K)')
plt.ylabel('Type C calibrated emf (mV)')

# Interesting. I cooled first to LN2 temperatures and then allowed the sample to heat up slowly by evaporating LN2
# The data agrees fairly well (within ~3 K) between the heating and cooling curves. I didn't heat all the way back up.

# Now lets fit the data to a polynowmial using least squares
fit_coeffs = np.polyfit(df['T'],df['TypeC_calib_mV'], deg = 10 , full=True)
# print(fit_coeffs)
fit_poly = np.poly1d(fit_coeffs[0])
print(fit_poly)


fig, ax = plt.subplots()
ax.plot(df['T'], df['TypeC_calib_mV'],'o',ms='0.5')
ax.plot(df['T'], fit_poly(df['T']) , 'o', ms='0.5')


"""
Explanation: The situation
Type C thermocouples are not NIST calibrated to below 273.15 K. For my research specific scenario, I need to cool my sample (Molybdenum) to cryogenic temperatures and also anneal to very high ~2000 K. There is no thermocouple with these properties. 
The solution
We know that Type K thermocouples are accurate down to cryogenic temperatures. So what I've done here is to read the Type K temperature and record the corresponding Type C mV to create a calibration table. Both thermocouples were spot welded to a large mass very close to one another to ensure the temperature readings will be accurate.
Then I will use a polynomial fit to get the low T calibration for the Type C thermocouple.
End of explanation
"""


# These mV values are also close ~0.5 degrees K of one another
print(fit_poly(273.15))  # fit
print(typeC.emf_mVK(273.15))  # NIST value  

"""
Explanation: The 10th degree polynomial appears to give the best fit overall. 
The lower order polynomials dont fit the curve exceedingly well below 100 K
Also, the polynomial tracks the heating curve (the slightly higher mV points from 80-150K) a little more closely than the cooling curve (295 to 80 K). Heating occurred much more slowly than cooling, so I expect it to me more accurate anyways.
End of explanation
"""


x = sp.symbols('x')
polynom = sp.Poly(fit_coeffs[0],x)
# print(fit_coeffs[0])
# find the second derivative of the polynomial
second_derivative = polynom.diff(x,x)
print(second_derivative)
sp.solve(second_derivative,x, domain= sp.S.Reals)

print(second_derivative.evalf(subs={x:77}))
print(second_derivative.evalf(subs={x:80}))
print('\n')
print(second_derivative.evalf(subs={x:120}))
print(second_derivative.evalf(subs={x:125}))
print('\n')
print(second_derivative.evalf(subs={x:135}))
print(second_derivative.evalf(subs={x:145}))
print('\n')
print(second_derivative.evalf(subs={x:283}))
print(second_derivative.evalf(subs={x:291}))

first_deriv = polynom.diff(x)
print(first_deriv)
sp.solve(first_deriv,x, domain= sp.S.Reals)

print(first_deriv.evalf(subs={x:80}))
print(first_deriv.evalf(subs={x:84}))

"""
Explanation: It's also a good idea to check that the polynomial does not have any inflection points, at least in the area we are interested in using the polynomial (77 K - 273.15 K). We can use the second derivative test to see if this will be important for our case.
End of explanation
"""


lowT_df = df.query('T<103')
# Now lets fit the data to an exponential
# print(np.min(lowT_df['TypeC_calib_mV']))
def func(x, a, b, c, d):
    return a * np.exp(b * x - c) + d

fit_coeffs = optimize.curve_fit(func, lowT_df['T'],lowT_df['TypeC_calib_mV'], p0=(1, 1, 90, -3))
print(fit_coeffs)
a = fit_coeffs[0][0]
b = fit_coeffs[0][1]
c = fit_coeffs[0][2]
d =  fit_coeffs[0][3]
expfunc = func(lowT_df['T'],a,b,c,d)

fig3, ax3 = plt.subplots()
# ax3.plot(lowT_df['T'], a*np.exp(b*lowT_df['TypeC_calib_mV']), 'o',ms='0.5')
ax3.plot(lowT_df['T'], lowT_df['TypeC_calib_mV'], 'o',ms='0.5')
ax3.plot(lowT_df['T'], expfunc, 'o',ms='0.5',color='r')

"""
Explanation: Well this is not optimal-- there exists a local minimum at 83.86 K in our polynomial fit. We can attempt to fit an exponential curve to this very low temperature data and append this to the polynomial function.
End of explanation
"""


# select data from 103 to 120 K just so we can see the point of intersection a little better
checkT_df = df.query('77<=T<=120') 
fig4, ax4 = plt.subplots()
ax4.plot(checkT_df['T'], fit_poly(checkT_df['T']), 'o', ms=0.5, label='polyfit', color='g')
ax4.plot(lowT_df['T'], expfunc, 'o', ms=0.5, label='expfunc', color='r')
ax4.plot(df['T'], df['TypeC_calib_mV'],'o',ms='0.5', label='Data', color='b')
ax4.set_xlim([80,110])
ax4.set_ylim([-1.88,-1.75])
ax4.legend()

"""
Explanation: This appears to be a better fit than the polynomial in this regime. Now lets concatenate these two functions and interpolate near the points around 100 K to smooth things out if necessary. Recall that the two functions are fit_poly and expfunc
End of explanation
"""


def switch_fcn(x, switchpoint, smooth):
    s = 0.5 + 0.5*np.tanh((x - switchpoint)/smooth)
    return s
sw = switch_fcn(df['T'], 103, 0.2)
expfunc2 = func(df['T'],a,b,c,d)
len(expfunc2)

fig, ax = plt.subplots()
ax.plot(df['T'], sw,'o', ms=0.5)


def combined(switch, low_f1, high_f2):
    comb = (1-switch)*low_f1 + switch*high_f2
    return comb
comb_fcn = combined(sw, expfunc2,fit_poly(df['T']))
fig, ax = plt.subplots()
ax.plot(df['T'], comb_fcn, 'o', ms=0.5)

fig5, ax5 = plt.subplots()
ax5.plot(df['T'],comb_fcn, 'o', ms=2, label='combined')
ax5.plot(checkT_df['T'], fit_poly(checkT_df['T']), 'o', ms=0.5, label='polyfit', color='g')
ax5.plot(lowT_df['T'], expfunc, 'o', ms=0.5, label='expfunc2', color='r')
ax5.set_xlim([80,110])
ax5.set_ylim([-1.88,-1.75])
ax5.legend()

"""
Explanation: The two fitted plots almost match near 103 K, but there is a little 'cusp'-like shape near the point of intersection. Let's smooth it out. Also, notice that the expfunc fit is a little better than the polyfit.
End of explanation
"""


# low temperature array
low_temp = np.arange(77.15,273.15, 0.1)
# low_temp_calib = fit_poly(low_temp)
low_temp_calib = combined(switch_fcn(low_temp, 103, 3), func(low_temp,a,b,c,d), fit_poly(low_temp))

# high temperature array
high_temp = np.arange(273.15,2588.15, 0.1)
high_temp_nist = typeC.emf_mVK(high_temp)

# concatentate and put into a dataframe and output to excel
Temperature = np.concatenate([low_temp, high_temp])
TypeC_mV = np.concatenate([low_temp_calib, high_temp_nist])

typeC_calibration = pd.DataFrame(data=TypeC_mV, index=Temperature, dtype='float32', columns = ['Type C (mV)'])
typeC_calibration.index.name = 'Temperature (Kelvin)'

print(typeC_calibration.head())
print(typeC_calibration.tail())

# Uncomment these lines and run the cell to output a calibration table
# write to excel

# xlwrite = pd.ExcelWriter('Type C calibration_low_res.xlsx')
# typeC_calibration.to_excel(xlwrite)
# xlwrite.save()

"""
Explanation: Now I will take the polynomial and take the values from 77 K to 273 K for calibration and append them to the NIST values
End of explanation
"""



low_calib = combined(switch_fcn(Temperature, 103, 3), func(Temperature,a,b,c,d), fit_poly(Temperature))

high_calib = pd.DataFrame(index=high_temp, data=high_temp_nist,columns=['mV'])
dummy_df = pd.DataFrame(index=low_temp, data=np.zeros(len(low_temp)),columns=['mV'])
concat_high_calib = dummy_df.append(high_calib)
print(concat_high_calib.loc[272.9:273.5])

freezept_calib = combined(switch_fcn(Temperature, 273.15, 0.45), low_calib, concat_high_calib['mV'] )
freezept_calib.index.name = 'T'

freezept_calib.loc[272.9:273.5]

"""
Explanation: But wait! Suppose we also want to fix that discontinuity at 273.15 K? We can apply the same procudure as before.
1. Apply a tanh(x) function:  $switch = 0.5 + 0.5np.tanh((x - switchpoint)/smooth)$
2. Combine both functions:  $comb = (1-switch)f1 + (switch)*f2 $
End of explanation
"""


fig, ax = plt.subplots()
freezept_calib.plot(ax=ax, label ='combined')
ax.plot(Temperature,low_calib, label = 'low calib')
ax.plot(Temperature,concat_high_calib, label= 'high_calib')
ax.set_ylim([-.04,0.04])
ax.set_xlim([268,277])
ax.legend()

print(signal.argrelmin(freezept_calib.values))
# print(signal.argrelextrema(freezept_calib.values,np.less))
# print(signal.argrelextrema(freezept_calib.values,np.greater))
# No local maxima or minima!

# How about candidates for inflection points?

df = np.gradient(freezept_calib,0.1,)
fig, ax = plt.subplots()
ax.plot(Temperature, df)

d2f = np.gradient(df, 0.1)
ax.plot(Temperature, d2f)

d2f[:5]

# Uncomment these lines and run the cell to output a calibration table
# write to excel

xlwrite = pd.ExcelWriter('Type C calibration_corrected_temp.xlsx')
# freezept_calib is a Series, not a Dataframe, so use the line below
freezept_calib.to_frame().to_excel(xlwrite)
xlwrite.save()

# 4212018
fig, ax = plt.subplots()
ax.plot(Temperature,TypeC_mV)

# np.fft.fft(TypeC_mV)

def scaled_data(data, low_range=0, high_range=1, standardize=False, print_=False):
    """
    scale data from input range to (low_range,high_range)
    assumes data is a np 1d array
    also allows capability for standardization (mean=0 and variance = 1)
    """
    _min_ = np.min(data)
    _max_ = np.max(data)
    
    if standardize is True:
        data = data - np.mean(data) #remove mean
        scaled_data = data/np.std(data)    #unit std
        if print_ is True:
            print('mean: '+str(np.mean(scaled_data)))
            print('std: '+str(np.std(scaled_data)))
    else:
        scaled_data = (high_range - low_range)*(data - _min_)/(_max_ - _min_) + low_range
    if print_ is True:
        print('max: = '+str(np.max(scaled_data)))
        print('min: = '+str(np.min(scaled_data)))
    
    return scaled_data

def revert_to_unscaled(scaled, original):
    """
    reverts normzlied data back to original scaling
    """
    scaled_min = np.min(scaled)
    scaled_max = np.max(scaled)
    orig_min = np.min(original)
    orig_max = np.max(original)
    
    data = (scaled - scaled_min)*(orig_max - orig_min)/(scaled_max - scaled_min) + orig_min
    
    return data

normmV = scaled_data(TypeC_mV, standardize=False, print_=True)
normT = scaled_data(Temperature, standardize=False, print_=True)
stdmV = scaled_data(TypeC_mV, standardize=True, print_=True)
stdT = scaled_data(Temperature, standardize=True, print_=True)

fig,ax = plt.subplots()
ax.plot(stdmV, stdT, label ='std')
plt.xlabel('stdmV')
plt.ylabel('stdT')

fig,ax = plt.subplots()
ax.plot(normmV, normT, label ='norm')
plt.xlabel('normmV')
plt.ylabel('normT')

def fit_fcn(mV, param):
    return param[0]*np.tanh(param[1]*mV+param[2]) +param[3]
fig, ax = plt.subplots()

ax.plot(normmV,normT, label='norm')
ax.plot(normmV,fit_fcn(normmV, param=[1,0.1,-1, 1]), label='start')
try:
#     popt, pcov = optimize.differential_evolution(fit_fcn, normmV, normT, p0=(2,1,-1,2), method='lm')
    bounds = [slice(0,100.,1), slice(-10.,100.,1), slice(-10.,100.,1),slice(-10.,100.,1)]
#     result = optimize.differential_evolution(lambda param:np.sum(fit_fcn(normmV, param) - normT)**2,bounds)
#     result = optimize.basinhopping(lambda param:np.sum(fit_fcn(normmV, param) - normT)**2,bounds)
    result = optimize.brute(lambda param:np.sum(fit_fcn(normmV, param) - normT)**2,bounds)
#     result = optimize.least_squares(lambda param:np.sum(fit_fcn(normmV, param) - normT)**2,[1,1,-1, 1])

    print(result)
    ax.plot(normmV,fit_fcn(normmV, result), label='fitted')
except RuntimeError:
    print("curve fit failed")

plt.legend()

"""
Explanation: The prior value at 273.15 K was -0.00867, when the actual value is 0. After the smoothing, the new value is -0.004336, about half of the prior value. Some of the values a little after 273.15 do not match exactly with the NIST table, but it is much better than the jump that we had before.
End of explanation
"""


