import numpy as np
import pandas as pd

# ---------------------------------------------------------
# 1) LOAD DATA
# ---------------------------------------------------------
df = pd.read_csv("training_data.csv", header=None,
                 names=["ts","temp","hum","leaf","wind"])

#sorts values by ts and resets index
df = df.sort_values("ts").reset_index(drop=True)

# Extract 4 features
data = df[["temp","hum","leaf","wind"]].values.astype(np.float32)

# ---------------------------------------------------------
# 2) BUILD TRAINING SET
# Input:  24 × 4 = 96 values
# Output: next 24 × 4 = 96 values
# ---------------------------------------------------------
WINDOW = 24
OUT = 24

X_list = []
Y_list = []

for i in range(len(data) - WINDOW - OUT):
    X_list.append(data[i:i+WINDOW].flatten())         # shape (96,)
    Y_list.append(data[i+WINDOW:i+WINDOW+OUT].flatten())  # shape (96,)

X = np.array(X_list)   # shape (N,96) input X
Y = np.array(Y_list)   # shape (N,96) output Y

print("Training samples:", X.shape[0])


# ---------------------------------------------------------
# 3) SOLVE LEAST SQUARES
# Y = XW + b
# Solve for W and b:
#    W = (Xᵀ X)^(-1) Xᵀ Y
#    b = mean(Y - XW)
# ---------------------------------------------------------

# Add bias trick: Xb = [X | 1]
Xb = np.hstack([X, np.ones((X.shape[0],1))])  # shape (N,97)

# Solve using np.linalg.lstsq
"""
residuals — sum of squared errors (Y-pred)^2, rank=97, s=singular values of Xb Large spread between max and min -> ill-conditioned

Good singular values --> stable solution set by rcond
Interpretation
Large spread the matrix is ill-conditioned
Condition number:
    κ= spread/MSE > >10 spread>10,000 ill conditioned
        This means:
        Some input features are highly correlated
        The least squares solution is sensitive to noise
        But since the MSE is reasonable, it still learned a usable mapping
        This is normal when using sliding windows from smooth seasonal data — the inputs are not independent.

"""
W_full, residuals, rank, s = np.linalg.lstsq(Xb, Y, rcond=None)

# W_full shape: (97,96)
W = W_full[:-1, :]     # (96 × 96)
b = W_full[-1, :]      # (96,)

print("W shape:", W.shape)
print("b shape:", b.shape)
print("Max singular value:", np.max(s))
print("Min singular value:", np.min(s))
print("Spread:", np.max(s) - np.min(s))
if residuals.size ==1:
    # residuals[0] = sum of squared errors
    mse = residuals[0] / X.shape[0]
else:
    # compute all manually if residuals empty
    Y_pred = Xb @ W_full
    mse = np.mean((Y - Y_pred)**2)

print("Mean of residuals (MSE):", mse)
# ---------------------------------------------------------
# 6) CONDITION NUMBER (κ) AFTER REGULARIZATION
# ---------------------------------------------------------
s = np.linalg.svd(Xb, compute_uv=False)
cond_number = np.max(s) / np.min(s)
print("Condition number (κ):", cond_number)
# ---------------------------------------------------------
# 4) SAVE MODEL
# ---------------------------------------------------------
np.savez("linear_model.npz", W=W, b=b)
print("Saved linear_model.npz")

# ---------------------------------------------------------
# 5) EXPORT AS C HEADER FILE (for ESP8266)
# ---------------------------------------------------------
with open("lin_model.h","w") as f:
    f.write("// Linear regression model for ESP8266\n")
    f.write("static const int INPUT_SIZE  = 96;\n")
    f.write("static const int OUTPUT_SIZE = 96;\n\n")

    # Write W
    f.write("static const float W_lin[96][96] = {\n")
    for i in range(96):
        row = ",".join([f"{v:.6f}" for v in W[i]])
        f.write("  {" + row + "},\n")
    f.write("};\n\n")

    # Write b
    f.write("static const float b_lin[96] = {")
    f.write(",".join([f"{v:.6f}" for v in b]))
    f.write("};\n")

print("Saved lin_model.h")