Pretrattamento iperspettrale

 ...ovvero tutto cio' che avrei voluto sapere in dottorato ma nessuno mi ha mai detto.

La domanda e': quanto conta il preprocessamento delle firme spettrali in funzione della successiva regressione? Tanto...di seguito un esempio

Partiamo da questo database di firme spettrali (1350-2500 nm) ed associate caratteristiche chimico fisiche dei suoli 
 

Sanderman, J., Smith, C., Safanelli, J. L., Mitu, S. M., Ge, Y., Murad, O., & Shepherd, K. (2023). Near-infrared (NIR) soil spectral library using the NeoSpectra Handheld NIR Analyzer by Si-Ware (1.0) [Data set]. Zenodo. https://doi.org/10.5281/zenodo.7586622 

 

Dal file Neospectra_WoodwellKSSL_avg_soil+site+NIR.csv ho estrapolata solo la colonna EC ed i valori di riflettanza NIR e proviamo senza pretrattamento ad effettuare un regressione Random Forest (2969 spettri di dataset)

==============================
Model Performance (Target: ec)
R² Score: 0.2661
RMSE:     4.0410
==============================
 

import pandas as pd

import numpy as np

import os

from sklearn.ensemble import RandomForestRegressor

from sklearn.model_selection import train_test_split

from sklearn.metrics import mean_squared_error, r2_score

# --- Settings ---

LABEL_FILE = 'label.csv'

SPECTRA_FOLDER = 'spectra'

ID_COL = 'id' # The column in label.csv

TARGET_COL = 'ec' # The lab value column

def load_data():

df_labels = pd.read_csv(LABEL_FILE)

df_labels = df_labels.dropna(subset=[ID_COL, TARGET_COL])

df_labels[ID_COL] = df_labels[ID_COL].astype(int).astype(str)

X_list = []

y_list = []

print(f"Searching for spectra in: {os.path.abspath(SPECTRA_FOLDER)}")

for _, row in df_labels.iterrows():

sample_id = row[ID_COL]

# Construct filename: e.g., '3748.csv'

file_path = os.path.join(SPECTRA_FOLDER, f"{sample_id}.csv")

if os.path.exists(file_path):

try:

# 2. Load spectrum and ensure wavelength alignment

spec_df = pd.read_csv(file_path)

spec_df = spec_df.sort_values('wavelength', ascending=False)

X_list.append(spec_df['reflectance'].values)

y_list.append(row[TARGET_COL])

except Exception as e:

print(f"Error reading {sample_id}.csv: {e}")

else:

print(f"File not found: {file_path}")

return np.array(X_list), np.array(y_list)

# --- Execution ---

X, y = load_data()

if X.size > 0:

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

rf_model = RandomForestRegressor(n_estimators=500, max_features='sqrt', n_jobs=-1, random_state=42)

print(f"\nTraining on {X_train.shape[0]} samples...")

rf_model.fit(X_train, y_train)

y_pred = rf_model.predict(X_test)

rmse = np.sqrt(mean_squared_error(y_test, y_pred))

r2 = r2_score(y_test, y_pred)

print("\n" + "="*30)

print(f"Model Performance (Target: {TARGET_COL})")

print(f"R² Score: {r2:.4f}")

print(f"RMSE: {rmse:.4f}")

print("="*30)

else:

print("\nNo data loaded. Check if the 'spectra' folder is in the same directory as this script.")

adesso stesso dataset ma con preprocessing cambiamento da riflettanza in assorbanza, SNV + Savitky Gorlay derivata prima ...a parita' di input il modello e' 

Training RF with 1st Derivative + SNV...

==============================
Derivative Results
R² Score: 0.6482
RMSE:     2.7978
==============================

import pandas as pd

import numpy as np

import os

from scipy.signal import savgol_filter

from sklearn.ensemble import RandomForestRegressor

from sklearn.model_selection import train_test_split

from sklearn.metrics import mean_squared_error, r2_score

# --- Settings ---

LABEL_FILE = 'label.csv'

SPECTRA_FOLDER = 'spectra'

ID_COL = 'id'

TARGET_COL = 'ec'

def load_and_preprocess_advanced():

df_labels = pd.read_csv(LABEL_FILE).dropna(subset=[ID_COL, TARGET_COL])

df_labels[ID_COL] = df_labels[ID_COL].astype(int).astype(str)

X_list, y_list = [], []

