# ---------------------------------------------------------------------------
# crop_yield.py
#
# Capstone project: predicting the yield of a farm plot (tonnes per hectare)
# from its crop, weather, soil and the fertiliser used.
#
# CBSE Class XII, Artificial Intelligence (843), Part C - Capstone Project.
# SDG 2 (Zero Hunger) and SDG 12 (Responsible Consumption and Production).
#
# The program follows the 10 steps of the Data Science Methodology:
#   1 Business Understanding   2 Analytic Approach   3 Data Requirements
#   4 Data Collection          5 Data Understanding  6 Data Preparation
#   7 Modelling                8 Evaluation          9 Deployment
#   10 Feedback
#
# Two regression models are compared - Linear Regression, as in the
# handbook, and a Random Forest - with a train-test split and with 5-fold
# cross validation, using MAE, MSE, RMSE and R squared.
#
# Libraries: pandas, numpy, matplotlib, scikit-learn, joblib.
# ---------------------------------------------------------------------------

# %% Step 1 - Business Understanding
print("=" * 64)
print("HOW MUCH WILL THIS FIELD YIELD?")
print("=" * 64)
print("Problem : farmers and planners cannot tell before harvest how much")
print("          a plot will yield, so input use and storage are guesswork.")
print("Goal    : predict yield (t/ha) from crop, weather, soil and inputs.")
print()

# %% Step 2 - Analytic Approach
# The answer is a number (tonnes per hectare), so this is PREDICTIVE
# analytics using REGRESSION.
print("Analytic approach : predictive - regression (the target is a number)")
print()

# %% Step 3 - Data Requirements
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split, KFold, cross_val_score
from sklearn.linear_model import LinearRegression
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
import joblib

pd.set_option("display.width", 130)
pd.set_option("display.max_columns", 20)

required = ["crop", "rainfall_mm", "temperature_c", "soil_ph", "nitrogen_kg_ha",
            "fertiliser_kg_ha", "irrigated", "yield_t_ha"]
print("Columns required :", required)

# %% Step 4 - Data Collection
data = pd.read_csv("crop_yield.csv")
print("Rows collected   :", len(data))
missing_columns = [c for c in required if c not in data.columns]
print("Required columns missing from the file :", missing_columns)
print()
print(data.head(8))
print()

# %% Step 5 - Data Understanding
print("--- Step 5: Data Understanding ---")
print("Missing values in each column:")
print(data.isnull().sum().to_string())
print()
print("Summary statistics:")
print(data.describe().round(2))
print()

# Chart 1: the distribution of yield.
plt.figure(figsize=(7, 4))
plt.hist(data["yield_t_ha"].clip(0, 10), bins=30, color="#59a14f", edgecolor="white")
plt.title("Distribution of yield (values above 10 shown at 10)")
plt.xlabel("Yield (tonnes per hectare)")
plt.ylabel("Number of plots")
plt.tight_layout()
plt.savefig("chart1_yield_distribution.png")
plt.close()

# Chart 2: fertiliser against yield.
ok = data[(data["yield_t_ha"] > 0) & (data["yield_t_ha"] < 15)]
plt.figure(figsize=(7, 5))
plt.scatter(ok["fertiliser_kg_ha"], ok["yield_t_ha"], s=12, alpha=0.5, color="#4e79a7")
plt.title("Fertiliser used and yield")
plt.xlabel("Fertiliser (kg per hectare)")
plt.ylabel("Yield (t/ha)")
plt.grid(True, linestyle="--", alpha=0.4)
plt.tight_layout()
plt.savefig("chart2_fertiliser_yield.png")
plt.close()

# Chart 3: yield of each crop, as box plots.
plt.figure(figsize=(7, 4.5))
crops_seen = sorted(ok["crop"].str.strip().str.title().unique())
plt.boxplot([ok[ok["crop"].str.strip().str.title() == c]["yield_t_ha"] for c in crops_seen], labels=crops_seen)
plt.title("Yield of each crop")
plt.ylabel("Yield (t/ha)")
plt.grid(True, axis="y", linestyle="--", alpha=0.4)
plt.tight_layout()
plt.savefig("chart3_yield_by_crop.png")
plt.close()
for c in crops_seen:
    part = ok[ok["crop"].str.strip().str.title() == c]["yield_t_ha"]
    print("Median yield of " + c.ljust(6), ":", round(part.median(), 2), "t/ha")
print()

# %% Step 6 - Data Preparation
print("--- Step 6: Data Preparation ---")
before = len(data)
data = data.drop_duplicates()
print("Duplicate rows removed          :", before - len(data))

data["crop"] = data["crop"].str.strip().str.title()
print("Crops after tidying the names   :", sorted(data["crop"].unique()))

