<gh_stars>0

from sklearn.feature_extraction.text import CountVectorizer
from sklearn.model_selection import ParameterGrid
from joblib import Parallel, delayed
from functools import partial
from scipy.stats import entropy
from scipy.stats import norm as Normal
import numpy as np
import pandas as pd
import math
import time
from functools import partial
import warnings
from pdb import set_trace as st

# ------------------------------------------------------------------------------
# Things to plot in global mode (when paralellizing the steps, once the
# hyperparameters were optimized)
max_probs = {}
min_probs = {}
max_itms = {}
min_itms = {}
# ------------------------------------------------------------------------------

class SetHashedDict:

    def __init__(self):
        self.x_keys = []
        self.y_keys = []
        self._values = {}
        self.coo_keys = []
        self.coo_mem = {}


    def __getitem__(self, ks):
        if isinstance(ks, list) or isinstance(ks, tuple):
            ka, kb = ks
            ix = self.x_keys.index(ka)
            iy = self.y_keys.index(kb)

            return self._values[(ix, iy)]
        else:
            ix = self.x_keys.index(ks)
            return self._values[ix]


    def __setitem__(self, ks, value):
        if isinstance(ks, list) or isinstance(ks, tuple):
            ka, kb = ks
            if not ka in self.x_keys:
                self.x_keys.append(ka)

            if not kb in self.y_keys:
                self.y_keys.append(kb)

            self._values[(
                self.x_keys.index(ka),
                self.y_keys.index(kb)
            )] = value
        else:
            if not ks in self.x_keys:
                self.x_keys.append(ks)

            self._values[self.x_keys.index(ks)] = value


    def __contains__(self, ks):
        if isinstance(ks, list) or isinstance(ks, tuple):
            ka, kb = ks
            return ka in self.x_keys and kb in self.y_keys
        else:
            return ks in self.x_keys


    def keys(self):
        if len(self.y_keys) > 0:
             for x, y in zip(self.x_keys, self.y_keys):
                yield x, y
        else:
            for k in self.x_keys:
                yield k

    def items(self):
        if len(self.y_keys) > 0:
            for it in self._values.items():
                yield (
                    self.x_keys[it[0][0]],
                    self.y_keys[it[0][1]], it[1])
        else:
            for it in self._values.items():
                yield (self.x_keys[it[0]], it[1])


    def values(self):
        return self._values.values()


    def set_coo_mem(self, ks):
        assert isinstance(ks, list) or isinstance(ks, tuple)
        k, value = ks
        if not k in self.coo_keys:
            self.coo_keys.append(k)

        self.coo_mem[self.coo_keys.index(k)] = value


    def get_coo_mem(self, ks):
        assert isinstance(ks, set)
        idx = self.coo_keys.index(ks)

        return self.coo_mem[idx]


class RandomSetDistributions(object):

    def __init__(self,
        gamma=None, bias=4, joint_method='conditional', kernel=None,
        n_jobs=-1, backend='processes', ngram_range=(1, 4), analyzer='char_wb'):
        self.gamma = gamma
        self.bias = bias
        self.analyzer = analyzer
        self.ngram_range = ngram_range
        self.projs = {}
        self.n_jobs = n_jobs
        self.backend = backend
        self.kernel = kernel
        self.set_rvs = {}
        self.distributions = {}
        self.it_metrics = {}
        self.Omega = {}
        self.check_np = lambda f : type(f).__module__ == np.__name__
        self.tokenizer = self.build_tokenizer()


    def build_tokenizer(self):
        tokenizer = CountVectorizer(
            analyzer=self.analyzer, ngram_range=self.ngram_range)
        return tokenizer.build_analyzer()


    def _build_vocab(self, rv):
        omega = []
        for x in self.set_rvs[rv]:
            if x in omega:
                continue
            else:
                omega.append(x)

        return omega


    def _build_set_outcomes(self, rv):
        sets = []
        for x in self.df[rv]:
            sets.append(set(self.tokenizer(str(x))))

        return sets


    def _kernel(self, x):
        # Pending to add more kernels from Hyndman and Scott
        try:
            if self.kernel is None:
                return np.array(x)
            elif self.kernel == 'gaussian':
                #return (1/np.sqrt(2 * math.pi)) * np.exp(-self.gamma * np.array(x) ** 2)
                return np.exp(-self.gamma * np.array(x) ** 2)
            elif self.kernel == 'gausset':
                return np.exp(-self.gamma * (np.array(x) + self.bias) ** 2)
            elif self.kernel == 'expset':
                return np.exp(self.gamma * (np.array(x) + self.bias))
            elif (isinstance(self.kernel, collections.abc.Callable)
                    and self.check_np(self.kernel)):
                return self.kernel(np.array(x))
            else:
                print("WARNING: Bad kernel function specified."
                        " Only intersect returned")
                return np.array(x)
        except OverflowError:
            print("OverflowError in kernel"
               " computations: f({}; {}) undetermined".format(x, self.gamma))
            return None


    def _projector(self, x, y):
        if self.kernel == 'gausset':
            # It would be great to include normalized distances
            # \cite{conci2018distance,app11041910}
            return len((x - y).union(y - x))
        elif self.kernel == 'expset':
            # # It would be great to include similarity metrics
            return len(x.intersection(y))
        else:
            return len(x.intersection(y))


    def _build_projections(self, rv):
        rv_projections = SetHashedDict()
        for x in self.Omega[rv]:
            projection = []
            for x_ in self.set_rvs[rv]:
                projection.append(self._projector(x, x_))

            rv_projections[x] = projection

        return rv_projections


    def _marginal(self, rv='X'):
        assert isinstance(rv, str)

        if rv in list(self.set_rvs.keys()):
            self.Omega[rv] = self._build_vocab(rv)
        else:
            self.set_rvs[rv] = self._build_set_outcomes(rv)
            self.Omega[rv] = self._build_vocab(rv)

        self.projs[rv] = self._build_projections(rv)

        mem = SetHashedDict()
        partition = 0
        for x, proj in self.projs[rv].items():
            f_X = self._kernel(proj).sum()
            mem[x] = f_X
            partition += f_X

        distribution = SetHashedDict()
        for x in self.Omega[rv]:
            lkhood = mem[x]
            if not 0.0 in [lkhood, partition]:
                distribution[x] = mem[x] / partition
            else:
                distribution[x] = 0.0

        self.distributions['P_' + rv] = distribution


    def _conditional(self, rvs=['X', 'Y']):
        assert (isinstance(rvs, list) or
                isinstance(rvs, tuple) and
                len(rvs) == 2), "Check 'rvs' argument for two RVs"
        try:
            omega_x = self.Omega[rvs[0]]
        except KeyError:
            self.Omega[rvs[0]] = self._build_vocab(rvs[0])
            omega_x = self.Omega[rvs[0]]
        try:
            omega_y = self.Omega[rvs[1]]
        except KeyError:
            self.Omega[rvs[1]] = self._build_vocab(rvs[1])
            omega_y = self.Omega[rvs[1]]
        try:
            projs_x = self.projs[rvs[0]]
        except KeyError:
            projs_x = self._build_projections(rvs[0])
        try:
            projs_y = self.projs[rvs[1]]
        except KeyError:
            projs_y = self._build_projections(rvs[1])
        # (1996) Estimating and Visualizing Conditional Densities [Hyndman - Bashtannyk - Grunwald]
        # (2005) From error probability to information theoretic (multi-modal) signal processing%0a[Butz - Thiran]
        # Estimating joint density using kernel product estimator (joint histogram, Butz)
        mem = SetHashedDict()
        for y, proj_y in projs_y.items():
            partition = 0
            for x, proj_x in projs_x.items():
                f_xy = self._kernel(proj_x).dot(self._kernel(proj_y))
                mem[x, y] = f_xy
                partition += f_xy

            mem.set_coo_mem((y, partition))

        distribution = SetHashedDict()
        for x in omega_x:
            for y in omega_y:
                lkhood = mem[x, y]
                partition = mem.get_coo_mem(y)
                if not 0.0 in [lkhood, partition]:
                    distribution[x, y] = mem[x, y] / mem.get_coo_mem(y)
                else:
                    distribution[x, y] = 0.0

        self.distributions['P_' + '|'.join(rvs)] = distribution


    def _entropy(self, rv):
        P_X = self.distributions['P_' + rv]
        P_X = np.array(list(P_X.values()))
        P_X = P_X[P_X != 0]
        if len(P_X) == 0:
            self.it_metrics['H(' + rv + ')'] = 0.0
            return 0.0

        H_X = entropy(P_X, base=2)
        self.it_metrics['H(' + rv + ')'] = H_X

        return H_X


    def _cond_entropy(self, rvs=['X', 'Y']):
        P_XgY = self.distributions['P_' + '|'.join(rvs)]
        P_Y = self.distributions['P_' + rvs[1]]
        omega_x = self.Omega[rvs[0]]
        omega_y = self.Omega[rvs[1]]
        entropi = lambda P: -sum([px * np.log2(px) for px in P if px > 0])
        H_Xgy = [
            P_Y[y] * entropi(
                [P_XgY[x, y] for x in omega_x])#, base=2)
                    for y in omega_y]
        H_XgY = sum(H_Xgy)
        self.it_metrics['H(' + '|'.join(rvs) + ')'] = H_XgY

        return H_XgY


    def _mutual_info(self, rvs=['X', 'Y']):
        try:
            H_X = self.it_metrics['H(' + rvs[0] + ')']
        except KeyError:
            H_X = self._entropy(rvs[0])
        try:
            self.it_metrics['H(' + '|'.join(rvs) + ')']
        except KeyError:
            H_XgY = self._cond_entropy(rvs)

        I_XY = H_X - H_XgY
        self.it_metrics['I(' + ';'.join(rvs) + ')'] = I_XY

        return I_XY


    def fit(self, df, it_rvs=['X,Y', 'Y,Z', 'Z,X']):
        # Verify wheter time can be shorten takin the idea of recursive computation
        # from Lodhi, Huma, et al. "Text classification using string kernels."
        #           Journal of Machine Learning Research 2.Feb (2002): 419-444.
        assert isinstance(df, pd.DataFrame) # Input data must a pd.DataFrame
        assert len(df.columns) >= 2 # Two or more columns are needed

        self.it_rvs = it_rvs
        self.df = df
        pairsRV = [rvs.split(',') for rvs in self.it_rvs]

        # Estimate Marginal distributions
        RVs = list(set(sum(pairsRV, [])))

        assert len(RVs)/len(set(df.columns).intersection(RVs)) == 1.0, ("One or"
            " more of the given columns in 'it_rvs' argument to fit() is not"
            "in the columns of the input dataframe.")

        for RV in RVs:
            self._marginal(RV)
        # Estimating pair-wise Information-Theoretic metrics

        for pair in pairsRV:
            self._conditional(pair)
            self._mutual_info(pair)
            # Change to this method if no problem when returnning the metrics:
            #key = 'I(' + ','.join(pair) + ')'
            #self.it_metrics[key] = self._mutual_info(pair)


def semantic_reward(df, cols, beta=1e3):
    line = {}
    mean_df = df[cols].mean().sort_values()
    a, b, c = mean_df.values
    sa, sb, sc = df[mean_df.index].std().values
    a_c = abs(a - c)
    #l = a_c / 2
    #dist = abs(b - l)
    # https://en.wikipedia.org/wiki/List_of_triangle_inequalities
    harmonic = 3 * a * b * c / (a * b + b * c + c * a)
    geometric = (a * b * c) ** (1. / 3.)
    arithmetic = np.mean([a, b, c])
    equilateral = (1. - ((geometric - harmonic) / a_c)) + (1. - ((arithmetic - geometric) / a_c))
    line["Reward"] = equilateral
    #try:
    #    z = abs(a - c) / (sb * math.sqrt(2 * math.pi))
    #    line["Reward"] = z * math.exp(-dist ** 2/(2 * sb ** 2))
    #except:
    #    line["Reward"] = 0.0

    pvals = (
                Normal(loc=a, scale=sa).pdf(b - 2 * sb),
                Normal(loc=b, scale=sb).pdf(c - 2 * sc),
                Normal(loc=c, scale=sc).pdf(a + 2 * sa)  # verificar (b + 2 * sb), etc. es decir en todas direcciones
            )
    line.update(
                zip(['ABpvalue',
                     'BCpvalue',
                     'CApvalue'], pvals))
    line["pReward"] = math.exp(-beta * max(pvals))
    print(line["pReward"])
    return line


def save_extrema(set_statistics):

    global max_probs
    global min_probs
    global max_itms
    global min_itms
    if max_probs == {}:
        max_probs = {k: [0, ''] for k in set_statistics.distributions.keys()}
    if min_probs == {}:
        min_probs = {k: [1.5, ''] for k in set_statistics.distributions.keys()}
    if max_itms == {}:
        max_itms = {k: [0, ''] for k in set_statistics.it_metrics.keys()}
    if min_itms == {}:
        min_itms = {k: [1000, ''] for k in set_statistics.it_metrics.keys()}

    for k, vs in set_statistics.distributions.items():
        mean_prob = np.mean(list(vs.values()))
        if mean_prob > max_probs[k][0]:
            max_probs[k] = (mean_prob, set_statistics)
        if mean_prob < min_probs[k][0]:
            min_probs[k] = (mean_prob, set_statistics)

    for k, vs in set_statistics.it_metrics.items():
        if vs > max_itms[k][0]:
            max_itms[k] = (vs, set_statistics)
        if vs < min_itms[k][0]:
            min_itms[k] = (vs, set_statistics)


def step_processing(
    df, it_rvs=['Y,X', 'Z,Y', 'X,Z'], kernel='gausset', gamma=1.0/50.0, bias=4.0):

    #start_time = time.time()

    df = df.dropna()
    random_sets = RandomSetDistributions(kernel=kernel, gamma=gamma, bias=bias)
    random_sets.fit(df, it_rvs=it_rvs)
    # --------------------------------------------------------------------------
#### Activate for plotting
    #save_extrema(random_sets)
    # --------------------------------------------------------------------------
    #end_time = time.time()
    #print("STEP time: {}".format(end_time-start_time))

    return random_sets.it_metrics


def episode_processing(chunk_size, kernel, gamma,bias, input_triplets, n_jobs=1):
    start_time = time.time()
    it_rvs = ['Y,X', 'Z,Y', 'X,Z']  # Order verified from De Marcken (1999)
    output_it = "results/it_{}_kernel-{}_gamma-{}_bias-{}_sample-{}.csv".format(
        'train' if '_train.' in input_triplets else 'test',
        kernel,
        gamma,
        bias,
        chunk_size)
    columns = [1, 2, 3]
    names = ['X', 'Y', 'Z']

    # Generate chunks from input csv triplets
    df_generator = pd.read_csv(
        input_triplets,
        sep='\t',
        names=names,
        chunksize=chunk_size,
        usecols=columns)

    step = partial(
        step_processing, it_rvs=it_rvs, kernel=kernel, gamma=gamma, bias=bias)
    result = Parallel(
        n_jobs=n_jobs,
####### Activate for plotting
	#require='sharedmem'
        prefer='processes'
        )(
            delayed(step)(df)
                for df in df_generator)

    end_time = time.time()
    print("Episode time: {}".format(end_time-start_time))
    print(output_it)
    return output_it, pd.DataFrame(result)
#    pd.DataFrame(results).to_csv(output_it)


def parse_params(strg):
    pairs = strg.split('/')[-1].replace('.csv', '').split('_')[2:]
    lpairs = dict([tuple(p.split('-')) if p.split('-')[0] == 'kernel'
                           else (p.split('-')[0], float(p.split('-')[1])) for p in pairs])
    return lpairs


if __name__ == "__main__":

    input_triplets = 'data/dis_train.txt.oie'
    beta = 1.0
    hyperpar = {
        'kernel': ['expset', 'gausset'],
        'gamma': [1., 1./10, 1./50., 1./100., 1./500., 1./1000.],
        'bias': [1, 5, 10, 15, 20, 25, 50],
        'chunk_size': [10, 50, 100, 200, 300],
    }


    grid = ParameterGrid(hyperpar)
    episode = partial(
        episode_processing, input_triplets=input_triplets)
    #param = [p for p in grid if p['kernel'] == 'gausset' and p['gamma'] == 1][-1]
    #episode(**param)

    results = Parallel(
        n_jobs=-1,
        prefer='processes'
        )(
            delayed(episode)(**params)
        for params in grid)
    its = ['I(Y;X)', 'I(Z;Y)', 'I(X;Z)']

    rewards = [{**semantic_reward(rdf, its, beta=beta), **parse_params(params)}
        for params, rdf in results]

    pd.DataFrame(rewards).to_csv('rewards.csv')

from sklearn.feature_extraction.text import CountVectorizer
from sklearn.model_selection import ParameterGrid
from joblib import Parallel, delayed
from functools import partial
from scipy.stats import entropy
from scipy.stats import norm as Normal
import numpy as np
import pandas as pd
import math
import time
from functools import partial
import warnings
from pdb import set_trace as st

# ------------------------------------------------------------------------------
# Things to plot in global mode (when paralellizing the steps, once the
# hyperparameters were optimized)
max_probs = {}
min_probs = {}
max_itms = {}
min_itms = {}
# ------------------------------------------------------------------------------

class SetHashedDict:

    def __init__(self):
        self.x_keys = []
        self.y_keys = []
        self._values = {}
        self.coo_keys = []
        self.coo_mem = {}


    def __getitem__(self, ks):
        if isinstance(ks, list) or isinstance(ks, tuple):
            ka, kb = ks
            ix = self.x_keys.index(ka)
            iy = self.y_keys.index(kb)

            return self._values[(ix, iy)]
        else:
            ix = self.x_keys.index(ks)
            return self._values[ix]


    def __setitem__(self, ks, value):
        if isinstance(ks, list) or isinstance(ks, tuple):
            ka, kb = ks
            if not ka in self.x_keys:
                self.x_keys.append(ka)

            if not kb in self.y_keys:
                self.y_keys.append(kb)

            self._values[(
                self.x_keys.index(ka),
                self.y_keys.index(kb)
            )] = value
        else:
            if not ks in self.x_keys:
                self.x_keys.append(ks)

            self._values[self.x_keys.index(ks)] = value


    def __contains__(self, ks):
        if isinstance(ks, list) or isinstance(ks, tuple):
            ka, kb = ks
            return ka in self.x_keys and kb in self.y_keys
        else:
            return ks in self.x_keys


    def keys(self):
        if len(self.y_keys) > 0:
             for x, y in zip(self.x_keys, self.y_keys):
                yield x, y
        else:
            for k in self.x_keys:
                yield k

    def items(self):
        if len(self.y_keys) > 0:
            for it in self._values.items():
                yield (
                    self.x_keys[it[0][0]],
                    self.y_keys[it[0][1]], it[1])
        else:
            for it in self._values.items():
                yield (self.x_keys[it[0]], it[1])


    def values(self):
        return self._values.values()


    def set_coo_mem(self, ks):
        assert isinstance(ks, list) or isinstance(ks, tuple)
        k, value = ks
        if not k in self.coo_keys:
            self.coo_keys.append(k)

        self.coo_mem[self.coo_keys.index(k)] = value


    def get_coo_mem(self, ks):
        assert isinstance(ks, set)
        idx = self.coo_keys.index(ks)

        return self.coo_mem[idx]


class RandomSetDistributions(object):

    def __init__(self,
        gamma=None, bias=4, joint_method='conditional', kernel=None,
        n_jobs=-1, backend='processes', ngram_range=(1, 4), analyzer='char_wb'):
        self.gamma = gamma
        self.bias = bias
        self.analyzer = analyzer
        self.ngram_range = ngram_range
        self.projs = {}
        self.n_jobs = n_jobs
        self.backend = backend
        self.kernel = kernel
        self.set_rvs = {}
        self.distributions = {}
        self.it_metrics = {}
        self.Omega = {}
        self.check_np = lambda f : type(f).__module__ == np.__name__
        self.tokenizer = self.build_tokenizer()


    def build_tokenizer(self):
        tokenizer = CountVectorizer(
            analyzer=self.analyzer, ngram_range=self.ngram_range)
        return tokenizer.build_analyzer()


    def _build_vocab(self, rv):
        omega = []
        for x in self.set_rvs[rv]:
            if x in omega:
                continue
            else:
                omega.append(x)

        return omega


    def _build_set_outcomes(self, rv):
        sets = []
        for x in self.df[rv]:
            sets.append(set(self.tokenizer(str(x))))

        return sets


    def _kernel(self, x):
        # Pending to add more kernels from Hyndman and Scott
        try:
            if self.kernel is None:
                return np.array(x)
            elif self.kernel == 'gaussian':
                #return (1/np.sqrt(2 * math.pi)) * np.exp(-self.gamma * np.array(x) ** 2)
                return np.exp(-self.gamma * np.array(x) ** 2)
            elif self.kernel == 'gausset':
                return np.exp(-self.gamma * (np.array(x) + self.bias) ** 2)
            elif self.kernel == 'expset':
                return np.exp(self.gamma * (np.array(x) + self.bias))
            elif (isinstance(self.kernel, collections.abc.Callable)
                    and self.check_np(self.kernel)):
                return self.kernel(np.array(x))
            else:
                print("WARNING: Bad kernel function specified."
                        " Only intersect returned")
                return np.array(x)
        except OverflowError:
            print("OverflowError in kernel"
               " computations: f({}; {}) undetermined".format(x, self.gamma))
            return None


    def _projector(self, x, y):
        if self.kernel == 'gausset':
            # It would be great to include normalized distances
            # \cite{conci2018distance,app11041910}
            return len((x - y).union(y - x))
        elif self.kernel == 'expset':
            # # It would be great to include similarity metrics
            return len(x.intersection(y))
        else:
            return len(x.intersection(y))


    def _build_projections(self, rv):
        rv_projections = SetHashedDict()
        for x in self.Omega[rv]:
            projection = []
            for x_ in self.set_rvs[rv]:
                projection.append(self._projector(x, x_))

            rv_projections[x] = projection

        return rv_projections


    def _marginal(self, rv='X'):
        assert isinstance(rv, str)

        if rv in list(self.set_rvs.keys()):
            self.Omega[rv] = self._build_vocab(rv)
        else:
            self.set_rvs[rv] = self._build_set_outcomes(rv)
            self.Omega[rv] = self._build_vocab(rv)

        self.projs[rv] = self._build_projections(rv)

        mem = SetHashedDict()
        partition = 0
        for x, proj in self.projs[rv].items():
            f_X = self._kernel(proj).sum()
            mem[x] = f_X
            partition += f_X

        distribution = SetHashedDict()
        for x in self.Omega[rv]:
            lkhood = mem[x]
            if not 0.0 in [lkhood, partition]:
                distribution[x] = mem[x] / partition
            else:
                distribution[x] = 0.0

        self.distributions['P_' + rv] = distribution


    def _conditional(self, rvs=['X', 'Y']):
        assert (isinstance(rvs, list) or
                isinstance(rvs, tuple) and
                len(rvs) == 2), "Check 'rvs' argument for two RVs"
        try:
            omega_x = self.Omega[rvs[0]]
        except KeyError:
            self.Omega[rvs[0]] = self._build_vocab(rvs[0])
            omega_x = self.Omega[rvs[0]]
        try:
            omega_y = self.Omega[rvs[1]]
        except KeyError:
            self.Omega[rvs[1]] = self._build_vocab(rvs[1])
            omega_y = self.Omega[rvs[1]]
        try:
            projs_x = self.projs[rvs[0]]
        except KeyError:
            projs_x = self._build_projections(rvs[0])
        try:
            projs_y = self.projs[rvs[1]]
        except KeyError:
            projs_y = self._build_projections(rvs[1])
        # (1996) Estimating and Visualizing Conditional Densities [Hyndman - Bashtannyk - Grunwald]
        # (2005) From error probability to information theoretic (multi-modal) signal processing%0a[Butz - Thiran]
        # Estimating joint density using kernel product estimator (joint histogram, Butz)
        mem = SetHashedDict()
        for y, proj_y in projs_y.items():
            partition = 0
            for x, proj_x in projs_x.items():
                f_xy = self._kernel(proj_x).dot(self._kernel(proj_y))
                mem[x, y] = f_xy
                partition += f_xy

            mem.set_coo_mem((y, partition))

        distribution = SetHashedDict()
        for x in omega_x:
            for y in omega_y:
                lkhood = mem[x, y]
                partition = mem.get_coo_mem(y)
                if not 0.0 in [lkhood, partition]:
                    distribution[x, y] = mem[x, y] / mem.get_coo_mem(y)
                else:
                    distribution[x, y] = 0.0

        self.distributions['P_' + '|'.join(rvs)] = distribution


    def _entropy(self, rv):
        P_X = self.distributions['P_' + rv]
        P_X = np.array(list(P_X.values()))
        P_X = P_X[P_X != 0]
        if len(P_X) == 0:
            self.it_metrics['H(' + rv + ')'] = 0.0
            return 0.0

        H_X = entropy(P_X, base=2)
        self.it_metrics['H(' + rv + ')'] = H_X

        return H_X


    def _cond_entropy(self, rvs=['X', 'Y']):
        P_XgY = self.distributions['P_' + '|'.join(rvs)]
        P_Y = self.distributions['P_' + rvs[1]]
        omega_x = self.Omega[rvs[0]]
        omega_y = self.Omega[rvs[1]]
        entropi = lambda P: -sum([px * np.log2(px) for px in P if px > 0])
        H_Xgy = [
            P_Y[y] * entropi(
                [P_XgY[x, y] for x in omega_x])#, base=2)
                    for y in omega_y]
        H_XgY = sum(H_Xgy)
        self.it_metrics['H(' + '|'.join(rvs) + ')'] = H_XgY

        return H_XgY


    def _mutual_info(self, rvs=['X', 'Y']):
        try:
            H_X = self.it_metrics['H(' + rvs[0] + ')']
        except KeyError:
            H_X = self._entropy(rvs[0])
        try:
            self.it_metrics['H(' + '|'.join(rvs) + ')']
        except KeyError:
            H_XgY = self._cond_entropy(rvs)

        I_XY = H_X - H_XgY
        self.it_metrics['I(' + ';'.join(rvs) + ')'] = I_XY

        return I_XY


    def fit(self, df, it_rvs=['X,Y', 'Y,Z', 'Z,X']):
        # Verify wheter time can be shorten takin the idea of recursive computation
        # from Lodhi, Huma, et al. "Text classification using string kernels."
        #           Journal of Machine Learning Research 2.Feb (2002): 419-444.
        assert isinstance(df, pd.DataFrame) # Input data must a pd.DataFrame
        assert len(df.columns) >= 2 # Two or more columns are needed

        self.it_rvs = it_rvs
        self.df = df
        pairsRV = [rvs.split(',') for rvs in self.it_rvs]

        # Estimate Marginal distributions
        RVs = list(set(sum(pairsRV, [])))

        assert len(RVs)/len(set(df.columns).intersection(RVs)) == 1.0, ("One or"
            " more of the given columns in 'it_rvs' argument to fit() is not"
            "in the columns of the input dataframe.")

        for RV in RVs:
            self._marginal(RV)
        # Estimating pair-wise Information-Theoretic metrics

        for pair in pairsRV:
            self._conditional(pair)
            self._mutual_info(pair)
            # Change to this method if no problem when returnning the metrics:
            #key = 'I(' + ','.join(pair) + ')'
            #self.it_metrics[key] = self._mutual_info(pair)


def semantic_reward(df, cols, beta=1e3):
    line = {}
    mean_df = df[cols].mean().sort_values()
    a, b, c = mean_df.values
    sa, sb, sc = df[mean_df.index].std().values
    a_c = abs(a - c)
    #l = a_c / 2
    #dist = abs(b - l)
    # https://en.wikipedia.org/wiki/List_of_triangle_inequalities
    harmonic = 3 * a * b * c / (a * b + b * c + c * a)
    geometric = (a * b * c) ** (1. / 3.)
    arithmetic = np.mean([a, b, c])
    equilateral = (1. - ((geometric - harmonic) / a_c)) + (1. - ((arithmetic - geometric) / a_c))
    line["Reward"] = equilateral
    #try:
    #    z = abs(a - c) / (sb * math.sqrt(2 * math.pi))
    #    line["Reward"] = z * math.exp(-dist ** 2/(2 * sb ** 2))
    #except:
    #    line["Reward"] = 0.0

    pvals = (
                Normal(loc=a, scale=sa).pdf(b - 2 * sb),
                Normal(loc=b, scale=sb).pdf(c - 2 * sc),
                Normal(loc=c, scale=sc).pdf(a + 2 * sa)  # verificar (b + 2 * sb), etc. es decir en todas direcciones
            )
    line.update(
                zip(['ABpvalue',
                     'BCpvalue',
                     'CApvalue'], pvals))
    line["pReward"] = math.exp(-beta * max(pvals))
    print(line["pReward"])
    return line


def save_extrema(set_statistics):

    global max_probs
    global min_probs
    global max_itms
    global min_itms
    if max_probs == {}:
        max_probs = {k: [0, ''] for k in set_statistics.distributions.keys()}
    if min_probs == {}:
        min_probs = {k: [1.5, ''] for k in set_statistics.distributions.keys()}
    if max_itms == {}:
        max_itms = {k: [0, ''] for k in set_statistics.it_metrics.keys()}
    if min_itms == {}:
        min_itms = {k: [1000, ''] for k in set_statistics.it_metrics.keys()}

    for k, vs in set_statistics.distributions.items():
        mean_prob = np.mean(list(vs.values()))
        if mean_prob > max_probs[k][0]:
            max_probs[k] = (mean_prob, set_statistics)
        if mean_prob < min_probs[k][0]:
            min_probs[k] = (mean_prob, set_statistics)

    for k, vs in set_statistics.it_metrics.items():
        if vs > max_itms[k][0]:
            max_itms[k] = (vs, set_statistics)
        if vs < min_itms[k][0]:
            min_itms[k] = (vs, set_statistics)


def step_processing(
    df, it_rvs=['Y,X', 'Z,Y', 'X,Z'], kernel='gausset', gamma=1.0/50.0, bias=4.0):

    #start_time = time.time()

    df = df.dropna()
    random_sets = RandomSetDistributions(kernel=kernel, gamma=gamma, bias=bias)
    random_sets.fit(df, it_rvs=it_rvs)
    # --------------------------------------------------------------------------
#### Activate for plotting
    #save_extrema(random_sets)
    # --------------------------------------------------------------------------
    #end_time = time.time()
    #print("STEP time: {}".format(end_time-start_time))

    return random_sets.it_metrics


def episode_processing(chunk_size, kernel, gamma,bias, input_triplets, n_jobs=1):
    start_time = time.time()
    it_rvs = ['Y,X', 'Z,Y', 'X,Z']  # Order verified from De Marcken (1999)
    output_it = "results/it_{}_kernel-{}_gamma-{}_bias-{}_sample-{}.csv".format(
        'train' if '_train.' in input_triplets else 'test',
        kernel,
        gamma,
        bias,
        chunk_size)
    columns = [1, 2, 3]
    names = ['X', 'Y', 'Z']

    # Generate chunks from input csv triplets
    df_generator = pd.read_csv(
        input_triplets,
        sep='\t',
        names=names,
        chunksize=chunk_size,
        usecols=columns)

    step = partial(
        step_processing, it_rvs=it_rvs, kernel=kernel, gamma=gamma, bias=bias)
    result = Parallel(
        n_jobs=n_jobs,
####### Activate for plotting
	#require='sharedmem'
        prefer='processes'
        )(
            delayed(step)(df)
                for df in df_generator)

    end_time = time.time()
    print("Episode time: {}".format(end_time-start_time))
    print(output_it)
    return output_it, pd.DataFrame(result)
#    pd.DataFrame(results).to_csv(output_it)


def parse_params(strg):
    pairs = strg.split('/')[-1].replace('.csv', '').split('_')[2:]
    lpairs = dict([tuple(p.split('-')) if p.split('-')[0] == 'kernel'
                           else (p.split('-')[0], float(p.split('-')[1])) for p in pairs])
    return lpairs


if __name__ == "__main__":

    input_triplets = 'data/dis_train.txt.oie'
    beta = 1.0
    hyperpar = {
        'kernel': ['expset', 'gausset'],
        'gamma': [1., 1./10, 1./50., 1./100., 1./500., 1./1000.],
        'bias': [1, 5, 10, 15, 20, 25, 50],
        'chunk_size': [10, 50, 100, 200, 300],
    }


    grid = ParameterGrid(hyperpar)
    episode = partial(
        episode_processing, input_triplets=input_triplets)
    #param = [p for p in grid if p['kernel'] == 'gausset' and p['gamma'] == 1][-1]
    #episode(**param)

    results = Parallel(
        n_jobs=-1,
        prefer='processes'
        )(
            delayed(episode)(**params)
        for params in grid)
    its = ['I(Y;X)', 'I(Z;Y)', 'I(X;Z)']

    rewards = [{**semantic_reward(rdf, its, beta=beta), **parse_params(params)}
        for params, rdf in results]

    pd.DataFrame(rewards).to_csv('rewards.csv')
# ------------------------------------------------------------------------------ # Things to plot in global mode (when paralellizing the steps, once the # hyperparameters were optimized) # ------------------------------------------------------------------------------ # Pending to add more kernels from Hyndman and Scott #return (1/np.sqrt(2 * math.pi)) * np.exp(-self.gamma * np.array(x) ** 2) # It would be great to include normalized distances # \cite{conci2018distance,app11041910} # # It would be great to include similarity metrics # (1996) Estimating and Visualizing Conditional Densities [Hyndman - Bashtannyk - Grunwald] # (2005) From error probability to information theoretic (multi-modal) signal processing%0a[Butz - Thiran] # Estimating joint density using kernel product estimator (joint histogram, Butz) #, base=2) # Verify wheter time can be shorten takin the idea of recursive computation # from Lodhi, Huma, et al. "Text classification using string kernels." #           Journal of Machine Learning Research 2.Feb (2002): 419-444. # Input data must a pd.DataFrame # Two or more columns are needed # Estimate Marginal distributions # Estimating pair-wise Information-Theoretic metrics # Change to this method if no problem when returnning the metrics: #key = 'I(' + ','.join(pair) + ')' #self.it_metrics[key] = self._mutual_info(pair) #l = a_c / 2 #dist = abs(b - l) # https://en.wikipedia.org/wiki/List_of_triangle_inequalities #try: #    z = abs(a - c) / (sb * math.sqrt(2 * math.pi)) #    line["Reward"] = z * math.exp(-dist ** 2/(2 * sb ** 2)) #except: #    line["Reward"] = 0.0 # verificar (b + 2 * sb), etc. es decir en todas direcciones #start_time = time.time() # -------------------------------------------------------------------------- #### Activate for plotting #save_extrema(random_sets) # -------------------------------------------------------------------------- #end_time = time.time() #print("STEP time: {}".format(end_time-start_time)) # Order verified from De Marcken (1999) # Generate chunks from input csv triplets ####### Activate for plotting #require='sharedmem' #    pd.DataFrame(results).to_csv(output_it) #param = [p for p in grid if p['kernel'] == 'gausset' and p['gamma'] == 1][-1] #episode(**param)