diff --git a/.gitignore b/.gitignore index 379f9da7..c4c55af7 100644 --- a/.gitignore +++ b/.gitignore @@ -52,6 +52,7 @@ regions/Norway/data/ logs/ venv/ +publish/ *.parquet *.pt diff --git a/.readthedocs.yaml b/.readthedocs.yaml index c8158f0a..590dfb58 100644 --- a/.readthedocs.yaml +++ b/.readthedocs.yaml @@ -16,7 +16,7 @@ build: - pip install poetry - poetry config virtualenvs.create false post_install: - - VIRTUAL_ENV=$READTHEDOCS_VIRTUALENV_PATH poetry install --with docs + - VIRTUAL_ENV=$READTHEDOCS_VIRTUALENV_PATH poetry install --with docs --without publish,viz,advanced # Build documentation in the "docs/" directory with Sphinx diff --git a/docs/conf.py b/docs/conf.py index 0cd4c431..d6dd5589 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -112,6 +112,9 @@ templates_path = ["_templates"] +# Print timing report +duration_show_summary = True + # -- Options for EPUB output epub_show_urls = "footnote" diff --git a/massbalancemachine/__init__.py b/massbalancemachine/__init__.py index c0588a68..85f1b242 100644 --- a/massbalancemachine/__init__.py +++ b/massbalancemachine/__init__.py @@ -2,6 +2,8 @@ sys.path.append(os.path.dirname(os.path.realpath(__file__))) +import importlib + __all__ = [ "dataloader", "models", @@ -24,6 +26,14 @@ import training import plots import metrics -import sampling + +# import sampling # Do not import by default since this is an advanced feature import utils from .config import * # Load config at the top level of the package + + +# Import only if the user asks for it +def __getattr__(name): + if name == "sampling": + return importlib.import_module(".sampling", __name__) + raise AttributeError(f"module {__name__} has no attribute {name}") diff --git a/massbalancemachine/data_preprocessing/wgms.py b/massbalancemachine/data_preprocessing/wgms.py index 87e06dd3..50dce20f 100644 --- a/massbalancemachine/data_preprocessing/wgms.py +++ b/massbalancemachine/data_preprocessing/wgms.py @@ -68,6 +68,9 @@ def build_monthly_data(data, cfg, rgi_region=None): ~((data.RGIId == "RGI60-01.23646") & (data.MONTH_DIFF == 24)) ] # West Yakutat + # Discard points before 1950 since ERA5 Land does not cover this period + data = data[data.YEAR > 1950] + region_name = get_region_name(rgi_region) dataset = Dataset(cfg, data=data, region_name=region_name, region_id=rgi_region) diff --git a/massbalancemachine/dataloader/DataLoader.py b/massbalancemachine/dataloader/DataLoader.py index a456452e..f781cbf7 100644 --- a/massbalancemachine/dataloader/DataLoader.py +++ b/massbalancemachine/dataloader/DataLoader.py @@ -69,7 +69,7 @@ def set_train_test_split( *, test_size: float = None, type_fold: str = "group-meas-id", - random_state: bool = False, + use_random_seed: bool = True, ) -> Tuple[Iterator[Any], Iterator[Any]]: """ Split the dataset into training and testing sets. @@ -77,6 +77,7 @@ def set_train_test_split( Args: test_size (float): Proportion of the dataset to include in the test split. type_fold (str): Type of splitting between train and test sets. Options are 'group-rgi','group-c_region' or 'group-meas-id'. + use_random_seed (bool): Whether the random seed should be used or not. This allows having reproducible splits. Returns: Tuple[Iterator[Any], Iterator[Any]]: Iterators for training and testing indices. @@ -96,13 +97,13 @@ def set_train_test_split( X, y, glacier_ids, stake_meas_id, regions = self._prepare_data_for_cv( self.data, self.meta_data_columns ) - if random_state == False: + if use_random_seed: gss = GroupShuffleSplit( n_splits=1, test_size=test_size, - random_state=self.random_seed, # commenting this improve randomness + random_state=self.random_seed, ) - elif random_state == True: + else: gss = GroupShuffleSplit( n_splits=1, test_size=test_size, diff --git a/massbalancemachine/models/TorchNeuralNetworkRegressor.py b/massbalancemachine/models/TorchNeuralNetworkRegressor.py index 6fcc9bfb..97db2067 100644 --- a/massbalancemachine/models/TorchNeuralNetworkRegressor.py +++ b/massbalancemachine/models/TorchNeuralNetworkRegressor.py @@ -7,13 +7,18 @@ def createModel(cfg, modelParams): nInp = len(cfg.featureColumns) + dropout = modelParams.get("dropout", 0.0) if modelParams["type"] == "sequential": assert len(modelParams["layers"]) > 0 l = [nn.Linear(nInp, modelParams["layers"][0])] for i in range(len(modelParams["layers"]) - 1): l.append(nn.ReLU()) + if dropout > 0: + l.append(nn.Dropout(dropout)) l.append(nn.Linear(modelParams["layers"][i], modelParams["layers"][i + 1])) l.append(nn.ReLU()) + if dropout > 0: + l.append(nn.Dropout(dropout)) l.append(nn.Linear(modelParams["layers"][-1], 1)) network = nn.Sequential(*l) return network @@ -178,12 +183,14 @@ def cumulative_pred(self): # TODO: implement pass - def evaluate_group_pred(self, geodataloader): + def evaluate_group_pred(self, geodataloader, val=False): grouped_ids = pd.DataFrame() with torch.no_grad(): - for g in geodataloader.glaciers(): + iterator = geodataloader.glaciersVal if val else geodataloader.glaciers + stakeMethod = geodataloader.stakesVal if val else geodataloader.stakes + for g in iterator(): # Get input features, metadata and ground truth - stakes, metadata, point_balance = geodataloader.stakes(g) + stakes, metadata, point_balance = stakeMethod(g) idAggr = metadata["ID"].values # Make prediction diff --git a/massbalancemachine/plots/__init__.py b/massbalancemachine/plots/__init__.py index 5a37f15f..3b8e9d61 100644 --- a/massbalancemachine/plots/__init__.py +++ b/massbalancemachine/plots/__init__.py @@ -8,5 +8,6 @@ from plots.profile_plots import profilePerGlacier from plots.temporal_plots import cumulatedMassChange from plots.input_plot import histogram_mb, scatterplot_mb +from plots.map_plots import mapGlacier from plots.train_plots import plot_training_history from plots.style import use_mbm_style, COLOR_ANNUAL, COLOR_WINTER diff --git a/massbalancemachine/plots/map_plots.py b/massbalancemachine/plots/map_plots.py index e69de29b..224f4b71 100644 --- a/massbalancemachine/plots/map_plots.py +++ b/massbalancemachine/plots/map_plots.py @@ -0,0 +1,68 @@ +import xarray as xr +import pyproj +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.colors as mcolors + +import data_processing + + +def mapGlacier(df, rgi_id, year, cfg, ax=None, max_abs=None, title=None, gdir=None): + + df_glacier_year = df[(df.RGIId == rgi_id) & (df.YEAR == year)] + + if gdir is None: + # Initialize the OGGM Config + data_processing.oggm_utils._initialize_oggm_config("") + gdir = data_processing.oggm_utils._initialize_glacier_directories( + [rgi_id], cfg + )[0] + + with xr.open_dataset(gdir.get_filepath("gridded_data")) as ds: + ds = ds.load() + + # Coordinate transformation from WGS84 to the projection of OGGM data + transf = pyproj.Transformer.from_proj( + pyproj.CRS.from_user_input("EPSG:4326"), + pyproj.CRS.from_user_input(ds.pyproj_srs), + always_xy=True, + ) + lon = df_glacier_year["POINT_LON"].to_numpy() + lat = df_glacier_year["POINT_LAT"].to_numpy() + x, y = transf.transform(lon, lat) + + # Convert projected coordinates to nearest grid indices + col = np.round((x - ds.x.values[0]) / (ds.x.values[1] - ds.x.values[0])).astype(int) + row = np.round((y - ds.y.values[0]) / (ds.y.values[1] - ds.y.values[0])).astype(int) + + # Make an empty grid matching OGGM's gridded_data + heat = xr.full_like(ds.topo, np.nan) + valid = (row >= 0) & (row < ds.sizes["y"]) & (col >= 0) & (col < ds.sizes["x"]) + assert all(valid), "Some of the projected points fall outside of the OGGM grid" + heat.values[row[valid], col[valid]] = df_glacier_year.loc[valid, "pred"].to_numpy() + + # Get background topography + smap = ds.salem.get_map(countries=False) + smap.set_shapefile(gdir.read_shapefile("outlines")) + smap.set_topography(ds.topo.data) + + # Build color normalization (white is MB=0) + max_abs = max_abs or df_glacier_year["pred"].abs().max() + norm = mcolors.TwoSlopeNorm(vmin=-max_abs, vcenter=0, vmax=max_abs) + + if ax is None: + fig, ax = plt.subplots(figsize=(9, 9)) + else: + fig = None + + # Plot annual MB + smap.set_cmap("RdBu") + smap.set_norm(norm) + smap.set_data(heat) + smap.plot(ax=ax) + smap.append_colorbar(ax=ax, label="Annual MB (m.w.e.)") + ax.set_title(title or f"{rgi_id} year {year}") + + plt.tight_layout() + + return fig diff --git a/massbalancemachine/plots/perf_plots.py b/massbalancemachine/plots/perf_plots.py index ff27733f..ece6dc56 100644 --- a/massbalancemachine/plots/perf_plots.py +++ b/massbalancemachine/plots/perf_plots.py @@ -291,9 +291,10 @@ def predVSTruthGlacierWide( geoErr, ax=None, title="Glacier wide MB", - ax_xlim=(-1.5, 0.5), + ax_xlim=(-1.5, 1.0), ax_ylim=(-1.5, 1.0), color="orange", + legend=False, ): if ax is None: @@ -305,7 +306,8 @@ def predVSTruthGlacierWide( ax.errorbar( geoTarget[g], geoPred[g], xerr=2 * geoErr[g], label=g, fmt="o", color=color ) - plt.text(geoTarget[g] + 0.02, geoPred[g] + 0.02, g, fontsize=10) + if legend: + plt.text(geoTarget[g] + 0.02, geoPred[g] + 0.02, g, fontsize=10) # Diagonal line pt = (0, 0) diff --git a/massbalancemachine/plots/temporal_plots.py b/massbalancemachine/plots/temporal_plots.py index a8dda1b3..b15bebc4 100644 --- a/massbalancemachine/plots/temporal_plots.py +++ b/massbalancemachine/plots/temporal_plots.py @@ -1,5 +1,6 @@ import numpy as np import matplotlib.pyplot as plt +from matplotlib.ticker import FormatStrFormatter from calendar import month_abbr import time @@ -68,7 +69,7 @@ def cumulatedMassChange( first_year = np.sort(monthly_df.YEAR.unique())[0] t = np.concatenate([[first_year], t]) y = np.concatenate([[0.0], y]) - ax.plot(t, np.cumsum(y), color=color_pred) + (line,) = ax.plot(t, np.cumsum(y), color=color_pred) nyear = monthly_df.YEAR.nunique() if geo is not None and test_gl in geo: @@ -95,6 +96,24 @@ def cumulatedMassChange( glacier_title = titles.get(test_gl) if titles is not None else None ax.set_title(glacier_title or test_gl.capitalize(), fontsize=20) + ax.tick_params(axis="x", labelsize=12) + ax.tick_params(axis="y", labelsize=12) + step_years_xticks = nyear // 10 + ax.set_xticks( + np.arange( + first_year, first_year + nyear + step_years_xticks, step_years_xticks + ) + ) + ax.xaxis.set_major_formatter(FormatStrFormatter("%.0f")) + + # Remove unused axes + for i in range(len(custom_order), len(axs)): + if isinstance(axs, list): + ax = axs[i] + else: + ax = axs.flatten()[i] + ax.set_visible(False) + # # Set axes limits # if ax_xlim is not None: # ax.set_xlim(ax_xlim) @@ -103,4 +122,4 @@ def cumulatedMassChange( plt.tight_layout() - return fig + return fig, line diff --git a/massbalancemachine/training/__init__.py b/massbalancemachine/training/__init__.py index 0eb11bf0..294a2b35 100644 --- a/massbalancemachine/training/__init__.py +++ b/massbalancemachine/training/__init__.py @@ -3,5 +3,6 @@ loadBestModel, compute_stake_loss, assessOnTest, + assessOnVal, eval_geodetic, ) diff --git a/massbalancemachine/training/training.py b/massbalancemachine/training/training.py index 32fefc0c..174e1960 100644 --- a/massbalancemachine/training/training.py +++ b/massbalancemachine/training/training.py @@ -226,9 +226,7 @@ def eval_geodetic(model, geo_dataloader, return_grid_pred=[]): # pbar = tqdm.tqdm(geo_dataloader.glaciersGeo(), total=geo_dataloader.lenGeo()) # for g in pbar: - pbar.set_description( - "Making geodetic pred for %s" % (current_g), refresh=True - ) + pbar.set_description("Geodetic pred for %s" % (current_g), refresh=True) pbar.update(1) # consume prefetched current batch if current_geo_future is not None: @@ -442,7 +440,7 @@ def assessOnTest(log_dir, model, geodataloader_test, light=False): ) = scores(predSummer, targetSummer) # Make plots - plot_pred_vs_obs(log_dir, targetAll, predAll, {"rmse": rmse, "mae": mae}) + plot_pred_vs_obs(log_dir, targetAll, predAll, {"rmse": rmse, "mae": mae, "r2": r2}) if not light: # Geodetic prediction @@ -483,6 +481,174 @@ def assessOnTest(log_dir, model, geodataloader_test, light=False): return stats +def assessOnVal(model, geodataloader, params, async_transfer=None): + wGeo = params["training"]["wGeo"] + scalingStakes = params["training"]["scalingStakes"] + statsVal = {} + + if async_transfer is None: + async_transfer = check_async_transfer_compatibility(geodataloader) + + model.eval() + with torch.no_grad(): + cntStake = 0 + cntGeo = 0 + lossStake = 0.0 + lossGeo = 0.0 + targetAll = torch.zeros(0) + predAll = torch.zeros(0) + periodAll = np.zeros( + 0, dtype=np.array(list(geodataloader.periodToInt.keys())).dtype + ) # Initialize with correct dtype + + glacier_iter = iter(geodataloader.glaciersVal()) + try: + current_g = next(glacier_iter) + except StopIteration: + current_g = None + + # geo future loaded at the first iteration + current_geo_future = None + if current_g is not None and wGeo > 0 and geodataloader.hasGeo(current_g): + current_geo_future = geodataloader.submit_geo(current_g) + + batch_idx = 0 + while current_g is not None: + # for g in geodataloader.glaciersVal(): + + # Look ahead and start loading geo for the next iteration + try: + next_g = next(glacier_iter) + except StopIteration: + next_g = None + + next_geo_future = None + if next_g is not None and wGeo > 0 and geodataloader.hasGeo(next_g): + next_geo_future = geodataloader.submit_geo(next_g) + + stakes, metadata, point_balance = geodataloader.stakesVal(current_g) + stakes = torch.tensor(stakes.astype(np.float32)).to(geodataloader.device) + point_balance = torch.tensor(point_balance.astype(np.float32)).to( + geodataloader.device + ) + l, ret, int_id = compute_stake_loss( + model, stakes, metadata, point_balance, returnPred=True + ) + target = ret["target"] + pred = ret["pred"] + targetAll = torch.concatenate((targetAll, target)) + predAll = torch.concatenate((predAll, pred)) + metadata = metadata.assign(ID_int=int_id) + grouped_ids = metadata.groupby("ID_int").agg({"PERIOD": "first"}) + periodAll = np.concatenate( + (periodAll, np.array(grouped_ids["PERIOD"].values)) + ) + + valScalingStakes = stakes.shape[0] if scalingStakes == "meas" else 1.0 + l = l * valScalingStakes + + lossStake += l + + if wGeo > 0 and geodataloader.hasGeo(current_g): + + # consume prefetched current batch + if current_geo_future is not None: + geoGrid, metadata, ygeo, errgeo, precomputed_meta = ( + current_geo_future.result() + ) + else: + geoGrid, metadata, ygeo, errgeo, precomputed_meta = ( + geodataloader.geo(current_g) + ) + + # geoGrid, metadata, ygeo, errgeo = geodataloader.geo(g) + # geoGrid = torch.tensor(geoGrid.astype(np.float32)).to( + # geodataloader.device + # ) + # ygeo = torch.tensor(ygeo.astype(np.float32)).to( + # geodataloader.device + # ) + # errgeo = torch.tensor(errgeo.astype(np.float32)).to( + # geodataloader.device + # ) + geoGrid = geoGrid.to(geodataloader.device, non_blocking=async_transfer) + ygeo = ygeo.to(geodataloader.device, non_blocking=async_transfer) + errgeo = errgeo.to(geodataloader.device, non_blocking=async_transfer) + geod_periods = geodataloader.periods_per_glacier[current_g] + lossGeo += compute_geo_loss( + model, + geoGrid, + metadata, + ygeo, + errgeo, + geod_periods, + precomputed_meta, + ) + cntGeo += 1 + + # Shift pipeline + current_g = next_g + current_geo_future = next_geo_future + batch_idx += 1 + + cntStake += valScalingStakes + lossStake /= cntStake + if wGeo > 0 and cntGeo > 0: + lossGeo /= cntGeo + loss = lossStake + wGeo * lossGeo + else: + lossGeo = torch.tensor(torch.nan) + loss = lossStake + mse, rmse, mae, pearson_corr, r2, bias = scores(predAll, targetAll) + + indAnnual = np.argwhere(periodAll == "annual")[:, 0] + indWinter = np.argwhere(periodAll == "winter")[:, 0] + predAnnual = predAll[indAnnual] + targetAnnual = targetAll[indAnnual] + predWinter = predAll[indWinter] + targetWinter = targetAll[indWinter] + ( + mse_annual, + rmse_annual, + mae_annual, + pearson_corr_annual, + r2_annual, + bias_annual, + ) = scores(predAnnual, targetAnnual) + ( + mse_winter, + rmse_winter, + mae_winter, + pearson_corr_winter, + r2_winter, + bias_winter, + ) = scores(predWinter, targetWinter) + + statsVal["lossValStake"] = lossStake.item() + statsVal["lossValGeo"] = lossGeo.item() + statsVal["lossVal"] = loss.item() + statsVal["mse"] = mse.item() + statsVal["rmse"] = rmse.item() + statsVal["mae"] = mae.item() + statsVal["pearson"] = pearson_corr.item() + statsVal["r2"] = r2.item() + statsVal["bias"] = bias.item() + statsVal["mse_annual"] = mse_annual.item() + statsVal["rmse_annual"] = rmse_annual.item() + statsVal["mae_annual"] = mae_annual.item() + statsVal["pearson_annual"] = pearson_corr_annual.item() + statsVal["r2_annual"] = r2_annual.item() + statsVal["bias_annual"] = bias_annual.item() + statsVal["mse_winter"] = mse_winter.item() + statsVal["rmse_winter"] = rmse_winter.item() + statsVal["mae_winter"] = mae_winter.item() + statsVal["pearson_winter"] = pearson_corr_winter.item() + statsVal["r2_winter"] = r2_winter.item() + statsVal["bias_winter"] = bias_winter.item() + + return statsVal + + def plot_pred_vs_obs(log_dir, target, pred, scores): grouped_ids = pd.DataFrame({"target": target.numpy(), "pred": pred.numpy()}) fig = predVSTruth( @@ -528,7 +694,7 @@ def loadBestModel(log_dir, model): if best is None: raise Exception("No model found.") model.load_state_dict(torch.load(files[best], weights_only=True)) - return files[best] + return files[best], bestVal def train_geo( @@ -564,17 +730,20 @@ def train_geo( scalingStakes = params["training"]["scalingStakes"] assert scalingStakes in ["meas", "glacier"] iterPerEpoch = len(geodataloader) - nColsProgressBar = 500 if _inJupyterNotebook else 100 + nColsProgressBar = 500 if _inJupyterNotebook else 85 # Setup logging run_name = datetime.datetime.now().strftime("%Y%m%d_%H%M%S") log_suffix = params["training"]["log_suffix"] if log_suffix != "": log_suffix = "_" + log_suffix + log_prefix = params["training"].get("log_prefix", "") + if log_prefix != "": + log_prefix = log_prefix + "_" if params["training"]["log_dir"] is None: log_dir = os.path.join( _default_log_dir, - f"geo_{run_name}{log_suffix}", + f"{log_prefix}geo_{run_name}{log_suffix}", ) else: log_dir = params["training"]["log_dir"] @@ -620,6 +789,7 @@ def train_geo( ): avg_training_loss = 0.0 + iter_counter = 0 model.train() with tqdm.tqdm( total=iterPerEpoch, @@ -699,6 +869,7 @@ def train_geo( torch.cuda.synchronize() stakesInferenceTime = time.time() - st + # TODO: is it scaled by the number of measurements properly? valScalingStakes = ( stakes.shape[0] if scalingStakes == "meas" else 1.0 ) @@ -786,7 +957,10 @@ def train_geo( statsTraining["lossGeo"].append(lossGeo.item()) statsTraining["loss"].append(loss.item()) statsTraining["lr"].append(lr) - avg_training_loss += statsTraining["loss"][-1] + if not torch.isnan(loss): + # This happens when there is no data to compute the loss (no geodetic term + no stakes) + avg_training_loss += loss.item() + iter_counter += 1 # Log to TensorBoard globalStep = iterPerEpoch * epoch + batch_idx @@ -820,279 +994,100 @@ def train_geo( batch_idx += 1 batch_bar.set_postfix( - loss=f"{statsTraining['loss'][-1]:.4f}", - lossStake=f"{statsTraining['lossStake'][-1]:.4f}", - lossGeo=f"{statsTraining['lossGeo'][-1]:.4f}", + l=f"{statsTraining['loss'][-1]:.3f}", + lStake=f"{statsTraining['lossStake'][-1]:.3f}", + lGeo=f"{statsTraining['lossGeo'][-1]:.3f}", ) batch_bar.update(1) - avg_training_loss /= iterPerEpoch + avg_training_loss /= iter_counter if scheduler is not None: scheduler.step() if freqVal and geodataloader.lenVal() > 0: - model.eval() - with torch.no_grad(): - cntStake = 0 - cntGeo = 0 - lossStake = 0.0 - lossGeo = 0.0 - targetAll = torch.zeros(0) - predAll = torch.zeros(0) - periodAll = np.zeros( - 0, dtype=np.array(list(geodataloader.periodToInt.keys())).dtype - ) # Initialize with correct dtype - - glacier_iter = iter(geodataloader.glaciersVal()) - try: - current_g = next(glacier_iter) - except StopIteration: - current_g = None - - # geo future loaded at the first iteration - current_geo_future = None - if ( - current_g is not None - and wGeo > 0 - and geodataloader.hasGeo(current_g) - ): - current_geo_future = geodataloader.submit_geo(current_g) - - batch_idx = 0 - while current_g is not None: - # for g in geodataloader.glaciersVal(): - - # Look ahead and start loading geo for the next iteration - try: - next_g = next(glacier_iter) - except StopIteration: - next_g = None - - next_geo_future = None - if ( - next_g is not None - and wGeo > 0 - and geodataloader.hasGeo(next_g) - ): - next_geo_future = geodataloader.submit_geo(next_g) - - stakes, metadata, point_balance = geodataloader.stakesVal( - current_g - ) - stakes = torch.tensor(stakes.astype(np.float32)).to( - geodataloader.device - ) - point_balance = torch.tensor( - point_balance.astype(np.float32) - ).to(geodataloader.device) - l, ret, int_id = compute_stake_loss( - model, stakes, metadata, point_balance, returnPred=True - ) - target = ret["target"] - pred = ret["pred"] - targetAll = torch.concatenate((targetAll, target)) - predAll = torch.concatenate((predAll, pred)) - metadata = metadata.assign(ID_int=int_id) - grouped_ids = metadata.groupby("ID_int").agg( - {"PERIOD": "first"} - ) - periodAll = np.concatenate( - (periodAll, np.array(grouped_ids["PERIOD"].values)) - ) - - valScalingStakes = ( - stakes.shape[0] if scalingStakes == "meas" else 1.0 - ) - l = l * valScalingStakes - - lossStake += l - - if wGeo > 0 and geodataloader.hasGeo(current_g): - - # consume prefetched current batch - if current_geo_future is not None: - geoGrid, metadata, ygeo, errgeo, precomputed_meta = ( - current_geo_future.result() - ) - else: - geoGrid, metadata, ygeo, errgeo, precomputed_meta = ( - geodataloader.geo(current_g) - ) + statsValEpoch = assessOnVal( + model, geodataloader, params, async_transfer=async_transfer + ) + metric2tbEntry = { + "lossValStake": "LossStake/val", + "lossValGeo": "LossGeo/val", + "lossVal": "Loss/val", + "rmse": "RMSE/val", + "rmse_annual": "RMSE_annual/val", + "rmse_winter": "RMSE_winter/val", + "mae": "MAE/val", + "mae_annual": "MAE_annual/val", + "mae_winter": "MAE_winter/val", + "pearson": "Pearson/val", + "pearson_annual": "Pearson_annual/val", + "pearson_winter": "Pearson_winter/val", + "r2": "R2/val", + "r2_annual": "R2_annual/val", + "r2_winter": "R2_winter/val", + "bias": "bias/val", + "bias_annual": "bias_annual/val", + "bias_winter": "bias_winter/val", + } + for k, v in statsValEpoch.items(): + statsVal[k].append(v) + # Log to TensorBoard + if k in metric2tbEntry: + writer.add_scalar(metric2tbEntry[k], statsVal[k][-1], epoch) + + # Check if current model is among the top 5 + scalarBestModelCriterion = ( + statsVal[bestModelCriterion][-1] + if len(statsVal[bestModelCriterion]) > 0 + else np.nan + ) + assert not np.isnan( + scalarBestModelCriterion + ), f"The statistics {bestModelCriterion} used for best model selection contains NaN." + higherBetterCriterion = ( + 1 + if any(crit in bestModelCriterion for crit in _maxCriterion) + else -1 + ) + model_path = os.path.join( + log_dir, + f"model_epoch{epoch}_{bestModelCriterion}{scalarBestModelCriterion:.4f}.pt", + ) - # geoGrid, metadata, ygeo, errgeo = geodataloader.geo(g) - # geoGrid = torch.tensor(geoGrid.astype(np.float32)).to( - # geodataloader.device - # ) - # ygeo = torch.tensor(ygeo.astype(np.float32)).to( - # geodataloader.device - # ) - # errgeo = torch.tensor(errgeo.astype(np.float32)).to( - # geodataloader.device - # ) - geoGrid = geoGrid.to( - geodataloader.device, non_blocking=async_transfer - ) - ygeo = ygeo.to( - geodataloader.device, non_blocking=async_transfer - ) - errgeo = errgeo.to( - geodataloader.device, non_blocking=async_transfer - ) - geod_periods = geodataloader.periods_per_glacier[current_g] - lossGeo += compute_geo_loss( - model, - geoGrid, - metadata, - ygeo, - errgeo, - geod_periods, - precomputed_meta, - ) - cntGeo += 1 - - # Shift pipeline - current_g = next_g - current_geo_future = next_geo_future - batch_idx += 1 - - cntStake += valScalingStakes - lossStake /= cntStake - if wGeo > 0 and cntGeo > 0: - lossGeo /= cntGeo - loss = lossStake + wGeo * lossGeo - else: - lossGeo = torch.tensor(torch.nan) - loss = lossStake - mse, rmse, mae, pearson_corr, r2, bias = scores(predAll, targetAll) - - indAnnual = np.argwhere(periodAll == "annual")[:, 0] - indWinter = np.argwhere(periodAll == "winter")[:, 0] - predAnnual = predAll[indAnnual] - targetAnnual = targetAll[indAnnual] - predWinter = predAll[indWinter] - targetWinter = targetAll[indWinter] - ( - mse_annual, - rmse_annual, - mae_annual, - pearson_corr_annual, - r2_annual, - bias_annual, - ) = scores(predAnnual, targetAnnual) - ( - mse_winter, - rmse_winter, - mae_winter, - pearson_corr_winter, - r2_winter, - bias_winter, - ) = scores(predWinter, targetWinter) - - statsVal["lossValStake"].append(lossStake.item()) - statsVal["lossValGeo"].append(lossGeo.item()) - statsVal["lossVal"].append(loss.item()) - statsVal["mse"].append(mse.item()) - statsVal["rmse"].append(rmse.item()) - statsVal["mae"].append(mae.item()) - statsVal["pearson"].append(pearson_corr.item()) - statsVal["r2"].append(r2.item()) - statsVal["bias"].append(bias.item()) - statsVal["mse_annual"].append(mse_annual.item()) - statsVal["rmse_annual"].append(rmse_annual.item()) - statsVal["mae_annual"].append(mae_annual.item()) - statsVal["pearson_annual"].append(pearson_corr_annual.item()) - statsVal["r2_annual"].append(r2_annual.item()) - statsVal["bias_annual"].append(bias_annual.item()) - statsVal["mse_winter"].append(mse_winter.item()) - statsVal["rmse_winter"].append(rmse_winter.item()) - statsVal["mae_winter"].append(mae_winter.item()) - statsVal["pearson_winter"].append(pearson_corr_winter.item()) - statsVal["r2_winter"].append(r2_winter.item()) - statsVal["bias_winter"].append(bias_winter.item()) + if ( + len(top_models) < 5 + or scalarBestModelCriterion + < higherBetterCriterion * top_models[0][0] + ): - # Log to TensorBoard - writer.add_scalar("LossStake/val", lossStake.item(), epoch) - writer.add_scalar("LossGeo/val", lossGeo.item(), epoch) - writer.add_scalar("Loss/val", loss.item(), epoch) - writer.add_scalar("RMSE/val", rmse.item(), epoch) - writer.add_scalar("RMSE_annual/val", rmse_annual.item(), epoch) - writer.add_scalar("RMSE_winter/val", rmse_winter.item(), epoch) - writer.add_scalar("MAE/val", mae.item(), epoch) - writer.add_scalar("MAE_annual/val", mae_annual.item(), epoch) - writer.add_scalar("MAE_winter/val", mae_winter.item(), epoch) - writer.add_scalar("Pearson/val", pearson_corr.item(), epoch) - writer.add_scalar( - "Pearson_annual/val", pearson_corr_annual.item(), epoch - ) - writer.add_scalar( - "Pearson_winter/val", pearson_corr_winter.item(), epoch - ) - writer.add_scalar("R2/val", r2.item(), epoch) - writer.add_scalar("R2_annual/val", r2_annual.item(), epoch) - writer.add_scalar("R2_winter/val", r2_winter.item(), epoch) - writer.add_scalar("bias/val", bias.item(), epoch) - writer.add_scalar("bias_annual/val", bias_annual.item(), epoch) - writer.add_scalar("bias_winter/val", bias_winter.item(), epoch) - - # Check if current model is among the top 5 - scalarBestModelCriterion = ( - statsVal[bestModelCriterion][-1] - if len(statsVal[bestModelCriterion]) > 0 - else np.nan - ) - assert not np.isnan( - scalarBestModelCriterion - ), f"The statistics {bestModelCriterion} used for best model selection contains NaN." - higherBetterCriterion = ( - 1 - if any(crit in bestModelCriterion for crit in _maxCriterion) - else -1 - ) - model_path = os.path.join( - log_dir, - f"model_epoch{epoch}_{bestModelCriterion}{scalarBestModelCriterion:.4f}.pt", + if len(top_models) >= 5: + # Pop the worst among top 5 + _, worst_path = heapq.heappop(top_models) + if os.path.exists(worst_path): + os.remove(worst_path) + + # Save the new model + # Use -scalarBestModelCriterion for max-heap + heapq.heappush( + top_models, + ( + higherBetterCriterion * scalarBestModelCriterion, + model_path, + ), ) + torch.save(model.state_dict(), model_path) - if ( - len(top_models) < 5 - or scalarBestModelCriterion - < higherBetterCriterion * top_models[0][0] - ): - - if len(top_models) >= 5: - # Pop the worst among top 5 - _, worst_path = heapq.heappop(top_models) - if os.path.exists(worst_path): - os.remove(worst_path) - - # Save the new model - # Use -scalarBestModelCriterion for max-heap - heapq.heappush( - top_models, - ( - higherBetterCriterion * scalarBestModelCriterion, - model_path, - ), - ) - torch.save(model.state_dict(), model_path) - - if ( - geodataloader_test is not None - and len(geodataloader_test) > 0 - ): - assessOnTest(log_dir, model, geodataloader_test, light=True) + if geodataloader_test is not None and len(geodataloader_test) > 0: + assessOnTest(log_dir, model, geodataloader_test, light=True) rmse = statsVal["rmse"][-1] if len(statsVal["rmse"]) > 0 else np.nan - pearson = ( - statsVal["pearson"][-1] if len(statsVal["pearson"]) > 0 else np.nan - ) + r2 = statsVal["r2"][-1] if len(statsVal["r2"]) > 0 else np.nan loss = statsVal["lossVal"][-1] if len(statsVal["lossVal"]) > 0 else np.nan geodataloader.onEpochEnd() # Show progress bar tqdm.tqdm.write( - f"[Epoch {epoch+1} / {Nepochs}] Avg training loss: {avg_training_loss:.3f}, Val loss: {loss:.3f}, RMSE: {rmse:.3f}, Pearson: {pearson:.3f}" + f"[Epoch {epoch+1} / {Nepochs}] Avg training loss: {avg_training_loss:.3f}, Loss: {loss:.3f}, RMSE: {rmse:.3f}, R2: {r2:.3f}" ) except KeyboardInterrupt: diff --git a/notebooks/data_processing_wgms.ipynb b/notebooks/data_processing_wgms.ipynb index 74de67de..bf4e0c83 100644 --- a/notebooks/data_processing_wgms.ipynb +++ b/notebooks/data_processing_wgms.ipynb @@ -86,9 +86,10 @@ "metadata": {}, "outputs": [], "source": [ - "# Get the RGI ID for each stake measurement for the region of interest\n", - "data = mbm.data_processing.utils.get_rgi(data=data, region=6)\n", - "data" + "# # Get the RGI ID for each stake measurement for the region of interest\n", + "# # The shape file is automatically retrieved with OGMM in the background\n", + "# data = mbm.data_processing.utils.get_rgi(data=data, region=6)\n", + "# data" ] }, { @@ -98,13 +99,14 @@ "metadata": {}, "outputs": [], "source": [ - "# # Specify the shape filename of the glaciers outline obtained from RGIv6\n", - "# glacier_outline_fname = 'path_to_your_own_shapefiles_folder/06_rgi60_Iceland.shp'\n", - "# # Load the glacier outlines\n", - "# glacier_outline = gpd.read_file(glacier_outline_fname)\n", - "# # Get the RGI ID for each stake measurement for the region of interest\n", - "# data = mbm.data_processing.utils.get_rgi(data=data, glacier_outlines=glacier_outline)\n", - "# data" + "# Specify the shape filename of the glaciers outline obtained from RGIv6\n", + "# The shape file is specified by the user\n", + "glacier_outline_fname = repoPath + 'notebooks/example_data/iceland/files/06_rgi62_Iceland.shp'\n", + "# Load the glacier outlines\n", + "glacier_outline = gpd.read_file(glacier_outline_fname)\n", + "# Get the RGI ID for each stake measurement for the region of interest\n", + "data = mbm.data_processing.utils.get_rgi(data=data, glacier_outlines=glacier_outline)\n", + "data" ] }, { diff --git a/notebooks/example_data/iceland/files/06_rgi62_Iceland.cpg b/notebooks/example_data/iceland/files/06_rgi62_Iceland.cpg new file mode 100644 index 00000000..cd89cb97 --- /dev/null +++ b/notebooks/example_data/iceland/files/06_rgi62_Iceland.cpg @@ -0,0 +1 @@ +ISO-8859-1 \ No newline at end of file diff --git a/notebooks/example_data/iceland/files/06_rgi62_Iceland.dbf b/notebooks/example_data/iceland/files/06_rgi62_Iceland.dbf new file mode 100644 index 00000000..3d3e075b Binary files /dev/null and b/notebooks/example_data/iceland/files/06_rgi62_Iceland.dbf differ diff --git a/notebooks/example_data/iceland/files/06_rgi62_Iceland.prj b/notebooks/example_data/iceland/files/06_rgi62_Iceland.prj new file mode 100644 index 00000000..a30c00a5 --- /dev/null +++ b/notebooks/example_data/iceland/files/06_rgi62_Iceland.prj @@ -0,0 +1 @@ +GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]] \ No newline at end of file diff --git a/notebooks/example_data/iceland/files/06_rgi62_Iceland.shp b/notebooks/example_data/iceland/files/06_rgi62_Iceland.shp new file mode 100644 index 00000000..905c00c1 Binary files /dev/null and b/notebooks/example_data/iceland/files/06_rgi62_Iceland.shp differ diff --git a/notebooks/example_data/iceland/files/06_rgi62_Iceland.shx b/notebooks/example_data/iceland/files/06_rgi62_Iceland.shx new file mode 100644 index 00000000..b734cf7f Binary files /dev/null and b/notebooks/example_data/iceland/files/06_rgi62_Iceland.shx differ diff --git a/notebooks/model_training_neuralnetwork.ipynb b/notebooks/model_training_neuralnetwork.ipynb index 4a8cff54..c9eae747 100644 --- a/notebooks/model_training_neuralnetwork.ipynb +++ b/notebooks/model_training_neuralnetwork.ipynb @@ -236,10 +236,10 @@ "metadata": {}, "outputs": [], "source": [ - "# For each of the XGBoost parameter, define the grid range\n", + "# For each of the neural network hyper-parameters, define the grid range\n", "parameters = {\n", - " 'lr': [0.005, 0.001, 0.0005],\n", - " 'max_epochs': [500, 1000],\n", + " 'lr': [0.005, 0.001],\n", + " 'max_epochs': [300, 600],\n", " 'optimizer': [torch.optim.Adam],\n", "}" ] diff --git a/notebooks/model_training_xgboost.ipynb b/notebooks/model_training_xgboost.ipynb index ddead639..f5674bbb 100644 --- a/notebooks/model_training_xgboost.ipynb +++ b/notebooks/model_training_xgboost.ipynb @@ -123,10 +123,9 @@ " 3,\n", " 4,\n", " 5,\n", - " 6,\n", " ],\n", - " 'learning_rate': [0.01, 0.1, 0.2, 0.3],\n", - " 'n_estimators': [100, 200, 300],\n", + " 'learning_rate': [0.01, 0.1],\n", + " 'n_estimators': [200, 300],\n", " 'gamma': [0, 1]\n", "}" ] @@ -172,7 +171,7 @@ "# versus runtime.\n", "custom_xgboost.randomsearch(\n", " parameters=parameters,\n", - " n_iter=20,\n", + " n_iter=5,\n", " splits=splits,\n", " features=df_X_train,\n", " targets=y_train,\n", diff --git a/pyproject.toml b/pyproject.toml index f80fe9a7..af588d88 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -25,7 +25,6 @@ xarray = "^2024.6.0" bottleneck = "^1.4.0" scikit-learn = "^1.5.1" pyproj = "^3.6.1" -cartopy = "^0.23.0" rasterio = "^1.3.10" geopandas = "^1.0.1" rioxarray = "^0.17.0" @@ -40,7 +39,6 @@ joblib = "^1.4.2" progressbar2 = "^4.4.2" oggm = "^1.6.1" xgboost = "^2.1.0" -contextily = "^1.6.0" dill = "^0.3.8" pdoc = "^14.6.1" tqdm = "^4.66.5" @@ -52,20 +50,18 @@ pytest = "^8.2.2" pytest-cov = "^6.1.1" zarr = "2.18.4" pyarrow = "^19.0.1" -tensorboard = "^2.19.0" torcheval = "^0.0.7" pip-tools = "7.5.0" pandoc = "2.4" nbconvert = "7.16.6" -sphinxemoji = "^0.3.0" numcodecs = "0.15.1" torchsummary = "^1.5" nbstripout = "^0.8.1" -joypy = "^0.2.6" -folium = "^0.20.0" -geomloss = "^0.2.6" -pykeops = "^2.3" setuptools = "81" +tensorboard = "^2.19.0" + +[tool.poetry.group.docs] +optional = false [tool.poetry.group.docs.dependencies] sphinx = "^7.4.6" @@ -75,10 +71,23 @@ sphinx-intl = "^2.2.0" sphinx-reredirects = "^0.1.5" sphinx-autobuild = "2025.8.25" sphinx-autorun = "^1.1.1" +sphinxemoji = "^0.3.0" myst-parser = "^2.0" nbsphinx = "^0.9" ipykernel = "^6.0" -pandoc = "2.4" +jupyterlab-execute-time = "^3.3.0" + +[tool.poetry.group.viz] +optional = false + +[tool.poetry.group.viz.dependencies] +cartopy = "^0.23.0" +contextily = "^1.6.0" +joypy = "^0.2.6" +folium = "^0.20.0" + +[tool.poetry.group.dev] +optional = false [tool.poetry.group.dev.dependencies] pytest = "^8.2.2" @@ -91,6 +100,22 @@ ruff = "^0.6.2" nbdime = "^4.0.1" nbdev = "^2.3.28" +[tool.poetry.group.advanced] +optional = false + +[tool.poetry.group.advanced.dependencies] +geomloss = "^0.2.6" +pykeops = "^2.3" +optuna = "^4.8.0" +plotly = "^6.7.0" +optuna-dashboard = "^0.20.0" + +[tool.poetry.group.publish] +optional = true + +[tool.poetry.group.publish.dependencies] +huggingface-hub = "^1.10.2" + [build-system] requires = ["poetry-core>=1.0.0"] build-backend = "poetry.core.masonry.api" diff --git a/scripts/common.py b/scripts/common.py index 134a2b5f..f0fbe6c1 100644 --- a/scripts/common.py +++ b/scripts/common.py @@ -6,6 +6,8 @@ import yaml import pandas as pd from sklearn.model_selection import train_test_split +from collections.abc import Mapping, Sequence +import math import massbalancemachine as mbm @@ -30,6 +32,7 @@ def parseParams(params): "model": { "type": params["model"]["type"], "layers": params["model"]["layers"], + "dropout": params["model"].get("dropout", 0.0), "inputs": inputs, "downscale": downscale, }, @@ -53,8 +56,10 @@ def parseParams(params): "bestModelCriterion": params["training"].get( "bestModelCriterion", "lossVal" ), + "splitVal": params["training"].get("splitVal", "group-meas-id"), "freqVal": params["training"].get("freqVal", 1), "log_suffix": params["training"].get("log_suffix", ""), + "log_prefix": params["training"].get("log_prefix", ""), "log_dir": params["training"].get("log_dir"), }, } @@ -68,3 +73,98 @@ def loadParams(modelType): print(exc) parsedParams = parseParams(params) return parsedParams + + +def default_glacier_name(rgi_id): + return { + # # Norway + # RGI60-08.00038; # Nigardsbreen + # RGI60-08.00087; # Jostedalsbreen + # RGI60-08.00147; # Folgefonna + # RGI60-08.00203; # Hardangerjøkulen + # Italy + "RGI60-11.00695": "Glatschiu dil segnas", + "RGI60-11.03005": "Miage", + "RGI60-11.03001": "Brenva", + "RGI60-11.01473": "Laaser Ferner", + "RGI60-11.00597": "Übeltalferner", + "RGI60-11.01776": "Langenferner/Vedretta Lunga", + # France, Mont Blanc + "RGI60-11.03643": "Mer de Glace/Geant", + "RGI60-11.03638": "Argentière", + "RGI60-11.03646": "Bossons", + "RGI60-11.03647": "Taconnaz", + "RGI60-11.03296": "Tricot", + "RGI60-11.03438": "Tete Rousse", + "RGI60-11.03648": "Bionnassay", + "RGI60-11.03601": "Armancette", + "RGI60-11.03650": "Covagnet", + "RGI60-11.03276": "Miage 1", + "RGI60-11.03388": "Miage 2", + "RGI60-11.03579": "Miage 3", + "RGI60-11.03649": "Miage 4", + "RGI60-11.03651": "Tré-la-Tête", + "RGI60-11.03339": "Glaciers", + # France, Belledonne + "RGI60-11.03674": "Saint Sorlin", + # France, Ecrins + "RGI60-11.03677": "Meije", + "RGI60-11.03684": "Blanc", + # France, Pyrénées + "RGI60-11.03232": "Ossoue", + "RGI60-11.03208": "Aneto", + # Austria + "RGI60-11.00897": "Hintereisferner", + "RGI60-11.00787": "Kesselwandferner", + # Switzerland + "RGI60-11.01270": "Grindelwald", + "RGI60-11.01450": "Aletsch", + "RGI60-11.01733": "Hangend", + "RGI60-11.01328": "Unteraar", + "RGI60-11.01238": "Rhone", + "RGI60-11.02249": "Tsanfleuron", + "RGI60-11.01702": "Kander", + "RGI60-11.00872": "Hüfi", + "RGI60-11.02774": "Giétro", + "RGI60-11.01876": "Gries", + "RGI60-11.02746": "Schwarzberg", + "RGI60-11.02810": "Arolla", + "RGI60-11.02775": "Orny", + "RGI60-11.02507": "Brunegg", + "RGI60-11.00804": "Silvretta", + "RGI60-11.00752": "Vorab", + "RGI60-11.02787": "Mont Collon", + "RGI60-11.01267": "Porchabella", + "RGI60-11.02634": "Prafleuri", + "RGI60-11.01946": "Morteratsch", + }.get(rgi_id) + + +def canonicalize(x): + """Convert params into a stable, comparable representation.""" + if isinstance(x, Mapping): + return tuple(sorted((str(k), canonicalize(v)) for k, v in x.items())) + if isinstance(x, tuple): + return tuple(canonicalize(v) for v in x) + if isinstance(x, list): + return tuple(canonicalize(v) for v in x) + if isinstance(x, float): + if math.isnan(x): + return "__nan__" + return x + return x + + +def already_completed_trial(study, candidate_params: dict): + from optuna.trial import TrialState + + candidate_key = canonicalize(candidate_params) + + for t in study.get_trials( + deepcopy=False, + states=(TrialState.COMPLETE,), + ): + if canonicalize(t.params) == candidate_key: + return True, t + + return False, None diff --git a/scripts/export_to_huggingface.py b/scripts/export_to_huggingface.py new file mode 100644 index 00000000..a8fe252b --- /dev/null +++ b/scripts/export_to_huggingface.py @@ -0,0 +1,46 @@ +import sys, os + +mbm_path = os.path.abspath(os.path.join(os.path.dirname(__file__), "../")) +sys.path.append(mbm_path) # Add root of repo to import MBM + +from pathlib import Path +import shutil +import argparse +from huggingface_hub import login, upload_folder, HfApi + +parser = argparse.ArgumentParser() +parser.add_argument("modelFolder", type=str, help="Folder of the model to export.") +parser.add_argument("tag", type=str, help="Tag to identify the model to publish.") +args = parser.parse_args() + +modelFolder = args.modelFolder +tag = args.tag + +repo_id = "MassBalanceMachine/MLP" + +api = HfApi() +login() + +source_folder = os.path.join(mbm_path, "logs", modelFolder) +source_model = os.path.join(source_folder, "best_model.json") +source_params = os.path.join(source_folder, "params.json") +if not os.path.isfile(source_model): + source_model = os.path.join(source_folder, "model.json") + assert os.path.isfile(source_model), f"The model file does not exist." +assert os.path.isfile(source_params), f"File {source_params} does not exist." +assert len(tag) > 0, "Tag cannot be empty" + +publish_folder = os.path.join(mbm_path, "publish", modelFolder) +print(f"{publish_folder=}") +if os.path.isdir(publish_folder): + shutil.rmtree(publish_folder) +os.makedirs(publish_folder, exist_ok=True) +shutil.copy2(source_model, os.path.join(publish_folder, "model.json")) +shutil.copy2(source_params, publish_folder) + + +commit_info = upload_folder( + folder_path=publish_folder, repo_id=repo_id, repo_type="model" +) +api.create_tag(repo_id=repo_id, tag=tag, revision=commit_info.oid) +print(f"Model was published under the tag {tag}") diff --git a/scripts/geo/compare.py b/scripts/geo/compare.py index 7696e6a6..e95ab8d2 100644 --- a/scripts/geo/compare.py +++ b/scripts/geo/compare.py @@ -25,13 +25,18 @@ parser = argparse.ArgumentParser("Compare two different models.") parser.add_argument("modelFolder1", type=str, help="Folder of the 1st model to load.") parser.add_argument("modelFolder2", type=str, help="Folder of the 2nd model to load.") -# parser.add_argument( -# "--cpu", -# dest="cpu", -# default=False, -# action="store_true", -# help="Force model to run on CPU, even if a GPU is available.", -# ) +parser.add_argument( + "--name1", + dest="name1", + default=None, + help="Optional name for the 1st model.", +) +parser.add_argument( + "--name2", + dest="name2", + default=None, + help="Optional name for the 2nd model.", +) parser.add_argument( "--plot", dest="plot", @@ -39,30 +44,41 @@ action="store_true", help="Display figures in addition to saving.", ) -# parser.add_argument( -# "--noTest", -# dest="noTest", -# default=False, -# action="store_true", -# help="Do not evaluate on test data.", -# ) -# parser.add_argument( -# "--onRegion", -# dest="onRegion", -# default=False, -# action="store_true", -# help="Evaluate prediction on the whole region in addition to classical plots.", -# ) +parser.add_argument( + "--noTrain", + dest="noTrain", + default=False, + action="store_true", + help="Do not compare on train data.", +) +parser.add_argument( + "--maps", + dest="maps", + default=False, + action="store_true", + help="Generate annual MB maps for test glaciers.", +) +parser.add_argument( + "--mapsTrain", + dest="mapsTrain", + default=[], + nargs="+", + help="Generate annual MB maps for specific train glaciers.", +) args = parser.parse_args() modelFolder1 = args.modelFolder1 modelFolder2 = args.modelFolder2 -# cpu = args.cpu +name1 = args.name1 +name2 = args.name2 plot = args.plot -# noTest = args.noTest -# onRegion = args.onRegion +noTrain = args.noTrain +maps = args.maps +mapsTrain = args.mapsTrain pathFolder1 = os.path.join("logs", modelFolder1) pathFolder2 = os.path.join("logs", modelFolder2) +name1 = name1 if name1 is not None else modelFolder1 +name2 = name2 if name2 is not None else modelFolder2 pathFolder = os.path.join("results/comp/", f"{modelFolder1}_{modelFolder2}") os.makedirs(pathFolder, exist_ok=True) @@ -72,36 +88,89 @@ with open(f"{pathFolder2}/params.json", "r") as f: params2 = json.load(f) +with open(f"{pathFolder1}/glacierNames.json", "r") as f: + glacierNames = json.load(f) -# Cumulated mass change on train data -df_gridded_monthly1 = pd.read_csv(f"{pathFolder1}/gridded_monthly_train.csv") -df_geo1 = pd.read_csv(f"{pathFolder1}/gridded_geodetic_train.csv") -geoTarget = df_geo1.set_index("RGIId").target.to_dict() -geoErr = df_geo1.set_index("RGIId").err.to_dict() -# Plot cumulated mass change -fig = mbm.plots.cumulatedMassChange( - df_gridded_monthly1, - geo={ - rgi_id: {"mean": geoTarget[rgi_id], "err": geoErr[rgi_id]} - for rgi_id in geoTarget - }, -) -del df_gridded_monthly1 -df_gridded_monthly2 = pd.read_csv(f"{pathFolder2}/gridded_monthly_train.csv") -mbm.plots.cumulatedMassChange( - df_gridded_monthly2, - geo=None, - axs=fig.axes, - color_pred="orange", -) +if not noTrain: + # Cumulated mass change on train data + df_gridded_monthly1 = pd.read_csv(f"{pathFolder1}/gridded_monthly_train.csv") + df_geo1 = pd.read_csv(f"{pathFolder1}/gridded_geodetic_train.csv") + geoTarget = df_geo1.set_index("RGIId").target.to_dict() + geoErr = df_geo1.set_index("RGIId").err.to_dict() + # Plot cumulated mass change + fig, l1 = mbm.plots.cumulatedMassChange( + df_gridded_monthly1, + geo={ + rgi_id: {"mean": geoTarget[rgi_id], "err": geoErr[rgi_id]} + for rgi_id in geoTarget + }, + ) + del df_gridded_monthly1 + df_gridded_monthly2 = pd.read_csv(f"{pathFolder2}/gridded_monthly_train.csv") + _, l2 = mbm.plots.cumulatedMassChange( + df_gridded_monthly2, + geo=None, + axs=fig.axes, + color_pred="orange", + titles={ + k: (f"{k} ({glacierNames[k]})" if glacierNames[k] is not None else None) + for k in glacierNames + }, + ) + del df_gridded_monthly2 + fig.legend([l1, l2], [name1, name2], loc="lower center", ncol=2) -fig.savefig(f"{pathFolder}/cumulated_mass_change_glaciers_train.pdf") -fig.savefig(f"{pathFolder}/cumulated_mass_change_glaciers_train.png", dpi=300) -if plot: - plt.show() -plt.close(fig) + fig.savefig(f"{pathFolder}/cumulated_mass_change_glaciers_train.pdf") + fig.savefig(f"{pathFolder}/cumulated_mass_change_glaciers_train.png", dpi=300) + if plot: + plt.show() + plt.close(fig) + + if len(mapsTrain) > 0: + df_gridded_annual1 = pd.read_csv(f"{pathFolder1}/gridded_annual_train.csv") + df_gridded_annual2 = pd.read_csv(f"{pathFolder2}/gridded_annual_train.csv") + + mapsFolder = f"{pathFolder}/maps" + os.makedirs(mapsFolder, exist_ok=True) + cfg = mbm.Config("11") # Fake cfg which is needed just for OGGM + assert set(mapsTrain).issubset(df_gridded_annual1.RGIId.unique()) + assert set(mapsTrain).issubset(df_gridded_annual2.RGIId.unique()) + for rgi_id in mapsTrain: + years = df_gridded_annual1[df_gridded_annual1.RGIId == rgi_id].YEAR.unique() + max1 = ( + df_gridded_annual1[df_gridded_annual1.RGIId == rgi_id].pred.abs().max() + ) + max2 = ( + df_gridded_annual2[df_gridded_annual2.RGIId == rgi_id].pred.abs().max() + ) + max_abs = max(max1, max2) + for year in years: + fig, axs = plt.subplots(1, 2, figsize=(12, 6)) + mbm.plots.mapGlacier( + df_gridded_annual1, + rgi_id, + year, + cfg, + ax=axs[0], + max_abs=max_abs, + title=name1, + ) + mbm.plots.mapGlacier( + df_gridded_annual2, + rgi_id, + year, + cfg, + ax=axs[1], + max_abs=max_abs, + title=name2, + ) + fig.suptitle(f"{rgi_id} year {year}") + plt.tight_layout() + fig.savefig(f"{mapsFolder}/{rgi_id}_{year}.pdf") + plt.close(fig) + del df_gridded_annual1, df_gridded_annual2 # Cumulated mass change on test data @@ -111,7 +180,7 @@ geoErr = df_geo1.set_index("RGIId").err.to_dict() # Plot cumulated mass change -fig = mbm.plots.cumulatedMassChange( +fig, l1 = mbm.plots.cumulatedMassChange( df_gridded_monthly1, geo={ rgi_id: {"mean": geoTarget[rgi_id], "err": geoErr[rgi_id]} @@ -120,12 +189,26 @@ ) del df_gridded_monthly1 df_gridded_monthly2 = pd.read_csv(f"{pathFolder2}/gridded_monthly_test.csv") -mbm.plots.cumulatedMassChange( +_, l2 = mbm.plots.cumulatedMassChange( df_gridded_monthly2, geo=None, axs=fig.axes, color_pred="orange", + titles={ + k: (f"{k} ({glacierNames[k]})" if glacierNames[k] is not None else None) + for k in glacierNames + }, +) +del df_gridded_monthly2 +fig.legend( + [l1, l2], + [name1, name2], + loc="lower center", + ncol=2, + fontsize=18, + bbox_to_anchor=(0.5, 0.02), ) +plt.tight_layout(rect=[0, 0.1, 1, 1]) fig.savefig(f"{pathFolder}/cumulated_mass_change_glaciers_test.pdf") @@ -133,3 +216,43 @@ if plot: plt.show() plt.close(fig) + +if maps: + df_gridded_annual1 = pd.read_csv(f"{pathFolder1}/gridded_annual_test.csv") + df_gridded_annual2 = pd.read_csv(f"{pathFolder2}/gridded_annual_test.csv") + + mapsFolder = f"{pathFolder}/maps" + os.makedirs(mapsFolder, exist_ok=True) + rgi_ids = df_gridded_annual1.RGIId.unique() + cfg = mbm.Config("11") # Fake cfg which is needed just for OGGM + assert set(rgi_ids) == set(df_gridded_annual2.RGIId.unique()) + for rgi_id in rgi_ids: + years = df_gridded_annual1[df_gridded_annual1.RGIId == rgi_id].YEAR.unique() + max1 = df_gridded_annual1[df_gridded_annual1.RGIId == rgi_id].pred.abs().max() + max2 = df_gridded_annual2[df_gridded_annual2.RGIId == rgi_id].pred.abs().max() + max_abs = max(max1, max2) + for year in years: + fig, axs = plt.subplots(1, 2, figsize=(12, 6)) + mbm.plots.mapGlacier( + df_gridded_annual1, + rgi_id, + year, + cfg, + ax=axs[0], + max_abs=max_abs, + title=name1, + ) + mbm.plots.mapGlacier( + df_gridded_annual2, + rgi_id, + year, + cfg, + ax=axs[1], + max_abs=max_abs, + title=name2, + ) + fig.suptitle(f"{rgi_id} year {year}") + plt.tight_layout() + fig.savefig(f"{mapsFolder}/{rgi_id}_{year}.pdf") + plt.close(fig) + del df_gridded_annual1, df_gridded_annual2 diff --git a/scripts/geo/eval.py b/scripts/geo/eval.py index a8e0a4a7..6c1fdb01 100644 --- a/scripts/geo/eval.py +++ b/scripts/geo/eval.py @@ -21,6 +21,7 @@ testData, setFeatures, ) +from scripts.common import default_glacier_name warnings.filterwarnings("ignore") @@ -62,6 +63,20 @@ action="store_true", help="Save predictions as CSV for further analysis or comparison.", ) +parser.add_argument( + "--maps", + dest="maps", + default=False, + action="store_true", + help="Generate annual MB maps for test glaciers.", +) +parser.add_argument( + "--mapsTrain", + dest="mapsTrain", + default=[], + nargs="+", + help="Generate annual MB maps for specific train glaciers.", +) args = parser.parse_args() modelFolder = args.modelFolder @@ -70,6 +85,8 @@ noTest = args.noTest onRegion = args.onRegion savePred = args.savePred +maps = args.maps +mapsTrain = args.mapsTrain pathFolder = os.path.join("logs", modelFolder) if not plot: @@ -167,7 +184,12 @@ data_test["y"] = test_set["y"] setFeatures(cfg, data_train, featuresInpModel) -df_X_train, y_train, df_X_val, y_val = trainValData(cfg, train_set, featuresInpModel) +df_X_train, y_train, df_X_val, y_val = trainValData( + cfg, + train_set, + featuresInpModel, + split_key=params["training"].get("splitVal", "group-meas-id"), +) df_X_test_subset = testData(cfg, test_set, featuresInpModel) @@ -184,7 +206,7 @@ # Load model and set to CPU -bestModelPath = mbm.training.loadBestModel(pathFolder, model) +bestModelPath, _ = mbm.training.loadBestModel(pathFolder, model) print(f"Loaded model {bestModelPath}") # loaded_model = mbm.models.CustomNeuralNetRegressor.load_model( # cfg, @@ -195,6 +217,7 @@ # model = model.to("cpu") +test_glacierNames = {} if len(df_X_test_subset) > 0 and not noTest: if sourceData == "switzerland": test_glaciers = list(data_test.GLACIER.unique()) @@ -202,6 +225,9 @@ test_glaciers = list(data_test.RGIId.unique()) elif "wgms" in sourceData: test_glaciers = list(data_test.RGIId.unique()) + test_glacierNames = mbm.data_processing.oggm_utils._glacier_name( + list(data_test.RGIId.unique()), cfg + ) assert set(df_X_test_subset.RGIId.unique()) == set(test_glaciers) @@ -249,6 +275,7 @@ scores_summer=scores_summer, ax_xlim=(-8, 6), ax_ylim=(-8, 6), + precLegend=2, ) fig.savefig(f"{pathFolder}/prediction_test_PMB.pdf") if plot: @@ -314,15 +341,46 @@ df_geo.to_csv(f"{pathFolder}/gridded_geodetic_test.csv") df_gridded_annual.to_csv(f"{pathFolder}/gridded_annual_test.csv") df_gridded_monthly.to_csv(f"{pathFolder}/gridded_monthly_test.csv") - del df_gridded_annual, df_gridded_monthly + if maps and len(test_glaciers) > 0: + mapsFolder = f"{pathFolder}/maps" + os.makedirs(mapsFolder, exist_ok=True) + # Initialize OGGM once for all to avoid repeated and useless computations + mbm.data_processing.oggm_utils._initialize_oggm_config("") + gdirs = mbm.data_processing.oggm_utils._initialize_glacier_directories( + test_glaciers, cfg + ) + for rgi_id, gdir in zip(test_glaciers, gdirs): + years = df_gridded_annual[df_gridded_annual.RGIId == rgi_id].YEAR.unique() + for year in years: + fig = mbm.plots.mapGlacier( + df_gridded_annual, rgi_id, year, cfg, gdir=gdir + ) + fig.savefig(f"{mapsFolder}/{rgi_id}_{year}.pdf") + plt.close(fig) + del df_gridded_annual, df_gridded_monthly + +else: + resTest = None if sourceData == "switzerland": - train_glaciers = list(data_train.GLACIER.unique()) + train_glaciers = list(df_X_train.GLACIER.unique()) + valid_glaciers = list(df_X_val.GLACIER.unique()) elif sourceData in ["iceland", "norway"]: - train_glaciers = list(data_train.RGIId.unique()) + train_glaciers = list(df_X_train.RGIId.unique()) + valid_glaciers = list(df_X_val.RGIId.unique()) elif "wgms" in sourceData: - train_glaciers = list(data_train.RGIId.unique()) + train_glaciers = list(df_X_train.RGIId.unique()) + valid_glaciers = list(df_X_val.RGIId.unique()) +train_glacierNames = mbm.data_processing.oggm_utils._glacier_name( + list(data_train.RGIId.unique()), cfg +) +glacierNames = train_glacierNames | test_glacierNames +for k in glacierNames: + if glacierNames[k] == "": + glacierNames[k] = default_glacier_name(k) +with open(os.path.join(pathFolder, "glacierNames.json"), "w") as f: + json.dump(glacierNames, f, indent=4, sort_keys=True) # Create dataloader train_gdl = mbm.dataloader.GeoDataLoader( @@ -330,12 +388,23 @@ train_glaciers, device=device, trainStakesDf=data_train, + glacierListVal=valid_glaciers, months_head_pad=months_head_pad, months_tail_pad=months_tail_pad, + valStakesDf=df_X_val, keyGlacierSel="GLACIER" if sourceData == "switzerland" else "RGIId", ) +with torch.no_grad(): + resVal = mbm.training.assessOnVal(model, train_gdl, params) + with open(os.path.join(pathFolder, "perf.json"), "w") as f: + json.dump({"test": resTest, "val": resVal}, f, indent=4) + +# PMB predictions grouped_ids_train = model.evaluate_group_pred(train_gdl) +grouped_ids_valid = model.evaluate_group_pred(train_gdl, val=True) + +# PMB train scores_train = mbm.metrics.seasonal_scores( grouped_ids_train, target_col="target", pred_col="pred" ) @@ -365,6 +434,7 @@ scores_summer=scores_summer, ax_xlim=(-14, 8), ax_ylim=(-14, 8), + precLegend=2, ) fig.savefig(f"{pathFolder}/prediction_train_PMB.pdf") if plot: @@ -404,8 +474,11 @@ scores[train_gl]["r2"]["s"] = scores_glacier["summer"]["r2"] scores[train_gl]["bias"]["s"] = scores_glacier["summer"]["bias"] +grouped_ids_train_valid = pd.concat( + [grouped_ids_train, grouped_ids_valid], ignore_index=True +) fig = mbm.plots.predVSTruthPerGlacier( - grouped_ids_train, + grouped_ids_train_valid, scores=scores, custom_order=train_gl_per_el, hue="PERIOD", @@ -416,6 +489,44 @@ plt.close(fig) +# PMB validation +scores_valid = mbm.metrics.seasonal_scores( + grouped_ids_valid, target_col="target", pred_col="pred" +) +scores_annual = { + "rmse": scores_valid["annual"]["rmse"], + "r2": scores_valid["annual"]["r2"], + "bias": scores_valid["annual"]["bias"], +} +scores_winter = { + "rmse": scores_valid["winter"]["rmse"], + "r2": scores_valid["winter"]["r2"], + "bias": scores_valid["winter"]["bias"], +} +if "summer" in scores_valid: + scores_summer = { + "rmse": scores_valid["summer"]["rmse"], + "r2": scores_valid["summer"]["r2"], + "bias": scores_valid["summer"]["bias"], + } +else: + scores_summer = None + +fig = mbm.plots.predVSTruthTimeSeries( + grouped_ids=grouped_ids_valid, + scores_annual=scores_annual, + scores_winter=scores_winter, + scores_summer=scores_summer, + ax_xlim=(-14, 8), + ax_ylim=(-14, 8), + precLegend=2, +) +fig.savefig(f"{pathFolder}/prediction_validation_PMB.pdf") +if plot: + plt.show() +plt.close(fig) + + geoPred, geoTarget, geoErr, dict_df_gridded = mbm.training.eval_geodetic( model, train_gdl, return_grid_pred=["annual", "monthly"] ) @@ -440,7 +551,13 @@ # Geodetic performance fig = mbm.plots.predVSTruthGlacierWide( - geoTarget, geoPred, geoErr, title="Glacier wide MB on train" + geoTarget, + geoPred, + geoErr, + title="Glacier wide MB on train", + ax_xlim=(-2.5, 2.0), + ax_ylim=(-2.5, 2.0), + legend=False, ) plt.savefig(os.path.join(pathFolder, "geodetic_train.png")) if plot: @@ -450,7 +567,12 @@ # Plot MB profile fig = mbm.plots.profilePerGlacier( - df_gridded_annual, custom_order=train_gl_per_el + df_gridded_annual, + custom_order=train_gl_per_el, + titles={ + k: (f"{k} ({glacierNames[k]})" if glacierNames[k] is not None else None) + for k in glacierNames + }, ) # , df_stakes=data_train) fig.savefig(f"{pathFolder}/PMB_profile_individual_glaciers_train.pdf") if plot: @@ -459,7 +581,7 @@ # Plot cumulated mass change -fig = mbm.plots.cumulatedMassChange( +fig, _ = mbm.plots.cumulatedMassChange( df_gridded_monthly, geo={ rgi_id: {"mean": geoTarget[rgi_id], "err": geoErr[rgi_id]} @@ -471,6 +593,22 @@ plt.show() plt.close(fig) +if len(mapsTrain) > 0: + mapsFolder = f"{pathFolder}/maps" + os.makedirs(mapsFolder, exist_ok=True) + # Initialize OGGM once for all to avoid repeated and useless computations + mbm.data_processing.oggm_utils._initialize_oggm_config("") + gdirs = mbm.data_processing.oggm_utils._initialize_glacier_directories( + mapsTrain, cfg + ) + for rgi_id, gdir in zip(mapsTrain, gdirs): + years = df_gridded_annual[df_gridded_annual.RGIId == rgi_id].YEAR.unique() + for year in years: + fig = mbm.plots.mapGlacier(df_gridded_annual, rgi_id, year, cfg, gdir=gdir) + fig.savefig(f"{mapsFolder}/{rgi_id}_{year}.pdf") + plt.close(fig) +del df_gridded_annual, df_gridded_monthly + # TODO: since we changed the iterator, is the evaluation consistent on train/test ? if onRegion: @@ -494,7 +632,11 @@ # Geodetic performance fig = mbm.plots.predVSTruthGlacierWide( - geoTarget, geoPred, geoErr, title="Glacier wide MB on the whole region" + geoTarget, + geoPred, + geoErr, + title="Glacier wide MB on the whole region", + legend=False, ) plt.savefig(os.path.join(pathFolder, "geodetic_region.png")) if plot: diff --git a/scripts/geo/train.py b/scripts/geo/train.py index 1a707705..95fe3869 100644 --- a/scripts/geo/train.py +++ b/scripts/geo/train.py @@ -11,6 +11,7 @@ from scripts.common import ( loadParams, + already_completed_trial, ) from scripts.nongeo.utils import getMetaData, setFeatures, trainValData, testData @@ -24,6 +25,13 @@ action="store_true", help="Force model to run on CPU, even if a GPU is available.", ) +parser.add_argument( + "-s", + "--suffix", + type=str, + default=None, + help="Suffix to add to the folder that contains the model once trained.", +) parser.add_argument( "--noTest", dest="noTest", @@ -51,20 +59,117 @@ default=None, help="Weight of the geodetic term.", ) +parser.add_argument( + "--gridsearch", + dest="gridsearch", + default=[], + nargs="+", + help="Grid search configuration file (name only) and grid search name.", +) args = parser.parse_args() params = loadParams(args.modelType) modelToLoad = args.load cpu = args.cpu +suffix = args.suffix noTest = args.noTest timeExec = args.time prof = args.prof wGeo = args.wGeo + +gridsearch = args.gridsearch +do_gridsearch = len(gridsearch) > 0 +if do_gridsearch: + assert len(gridsearch) == 2 + gridsearch_name = gridsearch[1] + gridsearch_config = gridsearch[0] + import yaml + import optuna + + with open("scripts/netcfg/" + gridsearch_config + ".yml") as stream: + try: + gridsearch_params = yaml.safe_load(stream) + except yaml.YAMLError as exc: + print(exc) + + def flatten_dict(d, parent_key="", sep="."): + result = {} + for key, value in d.items(): + new_key = f"{parent_key}{sep}{key}" if parent_key else key + if isinstance(value, dict): + result.update(flatten_dict(value, new_key, sep)) + else: + if isinstance(value, list): + if isinstance(value[0], (tuple, list)): + result[new_key] = tuple(tuple(v) for v in value) + else: + result[new_key] = tuple(value) + else: + result[new_key] = value + return result + + search_space = flatten_dict(gridsearch_params) + for k, v in search_space.items(): + if isinstance(v, (list, tuple)) and isinstance(v[0], (list, tuple)): + search_space[k] = tuple([",".join([str(e) for e in t]) for t in v]) + print(f"{search_space=}") + study = optuna.create_study( + study_name=gridsearch_name, + storage=optuna.storages.JournalStorage( + optuna.storages.journal.JournalFileBackend( + file_path="./journal_gridsearch.log" + ) + ), + sampler=optuna.samplers.GridSampler(search_space), + direction="minimize", + load_if_exists=True, + ) + trial = study.ask() + print(f"{trial.number=}") + params["training"]["log_prefix"] = gridsearch_name + f"_{trial.number}" + params["gridsearch"] = { + "study_name": trial.study.study_name, + "trial_number": trial.number, + "search_space": search_space, + } + + def recursive_update_from_flat(target, source, sep="."): + for flat_key, value in source.items(): + parts = flat_key.split(sep) + _set_recursive(target, parts, value) + + def _set_recursive(current, parts, value): + key = parts[0] + if len(parts) == 1: + current[key] = value + return + if key not in current or not isinstance(current[key], dict): + current[key] = {} + _set_recursive(current[key], parts[1:], value) + + candidate_params = { + k: trial.suggest_categorical(k, v) for k, v in search_space.items() + } + for k, v in candidate_params.items(): + if isinstance(v, str): + candidate_params[k] = [int(e) for e in v.split(",")] + print(f"{candidate_params=}") + exists, old_trial = already_completed_trial(study, candidate_params) + if exists: + print("This combination has already been tested.") + print("metric:", old_trial.value) + sys.exit(0) + recursive_update_from_flat(params, candidate_params) +else: + trial = None + if wGeo is not None: # Overwrite geodetic weight params["training"]["wGeo"] = wGeo wGeo = params["training"]["wGeo"] if params["training"]["log_suffix"] == "": params["training"]["log_suffix"] = f"wgeo={wGeo}" if wGeo > 0 else "" +if suffix is not None: + params["training"]["log_suffix"] += f"_{suffix}" featuresInpModel = params["model"]["inputs"] sourceData = params["training"]["source_data"] @@ -139,7 +244,12 @@ data_train["y"] = train_set["y"] setFeatures(cfg, data_train, featuresInpModel) -df_X_train, y_train, df_X_val, y_val = trainValData(cfg, train_set, featuresInpModel) +df_X_train, y_train, df_X_val, y_val = trainValData( + cfg, + train_set, + featuresInpModel, + split_key=params["training"].get("splitVal", "group-meas-id"), +) print( @@ -168,7 +278,7 @@ months_tail_pad=months_tail_pad, valStakesDf=df_X_val, keyGlacierSel="GLACIER" if sourceData == "switzerland" else "RGIId", - preloadGeodetic=True, + preloadGeodetic=wGeo > 0, ) @@ -178,7 +288,9 @@ model = model.to(device) if modelToLoad != "": - bestModelPath = mbm.training.loadBestModel(os.path.join("logs", modelToLoad), model) + bestModelPath, _ = mbm.training.loadBestModel( + os.path.join("logs", modelToLoad), model + ) print(f"Loaded model {bestModelPath}") optimType = params["training"]["optim"] @@ -226,7 +338,7 @@ months_head_pad=months_head_pad, months_tail_pad=months_tail_pad, keyGlacierSel="GLACIER" if sourceData == "switzerland" else "RGIId", - preloadGeodetic=True, + preloadGeodetic=wGeo > 0, ) else: gdl_test = None @@ -243,7 +355,7 @@ ) print() -bestModelPath = mbm.training.loadBestModel(ret["misc"]["log_dir"], model) +bestModelPath, bestVal = mbm.training.loadBestModel(ret["misc"]["log_dir"], model) print(f"Best model is {bestModelPath}") @@ -291,17 +403,33 @@ def default(self, obj): months_head_pad=months_head_pad, months_tail_pad=months_tail_pad, keyGlacierSel="GLACIER" if sourceData == "switzerland" else "RGIId", - preloadGeodetic=True, + preloadGeodetic=wGeo > 0, ) model.eval() with torch.no_grad(): print("Computing performance on test set") resTest = mbm.training.assessOnTest(ret["misc"]["log_dir"], model, gdl_test) + resVal = mbm.training.assessOnVal(model, gdl, params) + with open(os.path.join(ret["misc"]["log_dir"], "perf.json"), "w") as f: + json.dump({"test": resTest, "val": resVal}, f, indent=4) print("Performance:") print( json.dumps( - json.loads(json.dumps(resTest), parse_float=lambda x: round(float(x), 3)), + json.loads( + json.dumps({"test": resTest, "val": resVal}), + parse_float=lambda x: round(float(x), 3), + ), indent=2, ) ) + +if trial is not None: + trial.set_user_attr("log_dir", ret["misc"]["log_dir"]) + trial.set_user_attr("r2", resVal["r2"]) + trial.set_user_attr("bias", resVal["bias"]) + trial.set_user_attr("lossValStake", resVal["lossValStake"]) + trial.set_user_attr("lossValGeo", resVal["lossValGeo"]) + trial.set_user_attr("lossVal", resVal["lossVal"]) + trial.set_user_attr("rmse", resVal["rmse"]) + study.tell(trial, float(bestVal)) diff --git a/scripts/netcfg/stakes_wgms:11.yml b/scripts/netcfg/stakes_wgms:11.yml new file mode 100644 index 00000000..9f153fe2 --- /dev/null +++ b/scripts/netcfg/stakes_wgms:11.yml @@ -0,0 +1,100 @@ +model: + type: sequential + layers: [32, 16, 8] + inputs: + - ELEVATION_DIFFERENCE + - aspect + - fal + - slhf + - slope + - sshf + - ssrd + - str + - t2m + - tp + - svf +training: + source_data: wgms:11 + test_glaciers: + - RGI60-11.00597 + - RGI60-11.00918 + - RGI60-11.01238 + - RGI60-11.01776 + - RGI60-11.01946 + train_glaciers: + - RGI60-11.00002 + - RGI60-11.00006 + - RGI60-11.00116 + - RGI60-11.00135 + - RGI60-11.00190 + - RGI60-11.00251 + - RGI60-11.00267 + - RGI60-11.00289 + - RGI60-11.00300 + - RGI60-11.00603 + - RGI60-11.00638 + - RGI60-11.00647 + - RGI60-11.00719 + - RGI60-11.00746 + - RGI60-11.00781 + - RGI60-11.00787 + - RGI60-11.00804 + - RGI60-11.00817 + - RGI60-11.00819 + - RGI60-11.00833 + - RGI60-11.00843 + - RGI60-11.00878 + - RGI60-11.00892 + - RGI60-11.00897 + - RGI60-11.01198 + - RGI60-11.01280 + - RGI60-11.01367 + - RGI60-11.01376 + - RGI60-11.01450 + - RGI60-11.01509 + - RGI60-11.01662 + - RGI60-11.01704 + - RGI60-11.01806 + - RGI60-11.01834 + - RGI60-11.01876 + - RGI60-11.01930 + - RGI60-11.01962 + - RGI60-11.01987 + - RGI60-11.02024 + - RGI60-11.02072 + - RGI60-11.02244 + - RGI60-11.02245 + - RGI60-11.02249 + - RGI60-11.02282 + - RGI60-11.02285 + - RGI60-11.02299 + - RGI60-11.02448 + - RGI60-11.02540 + - RGI60-11.02600 + - RGI60-11.02624 + - RGI60-11.02648 + - RGI60-11.02674 + - RGI60-11.02679 + - RGI60-11.02692 + - RGI60-11.02704 + - RGI60-11.02740 + - RGI60-11.02745 + - RGI60-11.02746 + - RGI60-11.02764 + - RGI60-11.02766 + - RGI60-11.02770 + - RGI60-11.02773 + - RGI60-11.02774 + - RGI60-11.02775 + - RGI60-11.02787 + - RGI60-11.02801 + - RGI60-11.03135 + - RGI60-11.03166 + optim: ADAM + scheduler: StepLR + lr: 5e-3 + scheduler_step_size: 500 + Nepochs: 2000 + scalingStakes: glacier + bestModelCriterion: lossVal + wGeo: 0.0 diff --git a/scripts/nongeo/eval.py b/scripts/nongeo/eval.py index 23d9a1dc..391b9b67 100644 --- a/scripts/nongeo/eval.py +++ b/scripts/nongeo/eval.py @@ -108,7 +108,9 @@ data_train["y"] = train_set["y"] setFeatures(cfg, data_train, featuresInpModel) -df_X_train, y_train, df_X_val, y_val = trainValData(cfg, train_set, featuresInpModel) +df_X_train, y_train, df_X_val, y_val = trainValData( + cfg, train_set, featuresInpModel, split_key=params["training"]["splitVal"] +) df_X_test_subset = testData(cfg, test_set, featuresInpModel) diff --git a/scripts/nongeo/train.py b/scripts/nongeo/train.py index 25b7531b..f508b276 100644 --- a/scripts/nongeo/train.py +++ b/scripts/nongeo/train.py @@ -35,10 +35,11 @@ parser = argparse.ArgumentParser() parser.add_argument("modelType", type=str, help="Type of model to train.") parser.add_argument( - "--gpu", - type=bool, + "--cpu", + dest="cpu", default=False, - help="Train on GPU. By default training runs on CPU.", + action="store_true", + help="Force model to run on CPU, even if a GPU is available.", ) parser.add_argument( "-s", @@ -49,7 +50,7 @@ ) args = parser.parse_args() -runOnGpu = args.gpu +cpu = args.cpu suffix = args.suffix params = loadParams(args.modelType) featuresInpModel = params["model"]["inputs"] @@ -121,12 +122,14 @@ data_train["y"] = train_set["y"] setFeatures(cfg, data_train, featuresInpModel) -df_X_train, y_train, df_X_val, y_val = trainValData(cfg, train_set, featuresInpModel) +df_X_train, y_train, df_X_val, y_val = trainValData( + cfg, train_set, featuresInpModel, split_key=params["training"]["splitVal"] +) early_stop = EarlyStopping( monitor="valid_loss", - patience=20, + patience=50, threshold=1e-4, # Optional: stop only when improvement is very small ) lr_scheduler_cb = LRScheduler( @@ -134,7 +137,7 @@ monitor="valid_loss", mode="min", factor=0.5, - patience=5, + patience=20, threshold=0.01, threshold_mode="rel", verbose=True, @@ -147,7 +150,8 @@ def my_train_split(ds, y=None, **fit_params): return dataset, dataset_val -param_init = {"device": "cuda:0" if runOnGpu else "cpu"} +device = torch.device("cuda:0" if torch.cuda.is_available() and not cpu else "cpu") +param_init = {"device": device} model = mbm.models.buildModel(cfg, params=params) diff --git a/scripts/nongeo/utils.py b/scripts/nongeo/utils.py index fe64ff22..76e8e872 100644 --- a/scripts/nongeo/utils.py +++ b/scripts/nongeo/utils.py @@ -173,7 +173,7 @@ def forward(self, x): return self.model(x) -def trainValData(cfg, train_set, feature_columns): +def trainValData(cfg, train_set, feature_columns, split_key="group-meas-id"): """ Split training dataset into train and validation sets. @@ -187,7 +187,9 @@ def trainValData(cfg, train_set, feature_columns): data_train["y"] = train_set["y"] dataloader = mbm.dataloader.DataLoader(cfg, data=data_train) - train_itr, val_itr = dataloader.set_train_test_split(test_size=0.2) + train_itr, val_itr = dataloader.set_train_test_split( + test_size=0.2, type_fold=split_key + ) # Get all indices of the training and valing dataset at once from the iterators. Once called, the iterators are empty. train_indices, val_indices = list(train_itr), list(val_itr)