diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index f0da24d54dd183d0a4d40667ea0a9619b94f1e6e..9f1ba73f5a3b3fada8624dfd89fa6f12488a877b 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -32,7 +32,7 @@ tests:
     - django
   stage: test
   variables:
-    FAILURE_THRESHOLD: 90
+    FAILURE_THRESHOLD: 100
   before_script:
     - chmod +x ./CI/update_badge.sh
     - ./CI/update_badge.sh > /dev/null
diff --git a/requirements.txt b/requirements.txt
index 227da2b61976d6147902e05a54a3e414dc3f40cc..607563fd9042d68805c2be238ba982d62ba28974 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,26 +1,63 @@
+absl-py==0.9.0
+astor==0.8.1
+atomicwrites==1.3.0
+attrs==19.3.0
+Cartopy==0.17.0
+certifi==2019.11.28
+chardet==3.0.4
+cloudpickle==1.3.0
+coverage==5.0.3
+cycler==0.10.0
+Cython==0.29.15
+dask==2.11.0
+fsspec==0.6.2
+gast==0.3.3
+grpcio==1.27.2
+h5py==2.10.0
+idna==2.8
+importlib-metadata==1.5.0
 Keras==2.2.4
-numpy==1.15.4
+Keras-Applications==1.0.8
+Keras-Preprocessing==1.1.0
+kiwisolver==1.1.0
+locket==0.2.0
+Markdown==3.2.1
+matplotlib==3.2.0
+mock==4.0.1
+more-itertools==8.2.0
+numpy==1.18.1
+packaging==20.3
+pandas==1.0.1
+partd==1.1.0
+patsy==0.5.1
+Pillow==7.0.0
+pluggy==0.13.1
+protobuf==3.11.3
+py==1.8.1
+pydot==1.4.1
+pyparsing==2.4.6
+pyproj==2.5.0
+pyshp==2.1.0
+pytest==5.3.5
+pytest-cov==2.8.1
+pytest-html==2.0.1
+pytest-lazy-fixture==0.6.3
+pytest-metadata==1.8.0
+python-dateutil==2.8.1
+pytz==2019.3
+PyYAML==5.3
+requests==2.23.0
+scipy==1.4.1
+seaborn==0.10.0
+Shapely==1.7.0
+six==1.11.0
+statsmodels==0.11.1
+tensorboard==1.12.2
 tensorflow==1.12.0
-xarray==0.14.0
-pandas==0.25.1
-requests==2.22.0
-pytest==5.2.1
-pytest-lazy-fixture==0.6.1
-pytest-cov
-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
+termcolor==1.1.0
+toolz==0.10.0
+urllib3==1.25.8
+wcwidth==0.1.8
+Werkzeug==1.0.0
+xarray==0.15.0
+zipp==3.1.0
diff --git a/run.py b/run.py
index 8e4d9c46f4a39224335dd65b689f519943166d0f..556cd0b59ed023178fa7e6df1b5b03b9f6c5f157 100644
--- a/run.py
+++ b/run.py
@@ -17,7 +17,7 @@ def main(parser_args):
 
     with RunEnvironment():
         ExperimentSetup(parser_args, stations=['DEBW107', 'DEBY081', 'DEBW013', 'DEBW076', 'DEBW087', 'DEBW001'],
-                        station_type='background', trainable=True)
+                        station_type='background', trainable=False, create_new_model=False)
         PreProcessing()
 
         ModelSetup()
diff --git a/src/data_handling/bootstraps.py b/src/data_handling/bootstraps.py
new file mode 100644
index 0000000000000000000000000000000000000000..3e69950267f9c95ccc636e560a21731ade388432
--- /dev/null
+++ b/src/data_handling/bootstraps.py
@@ -0,0 +1,205 @@
+__author__ = 'Felix Kleinert, Lukas Leufen'
+__date__ = '2020-02-07'
+
+
+from src.run_modules.run_environment import RunEnvironment
+from src.data_handling.data_generator import DataGenerator
+import numpy as np
+import logging
+import dask.array as da
+import xarray as xr
+import os
+import re
+from src import helpers
+
+
+class BootStrapGenerator:
+
+    def __init__(self, orig_generator, boots, chunksize, bootstrap_path):
+        self.orig_generator: DataGenerator = orig_generator
+        self.stations = self.orig_generator.stations
+        self.variables = self.orig_generator.variables
+        self.boots = boots
+        self.chunksize = chunksize
+        self.bootstrap_path = bootstrap_path
+        self._iterator = 0
+        self.bootstrap_meta = []
+
+    def __len__(self):
+        """
+        display the number of stations
+        """
+        return len(self.orig_generator)*self.boots*len(self.variables)
+
+    def get_labels(self, key):
+        _, label = self.orig_generator[key]
+        for _ in range(self.boots):
+            yield label
+
+    def get_generator(self):
+        """
+        This is the implementation of the __next__ method of the iterator protocol. Get the data generator, and return
+        the history and label data of this generator.
+        :return:
+        """
+        while True:
+            for i, data in enumerate(self.orig_generator):
+                station = self.orig_generator.get_station_key(i)
+                logging.info(f"station: {station}")
+                hist, label = data
+                len_of_label = len(label)
+                shuffled_data = self.load_boot_data(station)
+                for var in self.variables:
+                    logging.info(f"  var: {var}")
+                    for boot in range(self.boots):
+                        logging.debug(f"boot: {boot}")
+                        boot_hist = hist.sel(variables=helpers.list_pop(self.variables, var))
+                        shuffled_var = shuffled_data.sel(variables=var, boots=boot).expand_dims("variables").drop("boots").transpose("datetime", "window", "Stations", "variables")
+                        boot_hist = boot_hist.combine_first(shuffled_var)
+                        boot_hist = boot_hist.sortby("variables")
+                        self.bootstrap_meta.extend([[var, station]]*len_of_label)
+                        yield boot_hist, label
+            return
+
+    def get_orig_prediction(self, path, file_name, prediction_name="CNN"):
+        file = os.path.join(path, file_name)
+        data = xr.open_dataarray(file)
+        for _ in range(self.boots):
+            yield data.sel(type=prediction_name).squeeze()
+
+    def load_boot_data(self, station):
+        files = os.listdir(self.bootstrap_path)
+        regex = re.compile(rf"{station}_\w*\.nc")
+        file_name = os.path.join(self.bootstrap_path, list(filter(regex.search, files))[0])
+        shuffled_data = xr.open_dataarray(file_name, chunks=100)
+        return shuffled_data
+
+
+class BootStraps(RunEnvironment):
+
+    def __init__(self, data, bootstrap_path, number_bootstraps=10):
+
+        super().__init__()
+        self.data: DataGenerator = data
+        self.number_bootstraps = number_bootstraps
+        self.bootstrap_path = bootstrap_path
+        self.chunks = self.get_chunk_size()
+        self.create_shuffled_data()
+        self._boot_strap_generator = BootStrapGenerator(self.data, self.number_bootstraps, self.chunks, self.bootstrap_path)
+
+    def get_boot_strap_meta(self):
+        return self._boot_strap_generator.bootstrap_meta
+
+    def boot_strap_generator(self):
+        return self._boot_strap_generator.get_generator()
+
+    def get_boot_strap_generator_length(self):
+        return self._boot_strap_generator.__len__()
+
+    def get_labels(self, key):
+        labels_list = []
+        chunks = None
+        for labels in self._boot_strap_generator.get_labels(key):
+            if len(labels_list) == 0:
+                chunks = (100, labels.data.shape[1])
+            labels_list.append(da.from_array(labels.data, chunks=chunks))
+        labels_out = da.concatenate(labels_list, axis=0)
+        return labels_out.compute()
+
+    def get_orig_prediction(self, path, name):
+        labels_list = []
+        chunks = None
+        for labels in self._boot_strap_generator.get_orig_prediction(path, name):
+            if len(labels_list) == 0:
+                chunks = (100, labels.data.shape[1])
+            labels_list.append(da.from_array(labels.data, chunks=chunks))
+        labels_out = da.concatenate(labels_list, axis=0)
+        labels_out = labels_out.compute()
+        return labels_out[~np.isnan(labels_out).any(axis=1), :]
+
+    def get_chunk_size(self):
+        hist, _ = self.data[0]
+        return (100, *hist.shape[1:], self.number_bootstraps)
+
+    def create_shuffled_data(self):
+        """
+        Create shuffled data. Use original test data, add dimension 'boots' with length number of bootstraps and insert
+        randomly selected variables. If there is a suitable local file for requested window size and number of
+        bootstraps, no additional file will be created inside this function.
+        """
+        logging.info("create shuffled bootstrap data")
+        variables_str = '_'.join(sorted(self.data.variables))
+        window = self.data.window_history_size
+        for station in self.data.stations:
+            valid, nboot = self.valid_bootstrap_file(station, variables_str, window)
+            if not valid:
+                logging.info(f'create bootstap data for {station}')
+                hist, _ = self.data[station]
+                data = hist.copy()
+                file_name = f"{station}_{variables_str}_hist{window}_nboots{nboot}_shuffled.nc"
+                file_path = os.path.join(self.bootstrap_path, file_name)
+                data = data.expand_dims({'boots': range(nboot)}, axis=-1)
+                shuffled_variable = []
+                for i, var in enumerate(data.coords['variables']):
+                    single_variable = data.sel(variables=var).values
+                    shuffled_variable.append(self.shuffle_single_variable(single_variable, chunks=(100, *data.shape[1:3], data.shape[-1])))
+                shuffled_variable_da = da.stack(shuffled_variable, axis=-2, ).rechunk("auto")
+                shuffled_data = xr.DataArray(shuffled_variable_da, coords=data.coords, dims=data.dims)
+                shuffled_data.to_netcdf(file_path)
+
+    def valid_bootstrap_file(self, station, variables, window):
+        """
+        Compare local bootstrap file with given settings for station, variables, window and number of bootstraps. If a
+        match was found, this method returns a tuple (True, None). In any other case, it returns (False, max_nboot),
+        where max_nboot is the highest boot number found in the local storage. A match is defined so that the window
+        length is ge than given window size form args and the number of boots is also ge than the given number of boots
+        from this class. Furthermore, this functions deletes local files, if the match the station pattern but don't fit
+        the window and bootstrap condition. This is performed, because it is assumed, that the corresponding file will
+        be created with a longer or at least same window size and numbers of bootstraps.
+        :param station:
+        :param variables:
+        :param window:
+        :return:
+        """
+        regex = re.compile(rf"{station}_{variables}_hist(\d+)_nboots(\d+)_shuffled")
+        max_nboot = self.number_bootstraps
+        for file in os.listdir(self.bootstrap_path):
+            match = regex.match(file)
+            if match:
+                window_file = int(match.group(1))
+                nboot_file = int(match.group(2))
+                max_nboot = max([max_nboot, nboot_file])
+                if (window_file >= window) and (nboot_file >= self.number_bootstraps):
+                    return True, None
+                else:
+                    os.remove(os.path.join(self.bootstrap_path, file))
+        return False, max_nboot
+
+    @staticmethod
+    def shuffle_single_variable(data: da.array, chunks) -> da.core.Array:
+        size = data.shape
+        return da.random.choice(data.reshape(-1,), size=size, chunks=chunks)
+
+
+if __name__ == "__main__":
+
+    from src.run_modules.experiment_setup import ExperimentSetup
+    from src.run_modules.run_environment import RunEnvironment
+    from src.run_modules.pre_processing import PreProcessing
+
+    formatter = '%(asctime)s - %(levelname)s: %(message)s  [%(filename)s:%(funcName)s:%(lineno)s]'
+    logging.basicConfig(format=formatter, level=logging.INFO)
+
+    with RunEnvironment() as run_env:
+        ExperimentSetup(stations=['DEBW107', 'DEBY081', 'DEBW013'],
+                        station_type='background', trainable=True, window_history_size=9)
+        PreProcessing()
+
+        data = run_env.data_store.get("generator", "general.test")
+        path = run_env.data_store.get("bootstrap_path", "general")
+        number_bootstraps = 10
+
+        boots = BootStraps(data, path, number_bootstraps)
+        for b in boots.boot_strap_generator():
+            a, c = b
+        logging.info(f"len is {len(boots.get_boot_strap_meta())}")
diff --git a/src/data_handling/data_generator.py b/src/data_handling/data_generator.py
index 79e1e7e72c1779d18a11652ab132c253e1dff806..24c9ada65b4bfd71de12785b2714cc5de94dc21f 100644
--- a/src/data_handling/data_generator.py
+++ b/src/data_handling/data_generator.py
@@ -80,8 +80,7 @@ class DataGenerator(keras.utils.Sequence):
             data = self.get_data_generator()
             self._iterator += 1
             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")
+                return data.get_transposed_history(), data.get_transposed_label()
             else:
                 self.__next__()  # pragma: no cover
         else:
@@ -94,7 +93,7 @@ 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.get_transposed_history(), data.label.squeeze("Stations").transpose("datetime", "window")
+        return data.get_transposed_history(), data.get_transposed_label()
 
     def setup_transformation(self, transformation):
         if transformation is None:
@@ -182,12 +181,13 @@ class DataGenerator(keras.utils.Sequence):
             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)
             if self.transformation is not None:
                 data.transform("datetime", **helpers.dict_pop(self.transformation, "scope"))
-            data.make_history_window(self.interpolate_dim, self.window_history_size)
+            data.interpolate(self.interpolate_dim, method=self.interpolate_method, limit=self.limit_nan_fill)
+            data.make_history_window(self.target_dim, self.window_history_size, self.interpolate_dim)
             data.make_labels(self.target_dim, self.target_var, self.interpolate_dim, self.window_lead_time)
-            data.history_label_nan_remove(self.interpolate_dim)
+            data.make_observation(self.target_dim, self.target_var, self.interpolate_dim)
+            data.remove_nan(self.interpolate_dim)
             if save_local_tmp_storage:
                 self._save_pickle_data(data)
         return data
diff --git a/src/data_handling/data_preparation.py b/src/data_handling/data_preparation.py
index 98b47a6df3825581564fa9aaef7be8698408760e..3fae09306ab65d18f19d770b525cdc2296215bcd 100644
--- a/src/data_handling/data_preparation.py
+++ b/src/data_handling/data_preparation.py
@@ -2,6 +2,7 @@ __author__ = 'Felix Kleinert, Lukas Leufen'
 __date__ = '2019-10-16'
 
 import datetime as dt
+from functools import reduce
 import logging
 import os
 from typing import Union, List, Iterable
@@ -15,6 +16,7 @@ from src import statistics
 
 # define a more general date type for type hinting
 date = Union[dt.date, dt.datetime]
+str_or_list = Union[str, List[str]]
 
 
 class DataPrep(object):
@@ -55,6 +57,7 @@ class DataPrep(object):
         self.std = None
         self.history = None
         self.label = None
+        self.observation = None
         self.kwargs = kwargs
         self.data = None
         self.meta = None
