Boulder_4hour.py

# coding: utf-8

# ### Exp 1.2 - RNN - Boulder- 2 HOUR summary:
# 1. import bird model ghi values
# 2. import actual data
# 3. merge the two based on the common timestamps
# 4. drop all rows if 'ghi' value is 0
# 5. ** DO NOT drop all rows with dw_solar < 0 **
# 6. Manage missing values: Drop any row if it has -9999.9
# 7. Make Kt column ('dw_solar'/'ghi') 
# 8. Group the data into hourly format (by finding mean)
# 9. assign hierarchical index
# 10. get the subdataframes for each day and shift the Kt values
# 11. drop the year, month, day, hour and minute columns
# 12. feed the data into RNN alrogithm
#     - RNN seq_lenth = 1
#     - learning rate = 0.001
#     - single RNN layer
#     - num of epochs = 5000
# 
# ** in this experiment we are not giving Kt as one of the features. We are using hierarchical index to make sub dataframe for every day and then shift the Kt values one cell up to make it as one hour forecast. **

# ### NOTE
# ** the test['Kt'].max() values after grouping into hourly data is 7279.907547 = this is a outlier. The reason might be having a dw_solar reading with QC=2 value. Still haven't checked that. But the second highest value of 'Kt' parameter in the testset is 672 (which is way lower) **

# In[1]:

import numpy as np
import pandas as pd
import datetime
import glob
import os.path
from pandas.compat import StringIO


# ### NREL Bird Model implementation: for obtaining clear sky GHI

# In[2]:

import itertools
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns


# In[3]:

#get_ipython().magic('matplotlib inline')
sns.set_color_codes()


# In[4]:

import pvlib
from pvlib import clearsky, atmosphere
from pvlib.location import Location


# In[5]:

bod = Location(40.0150,-105.2705, 'US/Mountain', 1655.064, 'Boulder')


# In[7]:

times2009 = pd.DatetimeIndex(start='2009-01-01', end='2010-01-01', freq='1min',
                        tz=bod.tz)   # 12 months of 2009 - For testing
times2010and2011 = pd.DatetimeIndex(start='2010-01-01', end='2012-01-01', freq='1min',
                        tz=bod.tz)   # 24 months of 2010 and 2011 - For training
#times2011 = pd.DatetimeIndex(start='2011-01-01', end='2012-01-01', freq='1min',
                       # tz=bod.tz)   # 12 months of 2011  - For training


# In[8]:

cs_2009 = bod.get_clearsky(times2009) 
cs_2010and2011 = bod.get_clearsky(times2010and2011) # ineichen with climatology table by default
#cs_2011 = bod.get_clearsky(times2011) 


# In[9]:

cs_2009.drop(['dni','dhi'],axis=1, inplace=True) #updating the same dataframe by dropping two columns
cs_2010and2011.drop(['dni','dhi'],axis=1, inplace=True) #updating the same dataframe by dropping two columns
#cs_2011.drop(['dni','dhi'],axis=1, inplace=True) #updating the same dataframe by dropping two columns


# In[10]:

cs_2009.reset_index(inplace=True)
cs_2010and2011.reset_index(inplace=True)
#cs_2011.reset_index(inplace=True)


# In[11]:

cs_2009['index']=cs_2009['index'].apply(lambda x:x.to_datetime())
cs_2009['year'] = cs_2009['index'].apply(lambda x:x.year)
cs_2009['month'] = cs_2009['index'].apply(lambda x:x.month)
cs_2009['day'] = cs_2009['index'].apply(lambda x:x.day)
cs_2009['hour'] = cs_2009['index'].apply(lambda x:x.hour)
cs_2009['min'] = cs_2009['index'].apply(lambda x:x.minute)


# In[12]:

cs_2010and2011['index']=cs_2010and2011['index'].apply(lambda x:x.to_datetime())
cs_2010and2011['year'] = cs_2010and2011['index'].apply(lambda x:x.year)
cs_2010and2011['month'] = cs_2010and2011['index'].apply(lambda x:x.month)
cs_2010and2011['day'] = cs_2010and2011['index'].apply(lambda x:x.day)
cs_2010and2011['hour'] = cs_2010and2011['index'].apply(lambda x:x.hour)
cs_2010and2011['min'] = cs_2010and2011['index'].apply(lambda x:x.minute)


# In[13]:

print(cs_2009.shape)
print(cs_2010and2011.shape)
#print(cs_2011.shape)


# In[15]:

cs_2009.drop(cs_2009.index[-1], inplace=True)
cs_2010and2011.drop(cs_2010and2011.index[-1], inplace=True)
#cs_2011.drop(cs_2011.index[-1], inplace=True)


# In[16]:

print(cs_2009.shape)
print(cs_2010and2011.shape)
#print(cs_2011.shape)


# In[18]:

cs_2010and2011.head()


# ### Import files from each year in a separate dataframe

# 
# - year            integer	 year, i.e., 1995
# - jday            integer	 Julian day (1 through 365 [or 366])
# - month           integer	 number of the month (1-12)
# - day             integer	 day of the month(1-31)
# - hour            integer	 hour of the day (0-23)
# - min             integer	 minute of the hour (0-59)
# - dt              real	 decimal time (hour.decimalminutes, e.g., 23.5 = 2330)
# - zen             real	 solar zenith angle (degrees)
# - dw_solar        real	 downwelling global solar (Watts m^-2)
# - uw_solar        real	 upwelling global solar (Watts m^-2)
# - direct_n        real	 direct-normal solar (Watts m^-2)
# - diffuse         real	 downwelling diffuse solar (Watts m^-2)
# - dw_ir           real	 downwelling thermal infrared (Watts m^-2)
# - dw_casetemp     real	 downwelling IR case temp. (K)
# - dw_dometemp     real	 downwelling IR dome temp. (K)
# - uw_ir           real	 upwelling thermal infrared (Watts m^-2)
# - uw_casetemp     real	 upwelling IR case temp. (K)
# - uw_dometemp     real	 upwelling IR dome temp. (K)
# - uvb             real	 global UVB (milliWatts m^-2)
# - par             real	 photosynthetically active radiation (Watts m^-2)
# - netsolar        real	 net solar (dw_solar - uw_solar) (Watts m^-2)
# - netir           real	 net infrared (dw_ir - uw_ir) (Watts m^-2)
# - totalnet        real	 net radiation (netsolar+netir) (Watts m^-2)
# - temp            real	 10-meter air temperature (?C)
# - rh              real	 relative humidity (%)
# - windspd         real	 wind speed (ms^-1)
# - winddir         real	 wind direction (degrees, clockwise from north)
# - pressure        real	 station pressure (mb)
# 

# In[19]:

cols = ['year', 'jday', 'month', 'day','hour','min','dt','zen','dw_solar','dw_solar_QC','uw_solar',
       'uw_solar_QC', 'direct_n','direct_n_QC','diffuse', 'diffuse_QC', 'dw_ir', 'dw_ir_QC', 'dw_casetemp',
       'dw_casetemp_QC', 'dw_dometemp','dw_dometemp_QC','uw_ir', 'uw_ir_QC', 'uw_casetemp','uw_casetemp_QC',
       'uw_dometemp','uw_dometemp_QC','uvb','uvb_QC','par','par_QC','netsolar','netsolar_QC','netir','netir_QC',
       'totalnet','totalnet_QC','temp','temp_QC','rh','rh_QC','windspd','windspd_QC','winddir','winddir_QC',
       'pressure','pressure_QC']


# In[20]:

path = r'./data/Boulder/Exp_1_2_train'
all_files = glob.glob(path + "/*.dat")
all_files.sort()

df_big_train = pd.concat([pd.read_csv(f, skipinitialspace = True, quotechar = '"',skiprows=(2),delimiter=' ', 
                 index_col=False,header=None, names=cols) for f in all_files],ignore_index=True)
df_big_train.shape


# In[22]:

path = r'./data/Boulder/Exp_1_2_test'
all_files = glob.glob(path + "/*.dat")
all_files.sort()

df_big_test = pd.concat((pd.read_csv(f, skipinitialspace = True, quotechar = '"',skiprows=(2),delimiter=' ', 
                 index_col=False,header=None, names=cols) for f in all_files),ignore_index=True)
df_big_test.shape


# In[25]:

df_big_test[df_big_test['dw_solar']==-9999.9].shape


# ### Merging Clear Sky GHI And the big dataframe

# In[26]:

df_train = pd.merge(df_big_train, cs_2010and2011, on=['year','month','day','hour','min'])
df_train.shape


# In[27]:

df_test = pd.merge(df_big_test, cs_2009, on=['year','month','day','hour','min'])
df_test.shape


# In[28]:

df_train.drop(['index'],axis=1, inplace=True) #updating the same dataframe by dropping the index columns from clear sky model
df_test.drop(['index'], axis=1, inplace=True)


# In[29]:

df_train.shape


# ### Managing missing values

# In[30]:

# Resetting index
df_train.reset_index(drop=True, inplace=True)
df_test.reset_index(drop=True, inplace=True)


# In[104]:

# Dropping rows with two or more -9999.9 values in columns


# In[31]:

# Step1: Get indices of all rows with 2 or more -999
missing_data_indices = np.where((df_train <=-9999.9).apply(sum, axis=1)>=2)[0]
# Step2: Drop those indices
df_train.drop(missing_data_indices, axis=0, inplace=True)
# Checking that the rows are dropped
df_train.shape


# In[32]:

missing_data_indices_test = np.where((df_test <= -9999.9).apply(sum, axis=1)>=2)[0]
df_test.drop(missing_data_indices_test, axis=0, inplace=True)
df_test.shape


# In[111]:

# For the rows with only one cell as -9999.9, replacing this cell with the mean of the column


# #### First resetting index after dropping rows in the previous part of the code

# In[33]:

# 2nd time - Reseting Index
df_train.reset_index(drop=True, inplace=True)
df_test.reset_index(drop=True, inplace=True)


# In[34]:

one_miss_train_idx = np.where((df_train <=-9999.9).apply(sum, axis=1)==1)[0]


# In[35]:

len(one_miss_train_idx)


# In[36]:

df_train.shape


# In[37]:

col_names = df_train.columns
from collections import defaultdict
stats = defaultdict(int)
total_single_missing_values = 0
for name in col_names:
    col_mean = df_train[~(df_train[name] == -9999.9)][name].mean()
    missing_indices = np.where((df_train[name] == -9999.9))
    stats[name] = len(missing_indices[0])
    df_train[name].loc[missing_indices] = col_mean
    total_single_missing_values += sum(df_train[name] == -9999.9)
    


# In[38]:

df_col_min = df_train.apply(min, axis=0)
df_col_max = df_train.apply(max, axis =0)
#print(df_col_min, df_col_max)


# In[39]:

train = np.where((df_train <=-9999.9).apply(sum, axis=1)==1)[0]


# In[40]:

len(train)


# In[115]:

# doing the same thing on test dataset


# In[41]:

one_miss_test_idx = np.where((df_test <=-9999.9).apply(sum, axis=1)==1)[0]
len(one_miss_test_idx)


# In[42]:

col_names_test = df_test.columns
from collections import defaultdict
stats_test = defaultdict(int)
total_single_missing_values_test = 0
for name in col_names_test:
    col_mean = df_test[~(df_test[name] == -9999.9)][name].mean()
    missing_indices = np.where((df_test[name] == -9999.9))
    stats_test[name] = len(missing_indices[0])
    df_test[name].loc[missing_indices] = col_mean
    total_single_missing_values_test += sum(df_test[name] == -9999.9)
    


# In[43]:

test = np.where((df_test <=-9999.9).apply(sum, axis=1)==1)[0]


# In[44]:

len(test)


# In[45]:

df_train.shape


# In[46]:

df_test.shape


# ### Exploratory Data Analysis

# In[47]:

dw_solar_everyday = df_test.groupby(['jday'])['dw_solar'].mean()
ghi_everyday = df_test.groupby(['jday'])['ghi'].mean()
j_day = df_test['jday'].unique()


# In[50]:

fig = plt.figure()

axes1 = fig.add_axes([0.1,0.1,0.8,0.8])
#axes2 = fig.add_axes([0.1,0.1,0.8,0.8])

axes1.scatter(j_day,dw_solar_everyday,label='Observed dw_solar',color='red')
axes1.scatter(j_day, ghi_everyday, label='Clear Sky GHI',color='green')

axes1.set_xlabel('Days')
axes1.set_ylabel('Solar Irradiance (Watts /m^2)')
axes1.set_title('Solar Irradiance - Test Year 2009')
axes1.legend(loc='best')

fig.savefig('./RNN Paper Results/Exp1_2/Boulder/4hour_Figure 2.jpg', bbox_inches = 'tight')


# In[51]:

sns.jointplot(x=dw_solar_everyday,y=ghi_everyday,kind='reg')
#plt.title('observed dw_solar vs clear sky ghi')
plt.xlabel('Observed global downwelling solar (Watts/m^2)')
plt.ylabel('Clear Sky GHI (Watts/m^2)')
plt.savefig('./RNN Paper Results/Exp1_2/Boulder/4hour_Figure 3', bbox_inches='tight')


# ### making the Kt (clear sky index at time t) column by first removing rows with ghi==0

# In[52]:

df_train = df_train[df_train['ghi']!=0]
df_test = df_test[df_test['ghi']!=0]
df_train['Kt'] = df_train['dw_solar']/df_train['ghi']
df_test['Kt'] = df_test['dw_solar']/df_test['ghi']


# In[56]:

df_train.reset_index(inplace=True)
df_test.reset_index(inplace=True)


# In[70]:

print("test Kt max: "+str(df_test['Kt'].max()))
print("test Kt min: "+str(df_test['Kt'].min()))
print("test Kt mean: "+str(df_test['Kt'].mean()))
print("\n")
print("train Kt max: "+str(df_train['Kt'].max()))
print("train Kt min: "+str(df_train['Kt'].min()))
print("train Kt mean: "+str(df_train['Kt'].mean()))


# In[73]:

plt.plot(df_train['Kt'])


# In[75]:

plt.plot(df_test['Kt'])


# In[72]:

df_train= df_train[df_train['Kt']< 5000]
df_train= df_train[df_train['Kt']> -1000]
df_test= df_test[df_test['Kt']< 5000]
df_test= df_test[df_test['Kt']> -1000]


# #### Group the data (train dataframe)

# In[101]:

zen = df_train.groupby(['year','month','day','hour'])['zen'].mean()
dw_solar = df_train.groupby(['year','month','day','hour'])['dw_solar'].mean()
uw_solar = df_train.groupby(['year','month','day','hour'])['uw_solar'].mean()
direct_n = df_train.groupby(['year','month','day','hour'])['direct_n'].mean()
diffuse = df_train.groupby(['year','month','day','hour'])['diffuse'].mean()
dw_ir = df_train.groupby(['year','month','day','hour'])['dw_ir'].mean()
dw_casetemp = df_train.groupby(['year','month','day','hour'])['dw_casetemp'].mean()
dw_dometemp = df_train.groupby(['year','month','day','hour'])['dw_dometemp'].mean()
uw_ir = df_train.groupby(['year','month','day','hour'])['uw_ir'].mean()
uw_casetemp = df_train.groupby(['year','month','day','hour'])['uw_casetemp'].mean()
uw_dometemp = df_train.groupby(['year','month','day','hour'])['uw_dometemp'].mean()
uvb = df_train.groupby(['year','month','day','hour'])['uvb'].mean()
par = df_train.groupby(['year','month','day','hour'])['par'].mean()
netsolar = df_train.groupby(['year','month','day','hour'])['netsolar'].mean()
netir = df_train.groupby(['year','month','day','hour'])['netir'].mean()
totalnet = df_train.groupby(['year','month','day','hour'])['totalnet'].mean()
temp = df_train.groupby(['year','month','day','hour'])['temp'].mean()
rh = df_train.groupby(['year','month','day','hour'])['rh'].mean()
windspd = df_train.groupby(['year','month','day','hour'])['windspd'].mean()
winddir = df_train.groupby(['year','month','day','hour'])['winddir'].mean()
pressure = df_train.groupby(['year','month','day','hour'])['pressure'].mean()
ghi = df_train.groupby(['year','month','day','hour'])['ghi'].mean()
Kt = df_train.groupby(['year','month','day','hour'])['Kt'].mean()


