import importlib
from random import Random
from statistics import mean, median
import sys
from typing import Counter, Union
from collections import defaultdict
import os
import pickle
from math import sqrt, ceil
from itertools import chain, islice, pairwise
import argparse
import json
import numpy as np
import pandas as pd
from matplotlib.mlab import detrend_linear
from scipy.stats import linregress, wasserstein_distance, ks_2samp
from slugify import slugify

from dhcollections.list_utils import deep_flatten
from prolothar_queue_mining.inference.queue.ras.ras import ArrivalLearningParameters, ServiceTimeLearningParameters

from prolothar_queue_mining.model.arrival_process.arrival_with_distribution import ArrivalWithDistribution
from prolothar_queue_mining.model.distribution.geometric import GeometricDistribution
from prolothar_queue_mining.model.population.list_with_replacement_population import ListWithReplacementPopulation

np.seterr(all='raise')
from prolothar_common.experiments.statistics import Statistics
import matplotlib
#not sure why this is necessary, but otherwise I sometimes observed a "Fail to allocate bitmap" exception
if matplotlib.get_backend() == 'TkAgg':
    matplotlib.use('Agg')
from matplotlib.ticker import FuncFormatter, FixedLocator
from matplotlib.colors import LinearSegmentedColormap

from statsmodels.tools.sm_exceptions import HessianInversionWarning
from statsmodels.tools.sm_exceptions import ConvergenceWarning
import warnings
warnings.filterwarnings("ignore", category=HessianInversionWarning)
warnings.filterwarnings("ignore", category=ConvergenceWarning)
from sklearn.tree import plot_tree
from sklearn.metrics import r2_score
from sklearn.feature_selection import mutual_info_regression

from tqdm import tqdm

from shsplots.matplotlib.plot_context import PlotContext, GridPlotContext
from shsplots.matplotlib.histogram import plot_histogram
from shsplots.matplotlib.density import plot_density

from experimentlogger.experiment_logger import ExperimentLogger

from prolothar_common.experiments.result_collector import ResultCollector

from prolothar_common.models.dataset import Dataset
from prolothar_common.models.dataset.transformer import OneHotEncoding
from prolothar_common.models.dataset.transformer import MinMaxScaling
from prolothar_common.models.dataset.instance import Instance
from prolothar_fdd.fodiscovery import GreedyFoDiscovery
from prolothar_subgroup_discovery import PySubgroup

from prolothar_queue_mining.model.queue import Queue
from prolothar_queue_mining.model.job import Job
from prolothar_queue_mining.model.distribution import PoissonDistribution
from prolothar_queue_mining.model.distribution import DiscreteDegenerateDistribution
from prolothar_queue_mining.model.distribution import PmfDefinedDistribution
from prolothar_queue_mining.model.distribution.plots import plot_pmf, plot_pmfs

from prolothar_queue_mining.inference.queue.times import infer_waiting_and_service_times
from prolothar_queue_mining.inference.queue import NaiveBruteForce
from prolothar_queue_mining.inference.queue import CueMin
from prolothar_queue_mining.inference.queue.kphf import Kphf
from prolothar_queue_mining.inference.sklearn.job_regression import train_random_forest_cv
from prolothar_queue_mining.inference.sklearn.train_arrival_window_departure_time_predictor import train_random_forest_cv as train_random_forest_cv_on_arrival_window
from prolothar_queue_mining.inference.sklearn.train_arrival_window_departure_time_predictor import train_linear_regression as train_linear_regression_on_arrival_window
from prolothar_queue_mining.inference.sklearn.train_arrival_window_departure_time_predictor import train_decision_tree_cv as train_decision_tree_cv_on_arrival_window
from prolothar_queue_mining.inference.queue.ras import RecurrentAdversarialServiceTimes
from prolothar_queue_mining.inference.queue.utils import count_nr_of_jobs_in_system

from prolothar_queue_mining.prediction.departure_time import DepartureTimePredictor
from prolothar_queue_mining.prediction.departure_time import OracleDepartureTimePredictor
from prolothar_queue_mining.prediction.departure_time import QueueDepartureTimePredictor
from prolothar_queue_mining.prediction.departure_time import FixedSojournTimeDepartureTimePredictor
from prolothar_queue_mining.prediction.departure_time import SklearnRegressionDepartureTimePredictor

from experiments.utils import load_dataset_with_waiting_times

AVG_SOJOURN_TIME_PREDICTOR_NAME = 'Avg Sojourn Time Predictor'
CUEMIN_NAME = 'CueMin'
KPHF_NAME = 'K-PHF'
NAIVE_BRUTE_FORCE_NAME = 'NBF'
RANDOM_FOREST_NAME = 'RF'
ORACLE_NAME = 'Oracle'
ARRIVAL_WINDOW_RANDOM_FOREST_NAME = 'AW+RF'
ARRIVAL_WINDOW_DECISION_TREE_NAME = 'AW+DT'
ARRIVAL_WINDOW_LINEAR_REGRESSION_NAME = 'AW+LR'
RAS_NAME = 'RAS'

class TimeUnit:
    def __init__(self, scaling_factor: float, unit_str: str, offset: float = 0):
        self.scaling_factor = scaling_factor
        self.unit_str = unit_str
        self.offset = offset
    def format(self, value: int) -> str:
        return f'{self.scale_value(value):.0f} {self.unit_str}'
    def divide_scaling_factor_by(self, dividend: int):
        self.scaling_factor /= dividend
    def scale_value(self, value: int) -> float:
        return (value + self.offset) / self.scaling_factor
    def inverse_scale_value(self, scaled_value: int) -> float:
        return scaled_value * self.scaling_factor - self.offset
    def to_matplotlib_formatter(self) -> FuncFormatter:
        return FuncFormatter(lambda value,pos: self.format(value))
    def to_matplotlib_major_locator(self, min_value, max_value):
        min_value_scaled = int(self.scale_value(min_value))
        max_value_scaled = ceil(self.scale_value(max_value))
        step_size = int((max_value_scaled - min_value_scaled) / 6)
        return FixedLocator([
            self.inverse_scale_value(scaled_value) for scaled_value
            in range(min_value_scaled, max_value_scaled, step_size)
        ])
    def relative(self, minimum: float) -> 'TimeUnit':
        return TimeUnit(self.scaling_factor, self.unit_str, offset=-minimum)
    def __repr__(self):
        return f'TimeUnit({self.scaling_factor}, {self.unit_str}, {self.offset})'
NoTime = TimeUnit(1, '')
Seconds = TimeUnit(1, 's')
Minutes = TimeUnit(60, 'm')
Hours = TimeUnit(3600, 'h')
Days = TimeUnit(3600*24, 'd')
Weeks = TimeUnit(3600*24*7, 'w')
TimeUnitMap = {
    time_unit.unit_str: time_unit
    for time_unit in [NoTime, Seconds, Minutes, Hours, Days, Weeks]
}

class SingleQueueInferenceExperiment:

    def __init__(
            self,
            dataset_name: str,
            categorical_features: list[str],
            numerical_features: list[str],
            nr_of_training_timesteps: int,
            test_until: int = -1,
            statistical_function_for_prediction: str = 'mean',
            methods = None,
            departure_time_unit: TimeUnit = NoTime,
            sojourn_time_unit: TimeUnit = NoTime,
            nr_of_simulation_runs_for_prediction: int = 1000,
            experiment_name: str = 'experiment',
            ground_truth_model: Queue = None,
            handcrafted_queues: list[tuple[str, Queue]] = [],
            error_analysis: bool = False,
            candidate_analysis: bool = False,
            plot_autocorrelation: bool = False,
            extrapolation_arrival_rates: list[float] = None,
            sojourn_time_outlier_threshold: int|None = None,
            record_cuemin_candidates: bool = True,
            large_train_test_gap: bool = False,
            should_evaluate_performance_spectrum: bool = False,
            plot_histograms: bool = False,
            cuemin_search_strategy: str = 'linear',
            cuemin_waiting_areas: list[str]|None = None,
            cuemin_patience: int = 30,
            limit_single_job_evaluation: int = -1,
            ras_num_steps: int = 100,
            ras_num_sequences: int = 100,
            ras_arrival_learning_parameters: ArrivalLearningParameters = ArrivalLearningParameters(),
            ras_service_time_learning_parameters: ServiceTimeLearningParameters = ServiceTimeLearningParameters(),
            ras_arrival_log_directory: str=None,
            ras_service_time_log_directory: str=None,
            seed: int = None):
        self.logger = ExperimentLogger(
            'experiments/evaluation/results/single_dataset_single_event',
            f'{experiment_name} - {dataset_name}')
        self.__seed = seed
        self.__nr_of_training_timesteps = nr_of_training_timesteps
        self.__test_until = test_until
        self.__nr_of_simulation_runs_for_prediction = nr_of_simulation_runs_for_prediction
        self.__result_collector = ResultCollector('Set')
        self.__categorical_features = categorical_features
        self.__numerical_features = numerical_features
        self.__metric_units = defaultdict(lambda: NoTime, {
            'sojourn time': sojourn_time_unit,
            'waiting time': sojourn_time_unit,
            'departure time': departure_time_unit
        })
        self.__ground_truth_model = ground_truth_model
        self.__handcrafted_queues = handcrafted_queues
        self.__error_analysis = error_analysis
        self.__candidate_analysis = candidate_analysis
        self.__methods = methods
        if extrapolation_arrival_rates is None:
            self.__extrapolation_arrival_rates = []
        else:
            self.__extrapolation_arrival_rates = extrapolation_arrival_rates
        self.__result_collector_extrapolation = ResultCollector(
            'arrival_rate_factor', mode=ResultCollector.modes.create_line_plot_mode(
                nr_of_decimals_in_text_file=3))
        self.__sojourn_time_outlier_threshold = sojourn_time_outlier_threshold
        self.__plot_autocorrelation = plot_autocorrelation
        self.__record_cuemin_candidates = record_cuemin_candidates
        self.__large_train_test_gap = large_train_test_gap
        self.__cuemin_search_strategy = cuemin_search_strategy
        self.__cuemin_waiting_areas = cuemin_waiting_areas
        self.__should_evaluate_performance_spectrum = should_evaluate_performance_spectrum
        self.__plot_histograms = plot_histograms
        self.__limit_single_job_evaluation = limit_single_job_evaluation
        self.__ras_num_steps = ras_num_steps
        self.__ras_num_sequences = ras_num_sequences
        self.__ras_arrival_learning_parameters = ras_arrival_learning_parameters
        self.__ras_service_time_learning_parameters = ras_service_time_learning_parameters
        self.__ras_arrival_log_directory = ras_arrival_log_directory
        self.__ras_service_time_log_directory = ras_service_time_log_directory
        self.__cuemin_patience = cuemin_patience
        if statistical_function_for_prediction == 'median':
            self.__statistical_function_for_prediction = median
        elif statistical_function_for_prediction == 'mean':
            self.__statistical_function_for_prediction = mean
        else:
            raise NotImplementedError(
                f'unknown statistical_function_for_prediction: {statistical_function_for_prediction}')

    def run(
            self, arrivals: list[tuple[Job, int]], departures: list[tuple[Job, int]],
            waiting_time_per_job: dict[Job, int]):
        try:
            arrivals, departures, waiting_time_per_job = self.__preprocess_data(arrivals, departures, waiting_time_per_job)
            if self.__plot_autocorrelation:
                self.__plot_autocorrelation_of_sojourn_times(arrivals, departures, 'actual')
                self.__plot_autocorrelation_of_nr_of_jobs_in_system(arrivals, departures, 'actual')
            if self.__error_analysis:
                self.__plot_arrivals_per_time_window(arrivals)
                self.__plot_inter_arrivaltimes_histogram(arrivals)
            train_arrivals, test_arrivals, train_departures, test_departures = \
                self.__train_test_split(arrivals, departures)
            arrivals = train_arrivals + test_arrivals

            if not self.__methods or ORACLE_NAME in self.__methods:
                self.logger.log('evaluate Oracle predictor')
                self.__evaluate_departure_time_predictor(
                    OracleDepartureTimePredictor({job: departure_time for job, departure_time in departures}),
                    ORACLE_NAME, arrivals, train_departures, test_departures, waiting_time_per_job)

            train_arrival_per_job = {job: time for job, time in train_arrivals}
            train_departure_per_job = {job: time for job, time in train_departures}
            self.logger.log('evaluate avg. sojourn time predictor baseline')
            self.__evaluate_departure_time_predictor(
                FixedSojournTimeDepartureTimePredictor(int(Statistics(
                    train_departure_per_job[job] - train_arrival_per_job[job]
                    for job in train_departure_per_job
                    if job in train_arrival_per_job
                ).mean())),
                AVG_SOJOURN_TIME_PREDICTOR_NAME, arrivals,
                train_departures, test_departures, waiting_time_per_job)

            self.__evaluate_ground_truth_model(
                arrivals, train_departures, test_departures, waiting_time_per_job)
            if not self.__methods or CUEMIN_NAME in self.__methods:
                self.__evaluate_cuemin(
                    arrivals, train_arrivals, train_departures, test_departures, waiting_time_per_job)
            if not self.__methods or KPHF_NAME in self.__methods:
                self.__evaluate_kphf(
                    arrivals, train_arrivals, train_departures, test_departures, waiting_time_per_job)
            if self.__methods and NAIVE_BRUTE_FORCE_NAME in self.__methods:
                self.__evaluate_naive_brute_force(
                    arrivals, train_arrivals, train_departures, test_departures, waiting_time_per_job)
            if not self.__methods or RANDOM_FOREST_NAME in self.__methods:
                self.__evaluate_random_forest(
                    arrivals, train_arrival_per_job, train_departures, test_departures, waiting_time_per_job)
            if not self.__methods or RAS_NAME in self.__methods:
                self.__evaluate_ras(
                    arrivals, train_arrivals, train_departures, test_departures, waiting_time_per_job)
            if not self.__methods or ARRIVAL_WINDOW_RANDOM_FOREST_NAME in self.__methods:
                self.__evaluate_arrival_window_random_forest(
                    arrivals, train_arrivals, train_departures, test_departures, waiting_time_per_job)
            if self.__methods and ARRIVAL_WINDOW_LINEAR_REGRESSION_NAME in self.__methods:
                self.__evaluate_arrival_window_linear_regression(
                    arrivals, train_arrivals, train_departures, test_departures, waiting_time_per_job)
            if self.__methods and ARRIVAL_WINDOW_DECISION_TREE_NAME in self.__methods:
                self.__evaluate_arrival_window_decision_tree(
                    arrivals, train_arrivals, train_departures, test_departures, waiting_time_per_job)

            self.__evaluate_handcrafted_queues(arrivals, waiting_time_per_job, train_departures, test_departures)

            self.__configure_markers_and_colors_for_result_collector()
        except Exception as e:
            self.logger.log(f'exception occured: {e}')
            raise e
        finally:
            self.logger.log('finished experiment')
            self.logger.close()

    def __preprocess_data(
            self, arrivals: list[tuple[Job, int]], departures: list[tuple[Job, int]],
            waiting_time_per_job: dict[Job, int]
    ) -> tuple[list[tuple[Job, int]], list[tuple[Job, int]], dict[Job, int]]:
        if self.__sojourn_time_outlier_threshold is not None:
            arrivals_per_job = dict(arrivals)
            sojourn_times = {
                job: departure_time - arrivals_per_job[job]
                for job, departure_time in departures
                if job in arrivals_per_job
            }
            arrivals = [
                (job, time) for job, time in arrivals
                if sojourn_times.get(job, 0) <= self.__sojourn_time_outlier_threshold
            ]
            departures = [
                (job, time) for job, time in departures
                if sojourn_times.get(job, 0) <= self.__sojourn_time_outlier_threshold
            ]
            waiting_time_per_job = dict(waiting_time_per_job)
            for job, sojourn_time in sojourn_times.items():
                if sojourn_time > self.__sojourn_time_outlier_threshold and job in waiting_time_per_job:
                    waiting_time_per_job.pop(job)
        return arrivals, departures, waiting_time_per_job

    def __plot_autocorrelation_of_sojourn_times(
            self, arrivals: list[tuple[Job, int]], departures: list[tuple[Job, int]], label: str):
        arrivals_per_job = dict(arrivals)
        actual_sojourn_times = [
            departure_time - arrivals_per_job[job]
            for job, departure_time in departures
            if job in arrivals_per_job
        ]
        with PlotContext(
                title='Autocorrelation of Sojourn Times', show=False,
                filepath=self.logger.get_absolute_path(f'autocorrelation_sojourn_times_{label}')) as plot_context:
            plot_context.get_axes().acorr(actual_sojourn_times, normed=True, detrend=detrend_linear, maxlags=10)
            plot_context.get_axes().set_xlim((-1,11))
            plot_context.get_axes().set_xlabel('lag')

    def __plot_autocorrelation_of_nr_of_jobs_in_system(
            self, arrivals: list[tuple[Job, int]], departures: list[tuple[Job, int]], label: str):
        _, y = count_nr_of_jobs_in_system(dict(arrivals), dict(departures))
        with PlotContext(
                title='Autocorrelation of Jobs in System', show=False,
                filepath=self.logger.get_absolute_path(
                    f'autocorrelation_nr_of_jobs_in_system_{label}')) as plot_context:
            plot_context.get_axes().acorr(y, normed=True, detrend=detrend_linear, maxlags=10)
            plot_context.get_axes().set_xlim((-1,11))
            plot_context.get_axes().set_xlabel('lag')

    def __evaluate_handcrafted_queues(self, arrivals, waiting_time_per_job, train_departures, test_departures):
        for name, queue in self.__handcrafted_queues:
            self.logger.log(f'evaluate manual model "{name}"')
            self.__evaluate_departure_time_predictor(
                QueueDepartureTimePredictor(
                    queue, seed=self.__seed,
                    repetitions=self.__nr_of_simulation_runs_for_prediction,
                    verbose=True),
                name, arrivals, train_departures, test_departures, waiting_time_per_job)
            self.__result_collector.set_color(name, ResultCollector.colors.GREY)
            self.__result_collector.set_marker(name, ResultCollector.markers.POINT)

    def __configure_markers_and_colors_for_result_collector(self):
        for result_collector in [self.__result_collector, self.__result_collector_extrapolation]:
            result_collector.set_color(CUEMIN_NAME, ResultCollector.colors.BLUISH_PURPLE)
            result_collector.set_marker(CUEMIN_NAME, ResultCollector.markers.PLUS)
            result_collector.set_color(AVG_SOJOURN_TIME_PREDICTOR_NAME, ResultCollector.colors.GREY)
            result_collector.set_marker(AVG_SOJOURN_TIME_PREDICTOR_NAME, ResultCollector.markers.TRIANGLE_UP)
            result_collector.set_color(ORACLE_NAME, ResultCollector.colors.GOLD)
            result_collector.set_marker(ORACLE_NAME, ResultCollector.markers.STAR)
            result_collector.set_color('Ground-Truth', ResultCollector.colors.BLACK)
            result_collector.set_marker('Ground-Truth', ResultCollector.markers.STAR)
            result_collector.set_color(NAIVE_BRUTE_FORCE_NAME, ResultCollector.colors.ORANGE)
            result_collector.set_marker(NAIVE_BRUTE_FORCE_NAME, ResultCollector.markers.X)
            result_collector.set_color(KPHF_NAME, ResultCollector.colors.RED)
            result_collector.set_marker(KPHF_NAME, ResultCollector.markers.SQUARE)
            result_collector.set_color(RANDOM_FOREST_NAME, ResultCollector.colors.GREEN)
            result_collector.set_marker(RANDOM_FOREST_NAME, ResultCollector.markers.DIAMOND)
            result_collector.set_color(ARRIVAL_WINDOW_RANDOM_FOREST_NAME, ResultCollector.colors.LIME)
            result_collector.set_marker(ARRIVAL_WINDOW_RANDOM_FOREST_NAME, ResultCollector.markers.PENTAGON)
            result_collector.set_color(RAS_NAME, ResultCollector.colors.RED)
            result_collector.set_marker(RAS_NAME, ResultCollector.markers.TRIANGLE_DOWN)
        self.__result_collector.save_plots(self.logger.get_path_to_experiment_directory())
        if self.__extrapolation_arrival_rates:
            self.__result_collector_extrapolation.save_plots(self.logger.get_path_to_experiment_directory())

    def __evaluate_ground_truth_model(
            self, arrivals, train_departures, test_departures, actual_waiting_times):
        if self.__ground_truth_model is not None:
            self.logger.log('evaluate ground-truth model')
            self.logger.log(str(self.__ground_truth_model))
            self.__evaluate_departure_time_predictor(
                QueueDepartureTimePredictor(
                    self.__ground_truth_model, seed=self.__seed,
                    repetitions=self.__nr_of_simulation_runs_for_prediction,
                    verbose=True),
                'Ground-Truth', arrivals,
                train_departures, test_departures, actual_waiting_times)

    def __evaluate_cuemin(
            self, arrivals, train_arrivals, train_departures, test_departures, actual_waiting_times):
        cuemin = CueMin(
            seed_for_distributions=self.__seed, patience=self.__cuemin_patience,
            categorical_attribute_names=self.__categorical_features,
            numerical_attribute_names=self.__numerical_features,
            record_candidates=self.__record_cuemin_candidates,
            recording_enabled=self.__record_cuemin_candidates,
            search_strategy_name=self.__cuemin_search_strategy,
            waiting_area_candidates=self.__cuemin_waiting_areas,
            verbose=True)
        self.logger.log('infer queue using CueMin')
        queue = cuemin.infer_queue(train_arrivals, train_departures)
        self.logger.log(str(queue))
        recording_dataframe = cuemin.get_recording_dataframe()
        self.__evaluate_departure_time_predictor(
            QueueDepartureTimePredictor(
                queue, seed=self.__seed,
                repetitions=self.__nr_of_simulation_runs_for_prediction,
                verbose=True),
            CUEMIN_NAME, arrivals, train_departures, test_departures, actual_waiting_times)

        if self.__error_analysis:
            assumed_service_times = list(infer_waiting_and_service_times(
                train_arrivals, train_departures, queue.get_waiting_area(), queue.get_nr_of_servers(),
                filter_delayed_jobs=True)[1].values())
            service_time_distribution_x_values = range(min(assumed_service_times), max(assumed_service_times) + 1)
            assumed_service_time_distribution = PmfDefinedDistribution.fit(assumed_service_times, laplace_smoothing=False)
            try:
                plot_pmfs(
                    [
                        assumed_service_time_distribution,
                        queue.get_servers()[0].get_service_time_definition().get_distribution()
                    ],
                    service_time_distribution_x_values,
                    show=False, filepath=self.logger.get_absolute_path('cuemin_service_time_distribution'),
                    legend_labels=['assumed', 'inferred'], mode='line')
            except AttributeError:
                plot_pmf(
                    assumed_service_time_distribution,
                    service_time_distribution_x_values,
                    show=False, filepath=self.logger.get_absolute_path('cuemin_service_time_distribution'))

        if self.__candidate_analysis:
            train_mae_list = []
            test_mae_list = []
            mdl_list = []
            self.logger.log('analyze correlation between MDL score and MAE on departure time for CueMin')
            for i, candidate in enumerate(tqdm(cuemin.get_recorded_candidates())):
                predictor = QueueDepartureTimePredictor(
                    candidate, seed=self.__seed,
                    repetitions=self.__nr_of_simulation_runs_for_prediction)
                predicted_departures = predictor.predict(arrivals)
                train_actual = {job: departure_time for job, departure_time in train_departures}
                test_actual = {job: departure_time for job, departure_time in test_departures}
                y_train_actual, y_train_predicted = zip(*[
                    (train_actual[job], predicted_departures[job])
                    for job in train_actual])
                y_test_actual, y_test_predicted = zip(*[
                    (test_actual[job], predicted_departures[job])
                    for job in test_actual])
                mdl_list.append(recording_dataframe.loc[i]['mdl_score'])
                train_mae_list.append(
                    Statistics(
                        abs(y_true - y_pred) for y_true, y_pred
                        in zip(y_train_actual, y_train_predicted)
                    ).mean()
                )
                test_mae_list.append(
                    Statistics(
                        abs(y_true - y_pred) for y_true, y_pred
                        in zip(y_test_actual, y_test_predicted)
                    ).mean()
                )
            if mdl_list:
                with PlotContext(
                        title='MDL vs MAE on departure times', show=False, use_tight_layout=True,
                        filepath=self.logger.get_absolute_path('cuemin_mdl_vs_mae.png')) as plot_context:
                    plot_context.get_axes().scatter(mdl_list, train_mae_list)
                    plot_context.get_axes().set_xlabel('L')
                    plot_context.get_axes().set_ylabel('MAE')
                    plot_context.get_axes().scatter(mdl_list, test_mae_list)
                    plot_context.get_axes().legend(['train', 'test'])
                recording_dataframe['train_MAE_departure_time'] = train_mae_list
                recording_dataframe['test_MAE_departure_time'] = test_mae_list
        recording_dataframe.to_excel(self.logger.get_absolute_path('cuemin.xlsx'))
        recording_dataframe.to_csv(self.logger.get_absolute_path('cuemin.csv'))

    def __evaluate_arrival_window_linear_regression(
            self, arrivals, train_arrivals, train_departures, test_departures, actual_waiting_times):
        self.logger.log('predict departure time with AW+LR (arrival window + linear regression)')
        departure_time_predictor = train_linear_regression_on_arrival_window(
            train_arrivals, train_departures,
            self.__categorical_features, self.__numerical_features,
            future_arrivals=0)
        self.__evaluate_departure_time_predictor(
            departure_time_predictor, ARRIVAL_WINDOW_LINEAR_REGRESSION_NAME,
            arrivals, train_departures, test_departures, actual_waiting_times)

    def __evaluate_arrival_window_decision_tree(
            self, arrivals, train_arrivals, train_departures, test_departures, actual_waiting_times):
        self.logger.log('predict departure time with AW+DT (arrival window + decision tree)')
        departure_time_predictor = train_decision_tree_cv_on_arrival_window(
            train_arrivals, train_departures,
            self.__categorical_features, self.__numerical_features,
            future_arrivals=0)
        self.logger.log('trained model')
        self.__evaluate_departure_time_predictor(
            departure_time_predictor, ARRIVAL_WINDOW_DECISION_TREE_NAME,
            arrivals, train_departures, test_departures, actual_waiting_times)

        with PlotContext(
                show=False,
                filepath=self.logger.get_absolute_path('decision_tree.pdf'),
                use_tight_layout=True) as plot_context:
            plot_context.get_figure().set_figheight(30)
            plot_context.get_figure().set_figwidth(30)
            decision_tree = departure_time_predictor.get_regression_model().best_estimator_
            plot_tree(
                decision_tree, ax=plot_context.get_axes(),
                feature_names=departure_time_predictor.get_feature_names_for_vector())

    def __evaluate_arrival_window_random_forest(
            self, arrivals, train_arrivals, train_departures, test_departures, actual_waiting_times):
        self.logger.log('predict departure time with AW+RF (arrival window + random forest)')
        departure_time_predictor = train_random_forest_cv_on_arrival_window(
            train_arrivals, train_departures,
            self.__categorical_features, self.__numerical_features,
            future_arrivals=0)
        self.logger.log('trained model')
        self.__evaluate_departure_time_predictor(
            departure_time_predictor, ARRIVAL_WINDOW_RANDOM_FOREST_NAME,
            arrivals, train_departures, test_departures, actual_waiting_times)

        random_forest = departure_time_predictor.get_regression_model().best_estimator_
        importances = random_forest.feature_importances_
        with PlotContext(
                show=False,
                filepath=self.logger.get_absolute_path('random_forest_feature_importance'),
                use_tight_layout=True) as plot_context:
            forest_importances = pd.Series(
                importances,
                index=departure_time_predictor.get_feature_names_for_vector()
            ).nlargest(n=20)
            forest_importances.plot.bar(ax=plot_context.get_axes())
            plot_context.get_axes().set_ylabel("Mean decrease in impurity")
            plot_context.get_axes().set_xticklabels(plot_context.get_axes().get_xticks(), rotation = 90)
            plot_context.get_axes().set_xticklabels(list(forest_importances.index), rotation = 90)
            plot_context.get_figure().set_figwidth(20)

    def __evaluate_kphf(
            self, arrivals, train_arrivals, train_departures, test_departures, actual_waiting_times):
        kphf = Kphf(k = 4, distribution_seed=self.__seed)
        self.logger.log('predict departure time with K-PHF')
        departure_time_predictor = kphf.infer_departure_time_predictor(train_arrivals, train_departures)
        self.logger.log('inferred K-PHF model')
        departure_time_predictor.set_show_progress_bar(True)
        departure_time_predictor.set_nr_of_repetitions(self.__nr_of_simulation_runs_for_prediction)
        self.__evaluate_departure_time_predictor(
            departure_time_predictor, KPHF_NAME,
            arrivals, train_departures, test_departures, actual_waiting_times)

    def __evaluate_ras(
            self, arrivals, train_arrivals, train_departures, test_departures, actual_waiting_times):
        ras = RecurrentAdversarialServiceTimes(
            random_seed=self.__seed,
            categorical_attribute_names=self.__categorical_features,
            numerical_attribute_names=self.__numerical_features,
            arrival_learning_parameters=self.__ras_arrival_learning_parameters,
            service_time_learning_parameters=self.__ras_service_time_learning_parameters,
            num_sequences=self.__ras_num_sequences,
            num_steps=self.__ras_num_steps,
            log_directory_arrivals=self.__ras_arrival_log_directory,
            log_directory_service_times=self.__ras_service_time_log_directory,
        )
        self.logger.log('predict departure time with RAS')
        departure_time_predictor = ras.train(train_arrivals, train_departures)
        self.logger.log('trained model')
        self.__evaluate_departure_time_predictor(
            departure_time_predictor, RAS_NAME,
            arrivals, train_departures, test_departures, actual_waiting_times)

    def __evaluate_random_forest(
            self, arrivals, train_arrival_per_job, train_departures, test_departures, actual_waiting_times):
        train_sojourn_times = {
            job: departure_time - train_arrival_per_job[job]
            for job, departure_time in train_departures
            if job in train_arrival_per_job
        }
        self.logger.log('infer random forest for departure time prediction')
        regressor, job_to_vector_transformer = train_random_forest_cv(
            train_sojourn_times, self.__categorical_features,
            self.__numerical_features)
        self.logger.log('trained model')
        predictor = SklearnRegressionDepartureTimePredictor(regressor, job_to_vector_transformer)

        self.__evaluate_departure_time_predictor(
            predictor, RANDOM_FOREST_NAME, arrivals, train_departures, test_departures, actual_waiting_times)

    def __evaluate_naive_brute_force(
            self, arrivals, train_arrivals, train_departures, test_departures, actual_waiting_times):
        naive_brute_force = NaiveBruteForce(
            list(range(1, 40)),
            batch_size_distribution_candidates=[
                PoissonDistribution(1, seed=self.__seed, shift=1),
                PoissonDistribution(3, seed=self.__seed, shift=1),
                PoissonDistribution(5, seed=self.__seed, shift=1),
                PoissonDistribution(7, seed=self.__seed, shift=1),
                PoissonDistribution(10, seed=self.__seed, shift=1),
                DiscreteDegenerateDistribution(1)
            ],
            seed=self.__seed, record_candidates=True
        )
        self.logger.log('infer queue using naive brute force search')
        queue = naive_brute_force.infer_queue(train_arrivals, train_departures)
        self.logger.log(str(queue))
        recording_dataframe = naive_brute_force.get_recording_dataframe()
        self.__evaluate_departure_time_predictor(
            QueueDepartureTimePredictor(
                queue, seed=self.__seed,
                repetitions=self.__nr_of_simulation_runs_for_prediction,
                verbose=True),
            NAIVE_BRUTE_FORCE_NAME, arrivals, train_departures, test_departures, actual_waiting_times)

        if self.__candidate_analysis:
            train_mae_list = []
            test_mae_list = []
            mdl_list = []
            self.logger.log('analyze correlation between naive brute force score and MAE on departure time')
            for i, candidate in enumerate(tqdm(naive_brute_force.get_recorded_candidates())):
                predictor = QueueDepartureTimePredictor(
                    candidate, seed=self.__seed,
                    repetitions=self.__nr_of_simulation_runs_for_prediction)
                predicted_departures = predictor.predict(arrivals)
                train_actual = {job: departure_time for job, departure_time in train_departures}
                test_actual = {job: departure_time for job, departure_time in test_departures}
                y_train_actual, y_train_predicted = zip(*[
                    (train_actual[job], predicted_departures[job])
                    for job in train_actual])
                y_test_actual, y_test_predicted = zip(*[
                    (test_actual[job], predicted_departures[job])
                    for job in test_actual])
                try:
                    mdl_list.append(recording_dataframe.loc[i]['wasserstein_distance'])
                    train_mae_list.append(
                        Statistics(
                            abs(y_true - y_pred) for y_true, y_pred
                            in zip(y_train_actual, y_train_predicted)
                        ).mean()
                    )
                    test_mae_list.append(
                        Statistics(
                            abs(y_true - y_pred) for y_true, y_pred
                            in zip(y_test_actual, y_test_predicted)
                        ).mean()
                    )
                except KeyError:
                    pass
            if mdl_list:
                with PlotContext(
                        title='Wasserstein distance vs MAE on departure times', show=False, use_tight_layout=True,
                        filepath=self.logger.get_absolute_path('naive_brute_force_wasserstein_distance_vs_mae.png')) as plot_context:
                    plot_context.get_axes().scatter(mdl_list, train_mae_list)
                    plot_context.get_axes().set_xlabel('wasserstein_distance')
                    plot_context.get_axes().set_ylabel('MAE')
                    plot_context.get_axes().scatter(mdl_list, test_mae_list)
                    plot_context.get_axes().legend(['train', 'test'])
                recording_dataframe['train_MAE_departure_time'] = train_mae_list
                recording_dataframe['test_MAE_departure_time'] = test_mae_list
        recording_dataframe.to_excel(self.logger.get_absolute_path('naive_brute_force.xlsx'))
        recording_dataframe.to_csv(self.logger.get_absolute_path('naive_brute_force.csv'))

    def __train_test_split(self, arrivals: list[tuple[Job, int]], departures: list[tuple[Job, int]]):
        train_arrivals = [(job, time) for job, time in arrivals if time <= self.__nr_of_training_timesteps]
        train_departures = [(job, time) for job, time in departures if time <= self.__nr_of_training_timesteps]
        if self.__test_until > 0:
            max_time_step = self.__test_until
        else:
            max_time_step = max(arrivals[-1][1], departures[-1][1])
        test_departures = [(job, time) for job, time in departures if self.__nr_of_training_timesteps < time <= max_time_step]
        test_arrivals = [(job, time) for job, time in arrivals if self.__nr_of_training_timesteps < time <= max_time_step]
        self.logger.log(f'training arrivals: {len(train_arrivals)}')
        self.logger.log(f'test arrivals: {len(test_arrivals)}')
        self.logger.log(f'training departures: {len(train_departures)}')
        self.logger.log(f'test departures: {len(test_departures)}')
        return train_arrivals, test_arrivals, train_departures, test_departures

    def __evaluate_departure_time_predictor(
            self, predictor: DepartureTimePredictor, predictor_name: str,
            arrivals: list[tuple[Job, int]],
            train_departures: list[tuple[Job, int]],
            test_departures: list[tuple[Job, int]],
            actual_waiting_times: dict[Job, int]):
        self.logger.log('start predicting departure times')
        if self.__large_train_test_gap:
            train_arrivals = [(job, arrival_time) for job, arrival_time in arrivals if arrival_time <= self.__nr_of_training_timesteps]
            test_arrivals = [(job, arrival_time) for job, arrival_time in arrivals if arrival_time > self.__nr_of_training_timesteps]
            predicted_waiting_times_distribution, predicted_departures_distribution = \
                predictor.predict_waiting_and_departure_times_distribution(train_arrivals)
            predicted_waiting_times_test, predicted_departures_test = predictor.predict_waiting_and_departure_times_distribution(test_arrivals)
            if predicted_waiting_times_distribution:
                predicted_waiting_times_distribution.update(predicted_waiting_times_test)
            predicted_departures_distribution.update(predicted_departures_test)
        else:
            predicted_waiting_times_distribution, predicted_departures_distribution = \
                predictor.predict_waiting_and_departure_times_distribution(arrivals)
        self.logger.log('finished predicting departure times')
        if predicted_waiting_times_distribution:
            predicted_waiting_times = {
                job: Statistics(predicted_time_list).mean() for job, predicted_time_list
                in predicted_waiting_times_distribution.items()
            }
        else:
            predicted_waiting_times = None
        predicted_departures = {
            job: self.__statistical_function_for_prediction(predicted_time_list) for job, predicted_time_list
            in predicted_departures_distribution.items()
        }
        arrivals_per_job = {job: arrival_time for job, arrival_time in arrivals}
        train_departures = {
            job: departure_time for job, departure_time in train_departures
            if job in arrivals_per_job}
        test_departures = {
            job: departure_time for job, departure_time in test_departures
            if job in arrivals_per_job}
        all_departures = {
            job: departure_time for job, departure_time in
            chain(train_departures.items(), test_departures.items())
        }
        self.__evaluate_job_predictions(
            train_departures, test_departures, predicted_departures, arrivals_per_job,
            'departure time', predictor_name)
        if actual_waiting_times and predicted_waiting_times:
            self.__evaluate_waiting_times_prediction(
                predictor_name, train_departures, test_departures,
                actual_waiting_times, predicted_waiting_times, predicted_departures,
                arrivals_per_job, predicted_waiting_times_distribution)

        self.__evaluate_nr_of_jobs_in_system(
            arrivals_per_job, all_departures, predicted_departures, predictor_name
        )

        if self.__should_evaluate_performance_spectrum:
            self.__evaluate_performance_spectrum(
                arrivals_per_job, train_departures, test_departures,
                predicted_departures, predictor_name)

        self.__evaluate_sojourn_time_predictions(
            arrivals_per_job, train_departures, test_departures,
            predicted_departures, predictor_name)

        self.__evaluate_sojourn_time_distribution_prediction(
            arrivals_per_job, train_departures, test_departures,
            predicted_departures_distribution, predictor_name)

        if self.__limit_single_job_evaluation != 0:
            self.__evaluate_single_job_predictions(
                arrivals_per_job, test_departures, actual_waiting_times,
                predicted_waiting_times_distribution, predicted_departures_distribution,
                predictor_name
            )

        if self.__plot_autocorrelation:
            if predictor_name != AVG_SOJOURN_TIME_PREDICTOR_NAME:
                self.__plot_autocorrelation_of_sojourn_times(
                    arrivals, list(predicted_departures.items()), predictor_name)
            self.__plot_autocorrelation_of_nr_of_jobs_in_system(
                arrivals, list(predicted_departures.items()), predictor_name)

        if self.__error_analysis \
        and not isinstance(predictor, OracleDepartureTimePredictor):
            self.__analyse_error_per_job({
                job: predicted_departure_time - all_departures[job]
                for job, predicted_departure_time in predicted_departures.items()
                if job in all_departures
            }, predictor_name, 'departure_time')

        if self.__extrapolation_arrival_rates:
            self.__evaluate_extrapolation_ability(predictor, predictor_name, arrivals)

    def __evaluate_single_job_predictions(
            self, arrivals_per_job: dict[Job, int], test_departures: dict[Job, int],
            actual_waiting_times: dict[Job, int]|None,
            predicted_waiting_times: dict[Job, list[int]]|None,
            predicted_departures: dict[Job, list[int]], predictor_name: str):
        x_scaler = self.__metric_units['sojourn time']
        test_departures = {
            job: departure_time for job, departure_time in test_departures.items()
            if job in arrivals_per_job and job in predicted_departures
        }
        if self.__limit_single_job_evaluation > 0:
            test_departures = {
                job: departure_time for job, departure_time in
                islice(test_departures.items(), self.__limit_single_job_evaluation)
            }
        prediction_analysis_directory = self.logger.get_absolute_path(f'predictions_{predictor_name}')
        os.mkdir(prediction_analysis_directory)
        if predicted_waiting_times:
            self.logger.log('test jobs with top-10 predicted waiting time:')
            for job, waiting_time_dist in sorted(
                    [(j,w) for j,w in predicted_waiting_times.items() if j in test_departures],
                    key=lambda x: max(x[1]), reverse=True)[:10]:
                self.logger.log(f'{job}: {round(median(waiting_time_dist))}')
        for job, actual_departure_time in test_departures.items():
            arrival_time = arrivals_per_job[job]
            actual_waiting_time = actual_waiting_times[job] if actual_waiting_times else None
            nr_of_rows = 1
            if predicted_waiting_times:
                predicted_waiting_time_distribution = Counter(predicted_waiting_times[job])
                total_count_waiting_times = sum(predicted_waiting_time_distribution.values())
                if total_count_waiting_times > 1:
                    nr_of_rows += 1
            predicted_departure_time_distribution = Counter(predicted_departures[job])
            total_count_departure_times = sum(predicted_departure_time_distribution.values())
            if total_count_departure_times > 1:
                nr_of_rows += 1
            with open(os.path.join(prediction_analysis_directory, f'{slugify(job.job_id)}.txt'), 'w') as f:
                with GridPlotContext(
                        nr_of_rows, 1, use_tight_layout=True,
                        title=f'Job {job.job_id}', show=False, sharex='all',
                        filepath=os.path.join(
                            prediction_analysis_directory,
                            f'{slugify(job.job_id)}.pdf')) as plot_context:
                    plot_context.get_axes(0, 0).vlines(
                        0, 0, 1,
                        label='a', color='green', alpha=0.5)
                    if actual_waiting_time:
                        plot_context.get_axes(0, 0).vlines(
                            x_scaler.scale_value(actual_waiting_time), 0, 1,
                            label='w', color='blue', alpha=0.5)
                        f.write(f'actual waiting time: {x_scaler.scale_value(actual_waiting_time)}\n')
                    plot_context.get_axes(0, 0).vlines(
                        x_scaler.scale_value(actual_departure_time - arrival_time), 0, 1,
                        label='d', color='red', alpha=0.5)
                    f.write(f'actual sojourn time: {x_scaler.scale_value(actual_departure_time - arrival_time)}\n')
                    if total_count_departure_times > 1:
                        x,y = zip(*[(x_scaler.scale_value(value - arrival_time), count/total_count_departure_times)
                                for value, count in predicted_departure_time_distribution.items()])
                        f.write(f'predicted sojourn time distribution x: {x}\n')
                        f.write(f'predicted sojourn time distribution y: {y}\n')
                        plot_context.get_axes(plot_context.get_nr_rows()-1, 0).vlines(
                            x, 0, y, label='__no_legend__',
                            color='red')
                        plot_context.get_axes(plot_context.get_nr_rows()-1, 0).set_ylabel('probability')
                        plot_context.get_axes(plot_context.get_nr_rows()-1, 0).set_title('predicted departure distribution')
                    plot_context.get_axes(0, 0).vlines(
                        x_scaler.scale_value(
                            round(self.__statistical_function_for_prediction(predicted_departure_time_distribution) - arrival_time)
                        ), 0, 1,
                        label='$\hat{d}$', color='red', linestyles='dashed')
                    f.write(f'predicted sojourn time: {x_scaler.scale_value(round(self.__statistical_function_for_prediction(predicted_departure_time_distribution) - arrival_time))}\n')
                    if predicted_waiting_times:
                        if total_count_waiting_times > 1:
                            x,y = zip(*[(x_scaler.scale_value(value), count/total_count_waiting_times)
                                    for value, count in predicted_waiting_time_distribution.items()])
                            f.write(f'predicted waiting time distribution x: {x}\n')
                            f.write(f'predicted waiting time distribution y: {y}\n')
                            plot_context.get_axes(1, 0).vlines(
                                x, 0, y, label='__no_legend__',
                                color='blue')
                            plot_context.get_axes(1, 0).set_ylabel('probability')
                            plot_context.get_axes(1, 0).set_title('predicted waiting time distribution')
                        plot_context.get_axes(0, 0).vlines(
                            x_scaler.scale_value(round(self.__statistical_function_for_prediction(predicted_waiting_time_distribution))), 0, 1,
                            label='$\hat{w}$', color='blue', linestyles='dashed')
                        f.write(f'predicted waiting time: {x_scaler.scale_value(round(self.__statistical_function_for_prediction(predicted_waiting_time_distribution)))}\n')
                    plot_context.get_axes(0, 0).legend(ncol=2)
                    plot_context.get_axes(0, 0).set_yticks([])
                    plot_context.get_axes(plot_context.get_nr_rows()-1, 0).set_xlabel(f'time [{x_scaler.unit_str}]')

    def __evaluate_extrapolation_ability(
            self, predictor: DepartureTimePredictor, predictor_name: str,
            arrivals: list[tuple[Job, int]]):
        if predictor_name == ORACLE_NAME:
            return
        actual_arrival_rate = len(arrivals) / (arrivals[-1][1] - arrivals[0][1])
        all_jobs = [job for job,_ in arrivals]
        seed_generator = Random(self.__seed)
        for factor in self.__extrapolation_arrival_rates:
            self.logger.log(f'run extrapolation experiment with arrival rate factor {factor}')
            arrival_process = ArrivalWithDistribution(
                ListWithReplacementPopulation(
                    all_jobs,
                    seed=seed_generator.randint(0, sys.maxsize)
                ),
                GeometricDistribution(
                    actual_arrival_rate * factor,
                    seed=seed_generator.randint(0, sys.maxsize)
                )
            )
            synthetic_arrivals = [
                (job, arrival_time)
                for arrival_time, job in [arrival_process.get_next_job()
                for _ in range(round(len(all_jobs) * factor))]
            ]
            predicted_departures = predictor.predict(synthetic_arrivals)
            self.__result_collector_extrapolation.add_result(
                predictor_name, 'mean_sojourn_time', factor,
                Statistics(
                    predicted_departures[job] - arrival_time
                    for job, arrival_time in synthetic_arrivals
                ).mean()
            )
            self.__result_collector_extrapolation.add_result(
                predictor_name, 'median_sojourn_time', factor,
                median(
                    predicted_departures[job] - arrival_time
                    for job, arrival_time in synthetic_arrivals
                )
            )
            self.__result_collector_extrapolation.add_result(
                predictor_name, 'median_inter_departure_time', factor,
                np.median(np.diff(sorted(predicted_departures.values())))
            )

    def __evaluate_waiting_times_prediction(
            self, predictor_name: str,
            train_departures: dict[Job, int], test_departures: dict[Job, int],
            actual_waiting_times: dict[Job, int], predicted_waiting_times: dict[Job, int],
            predicted_departures: dict[Job, int], arrivals_per_job: dict[Job, int],
            predicted_waiting_times_distribution: dict[Job, list[int]]):
        train_waiting_times = {
                job: waiting_time for job, waiting_time in actual_waiting_times.items()
                if job in train_departures
            }
        test_waiting_times = {
                job: waiting_time for job, waiting_time in actual_waiting_times.items()
                if job in test_departures
            }
        self.__evaluate_job_predictions(
                train_waiting_times, test_waiting_times, predicted_waiting_times,
                arrivals_per_job, 'waiting time', predictor_name)

        x_scaler = self.__metric_units['departure time']
        y_scaler = self.__metric_units['waiting time']
        actual_waiting_time_statistics_per_time_window = defaultdict(Statistics)
        predicted_waiting_time_statistics_per_time_window = defaultdict(Statistics)
        for job, departure_time in chain(train_departures.items(), test_departures.items()):
            time_window = int(x_scaler.scale_value(departure_time))
            actual_waiting_time_statistics_per_time_window[time_window].push(
                y_scaler.scale_value(actual_waiting_times[job]))
        for job, departure_time in predicted_departures.items():
            time_window = int(x_scaler.scale_value(departure_time))
            predicted_waiting_time_statistics_per_time_window[time_window].push(
                y_scaler.scale_value(predicted_waiting_times[job]))

        actual_x, actual_y = zip(*[
            (t, s.mean()) for t,s in actual_waiting_time_statistics_per_time_window.items()
        ])
        predicted_x, predicted_y = zip(*[
            (t, s.mean()) for t,s in predicted_waiting_time_statistics_per_time_window.items()
        ])

        with PlotContext(show=False, use_tight_layout=True, filepath=self.logger.get_absolute_path(
                        f'avg_waiting_times_per_time_window_{predictor_name}.pdf')) as plot_context:
            plot_context.get_axes().scatter(actual_x, actual_y)
            plot_context.get_axes().scatter(predicted_x, predicted_y)
            plot_context.get_axes().axvline(
                x_scaler.scale_value(self.__nr_of_training_timesteps), color='black')
            plot_context.get_axes().set_xlabel(f'departure time [{x_scaler.unit_str}]')
            plot_context.get_axes().set_ylabel(f'mean waiting time [{y_scaler.unit_str}]')
            plot_context.get_axes().legend(['actual', 'predicted'])
            plot_context.get_figure().set_figwidth(30)
        with PlotContext(
                show=False, use_tight_layout=True,
                log_filepath=self.logger.get_absolute_path(f'histogram_waiting_times_{predictor_name}.txt'),
                filepath=self.logger.get_absolute_path(f'histogram_waiting_times_{predictor_name}')) as plot_context:
            plot_histogram(
                list(chain(train_waiting_times.values(), test_waiting_times.values())),
                ax=plot_context.get_axes(), alpha=0.5, density=True)
            plot_histogram(
                [round(wt) for wt in deep_flatten(list(predicted_waiting_times_distribution.values()))],
                ax=plot_context.get_axes(), alpha=0.5, density=True)
            plot_context.get_axes().set_xlabel(f'waiting time')
            plot_context.get_axes().set_ylabel(f'density')
            plot_context.get_axes().legend(['actual', 'predicted'])

    def __analyse_error_per_job(
            self, error_per_job: dict[Job, float], predictor_name: str, metric_name: str):
        dataset = Dataset(self.__categorical_features, self.__numerical_features + ['error'])
        for job, error in error_per_job.items():
            features = dict(job.features)
            features['error'] = error
            dataset.add_instance(Instance(job.job_id, features))
        functional_dependency_list = GreedyFoDiscovery(beam_width=10, nr_of_results=3).run(dataset, 'error')
        self.logger.log(f'run functional dependency analysis for error on {metric_name}')
        for functional_dependency_pattern in functional_dependency_list:
            self.logger.log(f'{functional_dependency_pattern.attributes} => {functional_dependency_pattern.metrics}')

        try:
            subgroup_df = PySubgroup(result_set_size=3).discover_subgroups(dataset, 'error')
            self.logger.log(f'\n{subgroup_df.to_string()}')
        except FloatingPointError:
            self.logger.log('subgroup discovery failed due to FloatingPointError')

        self.__plot_prediction_error_vs_feature_values(error_per_job, predictor_name, metric_name)

        one_hot_encoded_dataset = OneHotEncoding().transform(dataset)
        MinMaxScaling().inplace_transform(one_hot_encoded_dataset)
        one_hot_feature_names = (
            one_hot_encoded_dataset.get_categorical_attribute_names() +
            self.__numerical_features)
        df = one_hot_encoded_dataset.to_dataframe()
        try:
            mutual_information = mutual_info_regression(
                df[one_hot_feature_names].to_numpy(),
                df['error'].to_numpy(),
                discrete_features=list(range(len(one_hot_encoded_dataset.get_categorical_attribute_names())))
            )
            self.logger.log('Top-3 features with highest mutual information on MAE:')
            for i in np.argsort(mutual_information)[-3:][::-1]:
                self.logger.log(f'{one_hot_feature_names[i]}: {mutual_information[i]}')
        except FloatingPointError:
            self.logger.log('mutual information computation failed due to FloatingPointError')

        for job, error in sorted(error_per_job.items(), key=lambda x: -abs(x[1]))[:10]:
            self.logger.log(f'{error}: {job}')

    def __plot_prediction_error_vs_feature_values(
            self, error_per_job: dict[Job, float], predictor_name: str, metric_name: str):
        for attribute in self.__numerical_features:
            x,y = zip(*[(job.features[attribute], error) for job, error in error_per_job.items()])
            with PlotContext(show=False, use_tight_layout=True, filepath=self.logger.get_absolute_path(
                    f'{metric_name}_error_{attribute}_{predictor_name}')) as plot_context:
                plot_context.get_axes().scatter(x,y)
                plot_context.get_axes().set_xlabel(attribute)
                plot_context.get_axes().set_ylabel(f'error on {metric_name}')
        for attribute in self.__categorical_features:
            grouped_error = defaultdict(list)
            for job, error in error_per_job.items():
                grouped_error[job.features[attribute]].append(error)
            categories = sorted(grouped_error.keys())
            with PlotContext(show=False, use_tight_layout=True, filepath=self.logger.get_absolute_path(
                    f'{metric_name}_error_{attribute}_{predictor_name}')) as plot_context:
                plot_context.get_axes().boxplot([grouped_error[c] for c in categories])
                plot_context.get_axes().set_xlabel(attribute)
                plot_context.get_axes().set_ylabel(f'error on {metric_name}')
                if attribute in self.__categorical_features:
                    plot_context.get_axes().set_xticklabels(
                        categories, rotation=90, ha='right')

    def __evaluate_nr_of_jobs_in_system(
            self, arrivals_per_job: dict[Job, int],
            all_departures: dict[Job, int],
            predicted_departures: dict[Job, int],
            predictor_name: str):
        x_actual, y_actual = count_nr_of_jobs_in_system(
            arrivals_per_job, all_departures)
        x_predicted, y_predicted = count_nr_of_jobs_in_system(
            arrivals_per_job, predicted_departures)

        with PlotContext(show=False, use_tight_layout=True, filepath=self.logger.get_absolute_path(
                        f'nr_of_jobs_in_system_{predictor_name}.pdf')) as plot_context:
            plot_context.get_axes().plot(x_actual, y_actual)
            plot_context.get_axes().plot(x_predicted, y_predicted)
            plot_context.get_axes().axvline(
                self.__nr_of_training_timesteps, color='black')
            plot_context.get_axes().set_xlabel(f'time')
            plot_context.get_axes().set_ylabel(f'jobs in system')
            plot_context.get_axes().legend(['actual', 'predicted', '_nolegend_'])
            plot_context.get_figure().set_figwidth(30)

        absolute_error = [abs(actual - predicted) for actual, predicted in zip(y_actual, y_predicted)]
        for error in absolute_error[:self.__nr_of_training_timesteps]:
            self.__result_collector.add_result(
                predictor_name, 'MAE jobs_in_system', 'Train', float(error))
        for error in absolute_error[self.__nr_of_training_timesteps:]:
            self.__result_collector.add_result(
                predictor_name, 'MAE jobs_in_system', 'Test', float(error))

    def __evaluate_sojourn_time_predictions(
            self, arrivals_per_job: dict[Job, int],
            train_departures: dict[Job, int],
            test_departures: dict[Job, int],
            predicted_departures: dict[Job, int],
            predictor_name: str):
        predicted_sojourn_times = {
            job: predicted_departures[job] - arrivals_per_job[job]
            for job in predicted_departures
            if predicted_departures[job] is not None
        }
        actual_sojourn_times_train = {
            job: train_departures[job] - arrivals_per_job[job]
            for job in train_departures
        }
        actual_sojourn_times_test = {
            job: test_departures[job] - arrivals_per_job[job]
            for job in test_departures
        }
        self.__evaluate_job_predictions(
            actual_sojourn_times_train, actual_sojourn_times_test, predicted_sojourn_times,
            arrivals_per_job, 'sojourn time', predictor_name)
        self.__evaluate_avg_sojourn_time_per_time_window(
            arrivals_per_job, train_departures, test_departures, predicted_departures,
            predicted_sojourn_times, predictor_name)
        self.__evaluate_median_sojourn_time_per_time_window(
            arrivals_per_job, train_departures, test_departures, predicted_departures,
            predicted_sojourn_times, predictor_name)
        with PlotContext(
                show=False, use_tight_layout=True,
                log_filepath=self.logger.get_absolute_path(f'histogram_sojourn_times_{predictor_name}.txt'),
                filepath=self.logger.get_absolute_path(f'histogram_sojourn_times_{predictor_name}')) as plot_context:
            plot_histogram(
                list(chain(actual_sojourn_times_train.values(), actual_sojourn_times_test.values())),
                ax=plot_context.get_axes(), alpha=0.5, density=True)
            plot_histogram(
                predicted_sojourn_times.values(),
                ax=plot_context.get_axes(), alpha=0.5, density=True)
            plot_context.get_axes().set_xlabel(f'sojourn time')
            plot_context.get_axes().set_ylabel(f'density')
            plot_context.get_axes().legend(['actual', 'predicted'])

    def __evaluate_sojourn_time_distribution_prediction(
            self, arrivals_per_job: dict[Job, int],
            train_departures: dict[Job, int],
            test_departures: dict[Job, int],
            predicted_departures_distribution: dict[Job, list[int]],
            predictor_name: str):
        actual_sojourn_times_train = {
            job: train_departures[job] - arrivals_per_job[job]
            for job in train_departures
        }
        actual_sojourn_times_test = {
            job: test_departures[job] - arrivals_per_job[job]
            for job in test_departures
        }
        predicted_sojourn_times = {
            job: [round(d) - arrivals_per_job[job] for d in dist]
            for job, dist in predicted_departures_distribution.items()
        }
        predicted_sojourn_times_train = {
            job: predicted_sojourn_times[job]
            for job in actual_sojourn_times_train
        }
        predicted_sojourn_times_test = {
            job: predicted_sojourn_times[job]
            for job in actual_sojourn_times_test
        }
        if predictor_name == AVG_SOJOURN_TIME_PREDICTOR_NAME:
            mean_actual_sojourn_time = Statistics(chain(
                actual_sojourn_times_train.values(),
                actual_sojourn_times_test.values())
            ).mean()
            mean_actual_sojourn_time_train = Statistics(
                actual_sojourn_times_train.values()
            ).mean()
            mean_actual_sojourn_time_test = Statistics(
                actual_sojourn_times_test.values()
            ).mean()
            self.logger.log(f'mean actual sojourn time: {mean_actual_sojourn_time}')
            self.logger.log(f'mean actual sojourn time train {mean_actual_sojourn_time_train}')
            self.logger.log(f'mean actual sojourn time test {mean_actual_sojourn_time_test}')
            median_actual_sojourn_time = median(chain(
                actual_sojourn_times_train.values(),
                actual_sojourn_times_test.values())
            )
            median_actual_sojourn_time_train = median(
                actual_sojourn_times_train.values()
            )
            median_actual_sojourn_time_test = median(
                actual_sojourn_times_test.values()
            )
            self.logger.log(f'median actual sojourn time: {median_actual_sojourn_time}')
            self.logger.log(f'median actual sojourn time train {median_actual_sojourn_time_train}')
            self.logger.log(f'median actual sojourn time test {median_actual_sojourn_time_test}')
        mean_predicted_sojourn_time = Statistics(chain(
            deep_flatten(list(predicted_sojourn_times_train.values())),
            deep_flatten(list(predicted_sojourn_times_test.values()))
        )).mean()
        mean_predicted_sojourn_time_train = Statistics(
            deep_flatten(list(predicted_sojourn_times_train.values()))
        ).mean()
        mean_predicted_sojourn_time_test = Statistics(
            deep_flatten(list(predicted_sojourn_times_test.values()))
        ).mean()
        self.logger.log(f'mean predicted sojourn time: {mean_predicted_sojourn_time}')
        self.logger.log(f'mean predicted sojourn time train {mean_predicted_sojourn_time_train}')
        self.logger.log(f'mean predicted sojourn time test {mean_predicted_sojourn_time_test}')
        median_predicted_sojourn_time = median(chain(
            deep_flatten(list(predicted_sojourn_times_train.values())),
            deep_flatten(list(predicted_sojourn_times_test.values()))
        ))
        median_predicted_sojourn_time_train = median(
            deep_flatten(list(predicted_sojourn_times_train.values()))
        )
        median_predicted_sojourn_time_test = median(
            deep_flatten(list(predicted_sojourn_times_test.values()))
        )
        self.logger.log(f'median predicted sojourn time: {median_predicted_sojourn_time}')
        self.logger.log(f'median predicted sojourn time train {median_predicted_sojourn_time_train}')
        self.logger.log(f'median predicted sojourn time test {median_predicted_sojourn_time_test}')
        self.__result_collector.add_result(
            predictor_name, 'wasserstein_distance', 'Train',
            wasserstein_distance(
                list(actual_sojourn_times_train.values()),
                deep_flatten(list(predicted_sojourn_times_train.values()))
            )
        )
        self.__result_collector.add_result(
            predictor_name, 'wasserstein_distance', 'Test',
            wasserstein_distance(
                list(actual_sojourn_times_test.values()),
                deep_flatten(list(predicted_sojourn_times_test.values()))
            )
        )
        self.__result_collector.add_result(
            predictor_name, 'KS', 'Train',
            ks_2samp(
                list(actual_sojourn_times_train.values()),
                deep_flatten(list(predicted_sojourn_times_train.values()))
            )[0]
        )
        self.__result_collector.add_result(
            predictor_name, 'KS', 'Test',
            ks_2samp(
                list(actual_sojourn_times_test.values()),
                deep_flatten(list(predicted_sojourn_times_test.values()))
            )[0]
        )

        if len(next(iter(predicted_departures_distribution.values()))) == 1 \
        or not self.__plot_histograms:
            #only one predicted value per job
            return
        with PlotContext(show=False, use_tight_layout=True, filepath=self.logger.get_absolute_path(
                        f'histogram_sojourn_times_{predictor_name}_dist')) as plot_context:
            plot_histogram(
                list(chain(actual_sojourn_times_train.values(), actual_sojourn_times_test.values())),
                ax=plot_context.get_axes(), alpha=0.5, density=True)
            plot_histogram(
                deep_flatten(list(predicted_sojourn_times.values())),
                ax=plot_context.get_axes(), alpha=0.5, density=True)
            plot_context.get_axes().set_xlabel(f'sojourn time')
            plot_context.get_axes().set_ylabel(f'density')
            plot_context.get_axes().legend(['actual', 'predicted'])
        with PlotContext(show=False, use_tight_layout=True, filepath=self.logger.get_absolute_path(
                        f'kde_sojourn_times_{predictor_name}_dist')) as plot_context:
            plot_density(
                list(chain(actual_sojourn_times_train.values(), actual_sojourn_times_test.values())),
                ax=plot_context.get_axes(), alpha=0.5, density=True)
            plot_density(
                deep_flatten(list(predicted_sojourn_times.values())),
                ax=plot_context.get_axes(), alpha=0.5, density=True, linestyle='dashed')
            plot_context.get_axes().set_xlabel(f'sojourn time')
            plot_context.get_axes().set_ylabel(f'density')
            plot_context.get_axes().legend(['actual', 'predicted'])

    def __evaluate_avg_sojourn_time_per_time_window(
            self, arrivals_per_job: dict[Job, int],
            train_departures: dict[Job, int],
            test_departures: dict[Job, int],
            predicted_departures: dict[Job, int],
            predicted_sojourn_times: dict[Job, int],
            predictor_name: str):
        x_scaler = self.__metric_units['departure time']
        y_scaler = self.__metric_units['sojourn time']
        actual_sojourn_time_statistics_per_time_window = defaultdict(Statistics)
        predicted_sojourn_time_statistics_per_time_window = defaultdict(Statistics)
        for job, departure_time in chain(train_departures.items(), test_departures.items()):
            time_window = int(x_scaler.scale_value(departure_time))
            actual_sojourn_time_statistics_per_time_window[time_window].push(
                y_scaler.scale_value(departure_time - arrivals_per_job[job]))
        for job, departure_time in predicted_departures.items():
            time_window = int(x_scaler.scale_value(departure_time))
            predicted_sojourn_time_statistics_per_time_window[time_window].push(
                y_scaler.scale_value(predicted_sojourn_times[job]))

        actual_x, actual_y, actual_c = zip(*[
            (t, s.mean(), len(s)) for t,s in actual_sojourn_time_statistics_per_time_window.items()
        ])
        predicted_x, predicted_y, predicted_c = zip(*[
            (t, s.mean(), len(s)) for t,s in predicted_sojourn_time_statistics_per_time_window.items()
        ])

        with PlotContext(show=False, use_tight_layout=True, filepath=self.logger.get_absolute_path(
                        f'avg_sojourn_times_per_time_window_{predictor_name}.pdf')) as plot_context:
            plot_context.get_axes().scatter(
                actual_x, actual_y, c=actual_c,
                cmap=LinearSegmentedColormap.from_list("", ["lightsteelblue","darkblue"]))
            plot_context.get_axes().scatter(
                predicted_x, predicted_y, c=predicted_c,
                cmap=LinearSegmentedColormap.from_list("", ["lightcoral","darkred"]))
            plot_context.get_axes().axvline(
                x_scaler.scale_value(self.__nr_of_training_timesteps), color='black')
            plot_context.get_axes().set_xlabel(f'time [{x_scaler.unit_str}]')
            plot_context.get_axes().set_ylabel(f'mean sojourn time [{y_scaler.unit_str}]')
            plot_context.get_axes().legend(['actual', 'predicted', '_nolegend_'])
            plot_context.get_figure().set_figwidth(30)

    def __evaluate_median_sojourn_time_per_time_window(
            self, arrivals_per_job: dict[Job, int],
            train_departures: dict[Job, int],
            test_departures: dict[Job, int],
            predicted_departures: dict[Job, int],
            predicted_sojourn_times: dict[Job, int],
            predictor_name: str):
        x_scaler = self.__metric_units['departure time']
        y_scaler = self.__metric_units['sojourn time']
        actual_sojourn_times_per_time_window = defaultdict(list)
        predicted_sojourn_times_per_time_window = defaultdict(list)
        for job, departure_time in chain(train_departures.items(), test_departures.items()):
            time_window = int(x_scaler.scale_value(departure_time))
            actual_sojourn_times_per_time_window[time_window].append(
                y_scaler.scale_value(departure_time - arrivals_per_job[job]))
        for job, departure_time in predicted_departures.items():
            time_window = int(x_scaler.scale_value(departure_time))
            predicted_sojourn_times_per_time_window[time_window].append(
                y_scaler.scale_value(predicted_sojourn_times[job]))

        actual_x, actual_y, actual_c = zip(*[
            (t, median(s), len(s)) for t,s in actual_sojourn_times_per_time_window.items()
        ])
        predicted_x, predicted_y, predicted_c = zip(*[
            (t, median(s), len(s)) for t,s in predicted_sojourn_times_per_time_window.items()
        ])

        with PlotContext(show=False, use_tight_layout=True, filepath=self.logger.get_absolute_path(
                        f'median_sojourn_times_per_time_window_{predictor_name}.pdf')) as plot_context:
            plot_context.get_axes().scatter(
                actual_x, actual_y, c=actual_c,
                cmap=LinearSegmentedColormap.from_list("", ["lightsteelblue","blue","darkblue"]))
            plot_context.get_axes().scatter(
                predicted_x, predicted_y, c=predicted_c,
                cmap=LinearSegmentedColormap.from_list("", ["lightcoral","red","darkred"]))
            plot_context.get_axes().axvline(
                x_scaler.scale_value(self.__nr_of_training_timesteps), color='black')
            plot_context.get_axes().set_xlabel(f'time [{x_scaler.unit_str}]')
            plot_context.get_axes().set_ylabel(f'median sojourn time [{y_scaler.unit_str}]')
            plot_context.get_axes().legend(['actual', 'predicted', '_nolegend_'])
            plot_context.get_figure().set_figwidth(30)

    def __plot_arrivals_per_time_window(self, arrivals: list[tuple[Job, int]]):
        x_scaler = self.__metric_units['departure time']
        arrivals_per_time_window = defaultdict(int)
        for job, arrival_time in arrivals:
            time_window = int(x_scaler.scale_value(arrival_time))
            arrivals_per_time_window[time_window] += 1

        actual_x, actual_y = zip(*[
            (time_window, arrival_count)
            for time_window,arrival_count in arrivals_per_time_window.items()
        ])

        with PlotContext(show=False, use_tight_layout=True, filepath=self.logger.get_absolute_path(
                        f'arrivals_per_time_window.pdf')) as plot_context:
            plot_context.get_axes().scatter(actual_x, actual_y)
            plot_context.get_axes().axvline(
                x_scaler.scale_value(self.__nr_of_training_timesteps), color='black')
            plot_context.get_axes().set_xlabel(f'arrival time [{x_scaler.unit_str}]')
            plot_context.get_axes().set_ylabel(f'arrivals')
            plot_context.get_figure().set_figwidth(30)

    def __plot_inter_arrivaltimes_histogram(self, arrivals: list[tuple[Job, int]]):
        arrival_times = [arrival_time for _,arrival_time in arrivals]
        interarrival_times = [b - a for a,b in pairwise(arrival_times)]
        with PlotContext(show=False, use_tight_layout=True, filepath=self.logger.get_absolute_path(
                        f'interarrival_times_histogram.png')) as plot_context:
            plot_histogram(interarrival_times, ax=plot_context.get_axes(), nr_of_bins=100)
            plot_context.get_axes().set_xlabel(f'interarrival time')
            plot_context.get_axes().set_ylabel(f'count')

    def __evaluate_performance_spectrum(
            self, arrivals_per_job: dict[Job, int],
            train_departures: dict[Job, int],
            test_departures: dict[Job, int],
            predicted_departures: dict[Job, int],
            predictor_name: str):
        train_departures = {
            job: departure_time for job, departure_time in train_departures.items()
            if departure_time is not None
        }
        test_departures = {
            job: departure_time for job, departure_time in test_departures.items()
            if departure_time is not None
        }
        predicted_departures = {
            job: departure_time for job, departure_time in predicted_departures.items()
            if departure_time is not None
        }
        arrivals_per_job = {
            job: arrival_time for job, arrival_time in arrivals_per_job.items()
            if job in train_departures or job in test_departures
        }

        min_arrival = min(arrivals_per_job.values())
        max_departure = max(max(test_departures.values()), max(predicted_departures.values()))
        y_value_arrival = 0
        y_value_actual = 1
        y_value_predicted = 2
        with PlotContext(
                title=f'performance spectrum actual vs predicted', show=False,
                use_tight_layout=True,
                filepath=self.logger.get_absolute_path(
                    f'performance_spectrum_actual_vs_predicted_{predictor_name}.png')) as plot_context:
            plot_context.get_axes().plot(
                [min_arrival, max_departure], [y_value_actual, y_value_actual],
                label='_nolegend_', color='black', zorder=1, linewidth=0.5)
            plot_context.get_axes().plot(
                [min_arrival, max_departure], [y_value_predicted, y_value_predicted],
                label='_nolegend_', color='black', zorder=1, linewidth=0.5)
            i = 0
            colors = ['red', 'blue', 'orange', 'purple', 'green', 'cyan', 'pink']
            for job, arrival_time in arrivals_per_job.items():
                if i == len(colors):
                    i = 0
                actual_departure_time = train_departures.get(job, test_departures.get(job, None))
                plot_context.get_axes().plot(
                    [arrival_time, actual_departure_time],
                    [y_value_arrival, y_value_actual], colors[i])
                plot_context.get_axes().plot(
                    [actual_departure_time, predicted_departures[job]],
                    [y_value_actual, y_value_predicted], colors[i])
                plot_context.get_axes().xaxis.set_major_formatter(
                    self.__metric_units['departure time'].relative(min_arrival).to_matplotlib_formatter())
                plot_context.get_axes().xaxis.set_major_locator(
                    self.__metric_units['departure time'].relative(min_arrival).to_matplotlib_major_locator(
                        min_arrival, max_departure))
                i += 1

    def __evaluate_job_predictions(
            self, train_actual: dict[Job, int], test_actual: dict[Job, int],
            prediction: dict[Job, int], arrivals_per_job: dict[Job, int],
            metric_name: str, predictor_name: str):
        y_train_actual, y_train_predicted = zip(*[
            (train_actual[job], prediction[job])
            for job in train_actual])
        y_test_actual, y_test_predicted = zip(*[
            (test_actual[job], prediction[job])
            for job in test_actual])
        minimal_value = min(min(y_train_actual), min(y_train_predicted))
        maximal_value = max(max(y_test_actual), max(y_test_predicted))

        with PlotContext(
                title=f'{metric_name} actual vs predicted', show=False,
                use_tight_layout=True,
                filepath=self.logger.get_absolute_path(
                    f'{metric_name}_actual_vs_predicted_{predictor_name}.png')) as plot_context:
            plot_context.get_axes().plot(
                [minimal_value, maximal_value],
                [minimal_value, maximal_value],
                label='_nolegend_', color='black', zorder=1, linewidth=0.5)
            plot_context.get_axes().scatter(y_train_actual, y_train_predicted, color='blue', marker='.', zorder=2)
            plot_context.get_axes().scatter(y_test_actual, y_test_predicted, color='red', marker='.', zorder=2)
            plot_context.get_axes().legend(['train', 'test'])
            plot_context.get_axes().set_xlabel('actual')
            plot_context.get_axes().set_ylabel('predicted')
            plot_context.get_axes().set_xbound(
                lower=minimal_value, upper=maximal_value)
            plot_context.get_axes().set_ybound(
                lower=minimal_value, upper=maximal_value)
            plot_context.get_axes().yaxis.set_major_formatter(
                self.__metric_units[metric_name].relative(minimal_value).to_matplotlib_formatter())
            plot_context.get_axes().yaxis.set_major_locator(
                self.__metric_units[metric_name].relative(minimal_value).to_matplotlib_major_locator(
                    minimal_value, maximal_value))
            plot_context.get_axes().xaxis.set_major_formatter(
                self.__metric_units[metric_name].relative(minimal_value).to_matplotlib_formatter())
            plot_context.get_axes().xaxis.set_major_locator(
                self.__metric_units[metric_name].relative(minimal_value).to_matplotlib_major_locator(
                    minimal_value, maximal_value))

        if self.__error_analysis and not metric_name.startswith('departure'):
            with PlotContext(
                    title=f'{metric_name} MAE dependent on arrival', show=False,
                    use_tight_layout=True, filepath=self.logger.get_absolute_path(
                        f'{metric_name}_MAE_dependent_on_arrival_{predictor_name}.png')) as plot_context:
                x = []
                y = []
                for job, arrival_time in arrivals_per_job.items():
                    try:
                        all_actual = train_actual | test_actual
                        y.append(abs(prediction[job] - all_actual[job]))
                        x.append(arrival_time)
                    except KeyError:
                        #no prediction available => do not plot
                        pass
                plot_context.get_axes().scatter(x, y)
                plot_context.get_axes().set_xlabel('arrival time')
                plot_context.get_axes().set_ylabel(f'MAE on {metric_name}')
                plot_context.get_axes().axvline(self.__nr_of_training_timesteps, color='black')
                slope, intercept, _, pvalue, _ = linregress(x, y)
                self.logger.log(f'linear regression for MAE on {metric_name} made by {predictor_name}: {slope}*T_A + {intercept}')
                self.logger.log(f'p-value for hypothesis slope=0: {pvalue}')
                self.logger.log((
                    f'f(0) = {intercept}, '
                    f'f({self.__nr_of_training_timesteps}) = {slope * self.__nr_of_training_timesteps + intercept}, '
                    f'f({max(x)}) = {slope * max(x) + intercept}'
                ))

        self.__result_collector.add_result(
            predictor_name, f'bias on {metric_name}', 'Train',
            Statistics(y_train_actual).mean() - Statistics(y_train_predicted).mean())
        self.__result_collector.add_result(
            predictor_name, f'bias on {metric_name}', 'Test',
            Statistics(y_test_actual).mean() - Statistics(y_test_predicted).mean())

        for y_true, y_pred in zip(y_train_actual, y_train_predicted):
            self.__result_collector.add_result(
                predictor_name, f'MAE {metric_name}', 'Train', abs(y_true - y_pred))
            self.__result_collector.add_result(
                predictor_name, f'MSE {metric_name}', 'Train', float((y_true - y_pred)**2))
        for y_true, y_pred in zip(y_test_actual, y_test_predicted):
            self.__result_collector.add_result(
                predictor_name, f'MAE {metric_name}', 'Test', abs(y_true - y_pred))
            self.__result_collector.add_result(
                predictor_name, f'MSE {metric_name}', 'Test', float((y_true - y_pred)**2))

        self.__result_collector.add_result(
            predictor_name, f'RMSE {metric_name}', 'Train',
            sqrt(Statistics((y_true - y_pred)**2 for y_true, y_pred in zip(y_train_actual, y_train_predicted)).mean()))
        self.__result_collector.add_result(
            predictor_name, f'RMSE {metric_name}', 'Test',
            sqrt(Statistics(abs(y_true - y_pred)**2 for y_true, y_pred in zip(y_test_actual, y_test_predicted)).mean()))

        self.__result_collector.add_result(
            predictor_name, f'R2 {metric_name}', 'Train',
           r2_score(y_train_actual, y_train_predicted))
        self.__result_collector.add_result(
            predictor_name, f'R2 {metric_name}', 'Test',
            r2_score(y_test_actual, y_test_predicted))

def load_ground_truth_model(data_folder: str) -> Union[Queue, None]:
    model_file = os.path.join(data_folder, 'model.pickle')
    if os.path.exists(model_file):
        with open(model_file, 'rb') as f:
            return pickle.load(f)

def load_handcrafted_queues(handcrafted_queues: list[str]) -> list[tuple[str, Queue]]:
    queue_list = []
    for file_name in handcrafted_queues:
        module = importlib.import_module(f'experiments.handcrafted_queues.{file_name}')
        queue_list.append((module.name, module.queue))
    return queue_list

def parse_ras_arrival_learning_parameters() -> ArrivalLearningParameters:
    return ArrivalLearningParameters(
        batch_size=arguments.ras_arrival_batch_size,
        num_epochs=arguments.ras_arrival_num_epochs,
        hidden_units=arguments.ras_arrival_hidden_units,
        learning_rate=arguments.ras_arrival_learning_rate
    )

def parse_ras_service_time_learning_parameters() -> ServiceTimeLearningParameters:
    return ServiceTimeLearningParameters(
        batch_size=arguments.ras_service_time_batch_size,
        num_epochs=arguments.ras_service_time_num_epochs,
        hidden_units=arguments.ras_service_time_hidden_units,
        learning_rate_gan=arguments.ras_learning_rate_gan,
        generator_dimensions=arguments.ras_generator_dimensions,
        n_discriminator_updates=arguments.ras_n_discriminator_updates
    )

if __name__ == '__main__':
    parser = argparse.ArgumentParser('queue mining experiments on single dataset')

    parser.add_argument(
        '--experiment_name', type=str,
        help='name of the experiment',
        required=False, default='experiment')

    parser.add_argument(
        '--data_folder', type=str,
        help='path to the folder containing the data',
        required=True)

    parser.add_argument(
        '--categorical_features', type=str,
        help='name of the attributes containing categorical values',
        nargs='*', default=[])

    parser.add_argument(
        '--numerical_features', type=str,
        help='name of the attributes containing numerical values',
        nargs='*', default=[])

    parser.add_argument(
        '--train_until', type=int,
        help=('number of timesteps for model inference. all jobs arrived after '
              'this are used for testset evaluation'),
        required=True)
    parser.add_argument(
        '--test_until', type=int,
        help=('all jobs arrived after this timestamp will not be used for testset evaluation'),
        required=False, default=-1)

    parser.add_argument(
        '--nr_of_simulation_runs_for_prediction', type=int,
        help=('number of simulations to average over for prediction'),
        required=False, default=10)

    parser.add_argument(
        '--departure_time_unit', type=str,
        help=('unit for the departure time, which is used in all plots'),
        required=False, default='', choices=['', 's', 'm', 'h', 'd', 'w'])

    parser.add_argument(
        '--sojourn_time_unit', type=str,
        help=('unit for the sojourn time, which is used in all plots'),
        required=False, default='', choices=['', 's', 'm', 'h', 'd', 'w'])

    parser.add_argument(
        '--time_factor', type=int,
        help='time_factor that was used to create this dataset => how many timesteps build 1 second',
        required=False, default=1)

    parser.add_argument(
        '--methods', type=str, nargs='*',
        help='list of methods for comparison',
        required=False, default=[])

    parser.add_argument(
        '--handcrafted_queues', type=str, nargs='*',
        help='a list of file names for handcrafted queue models',
        required=False, default=[])

    parser.add_argument(
        '--statistical_function_for_prediction', type=str,
        help='determines how to aggregate predictions from the predicted distribution',
        required=False, default='median', choices=['mean', 'median'])

    parser.add_argument(
        '--candidate_analysis', action='store_true',
        help=('add this flag for intensive analysis on candidates generated during search - '
              'e.g. for candidates generated in CueMin'),
        required=False, default=False)

    parser.add_argument(
        '--error_analysis', action='store_true',
        help=('add this flag for intensive analysis '
              '- e.g. functional dependency or subgroup discovery - on prediction errors'),
        required=False, default=False)

    parser.add_argument(
        '--plot_autocorrelation', action='store_true',
        help=('add this flag to plot autocorrelation of actual and predicted values'),
        required=False, default=False)

    parser.add_argument(
        '--record_cuemin_candidates', action='store_true',
        help=('add this flag to record generated candidates in CueMin'),
        required=False, default=False)

    parser.add_argument(
        '--large_train_test_gap', action='store_true',
        help=(
            'add this flag if there is a large time gap between train and test jobs, '
            'which is normally because train and test set are a completely '
            'indepdendent time range. if queue predictors expect a service batch size,'
            'they would predict very high sojourn times for jobs at the end of the '
            'training time range'
        ),
        required=False, default=False)

    parser.add_argument(
        '--evaluate_performance_spectrum', action='store_true',
        help=(
            'add this flag if you want to plot the performance spectrum for '
            'actual and predicted departure times'
        ),
        required=False, default=False)

    parser.add_argument(
        '--plot_histograms', action='store_true',
        help=(
            'add this flag if you want to plot histograms '
            'during evaluation (e.g. for sojourn time distribution)'
        ),
        required=False, default=False)

    parser.add_argument(
        '--extrapolation_arrival_rates', type=float,
        help='a list of arrival rates for extrapolation experiment',
        nargs='*', default=[])

    parser.add_argument(
        '--seed', type=int,
        help='random generator seed',
        required=False, default=None)

    parser.add_argument(
        '--limit_single_job_evaluation', type=int,
        help=(
            'specifies how many single job predictions are evaluated. '
            'set to 0 to turn off evaluation. '
            '-1 means no limit.'
        ),
        required=False, default=10)

    parser.add_argument(
        '--sojourn_time_outlier_threshold', type=int,
        help='if set, jobs with an actual sojourn time above this threshold will be filtered out',
        required=False, default=None)

    parser.add_argument(
        '--cuemin_search_strategy', type=str,
        help=(
            'use default (adaptive) unless results are not good. then use "linear", which can take a long time'
        ), required=False, default='adaptive')
    parser.add_argument(
        '--cuemin_waiting_areas', type=str, nargs='*',
        help='controls the list of possible waiting areas for cuemin',
        required=False, default=None, choices=['FCFS', 'LCFS', 'PQ-c', 'LR', 'FLIFO'])
    parser.add_argument(
        '--cuemin_cpus_sklearn', type=int,
        help='controls the number of cpus for sklearn algorithms used in cuemin',
        required=False, default=1)
    parser.add_argument(
        '--cuemin_patience', type=int,
        help='controls the patience parameter of cuemin',
        required=False, default=10)

    parser.add_argument(
        '--ras_num_sequences', type=int,
        help='num_sequences parameter in RAS',
        required=False, default=100)
    parser.add_argument(
        '--ras_num_steps', type=int,
        help='num_steps parameter in RAS',
        required=False, default=100)

    parser.add_argument(
        '--ras_arrival_batch_size', type=int,
        required=False, default=32)
    parser.add_argument(
        '--ras_arrival_num_epochs', type=int,
        required=False, default=100)
    parser.add_argument(
        '--ras_arrival_hidden_units', type=int,
        required=False, default=32)
    parser.add_argument(
        '--ras_arrival_learning_rate', type=float,
        required=False, default=0.00001)
    parser.add_argument(
        '--ras_arrival_log_directory', type=str,
        required=False, default=None)

    parser.add_argument(
        '--ras_service_time_batch_size', type=int,
        required=False, default=32)
    parser.add_argument(
        '--ras_service_time_num_epochs', type=int,
        required=False, default=100)
    parser.add_argument(
        '--ras_service_time_hidden_units', type=int,
        required=False, default=32)
    parser.add_argument(
        '--ras_learning_rate_gan', type=float,
        required=False, default=0.0001)
    parser.add_argument(
        '--ras_generator_dimensions', type=int,
        required=False, default=50)
    parser.add_argument(
        '--ras_n_discriminator_updates', type=int,
        required=False, default=10)
    parser.add_argument(
        '--ras_service_time_log_directory', type=str,
        required=False, default=None)

    arguments = parser.parse_args()

    for time_unit in TimeUnitMap.values():
        time_unit.divide_scaling_factor_by(arguments.time_factor)

    experiment = SingleQueueInferenceExperiment(
        os.path.basename(arguments.data_folder),
        arguments.categorical_features,
        arguments.numerical_features,
        arguments.train_until,
        test_until=arguments.test_until,
        methods=arguments.methods,
        experiment_name=arguments.experiment_name,
        seed=arguments.seed,
        departure_time_unit=TimeUnitMap[arguments.departure_time_unit],
        sojourn_time_unit=TimeUnitMap[arguments.sojourn_time_unit],
        statistical_function_for_prediction = arguments.statistical_function_for_prediction,
        nr_of_simulation_runs_for_prediction=arguments.nr_of_simulation_runs_for_prediction,
        ground_truth_model=load_ground_truth_model(arguments.data_folder),
        handcrafted_queues=load_handcrafted_queues(arguments.handcrafted_queues),
        candidate_analysis=arguments.candidate_analysis,
        plot_autocorrelation=arguments.plot_autocorrelation,
        error_analysis=arguments.error_analysis,
        limit_single_job_evaluation=arguments.limit_single_job_evaluation,
        extrapolation_arrival_rates=arguments.extrapolation_arrival_rates,
        sojourn_time_outlier_threshold=arguments.sojourn_time_outlier_threshold,
        large_train_test_gap=arguments.large_train_test_gap,
        should_evaluate_performance_spectrum=arguments.evaluate_performance_spectrum,
        plot_histograms=arguments.plot_histograms,
        record_cuemin_candidates=arguments.record_cuemin_candidates,
        cuemin_search_strategy=arguments.cuemin_search_strategy,
        cuemin_waiting_areas=arguments.cuemin_waiting_areas,
        cuemin_patience=arguments.cuemin_patience,
        ras_num_steps=arguments.ras_num_steps,
        ras_num_sequences=arguments.ras_num_sequences,
        ras_arrival_learning_parameters=parse_ras_arrival_learning_parameters(),
        ras_service_time_learning_parameters=parse_ras_service_time_learning_parameters(),
        ras_arrival_log_directory=arguments.ras_arrival_log_directory,
        ras_service_time_log_directory=arguments.ras_service_time_log_directory
    )
    with experiment.logger.open_text_file('arguments.json') as f:
        f.write(json.dumps(vars(arguments), indent=4, sort_keys=True))
    experiment.run(*load_dataset_with_waiting_times(arguments.data_folder))