#-*-coding:utf-8-*-
"""
importnumpyasnp
importpandasaspd
frompandasimport read_csv
frommatplotlibimport pyplotasplt
fromscipy.fftpackimportfft,rfft
importmath
importstatistics
fromscipy.signalimportwelch
fromsklearn.decompositionimportPCA
fromitertoolsimport *
fromsklearn.imputeimportSimpleImputer
importplotly.expressas px
defget_fft_values(y_values,T,N, f_s):
f_values=np.linspace(0.0,1.0/(2.0*T),N//2)
fft_values_= rfft(y_values)
fft_values= 2.0/N*np.abs(fft_values_[0:N//2])
returnf_values,fft_values
defget_psd_values(y_values, T,N,f_s):
f_values,psd_values =welch(y_values,fs=f_s)
returnf_values,psd_values
defmoving_avg(x,N):
cumsum=np.cumsum(np.insert(x, 0,0))
return(cumsum[N:]- cumsum[:-N])/N
defpad(seq, target_len,padding= None):
length=len(seq)
iflength> target_len:
raiseTooLongError("sequence toolong ({})for targetlength {}"
.format(length,target_length))
seq.extend([padding]*(target_len -length))
returnseq
cgm_ts = pd.read_csv("C:\\Users\\Shashank\\Documents\\IMP\\ASU\\Courses\\CSE
572\\DataFolder\\CGMDatenumLunchPat1.csv")
cgm = pd.read_csv("C:\\Users\\Shashank\\Documents\\IMP\\ASU\\Courses\\CSE
572\\DataFolder\\CGMSeriesLunchPat1.csv")
#foriin range(0,len(cgm3.index)):
# forjin range(0, len(cgm3.columns)):
# print(cgm3.values[i,j],cgm_ts3.values[i,j],i,j)
#Keep
#means=[]
#maxs=[]
#stdd=[]
#foriin range(33):
# means.append(cgm.loc[i,:].mean())
# maxs.append(cgm.loc[i,:].max())
# stdd.append(cgm.loc[i,:].std())
feature=[]
fl=[]
foriin range(len(cgm.values)):
feature_row =[]
#Feature Type1-RollingWindow
moving_average=moving_avg(cgm.loc[i,:].tolist(), 5).tolist()
moving_average=pad(moving_average, 27,np.nan)
feature_row.extend(moving_average)
plt.figure(i)
plt.plot(moving_avg(cgm_ts.loc[i,:].tolist(),5),moving_avg(cgm.loc[i,:].tolist(),5))
plt.plot(cgm_ts.loc[i,:].tolist(),cgm.loc[i,:].tolist())
#Feature Type2-Analyzingin theFrequencyDomain-FFT &PSD
y_values=cgm.loc[i,:].dropna().unique().tolist()
y_values_ts= cgm_ts.loc[i,:].dropna().unique().tolist()
num=len(y_values)
t_n=y_values[num-1]-y_values[0]
T=t_n/ num
fs=1 /T
f_values,fft_values =get_fft_values(y_values,T, num,fs)
#plt.figure(i+1)
#plt.plot(f_values,fft_values, linestyle= '-',color='blue')
#plt.xlabel('Frequency[Hz]',fontsize=16)
#plt.ylabel('Amplitude',fontsize=16)
fft_values= pad(fft_values.tolist(),15,np.nan)
feature_row.extend(fft_values)
#PSD(PowerSpectralDensity)
psdf_values,psd_values=get_psd_values(y_values, T,num,fs)
#plt.figure(i+2)
#plt.plot(psdf_values, psd_values,linestyle='-', color='blue')
#plt.xlabel('Frequency[Hz]')
#plt.ylabel('PSD[V**2/Hz]')
psd_values =pad(psd_values.tolist(),15,np.nan)
feature_row.extend(psd_values)
##FeatureType 3- StatisticalMethods
feature_row.extend([statistics.stdev(cgm.loc[i,:])])
feature_row.extend([max(cgm.loc[i,:])])
feature_row.append(np.sqrt(np.mean(cgm.loc[i,:]**2)))
#Feature Type4-RateoficreaseofCGMlevelsperunit time(or)Velocity
max_slope =[x- zforx, zinzip(y_values[:-1], y_values[1:])]
max_slope_time=[x-zforx,zinzip(y_values_ts[:-1],y_values_ts[1:])]
velocity= [y/w fory,w inzip(max_slope, max_slope_time)]
velocity_zero_crossings= np.where(np.diff(np.signbit(velocity)))[0]
velocity=pad(velocity,28,np.nan)
velocity_zero_crossings =pad(velocity_zero_crossings.tolist(),20,np.nan)
velocity_mean=np.mean(velocity)
velocity_max=np.max(velocity)
#plt.figure(i+3)
#plt.plot(y_values_ts,pad(velocity,int(len(y_values_ts)),0))
#plt.xlabel('time')
#plt.ylabel('velocity')
feature_row.extend(velocity)
feature_row.extend(velocity_zero_crossings)
feature_row.append(velocity_mean)
feature_row.append(velocity_max)
feature.append(feature_row)
fl.append(len(feature_row))
print(feature)
feature_2D_arr=np.asarray(feature).reshape(33,110)
print(feature_2D_arr,feature_2D_arr.shape)
##Replacingmissingvalueswith imputer
imp_mean= SimpleImputer(missing_values=np.nan,strategy='mean')
feature_ip= imp_mean.fit_transform(feature_2D_arr)
pca=PCA(n_components=5)
principalComps=pca.fit_transform(feature_ip)
print(pca.explained_variance_ratio_)
components=abs(pca.components_)
variances=pca.explained_variance_
x=[ifor iinrange(0,len(components[0]))]
foriin range(0,5):
plt.figure(figsize=(200,100))
plt.bar(x,components[i])
plt.xticks(np.arange(len(components[0])),x,rotation=90)
plt.show()
positives =np.array(np.argwhere(components[i] >0).flatten())
positive_sorted=np.argsort(components[i][:])
positives =positive_sorted
print(components[i])
# print(positives)
# print(columns)
print(variances[i])
#print(columns[positives])
print(components[i][positives])
print(components[i])
#Radarplotwithallthefeatures
radarPlotPCA=pd.DataFrame(dict(
r=pca.explained_variance_ratio_,
theta=['PC1','PC2','PC3',
'PC4','PC5']))
fig=px.line_polar(radarPlotPCA,r='r',theta='theta',line_close=True)
fig.show()
fig.write_image("fig1.png")