impossible = ((data["yield_t_ha"] <= 0) | (data["yield_t_ha"] > 15)).sum()
data = data[(data["yield_t_ha"] > 0) & (data["yield_t_ha"] <= 15)]
print("Impossible yields removed       :", impossible)

# Missing rainfall is filled with the median rainfall of the same crop:
# rice fields get far more rain than wheat fields, so one overall median
# would be wrong for both.
blank_rain = data["rainfall_mm"].isnull().sum()
data["rainfall_mm"] = data.groupby("crop")["rainfall_mm"].transform(lambda s: s.fillna(s.median()))
print("Blank rainfall values filled    :", blank_rain, "(median of the same crop)")
print("Rows ready for modelling        :", len(data))

# Text columns become numbers. irrigated: Yes = 1, No = 0. crop: one
# column per crop (one-hot encoding), so no crop is treated as "bigger".
data["irrigated"] = (data["irrigated"] == "Yes").astype(int)
data = pd.get_dummies(data, columns=["crop"], dtype=int)
features = ["rainfall_mm", "temperature_c", "soil_ph", "nitrogen_kg_ha", "fertiliser_kg_ha",
            "irrigated", "crop_Maize", "crop_Rice", "crop_Wheat"]
print("Features used                   :", len(features))
print()

# Chart 4: correlation of every feature with yield.
corr = data[features + ["yield_t_ha"]].corr()["yield_t_ha"].drop("yield_t_ha").sort_values()
print("Correlation of each feature with yield:")
print(corr.round(2).to_string())
print()
plt.figure(figsize=(7, 4.5))
colours = ["#e15759" if v < 0 else "#59a14f" for v in corr.values]
plt.barh(corr.index, corr.values, color=colours)
plt.axvline(0, color="#999999")
plt.title("Correlation with yield")
plt.xlabel("Correlation coefficient")
plt.tight_layout()
plt.savefig("chart4_correlation.png")
plt.close()

# %% Step 7 - Modelling
print("--- Step 7: Modelling ---")
X = data[features]
y = data["yield_t_ha"]
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, shuffle=True, random_state=7)
print("Training rows :", len(X_train))
print("Testing rows  :", len(X_test))

linear = LinearRegression()
linear.fit(X_train, y_train)
forest = RandomForestRegressor(n_estimators=200, min_samples_leaf=3, random_state=0)
forest.fit(X_train, y_train)

print("Linear Regression - intercept :", round(linear.intercept_, 3))
print("Linear Regression - fertiliser coefficient :",
      round(linear.coef_[features.index("fertiliser_kg_ha")], 4), "t/ha per kg")
print()

# %% Step 8 - Evaluation
print("--- Step 8: Evaluation ---")
results = []
for name, model in [["Linear Regression", linear], ["Random Forest", forest]]:
    pred = model.predict(X_test)
    mae = mean_absolute_error(y_test, pred)
    mse = mean_squared_error(y_test, pred)
    rmse = np.sqrt(mse)
    r2 = r2_score(y_test, pred)
    results.append([name, mae, mse, rmse, r2])
    print(name)
    print("  " + name + " MAE  :", round(mae, 3), "t/ha")
    print("  " + name + " MSE  :", round(mse, 3))
    print("  " + name + " RMSE :", round(rmse, 3), "t/ha")
    print("  " + name + " R2   :", round(r2, 3))
print()

# 5-fold cross validation: every row is in the test fold exactly once.
folds = KFold(n_splits=5, shuffle=True, random_state=1)
cv_lin = np.sqrt(-cross_val_score(LinearRegression(), X, y, cv=folds, scoring="neg_mean_squared_error"))
cv_rf = np.sqrt(-cross_val_score(RandomForestRegressor(n_estimators=200, min_samples_leaf=3, random_state=0),
                                 X, y, cv=folds, scoring="neg_mean_squared_error"))
print("RMSE in each of the 5 folds (t/ha):")
print("  Fold  Linear  Forest")
for k in range(5):
    print("  ", k + 1, "  ", round(cv_lin[k], 3), " ", round(cv_rf[k], 3))
print("Mean cross-validated RMSE, Linear Regression :", round(cv_lin.mean(), 3))
print("Mean cross-validated RMSE, Random Forest     :", round(cv_rf.mean(), 3))
print()

best_pred = forest.predict(X_test)

# Chart 5: predicted against actual, for both models side by side.
lin_pred = linear.predict(X_test)
plt.figure(figsize=(10, 4.8))
for k in range(2):
    plt.subplot(1, 2, k + 1)
    pred = [lin_pred, best_pred][k]
    plt.scatter(y_test, pred, s=14, alpha=0.6, color=["#9c9c9c", "#59a14f"][k])
    plt.plot([1, 9], [1, 9], color="#e15759", linestyle="--", label="Perfect prediction")
    plt.title(["Linear Regression", "Random Forest"][k])
    plt.xlabel("Actual yield (t/ha)")
    plt.ylabel("Predicted yield (t/ha)")
    plt.legend()
    plt.grid(True, linestyle="--", alpha=0.4)
plt.tight_layout()
plt.savefig("chart5_predicted_actual.png")
plt.close()

# Chart 6: the residuals (actual minus predicted) of the forest.
residuals = y_test - best_pred
print("Mean residual (Random Forest)        :", round(residuals.mean(), 3), "t/ha")
print("Residuals within +/- 0.5 t/ha        :", round((residuals.abs() <= 0.5).mean() * 100, 1), "%")
plt.figure(figsize=(7, 4))
plt.hist(residuals, bins=25, color="#59a14f", edgecolor="white")
plt.axvline(0, color="#e15759", linestyle="--")
plt.title("Random Forest: residuals on the test plots")
plt.xlabel("Actual minus predicted yield (t/ha)")
plt.ylabel("Number of plots")
plt.tight_layout()
plt.savefig("chart6_residuals.png")
plt.close()

# Chart 7: RMSE in each fold for both models.
x = np.arange(5)
plt.figure(figsize=(7, 4))
plt.bar(x - 0.2, cv_lin, width=0.4, label="Linear Regression", color="#9c9c9c")
plt.bar(x + 0.2, cv_rf, width=0.4, label="Random Forest", color="#59a14f")
plt.xticks(x, ["Fold 1", "Fold 2", "Fold 3", "Fold 4", "Fold 5"])
plt.ylabel("RMSE (t/ha) - lower is better")
plt.title("5-fold cross validation")
plt.legend()
plt.tight_layout()
plt.savefig("chart7_kfold.png")
plt.close()

# Chart 8: which features the forest relied on.
importance = pd.Series(forest.feature_importances_, index=features).sort_values()
print("Feature importance (Random Forest):")
print(importance.round(3).to_string())
print()
plt.figure(figsize=(7, 4.5))
plt.barh(importance.index, importance.values, color="#f28e2b")
plt.title("What the Random Forest relied on")
plt.xlabel("Importance")
plt.tight_layout()
plt.savefig("chart8_importance.png")
plt.close()

# %% Step 9 - Deployment
print("--- Step 9: Deployment ---")
# The trained model is saved to a file, so an app or website can load it
# and predict without training again.
joblib.dump(forest, "yield_model.pkl")
loaded = joblib.load("yield_model.pkl")
print("Model saved to yield_model.pkl and loaded back")

new_plots = pd.DataFrame({
    "rainfall_mm": [1180, 560, 640], "temperature_c": [26.5, 18.5, 27.0], "soil_ph": [6.7, 7.1, 7.9],
    "nitrogen_kg_ha": [85, 90, 60], "fertiliser_kg_ha": [130, 150, 40], "irrigated": [0, 1, 0],
    "crop_Maize": [0, 0, 1], "crop_Rice": [1, 0, 0], "crop_Wheat": [0, 1, 0],
})
names = ["Rice, rain-fed, good rain", "Wheat, irrigated", "Maize, rain-fed, little fertiliser"]
predicted = loaded.predict(new_plots)
for i in range(len(names)):
    print(names[i].ljust(36), "->", round(predicted[i], 2), "t/ha")
print()

# %% Step 10 - Feedback
print("--- Step 10: Feedback ---")
# After the next season's harvest, the real yields come back. Comparing
# them with what the model predicted shows whether it still works.
nxt = pd.read_csv("next_season.csv")
nxt["irrigated"] = (nxt["irrigated"] == "Yes").astype(int)
nxt = pd.get_dummies(nxt, columns=["crop"], dtype=int)
next_pred = loaded.predict(nxt[features])
next_rmse = np.sqrt(mean_squared_error(nxt["yield_t_ha"], next_pred))
print("Plots in the next season's feedback :", len(nxt))
print("RMSE on last season's test data     :", round(results[1][3], 3), "t/ha")
print("RMSE on the next season             :", round(next_rmse, 3), "t/ha")
print("Average rainfall, training season   :", round(data["rainfall_mm"].mean()), "mm")
print("Average rainfall, next season       :", round(nxt["rainfall_mm"].mean()), "mm")
if next_rmse > 1.2 * results[1][3]:
    print("Feedback: error has grown by more than 20% - retrain with the new season's data.")
else:
    print("Feedback: error is within 20% of last season - the model can stay in use.")