for _, row in df_labels.iterrows():

file_path = os.path.join(SPECTRA_FOLDER, f"{row[ID_COL]}.csv")

if os.path.exists(file_path):

spec_df = pd.read_csv(file_path).sort_values('wavelength', ascending=False)

# 1. Absorbance

ref = np.clip(spec_df['reflectance'].values, 0.01, 100)

absorbance = np.log10(100.0 / ref)

X_list.append(absorbance)

y_list.append(row[TARGET_COL])

X = np.array(X_list)

# 2. Savitzky-Golay 1st Derivative

# window_length: must be odd. 11-15 is usually good for NIR.

# polyorder: usually 2.

X_deriv = savgol_filter(X, window_length=11, polyorder=2, deriv=1)

# 3. SNV (applied to the derivative)

X_final = (X_deriv - np.mean(X_deriv, axis=1, keepdims=True)) / np.std(X_deriv, axis=1, keepdims=True)

return X_final, np.array(y_list)

# --- Execution ---

X, y = load_and_preprocess_advanced()

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

# Use the "Best" RF setup from your 0.42 run

rf = RandomForestRegressor(n_estimators=500, max_features='sqrt', n_jobs=-1, random_state=42)

print("Training RF with 1st Derivative + SNV...")

rf.fit(X_train, y_train)

y_pred = rf.predict(X_test)

print("\n" + "="*30)

print(f"Derivative Results")

print(f"R² Score: {r2_score(y_test, y_pred):.4f}")

print(f"RMSE: {np.sqrt(mean_squared_error(y_test, y_pred)):.4f}")

print("="*30)

Facciamo ancora un ultimo passo applicando la PCA (analisi delle componenti principali) al preprocessing

Training RF su componenti PCA...

==============================
PCA + Random Forest Results
R² Score: 0.7491
RMSE:     2.3630
==============================
 

import pandas as pd

import numpy as np

import os

from scipy.signal import savgol_filter

from sklearn.ensemble import RandomForestRegressor

from sklearn.decomposition import PCA

from sklearn.preprocessing import StandardScaler

from sklearn.model_selection import train_test_split

from sklearn.metrics import mean_squared_error, r2_score

# --- Settings ---

LABEL_FILE = 'label.csv'

SPECTRA_FOLDER = 'spectra'

ID_COL = 'id'

TARGET_COL = 'ec'

def load_and_preprocess_advanced():

df_labels = pd.read_csv(LABEL_FILE).dropna(subset=[ID_COL, TARGET_COL])

df_labels[ID_COL] = df_labels[ID_COL].astype(int).astype(str)

X_list, y_list = [], []

for _, row in df_labels.iterrows():

file_path = os.path.join(SPECTRA_FOLDER, f"{row[ID_COL]}.csv")

if os.path.exists(file_path):

spec_df = pd.read_csv(file_path).sort_values('wavelength', ascending=False)

ref = np.clip(spec_df['reflectance'].values, 0.01, 100)

absorbance = np.log10(100.0 / ref)

X_list.append(absorbance)

y_list.append(row[TARGET_COL])

X = np.array(X_list)

# Pre-processing classico (Derivata + SNV)

X_deriv = savgol_filter(X, window_length=11, polyorder=2, deriv=1)

X_snv = (X_deriv - np.mean(X_deriv, axis=1, keepdims=True)) / np.std(X_deriv, axis=1, keepdims=True)

return X_snv, np.array(y_list)

# --- 1. Caricamento ---

X, y = load_and_preprocess_advanced()

X_train_raw, X_test_raw, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# --- 2. PCA (Principal Component Analysis) ---

# n_components=0.99 significa che manteniamo abbastanza componenti per spiegare il 99% della varianza

pca = PCA(n_components=0.99, random_state=42)

# È fondamentale fittare la PCA solo sul train e trasformare il test

X_train_pca = pca.fit_transform(X_train_raw)

X_test_pca = pca.transform(X_test_raw)

print(f"Numero di Componenti Principali selezionate: {pca.n_components_}")

print(f"Varianza totale spiegata: {np.sum(pca.explained_variance_ratio_):.4f}")

# --- 3. Random Forest sulle Componenti Principali ---

# Con meno feature (PCA), possiamo spesso usare parametri più aggressivi

rf = RandomForestRegressor(n_estimators=500, max_features=None, n_jobs=-1, random_state=42)