# In[102]:

df_new_train = pd.concat([zen,dw_solar,uw_solar,direct_n,diffuse,dw_ir,dw_casetemp,dw_dometemp,uw_ir,uw_casetemp,uw_dometemp,
                    uvb,par,netsolar,netir,totalnet,temp,rh,windspd,winddir,pressure,ghi,Kt], axis=1)


# In[103]:

df_new_train.head()


# #### Groupdata - test dataframe

# In[78]:

test_zen = df_test.groupby(['month','day','hour'])['zen'].mean()
test_dw_solar = df_test.groupby(['month','day','hour'])['dw_solar'].mean()
test_uw_solar = df_test.groupby(['month','day','hour'])['uw_solar'].mean()
test_direct_n = df_test.groupby(['month','day','hour'])['direct_n'].mean()
test_diffuse = df_test.groupby(['month','day','hour'])['diffuse'].mean()
test_dw_ir = df_test.groupby(['month','day','hour'])['dw_ir'].mean()
test_dw_casetemp = df_test.groupby(['month','day','hour'])['dw_casetemp'].mean()
test_dw_dometemp = df_test.groupby(['month','day','hour'])['dw_dometemp'].mean()
test_uw_ir = df_test.groupby(['month','day','hour'])['uw_ir'].mean()
test_uw_casetemp = df_test.groupby(['month','day','hour'])['uw_casetemp'].mean()
test_uw_dometemp = df_test.groupby(['month','day','hour'])['uw_dometemp'].mean()
test_uvb = df_test.groupby(['month','day','hour'])['uvb'].mean()
test_par = df_test.groupby(['month','day','hour'])['par'].mean()
test_netsolar = df_test.groupby(['month','day','hour'])['netsolar'].mean()
test_netir = df_test.groupby(['month','day','hour'])['netir'].mean()
test_totalnet = df_test.groupby(['month','day','hour'])['totalnet'].mean()
test_temp = df_test.groupby(['month','day','hour'])['temp'].mean()
test_rh = df_test.groupby(['month','day','hour'])['rh'].mean()
test_windspd = df_test.groupby(['month','day','hour'])['windspd'].mean()
test_winddir = df_test.groupby(['month','day','hour'])['winddir'].mean()
test_pressure = df_test.groupby(['month','day','hour'])['pressure'].mean()
test_ghi = df_test.groupby(['month','day','hour'])['ghi'].mean()
test_Kt = df_test.groupby(['month','day','hour'])['Kt'].mean()


# In[79]:

df_new_test = pd.concat([test_zen,test_dw_solar,test_uw_solar,test_direct_n,test_diffuse,test_dw_ir,
                         test_dw_casetemp,test_dw_dometemp,test_uw_ir,test_uw_casetemp,test_uw_dometemp,
                    test_uvb,test_par,test_netsolar,test_netir,test_totalnet,test_temp,test_rh,
                         test_windspd,test_winddir,test_pressure,test_ghi,test_Kt], axis=1)


# In[178]:

df_new_test.loc[2].xs(17,level='day')


# ### Shifting Kt values to make 1 hour ahead forecast

# #### Train dataset

# In[104]:

levels_index= []
for m in df_new_train.index.levels:
    levels_index.append(m)


