<jupyter_start><jupyter_text>5. Estimation – Monte Carlo Augmented Data    About:    This notebook refers to the studies presented in Chapter 5.5 of the Ph.D. thesis [3].    We can not guarantee completeness or correctness of the code.    If you find bugs or if you have suggestions on how to improve the code, we encourage you to post your ideas as GitHub issue. Imports<jupyter_code>%load_ext autoreload
%autoreload 2
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import train_test_split
import sklearn.metrics as me
import numpy as np
import matplotlib.pyplot as plt
from tqdm.notebook import tqdm
import susi
import pandas as pd
import pickle

import utils
def fit_model(model, data):
    X_train, X_test, y_train, y_test = data
    
    # fit and predict
    model.fit(X_train, y_train)
    y_pred = model.predict(X_test)

    # evaluate
    r2 = me.r2_score(y_test, y_pred)
    rmse = np.sqrt(me.mean_squared_error(y_test, y_pred))
    mae = me.mean_absolute_error(y_test, y_pred)
    
    return r2, rmse, mae<jupyter_output><empty_output><jupyter_text>Regression over number of datapoints Generate estimations<jupyter_code># CHANGE
area = ["1", "3"]
# postfix = area[0][0]
postfix = "1+3"

test_sizes = np.round(np.arange(0.15, 0.35, 0.05), 2)
mc_args = {"std": 2, "n_new": 10, "area": area, "max_sm": 40}

# load data for checks only
_, y = utils.loadCSVData(area, max_sm=mc_args["max_sm"])
print("Area {0} with {1} datapoints and soil moisture of {2:.2f} ± {3:.2f} %.".format(
    area, len(y), np.mean(y), np.std(y)))

# for plots only
meas_error = 4.
meas_error_std = 2.2/2

# plot histo
fig, ax = plt.subplots(1, 1)
n, _, _ = plt.hist(y)
ax.set_ylim(0, max(n)*1.5)

ax.set_xlabel("Soil moisture in %")
ax.set_ylabel("Number of datapoints")
results = {
    "model": [],
    "r2": [],
    "rmse": [],
    "mae": [],
    "test_size": [],
    "random_state": []}

for i, test_size in enumerate(tqdm(test_sizes)):

    for random_state in range(40):
        data = utils.getMCSoilMoistureData(
            test_size=test_size, random_state=random_state, **mc_args)

        # --- SOM
        model = susi.SOMRegressor(n_rows=10, n_columns=10,
                                  n_iter_unsupervised=2000, n_iter_supervised=1000)
        r2, rmse, mae = fit_model(model, data)
        results["r2"].append(r2)
        results["rmse"].append(rmse)
        results["mae"].append(mae)
        results["test_size"].append(test_size)
        results["random_state"].append(random_state)
        results["model"].append("SOM")

        # --- RF
        model = RandomForestRegressor(n_estimators=100, n_jobs=-1)
        r2, rmse, mae = fit_model(model, data)
        results["r2"].append(r2)
        results["rmse"].append(rmse)
        results["mae"].append(mae)
        results["test_size"].append(test_size)
        results["random_state"].append(random_state)
        results["model"].append("RF")
        
