project1.Rmd

---
title: "Project1"
author: "Kumar Aiyer"
date: "01/15/2016"
output: html_document
---

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see <http://rmarkdown.rstudio.com>.

When you click the **Knit** button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

Load and transform the data for the analysis

```{r}
options(width=120)
#
options(warn=-1)

# write a function geteloaddatadf() - you will assign the return value to eloaddf
# in the function do the following
# 1. load the electric load data from elecloaddata.xlsx
# you are not allowed to convert to .csv. Find an appropriate R package that can read .xlsx files and load
# the data in a dataframe called eloaddf. The columns should be dates, kwh, kwVAR
#
# some notes on the data
# the file has 15 min interval data from 1/1/2014 to 12/30/2015
# when you import the data, the first column of your dataset should be dates in POSIXct format
# HINT: use the strptime and as.POSIX.ct functions to form the eloaddf$dates
#

geteloaddatadf<-function(){
library(xlsx)
eloaddf=read.xlsx('C:/Users/jyotika/Desktop/project1-master/elecloaddata.xlsx',sheetName = '15min intervals 2014 & 15')
library(stringr)
temp=paste(str_pad(eloaddf$DATE,6,pad="0"), str_pad(eloaddf$TIME,4,pad="0"), sep=" ")
eloaddf$dates= as.POSIXct(temp, format="%m%d%y %H%M",tz="GMT")
temp_time=eloaddf$TIME
temp_date=eloaddf$DATE
eloaddf$TIME<-NULL
eloaddf$DATE<-NULL
eloaddf<-eloaddf[,c("dates","kWh","kVARh")]
return(eloaddf)
}
# write a function getweatherdf() - you will assign the return value to weatherdf
# 2. Next load the weather data from NOAA into a data frame weatherdf. The data is in 1874606932872dat.txt
# This is 1 hour interval data for a specific weather station close to
# the location of the site from which electric load data was obtained
#

# you need to use fixed width parsing to read the data into a data frame.
# add a column called dates to the dataframe similar to #1 above
#
getweatherdf<-function()
{
weatherdf=read.fwf('C:/Users/jyotika/Desktop/project1-master/1874606932872dat.txt',sep="",widths=c(6,5,12,3,1,3,3,3,1,1,1,4,2,2,2,2,2,2,2,1,2,2,6,5,6,3,3,4,5,5,5,2))
colnames(weatherdf)<-c('USAF','WBAN','YR--MODAHRMN','DIR','SPD','GUS','CLG','SKC','L','M','H','VSB','MW','MW','MW','MW','AW','AW','AW','AW','W','TEMP','DEWP','SLP','ALT')
weatherdf = weatherdf[-1,]
weatherdf$DATES<-as.POSIXct(weatherdf$`YR--MODAHRMN`,format="%Y%m%d%H%M", tz="GMT")
rownames(weatherdf)=1:length(weatherdf$TEMP)
return(weatherdf)
}

# write a funcion getbillsdf() - you will assign the return value to billsdf
# 3. Next load the bill data from billdata.xlsx
# this data is monthly and carefully note the start and end date of each billing period. 
# name the fields of the dataframe as
# billdate, billstartdt, billenddt, kwh, mindemandkw, actualdemandkw, custcharge, 
# distchrgkw, mttkwh, tbckwh,nugckwh, sbckwh, rggieekwh, deliverykwh, 
# totdeliverychrg, supplychrg, totalchrg
#
getbillsdf<-function()
{
library(xlsx)
billsdf=read.xlsx2('C:/Users/jyotika/Desktop/project1-master/billdata.xlsx',1)
colnames(billsdf)<-c('billdate', 'billstartdt', 'billenddt', 'kwh', 'mindemandkw', 'actualdemandkw', 'custcharge', 'distchrgkw', 'mttkwh', 'tbckwh', 'nugckwh', 'sbckwh', 'rggieekwh', 'deliverykwh', 'totdeliverychrg', 'supplychrg', 'totalchrg')
return(billsdf)
}


eloaddf=geteloaddatadf()
billsdf=getbillsdf()
weatherdf=getweatherdf()
```

We now have 3 data sets

1. Electric load data in 15 min interval
2. Weather data in 60 min interval
3. Bill data monthly

Lets do some simple analysis

Display the monthly load profile

```{r}
# display a summary of the electric load data eloaddf$kwh by summarizing it by year, month and total kwh over each month
# your answer should display 24 rows without the header.

temp_kwh=0
indx=1
summary_elecload=matrix(nrow=24,ncol=3,byrow=FALSE)
for(i in 1:(length(eloaddf$dates)-1))
{
temp_kwh=temp_kwh+eloaddf$kWh[i]
if((i==length(eloaddf$dates)-1) || format(eloaddf$dates[i],"%m")!=format(eloaddf$dates[i+1],"%m"))
{
  summary_elecload[indx,1]=as.numeric(format(eloaddf$dates[i],"%Y"))
  summary_elecload[indx,2]=as.numeric(format(eloaddf$dates[i],"%m"))
  summary_elecload[indx,3]=temp_kwh+ eloaddf$kWh[i+1]
  temp_kwh=0
  indx=indx+1
  i=i+1
}
}
print(summary_elecload)

```