@@ -135,10 +138,12 @@ class DataPrep(object):
         return xarr, meta
 
     def _set_file_name(self):
-        return os.path.join(self.path, f"{''.join(self.station)}_{'_'.join(sorted(self.variables))}.nc")
+        all_vars = sorted(self.statistics_per_var.keys())
+        return os.path.join(self.path, f"{''.join(self.station)}_{'_'.join(all_vars)}.nc")
 
     def _set_meta_file_name(self):
-        return os.path.join(self.path, f"{''.join(self.station)}_{'_'.join(sorted(self.variables))}_meta.csv")
+        all_vars = sorted(self.statistics_per_var.keys())
+        return os.path.join(self.path, f"{''.join(self.station)}_{'_'.join(all_vars)}_meta.csv")
 
     def __repr__(self):
         return f"Dataprep(path='{self.path}', network='{self.network}', station={self.station}, " \
@@ -275,19 +280,20 @@ class DataPrep(object):
             std = None
         return mean, std, self._transform_method
 
-    def make_history_window(self, dim: str, window: int) -> None:
+    def make_history_window(self, dim_name_of_inputs: str, window: int, dim_name_of_shift: str) -> None:
         """
         This function uses shifts the data window+1 times and returns a xarray which has a new dimension 'window'
         containing the shifted data. This is used to represent history in the data. Results are stored in self.history .
 
-        :param dim: Dimension along shift will be applied
+        :param dim_name_of_inputs: Name of dimension which contains the input variables
         :param window: number of time steps to look back in history
                 Note: window will be treated as negative value. This should be in agreement with looking back on
                 a time line. Nonetheless positive values are allowed but they are converted to its negative
                 expression
+        :param dim_name_of_shift: Dimension along shift will be applied
         """
         window = -abs(window)
-        self.history = self.shift(dim, window)
+        self.history = self.shift(dim_name_of_shift, window).sel({dim_name_of_inputs: self.variables})
 
     def shift(self, dim: str, window: int) -> xr.DataArray:
         """
@@ -310,7 +316,7 @@ class DataPrep(object):
         res = xr.concat(res, dim=window_array)
         return res
 
-    def make_labels(self, dim_name_of_target: str, target_var: str, dim_name_of_shift: str, window: int) -> None:
+    def make_labels(self, dim_name_of_target: str, target_var: str_or_list, dim_name_of_shift: str, window: int) -> None:
         """
         This function creates a xarray.DataArray containing labels
 
@@ -322,7 +328,17 @@ class DataPrep(object):
         window = abs(window)
         self.label = self.shift(dim_name_of_shift, window).sel({dim_name_of_target: target_var})
 
-    def history_label_nan_remove(self, dim: str) -> None:
+    def make_observation(self, dim_name_of_target: str, target_var: str_or_list, dim_name_of_shift: str) -> None:
+        """
+        This function creates a xarray.DataArray containing labels
+
+        :param dim_name_of_target: Name of dimension which contains the target variable
+        :param target_var: Name of target variable(s) in 'dimension'
+        :param dim_name_of_shift: Name of dimension on which xarray.DataArray.shift will be applied
+        """
+        self.observation = self.shift(dim_name_of_shift, 0).sel({dim_name_of_target: target_var})
+
+    def remove_nan(self, dim: str) -> None:
         """
         All NAs slices in dim which contain nans in self.history or self.label are removed in both data sets.
         This is done to present only a full matrix to keras.fit.
@@ -334,14 +350,17 @@ 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)
+            non_nan_observation = self.observation.dropna(dim=dim)
+            intersect = reduce(np.intersect1d, (non_nan_history.coords[dim].values, non_nan_label.coords[dim].values, non_nan_observation.coords[dim].values))
 
         if len(intersect) == 0:
             self.history = None
             self.label = None
+            self.observation = None
         else:
             self.history = self.history.sel({dim: intersect})
             self.label = self.label.sel({dim: intersect})
+            self.observation = self.observation.sel({dim: intersect})
 
     @staticmethod
     def create_index_array(index_name: str, index_value: Iterable[int]) -> xr.DataArray:
@@ -395,7 +414,10 @@ class DataPrep(object):
         return data
 
     def get_transposed_history(self):
-        return self.history.transpose("datetime", "window", "Stations", "variables")
+        return self.history.transpose("datetime", "window", "Stations", "variables").copy()
+
+    def get_transposed_label(self):
+        return self.label.squeeze("Stations").transpose("datetime", "window").copy()
 
 
 if __name__ == "__main__":
diff --git a/src/helpers.py b/src/helpers.py
index e1496c3b232db29194878892647c59581b6a70a3..073a7bbf9ae3b7041591d48e4e5b7f3ef0efae42 100644
--- a/src/helpers.py
+++ b/src/helpers.py
@@ -49,9 +49,10 @@ class TimeTracking(object):
     method. Duration can always be shown by printing the time tracking object or calling get_current_duration.
     """
 
-    def __init__(self, start=True):
+    def __init__(self, start=True, name="undefined job"):
         self.start = None
         self.end = None
+        self._name = name
         if start:
             self._start()
 
@@ -93,7 +94,7 @@ class TimeTracking(object):
 
     def __exit__(self, exc_type, exc_val, exc_tb):
         self.stop()
-        logging.info(f"undefined job finished after {self}")
+        logging.info(f"{self._name} finished after {self}")
 
 
 def prepare_host(create_new=True, sampling="daily"):