# In[105]:

for i in levels_index[0]:
    for j in levels_index[1]:
        df_new_train.loc[i].loc[j]['Kt'] = df_new_train.loc[i].loc[j]['Kt'].shift(-4)


# In[106]:

df_new_train = df_new_train[~(df_new_train['Kt'].isnull())]


# #### Test dataset

# In[83]:

levels_index2= []
for m in df_new_test.index.levels:
    levels_index2.append(m)


# In[84]:

for i in levels_index2[0]:
    for j in levels_index2[1]:
        df_new_test.loc[i].loc[j]['Kt'] = df_new_test.loc[i].loc[j]['Kt'].shift(-4)


# In[85]:

df_new_test = df_new_test[~(df_new_test['Kt'].isnull())]


# In[86]:

df_new_test[df_new_test['Kt']==df_new_test['Kt'].max()]


# ### Normalize train and test dataframe

# In[110]:

train_norm = (df_new_train - df_new_train.mean()) / (df_new_train.max() - df_new_train.min())
test_norm =  (df_new_test - df_new_test.mean()) / (df_new_test.max() - df_new_test.min())


# In[111]:

train_norm.reset_index(inplace=True,drop=True)
test_norm.reset_index(inplace=True,drop=True)


# ### Making train and test sets with train_norm and test_norm

# #### finding the gcf (greatest common factor) of train and test dataset's length and chop off the extra rows to make it divisible with the batchsize

# In[89]:

from fractions import gcd
gcd(train_norm.shape[0],test_norm.shape[0])


# In[112]:

import math
def roundup(x):
    return int(math.ceil(x / 100.0)) * 100


# In[113]:

train_lim = roundup(train_norm.shape[0])
test_lim = roundup(test_norm.shape[0])

train_random = train_norm.sample(train_lim-train_norm.shape[0])
test_random = test_norm.sample(test_lim-test_norm.shape[0])

train_norm = train_norm.append(train_random)
test_norm = test_norm.append(test_random)


# In[114]:

X1 = train_norm.drop('Kt',axis=1)
y1 = train_norm['Kt']

X2 = test_norm.drop('Kt',axis=1)
y2 = test_norm['Kt']


# In[141]:

print("X1_train shape is {}".format(X1.shape))
print("y1_train shape is {}".format(y1.shape))
print("X2_test shape is {}".format(X2.shape))
print("y2_test shape is {}".format(y2.shape))


# In[142]:

X_train = np.array(X1)
y_train  = np.array(y1)
X_test = np.array(X2)
y_test = np.array(y2)


# ### start of RNN

# In[143]:

import torch
import torch.nn as nn
import torchvision.transforms as transforms
import torchvision.datasets as dsets
from torch.autograd import Variable


# In[144]:

class RNNModel(nn.Module):
    def __init__(self, input_dim, hidden_dim, layer_dim, output_dim):
        super(RNNModel, self).__init__()
        #Hidden Dimension
        self.hidden_dim = hidden_dim
        
        # Number of hidden layers
        self.layer_dim = layer_dim
        
        #Building the RNN
        self.rnn = nn.RNN(input_dim, hidden_dim, layer_dim, batch_first=True, nonlinearity='relu')
        
        # Readout layer
        self.fc = nn.Linear(hidden_dim, output_dim)
        
    def forward(self, x):
        # Initializing the hidden state with zeros
        # (layer_dim, batch_size, hidden_dim)
        h0 = Variable(torch.zeros(self.layer_dim, x.size(0), self.hidden_dim))
        
        #One time step (the last one perhaps?)
        out, hn = self.rnn(x, h0)
        
        # Indexing hidden state of the last time step
        # out.size() --> ??
        #out[:,-1,:] --> is it going to be 100,100
        out = self.fc(out[:,-1,:])
        # out.size() --> 100,1
        return out
        


# In[145]:

# Instantiating Model Class
input_dim = 22
hidden_dim = 15
layer_dim = 1
output_dim = 1
batch_size = 100