Now let us do some plotting of the load data

```{r}
# form a dataframe called eloadhrdf with two columns dates, kwh
# this dataframe sums the 15min kwh in the eloaddf to hourly data
# next create a plot frame with two panels side by side
# On the left panel show a heat map of kwh data for 2014 with x-axis as months and y-axis as hour of the day (1 to 24). use subsetting of the data frame rather than copying the data into yet another data frame
# On the right panel show a heat map of kwh data for 2015 with x-axis as months and y-axis as hour of the day (1 to 24). use subsetting of the data frame rather than copying the data into yet another data frame
arr={}
indx=1
for(i in seq(4,length(eloaddf$dates),4))
{
  arr[indx]= (eloaddf$kWh[i]+eloaddf$kWh[i-1]+eloaddf$kWh[i-2]+eloaddf$kWh[i-3])
  indx=indx+1
}
eloadhrdf=data.frame(eloaddf[seq(4,length(eloaddf$dates),4),1],arr)
colnames(eloadhrdf)=c("dates","kWh")



arr1_temp1=eloadhrdf$kWh[as.numeric(format(eloadhrdf$dates, '%Y'))==2014]
l1=length(arr1_temp1)
arr1_temp1=c(arr1_temp1,eloadhrdf$kWh[l1+1])
tempvar=c(1:23,0)
arr2_temp1=as.numeric(format(eloadhrdf$dates[1:(l1+1)],'%m'))
arr3_temp1=as.numeric(format(eloadhrdf$dates[1:(l1+1)],'%H'))

kwh2014_dm=matrix(nrow=24,ncol=12)
for(i in 1:12){
  for(j in 1:24){
    kwh2014_dm[j,i]=sum(arr1_temp1[arr3_temp1==tempvar[j] & arr2_temp1==i],na.rm=TRUE)
  }
}
arr1_temp2=eloadhrdf$kWh[as.numeric(format(eloadhrdf$dates, '%Y'))==2015]
l2=length(arr1_temp2)
arr1_temp2=arr1_temp2[-1]
arr2_temp2=as.numeric(format(eloadhrdf$dates[(l1+2):(l1+l2)],'%m'))
arr3_temp2=as.numeric(format(eloadhrdf$dates[(l1+2):(l1+l2)],'%H'))
kwh2015_dm=matrix(nrow=24,ncol=12)
for(i in 1:12){
  for(j in 1:24){
    kwh2015_dm[j,i]=sum(arr1_temp2[arr3_temp2==tempvar[j] & arr2_temp2==i],na.rm=TRUE)
  }
}

nf <- layout(matrix(c(1,2),2,2, byrow=T))
image(1:12,1:24,t(kwh2014_dm),col = heat.colors(50),ylab="Hours",xlab="Month")
axis(2, at = seq(1,24,1))
axis(1, at = seq(1,12,1))
title(main= "Heatmap kWh 2014")
image(1:12,1:24,t(kwh2015_dm),col = heat.colors(50),ylab="Hours",xlab="Month")
axis(2, at = seq(1,24,1))
axis(1, at = seq(1,12,1))
title(main= "Heatmap kWh 2015")
```

Usage increased in 2015 around midnight in jan and feb. But overall, 2014 usage seems higher for these 2 months.
May, June, July, August and Sept usage increased in 2015 compared to 2014.
Usage decreased in Oct, Nov and Dec in 2015 compared to 2014, throughout the day.)

Note: BOXPLOT done after next part, since data frames have to be used for this that I created in next part

We are now ready to build a simple predictive model.