results_df = pd.DataFrame(results)<jupyter_output><empty_output><jupyter_text>results_df.to_csv("results/mc_som_rf.csv")print(results_df.shape)results_df.head()<jupyter_code>## Plots<jupyter_output><empty_output><jupyter_text>color_list = ["tab:blue", "tab:orange"] fontsize = 18fig, ax = plt.subplots(1, 1, figsize=(6, 5))a_min = -0.5 * 100a_max = 1. * 100y_max = 60bins = np.arange(a_min, a_max, 0.1*100)for i, m in enumerate(np.unique(results_df.model)):    r2_list = results_df[results_df["model"]==m]["r2"].values    ax.hist(np.clip(r2_list*100, a_min=a_min, a_max=a_max), bins=bins, label=m, alpha=0.5,            color=color_list[i])    print("{0:10}  |\tOf {1} are {2} ({3:.1f} %) above 0. Median = {4:.1f} %. Mean = {5:.1f} %.".format(        m, len(np.ravel(r2_list)), sum(np.ravel(r2_list)>0.),        sum(np.ravel(r2_list)>0.)/len(np.ravel(r2_list))*100,        np.median(r2_list)*100, np.mean(r2_list)*100))    ax.axvline(np.median(r2_list)*100, label=m+" median",               color=color_list[i], linestyle="dashed")ax.set_xlim(a_min, a_max)ax.set_ylim(0, y_max)ax.set_xlabel("$R^2$ in %", fontsize=fontsize)ax.set_ylabel("Number of experiments", fontsize=fontsize, labelpad=10)for tick in ax.xaxis.get_major_ticks():    tick.label.set_fontsize(fontsize)for tick in ax.yaxis.get_major_ticks():    tick.label.set_fontsize(fontsize)    plt.legend(fontsize=fontsize*0.8)plt.tight_layout()plt.savefig("plots/mc_area"+postfix+"_hist_r2.pdf", bbox_inches="tight") fontsize = 18fig, ax = plt.subplots(1, 1, figsize=(6, 5))a_min = -0.5y_max = 120bins = np.arange(1, 12, 1)for i, m in enumerate(np.unique(results_df.model)):    mae_list = results_df[results_df["model"]==m]["mae"].values    ax.hist(mae_list, label=m, alpha=0.5, bins=bins, color=color_list[i])    print("{0:10}  |\t Median = {1:.1f} %. Mean = {2:.1f} %.".format(        m, np.median(mae_list), np.mean(mae_list)))    ax.axvline(np.median(mae_list), label=m+" median",               color=color_list[i], linestyle="dashed")ax.set_xlim(bins[0], bins[-1])ax.set_ylim(0, y_max)ax.set_xlabel("MAE in % soil moisture", fontsize=fontsize)ax.set_ylabel("Number of experiments", fontsize=fontsize, labelpad=10)for tick in ax.xaxis.get_major_ticks():    tick.label.set_fontsize(fontsize)for tick in ax.yaxis.get_major_ticks():    tick.label.set_fontsize(fontsize)    plt.legend(fontsize=fontsize*0.8)plt.tight_layout()plt.savefig("plots/mc_area"+postfix+"_hist_mae.pdf", bbox_inches="tight") fontsize = 18fig, ax = plt.subplots(1, 1, figsize=(6, 5))y_max = 120bins = np.arange(1, 13, 1)for i, m in enumerate(np.unique(results_df.model)):    rmse_list = results_df[results_df["model"]==m]["rmse"].values    ax.hist(rmse_list, label=m, alpha=0.5, bins=bins, color=color_list[i])    print("{0:10}  |\t Median = {1:.1f} %. Mean = {2:.1f} %.".format(        m, np.median(rmse_list), np.mean(rmse_list)))    ax.axvline(np.median(rmse_list), label=m+" median",               color=color_list[i], linestyle="dashed")ax.set_xlim(bins[0], bins[-1])ax.set_ylim(0, y_max)ax.set_xlabel("RMSE in % soil moisture", fontsize=fontsize)ax.set_ylabel("Number of experiments", fontsize=fontsize, labelpad=10)for tick in ax.xaxis.get_major_ticks():    tick.label.set_fontsize(fontsize)for tick in ax.yaxis.get_major_ticks():    tick.label.set_fontsize(fontsize)    plt.legend(fontsize=fontsize*0.8)plt.tight_layout()plt.savefig("plots/mc_area"+postfix+"_hist_rmse.pdf", bbox_inches="tight")<jupyter_code>## MC data histogram plot<jupyter_output><empty_output><jupyter_text>X, y = utils.loadCSVData(areas=area, max_sm=mc_args["max_sm"])X_new, y_new = utils.generateMCData(X=X, y=y, std=mc_args["std"], n_new=mc_args["n_new"], verbose=0)normed = Falsebins = np.arange(int(min(y_new)), int(max(y_new))+2, 2)fontsize = 15fig, ax = plt.subplots(1, 1, figsize=(7, 4))n, _, _ = ax.hist(y_new, density=normed, alpha=1.0, label="Monte Carlo", bins=bins)ax.hist(y, density=normed, alpha=1.0, label="Original", bins=bins, histtype="step", linewidth=3) ax.set_title("Area "+area[0][0], fontsize=fontsize)ax.legend(fontsize=fontsize*1.0, frameon=False)ax.set_ylim(0, max(n)*1.5)ax.set_xlabel("Soil moisture in %", fontsize=fontsize)ax.set_ylabel("Number of datapoints", fontsize=fontsize)for tick in ax.xaxis.get_major_ticks():    tick.label.set_fontsize(fontsize)for tick in ax.yaxis.get_major_ticks():    tick.label.set_fontsize(fontsize)    plt.tight_layout()plt.savefig("plots/mc_area"+postfix+"_hist.pdf", bbox_inches="tight")<jupyter_code>## Scatter plot<jupyter_output><empty_output><jupyter_text>results_df.iloc[309]   [(results_df["test_size"]==0.3) & (results_df["model"]=="SOM")]results_df["r2"].tail(40)<jupyter_code># CHANGE index for plot
i = 312

model_single = susi.SOMRegressor(n_rows=10, n_columns=10,
                          n_iter_unsupervised=2000, n_iter_supervised=1000)
X_train, X_test, y_train, y_test = utils.getMCSoilMoistureData(
    test_size=results_df.iloc[i]["test_size"],
    random_state=results_df.iloc[i]["random_state"], **mc_args)

# fit and predict
model_single.fit(X_train, y_train)
y_pred_single = model_single.predict(X_test)
print("R2 = {0:.1f} %".format(100*me.r2_score(y_test, y_pred_single)))

# plot
fontsize = 15
fig, ax = plt.subplots(1, 1, figsize=(6, 6))
factor = mc_args["n_new"]+1
ax.scatter(y_test, y_pred_single, alpha=0.4, label="Monte Carlo")
ax.scatter(y_test[::factor], y_pred_single[::factor], marker="x", s=100, label="Measurements")
ax.set_xlabel("Soil moisture (ground truth) in %", fontsize=fontsize, labelpad=10)
ax.set_ylabel("Soil moisture (estimated) in %", fontsize=fontsize, labelpad=10)
for tick in ax.xaxis.get_major_ticks():
    tick.label.set_fontsize(fontsize)
for tick in ax.yaxis.get_major_ticks():
    tick.label.set_fontsize(fontsize)

# min_plot = min(min(y_test), min(y_pred_single))
# max_plot = max(max(y_test), max(y_pred_single))
min_plot = 1
max_plot = 34
ax.set_xlim(min_plot-1, max_plot+1)
ax.set_ylim(min_plot-1, max_plot+1)
ax.plot([min_plot-2, max_plot+2], [min_plot-2, max_plot+2], color="grey", linewidth=1)
ax.legend(fontsize=fontsize)
plt.savefig("plots/mc_area"+postfix+"_scatter.pdf", bbox_inches="tight")<jupyter_output><empty_output><jupyter_text>SOM estimation map<jupyter_code>def plot_estimation_map_regression(estimation_map, sm_range=None,
                                   title="", fontsize=20):
    plt.figure(figsize=(7,5))
    if sm_range:
        plt.imshow(estimation_map, cmap="viridis_r",
                   vmin=sm_range[0], vmax=sm_range[1])
    else:
        plt.imshow(estimation_map, cmap="viridis_r")
    plt.xlabel("SOM columns", fontsize=fontsize)
    plt.ylabel("SOM rows", fontsize=fontsize)
    plt.xticks(fontsize=fontsize)
    plt.yticks(fontsize=fontsize)
    cbar = plt.colorbar()
    cbar.ax.tick_params(labelsize=fontsize)
    cbar.ax.set_ylabel('Soil moisture in %', fontsize=fontsize, labelpad=10)
    for label in cbar.ax.xaxis.get_ticklabels()[::2]:
        label.set_visible(False)
    plt.grid(b=False)
    plt.tight_layout()
    plt.savefig("plots/mc_area"+postfix+"_estimationmap"+title+".pdf",
                bbox_inches="tight")
esti_map = model_single.get_estimation_map()
plot_estimation_map_regression(np.squeeze(esti_map), sm_range=(9., 26.))<jupyter_output><empty_output><jupyter_text>Plot Estimation Map Area 1<jupyter_code>hyp_data_map = utils.predictSoilmoistureMap(
    area=area[0], model=model_single, dim_red_mode=None,
    # sm_range=(12., 29.),
    postfix="mc", verbose=1)
pickle.dump(hyp_data_map, open("estimations/hyp_data_map_mc_area1.p", "wb"))
hyp_data_map = pickle.load(open("estimations/hyp_data_map_mc_area1.p", "rb"))
print("Soil moisture between {0:.1f} and {1:.1f} %".format(
    np.min(hyp_data_map[hyp_data_map !=0]), np.max(hyp_data_map)))
utils.plotSoilMoistureMap(hyp_data_map, area=area[0][0], sm_range=(9., 26.), postfix="mc")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16219467061435974, 'fraction_numerical': 0.028589933382679497, 'mean_word_length': 5.099887133182844, 'pattern_counts': {'":': 10, '<': 33, '<?xml version=': 0, '>': 35, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '506419', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 4465, 'n_tokens_neox': 4219, 'n_words': 680, 'path': 'py/5_Estimation_MonteCarloAugmentedData.ipynb', 'repo_name': 'felixriese/alpaca-processing', 'size': '10808'}