import pandas as pd
import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error
from math import sqrt

import pprint
import matplotlib.pyplot as plt
import seaborn as sns

from sklearn.model_selection import cross_val_score
from sklearn import metrics

%matplotlib inline

"""
Explanation: Establishing a Baseline for the Problem
Using linear regression and its variants
With various feature combinations
Importing necessary libraries
End of explanation
"""


rice = pd.read_csv("/Users/macbook/Documents/BTP/Notebook/Rice.csv")
rice.head()

rice_haryana = rice[rice["State_Name"]=="Haryana"]
rice_haryana.head()

rainfall = pd.read_csv("/Users/macbook/Documents/BTP/Notebook/rainfall.csv")
rainfall.head()

rain_haryana = rainfall[rainfall["State"]=="Haryana"]
print(rain_haryana.head())
print(rain_haryana.describe())

X_hr = pd.read_csv("/Users/macbook/Documents/BTP/Notebook/haryana.csv")
X_hr.head()

"""
Explanation: Importing the required datasets
End of explanation
"""


X_finite = X_hr[np.isfinite(X_hr["X1"])]
X_finite = X_finite[np.isfinite(X_finite["X2"])]
X_finite = X_finite[np.isfinite(X_finite["X3"])]
X_finite = X_finite[np.isfinite(X_finite["X4"])]
X_finite = X_finite[np.isfinite(X_finite["Y"])]
X_finite.head()

Xn = X_finite
Xn.describe()

y = Xn["Y"]
X = Xn[["X1", "X2", "X3", "X4"]]

plt.figure(figsize=(9, 5))
plt.hist(y, bins=30)
plt.xlabel('Production Value',fontsize=15)
plt.ylabel('Occurences',fontsize=15)
plt.title('Distribution of the Rice Production Values',fontsize=18)

Xplot = Xn[["X1", "X2", "X3", "X4","Y"]]

var_name = "X1"
plt.figure(figsize=(10,6))
sns.regplot(x=var_name, y='Y', data=Xplot, scatter_kws={'alpha':0.6, 's':20})
plt.xlabel(var_name + " (Crop Produce of Last Year)", fontsize=15)
plt.ylabel('Y', fontsize=15)
plt.title("Distribution of y variable with feature "+var_name, fontsize=18)
plt.show()

var_name = "X2"
plt.figure(figsize=(10,6))
sns.regplot(x=var_name, y='Y', data=Xplot, scatter_kws={'alpha':0.6, 's':20})
plt.xlabel(var_name + " (Crop Produce of Last to Last Year)", fontsize=15)
plt.ylabel('Y', fontsize=15)
plt.title("Distribution of y variable with feature "+var_name, fontsize=18)
plt.show()

var_name = "X3"
plt.figure(figsize=(10,6))
sns.regplot(x=var_name, y='Y', data=Xplot, scatter_kws={'alpha':0.6, 's':20})
plt.xlabel(var_name + " (Rainfall of Present Year)", fontsize=15)
plt.ylabel('Y', fontsize=15)
plt.title("Distribution of y variable with feature "+var_name, fontsize=18)
plt.show()

var_name = "X4"
plt.figure(figsize=(10,6))
sns.regplot(x=var_name, y='Y', data=Xplot, scatter_kws={'alpha':0.6, 's':20})
plt.xlabel(var_name + " (Rainfall of Last Year)", fontsize=15)
plt.ylabel('Y', fontsize=15)
plt.title("Distribution of y variable with feature "+var_name, fontsize=18)
plt.show()

# Z-Score Normalization
cols = list(X.columns)

for col in cols:
    col_zscore = col + '_zscore'
    X[col_zscore] = (X[col] - X[col].mean())/X[col].std(ddof=0)
X = X[["X1_zscore", "X2_zscore", "X3_zscore", "X4_zscore"]]

X.head()

"""
Explanation: Cleaning and preparing the datasets
End of explanation
"""


X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=1)

"""
Explanation: Randomly splitting the dataset into training and testing sets
End of explanation
"""


alg = LinearRegression()
alg.fit(X_train, y_train)

coef = alg.coef_
coef = coef.round(decimals=2)
np.set_printoptions(suppress=True) #gem 
print("The coefficients for the linear regression model learnt are\n")
print(coef)
print()

y_predict = alg.predict(X_test)
rmse = sqrt(mean_squared_error(y_predict, y_test))
print(rmse)

"""
Explanation: First baseline using Linear Regression
End of explanation
"""


clf = LinearRegression()
scores = cross_val_score(clf, X, y, cv=5, scoring='neg_mean_squared_error')

for i in range(0,5):
    scores[i] = sqrt(-1*scores[i])

print(scores)
avg_rmse = scores.mean()
print("\n\nAvg RMSE is ",scores.mean())

# print(type(y_test))
# print(type(y_predict))
yt = y_test.as_matrix()
print(type(yt))

p = pd.DataFrame()
p["y_predicted"] = y_predict/1000
p["y_test"] = yt/1000

p["y_predicted"] = p["y_predicted"].round(decimals=1)
# p["y_test"] = p["y_test"].round(decimals=1)
p.describe()

"""
Explanation: Lets calculate the average RMSE (Cross Validation, 5-Fold)
End of explanation
"""


print(p)

"""
Explanation: Good enough results for the Haryana State
End of explanation
"""


rmse/1000

"""
Explanation: The Root Mean Square Error
It has the same unit as the data values.
With range of the test data set being [8,460], rmse as 23.3 is a decent one.
End of explanation
"""


rain1 = rainfall
rain2 = pd.read_csv("/Users/macbook/Documents/BTP/Notebook/rainfall_distt_2004-10_nax.csv")

print(rice.describe())
print(rain1.describe())
print(rain2.describe())

a = np.empty((rice.shape[0],1))*np.NAN
rice = rice.assign(X1 = a)
rice = rice.assign(X2 = a)
rice = rice.assign(X3 = a)
rice = rice.assign(X4 = a)
rice.head()

"""
Explanation: Constructing dataset for Whole India
Preparing columns for the features, using rice production and rainfall datasets
End of explanation
"""


l = rice.shape[0]
for row in range(0,l):
    if row-1<0 or rice.iloc[row,1] != rice.iloc[row-1,1]:
        continue
    else:
        rice.iloc[row,8] = rice.iloc[row-1,6]
        if row-2<0 or rice.iloc[row,1] != rice.iloc[row-2,1]:
            continue
        else:
            rice.iloc[row,9] = rice.iloc[row-2,6]

rice.head()

def func(s):
    x = s.strip()
    return x.lower()

rice['ind_district'] = rice['ind_district'].apply(func)
rice['Season'] = rice['Season'].apply(func)
rain1['ind_district'] = rain1['ind_district'].apply(func)
rain2['ind_district'] = rain2['ind_district'].apply(func)
rice.head()

rain1.head()

# can reduce the time by searching only one variable for some cases atleast
rice = rice[rice['Season'] == 'kharif']
l = rice.shape[0]

for row in range(0,l):
    
    dt = rice.iloc[row,1]
    yr = rice.iloc[row,2]
    
    if yr<=2002:
        
        # rainfall for the same year
        r = rain1.loc[(rain1.ind_district == dt) & (rain1.Year == yr)]
        if r.shape[0] == 1:
            rice.iloc[row,10] = r.iloc[0,3]
            
        # rainfall for the previous year
        r = rain1.loc[(rain1.ind_district == dt) & (rain1.Year == yr-1)]
        if r.shape[0] == 1:
            rice.iloc[row,11] = r.iloc[0,3]
            
    if yr>2004:
        
        # rainfall for the same year
        r = rain2.loc[(rain2.ind_district == dt) & (rain2.Year == yr)]
        if r.shape[0] == 1:
            rice.iloc[row,10] = r.iloc[0,3]
            
        # rainfall for the previous year
        r = rain2.loc[(rain2.ind_district == dt) & (rain2.Year == yr-1)]
        if r.shape[0] == 1:
            rice.iloc[row,11] = r.iloc[0,3]
        
    

# X1 = prod-1
# X2 = prod-2
# X3 = rain
# X4 = rain-1
rice.describe()

ricex = rice[np.isfinite(rice["Production"])]
ricex = ricex[np.isfinite(ricex["X1"])]
ricex = ricex[np.isfinite(ricex["X2"])]
ricex = ricex[np.isfinite(ricex["X3"])]
ricex = ricex[np.isfinite(ricex["X4"])]
ricex.head()

X = ricex[["X1","X2","X3","X4"]]
y = ricex[["Production"]]
ricex.describe()

plt.figure(figsize=(30, 10))
plt.hist(y, bins=250)
plt.xlabel('Production Value',fontsize=25)
plt.ylabel('Occurences',fontsize=25)
plt.title('Distribution of the Rice Production Values',fontsize=30)

Xplot = ricex[["X1", "X2", "X3", "X4","Production"]]

var_name = "X1"
plt.figure(figsize=(10,6))
sns.regplot(x=var_name, y='Production', data=Xplot, scatter_kws={'alpha':0.6, 's':20})
plt.xlabel(var_name + " (Crop Produce of Last Year)", fontsize=15)
plt.ylabel('Y', fontsize=15)
plt.title("Distribution of y variable with feature "+var_name, fontsize=18)
plt.show()

Xplot = ricex[["X1", "X2", "X3", "X4","Production"]]

var_name = "X2"
plt.figure(figsize=(10,6))
sns.regplot(x=var_name, y='Production', data=Xplot, scatter_kws={'alpha':0.6, 's':20})
plt.xlabel(var_name + " (Crop Produce of Last to Last Year)", fontsize=15)
plt.ylabel('Y', fontsize=15)
plt.title("Distribution of y variable with feature "+var_name, fontsize=18)
plt.show()

Xplot = ricex[["X1", "X2", "X3", "X4","Production"]]

var_name = "X3"
plt.figure(figsize=(10,6))
sns.regplot(x=var_name, y='Production', data=Xplot, scatter_kws={'alpha':0.6, 's':20})
plt.xlabel(var_name + " (Rainfall of Present Year)", fontsize=15)
plt.ylabel('Y', fontsize=15)
plt.title("Distribution of y variable with feature "+var_name, fontsize=18)
plt.show()

Xplot = ricex[["X1", "X2", "X3", "X4","Production"]]

var_name = "X4"
plt.figure(figsize=(10,6))
sns.regplot(x=var_name, y='Production', data=Xplot, scatter_kws={'alpha':0.6, 's':20})
plt.xlabel(var_name + " (Rainfall of Last Year)", fontsize=15)
plt.ylabel('Y', fontsize=15)
plt.title("Distribution of y variable with feature "+var_name, fontsize=18)
plt.show()

# Z-Score Normalization

cols = list(X.columns)
for col in cols:
    col_zscore = col + '_zscore'
    X[col_zscore] = (X[col] - X[col].mean())/X[col].std(ddof=0)

X = X[["X1_zscore", "X2_zscore", "X3_zscore", "X4_zscore"]]

X.head()

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=1)

alg = LinearRegression()
alg.fit(X_train, y_train)

coef = alg.coef_
intercept = alg.intercept_

coef = coef.round(decimals=2)
pp = pprint.PrettyPrinter()
pp.pprint(coef)
pp.pprint(intercept)

y_predict = alg.predict(X_test)

yp = y_predict
yt = y_test.as_matrix()
type(y_predict)

rmse = sqrt(mean_squared_error(y_predict, y_test))
print(rmse)

"""
Explanation: Constructing the features X1 and X2, the production for the last 2 years.
End of explanation
"""


clf = LinearRegression()
scores = cross_val_score(clf, X, y, cv=5, scoring='neg_mean_squared_error')

for i in range(0,5):
    scores[i] = sqrt(-1*scores[i])

print(scores)
avg_rmse = scores.mean()
print("\n\nAvg RMSE is ",scores.mean())

yt = yt/1000
yp = yp/1000
yt = yt.round(decimals=1)
yp = yp.round(decimals=1)

yo = np.concatenate((yp,yt),axis=1)

p = pd.DataFrame(data=yo,columns=['Predicted','Actual'])
p.describe()

"""
Explanation: Lets calculate the average RMSE (Cross Validation, 5-Fold)
End of explanation
"""


p

"""
Explanation: Result
This time the predictions are not as good as we got for just one state.
End of explanation
"""


rmse/1000

"""
Explanation: The range of the values is (0,1362) and rmse is 70.6
We can say the results are decent but not good enough as we got with just one state.
End of explanation
"""


# loading the dataset
ricep = pd.read_csv("/Users/macbook/Documents/BTP/Notebook/rice with soil(P).csv")

# Removing the rows with missing value of phosphorus
ricep = ricep[np.isfinite(ricep["phosphorus"])]

# Adding collumns for the other 4 features
a = np.empty((ricep.shape[0],1))*np.NAN
ricep = ricep.assign(X1 = a)
ricep = ricep.assign(X2 = a)
ricep = ricep.assign(X3 = a)
ricep = ricep.assign(X4 = a)
ricep.head()

# Constructing features X1 and X2
l = ricep.shape[0]
for row in range(0,l):
    if row-1<0 or ricep.iloc[row,1] != ricep.iloc[row-1,1]:
        continue
    else:
        ricep.iloc[row,8] = ricep.iloc[row-1,6]
        if row-2<0 or ricep.iloc[row,1] != ricep.iloc[row-2,1]:
            continue
        else:
            ricep.iloc[row,9] = ricep.iloc[row-2,6]

# Making the strings in the dataset uniform, with other datasets
ricep['ind_district'] = ricep['ind_district'].apply(func)
ricep['Season'] = ricep['Season'].apply(func)

# Constructing features X3 and X4
l = ricep.shape[0]

for row in range(0,l):
    
    dt = ricep.iloc[row,1]
    yr = ricep.iloc[row,2]
    
    if yr<=2002:
        
        # rainfall for the same year
        r = rain1.loc[(rain1.ind_district == dt) & (rain1.Year == yr)]
        if r.shape[0] == 1:
            ricep.iloc[row,10] = r.iloc[0,3]
            
        # rainfall for the previous year
        r = rain1.loc[(rain1.ind_district == dt) & (rain1.Year == yr-1)]
        if r.shape[0] == 1:
            ricep.iloc[row,11] = r.iloc[0,3]
            
    if yr>2004:
        
        # rainfall for the same year
        r = rain2.loc[(rain2.ind_district == dt) & (rain2.Year == yr)]
        if r.shape[0] == 1:
            ricep.iloc[row,10] = r.iloc[0,3]
            
        # rainfall for the previous year
        r = rain2.loc[(rain2.ind_district == dt) & (rain2.Year == yr-1)]
        if r.shape[0] == 1:
            ricep.iloc[row,11] = r.iloc[0,3]

# Removing rows with any missing values
ricep = ricep[np.isfinite(ricep["Production"])]
ricep = ricep[np.isfinite(ricep["X1"])]
ricep = ricep[np.isfinite(ricep["X2"])]
ricep = ricep[np.isfinite(ricep["X3"])]
ricep = ricep[np.isfinite(ricep["X4"])]
ricep.head()

ricep['phosphorus'] = ricep['phosphorus'] + 1
ricep.describe()

X = ricep[["X1","X2","X3","X4","phosphorus"]]
y = ricep[["Production"]]
ricep.to_csv("ricep.csv")

Xplot = ricep[["X1", "X2", "X3", "X4", "phosphorus", "Production"]]

var_name = "phosphorus"
plt.figure(figsize=(10,6))
sns.regplot(x=var_name, y='Production', data=Xplot, scatter_kws={'alpha':0.6, 's':20})
plt.xlabel(var_name + " (in Soil)", fontsize=15)
plt.ylabel('Y', fontsize=15)
plt.title("Distribution of y variable with feature "+var_name, fontsize=18)
plt.show()

# Z-Score Normalization

cols = list(X.columns)
for col in cols:
    col_zscore = col + '_zscore'
    X[col_zscore] = (X[col] - X[col].mean())/X[col].std(ddof=0)

X = X[["X1_zscore", "X2_zscore", "X3_zscore", "X4_zscore", "phosphorus_zscore"]]
X.head()

X.describe()

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=1)

alg = LinearRegression()
alg.fit(X_train, y_train)

coef = alg.coef_
intercept = alg.intercept_

coef = coef.round(decimals=2)
pp = pprint.PrettyPrinter()
pp.pprint(coef)
pp.pprint(intercept)

y_predict = alg.predict(X_test)

yp = y_predict
yt = y_test.as_matrix()
type(y_predict)

rmse = sqrt(mean_squared_error(y_predict, y_test))
print(rmse)

"""
Explanation: Adding a new feature
Quantity of Phosphorous in soil (low, medium, high)
This feaure is missing for some states in India, so little less number of samples
End of explanation
"""


clf = LinearRegression()
scores = cross_val_score(clf, X, y, cv=5, scoring='neg_mean_squared_error')

for i in range(0,5):
    scores[i] = sqrt(-1*scores[i])

print(scores)
avg_rmse_phos = scores.mean()
print("\n\nAvg RMSE is ",scores.mean())

yt = yt/1000
yp = yp/1000
yt = yt.round(decimals=1)
yp = yp.round(decimals=1)

yo = np.concatenate((yp,yt),axis=1)

p = pd.DataFrame(data=yo,columns=['Predicted','Actual'])
p.describe()

p

rmse/1000

"""
Explanation: Lets calculate the average RMSE (Cross Validation, 5-Fold)
End of explanation
"""


# Just the 4 original features (no soil data)
X_old = X[["X1_zscore", "X2_zscore", "X3_zscore", "X4_zscore"]]

# Seed is fixed, so the vector y_test is going to same as before
X_train, X_test, y_train, y_test = train_test_split(X_old, y, test_size=0.2, random_state=1)

alg = LinearRegression()
alg.fit(X_train, y_train)

coef = alg.coef_
intercept = alg.intercept_

coef = coef.round(decimals=2)
pp = pprint.PrettyPrinter()
pp.pprint(coef)
pp.pprint(intercept)

y_predict = alg.predict(X_test)
yp = y_predict
yt = y_test.as_matrix()
rmse = sqrt(mean_squared_error(y_predict, y_test))
print(rmse)

"""
Explanation: Now lets compare with other feature combinations
End of explanation
"""


clf = LinearRegression()
scores = cross_val_score(clf, X_old, y, cv=5, scoring='neg_mean_squared_error')

for i in range(0,5):
    scores[i] = sqrt(-1*scores[i])

print(scores)
avg_rmse_orig = scores.mean()
print("\n\nAvg RMSE is ",scores.mean())

"""
Explanation: Average RMSE (Cross Validation, 5-Fold)
End of explanation
"""


X_no_rain = X[["X1_zscore", "X2_zscore"]]

X_train, X_test, y_train, y_test = train_test_split(X_no_rain, y, test_size=0.2, random_state=1)

alg = LinearRegression()
alg.fit(X_train, y_train)

coef = alg.coef_
intercept = alg.intercept_

coef = coef.round(decimals=2)
pp = pprint.PrettyPrinter()
pp.pprint(coef)
pp.pprint(intercept)

y_predict = alg.predict(X_test)
yp = y_predict
yt = y_test.as_matrix()
rmse = sqrt(mean_squared_error(y_predict, y_test))
print(rmse)

"""
Explanation: Avg RMSE with original 4 features  :  51874.9                                             
Avg RMSE with 5 features(soil too) : 52272.6
So there was NO improvement in the RMSE, due to the new feature.
Now lets try by removing the rainfall features
End of explanation
"""


clf = LinearRegression()
scores = cross_val_score(clf, X_no_rain, y, cv=5, scoring='neg_mean_squared_error')

for i in range(0,5):
    scores[i] = sqrt(-1*scores[i])

print(scores)
avg_rmse_no_rain = scores.mean()
print("\n\nAvg RMSE is ",scores.mean())

"""
Explanation: Avg RMSE
End of explanation
"""


from sklearn import linear_model

reg = linear_model.RidgeCV(alphas=[1,2,3,4,5,6,7,7.1,7.2,7.3,8,9,10])
reg.fit(X_old, y)
reg.alpha_

X_train, X_test, y_train, y_test = train_test_split(X_old, y, test_size=0.2, random_state=1)

reg = linear_model.Ridge(alpha = 7.1)
reg.fit (X_train, y_train)
print(reg.coef_)

y_predict = reg.predict(X_test)
rmse = sqrt(mean_squared_error(y_predict, y_test))
print(rmse)

"""
Explanation: Avg RMSE with the original 4 features : 51874.9
Avg RMSE without the 2  rain features : 51921.9
Ridge Regression
End of explanation
"""


clf = linear_model.Ridge(alpha = 7.1)
scores = cross_val_score(clf, X_old, y, cv=5, scoring='neg_mean_squared_error')

for i in range(0,5):
    scores[i] = sqrt(-1*scores[i])
    
print(scores)
avg_rmse_ridge = scores.mean()
print("\n\nAvg RMSE is ",scores.mean())

"""
Explanation: Avg RMSE
End of explanation
"""


from sklearn import linear_model

reg = linear_model.LassoCV(alphas=[0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1])
reg.fit(X_old, y)
reg.alpha_

X_train, X_test, y_train, y_test = train_test_split(X_old, y, test_size=0.2, random_state=1)

las = linear_model.Lasso(alpha = 0.01)
las.fit (X_train, y_train)
print(las.coef_)

y_predict = las.predict(X_test)
rmse = sqrt(mean_squared_error(y_predict, y_test))
print(rmse)

"""
Explanation: Avg RMSE with Linear Regression : 51874.9
Avg RMSE with Ridge Regression : 51840.9
So Ridge is slightly better than Linear
Lasso Regression
End of explanation
"""


clf = linear_model.Lasso(alpha = 0.01)
scores = cross_val_score(clf, X_old, y, cv=5, scoring='neg_mean_squared_error')

for i in range(0,5):
    scores[i] = sqrt(-1*scores[i])
    
print(scores)
avg_rmse_las = scores.mean()
print("\n\nAvg RMSE is ",scores.mean())

"""
Explanation: Avg RMSE
End of explanation
"""