@@ -147,6 +148,13 @@ def set_experiment_name(experiment_date=None, experiment_path=None, sampling=Non
     return experiment_name, experiment_path
 
 
+def set_bootstrap_path(bootstrap_path, data_path, sampling):
+    if bootstrap_path is None:
+        bootstrap_path = os.path.join(data_path, "..", f"bootstrap_{sampling}")
+    check_path_and_create(bootstrap_path)
+    return bootstrap_path
+
+
 class PyTestRegex:
     """Assert that a given string meets some expectations."""
 
@@ -195,3 +203,18 @@ def float_round(number: float, decimals: int = 0, round_type: Callable = math.ce
 def dict_pop(dict: Dict, pop_keys):
     pop_keys = to_list(pop_keys)
     return {k: v for k, v in dict.items() if k not in pop_keys}
+
+
+def list_pop(list_full: list, pop_items):
+    pop_items = to_list(pop_items)
+    if len(pop_items) > 1:
+        return [e for e in list_full if e not in pop_items]
+    else:
+        list_pop = list_full.copy()
+        list_pop.remove(pop_items[0])
+        return list_pop
+
+
+def dict_pop(dict_orig: Dict, pop_keys):
+    pop_keys = to_list(pop_keys)
+    return {k: v for k, v in dict_orig.items() if k not in pop_keys}
diff --git a/src/model_modules/advanced_paddings.py b/src/model_modules/advanced_paddings.py
new file mode 100644
index 0000000000000000000000000000000000000000..1d48dfc0a87c05183fc5b8b7755f48efaf7b5428
--- /dev/null
+++ b/src/model_modules/advanced_paddings.py
@@ -0,0 +1,280 @@
+__author__ = 'Felix Kleinert'
+__date__ = '2020-03-02'
+
+import tensorflow as tf
+import numpy as np
+import keras.backend as K
+
+from keras.layers.convolutional import _ZeroPadding
+from keras.legacy import interfaces
+from keras.utils import conv_utils
+from keras.utils.generic_utils import transpose_shape
+from keras.backend.common import normalize_data_format
+
+
+class PadUtils:
+    """
+    Helper class for advanced paddings
+    """
+
+    @staticmethod
+    def get_padding_for_same(kernel_size, strides=1):
+        """
+        This methods calculates the padding size to keep input and output dimensions equal for a given kernel size
+        (STRIDES HAVE TO BE EQUAL TO ONE!)
+        :param kernel_size:
+        :return:
+        """
+        if strides != 1:
+            raise NotImplementedError("Strides other than 1 not implemented!")
+        if not all(isinstance(k, int) for k in kernel_size):
+            raise ValueError(f"The `kernel_size` argument must have a tuple of integers. Got: {kernel_size} "
+                             f"of type {[type(k) for k in kernel_size]}")
+
+        ks = np.array(kernel_size, dtype=np.int64)
+
+        if any(k <= 0 for k in ks):
+            raise ValueError(f"All values of kernel_size must be > 0. Got: {kernel_size} ")
+
+        if all(k % 2 == 1 for k in ks):  # (d & 0x1 for d in ks):
+            pad = ((ks - 1) / 2).astype(np.int64)
+            # convert numpy int to base int
+            pad = [np.asscalar(v) for v in pad]
+            return tuple(pad)
+            # return tuple(PadUtils.check_padding_format(pad))
+        else:
+            raise NotImplementedError(f"even kernel size not implemented. Got {kernel_size}")
+
+    @staticmethod
+    def spatial_2d_padding(padding=((1, 1), (1, 1)), data_format=None):
+        """Pads the 2nd and 3rd dimensions of a 4D tensor.
+
+        # Arguments
+            x: Tensor or variable.
+            padding: Tuple of 2 tuples, padding pattern.
+            data_format: string, `"channels_last"` or `"channels_first"`.
+
+        # Returns
+            A padded 4D tensor.
+
+        # Raises
+            ValueError: if `data_format` is neither `"channels_last"` or `"channels_first"`.
+        """
+        assert len(padding) == 2
+        assert len(padding[0]) == 2
+        assert len(padding[1]) == 2
+        data_format = normalize_data_format(data_format)
+
+        pattern = [[0, 0],
+                   list(padding[0]),
+                   list(padding[1]),
+                   [0, 0]]
+        pattern = transpose_shape(pattern, data_format, spatial_axes=(1, 2))
+        return pattern
+
+    @staticmethod
+    def check_padding_format(padding):
+        if isinstance(padding, int):
+            normalized_padding = ((padding, padding), (padding, padding))
+        elif hasattr(padding, '__len__'):
+            if len(padding) != 2:
+                raise ValueError('`padding` should have two elements. '
+                                 'Found: ' + str(padding))
+            for idx_pad, sub_pad in enumerate(padding):
+                if isinstance(sub_pad, str):
+                    raise ValueError(f'`padding[{idx_pad}]` is str but must be int')
+                if hasattr(sub_pad, '__len__'):
+                    if len(sub_pad) != 2:
+                        raise ValueError(f'`padding[{idx_pad}]` should have one or two elements. '
+                                         f'Found: {padding[idx_pad]}')
+                    if not all(isinstance(sub_k, int) for sub_k in padding[idx_pad]):
+                        raise ValueError(f'`padding[{idx_pad}]` should have one or two elements of type int. ' 
+                                         f"Found:{padding[idx_pad]} of type {[type(sub_k) for sub_k in padding[idx_pad]]}")
+            height_padding = conv_utils.normalize_tuple(padding[0], 2,
+                                                        '1st entry of padding')
+            if not all(k >= 0 for k in height_padding):
+                raise ValueError(f"The `1st entry of padding` argument must be >= 0. Received: {padding[0]} of type {type(padding[0])}")
+            width_padding = conv_utils.normalize_tuple(padding[1], 2,
+                                                       '2nd entry of padding')
+            if not all(k >= 0 for k in width_padding):
+                raise ValueError(f"The `2nd entry of padding` argument must be >= 0. Received: {padding[1]} of type {type(padding[1])}")
+            normalized_padding = (height_padding, width_padding)
+        else:
+            raise ValueError('`padding` should be either an int, '
+                             'a tuple of 2 ints '
+                             '(symmetric_height_pad, symmetric_width_pad), '
+                             'or a tuple of 2 tuples of 2 ints '
+                             '((top_pad, bottom_pad), (left_pad, right_pad)). '
+                             f'Found: {padding} of type {type(padding)}')
+        return normalized_padding
+
+
+class ReflectionPadding2D(_ZeroPadding):
+    """
+    Reflection padding layer for 2D input. This custum padding layer is built on keras' zero padding layers. Doc is copy
+    pasted from the original functions/methods:
+
+
+    This layer can add rows and columns of reflected values
+    at the top, bottom, left and right side of an image like tensor.
+
+    Example:
+                                                        6, 5,  4, 5, 6,  5, 4
+                                                              _________
+     1, 2, 3     RefPad(padding=[[1, 1,], [2, 2]])      3, 2,| 1, 2, 3,| 2, 1
+     4, 5, 6     =============================>>>>      6, 5,| 4, 5, 6,| 5, 4
+                                                              _________
+                                                        3, 2,  1, 2, 3,  2, 1
+
+
+
+    '# Arguments
+        padding: int, or tuple of 2 ints, or tuple of 2 tuples of 2 ints.
+            - If int: the same symmetric padding
+                is applied to height and width.
+            - If tuple of 2 ints:
+                interpreted as two different
+                symmetric padding values for height and width:
+                `(symmetric_height_pad, symmetric_width_pad)`.
+            - If tuple of 2 tuples of 2 ints:
+                interpreted as
+                `((top_pad, bottom_pad), (left_pad, right_pad))`
+        data_format: A string,
+            one of `"channels_last"` or `"channels_first"`.
+            The ordering of the dimensions in the inputs.
+            `"channels_last"` corresponds to inputs with shape
+            `(batch, height, width, channels)` while `"channels_first"`
+            corresponds to inputs with shape
+            `(batch, channels, height, width)`.
+            It defaults to the `image_data_format` value found in your
+            Keras config file at `~/.keras/keras.json`.
+            If you never set it, then it will be "channels_last".
+
+    # Input shape
+        4D tensor with shape:
+        - If `data_format` is `"channels_last"`:
+            `(batch, rows, cols, channels)`
+        - If `data_format` is `"channels_first"`:
+            `(batch, channels, rows, cols)`
+
+    # Output shape
+        4D tensor with shape:
+        - If `data_format` is `"channels_last"`:
+            `(batch, padded_rows, padded_cols, channels)`
+        - If `data_format` is `"channels_first"`:
+            `(batch, channels, padded_rows, padded_cols)`
+    '
+    """
+
+    @interfaces.legacy_zeropadding2d_support
+    def __init__(self,
+                 padding=(1, 1),
+                 data_format=None,
+                 **kwargs):
+        normalized_padding = PadUtils.check_padding_format(padding=padding)
+        super(ReflectionPadding2D, self).__init__(normalized_padding,
+                                                  data_format,
+                                                  **kwargs)
+
+    def call(self, inputs, mask=None):
+        pattern = PadUtils.spatial_2d_padding(padding=self.padding, data_format=self.data_format)
+        return tf.pad(inputs, pattern, 'REFLECT')
+
+
+class SymmetricPadding2D(_ZeroPadding):
+    """
+    Symmetric padding layer for 2D input. This custom padding layer is built on keras' zero padding layers. Doc is copy
+    pasted from the original functions/methods:
+
+
+    This layer can add rows and columns of symmetric values
+    at the top, bottom, left and right side of an image like tensor.
+
+        Example:
+                                                        2, 1,  1, 2, 3,  3, 2
+                                                              _________
+     1, 2, 3     SymPad(padding=[[1, 1,], [2, 2]])      2, 1,| 1, 2, 3,| 3, 2
+     4, 5, 6     =============================>>>>      5, 4,| 4, 5, 6,| 6, 5
+                                                              _________
+                                                        5, 4,  4, 5, 6,  6, 5
+
+
+    '# Arguments
+        padding: int, or tuple of 2 ints, or tuple of 2 tuples of 2 ints.
+            - If int: the same symmetric padding
+                is applied to height and width.
+            - If tuple of 2 ints:
+                interpreted as two different
+                symmetric padding values for height and width:
+                `(symmetric_height_pad, symmetric_width_pad)`.
+            - If tuple of 2 tuples of 2 ints:
+                interpreted as
+                `((top_pad, bottom_pad), (left_pad, right_pad))`
+        data_format: A string,
+            one of `"channels_last"` or `"channels_first"`.
+            The ordering of the dimensions in the inputs.
+            `"channels_last"` corresponds to inputs with shape
+            `(batch, height, width, channels)` while `"channels_first"`
+            corresponds to inputs with shape
+            `(batch, channels, height, width)`.
+            It defaults to the `image_data_format` value found in your
+            Keras config file at `~/.keras/keras.json`.
+            If you never set it, then it will be "channels_last".
+
+    # Input shape
+        4D tensor with shape:
+        - If `data_format` is `"channels_last"`:
+            `(batch, rows, cols, channels)`
+        - If `data_format` is `"channels_first"`:
+            `(batch, channels, rows, cols)`
+
+    # Output shape
+        4D tensor with shape:
+        - If `data_format` is `"channels_last"`:
+            `(batch, padded_rows, padded_cols, channels)`
+        - If `data_format` is `"channels_first"`:
+            `(batch, channels, padded_rows, padded_cols)`
+    '
+    """
+
+    @interfaces.legacy_zeropadding2d_support
+    def __init__(self,
+                 padding=(1, 1),
+                 data_format=None,
+                 **kwargs):
+        normalized_padding = PadUtils.check_padding_format(padding=padding)
+        super(SymmetricPadding2D, self).__init__(normalized_padding,
+                                                 data_format,
+                                                 **kwargs)
+
+    def call(self, inputs, mask=None):
+        pattern = PadUtils.spatial_2d_padding(padding=self.padding, data_format=self.data_format)
+        return tf.pad(inputs, pattern, 'SYMMETRIC')
+
+
+if __name__ == '__main__':
+    from keras.models import Model
+    from keras.layers import Conv2D, Flatten, Dense, Input
+
+    kernel_1 = (3, 3)
+    kernel_2 = (5, 5)
+    x = np.array(range(2000)).reshape(-1, 10, 10, 1)
+    y = x.mean(axis=(1, 2))
+
+    x_input = Input(shape=x.shape[1:])
+    pad1 = PadUtils.get_padding_for_same(kernel_size=kernel_1)
+    x_out = ReflectionPadding2D(padding=pad1, name="RefPAD")(x_input)
+    x_out = Conv2D(5, kernel_size=kernel_1, activation='relu')(x_out)
+
+    pad2 = PadUtils.get_padding_for_same(kernel_size=kernel_2)
+    x_out = SymmetricPadding2D(padding=pad2, name="SymPAD")(x_out)
+    x_out = Conv2D(2, kernel_size=kernel_2, activation='relu')(x_out)
+    x_out = Flatten()(x_out)
+    x_out = Dense(1, activation='linear')(x_out)
+
+    model = Model(inputs=x_input, outputs=x_out)
+    model.compile('adam', loss='mse')
+    model.summary()
+    model.fit(x, y, epochs=10)
+
+
diff --git a/src/model_modules/inception_model.py b/src/model_modules/inception_model.py
index 11093df56a4b262f75eae7a6f7c05e6e44e1435d..1cb7656335495f0261abb434e4a203cb4e63887e 100644
--- a/src/model_modules/inception_model.py
+++ b/src/model_modules/inception_model.py
@@ -5,6 +5,7 @@ import logging
 
 import keras
 import keras.layers as layers
+from src.model_modules.advanced_paddings import PadUtils, ReflectionPadding2D, SymmetricPadding2D
 
 
 class InceptionModelBase:
@@ -53,33 +54,39 @@ class InceptionModelBase:
         regularizer = kwargs.get('regularizer', keras.regularizers.l2(0.01))
         bn_settings = kwargs.get('bn_settings', {})
         act_settings = kwargs.get('act_settings', {})
+        padding = kwargs.get('padding', 'ZeroPad2D')
         logging.debug(f'Inception Block with activation: {activation}')
 
         block_name = f'Block_{self.number_of_blocks}{self.block_part_name()}_{tower_kernel[0]}x{tower_kernel[1]}'
+        padding_size = PadUtils.get_padding_for_same(tower_kernel)
 
         if tower_kernel == (1, 1):
             tower = layers.Conv2D(tower_filter,
                                   tower_kernel,
-                                  padding='same',
+                                  padding='valid',
                                   kernel_regularizer=regularizer,
                                   name=block_name)(input_x)
-            tower = self.act(tower, activation, **act_settings)
+            # tower = self.act(tower, activation, **act_settings)
         else:
             tower = layers.Conv2D(reduction_filter,
                                   (1, 1),
-                                  padding='same',
+                                  padding='valid',
                                   kernel_regularizer=regularizer,
                                   name=f'Block_{self.number_of_blocks}{self.block_part_name()}_1x1')(input_x)
             tower = self.act(tower, activation, **act_settings)
 
+            tower = self.padding_layer(padding)(padding=padding_size,
+                                                name=f'Block_{self.number_of_blocks}{self.block_part_name()}_Pad'
+                                                )(tower)
+
             tower = layers.Conv2D(tower_filter,
                                   tower_kernel,
-                                  padding='same',
+                                  padding='valid',
                                   kernel_regularizer=regularizer,
                                   name=block_name)(tower)
-            if batch_normalisation:
-                tower = self.batch_normalisation(tower, **bn_settings)
-            tower = self.act(tower, activation, **act_settings)
+        if batch_normalisation:
+            tower = self.batch_normalisation(tower, **bn_settings)
+        tower = self.act(tower, activation, **act_settings)
 
         return tower
 
@@ -101,6 +108,29 @@ class InceptionModelBase:
         else:
             return act_name.__name__
 
+    @staticmethod
+    def padding_layer(padding):
+        allowed_paddings = {
+            'RefPad2D': ReflectionPadding2D, 'ReflectionPadding2D': ReflectionPadding2D,
+            'SymPad2D': SymmetricPadding2D, 'SymmetricPadding2D': SymmetricPadding2D,
+            'ZeroPad2D': keras.layers.ZeroPadding2D, 'ZeroPadding2D': keras.layers.ZeroPadding2D
+        }
+        if isinstance(padding, str):
+            try:
+                pad2d = allowed_paddings[padding]
+            except KeyError as einfo:
+                raise NotImplementedError(
+                    f"`{einfo}' is not implemented as padding. " 
+                    "Use one of those: i) `RefPad2D', ii) `SymPad2D', iii) `ZeroPad2D'")
+        else:
+            if padding in allowed_paddings.values():
+                pad2d = padding
+            else:
+                raise TypeError(f"`{padding.__name__}' is not a valid padding layer type. "
+                                "Use one of those: "
+                                "i) ReflectionPadding2D, ii) SymmetricPadding2D, iii) ZeroPadding2D")
+        return pad2d
+
     def create_pool_tower(self, input_x, pool_kernel, tower_filter, activation='relu', max_pooling=True, **kwargs):
         """
         This function creates a "MaxPooling tower block"
@@ -114,6 +144,8 @@ class InceptionModelBase:
         self.part_of_block += 1
         self.act_number = 1
         act_settings = kwargs.get('act_settings', {})
+        padding = kwargs.get('padding', 'ZeroPad2D')
+        padding_size = PadUtils.get_padding_for_same(kernel_size=pool_kernel)
 
         # pooling block
         block_name = f"Block_{self.number_of_blocks}{self.block_part_name()}_"
@@ -123,10 +155,12 @@ class InceptionModelBase:
         else:
             block_type = "AvgPool"
             pooling = layers.AveragePooling2D
-        tower = pooling(pool_kernel, strides=(1, 1), padding='same', name=block_name+block_type)(input_x)
+
+        tower = self.padding_layer(padding)(padding=padding_size, name=block_name+'Pad')(input_x)
+        tower = pooling(pool_kernel, strides=(1, 1), padding='valid', name=block_name+block_type)(tower)
 
         # convolution block
-        tower = layers.Conv2D(tower_filter, (1, 1), padding='same', name=block_name+"1x1")(tower)
+        tower = layers.Conv2D(tower_filter, (1, 1), padding='valid', name=block_name+"1x1")(tower)
         tower = self.act(tower, activation, **act_settings)
 
         return tower
@@ -138,21 +172,28 @@ class InceptionModelBase:
         :param tower_conv_parts: dict containing settings for parts of inception block; Example:
                                  tower_conv_parts = {'tower_1': {'reduction_filter': 32,
                                                                  'tower_filter': 64,
-                                                                 'tower_kernel': (3, 1)},
+                                                                 'tower_kernel': (3, 1),
+                                                                 'activation' : 'relu',
+                                                                 'padding' : 'SymPad2D'}
                                                      'tower_2': {'reduction_filter': 32,
                                                                  'tower_filter': 64,
-                                                                 'tower_kernel': (5, 1)},
+                                                                 'tower_kernel': (5, 1),
+                                                                 'activation' : LeakyReLU,
+                                                                 'padding' : keras.layers.ZeroPadding2D}
                                                      'tower_3': {'reduction_filter': 32,
                                                                  'tower_filter': 64,
-                                                                 'tower_kernel': (1, 1)},
+                                                                 'tower_kernel': (1, 1),
+                                                                 'activation' : ELU,
+                                                                 'padding' : src.model_modules.advanced_paddings.ReflectionPadding2D}
                                                     }
         :param tower_pool_parts: dict containing settings for pool part of inception block; Example:
-                                 tower_pool_parts = {'pool_kernel': (3, 1), 'tower_filter': 64}
+                                 tower_pool_parts = {'pool_kernel': (3, 1), 'tower_filter': 64, 'padding': 'RefPad2D'}
         :return:
         """
         self.number_of_blocks += 1
         self.part_of_block = 0
         tower_build = {}
+        block_name = f"Block_{self.number_of_blocks}"
         for part, part_settings in tower_conv_parts.items():
             tower_build[part] = self.create_conv_tower(input_x, **part_settings, **kwargs)
         if 'max_pooling' in tower_pool_parts.keys():
@@ -165,16 +206,46 @@ class InceptionModelBase:
             tower_build['maxpool'] = self.create_pool_tower(input_x, **tower_pool_parts, **kwargs)
             tower_build['avgpool'] = self.create_pool_tower(input_x, **tower_pool_parts, **kwargs, max_pooling=False)
 
-        block = keras.layers.concatenate(list(tower_build.values()), axis=3)
+        block = keras.layers.concatenate(list(tower_build.values()), axis=3,
+                                         name=block_name+"_Co")
         return block
 
 
+# if __name__ == '__main__':
+#     from keras.models import Model
+#     from keras.layers import Conv2D, Flatten, Dense, Input
+#     import numpy as np
+#
+#
+#     kernel_1 = (3, 3)
+#     kernel_2 = (5, 5)
+#     x = np.array(range(2000)).reshape(-1, 10, 10, 1)
+#     y = x.mean(axis=(1, 2))
+#
+#     x_input = Input(shape=x.shape[1:])
+#     pad1 = PadUtils.get_padding_for_same(kernel_size=kernel_1)
+#     x_out = InceptionModelBase.padding_layer('RefPad2D')(padding=pad1, name="RefPAD1")(x_input)
+#     # x_out = ReflectionPadding2D(padding=pad1, name="RefPAD")(x_input)
+#     x_out = Conv2D(5, kernel_size=kernel_1, activation='relu')(x_out)
+#
+#     pad2 = PadUtils.get_padding_for_same(kernel_size=kernel_2)
+#     x_out = InceptionModelBase.padding_layer(SymmetricPadding2D)(padding=pad2, name="SymPAD1")(x_out)
+#     # x_out = SymmetricPadding2D(padding=pad2, name="SymPAD")(x_out)
+#     x_out = Conv2D(2, kernel_size=kernel_2, activation='relu')(x_out)
+#     x_out = Flatten()(x_out)
+#     x_out = Dense(1, activation='linear')(x_out)
+#
+#     model = Model(inputs=x_input, outputs=x_out)
+#     model.compile('adam', loss='mse')
+#     model.summary()
+#     # model.fit(x, y, epochs=10)
+
 if __name__ == '__main__':
     print(__name__)
     from keras.datasets import cifar10
     from keras.utils import np_utils
     from keras.layers import Input
-    from keras.layers.advanced_activations import LeakyReLU
+    from keras.layers.advanced_activations import LeakyReLU, ELU
     from keras.optimizers import SGD
     from keras.layers import Dense, Flatten, Conv2D, MaxPooling2D
     from keras.models import Model
@@ -183,11 +254,17 @@ if __name__ == '__main__':
     conv_settings_dict = {'tower_1': {'reduction_filter': 64,
                                       'tower_filter': 64,
                                       'tower_kernel': (3, 3),
-                                      'activation': LeakyReLU},
+                                      'activation': LeakyReLU,},
                           'tower_2': {'reduction_filter': 64,
                                       'tower_filter': 64,
                                       'tower_kernel': (5, 5),
-                                      'activation': 'relu'}
+                                      'activation': 'relu',
+                                      'padding': 'SymPad2D'},
+                          'tower_3': {'reduction_filter': 64,
+                                      'tower_filter': 64,
+                                      'tower_kernel': (1, 1),
+                                      'activation': ELU,
+                                      'padding': ReflectionPadding2D}
                           }
     pool_settings_dict = {'pool_kernel': (3, 3),
                           'tower_filter': 64,
@@ -205,17 +282,21 @@ if __name__ == '__main__':
 
     # create inception net
     inception_net = InceptionModelBase()
-    output = inception_net.inception_block(input_img, conv_settings_dict, pool_settings_dict)
+    output = inception_net.inception_block(input_img, conv_settings_dict, pool_settings_dict, batch_normalisation=True)
     output = Flatten()(output)
     output = Dense(10, activation='softmax')(output)
     model = Model(inputs=input_img, outputs=output)
     print(model.summary())
 
     # compile
-    epochs = 10
+    epochs = 1
     lrate = 0.01
     decay = lrate/epochs
     sgd = SGD(lr=lrate, momentum=0.9, decay=decay, nesterov=False)
     model.compile(loss='categorical_crossentropy', optimizer=sgd, metrics=['accuracy'])
     print(X_train.shape)
     keras.utils.plot_model(model, to_file='model.pdf', show_shapes=True, show_layer_names=True)
+    # model.fit(X_train, y_train, epochs=epochs, validation_data=(X_test, y_test))
+    print('test')
+
+
diff --git a/src/plotting/postprocessing_plotting.py b/src/plotting/postprocessing_plotting.py
index a41c636b5ab17d2039f7976fca625e9c8e11ce6e..854182613cdb63456dc8f62d2421560d829ee629 100644
--- a/src/plotting/postprocessing_plotting.py
+++ b/src/plotting/postprocessing_plotting.py
@@ -479,6 +479,76 @@ class PlotCompetitiveSkillScore(RunEnvironment):
         plt.close()
 
 
+class PlotBootstrapSkillScore(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 = ".", 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 model_setup: architecture type to specify plot name (default "CNN")
+        """
+        super().__init__()
+        self._labels = None
+        self._data = self._prepare_data(data)
+        self._plot(plot_folder, model_setup)
+
+    def _prepare_data(self, data: Dict) -> 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
+        :return: pre-processed data set
+        """
+        data = helpers.dict_to_xarray(data, "station")
+        self._labels = [str(i) + "d" for i in data.coords["ahead"].values]
+        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,  model_setup):
+        """
+        Main plot function to plot climatological skill score.
+        :param plot_folder: path to save the plot
+        :param model_setup: architecture type to specify plot name
+        """
+        fig, ax = plt.subplots()
+        sns.boxplot(x="boot_var", 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"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, model_setup)
+
+    @staticmethod
+    def _save(plot_folder, 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 model_setup: architecture type to specify plot name
+        """
+        plot_name = os.path.join(plot_folder, f"skill_score_bootstrap_{model_setup}.pdf")
+        logging.debug(f"... save plot to {plot_name}")
+        plt.savefig(plot_name, dpi=500)
+        plt.close('all')
+
+
 class PlotTimeSeries(RunEnvironment):
 
     def __init__(self, stations: List, data_path: str, name: str, window_lead_time: int = None, plot_folder: str = ".",
@@ -569,15 +639,16 @@ class PlotTimeSeries(RunEnvironment):
     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)
+            plot_data = data.sel(type="CNN", ahead=ahead).drop(["type", "ahead"]).squeeze().shift(index=ahead)
             label = f"{ahead}{self._sampling}"
-            ax.plot(index, plot_data.values, color=color[ahead-1], label=label)
+            ax.plot(plot_data, color=color[ahead-1], label=label)
 
     def _plot_obs(self, ax, data):
-        obs_data = data.sel(type="obs", ahead=1)
-        index = data.index + np.timedelta64(1, self._sampling)
-        ax.plot(index, obs_data.values, color=matplotlib.colors.cnames["green"], label="obs")
+        ahead = 1
+        obs_data = data.sel(type="obs", ahead=ahead).shift(index=ahead)
+        # index = data.index + np.timedelta64(1, self._sampling)
+        # ax.plot(index, obs_data.values, color=matplotlib.colors.cnames["green"], label="obs")
+        ax.plot(obs_data, color=matplotlib.colors.cnames["green"], label="obs")
 
     @staticmethod
     def _get_time_range(data):
diff --git a/src/run_modules/experiment_setup.py b/src/run_modules/experiment_setup.py
index 4c3b8872575ea9929f1d4ba3f5a42e222ac2fff4..56c22a81e48421438816855770b7477e84e3a8d8 100644
--- a/src/run_modules/experiment_setup.py
+++ b/src/run_modules/experiment_setup.py
@@ -28,13 +28,13 @@ class ExperimentSetup(RunEnvironment):
     trainable: Train new model if true, otherwise try to load existing model
     """
 
-    def __init__(self, parser_args=None, var_all_dict=None, stations=None, network=None, station_type=None, variables=None,
+    def __init__(self, parser_args=None, stations=None, network=None, station_type=None, variables=None,
                  statistics_per_var=None, start=None, end=None, window_history_size=None, target_var="o3", target_dim=None,
                  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=None, fraction_of_train=None,
                  experiment_path=None, plot_path=None, forecast_path=None, overwrite_local_data=None, sampling="daily",
-                 create_new_model=None, permute_data_on_training=None, transformation=None):
+                 create_new_model=None, bootstrap_path=None, permute_data_on_training=None, transformation=None):
 
         # create run framework
         super().__init__()
@@ -44,6 +44,9 @@ class ExperimentSetup(RunEnvironment):
         self._set_param("create_new_model", create_new_model, default=True)
         if self.data_store.get("create_new_model", "general"):
             trainable = True
+        data_path = self.data_store.get("data_path", "general")
+        bootstrap_path = helpers.set_bootstrap_path(bootstrap_path, data_path, sampling)
+        self._set_param("bootstrap_path", bootstrap_path)
         self._set_param("trainable", trainable, default=True)
         self._set_param("fraction_of_training", fraction_of_train, default=0.8)
         self._set_param("permute_data", permute_data_on_training, default=False, scope="general.train")
@@ -67,12 +70,11 @@ class ExperimentSetup(RunEnvironment):
         helpers.check_path_and_create(self.data_store.get("forecast_path", "general"))
 
         # setup for data
-        self._set_param("var_all_dict", var_all_dict, default=DEFAULT_VAR_ALL_DICT)
         self._set_param("stations", stations, default=DEFAULT_STATIONS)
         self._set_param("network", network, default="AIRBASE")
         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._set_param("statistics_per_var", statistics_per_var, default=DEFAULT_VAR_ALL_DICT)
+        self._set_param("variables", variables, default=list(self.data_store.get("statistics_per_var", "general").keys()))
         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")
@@ -84,6 +86,7 @@ class ExperimentSetup(RunEnvironment):
 
         # target
         self._set_param("target_var", target_var, default="o3")
+        self._check_target_var()
         self._set_param("target_dim", target_dim, default='variables')
         self._set_param("window_lead_time", window_lead_time, default=3)
 
@@ -133,16 +136,27 @@ class ExperimentSetup(RunEnvironment):
             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")
+        var = self.data_store.get("variables", "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}")
 
+    def _check_target_var(self):
+        target_var = helpers.to_list(self.data_store.get("target_var", "general"))
+        stat = self.data_store.get("statistics_per_var", "general")
+        var = self.data_store.get("variables", "general")
+        if not set(target_var).issubset(stat.keys()):
+            raise ValueError(f"Could not find target variable {target_var} in statistics_per_var.")
+        unused_vars = set(stat.keys()).difference(set(var).union(target_var))
+        if len(unused_vars) > 0:
+            logging.info(f"There are unused keys in statistics_per_var. Therefore remove keys: {unused_vars}")
+            stat_new = helpers.dict_pop(stat, list(unused_vars))
+            self._set_param("statistics_per_var", stat_new)
+
 
 if __name__ == "__main__":
 
diff --git a/src/run_modules/model_setup.py b/src/run_modules/model_setup.py
index e3945a542d60b09dc9855bd28be87cdba729ed72..d2c8d93fc957ecb2990e99000cbd3588e2d83eef 100644
--- a/src/run_modules/model_setup.py
+++ b/src/run_modules/model_setup.py
@@ -10,8 +10,8 @@ import tensorflow as tf
 
 from src.model_modules.keras_extensions import HistoryAdvanced, CallbackHandler
 # from src.model_modules.model_class import MyBranchedModel as MyModel
-# from src.model_modules.model_class import MyLittleModel as MyModel
-from src.model_modules.model_class import MyTowerModel as MyModel
+from src.model_modules.model_class import MyLittleModel as MyModel
+# from src.model_modules.model_class import MyTowerModel as MyModel
 from src.run_modules.run_environment import RunEnvironment
 
 
diff --git a/src/run_modules/post_processing.py b/src/run_modules/post_processing.py
index 06203c879872891f57c719040482fe052824c65e..0a61ee4f07d0c6eccf698aa16d3de9d7275e75f6 100644
--- a/src/run_modules/post_processing.py
+++ b/src/run_modules/post_processing.py
@@ -13,11 +13,12 @@ import xarray as xr
 from src import statistics
 from src.data_handling.data_distributor import Distributor
 from src.data_handling.data_generator import DataGenerator
+from src.data_handling.bootstraps import BootStraps
 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
+    PlotCompetitiveSkillScore, PlotTimeSeries, PlotBootstrapSkillScore
 from src.plotting.postprocessing_plotting import plot_conditional_quantiles
 from src.run_modules.run_environment import RunEnvironment
 
@@ -37,6 +38,7 @@ class PostProcessing(RunEnvironment):
         self.target_var = self.data_store.get("target_var", "general")
         self._sampling = self.data_store.get("sampling", "general")
         self.skill_scores = None
+        self.bootstrap_skill_scores = None
         self._run()
 
     def _run(self):
@@ -48,9 +50,67 @@ class PostProcessing(RunEnvironment):
             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.bootstrap_skill_scores = self.create_boot_straps()
         self.skill_scores = self.calculate_skill_scores()
         self.plot()
 
+    def create_boot_straps(self):
+
+        # forecast
+
+        bootstrap_path = self.data_store.get("bootstrap_path", "general")
+        forecast_path = self.data_store.get("forecast_path", "general")
+        window_lead_time = self.data_store.get("window_lead_time", "general")
+        bootstraps = BootStraps(self.test_data, bootstrap_path, 20)
+        with TimeTracking(name="boot predictions"):
+            bootstrap_predictions = self.model.predict_generator(generator=bootstraps.boot_strap_generator(),
+                                                                 steps=bootstraps.get_boot_strap_generator_length())
+        if isinstance(bootstrap_predictions, list):
+            bootstrap_predictions = bootstrap_predictions[-1]
+        bootstrap_meta = np.array(bootstraps.get_boot_strap_meta())
+        variables = np.unique(bootstrap_meta[:, 0])
+        for station in np.unique(bootstrap_meta[:, 1]):
+            coords = None
+            for boot in variables:
+                ind = np.all(bootstrap_meta == [boot, station], axis=1)
+                length = sum(ind)
+                sel = bootstrap_predictions[ind].reshape((length, window_lead_time, 1))
+                coords = (range(length), range(1, window_lead_time + 1))
+                tmp = xr.DataArray(sel, coords=(*coords, [boot]), dims=["index", "ahead", "type"])
+                file_name = os.path.join(forecast_path, f"bootstraps_{boot}_{station}.nc")
+                tmp.to_netcdf(file_name)
+            labels = bootstraps.get_labels(station).reshape((length, window_lead_time, 1))
+            file_name = os.path.join(forecast_path, f"bootstraps_labels_{station}.nc")
+            labels = xr.DataArray(labels, coords=(*coords, ["obs"]), dims=["index", "ahead", "type"])
+            labels.to_netcdf(file_name)
+
+        # file_name = os.path.join(forecast_path, f"bootstraps_orig.nc")
+        # orig = xr.open_dataarray(file_name)
+
+
+        # calc skill scores
+        skill_scores = statistics.SkillScores(None)
+        score = {}
+        for station in np.unique(bootstrap_meta[:, 1]):
+            file_name = os.path.join(forecast_path, f"bootstraps_labels_{station}.nc")
+            labels = xr.open_dataarray(file_name)
+            shape = labels.shape
+            orig = bootstraps.get_orig_prediction(forecast_path,  f"forecasts_norm_{station}_test.nc").reshape(shape)
+            orig = xr.DataArray(orig, coords=(range(shape[0]), range(1, shape[1] + 1), ["orig"]), dims=["index", "ahead", "type"])
+            skill = pd.DataFrame(columns=range(1, window_lead_time + 1))
+            for boot in variables:
+                file_name = os.path.join(forecast_path, f"bootstraps_{boot}_{station}.nc")
+                boot_data = xr.open_dataarray(file_name)
+                boot_data = boot_data.combine_first(labels)
+                boot_data = boot_data.combine_first(orig)
+                boot_scores = []
+                for iahead in range(window_lead_time):
+                    data = boot_data.sel(ahead=iahead + 1)
+                    boot_scores.append(skill_scores.general_skill_score(data, forecast_name=boot, reference_name="orig"))
+                skill.loc[boot] = np.array(boot_scores)
+            score[station] = xr.DataArray(skill, dims=["boot_var", "ahead"])
+        return score
+
     def _load_model(self):
         try:
             model = self.data_store.get("best_model", "general")
@@ -75,6 +135,7 @@ class PostProcessing(RunEnvironment):
         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")
+        PlotBootstrapSkillScore(self.bootstrap_skill_scores, 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):
@@ -96,64 +157,71 @@ class PostProcessing(RunEnvironment):
         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 = data.get_transposed_history()
 
             # get scaling parameters
             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)
+            for normalised in [True, False]:
+                # create empty arrays
+                nn_prediction, persistence_prediction, ols_prediction, observation = self._create_empty_prediction_arrays(data, count=4)
+
+                # nn forecast
+                nn_prediction = self._create_nn_forecast(input_data, nn_prediction, mean, std, transformation_method, normalised)
 
-            # persistence
-            persistence_prediction = self._create_persistence_forecast(input_data, persistence_prediction, mean, std, 
-                                                                       transformation_method)
+                # persistence
+                persistence_prediction = self._create_persistence_forecast(data, persistence_prediction, mean, std,
+                                                                           transformation_method, normalised)
 
-            # ols
-            ols_prediction = self._create_ols_forecast(input_data, ols_prediction, mean, std, transformation_method)
+                # ols
+                ols_prediction = self._create_ols_forecast(input_data, ols_prediction, mean, std, transformation_method, normalised)
 
-            # observation
-            observation = self._create_observation(data, None, mean, std, transformation_method)
+                # observation
+                observation = self._create_observation(data, observation, mean, std, transformation_method, normalised)
 
-            # merge all predictions
-            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,
-                                                          obs=observation,
-                                                          OLS=ols_prediction)
+                # merge all predictions
+                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,
+                                                              obs=observation,
+                                                              OLS=ols_prediction)
 
-            # save all forecasts locally
-            path = self.data_store.get("forecast_path", "general")
-            file = os.path.join(path, f"forecasts_{data.station[0]}_test.nc")
-            all_predictions.to_netcdf(file)
+                # save all forecasts locally
+                path = self.data_store.get("forecast_path", "general")
+                prefix = "forecasts_norm" if normalised else "forecasts"
+                file = os.path.join(path, f"{prefix}_{data.station[0]}_test.nc")
+                all_predictions.to_netcdf(file)
 
     def _get_frequency(self):
         getter = {"daily": "1D", "hourly": "1H"}
         return getter.get(self._sampling, None)
 
-    @staticmethod
-    def _create_observation(data, _, mean, std, transformation_method):
-        return statistics.apply_inverse_transformation(data.label.copy(), mean, std, transformation_method)
+    def _create_observation(self, data, _, mean, std, transformation_method, normalised):
+        obs = data.label.copy()
+        if not normalised:
+            obs = statistics.apply_inverse_transformation(obs, mean, std, transformation_method)
+        return obs
 
-    def _create_ols_forecast(self, input_data, ols_prediction, mean, std, transformation_method):
+    def _create_ols_forecast(self, input_data, ols_prediction, mean, std, transformation_method, normalised):
         tmp_ols = self.ols_model.predict(input_data)
-        tmp_ols = statistics.apply_inverse_transformation(tmp_ols, mean, std, transformation_method)
+        if not normalised:
+            tmp_ols = statistics.apply_inverse_transformation(tmp_ols, mean, std, transformation_method)
         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': self.target_var})
-        tmp_persi = statistics.apply_inverse_transformation(tmp_persi, mean, std, transformation_method)
+    def _create_persistence_forecast(self, data, persistence_prediction, mean, std, transformation_method, normalised):
+        tmp_persi = data.observation.copy().sel({'window': 0})
+        if not normalised:
+            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)),
                                                        axis=1)
         return persistence_prediction
 
-    def _create_nn_forecast(self, input_data, nn_prediction, mean, std, transformation_method):
+    def _create_nn_forecast(self, input_data, nn_prediction, mean, std, transformation_method, normalised):
         """
         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
@@ -166,8 +234,11 @@ class PostProcessing(RunEnvironment):
         :return:
         """
         tmp_nn = self.model.predict(input_data)
-        tmp_nn = statistics.apply_inverse_transformation(tmp_nn, mean, std, transformation_method)
-        if tmp_nn.ndim == 3:
+        if not normalised:
+            tmp_nn = statistics.apply_inverse_transformation(tmp_nn, mean, std, transformation_method)
+        if isinstance(tmp_nn, list):
+            nn_prediction.values = np.swapaxes(np.expand_dims(tmp_nn[-1], axis=1), 2, 0)
+        elif 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)
@@ -227,7 +298,7 @@ class PostProcessing(RunEnvironment):
         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_observation(data, None, mean, std, transformation_method)
+            external_data = self._create_observation(data, None, mean, std, transformation_method, normalised=False)
             external_data = external_data.squeeze("Stations").sel(window=1).drop(["window", "Stations", "variables"])
             return external_data.rename({'datetime': 'index'})
         except KeyError:
diff --git a/test/test_data_handling/test_bootstraps.py b/test/test_data_handling/test_bootstraps.py
new file mode 100644
index 0000000000000000000000000000000000000000..9dd23893ef903bfbd0595a482dceb32724c3b437
--- /dev/null
+++ b/test/test_data_handling/test_bootstraps.py
@@ -0,0 +1,64 @@
+
+from src.data_handling.bootstraps import BootStraps
+
+import pytest
+import os
+
+import numpy as np
+
+
+class TestBootstraps:
+
+    @pytest.fixture
+    def path(self):
+        path = os.path.join(os.path.dirname(__file__), "data")
+        if not os.path.exists(path):
+            os.makedirs(path)
+        return path
+
+    @pytest.fixture
+    def boot_no_init(self, path):
+        obj = object.__new__(BootStraps)
+        super(BootStraps, obj).__init__()
+        obj.number_bootstraps = 50
+        obj.bootstrap_path = path
+        return obj
+
+    def test_valid_bootstrap_file(self, path, boot_no_init):
+        station = "TESTSTATION"
+        variables = "var1_var2_var3"
+        window = 5
+        # empty case
+        assert len(os.listdir(path)) == 0
+        assert boot_no_init.valid_bootstrap_file(station, variables, window) == (False, 50)
+        # different cases, where files with bigger range are existing
+        os.mknod(os.path.join(path, f"{station}_{variables}_hist5_nboots50_shuffled.dat"))
+        assert boot_no_init.valid_bootstrap_file(station, variables, window) == (True, None)
+        os.mknod(os.path.join(path, f"{station}_{variables}_hist5_nboots100_shuffled.dat"))
+        assert boot_no_init.valid_bootstrap_file(station, variables, window) == (True, None)
+        os.mknod(os.path.join(path, f"{station}_{variables}_hist10_nboots50_shuffled.dat"))
+        os.mknod(os.path.join(path, f"{station}1_{variables}_hist10_nboots50_shuffled.dat"))
+        assert boot_no_init.valid_bootstrap_file(station, variables, window) == (True, None)
+        #  need to reload data and therefore remove not fitting files for this station
+        assert boot_no_init.valid_bootstrap_file(station, variables, 20) == (False, 100)
+        assert len(os.listdir(path)) == 1
+        # reload because expanded boot number
+        os.mknod(os.path.join(path, f"{station}_{variables}_hist5_nboots50_shuffled.dat"))
+        boot_no_init.number_bootstraps = 60
+        assert boot_no_init.valid_bootstrap_file(station, variables, window) == (False, 60)
+        assert len(os.listdir(path)) == 1
+        # reload because of expanded window size, but use maximum boot number from file names
+        os.mknod(os.path.join(path, f"{station}_{variables}_hist5_nboots60_shuffled.dat"))
+        boot_no_init.number_bootstraps = 50
+        assert boot_no_init.valid_bootstrap_file(station, variables, 20) == (False, 60)
+
+    def test_shuffle_single_variale(self, boot_no_init):
+        data = np.array([[1, 2, 3], [1, 2, 3], [1, 2, 3], [1, 2, 3]])
+        res = boot_no_init.shuffle_single_variable(data, chunks=(2, 3)).compute()
+        assert res.shape == data.shape
+        assert res.max() == data.max()
+        assert res.min() == data.min()
+        assert set(np.unique(res)).issubset({1, 2, 3})
+
+    def test_create_shuffled_data(self):
+        pass
\ No newline at end of file
diff --git a/test/test_data_handling/test_data_preparation.py b/test/test_data_handling/test_data_preparation.py
index ac449c4dc6d4c83a457eccc93a766ec4f17f58c9..91719f3dd16326ee6281c4db8ef3aa87e238d70f 100644
--- a/test/test_data_handling/test_data_preparation.py
+++ b/test/test_data_handling/test_data_preparation.py
@@ -27,6 +27,7 @@ class TestDataPrep:
         d.network = 'UBA'
         d.station = ['DEBW107']
         d.variables = ['o3', 'temp']
+        d.statistics_per_var = {'o3': 'dma8eu', 'temp': 'maximum'}
         d.station_type = "background"
         d.sampling = "daily"
         d.kwargs = None
@@ -39,7 +40,7 @@ class TestDataPrep:
         assert data.variables == ['o3', 'temp']
         assert data.station_type == "background"
         assert data.statistics_per_var == {'o3': 'dma8eu', 'temp': 'maximum'}
-        assert not all([data.mean, data.std, data.history, data.label, data.station_type])
+        assert not any([data.mean, data.std, data.history, data.label, data.observation])
         assert {'test': 'testKWARGS'}.items() <= data.kwargs.items()
 
     def test_init_no_stats(self):
@@ -125,6 +126,7 @@ class TestDataPrep:
         d.path = os.path.join(os.path.abspath(os.path.dirname(__file__)), "data")
         d.station = 'TESTSTATION'
         d.variables = ['a', 'bc']
+        d.statistics_per_var = {'a': 'dma8eu', 'bc': 'maximum'}
         assert d._set_file_name() == os.path.join(os.path.abspath(os.path.dirname(__file__)),
                                                   "data/TESTSTATION_a_bc.nc")
         assert d._set_meta_file_name() == os.path.join(os.path.abspath(os.path.dirname(__file__)),
@@ -258,29 +260,32 @@ class TestDataPrep:
         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
-        data.history_label_nan_remove('datetime')
-        assert data.history is None
-        assert data.label is None
-        data.make_history_window('datetime', 6)
+    def test_remove_nan_no_hist_or_label(self, data):
+        assert not any([data.history, data.label, data.observation])
+        data.remove_nan('datetime')
+        assert not any([data.history, data.label, data.observation])
+        data.make_history_window('variables', 6, 'datetime')
         assert data.history is not None
-        data.history_label_nan_remove('datetime')
+        data.remove_nan('datetime')
         assert data.history is None
         data.make_labels('variables', 'o3', 'datetime', 2)
-        assert data.label is not None
-        data.history_label_nan_remove('datetime')
-        assert data.label is None
+        data.make_observation('variables', 'o3', 'datetime')
+        assert all(map(lambda x: x is not None, [data.label, data.observation]))
+        data.remove_nan('datetime')
+        assert not any([data.history, data.label, data.observation])
 
-    def test_nan_remove(self, data):
-        data.make_history_window('datetime', -12)
+    def test_remove_nan(self, data):
+        data.make_history_window('variables', -12, 'datetime')
         data.make_labels('variables', 'o3', 'datetime', 3)
+        data.make_observation('variables', 'o3', 'datetime')
         shape = data.history.shape
-        data.history_label_nan_remove('datetime')
+        data.remove_nan('datetime')
         assert data.history.isnull().sum() == 0
         assert itemgetter(0, 1, 3)(shape) == itemgetter(0, 1, 3)(data.history.shape)
         assert shape[2] >= data.history.shape[2]
+        remaining_len = data.history.datetime.shape
+        assert remaining_len == data.label.datetime.shape
+        assert remaining_len == data.observation.datetime.shape
 
     def test_create_index_array(self, data):
         index_array = data.create_index_array('window', range(1, 4))
@@ -310,34 +315,52 @@ class TestDataPrep:
         res = data.shift('datetime', 4)
         window, orig = self.extract_window_data(res, data.data, 4)
         assert res.coords.dims == ('window', 'Stations', 'datetime', 'variables')
-        assert list(res.data.shape) == [4] + list(data.data.shape)
+        assert list(res.data.shape) == [4, *data.data.shape]
         assert np.testing.assert_array_equal(orig, window) is None
         res = data.shift('datetime', -3)
         window, orig = self.extract_window_data(res, data.data, -3)
-        assert list(res.data.shape) == [4] + list(data.data.shape)
+        assert list(res.data.shape) == [4, *data.data.shape]
         assert np.testing.assert_array_equal(orig, window) is None
         res = data.shift('datetime', 0)
         window, orig = self.extract_window_data(res, data.data, 0)
-        assert list(res.data.shape) == [1] + list(data.data.shape)
+        assert list(res.data.shape) == [1, *data.data.shape]
         assert np.testing.assert_array_equal(orig, window) is None
 
     def test_make_history_window(self, data):
         assert data.history is None
-        data.make_history_window('datetime', 5)
+        data.make_history_window("variables", 5, "datetime")
         assert data.history is not None
         save_history = data.history
-        data.make_history_window('datetime', -5)
+        data.make_history_window("variables", -5, "datetime")
         assert np.testing.assert_array_equal(data.history, save_history) is None
 
     def test_make_labels(self, data):
         assert data.label is None
         data.make_labels('variables', 'o3', 'datetime', 3)
         assert data.label.variables.data == 'o3'
-        assert list(data.label.shape) == [3] + list(data.data.shape)[:2]
-        save_label = data.label
+        assert list(data.label.shape) == [3, *data.data.shape[:2]]
+        save_label = data.label.copy()
         data.make_labels('variables', 'o3', 'datetime', -3)
         assert np.testing.assert_array_equal(data.label, save_label) is None
 
+    def test_make_labels_multiple(self, data):
+        assert data.label is None
+        data.make_labels("variables", ["o3", "temp"], "datetime", 4)
+        assert all(data.label.variables.data == ["o3", "temp"])
+        assert list(data.label.shape) == [4, *data.data.shape[:2], 2]
+
+    def test_make_observation(self, data):
+        assert data.observation is None
+        data.make_observation("variables", "o3", "datetime")
+        assert data.observation.variables.data == "o3"
+        assert list(data.observation.shape) == [1, 1, data.data.datetime.shape[0]]
+
+    def test_make_observation_multiple(self, data):
+        assert data.observation is None
+        data.make_observation("variables", ["o3", "temp"], "datetime")
+        assert all(data.observation.variables.data == ["o3", "temp"])
+        assert list(data.observation.shape) == [1, 1, data.data.datetime.shape[0], 2]
+
     def test_slice(self, data):
         res = data._slice(data.data, dt.date(1997, 1, 1), dt.date(1997, 1, 10), 'datetime')
         assert itemgetter(0, 2)(res.shape) == itemgetter(0, 2)(data.data.shape)
@@ -363,3 +386,12 @@ class TestDataPrep:
             data_new = DataPrep(os.path.join(os.path.dirname(__file__), 'data'), 'dummy', 'DEBW107', ['o3', 'temp'],
                                 station_type='traffic', statistics_per_var={'o3': 'dma8eu', 'temp': 'maximum'})
 
+    def test_get_transposed_history(self, data):
+        data.make_history_window("variables", 3, "datetime")
+        transposed = data.get_transposed_history()
+        assert transposed.coords.dims == ("datetime", "window", "Stations", "variables")
+
+    def test_get_transposed_label(self, data):
+        data.make_labels("variables", "o3", "datetime", 2)
+        transposed = data.get_transposed_label()
+        assert transposed.coords.dims == ("datetime", "window")
diff --git a/test/test_helpers.py b/test/test_helpers.py
index b807d2b8612b9ee006bff43f1ae4cfcfd2dd07e1..07ec244e078f977dca761274260275aab355c183 100644
--- a/test/test_helpers.py
+++ b/test/test_helpers.py
@@ -117,6 +117,14 @@ class TestTimeTracking:
         expression = PyTestRegex(r"undefined job finished after \d+:\d+:\d+ \(hh:mm:ss\)")
         assert caplog.record_tuples[-1] == ('root', 20, expression)
 
+    def test_name_enter_exit(self, caplog):
+        caplog.set_level(logging.INFO)
+        with TimeTracking(name="my job") as t:
+            assert t.start is not None
+            assert t.end is None
+        expression = PyTestRegex(r"my job finished after \d+:\d+:\d+ \(hh:mm:ss\)")
+        assert caplog.record_tuples[-1] == ('root', 20, expression)
+
 
 class TestPrepareHost:
 
diff --git a/test/test_model_modules/test_advanced_paddings.py b/test/test_model_modules/test_advanced_paddings.py
new file mode 100644
index 0000000000000000000000000000000000000000..5282eb6df34d4d395dbbdd1fd76fd71a95e9c8df
--- /dev/null
+++ b/test/test_model_modules/test_advanced_paddings.py
@@ -0,0 +1,419 @@
+import keras
+import pytest
+
+from src.model_modules.advanced_paddings import *
+
+
+class TestPadUtils:
+
+    def test_get_padding_for_same_negative_kernel_size(self):
+        print('In test_get_padding_for_same_negative_kernel_size')
+        with pytest.raises(ValueError) as einfo:
+            PadUtils.get_padding_for_same((-1, 2))
+        assert 'All values of kernel_size must be > 0. Got: (-1, 2) ' in str(einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            PadUtils.get_padding_for_same((1, -2))
+        assert 'All values of kernel_size must be > 0. Got: (1, -2) ' in str(einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            PadUtils.get_padding_for_same((-1, -2))
+        assert 'All values of kernel_size must be > 0. Got: (-1, -2) ' in str(einfo.value)
+
+    def test_get_padding_for_same_strides_greater_one(self):
+        with pytest.raises(NotImplementedError) as einfo:
+            PadUtils.get_padding_for_same((1, 1), strides=2)
+        assert 'Strides other than 1 not implemented!' in str(einfo.value)
+
+        with pytest.raises(NotImplementedError) as einfo:
+            PadUtils.get_padding_for_same((1, 1), strides=-1)
+        assert 'Strides other than 1 not implemented!' in str(einfo.value)
+
+    def test_get_padding_for_same_non_int_kernel(self):
+        with pytest.raises(ValueError) as einfo:
+            PadUtils.get_padding_for_same((1., 1))
+        assert "The `kernel_size` argument must have a tuple of integers. Got: (1.0, 1) " \
+               "of type [<class 'float'>, <class 'int'>]" in str(einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            PadUtils.get_padding_for_same((1, 1.))
+        assert "The `kernel_size` argument must have a tuple of integers. Got: (1, 1.0) " \
+               "of type [<class 'int'>, <class 'float'>]" in str(einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            PadUtils.get_padding_for_same((1, '1.'))
+        assert "The `kernel_size` argument must have a tuple of integers. Got: (1, '1.') " \
+               "of type [<class 'int'>, <class 'str'>]" in str(einfo.value)
+
+    def test_get_padding_for_same_stride_3d(self):
+        kernel = (3, 3, 3)
+        pad = PadUtils.get_padding_for_same(kernel)
+        assert pad == (1, 1, 1)
+        assert isinstance(pad, tuple)
+        assert isinstance(pad[0], int) and isinstance(pad[1], int)
+        assert not (isinstance(pad[0], np.int64) and isinstance(pad[1], np.int64) and isinstance(pad[2], np.int64))
+
+    def test_get_padding_for_same_even_pad(self):
+        with pytest.raises(NotImplementedError) as einfo:
+            PadUtils.get_padding_for_same((2, 1))
+        assert 'even kernel size not implemented. Got (2, 1)' in str(einfo.value)
+
+        with pytest.raises(NotImplementedError) as einfo:
+            PadUtils.get_padding_for_same((1, 4))
+        assert 'even kernel size not implemented. Got (1, 4)' in str(einfo.value)
+
+        with pytest.raises(NotImplementedError) as einfo:
+            PadUtils.get_padding_for_same((2, 4))
+        assert 'even kernel size not implemented. Got (2, 4)' in str(einfo.value)
+
+    ##################################################################################
+
+    def test_check_padding_format_negative_pads(self):
+
+        with pytest.raises(ValueError) as einfo:
+            PadUtils.check_padding_format((-2, 1))
+        assert "The `1st entry of padding` argument must be >= 0. Received: -2 of type <class 'int'>" in str(einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            PadUtils.check_padding_format((1, -1))
+        assert "The `2nd entry of padding` argument must be >= 0. Received: -1 of type <class 'int'>" in str(
+            einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            PadUtils.check_padding_format((-2, -1))
+        assert "The `1st entry of padding` argument must be >= 0. Received: -2 of type <class 'int'>" in str(
+            einfo.value)
+
+    def test_check_padding_format_len_of_pad_tuple(self):
+        with pytest.raises(ValueError) as einfo:
+            PadUtils.check_padding_format((1, 1, 2))
+        assert "`padding` should have two elements. Found: (1, 1, 2)" in str(einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            PadUtils.check_padding_format((1, 1, 2, 2))
+        assert "`padding` should have two elements. Found: (1, 1, 2, 2)" in str(einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            PadUtils.check_padding_format(((1, 1, 3), (2, 2, 4)))
+        assert "`padding[0]` should have one or two elements. Found: (1, 1, 3)" in str(einfo.value)
+
+        assert PadUtils.check_padding_format(((1, 1), (2, 2))) == ((1, 1), (2, 2))
+        assert PadUtils.check_padding_format((1, 2)) == ((1, 1), (2, 2))
+        assert PadUtils.check_padding_format(1) == ((1, 1), (1, 1))
+
+    def test_check_padding_format_tuple_of_none_integer(self):
+        with pytest.raises(ValueError) as einfo:
+            PadUtils.check_padding_format((1.2, 1))
+        assert "The `1st entry of padding` argument must be a tuple of 2 integers. Received: 1.2" in str(einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            PadUtils.check_padding_format((1, 1.))
+        assert "The `2nd entry of padding` argument must be a tuple of 2 integers. Received: 1.0" in str(einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            PadUtils.check_padding_format(1.2)
+        assert "`padding` should be either an int, a tuple of 2 ints (symmetric_height_pad, symmetric_width_pad), " \
+               "or a tuple of 2 tuples of 2 ints ((top_pad, bottom_pad), (left_pad, right_pad)). Found: 1.2 of type " \
+               "<class 'float'>" in str(einfo.value)
+
+    def test_check_padding_format_tuple_of_tuple_none_integer_first(self):
+        with pytest.raises(ValueError) as einfo:
+            PadUtils.check_padding_format(((1., 2), (3, 4)))
+        assert "`padding[0]` should have one or two elements of type int. Found:(1.0, 2) " \
+               "of type [<class 'float'>, <class 'int'>]" in str(einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            PadUtils.check_padding_format(((1, 2.), (3, 4)))
+        assert "`padding[0]` should have one or two elements of type int. Found:(1, 2.0) " \
+               "of type [<class 'int'>, <class 'float'>]" in str(einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            PadUtils.check_padding_format(((1, '2'), (3, 4)))
+        assert "`padding[0]` should have one or two elements of type int. Found:(1, '2') " \
+               "of type [<class 'int'>, <class 'str'>]" in str(einfo.value)
+
+    def test_check_padding_format_tuple_of_tuple_none_integer_second(self):
+        with pytest.raises(ValueError) as einfo:
+            PadUtils.check_padding_format(((1, 2), (3., 4)))
+        assert "`padding[1]` should have one or two elements of type int. Found:(3.0, 4) " \
+               "of type [<class 'float'>, <class 'int'>]" in str(einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            PadUtils.check_padding_format(((1, 2), (3, 4.)))
+        assert "`padding[1]` should have one or two elements of type int. Found:(3, 4.0) " \
+               "of type [<class 'int'>, <class 'float'>]" in str(einfo.value)
+
+    def test_check_padding_format_valid_mix_of_int_and_tuple(self):
+        assert PadUtils.check_padding_format(((1, 2), 3)) == ((1, 2), (3, 3))
+        assert PadUtils.check_padding_format((1, (2, 3))) == ((1, 1), (2, 3))
+
+    def test_check_padding_format_invalid_mixed_tuple_and_int(self):
+        with pytest.raises(ValueError) as einfo:
+            PadUtils.check_padding_format(((1., 2), 3))
+        assert "`padding[0]` should have one or two elements of type int. Found:(1.0, 2) " \
+               "of type [<class 'float'>, <class 'int'>]" in str(einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            PadUtils.check_padding_format(((1, 2), 3.))
+        assert "The `2nd entry of padding` argument must be a tuple of 2 integers. Received: 3.0" in str(einfo.value)
+
+    def test_check_padding_format_invalid_mixed_int_and_tuple(self):
+        with pytest.raises(ValueError) as einfo:
+            PadUtils.check_padding_format((1., (2, 3)))
+        assert "The `1st entry of padding` argument must be a tuple of 2 integers. Received: 1.0" in str(einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            PadUtils.check_padding_format((1, (2., 3)))
+        assert "`padding[1]` should have one or two elements of type int. Found:(2.0, 3) " \
+               "of type [<class 'float'>, <class 'int'>]" in str(einfo.value)
+
+
+class TestReflectionPadding2D:
+
+    @pytest.fixture
+    def input_x(self):
+        return keras.Input(shape=(10, 10, 3))
+
+    def test_init_tuple_of_valid_int(self):
+        pad = (1, 3)
+        layer_name = "RefPAD"
+        ref_pad = ReflectionPadding2D(padding=pad, name=layer_name)
+        assert ref_pad.padding == ((1, 1), (3, 3))
+        assert ref_pad.name == 'RefPAD'
+        assert ref_pad.data_format == 'channels_last'
+        assert ref_pad.rank == 2
+
+        pad = (0, 1)
+        ref_pad = ReflectionPadding2D(padding=pad, name=layer_name)
+        assert ref_pad.padding == ((0, 0), (1, 1))
+        assert ref_pad.name == 'RefPAD'
+        assert ref_pad.data_format == 'channels_last'
+        assert ref_pad.rank == 2
+
+        pad = (5, 3)
+        layer_name = "RefPAD_5x3"
+        ref_pad = ReflectionPadding2D(padding=pad, name=layer_name)
+        assert ref_pad.padding == ((5, 5), (3, 3))
+
+    def test_init_tuple_of_negative_int(self):
+        with pytest.raises(ValueError) as einfo:
+            ReflectionPadding2D(padding=(-1, 1))
+        assert "The `1st entry of padding` argument must be >= 0. Received: -1 of type <class 'int'>" in str(einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            ReflectionPadding2D(padding=(1, -2))
+        assert "The `2nd entry of padding` argument must be >= 0. Received: -2 of type <class 'int'>" in str(einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            ReflectionPadding2D(padding=(-1, -2))
+        assert "The `1st entry of padding` argument must be >= 0. Received: -1 of type <class 'int'>" in str(einfo.value)
+
+    def test_init_tuple_of_invalid_format_float(self):
+        with pytest.raises(ValueError) as einfo:
+            ReflectionPadding2D(padding=(1., 1))
+        assert 'The `1st entry of padding` argument must be a tuple of 2 integers. Received: 1.0' in str(einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            ReflectionPadding2D(padding=(1, 1.2))
+        assert 'The `2nd entry of padding` argument must be a tuple of 2 integers. Received: 1.2' in str(einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            ReflectionPadding2D(padding=(1., 1.2))
+        assert 'The `1st entry of padding` argument must be a tuple of 2 integers. Received: 1.0' in str(einfo.value)
+
+    def test_init_tuple_of_invalid_format_string(self):
+        with pytest.raises(ValueError) as einfo:
+            ReflectionPadding2D(padding=('1', 2))
+        # This error message is not the best as it is missing the type information.
+        # But it is raised by keras.utils.conv_utils which I will not touch.
+        assert "`padding[0]` is str but must be int" in str(einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            ReflectionPadding2D(padding=(1, '2'))
+        assert '`padding[1]` is str but must be int' in str(einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            ReflectionPadding2D(padding=('1', '2'))
+        assert '`padding[0]` is str but must be int' in str(einfo.value)
+
+    def test_init_int(self):
+        layer_name = "RefPAD"
+        ref_pad = ReflectionPadding2D(padding=1, name=layer_name)
+        assert ref_pad.padding == ((1, 1), (1, 1))
+        assert ref_pad.name == "RefPAD"
+
+    def test_init_tuple_of_tuple_of_valid_int(self):
+        ref_pad = ReflectionPadding2D(padding=((0, 1), (2, 3)), name="RefPAD")
+        assert ref_pad.padding == ((0, 1), (2, 3))
+        assert ref_pad.name == "RefPAD"
+
+    def test_init_tuple_of_tuple_of_invalid_int(self):
+        with pytest.raises(ValueError) as einfo:
+            ReflectionPadding2D(padding=((-4, 1), (2, 3)), name="RefPAD")
+        assert "The `1st entry of padding` argument must be >= 0. Received: (-4, 1) of type <class 'tuple'>" in str(
+            einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            ReflectionPadding2D(padding=((4, -1), (2, 3)), name="RefPAD")
+        assert "The `1st entry of padding` argument must be >= 0. Received: (4, -1) of type <class 'tuple'>" in str(
+            einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            ReflectionPadding2D(padding=((4, 1), (-2, 3)), name="RefPAD")
+        assert "The `2nd entry of padding` argument must be >= 0. Received: (-2, 3) of type <class 'tuple'>" in str(
+            einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            ReflectionPadding2D(padding=((4, 1), (2, -3)), name="RefPAD")
+        assert "The `2nd entry of padding` argument must be >= 0. Received: (2, -3) of type <class 'tuple'>" in str(
+            einfo.value)
+
+    def test_init_tuple_of_tuple_of_invalid_format(self):
+        with pytest.raises(ValueError) as einfo:
+            ReflectionPadding2D(padding=((0.1, 1), (2, 3)), name="RefPAD")
+        assert "`padding[0]` should have one or two elements of type int. Found:(0.1, 1) " \
+               "of type [<class 'float'>, <class 'int'>]" in str(einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            ReflectionPadding2D(padding=(1, 2.2))
+        assert "The `2nd entry of padding` argument must be a tuple of 2 integers. Received: 2.2" in str(einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            ReflectionPadding2D(padding=((0, 1), ('2', 3)), name="RefPAD")
+        assert "`padding[1]` should have one or two elements of type int. Found:('2', 3) " \
+               "of type [<class 'str'>, <class 'int'>]" in str(einfo.value)
+
+    def test_call(self, input_x):
+        # here it behaves like a "normal" keras layer, I don't know how to test those
+        pad = (1, 0)
+        layer_name = "RefPad_3x1"
+        ref_pad = ReflectionPadding2D(padding=pad, name=layer_name)(input_x)
+        assert ref_pad.get_shape().as_list() == [None, 12, 10, 3]
+        assert ref_pad.name == 'RefPad_3x1/MirrorPad:0'
+
+
+class TestSymmerticPadding2D:
+
+    @pytest.fixture
+    def input_x(self):
+        return keras.Input(shape=(10, 10, 3))
+
+    def test_init_tuple_of_valid_int(self):
+        pad = (1, 3)
+        layer_name = "SymPad"
+        sym_pad = SymmetricPadding2D(padding=pad, name=layer_name)
+        assert sym_pad.padding == ((1, 1), (3, 3))
+        assert sym_pad.name == 'SymPad'
+        assert sym_pad.data_format == 'channels_last'
+        assert sym_pad.rank == 2
+
+        pad = (0, 1)
+        sym_pad = SymmetricPadding2D(padding=pad, name=layer_name)
+        assert sym_pad.padding == ((0, 0), (1, 1))
+        assert sym_pad.name == 'SymPad'
+        assert sym_pad.data_format == 'channels_last'
+        assert sym_pad.rank == 2
+
+        pad = (5, 3)
+        layer_name = "SymPad_5x3"
+        sym_pad = SymmetricPadding2D(padding=pad, name=layer_name)
+        assert sym_pad.padding == ((5, 5), (3, 3))
+
+    def test_init_tuple_of_negative_int(self):
+        with pytest.raises(ValueError) as einfo:
+            SymmetricPadding2D(padding=(-1, 1))
+        assert "The `1st entry of padding` argument must be >= 0. Received: -1 of type <class 'int'>" in str(
+            einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            SymmetricPadding2D(padding=(1, -2))
+        assert "The `2nd entry of padding` argument must be >= 0. Received: -2 of type <class 'int'>" in str(
+            einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            SymmetricPadding2D(padding=(-1, -2))
+        assert "The `1st entry of padding` argument must be >= 0. Received: -1 of type <class 'int'>" in str(
+            einfo.value)
+
+    def test_init_tuple_of_invalid_format_float(self):
+        with pytest.raises(ValueError) as einfo:
+            SymmetricPadding2D(padding=(1., 1))
+        assert 'The `1st entry of padding` argument must be a tuple of 2 integers. Received: 1.0' in str(einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            SymmetricPadding2D(padding=(1, 1.2))
+        assert 'The `2nd entry of padding` argument must be a tuple of 2 integers. Received: 1.2' in str(einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            SymmetricPadding2D(padding=(1., 1.2))
+        assert 'The `1st entry of padding` argument must be a tuple of 2 integers. Received: 1.0' in str(einfo.value)
+
+    def test_init_tuple_of_invalid_format_string(self):
+        with pytest.raises(ValueError) as einfo:
+            SymmetricPadding2D(padding=('1', 2))
+        # This error message is not the best as it is missing the type information.
+        # But it is raised by keras.utils.conv_utils which I will not touch.
+        assert "`padding[0]` is str but must be int" in str(einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            SymmetricPadding2D(padding=(1, '2'))
+        assert '`padding[1]` is str but must be int' in str(einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            SymmetricPadding2D(padding=('1', '2'))
+        assert '`padding[0]` is str but must be int' in str(einfo.value)
+
+    def test_init_int(self):
+        layer_name = "SymPad"
+        sym_pad = SymmetricPadding2D(padding=1, name=layer_name)
+        assert sym_pad.padding == ((1, 1), (1, 1))
+        assert sym_pad.name == "SymPad"
+
+    def test_init_tuple_of_tuple_of_valid_int(self):
+        sym_pad = SymmetricPadding2D(padding=((0, 1), (2, 3)), name="SymPad")
+        assert sym_pad.padding == ((0, 1), (2, 3))
+        assert sym_pad.name == "SymPad"
+
+    def test_init_tuple_of_tuple_of_invalid_int(self):
+        with pytest.raises(ValueError) as einfo:
+            SymmetricPadding2D(padding=((-4, 1), (2, 3)), name="SymPad")
+        assert "The `1st entry of padding` argument must be >= 0. Received: (-4, 1) of type <class 'tuple'>" in str(
+            einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            SymmetricPadding2D(padding=((4, -1), (2, 3)), name="SymPad")
+        assert "The `1st entry of padding` argument must be >= 0. Received: (4, -1) of type <class 'tuple'>" in str(
+            einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            SymmetricPadding2D(padding=((4, 1), (-2, 3)), name="SymPad")
+        assert "The `2nd entry of padding` argument must be >= 0. Received: (-2, 3) of type <class 'tuple'>" in str(
+            einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            SymmetricPadding2D(padding=((4, 1), (2, -3)), name="SymPad")
+        assert "The `2nd entry of padding` argument must be >= 0. Received: (2, -3) of type <class 'tuple'>" in str(
+            einfo.value)
+
+    def test_init_tuple_of_tuple_of_invalid_format(self):
+        with pytest.raises(ValueError) as einfo:
+            SymmetricPadding2D(padding=((0.1, 1), (2, 3)), name="SymPad")
+        assert "`padding[0]` should have one or two elements of type int. Found:(0.1, 1) " \
+               "of type [<class 'float'>, <class 'int'>]" in str(einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            SymmetricPadding2D(padding=(1, 2.2))
+        assert "The `2nd entry of padding` argument must be a tuple of 2 integers. Received: 2.2" in str(einfo.value)
+
+        with pytest.raises(ValueError) as einfo:
+            SymmetricPadding2D(padding=((0, 1), ('2', 3)), name="SymPad")
+        assert "`padding[1]` should have one or two elements of type int. Found:('2', 3) " \
+               "of type [<class 'str'>, <class 'int'>]" in str(einfo.value)
+
+    def test_call(self, input_x):
+        # here it behaves like a "normal" keras layer, I don't know how to test those
+        pad = (1, 0)
+        layer_name = "SymPad_3x1"
+        sym_pad = SymmetricPadding2D(padding=pad, name=layer_name)(input_x)
+        assert sym_pad.get_shape().as_list() == [None, 12, 10, 3]
+        assert sym_pad.name == 'SymPad_3x1/MirrorPad:0'
diff --git a/test/test_model_modules/test_inception_model.py b/test/test_model_modules/test_inception_model.py
index aa5cb284ab196d733e04a9882fa4d5a4ef639a6d..9dee30788c34cd8d1a7572947ea2e568ac2006b7 100644
--- a/test/test_model_modules/test_inception_model.py
+++ b/test/test_model_modules/test_inception_model.py
@@ -2,6 +2,8 @@ import keras
 import pytest
 
 from src.model_modules.inception_model import InceptionModelBase
+from src.model_modules.advanced_paddings import ReflectionPadding2D, SymmetricPadding2D
+from src.helpers import PyTestRegex
 
 
 class TestInceptionModelBase:
@@ -32,7 +34,8 @@ class TestInceptionModelBase:
         assert base.block_part_name() == 'a'
 
     def test_create_conv_tower_3x3(self, base, input_x):
-        opts = {'input_x': input_x, 'reduction_filter': 64, 'tower_filter': 32, 'tower_kernel': (3, 3)}
+        opts = {'input_x': input_x, 'reduction_filter': 64, 'tower_filter': 32, 'tower_kernel': (3, 3),
+                'padding': 'SymPad2D'}
         tower = base.create_conv_tower(**opts)
         # check last element of tower (activation)
         assert base.part_of_block == 1
@@ -44,12 +47,17 @@ class TestInceptionModelBase:
         conv_layer = self.step_in(act_layer)
         assert isinstance(conv_layer, keras.layers.Conv2D)
         assert conv_layer.filters == 32
-        assert conv_layer.padding == 'same'
+        assert conv_layer.padding == 'valid'
         assert conv_layer.kernel_size == (3, 3)
         assert conv_layer.strides == (1, 1)
         assert conv_layer.name == "Block_0a_3x3"
+        # check previous element of tower (padding)
+        pad_layer = self.step_in(conv_layer)
+        assert isinstance(pad_layer, SymmetricPadding2D)
+        assert pad_layer.padding == ((1, 1), (1, 1))
+        assert pad_layer.name == 'Block_0a_Pad'
         # check previous element of tower (activation)
-        act_layer2 = self.step_in(conv_layer)
+        act_layer2 = self.step_in(pad_layer)
         assert isinstance(act_layer2, keras.layers.advanced_activations.ReLU)
         assert act_layer2.name == "Block_0a_act_1"
         # check previous element of tower (conv2D)
@@ -57,7 +65,49 @@ class TestInceptionModelBase:
         assert isinstance(conv_layer2, keras.layers.Conv2D)
         assert conv_layer2.filters == 64
         assert conv_layer2.kernel_size == (1, 1)
-        assert conv_layer2.padding == 'same'
+        assert conv_layer2.padding == 'valid'
+        assert conv_layer2.name == 'Block_0a_1x1'
+        assert conv_layer2.input._keras_shape == (None, 32, 32, 3)
+
+    def test_create_conv_tower_3x3_batch_norm(self, base, input_x):
+        # import keras
+        opts = {'input_x': input_x, 'reduction_filter': 64, 'tower_filter': 32, 'tower_kernel': (3, 3),
+                'padding': 'SymPad2D', 'batch_normalisation': True}
+        tower = base.create_conv_tower(**opts)
+        # check last element of tower (activation)
+        assert base.part_of_block == 1
+        # assert tower.name == 'Block_0a_act_2/Relu:0'
+        assert tower.name == 'Block_0a_act_2_1/Relu:0'
+        act_layer = tower._keras_history[0]
+        assert isinstance(act_layer, keras.layers.advanced_activations.ReLU)
+        assert act_layer.name == "Block_0a_act_2"
+        # check previous element of tower (batch_normal)
+        batch_layer = self.step_in(act_layer)
+        assert isinstance(batch_layer, keras.layers.BatchNormalization)
+        assert batch_layer.name == 'Block_0a_BN'
+        # check previous element of tower (conv2D)
+        conv_layer = self.step_in(batch_layer)
+        assert isinstance(conv_layer, keras.layers.Conv2D)
+        assert conv_layer.filters == 32
+        assert conv_layer.padding == 'valid'
+        assert conv_layer.kernel_size == (3, 3)
+        assert conv_layer.strides == (1, 1)
+        assert conv_layer.name == "Block_0a_3x3"
+        # check previous element of tower (padding)
+        pad_layer = self.step_in(conv_layer)
+        assert isinstance(pad_layer, SymmetricPadding2D)
+        assert pad_layer.padding == ((1, 1), (1, 1))
+        assert pad_layer.name == 'Block_0a_Pad'
+        # check previous element of tower (activation)
+        act_layer2 = self.step_in(pad_layer)
+        assert isinstance(act_layer2, keras.layers.advanced_activations.ReLU)
+        assert act_layer2.name == "Block_0a_act_1"
+        # check previous element of tower (conv2D)
+        conv_layer2 = self.step_in(act_layer2)
+        assert isinstance(conv_layer2, keras.layers.Conv2D)
+        assert conv_layer2.filters == 64
+        assert conv_layer2.kernel_size == (1, 1)
+        assert conv_layer2.padding == 'valid'
         assert conv_layer2.name == 'Block_0a_1x1'
         assert conv_layer2.input._keras_shape == (None, 32, 32, 3)
 
@@ -81,7 +131,7 @@ class TestInceptionModelBase:
         tower = base.create_conv_tower(**opts)
         # check last element of tower (activation)
         assert base.part_of_block == 1
-        assert tower.name == 'Block_0a_act_1_1/Relu:0'
+        assert tower.name == 'Block_0a_act_1_2/Relu:0'
         act_layer = tower._keras_history[0]
         assert isinstance(act_layer, keras.layers.advanced_activations.ReLU)
         assert act_layer.name == "Block_0a_act_1"
@@ -89,7 +139,7 @@ class TestInceptionModelBase:
         conv_layer = self.step_in(act_layer)
         assert isinstance(conv_layer, keras.layers.Conv2D)
         assert conv_layer.filters == 32
-        assert conv_layer.padding == 'same'
+        assert conv_layer.padding == 'valid'
         assert conv_layer.kernel_size == (1, 1)
         assert conv_layer.strides == (1, 1)
         assert conv_layer.name == "Block_0a_1x1"
@@ -107,7 +157,7 @@ class TestInceptionModelBase:
         tower = base.create_pool_tower(**opts)
         # check last element of tower (activation)
         assert base.part_of_block == 1
-        assert tower.name == 'Block_0a_act_1_3/Relu:0'
+        assert tower.name == 'Block_0a_act_1_4/Relu:0'
         act_layer = tower._keras_history[0]
         assert isinstance(act_layer, keras.layers.advanced_activations.ReLU)
         assert act_layer.name == "Block_0a_act_1"
@@ -115,7 +165,7 @@ class TestInceptionModelBase:
         conv_layer = self.step_in(act_layer)
         assert isinstance(conv_layer, keras.layers.Conv2D)
         assert conv_layer.filters == 32
-        assert conv_layer.padding == 'same'
+        assert conv_layer.padding == 'valid'
         assert conv_layer.kernel_size == (1, 1)
         assert conv_layer.strides == (1, 1)
         assert conv_layer.name == "Block_0a_1x1"
@@ -124,7 +174,12 @@ class TestInceptionModelBase:
         assert isinstance(pool_layer, keras.layers.pooling.MaxPooling2D)
         assert pool_layer.name == "Block_0a_MaxPool"
         assert pool_layer.pool_size == (3, 3)
-        assert pool_layer.padding == 'same'
+        assert pool_layer.padding == 'valid'
+        # check previous element of tower(padding)
+        pad_layer = self.step_in(pool_layer)
+        assert isinstance(pad_layer, keras.layers.convolutional.ZeroPadding2D)
+        assert pad_layer.name == "Block_0a_Pad"
+        assert pad_layer.padding == ((1, 1), (1, 1))
         # check avg pool tower
         opts = {'input_x': input_x, 'pool_kernel': (3, 3), 'tower_filter': 32}
         tower = base.create_pool_tower(max_pooling=False, **opts)
@@ -132,26 +187,48 @@ class TestInceptionModelBase:
         assert isinstance(pool_layer, keras.layers.pooling.AveragePooling2D)
         assert pool_layer.name == "Block_0b_AvgPool"
         assert pool_layer.pool_size == (3, 3)
-        assert pool_layer.padding == 'same'
+        assert pool_layer.padding == 'valid'
 
     def test_inception_block(self, base, input_x):
-        conv = {'tower_1': {'reduction_filter': 64, 'tower_kernel': (3, 3), 'tower_filter': 64},
-                'tower_2': {'reduction_filter': 64, 'tower_kernel': (5, 5), 'tower_filter': 64, 'activation': 'tanh'}}
-        pool = {'pool_kernel': (3, 3), 'tower_filter': 64}
+        conv = {'tower_1': {'reduction_filter': 64,
+                            'tower_kernel': (3, 3),
+                            'tower_filter': 64, },
+                'tower_2': {'reduction_filter': 64,
+                            'tower_kernel': (5, 5),
+                            'tower_filter': 64,
+                            'activation': 'tanh',
+                            'padding': 'SymPad2D', },
+                }
+        pool = {'pool_kernel': (3, 3), 'tower_filter': 64, 'padding': ReflectionPadding2D}
         opts = {'input_x': input_x, 'tower_conv_parts': conv, 'tower_pool_parts': pool}
         block = base.inception_block(**opts)
         assert base.number_of_blocks == 1
         concatenated = block._keras_history[0].input
         assert len(concatenated) == 4
         block_1a, block_1b, block_pool1, block_pool2 = concatenated
-        assert block_1a.name == 'Block_1a_act_2/Relu:0'  # <- sometimes keras changes given name (I don't know why yet)
-        assert block_1b.name == 'Block_1b_act_2_tanh/Tanh:0'
-        assert block_pool1.name == 'Block_1c_act_1/Relu:0'
-        assert block_pool2.name == 'Block_1d_act_1/Relu:0'
+        # keras_name_part_split
+        assert block_1a.name == PyTestRegex(r'Block_1a_act_2(_\d*)?/Relu:0')
+        assert block_1b.name == PyTestRegex(r'Block_1b_act_2_tanh(_\d*)?/Tanh:0')
+        assert block_pool1.name == PyTestRegex(r'Block_1c_act_1(_\d*)?/Relu:0')
+        assert block_pool2.name == PyTestRegex(r'Block_1d_act_1(_\d*)?/Relu:0')
         assert self.step_in(block_1a._keras_history[0]).name == "Block_1a_3x3"
         assert self.step_in(block_1b._keras_history[0]).name == "Block_1b_5x5"
+        assert self.step_in(block_1a._keras_history[0], depth=2).name == 'Block_1a_Pad'
+        assert isinstance(self.step_in(block_1a._keras_history[0], depth=2), keras.layers.ZeroPadding2D)
+        assert self.step_in(block_1b._keras_history[0], depth=2).name == 'Block_1b_Pad'
+        assert isinstance(self.step_in(block_1b._keras_history[0], depth=2), SymmetricPadding2D)
+        # pooling
         assert isinstance(self.step_in(block_pool1._keras_history[0], depth=2), keras.layers.pooling.MaxPooling2D)
+        assert self.step_in(block_pool1._keras_history[0], depth=3).name == 'Block_1c_Pad'
+        assert isinstance(self.step_in(block_pool1._keras_history[0], depth=3), ReflectionPadding2D)
+
         assert isinstance(self.step_in(block_pool2._keras_history[0], depth=2), keras.layers.pooling.AveragePooling2D)
+        assert self.step_in(block_pool2._keras_history[0], depth=3).name == 'Block_1d_Pad'
+        assert isinstance(self.step_in(block_pool2._keras_history[0], depth=3), ReflectionPadding2D)
+        # check naming of concat layer
+        assert block.name == PyTestRegex('Block_1_Co(_\d*)?/concat:0')
+        assert block._keras_history[0].name == 'Block_1_Co'
+        assert isinstance(block._keras_history[0], keras.layers.merge.Concatenate)
         # next block
         opts['input_x'] = block
         opts['tower_pool_parts']['max_pooling'] = True
@@ -160,15 +237,88 @@ class TestInceptionModelBase:
         concatenated = block._keras_history[0].input
         assert len(concatenated) == 3
         block_2a, block_2b, block_pool = concatenated
-        assert block_2a.name == 'Block_2a_act_2/Relu:0'
-        assert block_2b.name == 'Block_2b_act_2_tanh/Tanh:0'
-        assert block_pool.name == 'Block_2c_act_1/Relu:0'
+        assert block_2a.name == PyTestRegex(r'Block_2a_act_2(_\d*)?/Relu:0')
+        assert block_2b.name == PyTestRegex(r'Block_2b_act_2_tanh(_\d*)?/Tanh:0')
+        assert block_pool.name == PyTestRegex(r'Block_2c_act_1(_\d*)?/Relu:0')
         assert self.step_in(block_2a._keras_history[0]).name == "Block_2a_3x3"
+        assert self.step_in(block_2a._keras_history[0], depth=2).name == "Block_2a_Pad"
+        assert isinstance(self.step_in(block_2a._keras_history[0], depth=2), keras.layers.ZeroPadding2D)
+        # block 2b
         assert self.step_in(block_2b._keras_history[0]).name == "Block_2b_5x5"
+        assert self.step_in(block_2b._keras_history[0], depth=2).name == "Block_2b_Pad"
+        assert isinstance(self.step_in(block_2b._keras_history[0], depth=2), SymmetricPadding2D)
+        # block pool
         assert isinstance(self.step_in(block_pool._keras_history[0], depth=2), keras.layers.pooling.MaxPooling2D)
+        assert self.step_in(block_pool._keras_history[0], depth=3).name == 'Block_2c_Pad'
+        assert isinstance(self.step_in(block_pool._keras_history[0], depth=3), ReflectionPadding2D)
+        # check naming of concat layer
+        assert block.name == PyTestRegex(r'Block_2_Co(_\d*)?/concat:0')
+        assert block._keras_history[0].name == 'Block_2_Co'
+        assert isinstance(block._keras_history[0], keras.layers.merge.Concatenate)
+
+    def test_inception_block_invalid_batchnorm(self, base, input_x):
+        conv = {'tower_1': {'reduction_filter': 64,
+                            'tower_kernel': (3, 3),
+                            'tower_filter': 64, },
+                'tower_2': {'reduction_filter': 64,
+                            'tower_kernel': (5, 5),
+                            'tower_filter': 64,
+                            'activation': 'tanh',
+                            'padding': 'SymPad2D', },
+                }
+        pool = {'pool_kernel': (3, 3), 'tower_filter': 64, 'padding': ReflectionPadding2D, 'max_pooling': 'yes'}
+        opts = {'input_x': input_x, 'tower_conv_parts': conv, 'tower_pool_parts': pool, }
+        with pytest.raises(AttributeError) as einfo:
+            block = base.inception_block(**opts)
+        assert "max_pooling has to be either a bool or empty. Given was: yes" in str(einfo.value)
 
     def test_batch_normalisation(self, base, input_x):
         base.part_of_block += 1
         bn = base.batch_normalisation(input_x)._keras_history[0]
         assert isinstance(bn, keras.layers.normalization.BatchNormalization)
         assert bn.name == "Block_0a_BN"
+
+    def test_padding_layer_zero_padding(self, base, input_x):
+        padding_size = ((1, 1), (0, 0))
+        zp = base.padding_layer('ZeroPad2D')
+        assert zp == keras.layers.convolutional.ZeroPadding2D
+        assert base.padding_layer('ZeroPadding2D') == keras.layers.convolutional.ZeroPadding2D
+        assert base.padding_layer(keras.layers.ZeroPadding2D) == keras.layers.convolutional.ZeroPadding2D
+        assert zp.__name__ == 'ZeroPadding2D'
+        zp_ap = zp(padding=padding_size)(input_x)
+        assert zp_ap._keras_history[0].padding == ((1, 1), (0, 0))
+
+    def test_padding_layer_sym_padding(self, base, input_x):
+        padding_size = ((1, 1), (0, 0))
+        zp = base.padding_layer('SymPad2D')
+        assert zp == SymmetricPadding2D
+        assert base.padding_layer('SymmetricPadding2D') == SymmetricPadding2D
+        assert base.padding_layer(SymmetricPadding2D) == SymmetricPadding2D
+        assert zp.__name__ == 'SymmetricPadding2D'
+        zp_ap = zp(padding=padding_size)(input_x)
+        assert zp_ap._keras_history[0].padding == ((1, 1), (0, 0))
+
+    def test_padding_layer_ref_padding(self, base, input_x):
+        padding_size = ((1, 1), (0, 0))
+        zp = base.padding_layer('RefPad2D')
+        assert zp == ReflectionPadding2D
+        assert base.padding_layer('ReflectionPadding2D') == ReflectionPadding2D
+        assert base.padding_layer(ReflectionPadding2D) == ReflectionPadding2D
+        assert zp.__name__ == 'ReflectionPadding2D'
+        zp_ap = zp(padding=padding_size)(input_x)
+        assert zp_ap._keras_history[0].padding == ((1, 1), (0, 0))
+
+    def test_padding_layer_raises(self, base, input_x):
+        with pytest.raises(NotImplementedError) as einfo:
+            base.padding_layer('FalsePadding2D')
+        assert "`'FalsePadding2D'' is not implemented as padding. " \
+               "Use one of those: i) `RefPad2D', ii) `SymPad2D', iii) `ZeroPad2D'" in str(einfo.value)
+        with pytest.raises(TypeError) as einfo:
+            base.padding_layer(keras.layers.Conv2D)
+        assert "`Conv2D' is not a valid padding layer type. Use one of those: "\
+               "i) ReflectionPadding2D, ii) SymmetricPadding2D, iii) ZeroPadding2D" in str(einfo.value)
+
+
+
+
+
diff --git a/test/test_model_modules/test_keras_extensions.py b/test/test_model_modules/test_keras_extensions.py
index 17ab4f6d65c95a5a54c9d931818f889acadef532..35188933c476157a6ba8c244d3647f7d6f8bdc59 100644
--- a/test/test_model_modules/test_keras_extensions.py
+++ b/test/test_model_modules/test_keras_extensions.py
@@ -179,6 +179,13 @@ class TestCallbackHandler:
                                                             {"name": "other_clbk", "other_clbk": "callback_other",
                                                              "path": "otherpath"}]
 
+    def test_add_callback_raise(self, clbk_handler):
+        clbk_handler.editable = False
+        with pytest.raises(PermissionError) as einfo:
+            clbk_handler.add_callback("callback_new_instance", "this_path")
+        assert 'CallbackHandler is protected and cannot be edited.' in str(einfo.value)
+
+
     def test_get_callbacks_as_dict(self, clbk_handler_with_dummies):
         clbk = clbk_handler_with_dummies
         assert clbk.get_callbacks() == [{"callback": "callback_new_instance", "path": "this_path"},
diff --git a/test/test_modules/test_experiment_setup.py b/test/test_modules/test_experiment_setup.py
index 9e6d17627d1697a2150ea7f74a373a720d2f02ac..894e4b552af4231ccc12fb85aaaebf5bbc23edf3 100644
--- a/test/test_modules/test_experiment_setup.py
+++ b/test/test_modules/test_experiment_setup.py
@@ -54,11 +54,10 @@ class TestExperimentSetup:
         assert data_store.get("experiment_name", "general") == "TestExperiment"
         path = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", "TestExperiment"))
         assert data_store.get("experiment_path", "general") == path
-        default_var_all_dict = {'o3': 'dma8eu', 'relhum': 'average_values', 'temp': 'maximum', 'u': 'average_values',
-                                'v': 'average_values', 'no': 'dma8eu', 'no2': 'dma8eu', 'cloudcover': 'average_values',
-                                'pblheight': 'maximum'}
+        default_statistics_per_var = {'o3': 'dma8eu', 'relhum': 'average_values', 'temp': 'maximum',
+                                      'u': 'average_values', 'v': 'average_values', 'no': 'dma8eu', 'no2': 'dma8eu',
+                                      'cloudcover': 'average_values', 'pblheight': 'maximum'}
         # setup for data
-        assert data_store.get("var_all_dict", "general") == default_var_all_dict
         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',
@@ -69,8 +68,8 @@ class TestExperimentSetup:
         assert data_store.get("stations", "general") == default_stations
         assert data_store.get("network", "general") == "AIRBASE"
         assert data_store.get("station_type", "general") is None
-        assert data_store.get("variables", "general") == list(default_var_all_dict.keys())
-        assert data_store.get("statistics_per_var", "general") == default_var_all_dict
+        assert data_store.get("variables", "general") == list(default_statistics_per_var.keys())
+        assert data_store.get("statistics_per_var", "general") == default_statistics_per_var
         assert data_store.get("start", "general") == "1997-01-01"
         assert data_store.get("end", "general") == "2017-12-31"
         assert data_store.get("window_history_size", "general") == 13
@@ -98,11 +97,10 @@ class TestExperimentSetup:
     def test_init_no_default(self):
         experiment_path = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "data", "testExperimentFolder"))
         kwargs = dict(parser_args={"experiment_date": "TODAY"},
-                      var_all_dict={'o3': 'dma8eu', 'relhum': 'average_values', 'temp': 'maximum'},
+                      statistics_per_var={'o3': 'dma8eu', 'relhum': 'average_values', 'temp': 'maximum'},
                       stations=['DEBY053', 'DEBW059', 'DEBW027'], network="INTERNET", station_type="background",
-                      variables=["o3", "temp"],
-                      statistics_per_var=None, start="1999-01-01", end="2001-01-01", window_history_size=4,
-                      target_var="temp", target_dim="target", window_lead_time=10, dimensions="dim1",
+                      variables=["o3", "temp"], start="1999-01-01", end="2001-01-01", window_history_size=4,
+                      target_var="relhum", target_dim="target", window_lead_time=10, dimensions="dim1",
                       interpolate_dim="int_dim", interpolate_method="cubic", limit_nan_fill=5, train_start="2000-01-01",
                       train_end="2000-01-02", val_start="2000-01-03", val_end="2000-01-04", test_start="2000-01-05",
                       test_end="2000-01-06", use_all_stations_on_all_data_sets=False, trainable=False,
@@ -120,8 +118,6 @@ class TestExperimentSetup:
                                             "TODAY_network"))
         assert data_store.get("experiment_path", "general") == path
         # setup for data
-        assert data_store.get("var_all_dict", "general") == {'o3': 'dma8eu', 'relhum': 'average_values',
-                                                             'temp': 'maximum'}
         assert data_store.get("stations", "general") == ['DEBY053', 'DEBW059', 'DEBW027']
         assert data_store.get("network", "general") == "INTERNET"
         assert data_store.get("station_type", "general") == "background"
@@ -132,7 +128,7 @@ class TestExperimentSetup:
         assert data_store.get("end", "general") == "2001-01-01"
         assert data_store.get("window_history_size", "general") == 4
         # target
-        assert data_store.get("target_var", "general") == "temp"
+        assert data_store.get("target_var", "general") == "relhum"
         assert data_store.get("target_dim", "general") == "target"
         assert data_store.get("window_lead_time", "general") == 10
         # interpolation
@@ -173,8 +169,8 @@ class TestExperimentSetup:
     def test_compare_variables_and_statistics(self):
         experiment_path = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "data", "testExperimentFolder"))
         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,
+                      statistics_per_var={'o3': 'dma8eu', 'temp': 'maximum'},
+                      stations=['DEBY053', 'DEBW059', 'DEBW027'], variables=["o3", "relhum"],
                       experiment_path=experiment_path)
         with pytest.raises(ValueError) as e:
             ExperimentSetup(**kwargs)
diff --git a/test/test_modules/test_pre_processing.py b/test/test_modules/test_pre_processing.py
index 29172a1b8500b605859e925574535c6158c7d805..d58cbd41e2ce4f25f4cd79127256e313b4aac649 100644
--- a/test/test_modules/test_pre_processing.py
+++ b/test/test_modules/test_pre_processing.py
@@ -27,7 +27,7 @@ class TestPreProcessing:
     @pytest.fixture
     def obj_with_exp_setup(self):
         ExperimentSetup(parser_args={}, stations=['DEBW107', 'DEBY081', 'DEBW013', 'DEBW076', 'DEBW087', 'DEBW001'],
-                        var_all_dict={'o3': 'dma8eu', 'temp': 'maximum'}, station_type="background")
+                        statistics_per_var={'o3': 'dma8eu', 'temp': 'maximum'}, station_type="background")
         pre = object.__new__(PreProcessing)
         super(PreProcessing, pre).__init__()
         yield pre
@@ -35,7 +35,7 @@ class TestPreProcessing:
 
     def test_init(self, caplog):
         ExperimentSetup(parser_args={}, stations=['DEBW107', 'DEBY081', 'DEBW013', 'DEBW076', 'DEBW087'],
-                        var_all_dict={'o3': 'dma8eu', 'temp': 'maximum'})
+                        statistics_per_var={'o3': 'dma8eu', 'temp': 'maximum'})
         caplog.set_level(logging.INFO)
         with PreProcessing():
             assert caplog.record_tuples[0] == ('root', 20, 'PreProcessing started')