```{r}
#create a dataframe with hourly interval data inside your function by 
# combining selective columns from eloadhrdf and weatherdf
# your dataframe should be called modeldatadf and the columns should be dates, year, month, hrofday, temp, kwh
#
#
# write a simple function called predmodel. the model object should be the return parameter
# pass in the appropriate data frames.
# 
# you should fit a GLM with the following specification kwh ~ month + hrofday + temp
# your model should only use 2014 data for your prediction model
#
# use the summary function and display the results of the function

predmodel<-function(day,year,month,hour,temp,kwh){
  modeldatadf=data.frame(day,year,month,hour,temp,kwh)
  return(modeldatadf)
}


library(stats)
wdate=weatherdf$DATES
wtemp=as.numeric(as.character(weatherdf$TEMP))
indx_na=which(is.na(wtemp))
wdate=wdate[-indx_na]
wtemp=wtemp[-indx_na]
wmonth=as.numeric(format(wdate,"%m"))
wday=as.numeric(format(wdate,"%d"))
whour=as.numeric(format(wdate,"%H"))+1
wyear=as.numeric(format(wdate,"%Y"))
wspd=weatherdf$SPD[-indx_na]
wslp=weatherdf$SLP[-indx_na]
wdewp=weatherdf$DEWP[-indx_na]

hour=c(arr3_temp1,arr3_temp2)
day=as.numeric(format(eloadhrdf$dates,"%d"))
month=c(arr2_temp1,arr2_temp2)
year=c(rep(2014,length(arr2_temp1)),rep(2015,length(arr2_temp2)))
LOCS=(hour==0)
hour[LOCS]=24
kwh=c(arr1_temp1,arr1_temp2)
#which(is.na(hour))
for(i in seq(24,length(hour),24)){
    day[i]=day[i-1]
    month[i]=month[i-1]
    year[i]=year[i-1]
}
#new_df=data.frame(hour,day,month,year,kwh)
len_2014=l1+1
len_2015=l2-1

temp=rep(0,length(day))
slp={}
spd={}
dewp={}
temp={}
indx=1
for(yr in c(2014,2015)){
for(i0 in 1:12){
  if(i0==1||i0==3||i0==5||i0==7||i0==8||i0==10)
  daynum=31
  else if(i0==4||i0==6||i0==9||i0==11) 
  daynum=30
  else if(yr==2014 && i0==12)
  daynum=31
  else if(yr==2015 && i0==12)
  daynum=30
  else
  daynum=28
  for(i1 in 1:daynum){
    for(i2 in 1:24){
      temp_var=wtemp[wday==i1 & whour==i2 & wmonth==i0 & wyear==yr]
      tempspd=wspd[wday==i1 & whour==i2 & wmonth==i0 & wyear==yr]
      tempdewp=wdewp[wday==i1 & whour==i2 & wmonth==i0 & wyear==yr]
      tempslp=wslp[wday==i1 & whour==i2 & wmonth==i0 & wyear==yr]
      if(length(temp_var)==0)
      temp_var=0
      if(length(tempdewp)==0)
      tempdewp=0
      if(length(tempspd)==0)
      tempspd=0
      if(length(tempslp)==0)
      tempslp=0
      slp[indx]=tempslp
      spd[indx]=tempspd
      dewp[indx]=tempdewp
      temp[indx]=temp_var
      indx=indx+1
    }
  }
}
}

modeldatadf=predmodel(day,year,month,hour,temp,kwh)

new_df_2014=data.frame(modeldatadf[1:len_2014,])
new_df_2015=data.frame(modeldatadf[(len_2014+1):(len_2014+len_2015),])
rownames(new_df_2015)=1:len_2015
slp2014=slp[1:len_2014]
spd2014=spd[1:len_2014]
dewp2014=dewp[1:len_2014]
slp2015=slp[(len_2014+1):(len_2014+len_2015)]
spd2015=spd[(len_2014+1):(len_2014+len_2015)]
dewp2015=dewp[(len_2014+1):(len_2014+len_2015)]

#predictive model glm
model_glm=glm(new_df_2014$kwh ~ new_df_2014$month + new_df_2014$hour + new_df_2014$temp)
summary(model_glm,correlation=TRUE)

#COMMENTS: Temperature is the most significant feature with highest T value, followed by hour and month, respectively. All features are significant and not correlated with each other. 
```

BOXPLOT
```{r}
# plot the weather data. Use boxplots in ggplot2 with month on the x-axis and temperature in y-axis
library(ggplot2)
library(gridExtra)
library(cowplot)
g1=ggplot(new_df_2014, aes(factor(month), temp))
g2=ggplot(new_df_2015, aes(factor(month), temp))

plot_grid(g1+geom_boxplot(), g2+geom_boxplot(), labels=c("2014", "2015"), ncol = 2, nrow = 1)

#COMMENTS: This plot shows the variation between 2014 n 2015 temperature, different ranges of distribution. There is a hike in July, the lowest mean being in January-February. In general seems like the temp was higher in 2015.
```


Now show you skills in Machine Learning!

```{r}
#
# use the dataframe modeldatadf
# split it into training and testing data sets based on 2014 data for training and 2015 data for testing
# Use the GBM algorithm in the caret package in R to train and validate the model.
# You have free reign to display and explain your results graphically
#
#
library(caret)
library(plyr)
new_df_2014['slp']=slp2014
new_df_2014['spd']=spd2014
new_df_2014['dewp']=dewp2014
new_df_2015['slp']=slp2015
new_df_2015['spd']=spd2015
new_df_2015['dewp']=dewp2015
fitControl <- trainControl( method = "repeatedcv", number = 10, repeats = 10)
tr=train(new_df_2014[,c("hour","month","temp","dewp","day","slp")], new_df_2014[,c("kwh")], method="gbm",trControl=fitControl,verbose=FALSE)
plot(tr)
#filling some missing data before prediction
new_df_2015$temp[8734:8736]=new_df_2015$temp[8735]
p_res2015=predict(tr,new_df_2015[,c("hour","month","temp","dewp","day","slp")])
rmse_gbm=sqrt(mean((p_res2015-new_df_2015$kwh)^2))
print("rmse_gbm")
rmse_gbm
mae=mean(abs(p_res2015-new_df_2015$kwh))
print("mean_abs_error_gbm")
mae
summary(tr)

```

COMMENTS: The RMSE is around ~108, and the mean absolute error is ~83. Initially, the only features used were hour, month and temp. After that, I added more features like DEWP and SPD, which reduced the errors by a slight amount.
The relative inference of features explored was:
 var    rel.inf
hour   38.4390381
temp   25.6440059
month  22.7309565
dewp   6.0247762
day    3.3484539
slp    3.1590483
spd    0.6537213
This indicates, as expected, top features are hour, temp and month, followed by dewp, slp and spd~0.


Lets now compare the predicted model for 2015 with the bill data kwh!

```{r}
#
# run your machine learning model and create a data frame of dates, kwh for 1hr interval data for 2015. note you
# may need to include the last few days of 2014 in your dataset due to the billing dates in January (see billdata.xlsx)
# call your data frame pred2015df.
# now for each of the 12 rows (billing periods) in the billsdf, sum the kwh for the date range in each of the rows from pred2015df for the corresponding start and end of billing in billsdf 
# create a resultsdf which has billdate, predkwh (from pred2015df), actualkwh (from billsdf)
# display the results


startdate=as.Date(as.numeric(as.character(billsdf$billstartdt[-13])), origin='1899-12-30')
enddate=as.Date(as.numeric(as.character(billsdf$billenddt[-13])), origin='1899-12-30')
pred2015df=data.frame(modeldatadf[8473:17496,])
pred2015df$kwh=0
#Using GBM model
fillkwh=predict(tr,pred2015df[c('hour','day','temp')])
pred2015df$kwh=fillkwh
predkwh_bill=rep(0,12)
for(i in 1:12)
{
d1=as.numeric(format(startdate[i],"%d"))
m1=as.numeric(format(startdate[i],"%m"))
y1=as.numeric(format(startdate[i],"%Y"))
d2=as.numeric(format(enddate[i],"%d"))
m2=as.numeric(format(enddate[i],"%m"))
predkwh_bill[i]=sum(pred2015df$kwh[(pred2015df$day>=d1 & pred2015df$month==m1 & pred2015df$year==y1) | (pred2015df$year==2015 & pred2015df$day<=d2 & pred2015df$month==m2) ],na.rm=TRUE)
}
kwh_bill=as.numeric(as.character(billsdf$kwh[1:12]))

# Since bill date is missing for most of the entries, enddate of the bill is taken as bill date
resultsdf = data.frame(startdate,predkwh_bill,kwh_bill)

#displaying result
par(mfrow=c(1,2))
plot(predkwh_bill,type="b")
title(ylab="kwh",xlab="month",main="predicted kwh")
plot(kwh_bill, type="b")
title(ylab="kwh",xlab="month",main="bill kwh")

rmse_billkwh=sqrt(mean((kwh_bill-predkwh_bill)^2))
print("The RMSE for bill kwh: ")
print(rmse_billkwh)


```

COMMENTS: Since the bill data is not accurate, there is high difference in the prediction and bill kwh. Therefore, RMSE is high. Moreover, on looking at the plots of prediction vs bill kwh, Both have a peek around 7th bill, a dip and 2nd bill. Some of the pattern reflects some similarity.

This completes this little exploration of energy load data. Thank You!





CLUSTERING ASSIGNMENT

YEAR 2014
```{r}

#2014

spring14 = subset(new_df_2014, month>=3 & month<=5)
summer14 = subset(new_df_2014, month>=6 & month<=8)
fall14 = subset(new_df_2014, month>=9 & month<=11)
winter14 = subset(new_df_2014, month==12 | month<=2)

# TO GET INDEXING STARTING FROM 1 FOR ROWS
rownames(spring14)<-NULL
rownames(summer14)<-NULL
rownames(fall14)<-NULL
rownames(winter14)<-NULL

# BUILDING MATRIX TO GET HEATMAP OF TEMPERATURE
```

METHOD 1: FINDING SPLIT INTERVALS BY PLOT ANALYSIS

```{r}
# SPRING SEASON ANALYSIS
temp_spring = aggregate(spring14$temp, by = list(hour=spring14$hour), FUN = mean)
plot(temp_spring,ylab = "Avg. Temp")
lines(1:24,rep(46,24))
lines(1:24,rep(55,24))
image(1:24,1, matrix(temp_spring[,2]), xlab= "Hour")
```

FROM THE PLOT AND HEATMAP, I CHOOSE MY INTERVALS AS, (1). HOUR 4-12. (2). HOUR 1-3, 13-16, 24. (3). HOUR 17-23

```{r}
# FROM THE PLOT AND HEATMAP, I CHOOSE MY INTERVALS AS, (1). HOUR 4-12. (2). HOUR 1-3, 13-16, 24. (3). HOUR 17-23
cl14_1 = subset(spring14, hour>=4 & hour<=12)
cl14_2 = subset(spring14, (hour>=1 & hour<=3) | (hour>=13 & hour<=16) | hour==24)
cl14_3 = subset(spring14, hour>=17 & hour<=23)

# SUMMER SEASON ANALYSIS
temp_summer = aggregate(summer14$temp, by = list(hour=summer14$hour), FUN = mean)
plot(temp_summer,ylab = "Avg. Temp")
lines(1:24,rep(65,24))
lines(1:24,rep(73,24))
image(1:24,1, matrix(temp_summer[,2]), xlab= "Hour")
```

FROM THE PLOT AND HEATMAP, I CHOOSE MY INTERVALS AS, (1). HOUR 5-12. (2). HOUR 1-4, 13, 14. (3). HOUR 15-24

```{r}
# FROM THE PLOT AND HEATMAP, I CHOOSE MY INTERVALS AS, (1). HOUR 5-12. (2). HOUR 1-4, 13, 14. (3). HOUR 15-24
cl14_4 = subset(summer14, hour>=5 & hour<=12)
cl14_5 = subset(summer14, (hour>=1 & hour<=4) | hour==13 | hour==14)
cl14_6 = subset(summer14, hour>=15 & hour<=24)

# FALL SEASON ANALYSIS
temp_fall = aggregate(fall14$temp, by = list(hour=fall14$hour), FUN = mean)
plot(temp_fall,ylab = "Avg. Temp")
lines(1:24,rep(55,24))
lines(1:24,rep(60,24))
image(1:24,1, matrix(temp_fall[,2]), xlab= "Hour")
```

FROM THE PLOT AND HEATMAP, I CHOOSE MY INTERVALS AS, (1). HOUR 1-13. (2). HOUR 14, 15, 23, 24. (3). HOUR 16-22

```{r}
# FROM THE PLOT AND HEATMAP, I CHOOSE MY INTERVALS AS, (1). HOUR 1-13. (2). HOUR 14, 15, 23, 24. (3). HOUR 16-22
cl14_7 = subset(fall14, hour>=1 & hour<=13)
cl14_8 = subset(fall14, hour==14 | hour<=15 | hour==23 | hour==24)
cl14_9 = subset(fall14, hour>=16 & hour<=22)

# WINTER SEASON ANALYSIS
temp_winter = aggregate(winter14$temp, by = list(hour=winter14$hour), FUN = mean)
plot(temp_winter,ylab = "Avg. Temp")
lines(1:24,rep(31.5,24))
lines(1:24,rep(37.5,24))
image(1:24,1, matrix(temp_winter[,2]), xlab= "Hour")
```

FROM THE PLOT AND HEATMAP, I CHOOSE MY INTERVALS AS, (1). HOUR 4-13. (2). HOUR 1-3, 14-16, 22-24. (3). HOUR 17-21.

```{r}
# FROM THE PLOT AND HEATMAP, I CHOOSE MY INTERVALS AS, (1). HOUR 4-13. (2). HOUR 1-3, 14-16, 22-24. (3). HOUR 17-21.
cl14_10 = subset(winter14, hour>=4 & hour<=13)
cl14_11 = subset(winter14, (hour>=1 & hour<=3) | (hour>=14 & hour<=16) | (hour>=22 & hour<=24))
cl14_12 = subset(winter14, hour>=17 & hour<=21)

mean_kwh2014 = c(mean(cl14_1$kwh), mean(cl14_2$kwh), mean(cl14_3$kwh), mean(cl14_4$kwh), mean(cl14_5$kwh), mean(cl14_6$kwh), mean(cl14_7$kwh), mean(cl14_8$kwh), mean(cl14_9$kwh), mean(cl14_10$kwh), mean(cl14_11$kwh), mean(cl14_12$kwh))

plot(mean_kwh2014, xlab = 'cluster', main = ' Avg. kwh for the 12 clusters')
```


An analysis for this result can be demonstrated by the following plot

```{r}

plot(winter14$temp[winter14$hour==1])
```
This plot shows the variation of temperature for hour 1 in one season. It ranges between 5-58 degrees. Therefore clustering this way may not give significant distinction for each season based on average temperature. 