model = RNNModel(input_dim, hidden_dim, layer_dim, output_dim)

# Instantiating Loss Class
criterion = nn.MSELoss()

# Instantiate Optimizer Class
learning_rate = 0.001
optimizer = torch.optim.SGD(model.parameters(), lr=learning_rate)

# converting numpy array to torch tensor
X_train = torch.from_numpy(X_train)
y_train = torch.from_numpy(y_train)
X_test = torch.from_numpy(X_test)
y_test = torch.from_numpy(y_test)

# initializing lists to store losses over epochs:
train_loss = []
test_loss = []
train_iter = []
test_iter = []


# In[146]:

# Training the model
seq_dim = 1

n_iter =0
num_samples = len(X_train)
test_samples = len(X_test)
batch_size = 100
num_epochs = 1000
feat_dim = X_train.shape[1]

X_train = X_train.type(torch.FloatTensor)
y_train = y_train.type(torch.FloatTensor)
X_test = X_test.type(torch.FloatTensor)
y_test = y_test.type(torch.FloatTensor)

for epoch in range(num_epochs):
    for i in range(0, int(num_samples/batch_size -1)):
        
        
        features = Variable(X_train[i*batch_size:(i+1)*batch_size, :]).view(-1, seq_dim, feat_dim)
        Kt_value = Variable(y_train[i*batch_size:(i+1)*batch_size])
        
        #print("Kt_value={}".format(Kt_value))
        
        optimizer.zero_grad()
        
        outputs = model(features)
        #print("outputs ={}".format(outputs))
        
        loss = criterion(outputs, Kt_value)
        
        train_loss.append(loss.data[0])
        train_iter.append(n_iter)

        #print("loss = {}".format(loss))
        loss.backward()
        
        optimizer.step()
        
        
            
        if n_iter%200 == 0:
            for i in range(0,int(test_samples/batch_size -1)):
                features = Variable(X_test[i*batch_size:(i+1)*batch_size, :]).view(-1, seq_dim, feat_dim)
                Kt_test = Variable(y_test[i*batch_size:(i+1)*batch_size])
                
                outputs = model(features)
                
                mse = np.sqrt(np.mean((Kt_test.data.numpy() - outputs.data.numpy().squeeze())**2)/num_samples)
                
                test_iter.append(n_iter)
                test_loss.append(mse)
            print('Epoch: {} Iteration: {}. Train_MSE: {}. Test_MSE: {}'.format(epoch, n_iter, loss.data[0], mse))
                
        n_iter += 1  
                
               


# In[147]:

print(len(train_loss))
plt.plot(train_loss,'-')
plt.title('train_loss')


# In[148]:

print(len(test_loss))
plt.plot(test_loss,'r')
plt.title('test_loss')


# #### Demornamization

# In[ ]:

# normalized this way:
#test_norm =  (df_new_test - df_new_test.mean()) / (df_new_test.max() - df_new_test.min())


# In[150]:

rmse = np.sqrt(mse)


# In[151]:

rmse_denorm = (rmse * (df_new_test['Kt'].max() - df_new_test['Kt'].min()))+ df_new_test['Kt'].mean()


# In[152]:

print("rmse_denorm",rmse_denorm)


# In[153]:

#df_new_test['Kt'].sort_values()


# In[154]:

print(df_new_test['Kt'].describe())


# ### Saving train and test losses to a csv

# In[149]:

df_trainLoss = pd.DataFrame(data={'Train Loss':train_loss, 'iteration': train_iter}, columns=['Train Loss','iteration'])
df_trainLoss.to_csv('RNN Paper Results/Exp1_2/Boulder_4hour_1000Epoch_TrainLoss.csv')
df_testLoss = pd.DataFrame(data={'Test Loss':test_loss, 'iteration': test_iter}, columns=['Test Loss','iteration'])
df_testLoss.to_csv('RNN Paper Results/Exp1_2/Boulder_4hour_1000EpochTestLoss.csv')


# In[ ]:




# In[ ]:




# In[ ]: