/cmdp_collider_config.py

# -*- author: Dimitri Scheftelowitsch -*-
# -*- coding:utf-8 -*-

"""
A config file for collider

Generates a sequence of CMDPs and runs some not very sophisticated analysis on
them.
"""

import numpy as np
import pandas as pd
import matplotlib
from typing import Mapping, Any
from scipy.stats import t

matplotlib.use('PDF')
from matplotlib import pyplot as plt

experiment_name = 'CMDP'

values = {
    'scenarios': np.array([2, 3, 5]),
    'states': np.array([2, 5, 10, 20, 50, 100, 200, 300]),  # 10, 20, 50, 100, 200, 300
    'actions': np.array([2, 3, 5]),
    'run': np.arange(0, 30, 1), # np.arange(0, 30, 1),
    'gamma': np.array([0.9, 0.999]),
    'mode': ['deterministic', 'random'],
}

stages = [
    {
        'name': 'generate',
        'exe':  'python3',
        'args': ['/home/scheftel/projects/bmdp-in-python/cmdp_generator.py',
                 '--mode={mode}',
                 '--scenarios={scenarios}',
                 '--states={states}',
                 '--actions={actions}']
    },
    {
        'name': 'runIMGP',
        'exe':  '/home/scheftel/projects/BMDP/BMDPanalysis',
        'args': ['{generate}', '901', '{gamma}'],
        'timeout': 6050
    },
    {
        'name': 'runIMPP',
        'exe': '/home/scheftel/projects/BMDP/BMDPanalysis',
        'args': ['{generate}', '902', '{gamma}'],
        'timeout': 6050
    },
    {
        'name': 'runMIP',
        'exe': '/home/scheftel/projects/BMDP/BMDPanalysis',
        'args': ['{generate}', '904', '{gamma}'],
        'timeout': 6050
    },
    {
        'name': 'runNLP',
        'exe': '/home/scheftel/projects/BMDP/BMDPanalysis',
        'args': ['{generate}', '905', '{gamma}'],
        'timeout': 6050
    },
    {
        'name': 'runQCLP',
        'exe': '/home/scheftel/projects/BMDP/BMDPanalysis',
        'args': ['{generate}', '903', '{gamma}'],
        'timeout': 6010
    },
    {
        'name': 'debug',
        'exe': 'echo',
        'args': ['mode: {mode}, scenarios: {scenarios}, states: {states}, actions: {actions}, run: {run}']
    },
]


def postprocess(data: pd.DataFrame, values: Mapping[str, Any]):
    algorithms = ["IMGP", "IMPP", "QCLP", "MIP", "NLP"]
    major_tuples = ["scenarios", "states", "actions"]
    valid_run_filter = data['run'].map(lambda r: r in values['run'])
    data.loc[:, 'value'] = pd.Series(0.0, index=data.index)
    data.loc[:, 'relative_deviation'] = pd.Series(0.0, index=data.index)
    data.sort_values(by=major_tuples, inplace=True)

    models = values['mode']
    gammas = values['gamma']

    def filter_result(row):
        result = None
        run_stages = map("run{0}".format, algorithms)
        if row.Stage in run_stages:
            lines = row.output.decode('utf-8').split("\n")
            if len(lines) >= 2:
                try:
                    result = float(lines[-2])
                except ValueError:
                    result = None
        return result

    def compute_deviation(row):
        result = None
        run_stages = map("run{0}".format, algorithms)
        if (row.Stage in run_stages) and (not np.isnan(row.value)):
            query = (data.run == row.run) & (data.states == row.states) \
                & (data.actions == row.actions) & (data.scenarios == row.scenarios) \
                & (data['mode'] == row['mode']) & (data.gamma == row.gamma) \
                & ~np.isnan(data.value) & (data.value > 0.0)
            impp_value = data[query & (data.Stage == 'runIMPP')].value.iloc[0]
            cutoff = 2 * impp_value
            max_value = max(data[query & (data.value <= cutoff)].value)
            result = 1.0 - row.value / max_value
        return result

    data.loc[:, 'value'] = data.apply(filter_result, axis=1)

    # THIS CODE IS DUE TO A BUG IN GSL WHERE I HAD TO RERUN ONE SPECIFIC VALUE
    #data.loc[14051, 'value'] = 9.43569e+03
    #data.loc[14051, 'cpu_user'] = 184.74
    # END UGLY HACK

    data.loc[:, 'relative_deviation'] = data.apply(compute_deviation, axis=1)
    data_filtered = data[valid_run_filter]

    for model in models:
        for gamma in gammas:
            # tex_table_header = "\\begin{table}\\begin{tabular}"
            tex_code = ""
            tuples = data_filtered.loc[:, major_tuples].drop_duplicates()
            for _, tuple in tuples.iterrows():
                data_tuple = data_filtered[(data_filtered.scenarios == tuple.scenarios) &
                                           (data_filtered.states == tuple.states) &
                                           (data_filtered.actions == tuple.actions) &
                                           (data_filtered['mode'] == model) & (data_filtered.gamma == gamma)]
                avg_time = {}
                relative_deviation = {}
                errors = {}
                tex_code += "\\num{{{0}}} & \\num{{{1}}} & \\num{{{2}}} \n".format(tuple.scenarios, tuple.states,
                                                                                     tuple.actions)
                # find all runs that can be compared
                for algorithm in algorithms:
                    stage = "run{0}".format(algorithm)
                    invalid_value_filter = (np.isnan(data_tuple.value) | (data_tuple.value == 0.0))

                    successful_runs = data_tuple[~invalid_value_filter & (data_tuple.Stage == stage)]
                    errors[algorithm] = len(data_tuple[invalid_value_filter & (data_tuple.Stage == stage)])
                    avg_time[algorithm] = successful_runs.cpu_user.mean()

                    for run in values['run']:
                        data_run = data_tuple[data_tuple.run == run]
                        # this is an ugly hack to ensure valid values.
                        cutoff = data_run[data_run.Stage == 'runIMPP'].value.iloc[0]
                        cutoff = 2 * cutoff
                        validity_criterion = ~(np.isnan(data_run.value) | (data_run.value > cutoff))

                        try:
                            max_value = max(data_run[validity_criterion].value)
                        except ValueError:
                            max_value = None
                        data_algorithm_program = data_run[data_run.Stage == stage]
                        if len(data_algorithm_program) > 0:
                            program_value = data_algorithm_program.value.iloc[0]
                            if program_value > cutoff:
                                errors[algorithm] += 1
                            if not (np.isnan(program_value) or (program_value == 0.0) or (program_value > cutoff)):
                                deviation = 1.0 - program_value / max_value
                                if np.isnan(deviation):
                                    raise ArithmeticError
                                relative_deviation[algorithm] = relative_deviation.get(algorithm, 0) + deviation
                                data_query = (data.Stage == stage) & (data.run == run) \
                                    & (data.scenarios == tuple.scenarios) \
                                    & (data.states == tuple.states) \
                                    & (data.actions == tuple.actions) \
                                    & (data.gamma == gamma) \
                                    & (data['mode'] == model)
                                #data[data_query].relative_deviation = deviation

                    if algorithm in relative_deviation:
                        relative_deviation[algorithm] /= len(successful_runs)
                        if np.isnan(relative_deviation[algorithm]):
                            raise ArithmeticError
                    else:
                        relative_deviation[algorithm] = None

                    relative_deviation_string = "N/A"
                    if relative_deviation[algorithm] not in [np.nan, None]:
                        relative_deviation_string = "\\num{{{0:.2f}}}\%".format(relative_deviation[algorithm] * 100)

                    average_time_string = "N/A"
                    if not np.isnan(avg_time[algorithm]):
                        average_time_string = "\\num{{{0:.3f}}}\\si{{s}}".format(avg_time[algorithm])

                    tex_code += "& {0} & \\num{{{1}}} & {2} \n".format(average_time_string, errors[algorithm],
                                                                     relative_deviation_string)
                tex_code += "\\\\ \n"
            file = "CMDP_results_{0}_{1}.tex".format(model, gamma)
            f = open(file=file, mode='w+b')
            f.write(bytes(tex_code, 'utf-8'))
            f.close()

    pure_comparison = {
        'exact': 'MIP',
        'heuristic': 'IMPP',
        'exact_color': 'c',
        'heuristic_color': 'm',
        'filename': 'pure'
    }

    stationary_comparison = {
        'exact': 'NLP',
        'heuristic': 'IMGP',
        'exact_color': 'b',
        'heuristic_color': 'k',
        'filename': 'stationary'
    }

    plots = [pure_comparison, stationary_comparison]
    plot_axes = ['states', 'scenarios']
    subplots = ['cpu_user', 'relative_deviation']
    subplot_descriptions = {
        'cpu_user': 'CPU time',
        'relative_deviation': 'Deviation from optimum'
    }

    confidence_level = 0.99
    for plot in plots:
        for key in plot_axes:
            for subplot in subplots:
                filename = "plot_{0}_{1}_{2}.pdf".format(plot['filename'], key, subplot)
                x = values[key]
                heuristic_mean = np.zeros(values[key].shape)
                exact_mean = np.zeros(values[key].shape)
                heuristic_std = np.zeros(values[key].shape)
                exact_std = np.zeros(values[key].shape)
                for (i, value) in enumerate(values[key]):
                    query_heuristic = (data_filtered.Stage == "run{0}".format(plot['heuristic']))
                    query_exact = (data_filtered.Stage == "run{0}".format(plot['exact']))
                    query_value = (data_filtered[key] == value) & ~np.isnan(data_filtered[subplot])
                    data_key = data_filtered[query_value]
                    exact_mean[i] = data_key[query_exact][subplot].mean()
                    heuristic_mean[i] = data_key[query_heuristic][subplot].mean()
                    num_samples = len(data_key[query_heuristic][subplot])
                    t_star = t.isf((1 - confidence_level) / 2, num_samples - 1)
                    factor = t_star / np.sqrt(num_samples)
                    exact_std[i] = data_key[query_exact][subplot].std() * factor
                    heuristic_std[i] = data_key[query_heuristic][subplot].std() * factor

                style = '{0}--'
                error_style = '{0}-.'
                plt.clf()
                matplotlib.rcParams['lines.linewidth'] = 1.0
                if subplot == 'cpu_user':
                    plt.plot(x, exact_mean, style.format(plot['exact_color']))
                    (_, caps, _) = plt.errorbar(x, exact_mean, yerr=[exact_std, exact_std], fmt='+',
                                                color=plot['exact_color'], label=plot['exact'], marker='.', capsize=5)
                    #for cap in caps:
                    #    cap.set_markeredgewidth(5)
                    #eb[-1][0].set_linestyle('-.')
                plt.plot(x, heuristic_mean, style.format(plot['heuristic_color']))
                (_, caps, _) = plt.errorbar(x, heuristic_mean, yerr=[heuristic_std, heuristic_std], fmt='+',
                                            color=plot['heuristic_color'], label='heuristic', marker='.', capsize=5)
                #for cap in caps:
                #    cap.set_markeredgewidth(5)
                #eb[-1][0].set_linestyle('-.')
                plt.xlabel(key)
                plt.ylabel(subplot_descriptions[subplot])
                plt.legend()
                plt.savefig(filename)

if __name__ == "__main__":
    data = pd.read_pickle("CMDP_results.pickle")
    postprocess(data, values)