METHOD 2: USING K MEANS FOR CLUSTERING EACH SEASON INTO THE THREE TIME INTERVALS

```{r}
##### METHOD 2: USING K MEANS FOR CLUSTERING EACH SEASON INTO THE THREE TIME INTERVALS

# SPRING
kmSPR = kmeans(spring14$temp, 3, iter.max = 30)
spr1_km = spring14$kwh[kmSPR$cluster==1]
spr2_km = spring14$kwh[kmSPR$cluster==2]
spr3_km = spring14$kwh[kmSPR$cluster==3]

# SUMMER
kmSUM = kmeans(summer14$temp, 3, iter.max = 30)
sum1_km = summer14$kwh[kmSUM$cluster==1]
sum2_km = summer14$kwh[kmSUM$cluster==2]
sum3_km = summer14$kwh[kmSUM$cluster==3]

# FALL
kmFAL = kmeans(fall14$temp, 3, iter.max = 30)
fal1_km = fall14$kwh[kmFAL$cluster==1]
fal2_km = fall14$kwh[kmFAL$cluster==2]
fal3_km = fall14$kwh[kmFAL$cluster==3]

# WINTER
kmWIN = kmeans(winter14$temp, 3, iter.max = 30)
win1_km = winter14$kwh[kmWIN$cluster==1]
win2_km = winter14$kwh[kmWIN$cluster==2]
win3_km = winter14$kwh[kmWIN$cluster==3]

mean_kwh2014_km = c(mean(spr1_km), mean(spr2_km), mean(spr3_km), mean(sum1_km), mean(sum2_km), mean(sum3_km), mean(fal1_km), mean(fal2_km), mean(fal3_km), mean(win1_km), mean(win2_km), mean(win3_km))
plot(mean_kwh2014_km, xlab = 'cluster', main = ' Avg. kwh for the 12 clusters using KMeans')
```
As seen in this plot, the variation in average kwh per cluster is better than previously seen, however, still many clashes in cluster allotment would occur because different clusters for different seasons have some overlapping values.
To view this more quantitatively, lets run knn and random forest classifier and find out the 5 fold cross validation accuracy.


```{r}
TR_DATA = c(spr1_km, spr2_km, spr3_km, sum1_km, sum2_km, sum3_km, fal1_km, fal2_km, fal3_km, win1_km, win2_km, win3_km)
TR_OUT = c(rep(1, length(spr1_km)), rep(2, length(spr2_km)), rep(3, length(spr3_km)), rep(4, length(sum1_km)), rep(5, length(sum2_km)), rep(6, length(sum3_km)), rep(7, length(fal1_km)), rep(8, length(fal2_km)), rep(9, length(fal3_km)), rep(10, length(win1_km)), rep(11, length(win2_km)), rep(12, length(win3_km)))
library(KODAMA)
res = KNN.CV(matrix(TR_DATA), (TR_OUT), 1:length(TR_OUT),5)

table(res,TR_OUT)
print("ACCURACY FOR KNN:")
print(length(which(TR_OUT==as.numeric(res)))/length(TR_OUT))
```

The prediction accuracy is around ~20-21%. Depending on application, this maybe a good enough value.
For instance, if we have a group of kwh values which we know would belong to a particular cluster, we can look for maximum overlap in prediction.
 for example, lets take cluster no. 2.
 
```{r}
commlocs = {}
locsreal = which(TR_OUT==2)
for(i in 1:12){
locspred = which(as.numeric(res)==i)
commlocs[i] = 100*length(intersect(locsreal,locspred))/length(locsreal)
}

plot(commlocs, type="h", xlab="cluster", ylab="% overlap with clusters for cluster no.2")
#
```

As we can see, maximum overlap is with cluster 2

```{r}
#random forest
library(randomForest)
rf=randomForest(matrix(TR_DATA),factor(TR_OUT))
print("Accuracy for random forest:")
print(100*length(which(TR_OUT==as.numeric(rf$predicted)))/length(TR_OUT))
# Accuracy is ~18%. As done previously, lets check for cluster no.2
commlocs = {}
locsreal = which(TR_OUT==2)
for(i in 1:12){
locspred = which(rf$predicted==i)
commlocs[i] = 100*length(intersect(locsreal,locspred))/length(locsreal)
}

plot(commlocs, type="h", xlab="cluster", ylab="% overlap with clusters for cluster no.2")

```
Again, maxmimum no. of point belong to cluster 2 in out predcition. Therefore it could be useful in such applications.




YEAR 2015
```{r}
    
#2015

spring15 = subset(new_df_2015, month>=3 & month<=5)
summer15 = subset(new_df_2015, month>=6 & month<=8)
fall15 = subset(new_df_2015, month>=9 & month<=11)
winter15 = subset(new_df_2015, month==12 | month<=2)

# TO GET INDEXING STARTING FROM 1 FOR ROWS
rownames(spring15)<-NULL
rownames(summer15)<-NULL
rownames(fall15)<-NULL
rownames(winter15)<-NULL
```
METHOD 1: SPLIT BY PLOT ANALYSIS
```{r}
# SPRING SEASON ANALYSIS
temp_spring = aggregate(spring15$temp, by = list(hour=spring15$hour), FUN = mean)
plot(temp_spring,ylab = "Avg. Temp")
lines(1:24,rep(45.5,24))
lines(1:24,rep(57,24))
image(1:24,1, matrix(temp_spring[,2]), xlab= "Hour")
```

FROM THE PLOT AND HEATMAP, I CHOOSE MY INTERVALS AS, (1). HOUR 6-12. (2). HOUR 1-5, 13-16, 24. (3). HOUR 17-23

```{r}
# FROM THE PLOT AND HEATMAP, I CHOOSE MY INTERVALS AS, (1). HOUR 6-12. (2). HOUR 1-5, 13-16, 24. (3). HOUR 17-23
cl15_1 = subset(spring15, hour>=6 & hour<=12)
cl15_2 = subset(spring15, (hour>=1 & hour<=5) | (hour>=13 & hour<=16) | hour==24)
cl15_3 = subset(spring15, hour>=17 & hour<=23)

# SUMMER SEASON ANALYSIS
temp_summer = aggregate(summer15$temp, by = list(hour=summer15$hour), FUN = mean)
plot(temp_summer,ylab = "Avg. Temp")
lines(1:24,rep(68.5,24))
lines(1:24,rep(78,24))
image(1:24,1, matrix(temp_summer[,2]), xlab= "Hour")
```

FROM THE PLOT AND HEATMAP, I CHOOSE MY INTERVALS AS, (1). HOUR 6-12. (2). HOUR 1-5, 13-16, 24. (3). HOUR 17-23

```{r}
# FROM THE PLOT AND HEATMAP, I CHOOSE MY INTERVALS AS, (1). HOUR 6-12. (2). HOUR 1-5, 13-16, 24. (3). HOUR 17-23
cl15_4 = subset(summer15, hour>=6 & hour<=12)
cl15_5 = subset(summer15, (hour>=1 & hour<=4) | (hour>=13 & hour<=16) | hour==24)
cl15_6 = subset(summer15, hour>=17 & hour<=23)

# FALL SEASON ANALYSIS
temp_fall = aggregate(fall15$temp, by = list(hour=fall15$hour), FUN = mean)
plot(temp_fall,ylab = "Avg. Temp")
lines(1:24,rep(54,24))
lines(1:24,rep(65,24))
image(1:24,1, matrix(temp_fall[,2]), xlab= "Hour")
```

FROM THE PLOT AND HEATMAP, I CHOOSE MY INTERVALS AS, (1). HOUR 6-12. (2). HOUR 1-5, 13-15, 23, 24. (3). HOUR 16-22

```{r}
# FROM THE PLOT AND HEATMAP, I CHOOSE MY INTERVALS AS, (1). HOUR 6-12. (2). HOUR 1-5, 13-15, 23, 24. (3). HOUR 16-22
cl15_7 = subset(fall15, hour>=6 & hour<=12)
cl15_8 = subset(fall15, hour <=5 | (hour>=14 & hour<=15) | hour==23 | hour==24)
cl15_9 = subset(fall15, hour>=16 & hour<=22)

# WINTER SEASON ANALYSIS
temp_winter = aggregate(winter15$temp, by = list(hour=winter15$hour), FUN = mean)
plot(temp_winter,ylab = "Avg. Temp")
lines(1:24,rep(32,24))
lines(1:24,rep(38.5,24))
image(1:24,1, matrix(temp_winter[,2]), xlab= "Hour")
```

FROM THE PLOT AND HEATMAP, I CHOOSE MY INTERVALS AS, (1). HOUR 4-13. (2). HOUR 1-3, 14-16, 22-24. (3). HOUR 17-21.

```{r}
# FROM THE PLOT AND HEATMAP, I CHOOSE MY INTERVALS AS, (1). HOUR 4-13. (2). HOUR 1-3, 14-16, 22-24. (3). HOUR 17-21.
cl15_10 = subset(winter15, hour>=4 & hour<=13)
cl15_11 = subset(winter15, (hour>=1 & hour<=3) | (hour>=14 & hour<=16) | (hour>=22 & hour<=24))
cl15_12 = subset(winter15, hour>=17 & hour<=21)

mean_kwh2015 = c(mean(cl15_1$kwh), mean(cl15_2$kwh), mean(cl15_3$kwh), mean(cl15_4$kwh), mean(cl15_5$kwh), mean(cl15_6$kwh), mean(cl15_7$kwh), mean(cl15_8$kwh), mean(cl15_9$kwh), mean(cl15_10$kwh), mean(cl15_11$kwh), mean(cl15_12$kwh))

plot(mean_kwh2015, xlab = 'cluster', main = ' Avg. kwh for the 12 clusters')
```

An analysis for this result can be demonstrated by the following plot
```{r}

plot(winter15$temp[winter15$hour==1])
```
This plot shows the variation of temperature for hour 1 in one season. It ranges between 0-60 degrees. Therefore clustering this way may not give significant distinction for each season based on average temperature. 


METHOD 2: USING K MEANS FOR CLUSTERING EACH SEASON INTO THE THREE TIME INTERVALS
```{r}
# SPRING
kmSPR = kmeans(spring15$temp, 3, iter.max = 30)
spr1_km = spring15$kwh[kmSPR$cluster==1]
spr2_km = spring15$kwh[kmSPR$cluster==2]
spr3_km = spring15$kwh[kmSPR$cluster==3]

# SUMMER
kmSUM = kmeans(summer15$temp, 3, iter.max = 30)
sum1_km = summer15$kwh[kmSUM$cluster==1]
sum2_km = summer15$kwh[kmSUM$cluster==2]
sum3_km = summer15$kwh[kmSUM$cluster==3]

# FALL
kmFAL = kmeans(fall15$temp, 3, iter.max = 30)
fal1_km = fall15$kwh[kmFAL$cluster==1]
fal2_km = fall15$kwh[kmFAL$cluster==2]
fal3_km = fall15$kwh[kmFAL$cluster==3]

# WINTER
kmWIN = kmeans(winter15$temp, 3, iter.max = 30)
win1_km = winter15$kwh[kmWIN$cluster==1]
win2_km = winter15$kwh[kmWIN$cluster==2]
win3_km = winter15$kwh[kmWIN$cluster==3]

mean_kwh2015_km = c(mean(spr1_km), mean(spr2_km), mean(spr3_km), mean(sum1_km), mean(sum2_km), mean(sum3_km), mean(fal1_km), mean(fal2_km), mean(fal3_km), mean(win1_km), mean(win2_km), mean(win3_km))
plot(mean_kwh2015_km, xlab = 'cluster', main = ' Avg. kwh for the 12 clusters using KMeans')
```
As seen in this plot, the variation in average kwh per cluster is better than previously seen, however, still many clashes in cluster allotment would occur because different clusters for different seasons have some overlapping values.
To view this more quantitatively, lets run knn and random forest classifier and find out the 5 fold cross validation accuracy.
```{r}
TR_DATA = c(spr1_km, spr2_km, spr3_km, sum1_km, sum2_km, sum3_km, fal1_km, fal2_km, fal3_km, win1_km, win2_km, win3_km)
TR_OUT = c(rep(1, length(spr1_km)), rep(2, length(spr2_km)), rep(3, length(spr3_km)), rep(4, length(sum1_km)), rep(5, length(sum2_km)), rep(6, length(sum3_km)), rep(7, length(fal1_km)), rep(8, length(fal2_km)), rep(9, length(fal3_km)), rep(10, length(win1_km)), rep(11, length(win2_km)), rep(12, length(win3_km)))
library(KODAMA)
res = KNN.CV(matrix(TR_DATA), (TR_OUT), 1:length(TR_OUT),5)

table(res,TR_OUT)

print("ACCURACY FOR KNN:")
print(length(which(TR_OUT==as.numeric(res)))/length(TR_OUT))

```
The prediction accuracy is around ~20%. Depending on application, this maybe a good enough value.
for instance, if we have a group of kwh values which we know would belong to a particular cluster, we can look for maximum overlap in prediction.
for example, lets take cluster no. 1.
```{r}
commlocs = {}


locsreal = which(TR_OUT==1)
for(i in 1:12){
locspred = which(as.numeric(res)==i)
commlocs[i] = 100*length(intersect(locsreal,locspred))/length(locsreal)
}

plot(commlocs, type="h", xlab="cluster", ylab="% overlap with clusters for cluster no.1")
```
As we can see, maximum overlap is with cluster 1
```{r}
#random forest
rf=randomForest(matrix(TR_DATA),factor(TR_OUT))
print("Accuracy for random forest:")
print(100*length(which(TR_OUT==as.numeric(rf$predicted)))/length(TR_OUT))
# Accuracy is ~18%. As done previously, lets check for cluster no.1
commlocs = {}
locsreal = which(TR_OUT==1)
for(i in 1:12){
locspred = which(rf$predicted==i)
commlocs[i] = 100*length(intersect(locsreal,locspred))/length(locsreal)
}

plot(commlocs, type="h", xlab="cluster", ylab="% overlap with clusters for cluster no.1")

```
Again, maxmimum no. of point belong to cluster 1 in out predcition. Therefore it could be useful in such applications.