print("Training RF su componenti PCA...")

rf.fit(X_train_pca, y_train)

# --- 4. Risultati ---

y_pred = rf.predict(X_test_pca)

print("\n" + "="*30)

print(f"PCA + Random Forest Results")

print(f"R² Score: {r2_score(y_test, y_pred):.4f}")

print(f"RMSE: {np.sqrt(mean_squared_error(y_test, y_pred)):.4f}")

print("="*30)

 

in sintesi siamo partiti da R2= 0.27 a R2=0.75 solo di preprocessing 

Messo a posto il discorso preprocessing vediamo cosa succede a cambiare modello. Applicando PLSR si ha

Best number of components: 25 (CV R²: 0.2882)

==============================
Final PLS Model Results (n=25)
R² Score: 0.2732
RMSE:     4.0215
==============================
 

il drastico peggioramento dovrebbe essere dovuto al fatto della non linerita' della risposta spettrale nei confronti del parametro analizzato (PLSR lavora su una ipotesi di base di linearita')

facciamo un'ultimo tentativo con Boost Gradient

Riduzione dimensionalità: da 257 a 41 componenti PCA.
Training GBR su componenti PCA...

==============================
GBR + PCA Test Results:
R² Score: 0.7200
RMSE:     2.4963
==============================
 

import pandas as pd

import numpy as np

import os

from scipy.signal import savgol_filter

from sklearn.ensemble import GradientBoostingRegressor

from sklearn.model_selection import train_test_split

from sklearn.metrics import mean_squared_error, r2_score

# --- Settings ---

TRAIN_LABELS = 'label.csv'

TRAIN_DIR = 'spectra'

VAL_LABELS = 'valida.csv'

VAL_DIR = 'new_spectra'

ID_COL = 'id'

TARGET_COL = 'ec'

def preprocess_spectra(folder, id_list):

"""Absorbance -> SG Derivative -> SNV"""

features = []

processed_ids = []

for sid in id_list:

# Convert ID to string to avoid .0 issues

clean_id = str(int(float(sid)))

path = os.path.join(folder, f"{clean_id}.csv")

if os.path.exists(path):

spec_df = pd.read_csv(path).sort_values('wavelength', ascending=False)

# 1. Absorbance

ref = np.clip(spec_df['reflectance'].values, 0.01, 100)

absorbance = np.log10(100.0 / ref)

features.append(absorbance)

processed_ids.append(sid)

X = np.array(features)

# 2. 1st Derivative (Smooths and highlights peaks)

X = savgol_filter(X, window_length=11, polyorder=2, deriv=1)

# 3. SNV (Standardizes the scale)

X = (X - np.mean(X, axis=1, keepdims=True)) / np.std(X, axis=1, keepdims=True)

return X, processed_ids

# --- 1. Data Preparation ---

print("Loading and preprocessing training data...")

df_labels = pd.read_csv(TRAIN_LABELS).dropna(subset=[ID_COL, TARGET_COL])

X_train_full, train_ids = preprocess_spectra(TRAIN_DIR, df_labels[ID_COL].values)

# Align y with the IDs that actually had files

y_train_full = df_labels.set_index(ID_COL).loc[train_ids][TARGET_COL].values

# Split for internal testing

X_train, X_test, y_train, y_test = train_test_split(X_train_full, y_train_full, test_size=0.2, random_state=42)

# --- 2. Gradient Boosting Model ---

# learning_rate: smaller is usually better but requires more n_estimators

# max_depth: usually 3-5 for boosting to prevent overfitting

gbr = GradientBoostingRegressor(

n_estimators=1000,

learning_rate=0.05,

max_depth=4,

subsample=0.8, # Uses 80% of data for each tree to add randomness

random_state=42,

verbose=1

)

print("\nTraining Gradient Boosting Regressor...")

gbr.fit(X_train, y_train)

# Internal Evaluation

y_pred_internal = gbr.predict(X_test)

print("\n" + "="*30)

print(f"GBR Internal Test Results:")

print(f"R² Score: {r2_score(y_test, y_pred_internal):.4f}")

print(f"RMSE: {np.sqrt(mean_squared_error(y_test, y_pred_internal)):.4f}")

print("="*30)

# --- 3. Inference on Validation Set ---

print(f"\nProcessing inference for {VAL_DIR}...")

df_valida = pd.read_csv(VAL_LABELS)

X_val, val_ids = preprocess_spectra(VAL_DIR, df_valida[ID_COL].values)

if X_val.size > 0:

val_predictions = gbr.predict(X_val)

# Save Results

results = pd.DataFrame({ID_COL: val_ids, 'predicted_ec': val_predictions})

results.to_csv('valida_predictions_gbr.csv', index=False)

print("Inference results saved to 'valida_predictions_gbr.csv'")

import matplotlib.pyplot as plt

importances = gbr.feature_importances_

plt.plot(importances)

plt.title("Wavelength Importance for EC")

plt.show()

 

 

 

 

 

Lupi nei dintorni di Borselli

 29 marzo : Oratorio della Madonna dei Fossi

  

Hyperview Competition (2)

Riguardo al post precedente, la parte di CNN e' piuttosto pesante dal punto di vista di calcolo

Ho provato a vedere che risultati si ottengono passando direttamente alla seconda fase 

 

"""

HyperSoilNet Pipeline

=====================

Based on: "A Hybrid Framework for Soil Property Estimation from Hyperspectral Imaging"

https://www.mdpi.com/2072-4292/17/15/2568

Implements:

- Spectral feature engineering (derivatives, DWT, SVD, FFT)

- ML ensemble: Random Forest + XGBoost + KNN

- Property-specific weighted ensemble (Bayesian-optimised weights)

- 5-fold cross-validation with per-property R², RMSE, custom challenge score

Input layout

------------

labels.csv : sample_id, P, K, Mg, pH (header on first row)

<spectra_dir>/1.csv : one reflectance value per row, one column (no header)

<spectra_dir>/6.csv : same for sample 6, etc.

Usage

-----

# Train + cross-validate on labelled data

python hypersoilnet.py --labels label.csv --spectra_dir spectra/

# Train on all labelled data, then predict an unlabelled set

python hypersoilnet.py --labels label.csv --spectra_dir spectra/ \

--predict_dir unlabelled/ --output predictions.csv

"""

import argparse

import os

import warnings

from pathlib import Path

import numpy as np

import pandas as pd

from scipy.fft import fft

from scipy.signal import savgol_filter

try:

import pywt

HAS_PYWT = True

except ImportError:

HAS_PYWT = False

warnings.warn("PyWavelets (pywt) not found – wavelet features will be skipped. "

"Install with: pip install PyWavelets")

from sklearn.ensemble import RandomForestRegressor

from sklearn.neighbors import KNeighborsRegressor

from sklearn.preprocessing import StandardScaler

from sklearn.model_selection import KFold

from sklearn.metrics import mean_squared_error, r2_score

try:

from xgboost import XGBRegressor

HAS_XGB = True

except ImportError:

HAS_XGB = False

warnings.warn("XGBoost not found – ensemble will use RF + KNN only. "

"Install with: pip install xgboost")

try:

from scipy.optimize import minimize

HAS_SCIPY_OPT = True

except ImportError:

HAS_SCIPY_OPT = False

# ---------------------------------------------------------------------------

# Target column names (must match label CSV header)

# ---------------------------------------------------------------------------

TARGETS = ["P", "K", "Mg", "pH"]

# ---------------------------------------------------------------------------

# Property-specific ensemble weights from the paper (Table in Section 3.5)

# Order: [RF, XGB, KNN]

# ---------------------------------------------------------------------------

DEFAULT_WEIGHTS = {

"K": [0.45, 0.40, 0.15],

"P": [0.38, 0.45, 0.17], # paper calls it P2O5; mapped to "P" here

"Mg": [0.42, 0.38, 0.20],

"pH": [0.35, 0.50, 0.15],

}

# ===========================================================================

# 1. DATA LOADING

# ===========================================================================

def load_spectrum(filepath: str) -> np.ndarray:

"""Load a single-column CSV (no header) of reflectance values."""

df = pd.read_csv(filepath, header=None)

return df.iloc[:, 0].values.astype(float)

def load_dataset(labels_csv: str, spectra_dir: str):

"""

Returns

-------

spectra : np.ndarray shape (n_samples, n_bands)

labels : pd.DataFrame columns = TARGETS

sample_ids: list of sample ids

"""

labels_df = pd.read_csv(labels_csv)

# Strip whitespace from column names

labels_df.columns = labels_df.columns.str.strip()

spectra_list, valid_ids, valid_rows = [], [], []

for _, row in labels_df.iterrows():

