diff --git a/.gitignore b/.gitignore index cec17a77983bc017ef057b020600d15922e61f23..ff59ade5d38dac9c3cf2fecee6a676ee728a2162 100644 --- a/.gitignore +++ b/.gitignore @@ -58,4 +58,4 @@ htmlcov/ /test/test_modules/data/ report.html /TestExperiment/ -/testrun_network/ +/testrun_network*/ diff --git a/README.md b/README.md index 329b1f16ffc1f38ddc78408a5d04199b69bf7dc8..e49362e95e9a69c159a5b8d857ccb336cf58d3c6 100644 --- a/README.md +++ b/README.md @@ -6,4 +6,10 @@ This is a collection of all relevant functions used for ML stuff in the ESDE gro See a description [here](https://towardsdatascience.com/a-simple-guide-to-the-versions-of-the-inception-network-7fc52b863202) or take a look on the papers [Going Deeper with Convolutions (Szegedy et al., 2014)](https://arxiv.org/abs/1409.4842) -and [Network In Network (Lin et al., 2014)](https://arxiv.org/abs/1312.4400). \ No newline at end of file +and [Network In Network (Lin et al., 2014)](https://arxiv.org/abs/1312.4400). + + +# Installation + +* Install __proj__ on your machine using the console. E.g. for opensuse / leap `zypper install proj` +* c++ compiler required for cartopy installation \ No newline at end of file diff --git a/requirements.txt b/requirements.txt index e8ac7eb75d866b27d72af51fae31f417c471dd05..227da2b61976d6147902e05a54a3e414dc3f40cc 100644 --- a/requirements.txt +++ b/requirements.txt @@ -11,4 +11,16 @@ pytest-html pydot mock statsmodels +seaborn +dask==0.20.2 +toolz # for dask +cloudpickle # for dask +cython==0.29.14 +pyshp +six +pyproj +shapely +Cartopy==0.16.0 matplotlib +pillow +scipy \ No newline at end of file diff --git a/run.py b/run.py index 03eda04280a19e5c2bb9f1743c40f07e9e3fd2cc..8e4d9c46f4a39224335dd65b689f519943166d0f 100644 --- a/run.py +++ b/run.py @@ -2,15 +2,15 @@ __author__ = "Lukas Leufen" __date__ = '2019-11-14' -import logging import argparse +import logging from src.run_modules.experiment_setup import ExperimentSetup -from src.run_modules.run_environment import RunEnvironment -from src.run_modules.pre_processing import PreProcessing from src.run_modules.model_setup import ModelSetup -from src.run_modules.training import Training from src.run_modules.post_processing import PostProcessing +from src.run_modules.pre_processing import PreProcessing +from src.run_modules.run_environment import RunEnvironment +from src.run_modules.training import Training def main(parser_args): @@ -31,6 +31,7 @@ if __name__ == "__main__": formatter = '%(asctime)s - %(levelname)s: %(message)s [%(filename)s:%(funcName)s:%(lineno)s]' logging.basicConfig(format=formatter, level=logging.INFO) + # logging.basicConfig(format=formatter, level=logging.DEBUG) parser = argparse.ArgumentParser() parser.add_argument('--experiment_date', metavar='--exp_date', type=str, default=None, diff --git a/src/run_modules/modules.py b/run_hourly.py similarity index 51% rename from src/run_modules/modules.py rename to run_hourly.py index 5f0f12c19a87de7ba4ad1e0508906e75d1605563..af531aedbd275b133a087777334dba0ae24bd9c8 100644 --- a/src/run_modules/modules.py +++ b/run_hourly.py @@ -1,33 +1,41 @@ -import logging +__author__ = "Lukas Leufen" +__date__ = '2019-11-14' + + import argparse +import logging -from src.run_modules.run_environment import RunEnvironment from src.run_modules.experiment_setup import ExperimentSetup +from src.run_modules.model_setup import ModelSetup +from src.run_modules.post_processing import PostProcessing from src.run_modules.pre_processing import PreProcessing +from src.run_modules.run_environment import RunEnvironment +from src.run_modules.training import Training -class Training(RunEnvironment): +def main(parser_args): - def __init__(self): - super().__init__() + with RunEnvironment(): + ExperimentSetup(parser_args, stations=['DEBW107', 'DEBY081', 'DEBW013', 'DEBW076', 'DEBW087', 'DEBW001'], + station_type='background', trainable=True, sampling="hourly", window_history_size=48) + PreProcessing() + ModelSetup() -class PostProcessing(RunEnvironment): + Training() - def __init__(self): - super().__init__() + PostProcessing() if __name__ == "__main__": formatter = '%(asctime)s - %(levelname)s: %(message)s [%(filename)s:%(funcName)s:%(lineno)s]' - logging.basicConfig(format=formatter, level=logging.DEBUG) + logging.basicConfig(format=formatter, level=logging.INFO) + # logging.basicConfig(format=formatter, level=logging.DEBUG) parser = argparse.ArgumentParser() - parser.add_argument('--experiment_date', metavar='--exp_date', type=str, nargs=1, default=None, + parser.add_argument('--experiment_date', metavar='--exp_date', type=str, default=None, help="set experiment date as string") - parser_args = parser.parse_args() - with RunEnvironment(): - ExperimentSetup(parser_args, stations=['DEBW107', 'DEBY081', 'DEBW013', 'DEBW076', 'DEBW087', 'DEBW001'], - station_type='background') - PreProcessing() + args = parser.parse_args(["--experiment_date", "testrun"]) + + main(args) diff --git a/src/.gitignore b/src/.gitignore new file mode 100644 index 0000000000000000000000000000000000000000..8e2358dc56797578fe0de020aa827b1fef8663bf --- /dev/null +++ b/src/.gitignore @@ -0,0 +1 @@ +join_settings.py \ No newline at end of file diff --git a/src/data_handling/data_distributor.py b/src/data_handling/data_distributor.py index 74df5f6ac1c998e644fa7d89a688fc12dee21265..c6f38a6f0e70518956bcbbd51a6fdfc1a1e7849f 100644 --- a/src/data_handling/data_distributor.py +++ b/src/data_handling/data_distributor.py @@ -45,7 +45,7 @@ class Distributor(keras.utils.Sequence): for prev, curr in enumerate(range(1, num_mini_batches+1)): x = x_total[prev*self.batch_size:curr*self.batch_size, ...] y = [y_total[prev*self.batch_size:curr*self.batch_size, ...] for _ in range(mod_rank)] - if x is not None: + if x is not None: # pragma: no branch yield (x, y) if (k + 1) == len(self.generator) and curr == num_mini_batches and not fit_call: return diff --git a/src/data_handling/data_generator.py b/src/data_handling/data_generator.py index 1de0ab2092b46dfca6281963b260b1fc6bc65387..19a94fbb9dbbc8f382a225c852f34971a98395b8 100644 --- a/src/data_handling/data_generator.py +++ b/src/data_handling/data_generator.py @@ -1,12 +1,16 @@ __author__ = 'Felix Kleinert, Lukas Leufen' __date__ = '2019-11-07' +import os +from typing import Union, List, Tuple, Any + import keras +import xarray as xr +import pickle +import logging + from src import helpers from src.data_handling.data_preparation import DataPrep -import os -from typing import Union, List, Tuple -import xarray as xr class DataGenerator(keras.utils.Sequence): @@ -23,6 +27,9 @@ class DataGenerator(keras.utils.Sequence): interpolate_method: str = "linear", limit_nan_fill: int = 1, window_history_size: int = 7, window_lead_time: int = 4, transform_method: str = "standardise", **kwargs): self.data_path = os.path.abspath(data_path) + self.data_path_tmp = os.path.join(os.path.abspath(data_path), "tmp") + if not os.path.exists(self.data_path_tmp): + os.makedirs(self.data_path_tmp) self.network = network self.stations = helpers.to_list(stations) self.variables = variables @@ -70,11 +77,11 @@ class DataGenerator(keras.utils.Sequence): if self._iterator < self.__len__(): data = self.get_data_generator() self._iterator += 1 - if data.history is not None and data.label is not None: + if data.history is not None and data.label is not None: # pragma: no branch return data.history.transpose("datetime", "window", "Stations", "variables"), \ data.label.squeeze("Stations").transpose("datetime", "window") else: - self.__next__() + self.__next__() # pragma: no cover else: raise StopIteration @@ -85,24 +92,61 @@ class DataGenerator(keras.utils.Sequence): :return: The generator's time series of history data and its labels """ data = self.get_data_generator(key=item) - return data.history.transpose("datetime", "window", "Stations", "variables"), \ - data.label.squeeze("Stations").transpose("datetime", "window") + return data.get_transposed_history(), data.label.squeeze("Stations").transpose("datetime", "window") - def get_data_generator(self, key: Union[str, int] = None) -> DataPrep: + def get_data_generator(self, key: Union[str, int] = None, local_tmp_storage: bool = True) -> DataPrep: """ Select data for given key, create a DataPrep object and interpolate, transform, make history and labels and remove nans. :param key: station key to choose the data generator. + :param local_tmp_storage: say if data should be processed from scratch or loaded as already processed data from + tmp pickle file to save computational time (but of course more disk space required). :return: preprocessed data as a DataPrep instance """ station = self.get_station_key(key) - data = DataPrep(self.data_path, self.network, station, self.variables, station_type=self.station_type, - **self.kwargs) - data.interpolate(self.interpolate_dim, method=self.interpolate_method, limit=self.limit_nan_fill) - data.transform("datetime", method=self.transform_method) - data.make_history_window(self.interpolate_dim, self.window_history_size) - data.make_labels(self.target_dim, self.target_var, self.interpolate_dim, self.window_lead_time) - data.history_label_nan_remove(self.interpolate_dim) + try: + if not local_tmp_storage: + raise FileNotFoundError + data = self._load_pickle_data(station, self.variables) + except FileNotFoundError: + logging.info(f"load not pickle data for {station}") + data = DataPrep(self.data_path, self.network, station, self.variables, station_type=self.station_type, + **self.kwargs) + data.interpolate(self.interpolate_dim, method=self.interpolate_method, limit=self.limit_nan_fill) + data.transform("datetime", method=self.transform_method) + data.make_history_window(self.interpolate_dim, self.window_history_size) + data.make_labels(self.target_dim, self.target_var, self.interpolate_dim, self.window_lead_time) + data.history_label_nan_remove(self.interpolate_dim) + self._save_pickle_data(data) + return data + + def _save_pickle_data(self, data: Any): + """ + Save given data locally as .pickle in self.data_path_tmp with name '<station>_<var1>_<var2>_..._<varX>.pickle' + :param data: any data, that should be saved + """ + date = f"{self.kwargs.get('start')}_{self.kwargs.get('end')}" + vars = '_'.join(sorted(data.variables)) + station = ''.join(data.station) + file = os.path.join(self.data_path_tmp, f"{station}_{vars}_{date}_.pickle") + with open(file, "wb") as f: + pickle.dump(data, f) + logging.debug(f"save pickle data to {file}") + + def _load_pickle_data(self, station: Union[str, List[str]], variables: List[str]) -> Any: + """ + Load locally saved data from self.data_path_tmp and name '<station>_<var1>_<var2>_..._<varX>.pickle'. + :param station: station to load + :param variables: list of variables to load + :return: loaded data + """ + date = f"{self.kwargs.get('start')}_{self.kwargs.get('end')}" + vars = '_'.join(sorted(variables)) + station = ''.join(station) + file = os.path.join(self.data_path_tmp, f"{station}_{vars}_{date}_.pickle") + with open(file, "rb") as f: + data = pickle.load(f) + logging.debug(f"load pickle data from {file}") return data def get_station_key(self, key: Union[None, str, int, List[Union[None, str, int]]]) -> str: diff --git a/src/data_handling/data_preparation.py b/src/data_handling/data_preparation.py index db800f5e50146c4f73cde643f154bcb1a5047437..5bca71f52c9f136b5910d4e080491e0ff86484ae 100644 --- a/src/data_handling/data_preparation.py +++ b/src/data_handling/data_preparation.py @@ -1,17 +1,17 @@ __author__ = 'Felix Kleinert, Lukas Leufen' __date__ = '2019-10-16' - -import xarray as xr -import pandas as pd -import numpy as np +import datetime as dt import logging import os -from src import join, helpers -from src import statistics from typing import Union, List, Iterable -import datetime as dt +import numpy as np +import pandas as pd +import xarray as xr + +from src import join, helpers +from src import statistics # define a more general date type for type hinting date = Union[dt.date, dt.datetime] @@ -60,50 +60,56 @@ class DataPrep(object): self.meta = None self._transform_method = None self.statistics_per_var = kwargs.get("statistics_per_var", None) - if self.statistics_per_var is not None: + self.sampling = kwargs.get("sampling", "daily") + if self.statistics_per_var is not None or self.sampling == "hourly": self.load_data() else: - raise NotImplementedError # hourly data usage is not implemented yet - # self.data, self.meta = Fkf.read_hourly_data_from_csv_to_xarray(self.path, self.network, self.station, - # self.variables, **kwargs) + raise NotImplementedError("Either select hourly data or provide statistics_per_var.") def load_data(self): """ Load data and meta data either from local disk (preferred) or download new data from TOAR database if no local data is available. The latter case, store downloaded data locally if wished (default yes). """ - helpers.check_path_and_create(self.path) file_name = self._set_file_name() meta_file = self._set_meta_file_name() - try: - - logging.debug(f"try to load local data from: {file_name}") - data = self._slice_prep(xr.open_dataarray(file_name)) - self.data = self.check_for_negative_concentrations(data) - self.meta = pd.read_csv(meta_file, index_col=0) - if self.station_type is not None: - self.check_station_meta() - logging.debug("loading finished") - except FileNotFoundError as e: - logging.warning(e) - data, self.meta = self.download_data_from_join(file_name, meta_file) - data = self._slice_prep(data) - self.data = self.check_for_negative_concentrations(data) + if self.kwargs.get('overwrite_local_data', False): + logging.debug(f"overwrite_local_data is true, therefore reload {file_name} from JOIN") + if os.path.exists(file_name): + os.remove(file_name) + if os.path.exists(meta_file): + os.remove(meta_file) + self.download_data(file_name, meta_file) logging.debug("loaded new data from JOIN") + else: + try: + logging.debug(f"try to load local data from: {file_name}") + data = self._slice_prep(xr.open_dataarray(file_name)) + self.data = self.check_for_negative_concentrations(data) + self.meta = pd.read_csv(meta_file, index_col=0) + if self.station_type is not None: + self.check_station_meta() + logging.debug("loading finished") + except FileNotFoundError as e: + logging.warning(e) + self.download_data(file_name, meta_file) + logging.debug("loaded new data from JOIN") + + def download_data(self, file_name, meta_file): + data, self.meta = self.download_data_from_join(file_name, meta_file) + data = self._slice_prep(data) + self.data = self.check_for_negative_concentrations(data) def check_station_meta(self): """ Search for the entries in meta data and compare the value with the requested values. Raise a FileNotFoundError if the values mismatch. """ - check_dict = { - "station_type": self.station_type, - "network_name": self.network - } + check_dict = {"station_type": self.station_type, "network_name": self.network} for (k, v) in check_dict.items(): if self.meta.at[k, self.station[0]] != v: - logging.debug(f"meta data does not agree which given request for {k}: {v} (requested) != " + logging.debug(f"meta data does not agree with given request for {k}: {v} (requested) != " f"{self.meta.at[k, self.station[0]]} (local). Raise FileNotFoundError to trigger new " f"grapping from web.") raise FileNotFoundError @@ -117,7 +123,7 @@ class DataPrep(object): """ df_all = {} df, meta = join.download_join(station_name=self.station, stat_var=self.statistics_per_var, - station_type=self.station_type, network_name=self.network) + station_type=self.station_type, network_name=self.network, sampling=self.sampling) df_all[self.station[0]] = df # convert df_all to xarray xarr = {k: xr.DataArray(v, dims=['datetime', 'variables']) for k, v in df_all.items()} @@ -138,8 +144,8 @@ class DataPrep(object): return f"Dataprep(path='{self.path}', network='{self.network}', station={self.station}, " \ f"variables={self.variables}, station_type={self.station_type}, **{self.kwargs})" - def interpolate(self, dim: str, method: str = 'linear', limit: int = None, - use_coordinate: Union[bool, str] = True, **kwargs): + def interpolate(self, dim: str, method: str = 'linear', limit: int = None, use_coordinate: Union[bool, str] = True, + **kwargs): """ (Copy paste from dataarray.interpolate_na) Interpolate values according to different methods. @@ -193,6 +199,7 @@ class DataPrep(object): Perform inverse transformation :return: """ + def f_inverse(data, mean, std, method_inverse): if method_inverse == 'standardise': return statistics.standardise_inverse(data, mean, std), None, None @@ -319,8 +326,7 @@ class DataPrep(object): if (self.history is not None) and (self.label is not None): non_nan_history = self.history.dropna(dim=dim) non_nan_label = self.label.dropna(dim=dim) - intersect = np.intersect1d(non_nan_history.coords[dim].values, - non_nan_label.coords[dim].values) + intersect = np.intersect1d(non_nan_history.coords[dim].values, non_nan_label.coords[dim].values) if len(intersect) == 0: self.history = None @@ -380,8 +386,11 @@ class DataPrep(object): data.loc[..., used_chem_vars] = data.loc[..., used_chem_vars].clip(min=minimum) return data + def get_transposed_history(self): + if self.history is not None: + return self.history.transpose("datetime", "window", "Stations", "variables") -if __name__ == "__main__": +if __name__ == "__main__": dp = DataPrep('data/', 'dummy', 'DEBW107', ['o3', 'temp'], statistics_per_var={'o3': 'dma8eu', 'temp': 'maximum'}) print(dp) diff --git a/src/datastore.py b/src/datastore.py index 92623baf7c01f8199f653b7220db77a931986708..d9f844ff97acb3f5c6600205f91100219d9c53e6 100644 --- a/src/datastore.py +++ b/src/datastore.py @@ -2,8 +2,8 @@ __author__ = 'Lukas Leufen' __date__ = '2019-11-22' -from typing import Any, List, Tuple, Dict from abc import ABC +from typing import Any, List, Tuple, Dict class NameNotFoundInDataStore(Exception): @@ -144,6 +144,22 @@ class DataStoreByVariable(AbstractDataStore): """ return self._stride_through_scopes(name, scope)[2] + def get_default(self, name: str, scope: str, default: Any) -> Any: + """ + Same functionality like the standard get method. But this method adds a default argument that is returned if no + data was stored in the data store. Use this function with care, because it will not report any errors and just + return the given default value. Currently, there is no statement that reports, if the returned value comes from + the data store or the default value. + :param name: Name to look for + :param scope: scope to search the name for + :param default: default value that is return, if no data was found for given name and scope + :return: the stored object or the default value + """ + try: + return self._stride_through_scopes(name, scope)[2] + except (NameNotFoundInDataStore, NameNotFoundInScope): + return default + def _stride_through_scopes(self, name, scope, depth=0): if depth <= scope.count("."): local_scope = scope.rsplit(".", maxsplit=depth)[0] @@ -267,6 +283,22 @@ class DataStoreByScope(AbstractDataStore): """ return self._stride_through_scopes(name, scope)[2] + def get_default(self, name: str, scope: str, default: Any) -> Any: + """ + Same functionality like the standard get method. But this method adds a default argument that is returned if no + data was stored in the data store. Use this function with care, because it will not report any errors and just + return the given default value. Currently, there is no statement that reports, if the returned value comes from + the data store or the default value. + :param name: Name to look for + :param scope: scope to search the name for + :param default: default value that is return, if no data was found for given name and scope + :return: the stored object or the default value + """ + try: + return self._stride_through_scopes(name, scope)[2] + except (NameNotFoundInDataStore, NameNotFoundInScope): + return default + def _stride_through_scopes(self, name, scope, depth=0): if depth <= scope.count("."): local_scope = scope.rsplit(".", maxsplit=depth)[0] diff --git a/src/helpers.py b/src/helpers.py index 40a3f9762cd649651631e45d94b78c19562b9749..c33684508b7b36cbebc3cbf3ac826d1779f9df50 100644 --- a/src/helpers.py +++ b/src/helpers.py @@ -5,16 +5,17 @@ __date__ = '2019-10-21' import logging -import keras -import keras.backend as K import math -from typing import Union -import numpy as np import os import time import socket import datetime as dt +import keras.backend as K +import xarray as xr + +from typing import Dict, Callable + def to_list(arg): if not isinstance(arg, list): @@ -42,55 +43,6 @@ def l_p_loss(power: int): return loss -class LearningRateDecay(keras.callbacks.History): - """ - Decay learning rate during model training. Start with a base learning rate and lower this rate after every - n(=epochs_drop) epochs by drop value (0, 1], drop value = 1 means no decay in learning rate. - """ - - def __init__(self, base_lr: float = 0.01, drop: float = 0.96, epochs_drop: int = 8): - super().__init__() - self.lr = {'lr': []} - self.base_lr = self.check_param(base_lr, 'base_lr') - self.drop = self.check_param(drop, 'drop') - self.epochs_drop = self.check_param(epochs_drop, 'epochs_drop', upper=None) - - @staticmethod - def check_param(value: float, name: str, lower: Union[float, None] = 0, upper: Union[float, None] = 1): - """ - Check if given value is in interval. The left (lower) endpoint is open, right (upper) endpoint is closed. To - only one side of the interval, set the other endpoint to None. If both ends are set to None, just return the - value without any check. - :param value: value to check - :param name: name of the variable to display in error message - :param lower: left (lower) endpoint of interval, opened - :param upper: right (upper) endpoint of interval, closed - :return: unchanged value or raise ValueError - """ - if lower is None: - lower = -np.inf - if upper is None: - upper = np.inf - if lower < value <= upper: - return value - else: - raise ValueError(f"{name} is out of allowed range ({lower}, {upper}{')' if upper == np.inf else ']'}: " - f"{name}={value}") - - def on_epoch_begin(self, epoch: int, logs=None): - """ - Lower learning rate every epochs_drop epochs by factor drop. - :param epoch: current epoch - :param logs: ? - :return: update keras learning rate - """ - current_lr = self.base_lr * math.pow(self.drop, math.floor(epoch / self.epochs_drop)) - K.set_value(self.model.optimizer.lr, current_lr) - self.lr['lr'].append(current_lr) - logging.info(f"Set learning rate to {current_lr}") - return K.get_value(self.model.optimizer.lr) - - class TimeTracking(object): """ Track time to measure execution time. Time tracking automatically starts on initialisation and ends by calling stop @@ -136,25 +88,32 @@ class TimeTracking(object): def duration(self): return self._duration() + def __enter__(self): + return self + + def __exit__(self, exc_type, exc_val, exc_tb): + self.stop() + logging.info(f"undefined job finished after {self}") -def prepare_host(create_new=True): + +def prepare_host(create_new=True, sampling="daily"): hostname = socket.gethostname() try: user = os.getlogin() except OSError: user = "default" - if hostname == 'ZAM144': - path = f'/home/{user}/Data/toar_daily/' - elif hostname == 'zam347': - path = f'/home/{user}/Data/toar_daily/' - elif hostname == 'linux-gzsx': - path = f'/home/{user}/machinelearningtools/data/toar_daily/' - elif (len(hostname) > 2) and (hostname[:2] == 'jr'): - path = f'/p/project/cjjsc42/{user}/DATA/toar_daily/' - elif (len(hostname) > 2) and (hostname[:2] == 'jw'): - path = f'/p/home/jusers/{user}/juwels/intelliaq/DATA/toar_daily/' + if hostname == "ZAM144": + path = f"/home/{user}/Data/toar_{sampling}/" + elif hostname == "zam347": + path = f"/home/{user}/Data/toar_{sampling}/" + elif hostname == "linux-aa9b": + path = f"/home/{user}/machinelearningtools/data/toar_{sampling}/" + elif (len(hostname) > 2) and (hostname[:2] == "jr"): + path = f"/p/project/cjjsc42/{user}/DATA/toar_{sampling}/" + elif (len(hostname) > 2) and (hostname[:2] == "jw"): + path = f"/p/home/jusers/{user}/juwels/intelliaq/DATA/toar_{sampling}/" elif "runner-6HmDp9Qd-project-2411-concurrent" in hostname: - path = f'/home/{user}/machinelearningtools/data/toar_daily/' + path = f"/home/{user}/machinelearningtools/data/toar_{sampling}/" else: logging.error(f"unknown host '{hostname}'") raise OSError(f"unknown host '{hostname}'") @@ -173,12 +132,14 @@ def prepare_host(create_new=True): return path -def set_experiment_name(experiment_date=None, experiment_path=None): +def set_experiment_name(experiment_date=None, experiment_path=None, sampling=None): if experiment_date is None: experiment_name = "TestExperiment" else: experiment_name = f"{experiment_date}_network" + if sampling == "hourly": + experiment_name += f"_{sampling}" if experiment_path is None: experiment_path = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", experiment_name)) else: @@ -197,3 +158,35 @@ class PyTestRegex: def __repr__(self) -> str: return self._regex.pattern + + +def dict_to_xarray(d: Dict, coordinate_name: str) -> xr.DataArray: + """ + Convert a dictionary of 2D-xarrays to single 3D-xarray. The name of new coordinate axis follows <coordinate_name>. + :param d: dictionary with 2D-xarrays + :param coordinate_name: name of the new created axis (2D -> 3D) + :return: combined xarray + """ + xarray = None + for k, v in d.items(): + if xarray is None: + xarray = v + xarray.coords[coordinate_name] = k + else: + tmp_xarray = v + tmp_xarray.coords[coordinate_name] = k + xarray = xr.concat([xarray, tmp_xarray], coordinate_name) + return xarray + + +def float_round(number: float, decimals: int = 0, round_type: Callable = math.ceil) -> float: + """ + Perform given rounding operation on number with the precision of decimals. + :param number: the number to round + :param decimals: numbers of decimals of the rounding operations (default 0 -> round to next integer value) + :param round_type: the actual rounding operation. Can be any callable function like math.ceil, math.floor or python + built-in round operation. + :return: rounded number with desired precision + """ + multiplier = 10. ** decimals + return round_type(number * multiplier) / multiplier diff --git a/src/join.py b/src/join.py index 43271a7b0525b5d829ea761019176197c78c5468..351060f7bf4949801f94b04c13e3881f008389b6 100644 --- a/src/join.py +++ b/src/join.py @@ -2,14 +2,17 @@ __author__ = 'Felix Kleinert, Lukas Leufen' __date__ = '2019-10-16' -import requests -import logging -import pandas as pd import datetime as dt +import logging from typing import Iterator, Union, List, Dict + +import pandas as pd +import requests + from src import helpers +from src.join_settings import join_settings -join_url_base = 'https://join.fz-juelich.de/services/rest/surfacedata/' +# join_url_base = 'https://join.fz-juelich.de/services/rest/surfacedata/' str_or_none = Union[str, None] @@ -21,7 +24,7 @@ class EmptyQueryResult(Exception): def download_join(station_name: Union[str, List[str]], stat_var: dict, station_type: str = None, - network_name: str = None) -> [pd.DataFrame, pd.DataFrame]: + network_name: str = None, sampling: str = "daily") -> [pd.DataFrame, pd.DataFrame]: """ read data from JOIN/TOAR @@ -29,6 +32,7 @@ def download_join(station_name: Union[str, List[str]], stat_var: dict, station_t :param stat_var: key as variable like 'O3', values as statistics on keys like 'mean' :param station_type: set the station type like "traffic" or "background", can be none :param network_name: set the measurement network like "UBA" or "AIRBASE", can be none + :param sampling: sampling rate of the downloaded data, either set to daily or hourly (default daily) :returns: - df - data frame with all variables and statistics - meta - data frame with all meta information @@ -36,8 +40,15 @@ def download_join(station_name: Union[str, List[str]], stat_var: dict, station_t # make sure station_name parameter is a list station_name = helpers.to_list(station_name) + # get data connection settings + join_url_base, headers = join_settings(sampling) + # load series information - vars_dict = load_series_information(station_name, station_type, network_name) + vars_dict = load_series_information(station_name, station_type, network_name, join_url_base, headers) + + # correct stat_var values if data is not aggregated (hourly) + if sampling == "hourly": + [stat_var.update({k: "values"}) for k in stat_var.keys()] # download all variables with given statistic data = None @@ -49,10 +60,16 @@ def download_join(station_name: Union[str, List[str]], stat_var: dict, station_t # create data link opts = {'base': join_url_base, 'service': 'stats', 'id': vars_dict[var], 'statistics': stat_var[var], - 'sampling': 'daily', 'capture': 0, 'min_data_length': 1460} + 'sampling': sampling, 'capture': 0, 'min_data_length': 1460, 'format': 'json'} # load data - data = get_data(opts) + data = get_data(opts, headers) + + # adjust data format if given as list of list + # no branch cover because this just happens when downloading hourly data using a secret token, not available + # for CI testing. + if isinstance(data, list): # pragma: no branch + data = correct_data_format(data) # correct namespace of statistics stat = _correct_stat_name(stat_var[var]) @@ -70,30 +87,51 @@ def download_join(station_name: Union[str, List[str]], stat_var: dict, station_t raise EmptyQueryResult("No data found in JOIN.") -def get_data(opts: Dict) -> Union[Dict, List]: +def correct_data_format(data): + """ + Transform to the standard data format. For some cases (e.g. hourly data), the data is returned as list instead of + a dictionary with keys datetime, values and metadata. This functions addresses this issue and transforms the data + into the dictionary version. + :param data: data in hourly format + :return: the same data but formatted to fit with aggregated format + """ + formatted = {"datetime": [], + "values": [], + "metadata": data[-1]} + for d in data[:-1]: + for k, v in zip(["datetime", "values"], d): + formatted[k].append(v) + return formatted + + +def get_data(opts: Dict, headers: Dict) -> Union[Dict, List]: """ Download join data using requests framework. Data is returned as json like structure. Depending on the response structure, this can lead to a list or dictionary. :param opts: options to create the request url + :param headers: additional headers information like authorization, can be empty :return: requested data (either as list or dictionary) """ url = create_url(**opts) - response = requests.get(url) + response = requests.get(url, headers=headers) return response.json() -def load_series_information(station_name: List[str], station_type: str_or_none, network_name: str_or_none) -> Dict: +def load_series_information(station_name: List[str], station_type: str_or_none, network_name: str_or_none, + join_url_base: str, headers: Dict) -> Dict: """ List all series ids that are available for given station id and network name. :param station_name: Station name e.g. DEBW107 :param station_type: station type like "traffic" or "background" :param network_name: measurement network of the station like "UBA" or "AIRBASE" + :param join_url_base: base url name to download data from + :param headers: additional headers information like authorization, can be empty :return: all available series for requested station stored in an dictionary with parameter name (variable) as key and the series id as value. """ opts = {"base": join_url_base, "service": "series", "station_id": station_name[0], "station_type": station_type, "network_name": network_name} - station_vars = get_data(opts) + station_vars = get_data(opts, headers) vars_dict = {item[3].lower(): item[0] for item in station_vars} return vars_dict @@ -107,7 +145,11 @@ def _save_to_pandas(df: Union[pd.DataFrame, None], data: dict, stat: str, var: s :param var: variable the data is from (e.g. 'o3') :return: new created or concatenated data frame """ - index = map(lambda s: dt.datetime.strptime(s, "%Y-%m-%d %H:%M"), data['datetime']) + if len(data["datetime"][0]) == 19: + str_format = "%Y-%m-%d %H:%M:%S" + else: + str_format = "%Y-%m-%d %H:%M" + index = map(lambda s: dt.datetime.strptime(s, str_format), data['datetime']) if df is None: df = pd.DataFrame(data[stat], index=index, columns=[var]) else: @@ -151,8 +193,10 @@ def create_url(base: str, service: str, **kwargs: Union[str, int, float, None]) if __name__ == "__main__": + logging.basicConfig(level=logging.DEBUG) var_all_dic = {'o3': 'dma8eu', 'relhum': 'average_values', 'temp': 'maximum', 'u': 'average_values', 'v': 'average_values', 'no': 'dma8eu', 'no2': 'dma8eu', 'cloudcover': 'average_values', 'pblheight': 'maximum'} station = 'DEBW107' - download_join(station, var_all_dic) + # download_join(station, var_all_dic, sampling="daily") + download_join(station, var_all_dic, sampling="hourly") diff --git a/src/join_settings.py b/src/join_settings.py new file mode 100644 index 0000000000000000000000000000000000000000..365e8f39d25b28375eadf3b0dbda374feb5b158e --- /dev/null +++ b/src/join_settings.py @@ -0,0 +1,11 @@ + +def join_settings(sampling="daily"): + if sampling == "daily": # pragma: no branch + TOAR_SERVICE_URL = 'https://join.fz-juelich.de/services/rest/surfacedata/' + headers = {} + elif sampling == "hourly": + TOAR_SERVICE_URL = 'https://join.fz-juelich.de/services/rest/surfacedata/' + headers = {} + else: + raise NameError(f"Given sampling {sampling} is not supported, choose from either daily or hourly sampling.") + return TOAR_SERVICE_URL, headers diff --git a/src/model_modules/flatten.py b/src/model_modules/flatten.py index 707b0d4421fd506d16f93d7c56406abc878ab9f1..bbe92472ebb48e7486dede099dc098a161f51695 100644 --- a/src/model_modules/flatten.py +++ b/src/model_modules/flatten.py @@ -1,9 +1,10 @@ __author__ = "Felix Kleinert, Lukas Leufen" __date__ = '2019-12-02' -import keras from typing import Callable +import keras + def flatten_tail(input_X: keras.layers, name: str, bound_weight: bool = False, dropout_rate: float = 0.0, window_lead_time: int = 4, activation: Callable = keras.activations.relu, diff --git a/src/model_modules/inception_model.py b/src/model_modules/inception_model.py index 126dc15320f4b0ac7ff660b709601f083700d01f..11093df56a4b262f75eae7a6f7c05e6e44e1435d 100644 --- a/src/model_modules/inception_model.py +++ b/src/model_modules/inception_model.py @@ -1,9 +1,10 @@ __author__ = 'Felix Kleinert, Lukas Leufen' __date__ = '2019-10-22' +import logging + import keras import keras.layers as layers -import logging class InceptionModelBase: diff --git a/src/model_modules/keras_extensions.py b/src/model_modules/keras_extensions.py new file mode 100644 index 0000000000000000000000000000000000000000..e2a4b93219be2cbebfb35749560efa65c07226bb --- /dev/null +++ b/src/model_modules/keras_extensions.py @@ -0,0 +1,150 @@ +__author__ = 'Lukas Leufen, Felix Kleinert' +__date__ = '2020-01-31' + +import logging +import math +import pickle +from typing import Union + +import numpy as np +from keras import backend as K +from keras.callbacks import History, ModelCheckpoint + + +class HistoryAdvanced(History): + """ + This is almost an identical clone of the original History class. The only difference is that attributes epoch and + history are instantiated during the init phase and not during on_train_begin. This is required to resume an already + started but disrupted training from an saved state. This HistoryAdvanced callback needs to be added separately as + additional callback. To get the full history use this object for further steps instead of the default return of + training methods like fit_generator(). + + hist = HistoryAdvanced() + history = model.fit_generator(generator=.... , callbacks=[hist]) + history = hist + + If training was started from beginning this class is identical to the returned history class object. + """ + + def __init__(self): + self.epoch = [] + self.history = {} + super().__init__() + + def on_train_begin(self, logs=None): + pass + + +class LearningRateDecay(History): + """ + Decay learning rate during model training. Start with a base learning rate and lower this rate after every + n(=epochs_drop) epochs by drop value (0, 1], drop value = 1 means no decay in learning rate. + """ + + def __init__(self, base_lr: float = 0.01, drop: float = 0.96, epochs_drop: int = 8): + super().__init__() + self.lr = {'lr': []} + self.base_lr = self.check_param(base_lr, 'base_lr') + self.drop = self.check_param(drop, 'drop') + self.epochs_drop = self.check_param(epochs_drop, 'epochs_drop', upper=None) + self.epoch = [] + self.history = {} + + @staticmethod + def check_param(value: float, name: str, lower: Union[float, None] = 0, upper: Union[float, None] = 1): + """ + Check if given value is in interval. The left (lower) endpoint is open, right (upper) endpoint is closed. To + only one side of the interval, set the other endpoint to None. If both ends are set to None, just return the + value without any check. + :param value: value to check + :param name: name of the variable to display in error message + :param lower: left (lower) endpoint of interval, opened + :param upper: right (upper) endpoint of interval, closed + :return: unchanged value or raise ValueError + """ + if lower is None: + lower = -np.inf + if upper is None: + upper = np.inf + if lower < value <= upper: + return value + else: + raise ValueError(f"{name} is out of allowed range ({lower}, {upper}{')' if upper == np.inf else ']'}: " + f"{name}={value}") + + def on_train_begin(self, logs=None): + pass + + def on_epoch_begin(self, epoch: int, logs=None): + """ + Lower learning rate every epochs_drop epochs by factor drop. + :param epoch: current epoch + :param logs: ? + :return: update keras learning rate + """ + current_lr = self.base_lr * math.pow(self.drop, math.floor(epoch / self.epochs_drop)) + K.set_value(self.model.optimizer.lr, current_lr) + self.lr['lr'].append(current_lr) + logging.info(f"Set learning rate to {current_lr}") + return K.get_value(self.model.optimizer.lr) + + +class ModelCheckpointAdvanced(ModelCheckpoint): + """ + Enhance the standard ModelCheckpoint class by additional saves of given callbacks. Specify this callbacks as follow: + + lr = CustomLearningRate() + hist = CustomHistory() + callbacks_name = "your_custom_path_%s.pickle" + callbacks = [{"callback": lr, "path": callbacks_name % "lr"}, + {"callback": hist, "path": callbacks_name % "hist"}] + ckpt_callbacks = ModelCheckpointAdvanced(filepath=.... , callbacks=callbacks) + + Add this ckpt_callbacks as all other additional callbacks to the callback list. IMPORTANT: Always add ckpt_callbacks + as last callback to properly update all tracked callbacks, e.g. + + fit_generator(.... , callbacks=[lr, hist, ckpt_callbacks]) + + """ + def __init__(self, *args, **kwargs): + self.callbacks = kwargs.pop("callbacks") + super().__init__(*args, **kwargs) + + def update_best(self, hist): + """ + Update internal best on resuming a training process. Otherwise best is set to +/- inf depending on the + performance metric and the first trained model (first of the resuming training process) will always saved as + best model because its performance will be better than infinity. To prevent this behaviour and compare the + performance with the best model performance, call this method before resuming the training process. + :param hist: The History object from the previous (interrupted) training. + """ + self.best = hist.history.get(self.monitor)[-1] + + def update_callbacks(self, callbacks): + """ + Update all stored callback objects. The argument callbacks needs to follow the same convention like described + in the class description (list of dictionaries). Must be run before resuming a training process. + """ + self.callbacks = callbacks + + def on_epoch_end(self, epoch, logs=None): + """ + Save model as usual (see ModelCheckpoint class), but also save additional callbacks. + """ + super().on_epoch_end(epoch, logs) + + for callback in self.callbacks: + file_path = callback["path"] + if self.epochs_since_last_save == 0 and epoch != 0: + if self.save_best_only: + current = logs.get(self.monitor) + if current == self.best: + if self.verbose > 0: + print('\nEpoch %05d: save to %s' % (epoch + 1, file_path)) + with open(file_path, "wb") as f: + pickle.dump(callback["callback"], f) + else: + with open(file_path, "wb") as f: + if self.verbose > 0: + print('\nEpoch %05d: save to %s' % (epoch + 1, file_path)) + pickle.dump(callback["callback"], f) diff --git a/src/model_modules/linear_model.py b/src/model_modules/linear_model.py index 17a9b2326ab6ba1829ee4f65f0161de887e70778..3d5323e1b0303b497c1f26c4e84ee9b968380425 100644 --- a/src/model_modules/linear_model.py +++ b/src/model_modules/linear_model.py @@ -2,8 +2,8 @@ __author__ = "Felix Kleinert, Lukas Leufen" __date__ = '2019-12-11' -import statsmodels.api as sm import numpy as np +import statsmodels.api as sm class OrdinaryLeastSquaredModel: @@ -32,7 +32,7 @@ class OrdinaryLeastSquaredModel: def predict(self, data): data = sm.add_constant(self.reshape_xarray_to_numpy(data)) - return self.model.predict(data) + return np.atleast_2d(self.model.predict(data)) @staticmethod def reshape_xarray_to_numpy(data): diff --git a/src/model_modules/model_class.py b/src/model_modules/model_class.py index 1a8f7c4c400eaf75bdd1dc6af2e0993f662eac49..bb4801acfb5c3a643aecbcfad9cfdb758258d0ef 100644 --- a/src/model_modules/model_class.py +++ b/src/model_modules/model_class.py @@ -1,3 +1,5 @@ +import src.model_modules.keras_extensions + __author__ = "Lukas Leufen" __date__ = '2019-12-12' @@ -7,8 +9,6 @@ from typing import Any, Callable import keras -from src import helpers - class AbstractModelClass(ABC): @@ -109,8 +109,8 @@ class MyLittleModel(AbstractModelClass): self.regularizer = keras.regularizers.l2(0.1) self.initial_lr = 1e-2 self.optimizer = keras.optimizers.SGD(lr=self.initial_lr, momentum=0.9) - self.lr_decay = helpers.LearningRateDecay(base_lr=self.initial_lr, drop=.94, epochs_drop=10) - self.epochs = 2 + self.lr_decay = src.model_modules.keras_extensions.LearningRateDecay(base_lr=self.initial_lr, drop=.94, epochs_drop=10) + self.epochs = 20 self.batch_size = int(256) self.activation = keras.layers.PReLU @@ -154,3 +154,88 @@ class MyLittleModel(AbstractModelClass): """ self.loss = keras.losses.mean_squared_error + + +class MyBranchedModel(AbstractModelClass): + + """ + A customised model + + + with a 1x1 Conv, and 4 Dense layers (64, 32, 16, window_lead_time), where the last layer is the + output layer depending on the window_lead_time parameter. Dropout is used between the Convolution and the first + Dense layer. + """ + + def __init__(self, window_history_size, window_lead_time, channels): + + """ + Sets model and loss depending on the given arguments. + :param activation: activation function + :param window_history_size: number of historical time steps included in the input data + :param channels: number of variables used in input data + :param regularizer: <not used here> + :param dropout_rate: dropout rate used in the model [0, 1) + :param window_lead_time: number of time steps to forecast in the output layer + """ + + super().__init__() + + # settings + self.window_history_size = window_history_size + self.window_lead_time = window_lead_time + self.channels = channels + self.dropout_rate = 0.1 + self.regularizer = keras.regularizers.l2(0.1) + self.initial_lr = 1e-2 + self.optimizer = keras.optimizers.SGD(lr=self.initial_lr, momentum=0.9) + self.lr_decay = src.model_modules.keras_extensions.LearningRateDecay(base_lr=self.initial_lr, drop=.94, epochs_drop=10) + self.epochs = 20 + self.batch_size = int(256) + self.activation = keras.layers.PReLU + + # apply to model + self.set_model() + self.set_loss() + + def set_model(self): + + """ + Build the model. + :param activation: activation function + :param window_history_size: number of historical time steps included in the input data + :param channels: number of variables used in input data + :param dropout_rate: dropout rate used in the model [0, 1) + :param window_lead_time: number of time steps to forecast in the output layer + :return: built keras model + """ + + # add 1 to window_size to include current time step t0 + x_input = keras.layers.Input(shape=(self.window_history_size + 1, 1, self.channels)) + x_in = keras.layers.Conv2D(32, (1, 1), padding='same', name='{}_Conv_1x1'.format("major"))(x_input) + x_in = self.activation(name='{}_conv_act'.format("major"))(x_in) + x_in = keras.layers.Flatten(name='{}'.format("major"))(x_in) + x_in = keras.layers.Dropout(self.dropout_rate, name='{}_Dropout_1'.format("major"))(x_in) + x_in = keras.layers.Dense(64, name='{}_Dense_64'.format("major"))(x_in) + x_in = self.activation()(x_in) + out_minor_1 = keras.layers.Dense(self.window_lead_time, name='{}_Dense'.format("minor_1"))(x_in) + out_minor_1 = self.activation(name="minor_1")(out_minor_1) + x_in = keras.layers.Dense(32, name='{}_Dense_32'.format("major"))(x_in) + x_in = self.activation()(x_in) + out_minor_2 = keras.layers.Dense(self.window_lead_time, name='{}_Dense'.format("minor_2"))(x_in) + out_minor_2 = self.activation(name="minor_2")(out_minor_2) + x_in = keras.layers.Dense(16, name='{}_Dense_16'.format("major"))(x_in) + x_in = self.activation()(x_in) + x_in = keras.layers.Dense(self.window_lead_time, name='{}_Dense'.format("major"))(x_in) + out_main = self.activation(name="main")(x_in) + self.model = keras.Model(inputs=x_input, outputs=[out_minor_1, out_minor_2, out_main]) + + def set_loss(self): + + """ + Set the loss + :return: loss function + """ + + self.loss = [keras.losses.mean_absolute_error] + [keras.losses.mean_squared_error] + \ + [keras.losses.mean_squared_error] diff --git a/src/plotting/postprocessing_plotting.py b/src/plotting/postprocessing_plotting.py new file mode 100644 index 0000000000000000000000000000000000000000..01a565e3db284e56bf0b8c94420b71268fd21a80 --- /dev/null +++ b/src/plotting/postprocessing_plotting.py @@ -0,0 +1,596 @@ +__author__ = "Lukas Leufen, Felix Kleinert" +__date__ = '2019-12-17' + +import logging +import math +import os +import warnings +from typing import Dict, List, Tuple + +import cartopy.crs as ccrs +import cartopy.feature as cfeature +import matplotlib +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +import seaborn as sns +import xarray as xr +from matplotlib.backends.backend_pdf import PdfPages + +from src import helpers +from src.helpers import TimeTracking +from src.run_modules.run_environment import RunEnvironment + +logging.getLogger('matplotlib').setLevel(logging.WARNING) + + +class PlotMonthlySummary(RunEnvironment): + """ + Show a monthly summary over all stations for each lead time ("ahead") as box and whiskers plot. The plot is saved + in data_path with name monthly_summary_box_plot.pdf and 500dpi resolution. + """ + def __init__(self, stations: List, data_path: str, name: str, target_var: str, window_lead_time: int = None, + plot_folder: str = "."): + """ + Sets attributes and create plot + :param stations: all stations to plot + :param data_path: path, where the data is located + :param name: full name of the local files with a % as placeholder for the station name + :param target_var: display name of the target variable on plot's axis + :param window_lead_time: lead time to plot, if window_lead_time is higher than the available lead time or not given + the maximum lead time from data is used. (default None -> use maximum lead time from data). + :param plot_folder: path to save the plot (default: current directory) + """ + super().__init__() + self._data_path = data_path + self._data_name = name + self._data = self._prepare_data(stations) + self._window_lead_time = self._get_window_lead_time(window_lead_time) + self._plot(target_var, plot_folder) + + def _prepare_data(self, stations: List) -> xr.DataArray: + """ + Pre-process data required to plot. For each station, load locally saved predictions, extract the CNN and orig + prediction and group them into monthly bins (no aggregation, only sorting them). + :param stations: all stations to plot + :return: The entire data set, flagged with the corresponding month. + """ + forecasts = None + for station in stations: + logging.debug(f"... preprocess station {station}") + file_name = os.path.join(self._data_path, self._data_name % station) + data = xr.open_dataarray(file_name) + + data_cnn = data.sel(type="CNN").squeeze() + if len(data_cnn.shape) > 1: + data_cnn.coords["ahead"].values = [f"{days}d" for days in data_cnn.coords["ahead"].values] + + data_orig = data.sel(type="orig", ahead=1).squeeze() + data_orig.coords["ahead"] = "orig" + + data_concat = xr.concat([data_orig, data_cnn], dim="ahead") + data_concat = data_concat.drop("type") + + data_concat.index.values = data_concat.index.values.astype("datetime64[M]").astype(int) % 12 + 1 + data_concat = data_concat.clip(min=0) + + forecasts = xr.concat([forecasts, data_concat], 'index') if forecasts is not None else data_concat + return forecasts + + def _get_window_lead_time(self, window_lead_time: int): + """ + Extract the lead time from data and arguments. If window_lead_time is not given, extract this information from + data itself by the number of ahead dimensions. If given, check if data supports the give length. If the number + of ahead dimensions in data is lower than the given lead time, data's lead time is used. + :param window_lead_time: lead time from arguments to validate + :return: validated lead time, comes either from given argument or from data itself + """ + ahead_steps = len(self._data.ahead) + if window_lead_time is None: + window_lead_time = ahead_steps + return min(ahead_steps, window_lead_time) + + def _plot(self, target_var: str, plot_folder: str): + """ + Main plot function that creates a monthly grouped box plot over all stations but with separate boxes for each + lead time step. + :param target_var: display name of the target variable on plot's axis + :param plot_folder: path to save the plot + """ + data = self._data.to_dataset(name='values').to_dask_dataframe() + logging.debug("... start plotting") + color_palette = [matplotlib.colors.cnames["green"]] + sns.color_palette("Blues_d", self._window_lead_time).as_hex() + ax = sns.boxplot(x='index', y='values', hue='ahead', data=data.compute(), whis=1., palette=color_palette, + flierprops={'marker': '.', 'markersize': 1}, showmeans=True, + meanprops={'markersize': 1, 'markeredgecolor': 'k'}) + ax.set(xlabel='month', ylabel=f'{target_var}') + plt.tight_layout() + self._save(plot_folder) + + @staticmethod + def _save(plot_folder): + """ + Standard save method to store plot locally. The name of this plot is static. + :param plot_folder: path to save the plot + """ + plot_name = os.path.join(os.path.abspath(plot_folder), 'monthly_summary_box_plot.pdf') + logging.debug(f"... save plot to {plot_name}") + plt.savefig(plot_name, dpi=500) + plt.close('all') + + +class PlotStationMap(RunEnvironment): + """ + Plot geographical overview of all used stations as squares. Different data sets can be colorised by its key in the + input dictionary generators. The key represents the color to plot on the map. Currently, there is only a white + background, but this can be adjusted by loading locally stored topography data (not implemented yet). The plot is + saved under plot_path with the name station_map.pdf + """ + def __init__(self, generators: Dict, plot_folder: str = "."): + """ + Sets attributes and create plot + :param generators: dictionary with the plot color of each data set as key and the generator containing all stations + as value. + :param plot_folder: path to save the plot (default: current directory) + """ + super().__init__() + self._ax = None + self._plot(generators, plot_folder) + + def _draw_background(self): + """ + Draw coastline, lakes, ocean, rivers and country borders as background on the map. + """ + self._ax.add_feature(cfeature.COASTLINE.with_scale("50m"), edgecolor='black') + self._ax.add_feature(cfeature.LAKES.with_scale("50m")) + self._ax.add_feature(cfeature.OCEAN.with_scale("50m")) + self._ax.add_feature(cfeature.RIVERS.with_scale("50m")) + self._ax.add_feature(cfeature.BORDERS.with_scale("50m"), facecolor='none', edgecolor='black') + + def _plot_stations(self, generators): + """ + The actual plot function. Loops over all keys in generators dict and its containing stations and plots a square + and the stations's position on the map regarding the given color. + :param generators: dictionary with the plot color of each data set as key and the generator containing all + stations as value. + """ + if generators is not None: + for color, gen in generators.items(): + for k, v in enumerate(gen): + station_coords = gen.get_data_generator(k).meta.loc[['station_lon', 'station_lat']] + # station_names = gen.get_data_generator(k).meta.loc[['station_id']] + IDx, IDy = float(station_coords.loc['station_lon'].values), float( + station_coords.loc['station_lat'].values) + self._ax.plot(IDx, IDy, mfc=color, mec='k', marker='s', markersize=6, transform=ccrs.PlateCarree()) + + def _plot(self, generators: Dict, plot_folder: str): + """ + Main plot function to create the station map plot. Sets figure and calls all required sub-methods. + :param generators: dictionary with the plot color of each data set as key and the generator containing all + stations as value. + :param plot_folder: path to save the plot + """ + fig = plt.figure(figsize=(10, 5)) + self._ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree()) + self._ax.set_extent([0, 20, 42, 58], crs=ccrs.PlateCarree()) + self._draw_background() + self._plot_stations(generators) + self._save(plot_folder) + + @staticmethod + def _save(plot_folder): + """ + Standard save method to store plot locally. The name of this plot is static. + :param plot_folder: path to save the plot + """ + plot_name = os.path.join(os.path.abspath(plot_folder), 'station_map.pdf') + logging.debug(f"... save plot to {plot_name}") + plt.savefig(plot_name, dpi=500) + plt.close('all') + + +def plot_conditional_quantiles(stations: list, plot_folder: str = ".", rolling_window: int = 3, ref_name: str = 'orig', + pred_name: str = 'CNN', season: str = "", forecast_path: str = None, + plot_name_affix: str = "", units: str = "ppb"): + """ + This plot was originally taken from Murphy, Brown and Chen (1989): + https://journals.ametsoc.org/doi/pdf/10.1175/1520-0434%281989%29004%3C0485%3ADVOTF%3E2.0.CO%3B2 + + :param stations: stations to include in the plot (forecast data needs to be available already) + :param plot_folder: path to save the plot (default: current directory) + :param rolling_window: the rolling window mean will smooth the plot appearance (no smoothing in bin calculation, + this is only a cosmetic step, default: 3) + :param ref_name: name of the reference data series + :param pred_name: name of the investigated data series + :param season: season name to highlight if not empty + :param forecast_path: path to save the plot file + :param plot_name_affix: name to specify this plot (e.g. 'cali-ref', default: '') + :param units: units of the forecasted values (default: ppb) + """ + time = TimeTracking() + logging.debug(f"started plot_conditional_quantiles()") + # ignore warnings if nans appear in quantile grouping + warnings.filterwarnings("ignore", message="All-NaN slice encountered") + # ignore warnings if mean is calculated on nans + warnings.filterwarnings("ignore", message="Mean of empty slice") + # ignore warnings for y tick = 0 on log scale (instead of 0.00001 or similar) + warnings.filterwarnings("ignore", message="Attempted to set non-positive bottom ylim on a log-scaled axis.") + + def load_data(): + logging.debug("... load data") + data_collector = [] + for station in stations: + file = os.path.join(forecast_path, f"forecasts_{station}_test.nc") + data_tmp = xr.open_dataarray(file) + data_collector.append(data_tmp.loc[:, :, ['CNN', 'orig', 'OLS']].assign_coords(station=station)) + return xr.concat(data_collector, dim='station').transpose('index', 'type', 'ahead', 'station') + + def segment_data(data): + logging.debug("... segment data") + # combine index and station to multi index + data = data.stack(z=['index', 'station']) + # replace multi index by simple position index (order is not relevant anymore) + data.coords['z'] = range(len(data.coords['z'])) + # segment data of pred_name into bins + data.loc[pred_name, ...] = data.loc[pred_name, ...].to_pandas().T.apply(pd.cut, bins=bins, + labels=bins[1:]).T.values + return data + + def create_quantile_panel(data, q): + logging.debug("... create quantile panel") + # create empty xarray with dims: time steps ahead, quantiles, bin index (numbers create in previous step) + quantile_panel = xr.DataArray(np.full([data.ahead.shape[0], len(q), bins[1:].shape[0]], np.nan), + coords=[data.ahead, q, bins[1:]], dims=['ahead', 'quantiles', 'categories']) + # ensure that the coordinates are in the right order + quantile_panel = quantile_panel.transpose('ahead', 'quantiles', 'categories') + # calculate for each bin of the pred_name data the quantiles of the ref_name data + for bin in bins[1:]: + mask = (data.loc[pred_name, ...] == bin) + quantile_panel.loc[..., bin] = data.loc[ref_name, ...].where(mask).quantile(q, dim=['z']).T + + return quantile_panel + + def labels(plot_type, data_unit="ppb"): + names = (f"forecast concentration (in {data_unit})", f"observed concentration (in {data_unit})") + if plot_type == "orig": + return names + else: + return names[::-1] + + xlabel, ylabel = labels(ref_name, units) + + opts = {"q": [.1, .25, .5, .75, .9], "linetype": [':', '-.', '--', '-.', ':'], + "legend": ['.10th and .90th quantile', '.25th and .75th quantile', '.50th quantile', 'reference 1:1'], + "xlabel": xlabel, "ylabel": ylabel} + + # set name and path of the plot + base_name = "conditional_quantiles" + def add_affix(x): return f"_{x}" if len(x) > 0 else "" + plot_name = f"{base_name}{add_affix(season)}{add_affix(plot_name_affix)}_plot.pdf" + plot_path = os.path.join(os.path.abspath(plot_folder), plot_name) + + # check forecast path + if forecast_path is None: + raise ValueError("Forecast path is not given but required.") + + # load data and set data bins + orig_data = load_data() + bins = np.arange(0, math.ceil(orig_data.max().max()) + 1, 1).astype(int) + segmented_data = segment_data(orig_data) + quantile_panel = create_quantile_panel(segmented_data, q=opts["q"]) + + # init pdf output + pdf_pages = matplotlib.backends.backend_pdf.PdfPages(plot_path) + logging.debug(f"... plot path is {plot_path}") + + # create plot for each time step ahead + y2_max = 0 + for iteration, d in enumerate(segmented_data.ahead): + logging.debug(f"... plotting {d.values} time step(s) ahead") + # plot smoothed lines with rolling mean + smooth_data = quantile_panel.loc[d, ...].rolling(categories=rolling_window, center=True).mean().to_pandas().T + ax = smooth_data.plot(style=opts["linetype"], color='black', legend=False) + ax2 = ax.twinx() + # add reference line + ax.plot([0, bins.max()], [0, bins.max()], color='k', label='reference 1:1', linewidth=.8) + # add histogram of the segmented data (pred_name) + handles, labels = ax.get_legend_handles_labels() + segmented_data.loc[pred_name, d, :].to_pandas().hist(bins=bins, ax=ax2, color='k', alpha=.3, grid=False, + rwidth=1) + # add legend + plt.legend(handles[:3] + [handles[-1]], opts["legend"], loc='upper left', fontsize='large') + # adjust limits and set labels + ax.set(xlim=(0, bins.max()), ylim=(0, bins.max())) + ax.set_xlabel(opts["xlabel"], fontsize='x-large') + ax.tick_params(axis='x', which='major', labelsize=15) + ax.set_ylabel(opts["ylabel"], fontsize='x-large') + ax.tick_params(axis='y', which='major', labelsize=15) + ax2.yaxis.label.set_color('gray') + ax2.tick_params(axis='y', colors='gray') + ax2.yaxis.labelpad = -15 + ax2.set_yscale('log') + if iteration == 0: + y2_max = ax2.get_ylim()[1] + 100 + ax2.set(ylim=(0, y2_max * 10 ** 8), yticks=np.logspace(0, 4, 5)) + ax2.set_ylabel(' sample size', fontsize='x-large') + ax2.tick_params(axis='y', which='major', labelsize=15) + # set title and save current figure + title = f"{d.values} time step(s) ahead{f' ({season})' if len(season) > 0 else ''}" + plt.title(title) + pdf_pages.savefig() + # close all open figures / plots + pdf_pages.close() + plt.close('all') + logging.info(f"plot_conditional_quantiles() finished after {time}") + + +class PlotClimatologicalSkillScore(RunEnvironment): + """ + Create plot of climatological skill score after Murphy (1988) as box plot over all stations. A forecast time step + (called "ahead") is separately shown to highlight the differences for each prediction time step. Either each single + term is plotted (score_only=False) or only the resulting scores CASE I to IV are displayed (score_only=True, + default). Y-axis is adjusted following the data and not hard coded. The plot is saved under plot_folder path with + name skill_score_clim_{extra_name_tag}{model_setup}.pdf and resolution of 500dpi. + """ + def __init__(self, data: Dict, plot_folder: str = ".", score_only: bool = True, extra_name_tag: str = "", + model_setup: str = ""): + """ + Sets attributes and create plot + :param data: dictionary with station names as keys and 2D xarrays as values, consist on axis ahead and terms. + :param plot_folder: path to save the plot (default: current directory) + :param score_only: if true plot only scores of CASE I to IV, otherwise plot all single terms (default True) + :param extra_name_tag: additional tag that can be included in the plot name (default "") + :param model_setup: architecture type to specify plot name (default "CNN") + """ + super().__init__() + self._labels = None + self._data = self._prepare_data(data, score_only) + self._plot(plot_folder, score_only, extra_name_tag, model_setup) + + def _prepare_data(self, data: Dict, score_only: bool) -> pd.DataFrame: + """ + Shrink given data, if only scores are relevant. In any case, transform data to a plot friendly format. Also set + plot labels depending on the lead time dimensions. + :param data: dictionary with station names as keys and 2D xarrays as values + :param score_only: if true only scores of CASE I to IV are relevant + :return: pre-processed data set + """ + data = helpers.dict_to_xarray(data, "station") + self._labels = [str(i) + "d" for i in data.coords["ahead"].values] + if score_only: + data = data.loc[:, ["CASE I", "CASE II", "CASE III", "CASE IV"], :] + return data.to_dataframe("data").reset_index(level=[0, 1, 2]) + + def _label_add(self, score_only: bool): + """ + Adds the phrase "terms and " if score_only is disabled or empty string (if score_only=True). + :param score_only: if false all terms are relevant, otherwise only CASE I to IV + :return: additional label + """ + return "" if score_only else "terms and " + + def _plot(self, plot_folder, score_only, extra_name_tag, model_setup): + """ + Main plot function to plot climatological skill score. + :param plot_folder: path to save the plot + :param score_only: if true plot only scores of CASE I to IV, otherwise plot all single terms + :param extra_name_tag: additional tag that can be included in the plot name + :param model_setup: architecture type to specify plot name + """ + fig, ax = plt.subplots() + if not score_only: + fig.set_size_inches(11.7, 8.27) + sns.boxplot(x="terms", y="data", hue="ahead", data=self._data, ax=ax, whis=1., palette="Blues_d", + showmeans=True, meanprops={"markersize": 1, "markeredgecolor": "k"}, flierprops={"marker": "."}) + ax.axhline(y=0, color="grey", linewidth=.5) + ax.set(ylabel=f"{self._label_add(score_only)}skill score", xlabel="", title="summary of all stations") + handles, _ = ax.get_legend_handles_labels() + ax.legend(handles, self._labels) + plt.tight_layout() + self._save(plot_folder, extra_name_tag, model_setup) + + @staticmethod + def _save(plot_folder, extra_name_tag, model_setup): + """ + Standard save method to store plot locally. The name of this plot is dynamic. It includes the model setup like + 'CNN' and can additionally be adjusted using an extra name tag. + :param plot_folder: path to save the plot + :param extra_name_tag: additional tag that can be included in the plot name + :param model_setup: architecture type to specify plot name + """ + plot_name = os.path.join(plot_folder, f"skill_score_clim_{extra_name_tag}{model_setup}.pdf") + logging.debug(f"... save plot to {plot_name}") + plt.savefig(plot_name, dpi=500) + plt.close('all') + + +class PlotCompetitiveSkillScore(RunEnvironment): + """ + Create competitive skill score for the given model setup and the reference models ordinary least squared ("ols") and + the persistence forecast ("persi") for all lead times ("ahead"). The plot is saved under plot_folder with the name + skill_score_competitive_{model_setup}.pdf and resolution of 500dpi. + """ + def __init__(self, data: pd.DataFrame, plot_folder=".", model_setup="CNN"): + """ + :param data: data frame with index=['cnn-persi', 'ols-persi', 'cnn-ols'] and columns "ahead" containing the pre- + calculated comparisons for cnn, persistence and ols. + :param plot_folder: path to save the plot (default: current directory) + :param model_setup: architecture type (default "CNN") + """ + super().__init__() + self._labels = None + self._data = self._prepare_data(data) + self._plot(plot_folder, model_setup) + + def _prepare_data(self, data: pd.DataFrame) -> pd.DataFrame: + """ + Reformat given data and create plot labels. Introduces the dimensions stations and comparison + :param data: data frame with index=['cnn-persi', 'ols-persi', 'cnn-ols'] and columns "ahead" containing the pre- + calculated comparisons for cnn, persistence and ols. + :return: processed data + """ + data = pd.concat(data, axis=0) + data = xr.DataArray(data, dims=["stations", "ahead"]).unstack("stations") + data = data.rename({"stations_level_0": "stations", "stations_level_1": "comparison"}) + data = data.to_dataframe("data").unstack(level=1).swaplevel() + data.columns = data.columns.levels[1] + self._labels = [str(i) + "d" for i in data.index.levels[1].values] + return data.stack(level=0).reset_index(level=2, drop=True).reset_index(name="data") + + def _plot(self, plot_folder, model_setup): + """ + Main plot function to plot skill scores of the comparisons cnn-persi, ols-persi and cnn-ols. + :param plot_folder: path to save the plot + :param model_setup: + :return: architecture type to specify plot name + """ + fig, ax = plt.subplots() + sns.boxplot(x="comparison", y="data", hue="ahead", data=self._data, whis=1., ax=ax, palette="Blues_d", + showmeans=True, meanprops={"markersize": 3, "markeredgecolor": "k"}, flierprops={"marker": "."}, + order=["cnn-persi", "ols-persi", "cnn-ols"]) + ax.axhline(y=0, color="grey", linewidth=.5) + + ax.set(ylabel="skill score", xlabel="competing models", title="summary of all stations", ylim=self._ylim()) + handles, _ = ax.get_legend_handles_labels() + ax.legend(handles, self._labels) + plt.tight_layout() + self._save(plot_folder, model_setup) + + def _ylim(self) -> Tuple[float, float]: + """ + Calculate y-axis limits from data. Lower is the minimum of either 0 or data's minimum (reduced by small + subtrahend) and upper limit is data's maximum (increased by a small addend). + :return: + """ + lower = np.min([0, helpers.float_round(self._data.min()[2], 2) - 0.1]) + upper = helpers.float_round(self._data.max()[2], 2) + 0.1 + return lower, upper + + @staticmethod + def _save(plot_folder, model_setup): + """ + Standard save method to store plot locally. The name of this plot is dynamic by including the model setup. + :param plot_folder: path to save the plot + :param model_setup: architecture type to specify plot name + """ + plot_name = os.path.join(plot_folder, f"skill_score_competitive_{model_setup}.pdf") + logging.debug(f"... save plot to {plot_name}") + plt.savefig(plot_name, dpi=500) + plt.close() + + +class PlotTimeSeries(RunEnvironment): + + def __init__(self, stations: List, data_path: str, name: str, window_lead_time: int = None, plot_folder: str = ".", + sampling="daily"): + super().__init__() + self._data_path = data_path + self._data_name = name + self._stations = stations + self._window_lead_time = self._get_window_lead_time(window_lead_time) + self._sampling = self._get_sampling(sampling) + self._plot(plot_folder) + + @staticmethod + def _get_sampling(sampling): + if sampling == "daily": + return "D" + elif sampling == "hourly": + return "h" + + def _get_window_lead_time(self, window_lead_time: int): + """ + Extract the lead time from data and arguments. If window_lead_time is not given, extract this information from + data itself by the number of ahead dimensions. If given, check if data supports the give length. If the number + of ahead dimensions in data is lower than the given lead time, data's lead time is used. + :param window_lead_time: lead time from arguments to validate + :return: validated lead time, comes either from given argument or from data itself + """ + ahead_steps = len(self._load_data(self._stations[0]).ahead) + if window_lead_time is None: + window_lead_time = ahead_steps + return min(ahead_steps, window_lead_time) + + def _load_data(self, station): + logging.debug(f"... preprocess station {station}") + file_name = os.path.join(self._data_path, self._data_name % station) + data = xr.open_dataarray(file_name) + return data.sel(type=["CNN", "orig"]) + + def _plot(self, plot_folder): + pdf_pages = self._create_pdf_pages(plot_folder) + start, end = self._get_time_range(self._load_data(self._stations[0])) + for pos, station in enumerate(self._stations): + data = self._load_data(station) + fig, axes, factor = self._create_subplots(start, end) + nan_list = [] + for i_year in range(end - start + 1): + data_year = data.sel(index=f"{start + i_year}") + for i_half_of_year in range(factor): + pos = 2 * i_year + i_half_of_year + plot_data = self._create_plot_data(data_year, factor, i_half_of_year) + self._plot_orig(axes[pos], plot_data) + self._plot_ahead(axes[pos], plot_data) + if np.isnan(plot_data.values).all(): + nan_list.append(pos) + self._clean_up_axes(nan_list, axes, fig) + self._save_page(station, pdf_pages) + pdf_pages.close() + plt.close('all') + + @staticmethod + def _clean_up_axes(nan_list, axes, fig): + for i in reversed(nan_list): + fig.delaxes(axes[i]) + + @staticmethod + def _save_page(station, pdf_pages): + plt.suptitle(station) + plt.legend() + plt.tight_layout() + pdf_pages.savefig(dpi=500) + + @staticmethod + def _create_plot_data(data, factor, running_index): + if factor > 1: + if running_index == 0: + data = data.where(data["index.month"] < 7) + else: + data = data.where(data["index.month"] >= 7) + return data + + def _create_subplots(self, start, end): + factor = 1 + if self._sampling == "h": + factor = 2 + f, ax = plt.subplots((end - start + 1) * factor, sharey=True, figsize=(50, 30)) + return f, ax, factor + + def _plot_ahead(self, ax, data): + color = sns.color_palette("Blues_d", self._window_lead_time).as_hex() + for ahead in data.coords["ahead"].values: + plot_data = data.sel(type="CNN", ahead=ahead).drop(["type", "ahead"]).squeeze() + index = plot_data.index + np.timedelta64(int(ahead), self._sampling) + label = f"{ahead}{self._sampling}" + ax.plot(index, plot_data.values, color=color[ahead-1], label=label) + + def _plot_orig(self, ax, data): + orig_data = data.sel(type="orig", ahead=1) + index = data.index + np.timedelta64(1, self._sampling) + ax.plot(index, orig_data.values, color=matplotlib.colors.cnames["green"], label="orig") + + @staticmethod + def _get_time_range(data): + def f(x, f_x): + return pd.to_datetime(f_x(x.index.values)).year + return f(data, min), f(data, max) + + @staticmethod + def _create_pdf_pages(plot_folder): + """ + Standard save method to store plot locally. The name of this plot is static. + :param plot_folder: path to save the plot + """ + plot_name = os.path.join(os.path.abspath(plot_folder), 'timeseries_plot.pdf') + logging.debug(f"... save plot to {plot_name}") + return matplotlib.backends.backend_pdf.PdfPages(plot_name) diff --git a/src/plotting/training_monitoring.py b/src/plotting/training_monitoring.py index b18cce7a8993899621295644016f1e126d0dfac8..7e656895c5eecdabe1ef26869b68fb9494ed4c8c 100644 --- a/src/plotting/training_monitoring.py +++ b/src/plotting/training_monitoring.py @@ -5,12 +5,11 @@ __date__ = '2019-12-11' from typing import Union, Dict, List import keras -import pandas as pd import matplotlib import matplotlib.pyplot as plt +import pandas as pd -from src.helpers import LearningRateDecay - +from src.model_modules.keras_extensions import LearningRateDecay matplotlib.use('Agg') history_object = Union[Dict, keras.callbacks.History] @@ -19,46 +18,66 @@ lr_object = Union[Dict, LearningRateDecay] class PlotModelHistory: """ - Plots history of all losses for a training event. For default loss and val_loss are plotted. If further losses are - provided (name must somehow include the word `loss`), this additional information is added to the plot with an - separate y-axis scale on the right side (shared for all additional losses). The plot is saved locally. For a proper - saving behaviour, the parameter filename must include the absolute path for the plot. + Plots history of all plot_metrics (default: loss) for a training event. For default plot_metric and val_plot_metric + are plotted. If further metrics are provided (name must somehow include the word `<plot_metric>`), this additional + information is added to the plot with an separate y-axis scale on the right side (shared for all additional + metrics). The plot is saved locally. For a proper saving behaviour, the parameter filename must include the absolute + path for the plot. """ - def __init__(self, filename: str, history: history_object): + def __init__(self, filename: str, history: history_object, plot_metric: str = "loss", main_branch: bool = False): """ Sets attributes and create plot :param filename: saving name of the plot to create (preferably absolute path if possible), the filename needs a format ending like .pdf or .png to work. :param history: the history object (or a dict with at least 'loss' and 'val_loss' as keys) to plot loss from + :param plot_metric: the metric to plot (e.b. mean_squared_error, mse, mean_absolute_error, loss, default: loss) + :param main_branch: switch between only looking for metrics that go with 'main' or for all occurrences (default: + False -> look for losses from all branches, not only from main) """ if isinstance(history, keras.callbacks.History): history = history.history self._data = pd.DataFrame.from_dict(history) + self._plot_metric = self._get_plot_metric(history, plot_metric, main_branch) self._additional_columns = self._filter_columns(history) self._plot(filename) @staticmethod - def _filter_columns(history: Dict) -> List[str]: + def _get_plot_metric(history, plot_metric, main_branch): + if plot_metric.lower() == "mse": + plot_metric = "mean_squared_error" + elif plot_metric.lower() == "mae": + plot_metric = "mean_absolute_error" + available_keys = [k for k in history.keys() if plot_metric in k and ("main" in k.lower() if main_branch else True)] + available_keys.sort(key=len) + return available_keys[0] + + def _filter_columns(self, history: Dict) -> List[str]: """ - Select only columns named like %loss%. The default losses 'loss' and 'val_loss' are also removed. - :param history: a dict with at least 'loss' and 'val_loss' as keys (can be derived from keras History.history) - :return: filtered columns including all loss variations except loss and val_loss. + Select only columns named like %<plot_metric>%. The default metrics '<plot_metric>' and 'val_<plot_metric>' are + also removed. + :param history: a dict with at least '<plot_metric>' and 'val_<plot_metric>' as keys (can be derived from keras + History.history) + :return: filtered columns including all plot_metric variations except <plot_metric> and val_<plot_metric>. """ - cols = list(filter(lambda x: "loss" in x, history.keys())) - cols.remove("val_loss") - cols.remove("loss") + cols = list(filter(lambda x: self._plot_metric in x, history.keys())) + # heuristic: there is always val_<plot_metric> and <plot_metric> available in cols, because this is generated by + # the keras framework. If this metric isn't available the self._get_plot_metric() will fail before (but only + # because it is executed before) + cols.remove(f"val_{self._plot_metric}") + cols.remove(self._plot_metric) return cols def _plot(self, filename: str) -> None: """ - Actual plot routine. Plots loss and val_loss as default. If more losses are provided, they will be added with - an additional yaxis on the right side. The plot is saved in filename. + Actual plot routine. Plots <plot_metric> and val_<plot_metric> as default. If more plot_metrics are provided, + they will be added with an additional yaxis on the right side. The plot is saved in filename. :param filename: name (including total path) of the plot to save. """ - ax = self._data[["loss", "val_loss"]].plot(linewidth=0.7) + ax = self._data[[self._plot_metric, f"val_{self._plot_metric}"]].plot(linewidth=0.7) if len(self._additional_columns) > 0: self._data[self._additional_columns].plot(linewidth=0.7, secondary_y=True, ax=ax) - ax.set(xlabel="epoch", ylabel="loss", title=f"Model loss: best = {self._data[['val_loss']].min().values}") + title = f"Model {self._plot_metric}: best = {self._data[[f'val_{self._plot_metric}']].min().values}" + ax.set(xlabel="epoch", ylabel=self._plot_metric, title=title) ax.axhline(y=0, color="gray", linewidth=0.5) plt.tight_layout() plt.savefig(filename) diff --git a/src/run_modules/README.md b/src/run_modules/README.md new file mode 100644 index 0000000000000000000000000000000000000000..581811f1d35b3abd2c7d04a60d1027237bde0a91 --- /dev/null +++ b/src/run_modules/README.md @@ -0,0 +1,99 @@ +# On data handling + +This readme declares which function loads which data and where it is stored. + +## experiment setup + +**data_path** is the destination where all downloaded data is locally stored. Data is downloaded from TOARDB either using +the JOIN interface or a direct connection to the underlying PostgreSQL DB. If data was already downloaded, no new +download will be started. Missing data will be downloaded on the fly and saved in data_path. + +`data_path = src.helpers.prepare_host()` + + Current implementation leads to following paths: + + | hostname | path | comment | + | --- | --- | --- | + | ZAM144 | `/home/{user}/Data/toar_daily/` | notebook Felix | + | zam347 | `/home/{user}/Data/toar_daily/` | ESDE server | + | linux-gzsx | `/home/{user}/machinelearningtools/data/toar_daily/` | notebook Lukas | + | jureca | `/p/project/cjjsc42/{user}/DATA/toar_daily/` | JURECA | + | juwels | `/p/home/jusers/{user}/juwels/intelliaq/DATA/toar_daily/` | JUWELS | + | runner-6HmDp9Qd-project-2411-concurrent | `/home/{user}/machinelearningtools/data/toar_daily/` | gitlab-runner | + +**experiment_path** is the root folder in that all results from the experiment are saved. For each experiment there should +be distinct folder. Experiment path is can be set in ExperimentSetup. `experiment_date` can be set by parser_args and +`experiment_path` (this argument is not the same as the internal stored experiment_path!) as args. The *experiment_path* +is the combination of both given arguments `os.path.join(experiment_path, f"{experiment_date}_network")`. Inside this +folder, several subfolders are created in the course of the program. + +```bash +data_path + <station1>_<var1>_<var2>_..._<varx>.nc + <station1>_<var1>_<var2>_..._<varx>_meta.csv + <station2>_<var1>_<var2>_..._<varx>.nc + <station2>_<var1>_<var2>_..._<varx>_meta.csv +------ +experiment_path +| history.json +| history_lr.json +| <experiment_name>_model.pdf +| <experiment_name>_model-best.h5 +| <experiment_name>_my_model.h5 +├─── forecasts +| forecasts_<station1>_test.nc +| forecasts_<station2>_test.nc +| ... +└─── plots + conditional_quantiles_cali-ref_plot.pdf + conditional_quantiles_like-bas_plot.pdf + monthly_summary_box_plot.pdf + skill_score_clim_all_terms_<architecture>.pdf + skill_score_clim_<architecture>.pdf + skill_score_competitive_<architecture>.pdf + station_map.pdf + <experiment_name>_history_learning_rate.pdf + <experiment_name>_history_loss.pdf + <experiment_name>_history_main_loss.pdf + <experiment_name>_history_main_mse.pdf + ... + +``` + +**plot_path** includes all created plots. If not given, this is create into the experiment_path by default (as shown in +the folder structure above). Can be customised by `ExperimentSetup(plot_path=<path>)`. + +**forecast_path** is the place, where all forecasts are stored as netcdf file. Each file consists exactly one single +station. If not given, this is create into the experiment_path by default (as shown in the folder structure above). Can +be customised by `ExperimentSetup(forecast_path=<path>)`. + +## pre-processing + +Each requested station is check whether it is already included in *data_path*. The files all following the naming +convention `<station_name>_<sorted_list_of_all_variables_split_by_underscore>.nc`. E.g. the station *DEBW013* with the +variables cloudcover, NO, NO2, O3 and temp (all TOARDB short names) is saved as `DEBW013_cloudcover_no_no2_o3_temp.nc`, +whereas the same station with only O3 and temperature becomes `DEBW013_o3_temp.nc`. Although all data of the latter file +is potentially also included in the former file, the program will always download the data specification for new and +save this data into a new file. Only if the exactly fitting file is available locally, no data is downloaded. **NOTE**: +There is no check on data time range, only the name is compared. Set `overwrite_local_data=True` +in `experiment_setup.py` to overwrite local data by downloading new data. + +## model setup + +**checkpoint** is created inside *experiment_path* as `<experiment_name>_model-best.h5`. + +The architecture of the model is plotted into *experiment_path* as `<experiment_name>_model.pdf` + +## training + +Training metrics are saved in `history.json` and `history_lr.json`. + +Best model is saved in `<experiment_name>_my_model.h5`. + +## post-processing + +During the *make_forecast* method, all calculated forecasts of the neural network, persistence, ordinary least squared +and the target values with the regarding lead time are saved locally inside *forecast_path* as +`forecasts_<station>_test.nc`. + +All plots are created inside *plot_path*. \ No newline at end of file diff --git a/src/run_modules/experiment_setup.py b/src/run_modules/experiment_setup.py index b0c628b8966aa7cfedbf389552a9893acca297fa..44f1a3821dfacba146a0aabe8fb0254068d9e6d3 100644 --- a/src/run_modules/experiment_setup.py +++ b/src/run_modules/experiment_setup.py @@ -2,15 +2,14 @@ __author__ = "Lukas Leufen, Felix Kleinert" __date__ = '2019-11-15' -import logging import argparse -from typing import Union, Dict, Any +import logging import os +from typing import Union, Dict, Any from src import helpers from src.run_modules.run_environment import RunEnvironment - DEFAULT_STATIONS = ['DEBW107', 'DEBY081', 'DEBW013', 'DEBW076', 'DEBW087', 'DEBY052', 'DEBY032', 'DEBW022', 'DEBY004', 'DEBY020', 'DEBW030', 'DEBW037', 'DEBW031', 'DEBW015', 'DEBW073', 'DEBY039', 'DEBW038', 'DEBW081', 'DEBY075', 'DEBW040', 'DEBY053', 'DEBW059', 'DEBW027', 'DEBY072', 'DEBW042', 'DEBW039', 'DEBY001', @@ -33,19 +32,20 @@ class ExperimentSetup(RunEnvironment): window_lead_time=None, dimensions=None, interpolate_dim=None, interpolate_method=None, limit_nan_fill=None, train_start=None, train_end=None, val_start=None, val_end=None, test_start=None, test_end=None, use_all_stations_on_all_data_sets=True, trainable=False, fraction_of_train=None, - experiment_path=None, plot_path=None, forecast_path=None): + experiment_path=None, plot_path=None, forecast_path=None, overwrite_local_data=None, sampling="daily"): # create run framework super().__init__() # experiment setup - self._set_param("data_path", helpers.prepare_host()) + self._set_param("data_path", helpers.prepare_host(sampling=sampling)) self._set_param("trainable", trainable, default=False) self._set_param("fraction_of_training", fraction_of_train, default=0.8) # set experiment name exp_date = self._get_parser_args(parser_args).get("experiment_date") - exp_name, exp_path = helpers.set_experiment_name(experiment_date=exp_date, experiment_path=experiment_path) + exp_name, exp_path = helpers.set_experiment_name(experiment_date=exp_date, experiment_path=experiment_path, + sampling=sampling) self._set_param("experiment_name", exp_name) self._set_param("experiment_path", exp_path) helpers.check_path_and_create(self.data_store.get("experiment_path", "general")) @@ -67,9 +67,12 @@ class ExperimentSetup(RunEnvironment): self._set_param("station_type", station_type, default=None) self._set_param("variables", variables, default=list(self.data_store.get("var_all_dict", "general").keys())) self._set_param("statistics_per_var", statistics_per_var, default=self.data_store.get("var_all_dict", "general")) + self._compare_variables_and_statistics() self._set_param("start", start, default="1997-01-01", scope="general") self._set_param("end", end, default="2017-12-31", scope="general") self._set_param("window_history_size", window_history_size, default=13) + self._set_param("overwrite_local_data", overwrite_local_data, default=False, scope="general.preprocessing") + self._set_param("sampling", sampling) # target self._set_param("target_var", target_var, default="o3") @@ -94,6 +97,10 @@ class ExperimentSetup(RunEnvironment): self._set_param("start", test_start, default="2010-01-01", scope="general.test") self._set_param("end", test_end, default="2017-12-31", scope="general.test") + # train_val parameters + self._set_param("start", self.data_store.get("start", "general.train"), scope="general.train_val") + self._set_param("end", self.data_store.get("end", "general.val"), scope="general.train_val") + # use all stations on all data sets (train, val, test) self._set_param("use_all_stations_on_all_data_sets", use_all_stations_on_all_data_sets, default=True) @@ -117,6 +124,17 @@ class ExperimentSetup(RunEnvironment): else: return {} + def _compare_variables_and_statistics(self): + + logging.debug("check if all variables are included in statistics_per_var") + var = self.data_store.get("variables", "general") + stat = self.data_store.get("statistics_per_var", "general") + if not set(var).issubset(stat.keys()): + missing = set(var).difference(stat.keys()) + raise ValueError(f"Comparison of given variables and statistics_per_var show that not all requested " + f"variables are part of statistics_per_var. Please add also information on the missing " + f"statistics for the variables: {missing}") + if __name__ == "__main__": diff --git a/src/run_modules/model_setup.py b/src/run_modules/model_setup.py index a47ef67ad5781ff37ce812aa931dbd195d4513dc..7049af591cdf73434459d4c3f5f6c11e80ab64c0 100644 --- a/src/run_modules/model_setup.py +++ b/src/run_modules/model_setup.py @@ -2,20 +2,20 @@ __author__ = "Lukas Leufen, Felix Kleinert" __date__ = '2019-12-02' -import keras -from keras import losses -from keras.callbacks import ModelCheckpoint -from keras.regularizers import l2 -from keras.optimizers import SGD -import tensorflow as tf import logging import os -from src.run_modules.run_environment import RunEnvironment -from src.helpers import l_p_loss, LearningRateDecay -from src.model_modules.inception_model import InceptionModelBase +import keras +import tensorflow as tf +from keras import losses + +from src.helpers import l_p_loss from src.model_modules.flatten import flatten_tail -from src.model_modules.model_class import MyLittleModel +from src.model_modules.inception_model import InceptionModelBase +from src.model_modules.keras_extensions import HistoryAdvanced, ModelCheckpointAdvanced +# from src.model_modules.model_class import MyBranchedModel as MyModel +from src.model_modules.model_class import MyLittleModel as MyModel +from src.run_modules.run_environment import RunEnvironment class ModelSetup(RunEnvironment): @@ -29,13 +29,11 @@ class ModelSetup(RunEnvironment): exp_name = self.data_store.get("experiment_name", "general") self.scope = "general.model" self.checkpoint_name = os.path.join(path, f"{exp_name}_model-best.h5") + self.callbacks_name = os.path.join(path, f"{exp_name}_model-best-callbacks-%s.pickle") self._run() def _run(self): - # create checkpoint - self._set_checkpoint() - # set channels depending on inputs self._set_channels() @@ -49,6 +47,9 @@ class ModelSetup(RunEnvironment): if self.data_store.get("trainable", self.scope) is False: self.load_weights() + # create checkpoint + self._set_checkpoint() + # compile model self.compile_model() @@ -63,7 +64,17 @@ class ModelSetup(RunEnvironment): self.data_store.set("model", self.model, self.scope) def _set_checkpoint(self): - checkpoint = ModelCheckpoint(self.checkpoint_name, verbose=1, monitor='val_loss', save_best_only=True, mode='auto') + """ + Must be run after all callback functions that shall be tracked during training have been created (currently this + affects the learning rate decay and the advanced history [actually created in this method]). + """ + lr = self.data_store.get("lr_decay", scope="general.model") + hist = HistoryAdvanced() + self.data_store.set("hist", hist, scope="general.model") + callbacks = [{"callback": lr, "path": self.callbacks_name % "lr"}, + {"callback": hist, "path": self.callbacks_name % "hist"}] + checkpoint = ModelCheckpointAdvanced(filepath=self.checkpoint_name, verbose=1, monitor='val_loss', + save_best_only=True, mode='auto', callbacks=callbacks) self.data_store.set("checkpoint", checkpoint, self.scope) def load_weights(self): @@ -76,7 +87,7 @@ class ModelSetup(RunEnvironment): def build_model(self): args_list = ["window_history_size", "window_lead_time", "channels"] args = self.data_store.create_args_dict(args_list, self.scope) - self.model = MyLittleModel(**args) + self.model = MyModel(**args) self.get_model_settings() def get_model_settings(self): diff --git a/src/run_modules/post_processing.py b/src/run_modules/post_processing.py index 35d93dcbd932d1c298c0744fcd0205697576bb4c..03d2e36e8662a573b96c970747e9fe4445244e9b 100644 --- a/src/run_modules/post_processing.py +++ b/src/run_modules/post_processing.py @@ -5,38 +5,83 @@ __date__ = '2019-12-11' import logging import os +import keras import numpy as np import pandas as pd import xarray as xr -import statsmodels.api as sm -from src.run_modules.run_environment import RunEnvironment -from src.data_handling.data_distributor import Distributor -from src.model_modules.linear_model import OrdinaryLeastSquaredModel from src import statistics -from src import helpers +from src.data_handling.data_distributor import Distributor +from src.data_handling.data_generator import DataGenerator +from src.datastore import NameNotFoundInDataStore from src.helpers import TimeTracking +from src.model_modules.linear_model import OrdinaryLeastSquaredModel +from src.plotting.postprocessing_plotting import PlotMonthlySummary, PlotStationMap, PlotClimatologicalSkillScore, \ + PlotCompetitiveSkillScore, PlotTimeSeries +from src.plotting.postprocessing_plotting import plot_conditional_quantiles +from src.run_modules.run_environment import RunEnvironment class PostProcessing(RunEnvironment): def __init__(self): super().__init__() - self.model = self.data_store.get("best_model", "general") + self.model: keras.Model = self._load_model() self.ols_model = None - self.batch_size = self.data_store.get("batch_size", "general.model") - self.test_data = self.data_store.get("generator", "general.test") + self.batch_size: int = self.data_store.get_default("batch_size", "general.model", 64) + self.test_data: DataGenerator = self.data_store.get("generator", "general.test") self.test_data_distributed = Distributor(self.test_data, self.model, self.batch_size) - self.train_data = self.data_store.get("generator", "general.train") + self.train_data: DataGenerator = self.data_store.get("generator", "general.train") + self.train_val_data: DataGenerator = self.data_store.get("generator", "general.train_val") + self.plot_path: str = self.data_store.get("plot_path", "general") + self.target_var = self.data_store.get("target_var", "general") + self._sampling = self.data_store.get("sampling", "general") + self.skill_scores = None self._run() def _run(self): - self.train_ols_model() - preds_for_all_stations = self.make_prediction() + with TimeTracking(): + self.train_ols_model() + logging.info("take a look on the next reported time measure. If this increases a lot, one should think to " + "skip make_prediction() whenever it is possible to save time.") + with TimeTracking(): + self.make_prediction() + logging.info("take a look on the next reported time measure. If this increases a lot, one should think to " + "skip make_prediction() whenever it is possible to save time.") + self.skill_scores = self.calculate_skill_scores() + self.plot() + + def _load_model(self): + try: + model = self.data_store.get("best_model", "general") + except NameNotFoundInDataStore: + logging.info("no model saved in data store. trying to load model from experiment") + path = self.data_store.get("experiment_path", "general") + name = f"{self.data_store.get('experiment_name', 'general')}_my_model.h5" + model_name = os.path.join(path, name) + model = keras.models.load_model(model_name) + return model + + def plot(self): + logging.debug("Run plotting routines...") + path = self.data_store.get("forecast_path", "general") + + plot_conditional_quantiles(self.test_data.stations, pred_name="CNN", ref_name="orig", + forecast_path=path, plot_name_affix="cali-ref", plot_folder=self.plot_path) + plot_conditional_quantiles(self.test_data.stations, pred_name="orig", ref_name="CNN", + forecast_path=path, plot_name_affix="like-bas", plot_folder=self.plot_path) + PlotStationMap(generators={'b': self.test_data}, plot_folder=self.plot_path) + PlotMonthlySummary(self.test_data.stations, path, r"forecasts_%s_test.nc", self.target_var, + plot_folder=self.plot_path) + PlotClimatologicalSkillScore(self.skill_scores[1], plot_folder=self.plot_path, model_setup="CNN") + PlotClimatologicalSkillScore(self.skill_scores[1], plot_folder=self.plot_path, score_only=False, + extra_name_tag="all_terms_", model_setup="CNN") + PlotCompetitiveSkillScore(self.skill_scores[0], plot_folder=self.plot_path, model_setup="CNN") + PlotTimeSeries(self.test_data.stations, path, r"forecasts_%s_test.nc", plot_folder=self.plot_path, sampling=self._sampling) def calculate_test_score(self): - test_score = self.model.evaluate(generator=self.test_data_distributed.distribute_on_batches(), - use_multiprocessing=False, verbose=0, steps=1) + test_score = self.model.evaluate_generator(generator=self.test_data_distributed.distribute_on_batches(), + use_multiprocessing=False, verbose=0, steps=1) logging.info(f"test score = {test_score}") self._save_test_score(test_score) @@ -49,16 +94,16 @@ class PostProcessing(RunEnvironment): def train_ols_model(self): self.ols_model = OrdinaryLeastSquaredModel(self.train_data) - def make_prediction(self, freq="1D"): - nn_prediction_all_stations = [] - for i, v in enumerate(self.test_data): + def make_prediction(self): + logging.debug("start make_prediction") + for i, _ in enumerate(self.test_data): data = self.test_data.get_data_generator(i) nn_prediction, persistence_prediction, ols_prediction = self._create_empty_prediction_arrays(data, count=3) - input_data = self.test_data[i][0] + input_data = data.get_transposed_history() # get scaling parameters - mean, std, transformation_method = data.get_transformation_information(variable='o3') + mean, std, transformation_method = data.get_transformation_information(variable=self.target_var) # nn forecast nn_prediction = self._create_nn_forecast(input_data, nn_prediction, mean, std, transformation_method) @@ -74,7 +119,7 @@ class PostProcessing(RunEnvironment): orig_pred = self._create_orig_forecast(data, None, mean, std, transformation_method) # merge all predictions - full_index = self.create_fullindex(data.data.indexes['datetime'], freq) + full_index = self.create_fullindex(data.data.indexes['datetime'], self._get_frequency()) all_predictions = self.create_forecast_arrays(full_index, list(data.label.indexes['window']), CNN=nn_prediction, persi=persistence_prediction, @@ -86,22 +131,24 @@ class PostProcessing(RunEnvironment): file = os.path.join(path, f"forecasts_{data.station[0]}_test.nc") all_predictions.to_netcdf(file) - # save nn forecast to return variable - nn_prediction_all_stations.append(nn_prediction) - return nn_prediction_all_stations + def _get_frequency(self): + getter = {"daily": "1D", "hourly": "1H"} + return getter.get(self._sampling, None) @staticmethod - def _create_orig_forecast(data, placeholder, mean, std, transformation_method): - return statistics.apply_inverse_transformation(data.label, mean, std, transformation_method) + def _create_orig_forecast(data, _, mean, std, transformation_method): + return statistics.apply_inverse_transformation(data.label.copy(), mean, std, transformation_method) def _create_ols_forecast(self, input_data, ols_prediction, mean, std, transformation_method): tmp_ols = self.ols_model.predict(input_data) tmp_ols = statistics.apply_inverse_transformation(tmp_ols, mean, std, transformation_method) - ols_prediction.values = np.swapaxes(np.expand_dims(tmp_ols, axis=1), 2, 0) + tmp_ols = np.expand_dims(tmp_ols, axis=1) + target_shape = ols_prediction.values.shape + ols_prediction.values = np.swapaxes(tmp_ols, 2, 0) if target_shape != tmp_ols.shape else tmp_ols return ols_prediction def _create_persistence_forecast(self, input_data, persistence_prediction, mean, std, transformation_method): - tmp_persi = input_data.sel({'window': 0, 'variables': 'o3'}) + tmp_persi = input_data.sel({'window': 0, 'variables': self.target_var}) tmp_persi = statistics.apply_inverse_transformation(tmp_persi, mean, std, transformation_method) window_lead_time = self.data_store.get("window_lead_time", "general") persistence_prediction.values = np.expand_dims(np.tile(tmp_persi.squeeze('Stations'), (window_lead_time, 1)), @@ -109,14 +156,30 @@ class PostProcessing(RunEnvironment): return persistence_prediction def _create_nn_forecast(self, input_data, nn_prediction, mean, std, transformation_method): + """ + create the nn forecast for given input data. Inverse transformation is applied to the forecast to get the output + in the original space. Furthermore, only the output of the main branch is returned (not all minor branches, if + the network has multiple output branches). The main branch is defined to be the last entry of all outputs. + :param input_data: + :param nn_prediction: + :param mean: + :param std: + :param transformation_method: + :return: + """ tmp_nn = self.model.predict(input_data) tmp_nn = statistics.apply_inverse_transformation(tmp_nn, mean, std, transformation_method) - nn_prediction.values = np.swapaxes(np.expand_dims(tmp_nn, axis=1), 2, 0) + if tmp_nn.ndim == 3: + nn_prediction.values = np.swapaxes(np.expand_dims(tmp_nn[-1, ...], axis=1), 2, 0) + elif tmp_nn.ndim == 2: + nn_prediction.values = np.swapaxes(np.expand_dims(tmp_nn, axis=1), 2, 0) + else: + raise NotImplementedError(f"Number of dimension of model output must be 2 or 3, but not {tmp_nn.dims}.") return nn_prediction @staticmethod def _create_empty_prediction_arrays(generator, count=1): - return [generator.label.copy()] * count + return [generator.label.copy() for _ in range(count)] @staticmethod def create_fullindex(df, freq): @@ -130,7 +193,7 @@ class PostProcessing(RunEnvironment): elif isinstance(df, xr.DataArray): earliest = df.index[0].values latest = df.index[-1].values - elif isinstance(df, pd.core.indexes.datetimes.DatetimeIndex): + elif isinstance(df, pd.DatetimeIndex): earliest = df[0] latest = df[-1] else: @@ -162,4 +225,27 @@ class PostProcessing(RunEnvironment): res.loc[match_index, :, k] = v.sel({'datetime': match_index}).squeeze('Stations').transpose() return res - + def _get_external_data(self, station): + try: + data = self.train_val_data.get_data_generator(station) + mean, std, transformation_method = data.get_transformation_information(variable=self.target_var) + external_data = self._create_orig_forecast(data, None, mean, std, transformation_method) + external_data = external_data.squeeze("Stations").sel(window=1).drop(["window", "Stations", "variables"]) + return external_data.rename({'datetime': 'index'}) + except KeyError: + return None + + def calculate_skill_scores(self): + path = self.data_store.get("forecast_path", "general") + window_lead_time = self.data_store.get("window_lead_time", "general") + skill_score_competitive = {} + skill_score_climatological = {} + for station in self.test_data.stations: + file = os.path.join(path, f"forecasts_{station}_test.nc") + data = xr.open_dataarray(file) + skill_score = statistics.SkillScores(data) + external_data = self._get_external_data(station) + skill_score_competitive[station] = skill_score.skill_scores(window_lead_time) + skill_score_climatological[station] = skill_score.climatological_skill_scores(external_data, + window_lead_time) + return skill_score_competitive, skill_score_climatological diff --git a/src/run_modules/pre_processing.py b/src/run_modules/pre_processing.py index 9faae7943eb4cd4adb6572c524f5d1795aaa65fe..4660a8116b6d0b860a7d0d50b92cee5e0deb77d8 100644 --- a/src/run_modules/pre_processing.py +++ b/src/run_modules/pre_processing.py @@ -7,12 +7,12 @@ from typing import Tuple, Dict, List from src.data_handling.data_generator import DataGenerator from src.helpers import TimeTracking -from src.run_modules.run_environment import RunEnvironment from src.join import EmptyQueryResult - +from src.run_modules.run_environment import RunEnvironment DEFAULT_ARGS_LIST = ["data_path", "network", "stations", "variables", "interpolate_dim", "target_dim", "target_var"] -DEFAULT_KWARGS_LIST = ["limit_nan_fill", "window_history_size", "window_lead_time", "statistics_per_var", "station_type"] +DEFAULT_KWARGS_LIST = ["limit_nan_fill", "window_history_size", "window_lead_time", "statistics_per_var", + "station_type", "overwrite_local_data", "start", "end", "sampling"] class PreProcessing(RunEnvironment): @@ -33,48 +33,51 @@ class PreProcessing(RunEnvironment): self._run() def _run(self): - args = self.data_store.create_args_dict(DEFAULT_ARGS_LIST) - kwargs = self.data_store.create_args_dict(DEFAULT_KWARGS_LIST) - valid_stations = self.check_valid_stations(args, kwargs, self.data_store.get("stations", "general")) + args = self.data_store.create_args_dict(DEFAULT_ARGS_LIST, scope="general.preprocessing") + kwargs = self.data_store.create_args_dict(DEFAULT_KWARGS_LIST, scope="general.preprocessing") + valid_stations = self.check_valid_stations(args, kwargs, self.data_store.get("stations", "general"), load_tmp=False) self.data_store.set("stations", valid_stations, "general") self.split_train_val_test() + self.report_pre_processing() def report_pre_processing(self): - logging.info(20 * '##') + logging.debug(20 * '##') n_train = len(self.data_store.get('generator', 'general.train')) n_val = len(self.data_store.get('generator', 'general.val')) n_test = len(self.data_store.get('generator', 'general.test')) n_total = n_train + n_val + n_test - logging.info(f"Number of all stations: {n_total}") - logging.info(f"Number of training stations: {n_train}") - logging.info(f"Number of val stations: {n_val}") - logging.info(f"Number of test stations: {n_test}") - logging.info(f"TEST SHAPE OF GENERATOR CALL: {self.data_store.get('generator', 'general.test')[0][0].shape}" - f"{self.data_store.get('generator', 'general.test')[0][1].shape}") + logging.debug(f"Number of all stations: {n_total}") + logging.debug(f"Number of training stations: {n_train}") + logging.debug(f"Number of val stations: {n_val}") + logging.debug(f"Number of test stations: {n_test}") + logging.debug(f"TEST SHAPE OF GENERATOR CALL: {self.data_store.get('generator', 'general.test')[0][0].shape}" + f"{self.data_store.get('generator', 'general.test')[0][1].shape}") def split_train_val_test(self): fraction_of_training = self.data_store.get("fraction_of_training", "general") stations = self.data_store.get("stations", "general") - train_index, val_index, test_index = self.split_set_indices(len(stations), fraction_of_training) - for (ind, scope) in zip([train_index, val_index, test_index], ["train", "val", "test"]): + train_index, val_index, test_index, train_val_index = self.split_set_indices(len(stations), fraction_of_training) + subset_names = ["train", "val", "test", "train_val"] + for (ind, scope) in zip([train_index, val_index, test_index, train_val_index], subset_names): self.create_set_split(ind, scope) @staticmethod - def split_set_indices(total_length: int, fraction: float) -> Tuple[slice, slice, slice]: + def split_set_indices(total_length: int, fraction: float) -> Tuple[slice, slice, slice, slice]: """ create the training, validation and test subset slice indices for given total_length. The test data consists on (1-fraction) of total_length (fraction*len:end). Train and validation data therefore are made from fraction of total_length (0:fraction*len). Train and validation data is split by the factor 0.8 for train and 0.2 for - validation. + validation. In addition, split_set_indices returns also the combination of training and validation subset. :param total_length: list with all objects to split :param fraction: ratio between test and union of train/val data - :return: slices for each subset in the order: train, val, test + :return: slices for each subset in the order: train, val, test, train_val """ pos_test_split = int(total_length * fraction) train_index = slice(0, int(pos_test_split * 0.8)) val_index = slice(int(pos_test_split * 0.8), pos_test_split) test_index = slice(pos_test_split, total_length) - return train_index, val_index, test_index + train_val_index = slice(0, pos_test_split) + return train_index, val_index, test_index, train_val_index def create_set_split(self, index_list, set_name): scope = f"general.{set_name}" @@ -93,7 +96,7 @@ class PreProcessing(RunEnvironment): self.data_store.set("generator", data_set, scope) @staticmethod - def check_valid_stations(args: Dict, kwargs: Dict, all_stations: List[str]): + def check_valid_stations(args: Dict, kwargs: Dict, all_stations: List[str], load_tmp=True): """ Check if all given stations in `all_stations` are valid. Valid means, that there is data available for the given time range (is included in `kwargs`). The shape and the loading time are logged in debug mode. @@ -114,9 +117,10 @@ class PreProcessing(RunEnvironment): for station in all_stations: t_inner.run() try: - (history, label) = data_gen[station] + # (history, label) = data_gen[station] + data = data_gen.get_data_generator(key=station, local_tmp_storage=load_tmp) valid_stations.append(station) - logging.debug(f"{station}: history_shape = {history.shape}") + logging.debug(f'{station}: history_shape = {data.history.transpose("datetime", "window", "Stations", "variables").shape}') logging.debug(f"{station}: loading time = {t_inner}") except (AttributeError, EmptyQueryResult): continue diff --git a/src/run_modules/training.py b/src/run_modules/training.py index 272609a31a3e3c91d6857ed841d5dd2783c66f35..e2a98f27c65e6050b0edae2bbc178abbf97ab646 100644 --- a/src/run_modules/training.py +++ b/src/run_modules/training.py @@ -1,29 +1,32 @@ __author__ = "Lukas Leufen, Felix Kleinert" __date__ = '2019-12-05' +import json import logging import os -import json +import pickle + import keras -from src.run_modules.run_environment import RunEnvironment from src.data_handling.data_distributor import Distributor +from src.model_modules.keras_extensions import LearningRateDecay, ModelCheckpointAdvanced from src.plotting.training_monitoring import PlotModelHistory, PlotModelLearningRate -from src.helpers import LearningRateDecay +from src.run_modules.run_environment import RunEnvironment class Training(RunEnvironment): def __init__(self): super().__init__() - self.model = self.data_store.get("model", "general.model") + self.model: keras.Model = self.data_store.get("model", "general.model") self.train_set = None self.val_set = None self.test_set = None self.batch_size = self.data_store.get("batch_size", "general.model") self.epochs = self.data_store.get("epochs", "general.model") - self.checkpoint = self.data_store.get("checkpoint", "general.model") + self.checkpoint: ModelCheckpointAdvanced = self.data_store.get("checkpoint", "general.model") self.lr_sc = self.data_store.get("lr_decay", "general.model") + self.hist = self.data_store.get("hist", "general.model") self.experiment_name = self.data_store.get("experiment_name", "general") self._run() @@ -35,7 +38,7 @@ class Training(RunEnvironment): 2) make_predict_function(): create predict function before distribution on multiple nodes (detailed information in method description) 3) train(): - train model and save callbacks + start or resume training of model and save callbacks 4) save_model(): save best model from training as final model """ @@ -73,17 +76,42 @@ class Training(RunEnvironment): def train(self) -> None: """ Perform training using keras fit_generator(). Callbacks are stored locally in the experiment directory. Best - model from training is saved for class variable model. + model from training is saved for class variable model. If the file path of checkpoint is not empty, this method + assumes, that this is not a new training starting from the very beginning, but a resumption from a previous + started but interrupted training (or a stopped and now continued training). Train will automatically load the + locally stored information and the corresponding model and proceed with the already started training. """ logging.info(f"Train with {len(self.train_set)} mini batches.") - history = self.model.fit_generator(generator=self.train_set.distribute_on_batches(), - steps_per_epoch=len(self.train_set), - epochs=self.epochs, - verbose=2, - validation_data=self.val_set.distribute_on_batches(), - validation_steps=len(self.val_set), - callbacks=[self.checkpoint, self.lr_sc]) - self.save_callbacks(history) + if not os.path.exists(self.checkpoint.filepath): + history = self.model.fit_generator(generator=self.train_set.distribute_on_batches(), + steps_per_epoch=len(self.train_set), + epochs=self.epochs, + verbose=2, + validation_data=self.val_set.distribute_on_batches(), + validation_steps=len(self.val_set), + callbacks=[self.lr_sc, self.hist, self.checkpoint]) + else: + logging.info("Found locally stored model and checkpoints. Training is resumed from the last checkpoint.") + lr_filepath = self.checkpoint.callbacks[0]["path"] + hist_filepath = self.checkpoint.callbacks[1]["path"] + self.lr_sc = pickle.load(open(lr_filepath, "rb")) + self.hist = pickle.load(open(hist_filepath, "rb")) + self.model = keras.models.load_model(self.checkpoint.filepath) + initial_epoch = max(self.hist.epoch) + 1 + callbacks = [{"callback": self.lr_sc, "path": lr_filepath}, + {"callback": self.hist, "path": hist_filepath}] + self.checkpoint.update_callbacks(callbacks) + self.checkpoint.update_best(self.hist) + _ = self.model.fit_generator(generator=self.train_set.distribute_on_batches(), + steps_per_epoch=len(self.train_set), + epochs=self.epochs, + verbose=2, + validation_data=self.val_set.distribute_on_batches(), + validation_steps=len(self.val_set), + callbacks=[self.lr_sc, self.hist, self.checkpoint], + initial_epoch=initial_epoch) + history = self.hist + self.save_callbacks_as_json(history) self.load_best_model(self.checkpoint.filepath) self.create_monitoring_plots(history, self.lr_sc) @@ -110,7 +138,7 @@ class Training(RunEnvironment): except OSError: logging.info('no weights to reload...') - def save_callbacks(self, history: keras.callbacks.History) -> None: + def save_callbacks_as_json(self, history: keras.callbacks.History) -> None: """ Save callbacks (history, learning rate) of training. * history.history -> history.json @@ -134,5 +162,17 @@ class Training(RunEnvironment): """ path = self.data_store.get("plot_path", "general") name = self.data_store.get("experiment_name", "general") - PlotModelHistory(filename=os.path.join(path, f"{name}_history_loss_val_loss.pdf"), history=history) + + # plot history of loss and mse (if available) + filename = os.path.join(path, f"{name}_history_loss.pdf") + PlotModelHistory(filename=filename, history=history) + multiple_branches_used = len(history.model.output_names) > 1 # means that there are multiple output branches + if multiple_branches_used: + filename = os.path.join(path, f"{name}_history_main_loss.pdf") + PlotModelHistory(filename=filename, history=history, main_branch=True) + if "mean_squared_error" in history.model.metrics_names: + filename = os.path.join(path, f"{name}_history_main_mse.pdf") + PlotModelHistory(filename=filename, history=history, plot_metric="mse", main_branch=multiple_branches_used) + + # plot learning rate PlotModelLearningRate(filename=os.path.join(path, f"{name}_history_learning_rate.pdf"), lr_sc=lr_sc) diff --git a/src/statistics.py b/src/statistics.py index 6f34187e32949910df5762a45d868701920b610f..df73784df830d5f7b96bf0fcd18a65d362516f12 100644 --- a/src/statistics.py +++ b/src/statistics.py @@ -1,6 +1,11 @@ -__author__ = 'Lukas Leufen' +from scipy import stats + +from src.run_modules.run_environment import RunEnvironment + +__author__ = 'Lukas Leufen, Felix Kleinert' __date__ = '2019-10-23' +import numpy as np import xarray as xr import pandas as pd from typing import Union, Tuple @@ -70,3 +75,145 @@ def centre_inverse(data: Data, mean: Data) -> Data: :return: """ return data + mean + + +def mean_squared_error(a, b): + return np.square(a - b).mean() + + +class SkillScores(RunEnvironment): + + def __init__(self, internal_data): + super().__init__() + self.internal_data = internal_data + + def skill_scores(self, window_lead_time): + ahead_names = list(range(1, window_lead_time + 1)) + skill_score = pd.DataFrame(index=['cnn-persi', 'ols-persi', 'cnn-ols']) + for iahead in ahead_names: + data = self.internal_data.sel(ahead=iahead) + skill_score[iahead] = [self.general_skill_score(data, forecast_name="CNN", reference_name="persi"), + self.general_skill_score(data, forecast_name="OLS", reference_name="persi"), + self.general_skill_score(data, forecast_name="CNN", reference_name="OLS")] + return skill_score + + def climatological_skill_scores(self, external_data, window_lead_time): + ahead_names = list(range(1, window_lead_time + 1)) + + all_terms = ['AI', 'AII', 'AIII', 'AIV', 'BI', 'BII', 'BIV', 'CI', 'CIV', 'CASE I', 'CASE II', 'CASE III', + 'CASE IV'] + skill_score = xr.DataArray(np.full((len(all_terms), len(ahead_names)), np.nan), coords=[all_terms, ahead_names], + dims=['terms', 'ahead']) + + for iahead in ahead_names: + + data = self.internal_data.sel(ahead=iahead) + + skill_score.loc[["CASE I", "AI", "BI", "CI"], iahead] = np.stack(self._climatological_skill_score( + data, mu_type=1, forecast_name="CNN").values.flatten()) + + skill_score.loc[["CASE II", "AII", "BII"], iahead] = np.stack(self._climatological_skill_score( + data, mu_type=2, forecast_name="CNN").values.flatten()) + + if external_data is not None: + skill_score.loc[["CASE III", "AIII"], iahead] = np.stack(self._climatological_skill_score( + data, mu_type=3, forecast_name="CNN", + external_data=external_data).values.flatten()) + + skill_score.loc[["CASE IV", "AIV", "BIV", "CIV"], iahead] = np.stack(self._climatological_skill_score( + data, mu_type=4, forecast_name="CNN", + external_data=external_data).values.flatten()) + + return skill_score + + def _climatological_skill_score(self, data, mu_type=1, observation_name="orig", forecast_name="CNN", external_data=None): + kwargs = {"external_data": external_data} if external_data is not None else {} + return self.__getattribute__(f"skill_score_mu_case_{mu_type}")(data, observation_name, forecast_name, **kwargs) + + @staticmethod + def general_skill_score(data, observation_name="orig", forecast_name="CNN", reference_name="persi"): + data = data.dropna("index") + observation = data.sel(type=observation_name) + forecast = data.sel(type=forecast_name) + reference = data.sel(type=reference_name) + mse = mean_squared_error + skill_score = 1 - mse(observation, forecast) / mse(observation, reference) + return skill_score.values + + @staticmethod + def skill_score_pre_calculations(data, observation_name, forecast_name): + + data = data.loc[..., [observation_name, forecast_name]].drop("ahead") + data = data.dropna("index") + + mean = data.mean("index") + sigma = np.sqrt(data.var("index")) + # r, p = stats.spearmanr(data.loc[..., [forecast_name, observation_name]]) + r, p = stats.pearsonr(data.loc[..., forecast_name], data.loc[..., observation_name]) + + AI = np.array(r ** 2) + BI = ((r - (sigma.loc[..., forecast_name] / sigma.loc[..., observation_name])) ** 2).values + CI = (((mean.loc[..., forecast_name] - mean.loc[..., observation_name]) / sigma.loc[ + ..., observation_name]) ** 2).values + + suffix = {"mean": mean, "sigma": sigma, "r": r, "p": p} + return AI, BI, CI, data, suffix + + def skill_score_mu_case_1(self, data, observation_name="orig", forecast_name="CNN"): + AI, BI, CI, data, _ = self.skill_score_pre_calculations(data, observation_name, forecast_name) + skill_score = np.array(AI - BI - CI) + return pd.DataFrame({"skill_score": [skill_score], "AI": [AI], "BI": [BI], "CI": [CI]}).to_xarray().to_array() + + def skill_score_mu_case_2(self, data, observation_name="orig", forecast_name="CNN"): + AI, BI, CI, data, suffix = self.skill_score_pre_calculations(data, observation_name, forecast_name) + monthly_mean = self.create_monthly_mean_from_daily_data(data) + data = xr.concat([data, monthly_mean], dim="type") + sigma = suffix["sigma"] + sigma_monthly = np.sqrt(monthly_mean.var()) + # r, p = stats.spearmanr(data.loc[..., [observation_name, observation_name + "X"]]) + r, p = stats.pearsonr(data.loc[..., observation_name], data.loc[..., observation_name + "X"]) + AII = np.array(r ** 2) + BII = ((r - sigma_monthly / sigma.loc[observation_name]) ** 2).values + skill_score = np.array((AI - BI - CI - AII + BII) / (1 - AII + BII)) + return pd.DataFrame({"skill_score": [skill_score], "AII": [AII], "BII": [BII]}).to_xarray().to_array() + + def skill_score_mu_case_3(self, data, observation_name="orig", forecast_name="CNN", external_data=None): + AI, BI, CI, data, suffix = self.skill_score_pre_calculations(data, observation_name, forecast_name) + mean, sigma = suffix["mean"], suffix["sigma"] + AIII = (((external_data.mean().values - mean.loc[observation_name]) / sigma.loc[observation_name])**2).values + skill_score = np.array((AI - BI - CI + AIII) / 1 + AIII) + return pd.DataFrame({"skill_score": [skill_score], "AIII": [AIII]}).to_xarray().to_array() + + def skill_score_mu_case_4(self, data, observation_name="orig", forecast_name="CNN", external_data=None): + AI, BI, CI, data, suffix = self.skill_score_pre_calculations(data, observation_name, forecast_name) + monthly_mean_external = self.create_monthly_mean_from_daily_data(external_data, columns=data.type.values, index=data.index) + data = xr.concat([data, monthly_mean_external], dim="type") + mean, sigma = suffix["mean"], suffix["sigma"] + monthly_mean_external = self.create_monthly_mean_from_daily_data(external_data, columns=data.type.values) + mean_external = monthly_mean_external.mean() + sigma_external = np.sqrt(monthly_mean_external.var()) + + # r_mu, p_mu = stats.spearmanr(data.loc[..., [observation_name, observation_name+'X']]) + r_mu, p_mu = stats.pearsonr(data.loc[..., observation_name], data.loc[..., observation_name + "X"]) + + AIV = np.array(r_mu**2) + BIV = ((r_mu - sigma_external / sigma.loc[observation_name])**2).values + CIV = (((mean_external - mean.loc[observation_name]) / sigma.loc[observation_name])**2).values + skill_score = np.array((AI - BI - CI - AIV + BIV + CIV) / (1 - AIV + BIV + CIV)) + return pd.DataFrame({"skill_score": [skill_score], "AIV": [AIV], "BIV": [BIV], "CIV": CIV}).to_xarray().to_array() + + @staticmethod + def create_monthly_mean_from_daily_data(data, columns=None, index=None): + if columns is None: + columns = data.type.values + if index is None: + index = data.index + coordinates = [index, [v + "X" for v in list(columns)]] + empty_data = np.full((len(index), len(columns)), np.nan) + monthly_mean = xr.DataArray(empty_data, coords=coordinates, dims=["index", "type"]) + mu = data.groupby("index.month").mean() + + for month in mu.month: + monthly_mean[monthly_mean.index.dt.month == month, :] = mu[mu.month == month].values + + return monthly_mean diff --git a/test/test_data_handling/test_data_distributor.py b/test/test_data_handling/test_data_distributor.py index cb51f20c8771ec49116731f02c7b462a62405394..4c6dbb1c38f2e4a49e53883fbe3cb33cb565118a 100644 --- a/test/test_data_handling/test_data_distributor.py +++ b/test/test_data_handling/test_data_distributor.py @@ -1,6 +1,5 @@ import math import os -import shutil import keras import numpy as np diff --git a/test/test_data_handling/test_data_generator.py b/test/test_data_handling/test_data_generator.py index 879436afddb8da8d11d6cc585da7c703aa12ef8a..142acd166604951352ad6686548c2cb76f609ce0 100644 --- a/test/test_data_handling/test_data_generator.py +++ b/test/test_data_handling/test_data_generator.py @@ -1,7 +1,12 @@ -import pytest import os + +import pytest + import shutil +import numpy as np +import pickle from src.data_handling.data_generator import DataGenerator +from src.data_handling.data_preparation import DataPrep class TestDataGenerator: @@ -15,7 +20,13 @@ class TestDataGenerator: @pytest.fixture def gen(self): return DataGenerator(os.path.join(os.path.dirname(__file__), 'data'), 'AIRBASE', 'DEBW107', ['o3', 'temp'], - 'datetime', 'variables', 'o3') + 'datetime', 'variables', 'o3', start=2010, end=2014) + + class DummyDataPrep: + def __init__(self, data): + self.station = "DEBW107" + self.variables = ["o3", "temp"] + self.data = data def test_init(self, gen): assert gen.data_path == os.path.join(os.path.dirname(__file__), 'data') @@ -31,28 +42,20 @@ class TestDataGenerator: assert gen.window_history_size == 7 assert gen.window_lead_time == 4 assert gen.transform_method == "standardise" - assert gen.kwargs == {} + assert gen.kwargs == {"start": 2010, "end": 2014} def test_repr(self, gen): path = os.path.join(os.path.dirname(__file__), 'data') assert gen.__repr__().rstrip() == f"DataGenerator(path='{path}', network='AIRBASE', stations=['DEBW107'], "\ f"variables=['o3', 'temp'], station_type=None, interpolate_dim='datetime', " \ - f"target_dim='variables', target_var='o3', **{{}})".rstrip() + f"target_dim='variables', target_var='o3', **{{'start': 2010, 'end': 2014}})"\ + .rstrip() def test_len(self, gen): assert len(gen) == 1 gen.stations = ['station1', 'station2', 'station3'] assert len(gen) == 3 - def test_getitem(self, gen): - gen.kwargs = {'statistics_per_var': {'o3': 'dma8eu', 'temp': 'maximum'}} - station = gen["DEBW107"] - assert len(station) == 2 - assert station[0].Stations.data == "DEBW107" - assert station[0].data.shape[1:] == (8, 1, 2) - assert station[1].data.shape[-1] == gen.window_lead_time - assert station[0].data.shape[1] == gen.window_history_size + 1 - def test_iter(self, gen): assert hasattr(gen, '_iterator') is False iter(gen) @@ -64,6 +67,15 @@ class TestDataGenerator: for i, d in enumerate(gen, start=1): assert i == gen._iterator + def test_getitem(self, gen): + gen.kwargs = {'statistics_per_var': {'o3': 'dma8eu', 'temp': 'maximum'}} + station = gen["DEBW107"] + assert len(station) == 2 + assert station[0].Stations.data == "DEBW107" + assert station[0].data.shape[1:] == (8, 1, 2) + assert station[1].data.shape[-1] == gen.window_lead_time + assert station[0].data.shape[1] == gen.window_history_size + 1 + def test_get_station_key(self, gen): gen.stations.append("DEBW108") f = gen.get_station_key @@ -85,3 +97,38 @@ class TestDataGenerator: with pytest.raises(KeyError) as e: f(6.5) assert "key has to be from Union[str, int]. Given was 6.5 (float)" + + def test_get_data_generator(self, gen): + gen.kwargs = {"statistics_per_var": {'o3': 'dma8eu', 'temp': 'maximum'}} + file = os.path.join(gen.data_path_tmp, f"DEBW107_{'_'.join(sorted(gen.variables))}_None_None_.pickle") + if os.path.exists(file): + os.remove(file) + assert not os.path.exists(file) + assert isinstance(gen.get_data_generator("DEBW107", local_tmp_storage=False), DataPrep) + t = os.stat(file).st_ctime + assert os.path.exists(file) + assert isinstance(gen.get_data_generator("DEBW107"), DataPrep) + assert os.stat(file).st_mtime == t + os.remove(file) + assert isinstance(gen.get_data_generator("DEBW107"), DataPrep) + assert os.stat(file).st_ctime > t + + def test_save_pickle_data(self, gen): + file = os.path.join(gen.data_path_tmp, f"DEBW107_{'_'.join(sorted(gen.variables))}_2010_2014_.pickle") + if os.path.exists(file): + os.remove(file) + assert not os.path.exists(file) + data = self.DummyDataPrep(np.ones((10, 2))) + gen._save_pickle_data(data) + assert os.path.exists(file) + os.remove(file) + + def test_load_pickle_data(self, gen): + file = os.path.join(gen.data_path_tmp, f"DEBW107_{'_'.join(sorted(gen.variables))}_2010_2014_.pickle") + data = self.DummyDataPrep(np.ones((10, 2))) + with open(file, "wb") as f: + pickle.dump(data, f) + assert os.path.exists(file) + res = gen._load_pickle_data("DEBW107", ["o3", "temp"]).data + assert np.testing.assert_almost_equal(res, np.ones((10, 2))) is None + os.remove(file) diff --git a/test/test_data_handling/test_data_preparation.py b/test/test_data_handling/test_data_preparation.py index 12b619d9e31990f6cc24216ff84ad9d030265e36..72bacaf9cc1e5a9dc9736b8e8eb7161f35d8ea69 100644 --- a/test/test_data_handling/test_data_preparation.py +++ b/test/test_data_handling/test_data_preparation.py @@ -1,12 +1,15 @@ -import pytest +import datetime as dt import os -from src.data_handling.data_preparation import DataPrep -from src.join import EmptyQueryResult +from operator import itemgetter +import logging + import numpy as np -import xarray as xr -import datetime as dt import pandas as pd -from operator import itemgetter +import pytest +import xarray as xr + +from src.data_handling.data_preparation import DataPrep +from src.join import EmptyQueryResult class TestDataPrep: @@ -17,6 +20,18 @@ class TestDataPrep: station_type='background', test='testKWARGS', statistics_per_var={'o3': 'dma8eu', 'temp': 'maximum'}) + @pytest.fixture + def data_prep_no_init(self): + d = object.__new__(DataPrep) + d.path = os.path.join(os.path.abspath(os.path.dirname(__file__)), 'data') + d.network = 'UBA' + d.station = ['DEBW107'] + d.variables = ['o3', 'temp'] + d.station_type = "background" + d.sampling = "daily" + d.kwargs = None + return d + def test_init(self, data): assert data.path == os.path.join(os.path.abspath(os.path.dirname(__file__)), 'data') assert data.network == 'AIRBASE' @@ -31,16 +46,79 @@ class TestDataPrep: with pytest.raises(NotImplementedError): DataPrep('data/', 'dummy', 'DEBW107', ['o3', 'temp']) - def test_repr(self): - d = object.__new__(DataPrep) - d.path = 'data/test' - d.network = 'dummy' - d.station = ['DEBW107'] - d.variables = ['o3', 'temp'] - d.station_type = "traffic" - d.kwargs = None - assert d.__repr__().rstrip() == "Dataprep(path='data/test', network='dummy', station=['DEBW107'], "\ - "variables=['o3', 'temp'], station_type=traffic, **None)".rstrip() + def test_download_data(self, data_prep_no_init): + file_name = data_prep_no_init._set_file_name() + meta_file = data_prep_no_init._set_meta_file_name() + data_prep_no_init.kwargs = {"store_data_locally": False} + data_prep_no_init.statistics_per_var = {'o3': 'dma8eu', 'temp': 'maximum'} + data_prep_no_init.download_data(file_name, meta_file) + assert isinstance(data_prep_no_init.data, xr.DataArray) + + def test_download_data_from_join(self, data_prep_no_init): + file_name = data_prep_no_init._set_file_name() + meta_file = data_prep_no_init._set_meta_file_name() + data_prep_no_init.kwargs = {"store_data_locally": False} + data_prep_no_init.statistics_per_var = {'o3': 'dma8eu', 'temp': 'maximum'} + xarr, meta = data_prep_no_init.download_data_from_join(file_name, meta_file) + assert isinstance(xarr, xr.DataArray) + assert isinstance(meta, pd.DataFrame) + + def test_check_station_meta(self, caplog, data_prep_no_init): + caplog.set_level(logging.DEBUG) + file_name = data_prep_no_init._set_file_name() + meta_file = data_prep_no_init._set_meta_file_name() + data_prep_no_init.kwargs = {"store_data_locally": False} + data_prep_no_init.statistics_per_var = {'o3': 'dma8eu', 'temp': 'maximum'} + data_prep_no_init.download_data(file_name, meta_file) + assert data_prep_no_init.check_station_meta() is None + data_prep_no_init.station_type = "traffic" + with pytest.raises(FileNotFoundError) as e: + data_prep_no_init.check_station_meta() + msg = "meta data does not agree with given request for station_type: traffic (requested) != background (local)" + assert caplog.record_tuples[-1][:-1] == ('root', 10) + assert msg in caplog.record_tuples[-1][-1] + + def test_load_data_overwrite_local_data(self, data_prep_no_init): + data_prep_no_init.statistics_per_var = {'o3': 'dma8eu', 'temp': 'maximum'} + file_path = data_prep_no_init._set_file_name() + meta_file_path = data_prep_no_init._set_meta_file_name() + os.remove(file_path) + os.remove(meta_file_path) + assert not os.path.exists(file_path) + assert not os.path.exists(meta_file_path) + data_prep_no_init.kwargs = {"overwrite_local_data": True} + data_prep_no_init.load_data() + assert os.path.exists(file_path) + assert os.path.exists(meta_file_path) + t = os.stat(file_path).st_ctime + tm = os.stat(meta_file_path).st_ctime + data_prep_no_init.load_data() + assert os.path.exists(file_path) + assert os.path.exists(meta_file_path) + assert os.stat(file_path).st_ctime > t + assert os.stat(meta_file_path).st_ctime > tm + assert isinstance(data_prep_no_init.data, xr.DataArray) + assert isinstance(data_prep_no_init.meta, pd.DataFrame) + + def test_load_data_keep_local_data(self, data_prep_no_init): + data_prep_no_init.statistics_per_var = {'o3': 'dma8eu', 'temp': 'maximum'} + data_prep_no_init.station_type = None + data_prep_no_init.kwargs = {} + file_path = data_prep_no_init._set_file_name() + data_prep_no_init.load_data() + assert os.path.exists(file_path) + t = os.stat(file_path).st_ctime + data_prep_no_init.load_data() + assert os.path.exists(data_prep_no_init._set_file_name()) + assert os.stat(file_path).st_ctime == t + assert isinstance(data_prep_no_init.data, xr.DataArray) + assert isinstance(data_prep_no_init.meta, pd.DataFrame) + + def test_repr(self, data_prep_no_init): + path = os.path.join(os.path.abspath(os.path.dirname(__file__)), 'data') + assert data_prep_no_init.__repr__().rstrip() == f"Dataprep(path='{path}', network='UBA', " \ + f"station=['DEBW107'], variables=['o3', 'temp'], " \ + f"station_type=background, **None)".rstrip() def test_set_file_name_and_meta(self): d = object.__new__(DataPrep) @@ -133,6 +211,16 @@ class TestDataPrep: with pytest.raises(NotImplementedError): data.inverse_transform() + def test_get_transformation_information(self, data): + assert (None, None, None) == data.get_transformation_information("o3") + mean_test = data.data.mean("datetime").sel(variables='o3').values + std_test = data.data.std("datetime").sel(variables='o3').values + data.transform('datetime') + mean, std, info = data.get_transformation_information("o3") + assert np.testing.assert_almost_equal(mean, mean_test) is None + assert np.testing.assert_almost_equal(std, std_test) is None + assert info == "standardise" + def test_nan_remove_no_hist_or_label(self, data): assert data.history is None assert data.label is None diff --git a/test/test_datastore.py b/test/test_datastore.py index 3f61c227be3a05d78f825eca77b0d6cbbc617ce1..95a58deafc915dd6193960e77bb99cc8ab8d85cb 100644 --- a/test/test_datastore.py +++ b/test/test_datastore.py @@ -2,9 +2,10 @@ __author__ = 'Lukas Leufen' __date__ = '2019-11-22' +import pytest + from src.datastore import AbstractDataStore, DataStoreByVariable, DataStoreByScope from src.datastore import NameNotFoundInDataStore, NameNotFoundInScope, EmptyScope -import pytest class TestAbstractDataStore: @@ -16,6 +17,12 @@ class TestAbstractDataStore: def test_init(self, ds): assert ds._store == {} + def test_clear_data_store(self, ds): + ds._store["test"] = "test2" + assert len(ds._store.keys()) == 1 + ds.clear_data_store() + assert len(ds._store.keys()) == 0 + class TestDataStoreByVariable: @@ -49,6 +56,12 @@ class TestDataStoreByVariable: ds.get("number3", "general") assert "Couldn't find number3 in data store" in e.value.args[0] + def test_get_default(self, ds): + ds.set("number", 3, "general") + assert ds.get_default("number", "general", 45) == 3 + assert ds.get_default("number", "general.sub", 45) == 3 + assert ds.get_default("number", "other", 45) == 45 + def test_search(self, ds): ds.set("number", 22, "general") ds.set("number", 22, "general2") @@ -118,6 +131,25 @@ class TestDataStoreByVariable: assert ds.search_scope("general.sub.sub", current_scope_only=False, return_all=True) == \ [("number", "general.sub.sub", "ABC"), ("number1", "general.sub", 22), ("number2", "general.sub.sub", 3)] + def test_create_args_dict(self, ds): + ds.set("tester1", 1, "general") + ds.set("tester2", 11, "general") + ds.set("tester2", 10, "general.sub") + ds.set("tester3", 21, "general") + args = ["tester1", "tester2", "tester3", "tester4"] + assert ds.create_args_dict(args) == {"tester1": 1, "tester2": 11, "tester3": 21} + assert ds.create_args_dict(args, "general.sub") == {"tester1": 1, "tester2": 10, "tester3": 21} + assert ds.create_args_dict(["notAvail", "alsonot"]) == {} + + def test_set_args_from_dict(self, ds): + ds.set_args_from_dict({"tester1": 1, "tester2": 10, "tester3": 21}) + assert ds.get("tester1", "general") == 1 + assert ds.get("tester2", "general") == 10 + assert ds.get("tester3", "general") == 21 + ds.set_args_from_dict({"tester1": 111}, "general.sub") + assert ds.get("tester1", "general.sub") == 111 + assert ds.get("tester3", "general.sub") == 21 + class TestDataStoreByScope: @@ -151,6 +183,12 @@ class TestDataStoreByScope: ds.get("number3", "general") assert "Couldn't find number3 in data store" in e.value.args[0] + def test_get_default(self, ds): + ds.set("number", 3, "general") + assert ds.get_default("number", "general", 45) == 3 + assert ds.get_default("number", "general.sub", 45) == 3 + assert ds.get_default("number", "other", 45) == 45 + def test_search(self, ds): ds.set("number", 22, "general") ds.set("number", 22, "general2") @@ -220,3 +258,21 @@ class TestDataStoreByScope: assert ds.search_scope("general.sub.sub", current_scope_only=False, return_all=True) == \ [("number", "general.sub.sub", "ABC"), ("number1", "general.sub", 22), ("number2", "general.sub.sub", 3)] + def test_create_args_dict(self, ds): + ds.set("tester1", 1, "general") + ds.set("tester2", 11, "general") + ds.set("tester2", 10, "general.sub") + ds.set("tester3", 21, "general") + args = ["tester1", "tester2", "tester3", "tester4"] + assert ds.create_args_dict(args) == {"tester1": 1, "tester2": 11, "tester3": 21} + assert ds.create_args_dict(args, "general.sub") == {"tester1": 1, "tester2": 10, "tester3": 21} + assert ds.create_args_dict(["notAvail", "alsonot"]) == {} + + def test_set_args_from_dict(self, ds): + ds.set_args_from_dict({"tester1": 1, "tester2": 10, "tester3": 21}) + assert ds.get("tester1", "general") == 1 + assert ds.get("tester2", "general") == 10 + assert ds.get("tester3", "general") == 21 + ds.set_args_from_dict({"tester1": 111}, "general.sub") + assert ds.get("tester1", "general.sub") == 111 + assert ds.get("tester3", "general.sub") == 21 \ No newline at end of file diff --git a/test/test_helpers.py b/test/test_helpers.py index ce5d28a63d63dc4a793e6e07c60f95cb411ae97e..b807d2b8612b9ee006bff43f1ae4cfcfd2dd07e1 100644 --- a/test/test_helpers.py +++ b/test/test_helpers.py @@ -1,11 +1,13 @@ -import pytest -from src.helpers import * import logging import os +import platform + import keras -import numpy as np import mock -import platform +import numpy as np +import pytest + +from src.helpers import * class TestToList: @@ -44,44 +46,6 @@ class TestLoss: assert hist.history['loss'][0] == 2.25 -class TestLearningRateDecay: - - def test_init(self): - lr_decay = LearningRateDecay() - assert lr_decay.lr == {'lr': []} - assert lr_decay.base_lr == 0.01 - assert lr_decay.drop == 0.96 - assert lr_decay.epochs_drop == 8 - - def test_check_param(self): - lr_decay = object.__new__(LearningRateDecay) - assert lr_decay.check_param(1, "tester") == 1 - assert lr_decay.check_param(0.5, "tester") == 0.5 - with pytest.raises(ValueError) as e: - lr_decay.check_param(0, "tester") - assert "tester is out of allowed range (0, 1]: tester=0" in e.value.args[0] - with pytest.raises(ValueError) as e: - lr_decay.check_param(1.5, "tester") - assert "tester is out of allowed range (0, 1]: tester=1.5" in e.value.args[0] - assert lr_decay.check_param(1.5, "tester", upper=None) == 1.5 - with pytest.raises(ValueError) as e: - lr_decay.check_param(0, "tester", upper=None) - assert "tester is out of allowed range (0, inf): tester=0" in e.value.args[0] - assert lr_decay.check_param(0.5, "tester", lower=None) == 0.5 - with pytest.raises(ValueError) as e: - lr_decay.check_param(0.5, "tester", lower=None, upper=0.2) - assert "tester is out of allowed range (-inf, 0.2]: tester=0.5" in e.value.args[0] - assert lr_decay.check_param(10, "tester", upper=None, lower=None) - - def test_on_epoch_begin(self): - lr_decay = LearningRateDecay(base_lr=0.02, drop=0.95, epochs_drop=2) - model = keras.Sequential() - model.add(keras.layers.Dense(1, input_dim=1)) - model.compile(optimizer=keras.optimizers.Adam(), loss=l_p_loss(2)) - model.fit(np.array([1, 0, 2, 0.5]), np.array([1, 1, 0, 0.5]), epochs=5, callbacks=[lr_decay]) - assert lr_decay.lr['lr'] == [0.02, 0.02, 0.02 * 0.95, 0.02 * 0.95, 0.02 * 0.95 * 0.95] - - class TestTimeTracking: def test_init(self): @@ -145,10 +109,18 @@ class TestTimeTracking: duration = t.stop(get_duration=True) assert duration == t.duration() + def test_enter_exit(self, caplog): + caplog.set_level(logging.INFO) + with TimeTracking() as t: + assert t.start is not None + assert t.end is None + expression = PyTestRegex(r"undefined job finished after \d+:\d+:\d+ \(hh:mm:ss\)") + assert caplog.record_tuples[-1] == ('root', 20, expression) + class TestPrepareHost: - @mock.patch("socket.gethostname", side_effect=["linux-gzsx", "ZAM144", "zam347", "jrtest", "jwtest"]) + @mock.patch("socket.gethostname", side_effect=["linux-aa9b", "ZAM144", "zam347", "jrtest", "jwtest"]) @mock.patch("os.getlogin", return_value="testUser") @mock.patch("os.path.exists", return_value=True) def test_prepare_host(self, mock_host, mock_user, mock_path): @@ -170,10 +142,10 @@ class TestPrepareHost: prepare_host() assert "unknown host 'NotExistingHostName'" in e.value.args[0] if "runner-6HmDp9Qd-project-2411-concurrent" not in platform.node(): - mock_host.return_value = "linux-gzsx" + mock_host.return_value = "linux-aa9b" with pytest.raises(NotADirectoryError) as e: prepare_host() - assert "does not exist for host 'linux-gzsx'" in e.value.args[0] + assert "does not exist for host 'linux-aa9b'" in e.value.args[0] class TestSetExperimentName: @@ -189,3 +161,55 @@ class TestSetExperimentName: def test_set_experiment_from_sys(self): exp_name, _ = set_experiment_name(experiment_date="2019-11-14") assert exp_name == "2019-11-14_network" + + +class TestPytestRegex: + + @pytest.fixture + def regex(self): + return PyTestRegex("teststring") + + def test_pytest_regex_init(self, regex): + assert regex._regex.pattern == "teststring" + + def test_pytest_regex_eq(self, regex): + assert regex == "teststringabcd" + assert regex != "teststgabcd" + + def test_pytest_regex_repr(self, regex): + assert regex.__repr__() == "teststring" + + +class TestDictToXarray: + + def test_dict_to_xarray(self): + array1 = xr.DataArray(np.random.randn(2, 3), dims=('x', 'y'), coords={'x': [10, 20]}) + array2 = xr.DataArray(np.random.randn(2, 3), dims=('x', 'y'), coords={'x': [10, 20]}) + d = {"number1": array1, "number2": array2} + res = dict_to_xarray(d, "merge_dim") + assert type(res) == xr.DataArray + assert sorted(list(res.coords)) == ["merge_dim", "x"] + assert res.shape == (2, 2, 3) + + +class TestFloatRound: + + def test_float_round_ceil(self): + assert float_round(4.6) == 5 + assert float_round(239.3992) == 240 + + def test_float_round_decimals(self): + assert float_round(23.0091, 2) == 23.01 + assert float_round(23.1091, 3) == 23.11 + + def test_float_round_type(self): + assert float_round(34.9221, 2, math.floor) == 34.92 + assert float_round(34.9221, 0, math.floor) == 34. + assert float_round(34.9221, 2, round) == 34.92 + assert float_round(34.9221, 0, round) == 35. + + def test_float_round_negative(self): + assert float_round(-34.9221, 2, math.floor) == -34.93 + assert float_round(-34.9221, 0, math.floor) == -35. + assert float_round(-34.9221, 2) == -34.92 + assert float_round(-34.9221, 0) == -34. diff --git a/test/test_join.py b/test/test_join.py index 865ae80dfaaa0244eb7592e65ef134a23b36634c..fe3d33d6296c16bfc72675bc1737aad12ee3c8b9 100644 --- a/test/test_join.py +++ b/test/test_join.py @@ -1,15 +1,18 @@ from typing import Iterable -import datetime as dt + import pytest from src.join import * from src.join import _save_to_pandas, _correct_stat_name, _lower_list +from src.join_settings import join_settings class TestJoinUrlBase: def test_url(self): - assert join_url_base == 'https://join.fz-juelich.de/services/rest/surfacedata/' + url, headers = join_settings() + assert url == 'https://join.fz-juelich.de/services/rest/surfacedata/' + assert headers == {} class TestDownloadJoin: @@ -25,22 +28,35 @@ class TestDownloadJoin: assert e.value.args[-1] == "No data found in JOIN." +class TestCorrectDataFormat: + + def test_correct_data_format(self): + list_data = [["2020-01-01 06:00:01", 23.], ["2020-01-01 06:00:11", 24.], ["2020-01-01 06:00:21", 25.], + ["2020-01-01 06:00:31", 26.], ["2020-01-01 06:00:41", 27.], ["2020-01-01 06:00:51", 23.], + {"station": "test_station_001", "author": "ME", "success": True}] + dict_data = correct_data_format(list_data) + assert dict_data == {"datetime": ["2020-01-01 06:00:01", "2020-01-01 06:00:11", "2020-01-01 06:00:21", + "2020-01-01 06:00:31", "2020-01-01 06:00:41", "2020-01-01 06:00:51"], + "values": [23., 24., 25., 26., 27., 23.], + "metadata": {"station": "test_station_001", "author": "ME", "success": True}} + + class TestGetData: def test(self): - opts = {"base": join_url_base, "service": "series", "station_id": 'DEBW107', "network_name": "UBA", + opts = {"base": join_settings()[0], "service": "series", "station_id": 'DEBW107', "network_name": "UBA", "parameter_name": "o3,no2"} - assert get_data(opts) == [[17057, 'UBA', 'DEBW107', 'O3'], [17058, 'UBA', 'DEBW107', 'NO2']] + assert get_data(opts, headers={}) == [[17057, 'UBA', 'DEBW107', 'O3'], [17058, 'UBA', 'DEBW107', 'NO2']] class TestLoadSeriesInformation: def test_standard_query(self): expected_subset = {'o3': 23031, 'no2': 39002, 'temp--lubw': 17059, 'wspeed': 17060} - assert expected_subset.items() <= load_series_information(['DEBW107'], None, None).items() + assert expected_subset.items() <= load_series_information(['DEBW107'], None, None, join_settings()[0], {}).items() def test_empty_result(self): - assert load_series_information(['DEBW107'], "traffic", None) == {} + assert load_series_information(['DEBW107'], "traffic", None, join_settings()[0], {}) == {} class TestSaveToPandas: @@ -53,6 +69,10 @@ class TestSaveToPandas: def date(self): return ['1997-01-01 00:00', '1997-01-02 00:00', '1997-01-03 00:00', '1997-01-04 00:00'] + @pytest.fixture + def date_len19(self): + return ['1997-01-01 00:00:00', '1997-01-02 00:00:00', '1997-01-03 00:00:00', '1997-01-04 00:00:00'] + @pytest.fixture def values(self): return [86.21, 94.76, 76.96, 99.89] @@ -75,6 +95,10 @@ class TestSaveToPandas: df_concat = pd.concat([create_df, next_df], axis=1) assert pd.testing.assert_frame_equal(df_concat, _save_to_pandas(create_df, data, 'max', 'temperature')) is None + def test_alternative_date_format(self, date_len19, values, create_df): + data = {'datetime': date_len19, 'mean': values, 'metadata': None} + assert pd.testing.assert_frame_equal(create_df, _save_to_pandas(None, data, 'mean', 'cloudcover')) is None + class TestCorrectStatName: diff --git a/test/test_model_modules/test_inception_model.py b/test/test_model_modules/test_inception_model.py index b2a923c47c4351ed4c2f543a4e30d25fdbaa58ea..aa5cb284ab196d733e04a9882fa4d5a4ef639a6d 100644 --- a/test/test_model_modules/test_inception_model.py +++ b/test/test_model_modules/test_inception_model.py @@ -1,6 +1,7 @@ +import keras import pytest + from src.model_modules.inception_model import InceptionModelBase -import keras class TestInceptionModelBase: diff --git a/test/test_model_modules/test_keras_extensions.py b/test/test_model_modules/test_keras_extensions.py new file mode 100644 index 0000000000000000000000000000000000000000..2f6565b4cabe295169047a6582d2b89cbf387062 --- /dev/null +++ b/test/test_model_modules/test_keras_extensions.py @@ -0,0 +1,62 @@ +import keras +import numpy as np +import pytest + +from src.helpers import l_p_loss +from src.model_modules.keras_extensions import * + + +class TestHistoryAdvanced: + + def test_init(self): + hist = HistoryAdvanced() + assert hist.validation_data is None + assert hist.model is None + assert isinstance(hist.epoch, list) and len(hist.epoch) == 0 + assert isinstance(hist.history, dict) and len(hist.history.keys()) == 0 + + def test_on_train_begin(self): + hist = HistoryAdvanced() + hist.epoch = [1, 2, 3] + hist.history = {"mse": [10, 7, 4]} + hist.on_train_begin() + assert hist.epoch == [1, 2, 3] + assert hist.history == {"mse": [10, 7, 4]} + + +class TestLearningRateDecay: + + def test_init(self): + lr_decay = LearningRateDecay() + assert lr_decay.lr == {'lr': []} + assert lr_decay.base_lr == 0.01 + assert lr_decay.drop == 0.96 + assert lr_decay.epochs_drop == 8 + + def test_check_param(self): + lr_decay = object.__new__(LearningRateDecay) + assert lr_decay.check_param(1, "tester") == 1 + assert lr_decay.check_param(0.5, "tester") == 0.5 + with pytest.raises(ValueError) as e: + lr_decay.check_param(0, "tester") + assert "tester is out of allowed range (0, 1]: tester=0" in e.value.args[0] + with pytest.raises(ValueError) as e: + lr_decay.check_param(1.5, "tester") + assert "tester is out of allowed range (0, 1]: tester=1.5" in e.value.args[0] + assert lr_decay.check_param(1.5, "tester", upper=None) == 1.5 + with pytest.raises(ValueError) as e: + lr_decay.check_param(0, "tester", upper=None) + assert "tester is out of allowed range (0, inf): tester=0" in e.value.args[0] + assert lr_decay.check_param(0.5, "tester", lower=None) == 0.5 + with pytest.raises(ValueError) as e: + lr_decay.check_param(0.5, "tester", lower=None, upper=0.2) + assert "tester is out of allowed range (-inf, 0.2]: tester=0.5" in e.value.args[0] + assert lr_decay.check_param(10, "tester", upper=None, lower=None) + + def test_on_epoch_begin(self): + lr_decay = LearningRateDecay(base_lr=0.02, drop=0.95, epochs_drop=2) + model = keras.Sequential() + model.add(keras.layers.Dense(1, input_dim=1)) + model.compile(optimizer=keras.optimizers.Adam(), loss=l_p_loss(2)) + model.fit(np.array([1, 0, 2, 0.5]), np.array([1, 1, 0, 0.5]), epochs=5, callbacks=[lr_decay]) + assert lr_decay.lr['lr'] == [0.02, 0.02, 0.02 * 0.95, 0.02 * 0.95, 0.02 * 0.95 * 0.95] diff --git a/test/test_model_modules/test_model_class.py b/test/test_model_modules/test_model_class.py index 0af16012336c9dbcc3133d7dac1f365e276d11bc..0dbd2d9b67a0748bf09eb4f59e1888aae1ea405d 100644 --- a/test/test_model_modules/test_model_class.py +++ b/test/test_model_modules/test_model_class.py @@ -1,5 +1,5 @@ -import pytest import keras +import pytest from src.model_modules.model_class import AbstractModelClass diff --git a/test/test_modules/test_experiment_setup.py b/test/test_modules/test_experiment_setup.py index 4a2ba9f9e641612a7f9ffbc81b93b68fd5dbccd8..7a4f16fd055e0af0aea95181d475d061c185ca92 100644 --- a/test/test_modules/test_experiment_setup.py +++ b/test/test_modules/test_experiment_setup.py @@ -1,11 +1,11 @@ -import pytest -import logging import argparse +import logging import os -from src.run_modules.experiment_setup import ExperimentSetup +import pytest + from src.helpers import TimeTracking, prepare_host -from src.datastore import NameNotFoundInScope, NameNotFoundInDataStore +from src.run_modules.experiment_setup import ExperimentSetup class TestExperimentSetup: @@ -149,3 +149,14 @@ class TestExperimentSetup: assert data_store.get("end", "general.test") == "2000-01-06" # use all stations on all data sets (train, val, test) assert data_store.get("use_all_stations_on_all_data_sets", "general.test") is False + + def test_compare_variables_and_statistics(self): + kwargs = dict(parser_args={"experiment_date": "TODAY"}, + var_all_dict={'o3': 'dma8eu', 'temp': 'maximum'}, + stations=['DEBY053', 'DEBW059', 'DEBW027'], variables=["o3", "relhum"], statistics_per_var=None) + with pytest.raises(ValueError) as e: + ExperimentSetup(**kwargs) + assert "for the variables: {'relhum'}" in e.value.args[0] + + kwargs["variables"] = ["o3", "temp"] + assert ExperimentSetup(**kwargs) is not None diff --git a/test/test_modules/test_model_setup.py b/test/test_modules/test_model_setup.py index d604d7474af84740a4b7a1cc51e5e94f1c94533b..2864ae45bcd7d3c6109d6d84fe5ea152a7d86384 100644 --- a/test/test_modules/test_model_setup.py +++ b/test/test_modules/test_model_setup.py @@ -1,13 +1,12 @@ -import pytest import os -import keras -import mock -from src.run_modules.model_setup import ModelSetup -from src.run_modules.run_environment import RunEnvironment +import pytest + from src.data_handling.data_generator import DataGenerator -from src.model_modules.model_class import AbstractModelClass from src.datastore import EmptyScope +from src.model_modules.model_class import AbstractModelClass +from src.run_modules.model_setup import ModelSetup +from src.run_modules.run_environment import RunEnvironment class TestModelSetup: @@ -18,6 +17,9 @@ class TestModelSetup: super(ModelSetup, obj).__init__() obj.scope = "general.modeltest" obj.model = None + obj.callbacks_name = "placeholder_%s_str.pickle" + obj.data_store.set("lr_decay", "dummy_str", "general.model") + obj.data_store.set("hist", "dummy_str", "general.model") yield obj RunEnvironment().__del__() diff --git a/test/test_modules/test_post_processing.py b/test/test_modules/test_post_processing.py index e6a245ccb27807d59114d0a7c69e472d98bc1e58..67897451eaa5de065b644d1f9868a846cb57d84e 100644 --- a/test/test_modules/test_post_processing.py +++ b/test/test_modules/test_post_processing.py @@ -1,8 +1,3 @@ -import keras - -from src.run_modules.post_processing import PostProcessing - - class TestPostProcessing: def test_init(self): diff --git a/test/test_modules/test_pre_processing.py b/test/test_modules/test_pre_processing.py index a562e7b05a79f0068b10e9e36771669fe47d4ce8..9d1feb03ac71b980f6a4cd1b0e6cac2a52d9625b 100644 --- a/test/test_modules/test_pre_processing.py +++ b/test/test_modules/test_pre_processing.py @@ -1,11 +1,12 @@ import logging + import pytest +from src.data_handling.data_generator import DataGenerator +from src.datastore import NameNotFoundInScope from src.helpers import PyTestRegex from src.run_modules.experiment_setup import ExperimentSetup from src.run_modules.pre_processing import PreProcessing, DEFAULT_ARGS_LIST, DEFAULT_KWARGS_LIST -from src.data_handling.data_generator import DataGenerator -from src.datastore import NameNotFoundInScope from src.run_modules.run_environment import RunEnvironment @@ -47,14 +48,15 @@ class TestPreProcessing: assert obj_with_exp_setup.data_store.search_name("generator") == [] assert obj_with_exp_setup._run() is None assert obj_with_exp_setup.data_store.search_name("generator") == sorted(["general.train", "general.val", - "general.test"]) + "general.train_val", "general.test"]) def test_split_train_val_test(self, obj_with_exp_setup): assert obj_with_exp_setup.data_store.search_name("generator") == [] obj_with_exp_setup.split_train_val_test() data_store = obj_with_exp_setup.data_store assert data_store.search_scope("general.train") == sorted(["generator", "start", "end", "stations"]) - assert data_store.search_name("generator") == sorted(["general.train", "general.val", "general.test"]) + assert data_store.search_name("generator") == sorted(["general.train", "general.val", "general.test", + "general.train_val"]) def test_create_set_split_not_all_stations(self, caplog, obj_with_exp_setup): caplog.set_level(logging.DEBUG) @@ -93,10 +95,11 @@ class TestPreProcessing: def test_split_set_indices(self, obj_super_init): dummy_list = list(range(0, 15)) - train, val, test = obj_super_init.split_set_indices(len(dummy_list), 0.9) + train, val, test, train_val = obj_super_init.split_set_indices(len(dummy_list), 0.9) assert dummy_list[train] == list(range(0, 10)) assert dummy_list[val] == list(range(10, 13)) assert dummy_list[test] == list(range(13, 15)) + assert dummy_list[train_val] == list(range(0, 13)) def test_create_args_dict_default_scope(self, obj_super_init): assert obj_super_init.data_store.create_args_dict(["NAME1", "NAME2"]) == {"NAME1": 1, "NAME2": 2} diff --git a/test/test_modules/test_training.py b/test/test_modules/test_training.py index 3426eb1d355a6690ee57c3ee45e5088d7df9c249..485348ceca740d8263394fca36efbfbde6dd2d0d 100644 --- a/test/test_modules/test_training.py +++ b/test/test_modules/test_training.py @@ -1,20 +1,22 @@ +import glob +import json +import logging +import os +import shutil + import keras +import mock import pytest from keras.callbacks import ModelCheckpoint, History -import mock -import os -import json -import shutil -import logging -import glob -from src.model_modules.inception_model import InceptionModelBase -from src.model_modules.flatten import flatten_tail -from src.run_modules.training import Training -from src.run_modules.run_environment import RunEnvironment from src.data_handling.data_distributor import Distributor from src.data_handling.data_generator import DataGenerator -from src.helpers import LearningRateDecay, PyTestRegex +from src.helpers import PyTestRegex +from src.model_modules.flatten import flatten_tail +from src.model_modules.inception_model import InceptionModelBase +from src.model_modules.keras_extensions import LearningRateDecay, HistoryAdvanced +from src.run_modules.run_environment import RunEnvironment +from src.run_modules.training import Training def my_test_model(activation, window_history_size, channels, dropout_rate, add_minor_branch=False): @@ -37,7 +39,7 @@ def my_test_model(activation, window_history_size, channels, dropout_rate, add_m class TestTraining: @pytest.fixture - def init_without_run(self, path, model, checkpoint): + def init_without_run(self, path: str, model: keras.Model, checkpoint: ModelCheckpoint): obj = object.__new__(Training) super(Training, obj).__init__() obj.model = model @@ -48,6 +50,7 @@ class TestTraining: obj.epochs = 2 obj.checkpoint = checkpoint obj.lr_sc = LearningRateDecay() + obj.hist = HistoryAdvanced() obj.experiment_name = "TestExperiment" obj.data_store.set("generator", mock.MagicMock(return_value="mock_train_gen"), "general.train") obj.data_store.set("generator", mock.MagicMock(return_value="mock_val_gen"), "general.val") @@ -82,6 +85,7 @@ class TestTraining: 'loss': [0.6795708956961347, 0.45963566494176616], 'mean_squared_error': [0.6795708956961347, 0.45963566494176616], 'mean_absolute_error': [0.6523177288928538, 0.5363963260296364]} + h.model = mock.MagicMock() return h @pytest.fixture @@ -103,7 +107,7 @@ class TestTraining: return ModelCheckpoint(os.path.join(path, "model_checkpoint"), monitor='val_loss', save_best_only=True) @pytest.fixture - def ready_to_train(self, generator, init_without_run): + def ready_to_train(self, generator: DataGenerator, init_without_run: Training): init_without_run.train_set = Distributor(generator, init_without_run.model, init_without_run.batch_size) init_without_run.val_set = Distributor(generator, init_without_run.model, init_without_run.batch_size) init_without_run.model.compile(optimizer=keras.optimizers.SGD(), loss=keras.losses.mean_absolute_error) @@ -131,6 +135,7 @@ class TestTraining: obj.data_store.set("epochs", 2, "general.model") obj.data_store.set("checkpoint", checkpoint, "general.model") obj.data_store.set("lr_decay", LearningRateDecay(), "general.model") + obj.data_store.set("hist", HistoryAdvanced(), "general.model") obj.data_store.set("experiment_name", "TestExperiment", "general") obj.data_store.set("experiment_path", path, "general") path_plot = os.path.join(path, "plots") @@ -187,26 +192,28 @@ class TestTraining: assert caplog.record_tuples[1] == ("root", 20, PyTestRegex("no weights to reload...")) def test_save_callbacks_history_created(self, init_without_run, history, path): - init_without_run.save_callbacks(history) + init_without_run.save_callbacks_as_json(history) assert "history.json" in os.listdir(path) def test_save_callbacks_lr_created(self, init_with_lr, history, path): - init_with_lr.save_callbacks(history) + init_with_lr.save_callbacks_as_json(history) assert "history_lr.json" in os.listdir(path) def test_save_callbacks_inspect_history(self, init_without_run, history, path): - init_without_run.save_callbacks(history) + init_without_run.save_callbacks_as_json(history) with open(os.path.join(path, "history.json")) as jfile: hist = json.load(jfile) assert hist == history.history def test_save_callbacks_inspect_lr(self, init_with_lr, history, path): - init_with_lr.save_callbacks(history) + init_with_lr.save_callbacks_as_json(history) with open(os.path.join(path, "history_lr.json")) as jfile: lr = json.load(jfile) assert lr == init_with_lr.lr_sc.lr def test_create_monitoring_plots(self, init_without_run, learning_rate, history, path): assert len(glob.glob(os.path.join(path, "plots", "TestExperiment_history_*.pdf"))) == 0 + history.model.output_names = mock.MagicMock(return_value=["Main"]) + history.model.metrics_names = mock.MagicMock(return_value=["loss", "mean_squared_error"]) init_without_run.create_monitoring_plots(history, learning_rate) assert len(glob.glob(os.path.join(path, "plots", "TestExperiment_history_*.pdf"))) == 2 diff --git a/test/test_plotting/test_training_monitoring.py b/test/test_plotting/test_training_monitoring.py index d38fc623b7fd2f5f8c779ea03757d7bf0794089a..7e4e21c1a28b35bef4aa6e613756378fe41611b5 100644 --- a/test/test_plotting/test_training_monitoring.py +++ b/test/test_plotting/test_training_monitoring.py @@ -1,9 +1,10 @@ +import os + import keras import pytest -import os +from src.model_modules.keras_extensions import LearningRateDecay from src.plotting.training_monitoring import PlotModelLearningRate, PlotModelHistory -from src.helpers import LearningRateDecay @pytest.fixture @@ -17,50 +18,81 @@ def path(): class TestPlotModelHistory: @pytest.fixture - def history(self): + def default_history(self): hist = keras.callbacks.History() hist.epoch = [0, 1] hist.history = {'val_loss': [0.5586272982587484, 0.45712877659670287], - 'val_mean_squared_error': [0.5586272982587484, 0.45712877659670287], - 'val_mean_absolute_error': [0.595368885413389, 0.530547587585537], - 'loss': [0.6795708956961347, 0.45963566494176616], - 'mean_squared_error': [0.6795708956961347, 0.45963566494176616], - 'mean_absolute_error': [0.6523177288928538, 0.5363963260296364]} + 'val_mean_squared_error': [0.5586272982587484, 0.45712877659670287], + 'val_mean_absolute_error': [0.595368885413389, 0.530547587585537], + 'loss': [0.6795708956961347, 0.45963566494176616], + 'mean_squared_error': [0.6795708956961347, 0.45963566494176616], + 'mean_absolute_error': [0.6523177288928538, 0.5363963260296364]} return hist @pytest.fixture - def history_var(self): + def history_additional_loss(self): hist = keras.callbacks.History() hist.epoch = [0, 1] hist.history = {'val_loss': [0.5586272982587484, 0.45712877659670287], - 'test_loss': [0.595368885413389, 0.530547587585537], - 'loss': [0.6795708956961347, 0.45963566494176616], - 'mean_squared_error': [0.6795708956961347, 0.45963566494176616], - 'mean_absolute_error': [0.6523177288928538, 0.5363963260296364]} + 'test_loss': [0.595368885413389, 0.530547587585537], + 'loss': [0.6795708956961347, 0.45963566494176616], + 'mean_squared_error': [0.6795708956961347, 0.45963566494176616], + 'mean_absolute_error': [0.6523177288928538, 0.5363963260296364]} return hist + @pytest.fixture + def history_with_main(self, default_history): + default_history.history["main_val_loss"] = [0.5586272982587484, 0.45712877659670287] + default_history.history["main_loss"] = [0.6795708956961347, 0.45963566494176616] + return default_history + @pytest.fixture def no_init(self): return object.__new__(PlotModelHistory) - def test_plot_from_hist_obj(self, history, path): + def test_get_plot_metric(self, no_init, default_history): + history = default_history.history + metric = no_init._get_plot_metric(history, plot_metric="loss", main_branch=False) + assert metric == "loss" + metric = no_init._get_plot_metric(history, plot_metric="mean_squared_error", main_branch=False) + assert metric == "mean_squared_error" + + def test_get_plot_metric_short_metric(self, no_init, default_history): + history = default_history.history + metric = no_init._get_plot_metric(history, plot_metric="mse", main_branch=False) + assert metric == "mean_squared_error" + metric = no_init._get_plot_metric(history, plot_metric="mae", main_branch=False) + assert metric == "mean_absolute_error" + + def test_get_plot_metric_main_branch(self, no_init, history_with_main): + history = history_with_main.history + metric = no_init._get_plot_metric(history, plot_metric="loss", main_branch=True) + assert metric == "main_loss" + + def test_filter_columns(self, no_init): + no_init._plot_metric = "loss" + res = no_init._filter_columns({'loss': None, 'another_loss': None, 'val_loss': None, 'wrong': None}) + assert res == ['another_loss'] + no_init._plot_metric = "mean_squared_error" + res = no_init._filter_columns({'mean_squared_error': None, 'another_loss': None, 'val_mean_squared_error': None, + 'wrong': None}) + assert res == [] + + def test_plot_from_hist_obj(self, default_history, path): assert "hist_obj.pdf" not in os.listdir(path) - PlotModelHistory(os.path.join(path, "hist_obj.pdf"), history) + PlotModelHistory(os.path.join(path, "hist_obj.pdf"), default_history) assert "hist_obj.pdf" in os.listdir(path) - def test_plot_from_hist_dict(self, history, path): + def test_plot_from_hist_dict(self, default_history, path): assert "hist_dict.pdf" not in os.listdir(path) - PlotModelHistory(os.path.join(path, "hist_dict.pdf"), history.history) + PlotModelHistory(os.path.join(path, "hist_dict.pdf"), default_history.history) assert "hist_dict.pdf" in os.listdir(path) - def test_plot_additional_loss(self, history_var, path): + def test_plot_additional_loss(self, history_additional_loss, path): assert "hist_additional.pdf" not in os.listdir(path) - PlotModelHistory(os.path.join(path, "hist_additional.pdf"), history_var) + PlotModelHistory(os.path.join(path, "hist_additional.pdf"), history_additional_loss) assert "hist_additional.pdf" in os.listdir(path) - def test_filter_list(self, no_init): - res = no_init._filter_columns({'loss': None, 'another_loss': None, 'val_loss': None, 'wrong': None}) - assert res == ['another_loss'] class TestPlotModelLearningRate: diff --git a/test/test_statistics.py b/test/test_statistics.py index fa7de6022f0957031f90fbf951679594583eccab..308ac655787e69f90b45e65e7e7df8f35875f652 100644 --- a/test/test_statistics.py +++ b/test/test_statistics.py @@ -1,7 +1,8 @@ +import numpy as np +import pandas as pd import pytest import xarray as xr -import pandas as pd -import numpy as np + from src.statistics import standardise, standardise_inverse, centre, centre_inverse