sid = int(row["sample_id"])

fpath = os.path.join(spectra_dir, f"{sid}.csv")

if not os.path.exists(fpath):

warnings.warn(f"Spectrum file not found, skipping: {fpath}")

continue

spectra_list.append(load_spectrum(fpath))

valid_ids.append(sid)

valid_rows.append(row)

spectra = np.vstack(spectra_list)

labels = pd.DataFrame(valid_rows)[TARGETS].reset_index(drop=True)

return spectra, labels, valid_ids

def load_unlabelled(spectra_dir: str):

"""Load all CSVs in a directory that have no matching label."""

paths = sorted(Path(spectra_dir).glob("*.csv"))

spectra_list, ids = [], []

for p in paths:

spectra_list.append(load_spectrum(str(p)))

ids.append(p.stem)

return np.vstack(spectra_list), ids

# ===========================================================================

# 2. FEATURE ENGINEERING (Section 3.4 of the paper)

# ===========================================================================

def spectral_derivatives(spectrum: np.ndarray, orders=(1, 2, 3)) -> np.ndarray:

"""1st, 2nd, 3rd order Savitzky-Golay derivatives."""

feats = []

n = len(spectrum)

for order in orders:

polyorder = order + 1

# window must be odd, > polyorder, and <= len(spectrum)

window = min(11, n)

if window % 2 == 0:

window -= 1

window = max(window, polyorder + 1)

if (polyorder + 1) % 2 == 0:

window = polyorder + 2 # next odd number above polyorder

if window > n:

# spectrum too short for this derivative order – return zeros

feats.append(np.zeros(n))

continue

d = savgol_filter(spectrum, window_length=window, polyorder=polyorder,

deriv=order)

feats.append(d)

return np.concatenate(feats)

def wavelet_features(spectrum: np.ndarray, wavelet="dmey", level=4) -> np.ndarray:

"""DWT approximation + detail coefficient statistics (paper eq. 5)."""

if not HAS_PYWT:

return np.array([])

# Clamp level to what the signal length allows

max_level = pywt.dwt_max_level(len(spectrum), wavelet)

level = min(level, max_level)

coeffs = pywt.wavedec(spectrum, wavelet, level=level)

feats = []

for c in coeffs:

feats.extend([c.mean(), c.std(), np.abs(c).max()])

return np.array(feats)

def svd_features(spectrum: np.ndarray, n_singular=5) -> np.ndarray:

"""

Treat spectrum as a 1-D signal; reshape into a 2-D matrix and compute SVD.

Top singular values + their ratios (paper eq. 6).

"""

n = len(spectrum)

# Build a Hankel-like matrix

n_rows = max(2, n // 10)

n_cols = n - n_rows + 1

mat = np.array([spectrum[i:i + n_cols] for i in range(n_rows)])

_, s, _ = np.linalg.svd(mat, full_matrices=False)

s = s[:n_singular] if len(s) >= n_singular else np.pad(s, (0, n_singular - len(s)))

ratios = s[:-1] / (s[1:] + 1e-10)

return np.concatenate([s, ratios])

def fft_features(spectrum: np.ndarray, n_components=20) -> np.ndarray:

"""Real + imaginary FFT components (paper eq. 7)."""

f = fft(spectrum)

f = f[:n_components]

return np.concatenate([f.real, f.imag])

def statistical_features(spectrum: np.ndarray) -> np.ndarray:

"""Basic statistics of the raw spectrum."""

return np.array([

spectrum.mean(),

spectrum.std(),

np.percentile(spectrum, 25),

np.percentile(spectrum, 75),

spectrum.max() - spectrum.min(),

np.argmax(spectrum),

np.argmin(spectrum),

])

def extract_features(spectrum: np.ndarray) -> np.ndarray:

"""Combine all feature groups into one vector."""

parts = [

spectrum, # raw reflectance

statistical_features(spectrum),

spectral_derivatives(spectrum),

wavelet_features(spectrum),

svd_features(spectrum),

fft_features(spectrum),

]

return np.concatenate([p for p in parts if len(p) > 0])

def build_feature_matrix(spectra: np.ndarray) -> np.ndarray:

rows = [extract_features(s) for s in spectra]

return np.array(rows)

# ===========================================================================

# 3. MODEL BUILDING (Section 3.5)

# ===========================================================================

def build_models():

models = {

"RF": RandomForestRegressor(

n_estimators=150,

max_depth=30,

min_samples_leaf=5,

criterion="squared_error",

bootstrap=True,

n_jobs=-1,

random_state=42,

),

"KNN": KNeighborsRegressor(

n_neighbors=7,

weights="distance",

metric="euclidean",

),

}

if HAS_XGB:

models["XGB"] = XGBRegressor(

n_estimators=100,

learning_rate=0.1,

max_depth=5,

reg_alpha=0.01,

reg_lambda=1.0,

early_stopping_rounds=15,

eval_metric="rmse",

random_state=42,

verbosity=0,

)

return models

# ===========================================================================

# 4. WEIGHTED ENSEMBLE PREDICTION (paper eq. 8)

# ===========================================================================

def get_weights(target: str, models: dict) -> list:

"""

Return per-model weights for a target property.

If XGB is unavailable, re-normalise between RF and KNN.

"""

paper_w = DEFAULT_WEIGHTS.get(target, [0.40, 0.40, 0.20])

rf_w, xgb_w, knn_w = paper_w

if not HAS_XGB:

total = rf_w + knn_w

return [rf_w / total, knn_w / total]

return [rf_w, xgb_w, knn_w]

def ensemble_predict(fitted_models: dict, X: np.ndarray, target: str) -> np.ndarray:

keys = list(fitted_models.keys())

weights = get_weights(target, fitted_models)

preds = np.stack([fitted_models[k].predict(X) for k in keys], axis=1)

w = np.array(weights[:len(keys)])

w = w / w.sum()

return preds @ w

# ===========================================================================

# 5. TRAINING HELPERS

# ===========================================================================

def fit_ensemble(X_train: np.ndarray, y_train: np.ndarray,

X_val: np.ndarray = None, y_val: np.ndarray = None,

target: str = "P") -> dict:

"""Fit all regressors for one target property."""

scaler = StandardScaler()

X_tr_s = scaler.fit_transform(X_train)

X_val_s = scaler.transform(X_val) if X_val is not None else None

models = build_models()

fitted = {}

for name, mdl in models.items():

if name == "XGB" and X_val_s is not None:

mdl.fit(X_tr_s, y_train,

eval_set=[(X_val_s, y_val)],

verbose=False)

elif name == "XGB":

# No validation set → disable early stopping

mdl.set_params(early_stopping_rounds=None)

mdl.fit(X_tr_s, y_train)

else:

mdl.fit(X_tr_s, y_train)

fitted[name] = mdl

return {"scaler": scaler, "models": fitted}

# ===========================================================================

# 6. CROSS-VALIDATION (Section 3.6 — 5-fold CV)

# ===========================================================================

def cross_validate(X: np.ndarray, labels: pd.DataFrame, n_splits: int = 5):

kf = KFold(n_splits=n_splits, shuffle=True, random_state=42)

results = {t: {"r2": [], "rmse": [], "mse": []} for t in TARGETS}

for fold, (tr_idx, val_idx) in enumerate(kf.split(X), 1):

X_tr, X_val = X[tr_idx], X[val_idx]

print(f"\n Fold {fold}/{n_splits}")

for target in TARGETS:

y_tr = labels[target].values[tr_idx]

y_val = labels[target].values[val_idx]

bundle = fit_ensemble(X_tr, y_tr, X_val, y_val, target)

scaler = bundle["scaler"]

X_tr_s = scaler.transform(X_tr)

X_val_s = scaler.transform(X_val)

preds = ensemble_predict(bundle["models"], X_val_s, target)

mse = mean_squared_error(y_val, preds)

r2 = r2_score(y_val, preds)

rmse = np.sqrt(mse)

results[target]["r2"].append(r2)

results[target]["rmse"].append(rmse)

results[target]["mse"].append(mse)

print(f" {target:4s} R²={r2:.3f} RMSE={rmse:.3f}")

print("\n" + "=" * 55)

print(f"{'Target':<6} {'R² mean':>9} {'R² std':>8} {'RMSE mean':>10}")

print("-" * 55)

baseline_mse = {}

for target in TARGETS:

r2_m = np.mean(results[target]["r2"])

r2_s = np.std(results[target]["r2"])

rmse_m = np.mean(results[target]["rmse"])

print(f"{target:<6} {r2_m:>9.3f} {r2_s:>8.3f} {rmse_m:>10.3f}")

baseline_mse[target] = np.var(labels[target].values) # naive baseline

# Custom challenge score (paper eq. 11) – mean ratio algo_MSE / baseline_MSE

mean_mses = {t: np.mean(results[t]["mse"]) for t in TARGETS}

score = np.mean([mean_mses[t] / baseline_mse[t] for t in TARGETS])

print(f"\n Custom challenge score (lower is better): {score:.4f}")

print("=" * 55)

return results

# ===========================================================================

# 7. FULL TRAINING + PREDICTION

# ===========================================================================

def train_full(X: np.ndarray, labels: pd.DataFrame):

"""Train on the entire labelled dataset."""

# Hold out 20 % for ensemble weight validation (paper Section 3.6)

n = len(X)

val_n = max(1, int(0.2 * n))

rng = np.random.default_rng(42)

val_idx = rng.choice(n, val_n, replace=False)

tr_idx = np.setdiff1d(np.arange(n), val_idx)

bundles = {}

for target in TARGETS:

y = labels[target].values

bundle = fit_ensemble(X[tr_idx], y[tr_idx],

X[val_idx], y[val_idx], target)

bundles[target] = bundle

print(f" Trained ensemble for {target}")

return bundles

def predict(bundles: dict, X_new: np.ndarray) -> pd.DataFrame:

preds = {}

for target in TARGETS:

scaler = bundles[target]["scaler"]

X_s = scaler.transform(X_new)

preds[target] = ensemble_predict(bundles[target]["models"], X_s, target)

return pd.DataFrame(preds)

# ===========================================================================

# 8. MAIN

# ===========================================================================

def parse_args():

p = argparse.ArgumentParser(description="HyperSoilNet – soil property estimation")

p.add_argument("--labels", required=True, help="Path to label CSV")

p.add_argument("--spectra_dir", required=True, help="Dir containing <sample_id>.csv spectra")

p.add_argument("--predict_dir", default=None, help="Dir with unlabelled spectra to predict")

p.add_argument("--output", default="predictions.csv",

help="Output CSV for predictions (default: predictions.csv)")

p.add_argument("--no_cv", action="store_true",

help="Skip cross-validation (faster)")

p.add_argument("--cv_folds", type=int, default=5)

return p.parse_args()

def main():

args = parse_args()

# ---- Load & featurise labelled data ------------------------------------

print("\n[1/4] Loading labelled data …")

spectra, labels, ids = load_dataset(args.labels, args.spectra_dir)

print(f" {len(ids)} samples loaded, {spectra.shape[1]} bands each")

print("\n[2/4] Extracting features …")

X = build_feature_matrix(spectra)

print(f" Feature matrix: {X.shape}")

# ---- Cross-validation --------------------------------------------------

if not args.no_cv:

print(f"\n[3/4] {args.cv_folds}-fold cross-validation …")

cross_validate(X, labels, n_splits=args.cv_folds)

else:

print("\n[3/4] Cross-validation skipped (--no_cv)")

# ---- Train on full dataset ---------------------------------------------

print("\n[4/4] Training on full dataset …")

bundles = train_full(X, labels)

# ---- Predict unlabelled data if requested ------------------------------

if args.predict_dir:

print(f"\n[+] Predicting unlabelled spectra in: {args.predict_dir}")

spectra_new, new_ids = load_unlabelled(args.predict_dir)

X_new = build_feature_matrix(spectra_new)

pred_df = predict(bundles, X_new)

pred_df.insert(0, "sample_id", new_ids)

pred_df.to_csv(args.output, index=False)

print(f" Predictions saved to: {args.output}")

print(pred_df.to_string(index=False))

else:

# Run predictions on the training set itself as a sanity check

print("\n[+] Running predictions on training set (sanity check) …")

all_preds = predict(bundles, X)

all_preds.insert(0, "sample_id", ids)

out_path = "train_predictions.csv"

all_preds.to_csv(out_path, index=False)

print(f" Saved to: {out_path}")

print("\nDone.\n")

if __name__ == "__main__":

main()

Training 
python hypersoilnet.py --labels label.csv --spectra_dir spectra/

Validazione 

python hypersoilnet.py --labels label.csv --spectra_dir spectra/ \
                       --predict_dir new_spectra/ --output predictions.csv

 

Come si vede dai grafici sottostanti i risultati sono in generali peggiori ma non cosi' tanto da essere scartati a priori ....potrebbe essere una valida pre analisi