text:
import warnings
import numpy as np
import theano
import theano.tensor as tt

import pymc3 as pm
from .updates import adam
from ..distributions.dist_math import rho2sd, log_normal
from ..model import modelcontext, ArrayOrdering
from ..theanof import tt_rng, memoize, change_flags, GradScale

__all__ = [
    'ObjectiveFunction',
    'Operator',
    'TestFunction',
    'Approximation'
]

class ObjectiveUpdates(theano.OrderedUpdates):
    """
    OrderedUpdates extension for storing loss
    """
    loss = None

def _warn_not_used(smth, where):
    warnings.warn('`%s` is not used for %s and ignored' % (smth, where))

class ObjectiveFunction(object):
    """
    Helper class for construction loss and updates for variational inference

    Parameters
    ----------
    op : Operator
    tf : TestFunction
    """
    def __init__(self, op, tf):
        self.op = op
        self.tf = tf

    obj_params = property(lambda self: self.op.approx.params)
    test_params = property(lambda self: self.tf.params)
    approx = property(lambda self: self.op.approx)

    def random(self, size=None):
        """
        Posterior distribution from initial latent space

        Parameters
        ----------
        size : int
            number of samples from distribution

        Returns
        -------
        posterior space (theano)
        """
        return self.op.approx.random(size)

    def updates(self, obj_n_mc=None, tf_n_mc=None, obj_optimizer=adam, test_optimizer=adam,
                more_obj_params=None, more_tf_params=None, more_updates=None, more_replacements=None):
        """
        Calculates gradients for objective function, test function and then
        constructs updates for optimization step

        Parameters
        ----------
        obj_n_mc : int
            Number of monte carlo samples used for approximation of objective gradients
        tf_n_mc : int
            Number of monte carlo samples used for approximation of test function gradients
        obj_optimizer : function (loss, params) -> updates
            Optimizer that is used for objective params
        test_optimizer : function (loss, params) -> updates
            Optimizer that is used for test function params
        more_obj_params : list
            Add custom params for objective optimizer
        more_tf_params : list
            Add custom params for test function optimizer
        more_updates : dict
            Add custom updates to resulting updates
        more_replacements : dict
            Apply custom replacements before calculating gradients

        Returns
        -------
        ObjectiveUpdates
        """
        if more_obj_params is None:
            more_obj_params = []
        if more_tf_params is None:
            more_tf_params = []
        if more_updates is None:
            more_updates = dict()
        if more_replacements is None:
            more_replacements = dict()
        if not self.op.RETURNS_LOSS and (more_obj_params or more_tf_params):
            raise ValueError('%s does not support, differentiation with '
                             'additional params, try passing your updates '
                             'via more_updates kwarg' % self.op.__class__)
        resulting_updates = ObjectiveUpdates()
        if self.test_params:
            self.add_test_updates(
                resulting_updates,
                tf_n_mc=tf_n_mc,
                test_optimizer=test_optimizer,
                more_tf_params=more_tf_params,
                more_replacements=more_replacements
            )
        else:
            if tf_n_mc is not None:
                _warn_not_used('tf_n_mc', self.op)
            if more_tf_params:
                _warn_not_used('more_tf_params', self.op)
        self.add_obj_updates(
            resulting_updates,
            obj_n_mc=obj_n_mc,
            obj_optimizer=obj_optimizer,
            more_obj_params=more_obj_params,
            more_replacements=more_replacements
        )
        resulting_updates.update(more_updates)
        return resulting_updates

    def add_test_updates(self, updates, tf_n_mc=None, test_optimizer=adam, more_tf_params=None, more_replacements=None):
        tf_z = self.random(tf_n_mc)
        tf_target = -self(tf_z)
        tf_target = theano.clone(tf_target, more_replacements, strict=False)
        updates.update(test_optimizer(tf_target, self.test_params + more_tf_params))

    def add_obj_updates(self, updates, obj_n_mc=None, obj_optimizer=adam, more_obj_params=None, more_replacements=None):
        obj_z = self.random(obj_n_mc)
        obj_target = self(obj_z)
        obj_target = theano.clone(obj_target, more_replacements, strict=False)
        updates.update(obj_optimizer(obj_target, self.obj_params + more_obj_params))
        if self.op.RETURNS_LOSS:
            updates.loss = obj_target

    @memoize
    def step_function(self, obj_n_mc=None, tf_n_mc=None,
                      obj_optimizer=adam, test_optimizer=adam,
                      more_obj_params=None, more_tf_params=None,
                      more_updates=None, more_replacements=None, score=False,
                      fn_kwargs=None):
        """
        Step function that should be called on each optimization step.

        Generally it solves the following problem:
        .. math::

                \textbf{\lambda^{*}} = \inf_{\lambda} \sup_{\theta} t(\mathbb{E}_{\lambda}[(O^{p,q}f_{\theta})(z)])

        Parameters
        ----------
        obj_n_mc : int
            Number of monte carlo samples used for approximation of objective gradients
        tf_n_mc : int
            Number of monte carlo samples used for approximation of test function gradients
        obj_optimizer : function (loss, params) -> updates
            Optimizer that is used for objective params
        test_optimizer : function (loss, params) -> updates
            Optimizer that is used for test function params
        more_obj_params : list
            Add custom params for objective optimizer
        more_tf_params : list
            Add custom params for test function optimizer
        more_updates : dict
            Add custom updates to resulting updates
        score : bool
            calculate loss on each step? Defaults to False for speed
        fn_kwargs : dict
            Add kwargs to theano.function (e.g. `{'profile': True}`)
        more_replacements : dict
            Apply custom replacements before calculating gradients

        Returns
        -------
        theano.function
        """
        if fn_kwargs is None:
            fn_kwargs = {}
        if score and not self.op.RETURNS_LOSS:
            raise NotImplementedError('%s does not have loss' % self.op)
        updates = self.updates(obj_n_mc=obj_n_mc, tf_n_mc=tf_n_mc,
                               obj_optimizer=obj_optimizer,
                               test_optimizer=test_optimizer,
                               more_obj_params=more_obj_params,
                               more_tf_params=more_tf_params,
                               more_updates=more_updates,
                               more_replacements=more_replacements)
        if score:
            step_fn = theano.function([], updates.loss, updates=updates, **fn_kwargs)
        else:
            step_fn = theano.function([], None, updates=updates, **fn_kwargs)
        return step_fn

    @memoize
    def score_function(self, sc_n_mc=None, more_replacements=None, fn_kwargs=None):   # pragma: no cover
        """
        Compiles scoring function that operates which takes no inputs and returns Loss

        Parameters
        ----------
        sc_n_mc : int
            number of scoring MC samples
        more_replacements:
            Apply custom replacements before compiling a function
        fn_kwargs:
            arbitrary kwargs passed to theano.function

        Returns
        -------
        theano.function
        """
        if fn_kwargs is None:
            fn_kwargs = {}
        if not self.op.RETURNS_LOSS:
            raise NotImplementedError('%s does not have loss' % self.op)
        if more_replacements is None:
            more_replacements = {}
        loss = theano.clone(self(self.random(sc_n_mc)), more_replacements, strict=False)
        return theano.function([], loss, **fn_kwargs)

    def __getstate__(self):
        return self.op, self.tf

    def __setstate__(self, state):
        self.__init__(*state)

    def __call__(self, z):
        if z.ndim > 1:
            a = theano.scan(
                lambda z_: theano.clone(self.op.apply(self.tf), {self.op.input: z_}, strict=False),
                sequences=z, n_steps=z.shape[0])[0].mean()
        else:
            a = theano.clone(self.op.apply(self.tf), {self.op.input: z}, strict=False)
        return tt.abs_(a)

class Operator(object):
    """
    Base class for Operator

    Parameters
    ----------
    approx : Approximation

    Subclassing
    -----------
    For implementing Custom operator it is needed to define `.apply(f)` method
    """

    HAS_TEST_FUNCTION = False
    RETURNS_LOSS = True
    SUPPORT_AEVB = True
    OBJECTIVE = ObjectiveFunction

    def __init__(self, approx):
        if not self.SUPPORT_AEVB and approx.local_vars:
            raise ValueError('%s does not support AEVB, '
                             'please change inference method' % type(self))
        self.model = approx.model
        self.approx = approx

    flat_view = property(lambda self: self.approx.flat_view)
    input = property(lambda self: self.approx.flat_view.input)

    def logp(self, z):
        factors = ([tt.sum(var.logpt)for var in self.model.basic_RVs] +
                   [tt.sum(var) for var in self.model.potentials])

        p = self.approx.to_flat_input(tt.add(*factors))
        p = theano.clone(p, {self.input: z})
        return p

    def logp_norm(self, z):
        t = self.approx.normalizing_constant
        factors = ([tt.sum(var.logpt) / t for var in self.model.basic_RVs] +
                   [tt.sum(var) / t for var in self.model.potentials])
        logpt = tt.add(*factors)
        p = self.approx.to_flat_input(logpt)
        p = theano.clone(p, {self.input: z})
        return p

    def logq(self, z):
        return self.approx.logq(z)

    def logq_norm(self, z):
        return self.approx.logq_norm(z)

    def apply(self, f):   # pragma: no cover
        """
        Operator itself
        .. math::

            (O^{p,q}f_{\theta})(z)

        Parameters
        ----------
        f : function or None
            function that takes `z = self.input` and returns
            same dimension output

        Returns
        -------
        symbolically applied operator
        """
        raise NotImplementedError

    def __call__(self, f=None):
        if self.HAS_TEST_FUNCTION:
            if f is None:
                raise ValueError('Operator %s requires TestFunction' % self)
            else:
                if not isinstance(f, TestFunction):
                    f = TestFunction.from_function(f)
        else:
            if f is not None:
                warnings.warn('TestFunction for %s is redundant and removed' % self)
            else:
                pass
            f = TestFunction()
        f.setup(self.approx.total_size)
        return self.OBJECTIVE(self, f)

    def __getstate__(self):
        # pickle only important parts
        return self.approx

    def __setstate__(self, approx):
        self.__init__(approx)

    def __str__(self):    # pragma: no cover
        return '%(op)s[%(ap)s]' % dict(op=self.__class__.__name__,
                                       ap=self.approx.__class__.__name__)

def cast_to_list(params):
    """
    Helper function for getting a list from
    usable representation of parameters

    Parameters
    ----------
    params : {list|tuple|dict|theano.shared|None}

    Returns
    -------
    list
    """
    if isinstance(params, list):
        return params
    elif isinstance(params, tuple):
        return list(params)
    elif isinstance(params, dict):
        return list(params.values())
    elif isinstance(params, theano.compile.SharedVariable):
        return [params]
    elif params is None:
        return []
    else:
        raise TypeError('Unknown type %s for %r, need list, dict or shared variable')

class TestFunction(object):
    def __init__(self):
        self._inited = False
        self.shared_params = None

    def create_shared_params(self, dim):
        """
        Returns
        -------
        {dict|list|theano.shared}
        """
        pass

    @property
    def params(self):
        return cast_to_list(self.shared_params)

    def __call__(self, z):
        raise NotImplementedError

    def setup(self, dim):
        if not self._inited:
            self._setup(dim)
            self.shared_params = self.create_shared_params(dim)
            self._inited = True

    def _setup(self, dim):
        """
        Does some preparation stuff before calling `.create_shared_params()`

        Parameters
        ----------
        dim : int dimension of posterior distribution
        """
        pass

    @classmethod
    def from_function(cls, f):
        if not callable(f):
            raise ValueError('Need callable, got %r' % f)
        obj = TestFunction()
        obj.__call__ = f
        return obj

class Approximation(object):
    """
    Base class for approximations.

    Parameters
    ----------
    local_rv : dict
        mapping {model_variable -> local_variable}
        Local Vars are used for Autoencoding Variational Bayes
        See (AEVB; Kingma and Welling, 2014) for details

    model : PyMC3 model for inference

    cost_part_grad_scale : float or scalar tensor
        Scaling score part of gradient can be useful near optimum for
        archiving better convergence properties. Common schedule is
        1 at the start and 0 in the end. So slow decay will be ok.
        See (Sticking the Landing; Geoffrey Roeder,
        Yuhuai Wu, David Duvenaud, 2016) for details

    Subclassing
    -----------
    Defining an approximation needs
    custom implementation of the following methods:
        - `.create_shared_params()`
            Returns {dict|list|theano.shared}

        - `.random_global(size=None, no_rand=False)`
            Generate samples from posterior. If `no_rand==False`:
            sample from MAP of initial distribution.
            Returns TensorVariable

        - `.log_q_W_global(z)`
            It is needed only if used with operator
            that requires :math:`logq` of an approximation
            Returns Scalar

    Notes
    -----
    There are some defaults for approximation classes that can be
    optionally overriden.
        - `initial_dist_name`
            string that represents name of the initial distribution.
            In most cases if will be `uniform` or `normal`
        - `initial_dist_map`
            float where initial distribution has maximum density

    References
    ----------
    - Geoffrey Roeder, Yuhuai Wu, David Duvenaud, 2016
        Sticking the Landing: A Simple Reduced-Variance Gradient for ADVI
        approximateinference.org/accepted/RoederEtAl2016.pdf

    - Kingma, D. P., & Welling, M. (2014).
      Auto-Encoding Variational Bayes. stat, 1050, 1.
    """
    initial_dist_name = 'normal'
    initial_dist_map = 0.

    def __init__(self, local_rv=None, model=None, cost_part_grad_scale=1, **kwargs):
        model = modelcontext(model)
        self.model = model
        self.check_model(model, **kwargs)
        if local_rv is None:
            local_rv = {}

        def get_transformed(v):
            if hasattr(v, 'transformed'):
                return v.transformed
            return v

        known = {get_transformed(k): v for k, v in local_rv.items()}
        self.known = known
        self.local_vars = self.get_local_vars(**kwargs)
        self.global_vars = self.get_global_vars(**kwargs)
        self.order = ArrayOrdering(self.local_vars + self.global_vars)
        self.flat_view = model.flatten(
            vars=self.local_vars + self.global_vars
        )
        self.grad_scale_op = GradScale(cost_part_grad_scale)
        self._setup(**kwargs)
        self.shared_params = self.create_shared_params(**kwargs)

    @property
    def normalizing_constant(self):
        t = self.to_flat_input(tt.max([v.scaling for v in self.model.basic_RVs]))
        t = theano.clone(t, {self.input: tt.zeros(self.total_size)})
        return t

    def _setup(self, **kwargs):
        pass

    def get_global_vars(self, **kwargs):
        return [v for v in self.model.free_RVs if v not in self.known]

    def get_local_vars(self, **kwargs):
        return [v for v in self.model.free_RVs if v in self.known]

    def __getstate__(self):
        state = self.__dict__.copy()
        # can be inferred from the rest parts
        state.pop('flat_view')
        state.pop('order')
        return state

    def __setstate__(self, state):
        self.__dict__.update(state)
        self.order = ArrayOrdering(self.local_vars + self.global_vars)
        self.flat_view = self.model.flatten(
            vars=self.local_vars + self.global_vars
        )

    _view = property(lambda self: self.flat_view.view)
    input = property(lambda self: self.flat_view.input)

    def check_model(self, model, **kwargs):
        """
        Checks that model is valid for variational inference
        """
        vars_ = [var for var in model.vars if not isinstance(var, pm.model.ObservedRV)]
        if any([var.dtype in pm.discrete_types for var in vars_]):  # pragma: no cover
            raise ValueError('Model should not include discrete RVs')

    def create_shared_params(self, **kwargs):
        """
        Returns
        -------
        {dict|list|theano.shared}
        """
        pass

    def _local_mu_rho(self):
        mu = []
        rho = []
        for var in self.local_vars:
            mu.append(self.known[var][0].ravel())
            rho.append(self.known[var][1].ravel())
        mu = tt.concatenate(mu)
        rho = tt.concatenate(rho)
        return mu, rho

    def construct_replacements(self, include=None, exclude=None,
                               more_replacements=None):
        """
        Construct replacements with given conditions

        Parameters
        ----------
        include : list
            latent variables to be replaced
        exclude : list
            latent variables to be excluded for replacements
        more_replacements : dict
            add custom replacements to graph, e.g. change input source

        Returns
        -------
        dict
            Replacements
        """
        if include is not None and exclude is not None:
            raise ValueError('Only one parameter is supported {include|exclude}, got two')
        if include is not None:    # pragma: no cover
            replacements = {k: v for k, v
                            in self.flat_view.replacements.items() if k in include}
        elif exclude is not None:  # pragma: no cover
            replacements = {k: v for k, v
                            in self.flat_view.replacements.items() if k not in exclude}
        else:
            replacements = self.flat_view.replacements.copy()
        if more_replacements is not None:   # pragma: no cover
            replacements.update(more_replacements)
        return replacements

    def apply_replacements(self, node, deterministic=False,
                           include=None, exclude=None,
                           more_replacements=None):
        """
        Replace variables in graph with variational approximation. By default, replaces all variables

        Parameters
        ----------
        node : Theano Variables (or Theano expressions)
            node or nodes for replacements
        deterministic : bool
            whether to use zeros as initial distribution
            if True - zero initial point will produce constant latent variables
        include : list
            latent variables to be replaced
        exclude : list
            latent variables to be excluded for replacements
        more_replacements : dict
            add custom replacements to graph, e.g. change input source

        Returns
        -------
        node(s) with replacements
        """
        replacements = self.construct_replacements(
            include, exclude, more_replacements
        )
        node = theano.clone(node, replacements, strict=False)
        posterior = self.random(no_rand=deterministic)
        return theano.clone(node, {self.input: posterior}, strict=False)

    def sample_node(self, node, size=100,
                    more_replacements=None):
        """
        Samples given node or nodes over shared posterior

        Parameters
        ----------
        node : Theano Variables (or Theano expressions)
        size : scalar
            number of samples
        more_replacements : dict
            add custom replacements to graph, e.g. change input source

        Returns
        -------
        sampled node(s) with replacements
        """
        if more_replacements is not None:   # pragma: no cover
            node = theano.clone(node, more_replacements, strict=False)
        posterior = self.random(size)
        node = self.to_flat_input(node)

        def sample(z): return theano.clone(node, {self.input: z}, strict=False)
        nodes, _ = theano.scan(sample, posterior, n_steps=size)
        return nodes

    def scale_grad(self, inp):
        """
        Rescale gradient of input

        References
        ----------
        - Geoffrey Roeder, Yuhuai Wu, David Duvenaud, 2016
            Sticking the Landing: A Simple Reduced-Variance Gradient for ADVI
            approximateinference.org/accepted/RoederEtAl2016.pdf
        """
        return self.grad_scale_op(inp)

    def to_flat_input(self, node):
        """
        Replaces vars with flattened view stored in self.input
        """
        return theano.clone(node, self.flat_view.replacements, strict=False)

    @property
    def params(self):
        return cast_to_list(self.shared_params)

    def initial(self, size, no_rand=False, l=None):
        """
        Initial distribution for constructing posterior

        Parameters
        ----------
        size : int - number of samples
        no_rand : bool - return zeros if True
        l : length of sample, defaults to latent space dim

        Returns
        -------
        Tensor
            sampled latent space shape == size + latent_dim
        """

        theano_condition_is_here = isinstance(no_rand, tt.Variable)
        if l is None:   # pragma: no cover
            l = self.total_size
        if size is None:
            shape = (l, )
        else:
            shape = (size, l)
        shape = tt.stack(*shape)
        if theano_condition_is_here:
            no_rand = tt.as_tensor(no_rand)
            sample = getattr(tt_rng(), self.initial_dist_name)(shape)
            space = tt.switch(
                no_rand,
                tt.ones_like(sample) * self.initial_dist_map,
                sample
            )
        else:
            if no_rand:
                return tt.ones(shape) * self.initial_dist_map
            else:
                return getattr(tt_rng(), self.initial_dist_name)(shape)
        return space

    def random_local(self, size=None, no_rand=False):
        """
        Implements posterior distribution from initial latent space

        Parameters
        ----------
        size : number of samples from distribution
        no_rand : whether use deterministic distribution

        Returns
        -------
        local posterior space
        """

        mu, rho = self._local_mu_rho()
        e = self.initial(size, no_rand, self.local_size)
        return e * rho2sd(rho) + mu

    def random_global(self, size=None, no_rand=False):  # pragma: no cover
        """
        Implements posterior distribution from initial latent space

        Parameters
        ----------
        size : number of samples from distribution
        no_rand : whether use deterministic distribution

        Returns
        -------
        global posterior space
        """
        raise NotImplementedError

    def random(self, size=None, no_rand=False):
        """
        Implements posterior distribution from initial latent space

        Parameters
        ----------
        size : number of samples from distribution
        no_rand : whether use deterministic distribution

        Returns
        -------
        posterior space (theano)
        """
        if size is None:
            ax = 0
        else:
            ax = 1
        if self.local_vars and self.global_vars:
            return tt.concatenate([
                self.random_local(size, no_rand),
                self.random_global(size, no_rand)
            ], axis=ax)
        elif self.local_vars:   # pragma: no cover
            return self.random_local(size, no_rand)
        elif self.global_vars:
            return self.random_global(size, no_rand)
        else:   # pragma: no cover
            raise ValueError('No FreeVARs in model')

    @property
    @memoize
    @change_flags(compute_test_value='off')
    def random_fn(self):
        """
        Implements posterior distribution from initial latent space

        Parameters
        ----------
        size : number of samples from distribution
        no_rand : whether use deterministic distribution

        Returns
        -------
        posterior space (numpy)
        """
        In = theano.In
        size = tt.iscalar('size')
        no_rand = tt.bscalar('no_rand')
        posterior = self.random(size, no_rand)
        fn = theano.function([In(size, 'size', 1, allow_downcast=True),
                              In(no_rand, 'no_rand', 0, allow_downcast=True)],
                             posterior)

        def inner(size=None, no_rand=False):
            if size is None:
                return fn(1, int(no_rand))[0]
            else:
                return fn(size, int(no_rand))

        return inner

    def sample_vp(self, draws=1, hide_transformed=False):
        """
        Draw samples from variational posterior.

        Parameters
        ----------
        draws : int
            Number of random samples.
        hide_transformed : bool
            If False, transformed variables are also sampled. Default is True.

        Returns
        -------
        trace : pymc3.backends.base.MultiTrace
            Samples drawn from the variational posterior.
        """
        if hide_transformed:
            vars_sampled = [v_ for v_ in self.model.unobserved_RVs
                            if not str(v_).endswith('_')]
        else:
            vars_sampled = [v_ for v_ in self.model.unobserved_RVs]
        posterior = self.random_fn(draws)
        names = [var.name for var in self.local_vars + self.global_vars]
        samples = {name: self.view(posterior, name)
                   for name in names}

        def points():
            for i in range(draws):
                yield {name: samples[name][i]
                       for name in names}

        trace = pm.sampling.NDArray(model=self.model, vars=vars_sampled)
        try:
            trace.setup(draws=draws, chain=0)
            for point in points():
                trace.record(point)
        finally:
            trace.close()
        return pm.sampling.MultiTrace([trace])

    def log_q_W_local(self, z):
        """
        log_q_W samples over q for local vars
        Gradient wrt mu, rho in density parametrization
        is set to zero to lower variance of ELBO
        """
        if not self.local_vars:
            return tt.constant(0)
        mu, rho = self._local_mu_rho()
        mu = self.scale_grad(mu)
        rho = self.scale_grad(rho)
        logp = log_normal(z[self.local_slc], mu, rho=rho)
        scaling = []
        for var in self.local_vars:
            scaling.append(tt.ones(var.dsize)*var.scaling)
        scaling = tt.concatenate(scaling)
        if z.ndim > 1:  # pragma: no cover
            # rare case when logq(z) is called directly
            logp *= scaling[None]
        else:
            logp *= scaling
        return self.to_flat_input(tt.sum(logp))

    def log_q_W_global(self, z):    # pragma: no cover
        """
        log_q_W samples over q for global vars
        """
        raise NotImplementedError

    def logq(self, z):
        """
        Total logq for approximation
        """
        return self.log_q_W_global(z) + self.log_q_W_local(z)

    def logq_norm(self, z):
        return self.logq(z) / self.normalizing_constant

    def view(self, space, name, reshape=True):
        """
        Construct view on a variable from flattened `space`

        Parameters
        ----------
        space : space to take view of variable from
        name : str
            name of variable
        reshape : bool
            whether to reshape variable from vectorized view

        Returns
        -------
        variable view
        """
        theano_is_here = isinstance(space, tt.TensorVariable)
        slc = self._view[name].slc
        _, _, _shape, dtype = self._view[name]
        if space.ndim == 2:
            view = space[:, slc]
        elif space.ndim < 2:
            view = space[slc]
        else:   # pragma: no cover
            raise ValueError('Space should have no more than 2 dims, got %d' % space.ndim)
        if reshape:
            if len(_shape) > 0:
                if theano_is_here:
                    shape = tt.concatenate([space.shape[:-1],
                                            tt.as_tensor(_shape)])
                else:
                    shape = np.concatenate([space.shape[:-1],
                                            _shape]).astype(int)

            else:
                shape = space.shape[:-1]
            if theano_is_here:
                view = view.reshape(shape, ndim=space.ndim + len(_shape) - 1)
            else:
                view = view.reshape(shape)
        return view.astype(dtype)

    @property
    def total_size(self):
        return self.order.dimensions

    @property
    def local_size(self):
        size = sum([0] + [v.dsize for v in self.local_vars])
        return size

    @property
    def global_size(self):
        return self.total_size - self.local_size

    @property
    def local_slc(self):
        return slice(0, self.local_size)

    @property
    def global_slc(self):
        return slice(self.local_size, self.total_size)

reason: import numpy,import theano,import pymc3


text:
# MIT License (check the LICENSE file)
# Copyright (c) Author of "Balancing detectability and performance
# of attacks on the control channel of Markov Decision Processes",
# Alessio Russo, alessior@kth.se, 2021.
#

import cvxpy as cp
import numpy as np
from scipy.linalg import solve_discrete_lyapunov, solve_discrete_are
from collections import namedtuple
System = namedtuple('System', ['A', 'B', 'K', 'Kbar', 'R', 'theta', 'Q'])

def dlqr(A: np.ndarray, B: np.ndarray, Q: np.ndarray, R: np.ndarray) -> np.ndarray:
    ''' Compute optimal feedback gain for a discrete-time system '''
    P = solve_discrete_are(A, B, Q, R)
    return np.linalg.inv(B.T @ P @ B + R) * (B.T @ P @ A)

def compute_power_series(X: np.ndarray, Z: np.ndarray = None, eps: float = 1e-5):
    ''' Function used to compute Y=\\sum_{i=0}^\\infty (X^T)^i Z X^i
        eps tunes the approximated value of the series.
        Returns Y and the value of i for which the series is truncated.
    '''
    L_bar = np.identity(X.shape[0])
    L_pwr = np.matrix(np.identity(X.shape[0]))
    n = 0
    if Z is None:
        Z = np.matrix(np.identity(X.shape[0]))
    elif np.all(np.isclose(Z, 0.)):
        return 0*X, 0, 0
    while np.linalg.norm(L_pwr) > eps:
        L_pwr = L_pwr @ (X.T @ Z @ X)
        L_bar += L_pwr
        n += 1
    return L_bar, L_pwr, n

def compute_R(sys: System, beta: float, Ebar: np.ndarray) -> np.ndarray:
    ''' Computes the optimal attack covariance matrix according to the
        value of beta, Ebar, Kbar
    '''
    Binv = np.linalg.inv(sys.B.T @ sys.B) @ sys.B.T @ sys.Q
    Qinv = np.linalg.inv(sys.Q)
    E1 = 0.5 * beta * compute_power_series(Ebar, sys.Kbar.T @ sys.B.T @ Qinv @ sys.B @ sys.Kbar)[0]
    E2 = compute_power_series(Ebar)[0]

    Ebarinv = np.linalg.inv(E1 -E2)
    R = - Binv @ np.linalg.inv(0.5 * beta * Ebarinv + sys.Q) @ Binv.T
    return R

def compute_llr(x: float,  next_x: float, sys: System) -> float:
    ''' Computes log-likelihood ratio '''
    n,m = sys.B.shape
    mu0 = (sys.A + sys.B @ sys.K) @ x
    mu1 = mu0 + sys.B @ sys.Kbar @ x
    Sigma1 = sys.Q + sys.B @ sys.R @ sys.B.T
    Sigma0 = sys.Q
    z = next_x[None, :]

    t1 = -0.5 * (z- mu1) @ np.linalg.inv(Sigma1) @ (z - mu1).T
    t1 -= 0.5 * np.log(np.linalg.det(Sigma1))
    t0 = -0.5 * (z - mu0) @ np.linalg.inv(Sigma0) @ (z - mu0).T
    t0 -= 0.5 * np.log(np.linalg.det(Sigma0))

    return (t1 - t0).item()

def simulate(N: int, T: int, sys: System, attack_time: int = 0) -> (np.ndarray, np.ndarray):
    ''' Performs N simulations of the given system, each of T steps
        returns the states and actions
    '''
    n, m = sys.B.shape
    x = np.zeros((N, T+1, n))
    u = np.zeros((N, T, m))
    llr =  np.zeros((N, T))
    cusum = np.zeros((N, T+1))
    for i in range(N):
        x[i, 0] = np.random.normal(size=(n,))
        for t in range(T):
            u[i, t] = sys.K @ x[i, t]
            if t >= attack_time:
                u[i, t] += (sys.Kbar @ x[i, t] + np.random.multivariate_normal(np.zeros(m), sys.R)).tolist()[0]
            # u[i, t, :] = ((sys.K + sys.Kbar) @ x[i, t,:] + np.random.multivariate_normal(np.zeros(m), sys.R)).flatten()
            x[i, t+1,:] = sys.A@x[i, t] + sys.B@u[i, t] + np.random.multivariate_normal(np.zeros(n), sys.Q)
            llr[i, t] = compute_llr(x[i, t], x[i, t+1], sys)
            cusum[i, t+1] = max(0,  cusum[i, t]+llr[i, t])
    return x, u, llr, cusum

def compute_beta_theta(sys: System, Linv: np.ndarray) -> float:
    ''' Compute the minimum value of beta for theta '''
    I = np.identity(sys.A.shape[0])
    Qinv = np.linalg.inv(sys.Q)
    x = sys.B.T @ Linv.T @ Linv @ sys.B
    y = sys.B.T @ (I + sys.B @ sys.Kbar @ Linv).T @ Qinv @ (I + sys.B @ sys.Kbar @ Linv) @ sys.B
    return 2 * np.linalg.eig(x)[0].max() / np.linalg.eig(y)[0].min()

def compute_beta_R(sys: System, Enum: np.ndarray, Eden: np.ndarray) -> (float, float):
    ''' Compute the values beta_min, beta_max for the optimal attack
        covariance matrix R
    '''
    beta = cp.Variable(1)
    I = np.identity(sys.A.shape[0])
    Ebeta = (0.5 * beta * Eden - Enum)

    constraints = [-Ebeta >> 1e-9 * I, beta >= 1e-5, 0.5 * beta * I + sys.Q @ Ebeta >> 1e-9*I]
    prob = cp.Problem(cp.Minimize(beta), constraints)
    prob.solve()

    beta_min = beta.value.item()
    prob = cp.Problem(cp.Maximize(beta), constraints)
    prob.solve()
    if beta.value is None:
        beta_max = np.infty
    else:
        beta_max = beta.value.item()
    return beta_min, beta_max

def compute_I(sys: System, Linv: np.ndarray) -> float:
    ''' Compute information rate '''
    I = np.identity(sys.A.shape[0])
    Qinv = np.linalg.inv(sys.Q)
    term1 = np.trace(Qinv @ sys.B @ sys.R @ sys.B.T)
    term2 = np.log(np.linalg.det(sys.Q) / np.linalg.det(sys.Q + sys.B @ sys.R @ sys.B.T))
    x = sys.B.T @ (I +  sys.B @ sys.Kbar @ Linv).T @ Qinv @ (I +  sys.B @ sys.Kbar @ Linv) @ sys.B
    term3 = sys.theta.T @ x @ sys.theta
    Sigma = solve_discrete_lyapunov(sys.A + sys.B @ sys.K + sys.B @ sys.Kbar, sys.Q + sys.B @ sys.R @ sys.B.T)
    term4 = np.trace(sys.Kbar.T @ sys.B.T @ Qinv @ sys.B @ sys.Kbar @ Sigma)
    return 0.5 * (term1 + term2 + term3 + term4)

def compute_attack(sys: System, debug: bool = False):
    ''' Computes R, beta_min, beta_max '''
    # Eigenvalues
    I = np.identity(sys.A.shape[0])
    E = sys.A + sys.B @ sys.K
    Ebar = sys.A + sys.B @ (sys.K + sys.Kbar)
    L = I - Ebar
    Linv = np.linalg.inv(L)
    Qinv = np.linalg.inv(sys.Q)
    eigE = np.linalg.eig(E)[0]
    eigEbar = np.linalg.eig(Ebar)[0]
    if debug:
        print('Eigenvalues of A+BK: {}'.format(eigE))
        print('Eigenvalues of A+BK+BKbar: {}'.format(eigEbar))

    # Check condition on R
    Epower_num = compute_power_series(Ebar)[0]
    Epower_den = compute_power_series(Ebar, sys.Kbar.T @ sys.B.T @ Qinv @ sys.B @ sys.Kbar)[0]

    beta_min, beta_max = compute_beta_R(sys, Epower_num, Epower_den)
    beta_theta = compute_beta_theta(sys, Linv)

    # Compute R

    beta = beta_min * (1 + 1e-3)
    R = compute_R(sys, beta, Ebar)
    if debug:
        print('Chosen beta: {} - {} - beta_theta: {}'.format(beta_min, beta_max, beta_theta))
        print('Matrix R: \n{}'.format(R))
        print('Eigenvalues of R: {}'.format(np.linalg.eig(R)[0]))

    if beta_max < np.infty:
        beta_min = beta_min + (beta_max-beta_min)/50
    return R, (beta_min, beta_max), beta, Linv, Ebar

def compute_ergodic_reward(sys: System, Linv: np.ndarray) -> float:
    ''' Computes the ergodic reward of a linear system '''
    term1 = sys.theta.T @ (sys.B.T @ Linv.T @ Linv @ sys.B) @ sys.theta
    term2 = np.trace(
            solve_discrete_lyapunov(
                sys.A + sys.B @ sys.K + sys.B @ sys.Kbar, sys.Q + sys.B @ sys.R @ sys.B.T)
            )
    return term1 + term2

ext: py

reason: import numpy,from scipy,import cvxpy


text:
from tensorflow.keras.utils import Sequence
import math
import numpy as np
import cv2
import gc
#from skimage.io import imread
import tensorflow as tf
from keras import backend

class Dataloader(Sequence):
    def __init__(self,x_set,y_set,batch_size):
        self.x, self.y = x_set, y_set
        self.batch_size = batch_size
        self.indices = np.arange(len(self.x))
        self.längd = len(self.x)
        self.antalPositiva = 0
    def __len__(self):
        return math.ceil(len(self.x) / self.batch_size)

    def createStacks(self,batch_x,batch_y,j,index):
        train_X = list()
        train_y = list()
        sista = len(batch_x)
        if j == 0:
            for i in range(self.batch_size):
                train_X.append(batch_x[i:i + self.batch_size])
                train_y.append(batch_y[i + self.batch_size])
        else:
            for i in range(self.batch_size):
                if(i < (32 - j - 16)):
                    train_X.append(batch_x[i:i + self.batch_size])
                    train_y.append(batch_y[i + self.batch_size])
                else:
                    clip_x = list()
                    for k in range(self.batch_size):
                        clip_x.append(batch_x[sista - 1])
                    train_X.append(clip_x)
                    train_y.append(batch_y[sista - 1])
        train_X = np.array(train_X)
        #if train_X.shape != (16,16,240,320,3):
        #    print(train_X.shape)
        train_X = train_X.reshape(train_X.shape[0],train_X.shape[2],train_X.shape[3],train_X.shape[4],train_X.shape[1])
        train_X = train_X.astype("float32") / 255
        train_y = np.array(train_y)
        #if train_y.max == 1:
        #    print(1)
        return train_X, train_y

    def __getitem__(self, index):
        #längd =len(self) - 2
        #index = längd
        #längd2 =
        i = (index + 2) * self.batch_size
        #l = (längd2 + 2) * self.batch_size
        #k = längd * self.batch_size
        j = 0
        if i > (self.längd -2):
            j = i - len(self.x)
            #print(j)
            #print(index)
        batch_x = self.x[index * self.batch_size:(index + 2) * self.batch_size - j]
        batch_y = self.y[index * self.batch_size:(index + 2) * self.batch_size - j]
        batch_x = np.array([cv2.imread(file_name)
               for file_name in batch_x])
        batch_y = np.array(batch_y)
        #if(batch_y.max == 1):
        #    self.antalPositiva += 1
        batch_x,batch_y = self.createStacks(batch_x,batch_y,j,i)
        return batch_x, batch_y

    #def on_epoch_end(self) -> None:
        #gc.collect()
        #np.random.shuffle(self.indices)

max_stars_repo_name: evinus/My-appproch-One

max_issues_repo_name: evinus/My-appproch-One

max_forks_repo_name: evinus/My-appproch-One

reason: import numpy


text:
import sys
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import re
from matplotlib.ticker import MultipleLocator
from matplotlib import gridspec

#majorLocatorX = MultipleLocator(2)
#minorLocatorX = MultipleLocator(1)
#majorLocatorY = MultipleLocator(0.05)
#minorLocatorY = MultipleLocator(0.025)

filename1 = '/home/sam/Documents/thesis/data/15N_levels.dat'
filename2 = '/home/sam/Documents/thesis/data/PR_final_7_8.dat'

e1 = []
jp1 = []
j1 = []
p1 = []

e2 = [[] for num in range(3)]
jp2 = [[] for num in range(3)]
j2 = [[] for num in range(3)]
p2 = [[] for num in range(3)]

with open(filename1) as f1:
    data1 = f1.read()
data1 = data1.split('\n')

with open(filename2) as f2:
    data2 = f2.read()
data2 = data2.split('\n')

num_states = 15
num = 1
#while(num < num_states and num < len(data1)-1):
for num in [1,4]:
    line = data1[num].split("\t")
    e1.append([float(line[4])/1000.0,float(line[4])/1000.0])
    jp = str(line[6])
    jp = jp.replace('(','').replace(')','').replace(' ','')

    if( len(jp) != 0 ):
        if(jp not in jp1):
            jp1.append(jp)
        else:
            jp1.append('')
    else:
        if('J^{\pi}' not in jp1):
            jp1.append('J^{\pi}')
        else:
            jp1.append('')

    if( len(jp) == 0 ):
        p1.append('-')
    elif( jp[-1] == '-' ):
        p1.append('--')
    else:
        p1.append('-')

    if( len(jp) == 0 ):
        j1.append('k')
    else:
        jp = jp[:-1]

    if( jp == '1/2' ):
        j1.append('b')
    elif( jp == '3/2' ):
        j1.append('r')
    elif( jp == '5/2' ):
        j1.append('g')
    elif( jp == '7/2' ):
        j1.append('m')
    elif( jp == '9/2' ):
        j1.append('c')

    num = num + 1

indices = {0:0, 1:1, 4:2}

num1 = [1, 1, 1]
num2 = 1
while((num1[0] < num_states or num1[1] < num_states or num1[2] < num_states) and num2 < len(data2)-1):
    line = data2[num2].split()
    if(  float(line[0]) != 16.0 ):
        num2 = num2 + 1
        continue
    index = indices[int(line[2])]
    e2[index].append([float(line[11]),float(line[11])])
    jp = str(line[6])[1:-1]

    if( jp[0] == 'p' or jp[0] == 'f'):
        p2[index].append('--')
        jp = jp + '-'
    else:
        p2[index].append('-')
        jp = jp + '+'

    jp = jp[1:]

    if(jp not in jp2[0] and jp not in jp2[1] and jp not in jp2[2] and jp not in jp1):
        jp2[index].append(jp)
    else:
        jp2[index].append('')

    if( jp[:-1] == '1/2' ):
        j2[index].append('b')
    elif( jp[:-1] == '3/2' ):
        j2[index].append('r')
    elif( jp[:-1] == '5/2' ):
        j2[index].append('g')
    elif( jp[:-1] == '7/2' ):
        j2[index].append('m')
    elif( jp[:-1] == '9/2' ):
        j2[index].append('c')

    num1[index] = num1[index] + 1
    num2 = num2 + 1

filename3 = '/home/sam/Documents/thesis/data/15O_levels.dat'
filename4 = '/home/sam/Documents/thesis/data/PR_final_8_7.dat'

e3 = []
jp3 = []
j3 = []
p3 = []

e4 = [[] for num in range(3)]
jp4 = [[] for num in range(3)]
j4 = [[] for num in range(3)]
p4 = [[] for num in range(3)]

with open(filename3) as f3:
    data3 = f3.read()
data3 = data3.split('\n')

with open(filename4) as f4:
    data4 = f4.read()
data4 = data4.split('\n')

num_states = 15
num = 1
#while(num < num_states and num < len(data3)-1):
for num in [1,4]:
    line = data3[num].split("\t")
    e3.append([float(line[4])/1000.0,float(line[4])/1000.0])
    jp = str(line[6])
    jp = jp.replace('(','').replace(')','').replace(' ','')

    if( len(jp) != 0 ):
        if(jp not in jp3):
            jp3.append(jp)
        else:
            jp3.append('')
    else:
        if('J^{\pi}' not in jp3):
            jp3.append('J^{\pi}')
        else:
            jp3.append('')

    if( len(jp) == 0 ):
        p3.append('-')
    elif( jp[-1] == '-' ):
        p3.append('--')
    else:
        p3.append('-')

    if( len(jp) == 0 ):
        j3.append('k')
    else:
        jp = jp[:-1]

    if( jp == '1/2' ):
        j3.append('b')
    elif( jp == '3/2' ):
        j3.append('r')
    elif( jp == '5/2' ):
        j3.append('g')
    elif( jp == '7/2' ):
        j3.append('m')
    elif( jp == '9/2' ):
        j3.append('c')

    num = num + 1

indices = {0:0, 1:1, 4:2}

num1 = [1, 1, 1]
num2 = 1
while((num1[0] < num_states or num1[1] < num_states or num1[2] < num_states) and num2 < len(data4)-1):
    line = data4[num2].split()
    if(  float(line[0]) != 16.0 ):
        num2 = num2 + 1
        continue
    index = indices[int(line[2])]
    e4[index].append([float(line[11]),float(line[11])])
    jp = str(line[6])[1:-1]

    if( jp[0] == 'p' or jp[0] == 'f'):
        p4[index].append('--')
        jp = jp + '-'
    else:
        p4[index].append('-')
        jp = jp + '+'

    jp = jp[1:]

    if(jp not in jp4[0] and jp not in jp4[1] and jp not in jp4[2] and jp not in jp3):
        jp4[index].append(jp)
    else:
        jp4[index].append('')

    if( jp[:-1] == '1/2' ):
        j4[index].append('b')
    elif( jp[:-1] == '3/2' ):
        j4[index].append('r')
    elif( jp[:-1] == '5/2' ):
        j4[index].append('g')
    elif( jp[:-1] == '7/2' ):
        j4[index].append('m')
    elif( jp[:-1] == '9/2' ):
        j4[index].append('c')

    num1[index] = num1[index] + 1
    num2 = num2 + 1

#xa = [[-14.0,-10.0], [-6.0,-2.0], [2.0,6.0], [10.0,14.0]]
xa = [[-12.0,-6.0], [-3.0,3.0], [2.0,6.0], [6.0,12.0]]

plt.rc('font', family='serif')

fig = plt.figure(figsize=(16, 6))
gs = gridspec.GridSpec(1, 2)

ax1 = plt.subplot(gs[0])
for num in range(len(e1)):
    if( jp1[num] == '' ):
        labelstring = ""
    else:
        labelstring = r'$\mathrm{' + jp1[num] + '}$'

    if(p1[num] == '--'):
        ax1.plot(xa[3], e1[num], linestyle=p1[num], dashes=(7,1), color=j1[num], linewidth=2.5, label=labelstring)
    else:
        ax1.plot(xa[3], e1[num], linestyle=p1[num], color=j1[num], linewidth=2.5, label=labelstring)

for index in range(3):
    offset = e2[index][0][0]
    for num in range(len(e2[index])):
        e2[index][num][0] -= offset
        e2[index][num][1] -= offset

for index in [0,1]:
    for num in range(len(e2[index])):
        if( jp2[index][num] == '' ):
            labelstring = ""
        else:
            labelstring = r'$\mathrm{' + jp2[index][num] + '}$'

        if(p2[index][num] == '--'):
            ax1.plot(xa[index], e2[index][num], linestyle=p2[index][num], dashes=(7,1), color=j2[index][num], linewidth=2.5, label=labelstring)
        else:
            ax1.plot(xa[index], e2[index][num], linestyle=p2[index][num], color=j2[index][num], linewidth=2.5, label=labelstring)

ax1.plot([xa[0][1], xa[1][0]], [e2[0][0], e2[1][0]], linestyle=':', color='b', linewidth=2.5)
ax1.plot([xa[1][1], xa[3][0]], [e2[1][0], e1[0]], linestyle=':', color='b', linewidth=2.5)

ax1.plot([xa[0][1], xa[1][0]], [e2[0][2], e2[1][1]], linestyle=':', color='r', linewidth=2.5)
ax1.plot([xa[1][1], xa[3][0]], [e2[1][1], e1[1]], linestyle=':', color='r', linewidth=2.5)

#ax1.plot([xa[0][1], xa[1][0]], [e2[0][4], e2[1][2]], linestyle=':', color='r', linewidth=2.5)
#ax1.plot([xa[1][1], xa[3][0]], [e2[1][2], e1[2]], linestyle=':', color='r', linewidth=2.5)

ax1.axis([-13.0, 15.0, -0.5, 12.0], fontsize=22)
ax1.tick_params(axis='x',labelsize=20)
ax1.tick_params(axis='y',labelsize=16)
ax1.set_ylabel(r'$\mathrm{E\ (MeV)}$', fontsize=20)
annotation_string = r'$\mathrm{^{15}N}$'
ax1.annotate(annotation_string, fontsize=24, xy=(0.9, 0.05), xycoords='axes fraction')
plt.setp(ax1, xticks=[-9.0,0.0,9.0], xticklabels=[r'$\beta=0$',r'$\beta=1$',r'$\mathrm{Exp.}$'])

ax1.legend(loc='upper right', frameon=False, fontsize=18, labelspacing=0.0)

ax2 = plt.subplot(gs[1])
for num in range(len(e3)):
    if( jp3[num] == '' ):
        labelstring = ""
    else:
        labelstring = r'$\mathrm{' + jp3[num] + '}$'

    if(p3[num] == '--'):
        ax2.plot(xa[3], e3[num], linestyle=p3[num], dashes=(7,1), color=j3[num], linewidth=2.5, label=labelstring)
    else:
        ax2.plot(xa[3], e3[num], linestyle=p3[num], color=j3[num], linewidth=2.5, label=labelstring)

for index in range(3):
    offset = e4[index][0][0]
    for num in range(len(e4[index])):
        e4[index][num][0] -= offset
        e4[index][num][1] -= offset

for index in [0,1]:
    for num in range(len(e4[index])):
        if( jp4[index][num] == '' ):
            labelstring = ""
        else:
            labelstring = r'$\mathrm{' + jp4[index][num] + '}$'

        if(p4[index][num] == '--'):
            ax2.plot(xa[index], e4[index][num], linestyle=p4[index][num], dashes=(7,1), color=j4[index][num], linewidth=2.5, label=labelstring)
        else:
            ax2.plot(xa[index], e4[index][num], linestyle=p4[index][num], color=j4[index][num], linewidth=2.5, label=labelstring)

ax2.plot([xa[0][1], xa[1][0]], [e4[0][0], e4[1][0]], linestyle=':', color='b', linewidth=2.5)
ax2.plot([xa[1][1], xa[3][0]], [e4[1][0], e3[0]], linestyle=':', color='b', linewidth=2.5)

ax2.plot([xa[0][1], xa[1][0]], [e4[0][2], e4[1][1]], linestyle=':', color='r', linewidth=2.5)
ax2.plot([xa[1][1], xa[3][0]], [e4[1][1], e3[1]], linestyle=':', color='r', linewidth=2.5)

#ax2.plot([xa[0][1], xa[1][0]], [e4[0][4], e4[1][2]], linestyle=':', color='r', linewidth=2.5)
#ax2.plot([xa[1][1], xa[3][0]], [e4[1][2], e3[2]], linestyle=':', color='r', linewidth=2.5)

ax2.axis([-13.0, 15.0, -0.5, 12.0], fontsize=22)
ax2.tick_params(axis='x',labelsize=20)
plt.setp(ax2.get_yticklabels(), visible=False)
annotation_string = r'$\mathrm{^{15}O}$'
ax2.annotate(annotation_string, fontsize=24, xy=(0.9, 0.05), xycoords='axes fraction')
plt.setp(ax2, xticks=[-9.0,0.0,9.0], xticklabels=[r'$\beta=0$',r'$\beta=1$',r'$\mathrm{Exp.}$'])

ax2.legend(loc='upper right', frameon=False, fontsize=18, labelspacing=0.0)

plt.tight_layout()
plt.savefig('N15_O15_EOM_2.pdf', format='pdf', bbox_inches='tight')
plt.show()

max_stars_repo_name: novarios/Thesis

max_issues_repo_name: novarios/Thesis

max_forks_repo_name: novarios/Thesis


text:
# Copyright 2018 Google LLC
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#      http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License

"""Tests for tensorflow_data_validation.arrow.decoded_examples_to_arrow."""

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

from absl.testing import absltest
from absl.testing import parameterized
import numpy as np
import pyarrow as pa
import six
from tensorflow_data_validation.arrow import arrow_util
from tensorflow_data_validation.arrow import decoded_examples_to_arrow

_INVALID_INPUT_TEST_CASES = [
    dict(
        testcase_name="list_of_non_dict",
        test_input=[1, 2],
        expected_error=ValueError,
        expected_error_regexp="Unexpected Arrow type created from input"),
    dict(
        testcase_name="list_of_dict_of_non_str_key",
        test_input=[{
            1: None
        }],
        expected_error=pa.ArrowTypeError,
        expected_error_regexp="Expected dict key of type str or bytes"),
    dict(
        testcase_name="unsupported_ndarray_type",
        test_input=[{
            "a": np.array([1j, 2j, 3j], dtype=np.complex64)
        }],
        expected_error=RuntimeError,
        expected_error_regexp="Unsupported numpy type"),
]

_CONVERSION_TEST_CASES = [
    dict(
        testcase_name="unicode_feature_name",
        input_examples=[{
            u"\U0001f951": np.array([1, 2, 3], dtype=np.int64),
        }],
        expected_output={
            u"\U0001f951": pa.array([[1, 2, 3]], type=pa.list_(pa.int64())),
        }),
    dict(
        testcase_name="supported_ndarray_types",
        input_examples=[
            {
                "int64_feature": np.array([1, 2, 3], dtype=np.int64),
                "uint64_feature": np.array([1, 2, 3], dtype=np.uint64),
                "int32_feature": np.array([1, 2, 3], dtype=np.int32),
                "uint32_feature": np.array([1, 2, 3], dtype=np.uint32),
                "float_feature": np.array([1.], dtype=np.float32),
                "double_feature": np.array([1.], dtype=np.float64),
                "bytes_feature": np.array([b"abc", b"def"], dtype=np.object),
                "unicode_feature": np.array([u"abc", u"def"], dtype=np.object),
            },
            {
                "int64_feature": np.array([4], dtype=np.int64),
                "int32_feature": np.array([4], dtype=np.int32),
                "float_feature": np.array([2., 3., 4.], dtype=np.float32),
                "double_feature": np.array([2., 3., 4.], dtype=np.float64),
                "bytes_feature": np.array([b"ghi"], dtype=np.object),
                "unicode_feature": np.array([u"ghi"], dtype=np.object),
            },
        ],
        expected_output={
            "int64_feature":
                pa.array([[1, 2, 3], [4]], type=pa.list_(pa.int64())),
            "uint64_feature":
                pa.array([[1, 2, 3], None], type=pa.list_(pa.uint64())),
            "int32_feature":
                pa.array([[1, 2, 3], [4]], type=pa.list_(pa.int32())),
            "uint32_feature":
                pa.array([[1, 2, 3], None], type=pa.list_(pa.uint32())),
            "float_feature":
                pa.array([[1.], [2., 3., 4.]], type=pa.list_(pa.float32())),
            "double_feature":
                pa.array([[1.], [2., 3., 4.]], type=pa.list_(pa.float64())),
            "bytes_feature":
                pa.array([[b"abc", b"def"], [b"ghi"]],
                         type=pa.list_(pa.binary())),
            "unicode_feature":
                pa.array([[b"abc", b"def"], [b"ghi"]],
                         type=pa.list_(pa.string())),
        }),
    dict(
        testcase_name="mixed_unicode_and_bytes",
        input_examples=[
            {
                "a": np.array([b"abc"], dtype=np.object),
            },
            {
                "a": np.array([u"def"], dtype=np.object),
            },
        ],
        expected_output={
            "a":
                pa.array([[b"abc"], [b"def"]], type=pa.list_(pa.binary()))
        }),
    dict(
        testcase_name="none_feature_value",
        input_examples=[{
            "a": np.array([1, 2, 3], dtype=np.int64),
        }, {
            "a": None,
        }, {
            "a": None,
        }, {
            "a": np.array([4], dtype=np.int64),
        }],
        expected_output={
            "a":
                pa.array([[1, 2, 3], None, None, [4]],
                         type=pa.list_(pa.int64())),
        }),
    dict(
        testcase_name="empty_feature_value",
        input_examples=[{
            "a": np.array([], dtype=np.int64),
        }],
        expected_output={
            "a": pa.array([[]], type=pa.list_(pa.int64())),
        }),
    dict(
        testcase_name="missing_feature",
        input_examples=[{
            "f1": np.array([1, 2, 3], dtype=np.int64),
        }, {
            "f2": np.array([1., 2., 3.], dtype=np.float32),
        }, {
            "f3": np.array([b"abc", b"def"], dtype=np.object),
        }, {
            "f1": np.array([4, 5, 6], dtype=np.int64),
            "f4": np.array([8], dtype=np.int64),
        }],
        expected_output={
            "f1":
                pa.array([[1, 2, 3], None, None, [4, 5, 6]],
                         pa.list_(pa.int64())),
            "f2":
                pa.array([None, [1., 2., 3.], None, None],
                         pa.list_(pa.float32())),
            "f3":
                pa.array([None, None, [b"abc", b"def"], None],
                         pa.list_(pa.binary())),
            "f4":
                pa.array([None, None, None, [8]], pa.list_(pa.int64())),
        }),
    dict(
        testcase_name="null_array",
        input_examples=[{
            "a": None,
        }, {
            "a": None,
        }],
        expected_output={
            "a": pa.array([None, None], type=pa.null()),
        })
]

class DecodedExamplesToArrowPyTest(parameterized.TestCase):

  @parameterized.named_parameters(*_INVALID_INPUT_TEST_CASES)
  def test_invalid_input(self, test_input, expected_error,
                         expected_error_regexp):
    with self.assertRaisesRegex(expected_error, expected_error_regexp):
      decoded_examples_to_arrow.DecodedExamplesToRecordBatch(test_input)

  @parameterized.named_parameters(*_CONVERSION_TEST_CASES)
  def test_conversion(self, input_examples, expected_output):
    record_batch = decoded_examples_to_arrow.DecodedExamplesToRecordBatch(
        input_examples)
    self.assertLen(expected_output, record_batch.num_columns)
    for feature_name, expected_arrow_array in six.iteritems(expected_output):
      actual = arrow_util.get_column(record_batch, feature_name)
      self.assertTrue(
          expected_arrow_array.equals(actual),
          "{} vs {}".format(expected_arrow_array, actual))

  def test_conversion_empty_input(self):
    record_batch = decoded_examples_to_arrow.DecodedExamplesToRecordBatch([])
    self.assertEqual(record_batch.num_columns, 0)
    self.assertEqual(record_batch.num_rows, 0)

  def test_conversion_empty_examples(self):
    input_examples = [{}] * 10
    record_batch = decoded_examples_to_arrow.DecodedExamplesToRecordBatch(
        input_examples)
    self.assertEqual(record_batch.num_rows, 10)
    self.assertEqual(record_batch.num_columns, 0)

if __name__ == "__main__":
  absltest.main()

max_stars_repo_name: anukaal/data-validation

max_issues_repo_name: anukaal/data-validation

max_forks_repo_name: anukaal/data-validation


text:
import os
import cv2
import tqdm
import torch
import numpy as np
import seaborn as sns
import pandas as pd
import codecs
import json
import matplotlib.pyplot as plt
plt.switch_backend('agg')
from solver import Solver
from models.model import Model
from datasets.steel_dataset import provider
from utils.set_seed import seed_torch
from utils.cal_dice_iou import compute_dice_class
from config import get_seg_config

class ChooseThresholdMinArea():
    ''' 选择每一类的像素阈值和最小连通域
    '''
    def __init__(self, model, model_name, valid_loader, fold, save_path, class_num=4):
        ''' 模型初始化

        Args:
            model: 使用的模型
            model_name: 当前模型的名称
            valid_loader: 验证数据的Dataloader
            fold: 当前为多少折
            save_path: 保存结果的路径
            class_num: 有多少个类别
        '''
        self.model = model
        self.model_name = model_name
        self.valid_loader = valid_loader
        self.fold = fold
        self.save_path = save_path
        self.class_num = class_num

        self.model.eval()
        self.solver = Solver(model)

    def choose_threshold_minarea(self):
        ''' 采用网格法搜索各个类别最优像素阈值和最优最小连通域，并画出各个类别搜索过程中的热力图

        Return:
            best_thresholds_little: 每一个类别的最优阈值
            best_minareas_little: 每一个类别的最优最小连通取余
            max_dices_little: 每一个类别的最大dice值
        '''
        init_thresholds_range, init_minarea_range = np.arange(0.50, 0.71, 0.03), np.arange(768, 2305, 256)

        # 阈值列表和最小连通域列表，大小为 Nx4
        thresholds_table_big = np.array([init_thresholds_range, init_thresholds_range, \
                                         init_thresholds_range, init_thresholds_range])  # 阈值列表
        minareas_table_big = np.array([init_minarea_range, init_minarea_range, \
                                       init_minarea_range, init_minarea_range])  # 最小连通域列表

        f, axes = plt.subplots(figsize=(28.8, 18.4), nrows=2, ncols=self.class_num)
        cmap = sns.cubehelix_palette(start=1.5, rot=3, gamma=0.8, as_cmap=True)

        best_thresholds_big, best_minareas_big, max_dices_big = self.grid_search(thresholds_table_big, minareas_table_big, axes[0,:], cmap)
        print('best_thresholds_big:{}, best_minareas_big:{}, max_dices_big:{}'.format(best_thresholds_big, best_minareas_big, max_dices_big))

        # 开始细分类
        thresholds_table_little, minareas_table_little = list(), list()
        for best_threshold_big, best_minarea_big in zip(best_thresholds_big, best_minareas_big):
            thresholds_table_little.append(np.arange(best_threshold_big-0.03, best_threshold_big+0.03, 0.015))  # 阈值列表
            minareas_table_little.append(np.arange(best_minarea_big-256, best_minarea_big+257, 128))  # 像素阈值列表
        thresholds_table_little, minareas_table_little = np.array(thresholds_table_little), np.array(minareas_table_little)

        best_thresholds_little, best_minareas_little, max_dices_little = self.grid_search(thresholds_table_little, minareas_table_little, axes[1,:], cmap)
        print('best_thresholds_little:{}, best_minareas_little:{}, max_dices_little:{}'.format(best_thresholds_little, best_minareas_little, max_dices_little))

        f.savefig(os.path.join(self.save_path, self.model_name + '_fold'+str(self.fold)), bbox_inches='tight')
        # plt.show()
        plt.close()

        return best_thresholds_little, [float(x) for x in best_minareas_little], max_dices_little

    def grid_search(self, thresholds_table, minareas_table, axes, cmap):
        ''' 给定包含各个类别搜索区间的thresholds_table和minareas_table，求的各个类别的最优像素阈值，最优最小连通域，最高dice；
            并画出各个类别搜索过程中的热力图

        Args:
            thresholds_table: 待搜索的阈值范围，维度为[4, N]，numpy类型
            minareas_table: 待搜索的最小连通域范围，维度为[4, N]，numpy类型
            axes: 画各个类别搜索热力图时所需要的画柄，尺寸为[class_num]
            cmap: 画图时所需要的cmap

        return:
            best_thresholds: 各个类别的最优像素阈值，尺寸为[class_num]
            best_minareas: 各个类别的最优最小连通域，尺寸为[class_num]
            max_dices: 各个类别的最大dice，尺寸为[class_num]
        '''
        dices_table = np.zeros((self.class_num, np.shape(thresholds_table)[1], np.shape(minareas_table)[1]))
        tbar = tqdm.tqdm(self.valid_loader)
        with torch.no_grad():
            for i, samples in enumerate(tbar):
                if len(samples) == 0:
                    continue
                images, masks = samples[0], samples[1]
                # 完成网络的前向传播
                masks_predict_allclasses = self.solver.forward(images)
                dices_table += self.grid_search_batch(thresholds_table, minareas_table, masks_predict_allclasses, masks)

        dices_table = dices_table/len(tbar)
        best_thresholds, best_minareas, max_dices = list(), list(), list()
        # 处理每一类的预测结果
        for each_class, dices_oneclass_table in enumerate(dices_table):
            max_dice = np.max(dices_oneclass_table)
            max_location = np.unravel_index(np.argmax(dices_oneclass_table, axis=None),
                                            dices_oneclass_table.shape)
            best_thresholds.append(thresholds_table[each_class, max_location[0]])
            best_minareas.append(minareas_table[each_class, max_location[1]])
            max_dices.append(max_dice)

            data = pd.DataFrame(data=dices_oneclass_table, index=np.around(thresholds_table[each_class,:], 3), columns=minareas_table[each_class,:])
            sns.heatmap(data, linewidths=0.05, ax=axes[each_class], vmax=np.max(dices_oneclass_table), vmin=np.min(dices_oneclass_table), cmap=cmap,
                        annot=True, fmt='.4f')
            axes[each_class].set_title('search result')
        return best_thresholds, best_minareas, max_dices

    def grid_search_batch(self, thresholds_table, minareas_table, masks_predict_allclasses, masks_allclasses):
        '''给定thresholds、minareas矩阵、一个batch的预测结果和真实标签，遍历每个类的每一个组合得到对应的dice值

        Args:
            thresholds_table: 待搜索的阈值范围，维度为[4, N]，numpy类型
            minareas_table: 待搜索的最小连通域范围，维度为[4, N]，numpy类型
            masks_predict_allclasses: 所有类别的某个batch的预测结果且未经过sigmoid，维度为[batch_size, class_num, height, width]
            masks_allclasses: 所有类别的某个batch的真实类标，维度为[batch_size, class_num, height, width]

        Return:
            dices_table: 各个类别在其各自的所有搜索组合中所得到的dice值，维度为[4, M, N]
        '''

        # 得到每一个类别的搜索阈值区间和最小连通域搜索区间
        dices_table = list()
        for each_class, (thresholds_range, minareas_range) in enumerate(zip(thresholds_table, minareas_table)):
            # 得到每一类的预测结果和真实类标
            masks_predict_oneclass = masks_predict_allclasses[:, each_class, ...]
            masks_oneclasses = masks_allclasses[:, each_class, ...]

            dices_range = self.post_process(thresholds_range, minareas_range, masks_predict_oneclass, masks_oneclasses)
            dices_table.append(dices_range)
        # 得到的大小为 4 x len(thresholds_range) x len(minareas_range)
        return np.array(dices_table)

    def post_process(self, thresholds_range, minareas_range, masks_predict_oneclass, masks_oneclasses):
        '''给定某个类别的某个batch的数据，遍历所有搜索组合，得到每个组合的dice值

        Args:
            thresholds_range: 具体某个类别的像素阈值搜索区间，尺寸为[M]
            minareas_range: 具体某个类别的最小连通域搜索区间，尺寸为[N]
            masks_predict_oneclass: 预测出的某个类别的该batch的tensor向量且未经过sigmoid，维度为[batch_size, height, width]
            masks_oneclasses: 某个类别的该batch的真实类标，维度为[batch_size, height, width]

        Return:
            dices_range: 某个类别的该batch的所有搜索组合得到dice矩阵，维度为[M, N]
        '''

        # 注意，损失函数中包含sigmoid函数，一般情况下需要手动经过sigmoid函数
        masks_predict_oneclass = torch.sigmoid(masks_predict_oneclass).detach().cpu().numpy()
        dices_range = np.zeros((len(thresholds_range), len(minareas_range)))

        # 遍历每一个像素阈值和最小连通域
        for index_threshold, threshold in enumerate(thresholds_range):
            for index_minarea, minarea in enumerate(minareas_range):
                batch_preds = list()
                # 遍历每一张图片
                for pred in masks_predict_oneclass:
                    mask = cv2.threshold(pred, threshold, 1, cv2.THRESH_BINARY)[1]
                    # 将背景标记为 0，其他的块从 1 开始的正整数标记
                    num_component, component = cv2.connectedComponents(mask.astype(np.uint8))
                    predictions = np.zeros((256, 1600), np.float32)
                    num = 0
                    for c in range(1, num_component):
                        p = (component == c)
                        if p.sum() > minarea:
                            predictions[p] = 1
                            num += 1
                    batch_preds.append(predictions)
                dice = compute_dice_class(torch.from_numpy(np.array(batch_preds)), masks_oneclasses)
                dices_range[index_threshold, index_minarea] = dice
        return dices_range

def get_model(model_name, load_path):
    ''' 加载网络模型并加载对应的权重

    Args:
        model_name: 当前模型的名称
        load_path: 当前模型的权重路径

    Return:
        model: 加载出来的模型
    '''
    model = Model(model_name).create_model()
    Solver(model).load_checkpoint(load_path)
    return model

if __name__ == "__main__":
    config = get_seg_config()
    mean=(0.485, 0.456, 0.406)
    std=(0.229, 0.224, 0.225)
    mask_only = True

    dataloaders = provider(config.dataset_root, os.path.join(config.dataset_root, 'train.csv'), mean, std, config.batch_size, config.num_workers, config.n_splits, mask_only)
    results = {}
    # 存放权重的路径
    model_path = os.path.join(config.save_path, config.model_name)
    best_thresholds_sum, best_minareas_sum, max_dices_sum = [0 for x in range(len(dataloaders))], \
                                                            [0 for x in range(len(dataloaders))], [0 for x in range(len(dataloaders))]
    for fold_index, [train_loader, valid_loader] in enumerate(dataloaders):
        if fold_index != 1:
            continue
        # 存放权重的路径+文件名
        load_path = os.path.join(model_path, '%s_fold%d_best.pth' % (config.model_name, fold_index))
        # 加载模型
        model = get_model(config.model_name, load_path)
        mychoose_threshold_minarea = ChooseThresholdMinArea(model, config.model_name, valid_loader, fold_index, model_path)
        best_thresholds, best_minareas, max_dices = mychoose_threshold_minarea.choose_threshold_minarea()
        result = {'best_thresholds': best_thresholds, 'best_minareas': best_minareas, 'max_dices': max_dices}
        results[str(fold_index)] = result

        best_thresholds_sum = [x+y for x,y in zip(best_thresholds_sum, best_thresholds)]
        best_minareas_sum = [x+y for x,y in zip(best_minareas_sum, best_minareas)]
        max_dices_sum = [x+y for x,y in zip(max_dices_sum, max_dices)]
    best_thresholds_average, best_minareas_average, max_dices_average = [x/len(dataloaders) for x in best_thresholds_sum], \
                                                                        [x/len(dataloaders) for x in best_minareas_sum], [x/len(dataloaders) for x in max_dices_sum]
    results['mean'] = {'best_thresholds': best_thresholds_average, 'best_minareas': best_minareas_average, 'max_dices': max_dices_average}
    with codecs.open(model_path + '/%s_result.json' % (config.model_name), 'w', "utf-8") as json_file:
        json.dump(results, json_file, ensure_ascii=False)

    print('save the result')

max_stars_repo_name: zdaiot/Kaggle-Steel-Defect-Detection

max_issues_repo_name: DICLab/Kaggle-Steel-Defect-Detection

max_forks_repo_name: DICLab/Kaggle-Steel-Defect-Detection


text:
#    Copyright 2015-2016 ARM Limited
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#

"""Scheduler specific Functionality for the
stats framework

The Scheduler stats aggregation is based on a signal
which is generated by the combination of two triggers
from the events with the following parameters

========================= ============ =============
EVENT                       VALUE          FILTERS
========================= ============ =============
:func:`sched_switch`        1           next_pid
:func:`sched_switch`      -1           prev_pid
========================= ============ =============

Both these Triggers are provided by the event
:mod:`trappy.sched.SchedSwitch` which correspond to
the :code:`sched_switch` unique word in the trace

.. seealso:: :mod:`trappy.stats.Trigger.Trigger`

Using the above information the following signals are
generated.

**EVENT SERIES**

This is a combination of the two triggers as specified
above and has alternating +/- 1 values and is merely
a representation of the position in time when the process
started or stopped running on a CPU

**RESIDENCY SERIES**

This series is a cumulative sum of the event series and
is a representation of the continuous residency of the
process on a CPU

The pivot for the aggregators is the CPU on which the
event occurred on. If N is the number of CPUs in the
system, N signal for each CPU are generated. These signals
can then be aggregated by specifying a Topology

.. seealso:: :mod:`trappy.stats.Topology.Topology`
"""
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function

import numpy as np
from trappy.stats.Trigger import Trigger

WINDOW_SIZE = 0.0001
"""A control config for filter events. Some analyses
may require ignoring of small interruptions"""

# Trigger Values
SCHED_SWITCH_IN = 1
"""Value of the event when a task is **switch in**
or scheduled on a CPU"""
SCHED_SWITCH_OUT = -1
"""Value of the event when a task is **switched out**
or relinquishes a CPU"""
NO_EVENT = 0
"""Signifies no event on an event trace"""

# Field Names
CPU_FIELD = "__cpu"
"""The column in the sched_switch event that
indicates the CPU on which the event occurred
"""
NEXT_PID_FIELD = "next_pid"
"""The column in the sched_switch event that
indicates the PID of the next process to be scheduled
"""
PREV_PID_FIELD = "prev_pid"
"""The column in the sched_switch event that
indicates the PID of the process that was scheduled
in
"""
TASK_RUNNING = 1
"""The column in the sched_switch event that
indicates the CPU on which the event occurred
"""
TASK_NOT_RUNNING = 0
"""In a residency series, a zero indicates
that the task is not running
"""
TIME_INVAL = -1
"""Standard Value to indicate invalid time data"""
SERIES_SANTIZED = "_sched_sanitized"
"""A memoized flag which is set when an event series
is checked for boundary conditions
"""

def sanitize_asymmetry(series, window=None):
    """Sanitize the cases when a :code:`SWITCH_OUT`
    happens before a :code:`SWITCH_IN`. (The case when
    a process is already running before the trace started)

    :param series: Input Time Series data
    :type series: :mod:`pandas.Series`

    :param window: A tuple indicating a time window
    :type window: tuple
    """

    if not hasattr(series, SERIES_SANTIZED):

        events = series[series != 0]
        if len(series) >= 2 and len(events):
            if series.values[0] == SCHED_SWITCH_OUT:
                series.values[0] = TASK_NOT_RUNNING

            elif events.values[0] == SCHED_SWITCH_OUT:
                series.values[0] = SCHED_SWITCH_IN
                if window:
                    series.index.values[0] = window[0]

            if series.values[-1] == SCHED_SWITCH_IN:
                series.values[-1] = TASK_NOT_RUNNING

            elif events.values[-1] == SCHED_SWITCH_IN:
                series.values[-1] = SCHED_SWITCH_OUT
                if window:
                    series.index.values[-1] = window[1]

        # No point if the series just has one value and
        # one event. We do not have sufficient data points
        # for any calculation. We should Ideally never reach
        # here.
        elif len(series) == 1:
            series.values[0] = 0

        setattr(series, SERIES_SANTIZED, True)

    return series

def csum(series, window=None, filter_gaps=False):
    """:func:`aggfunc` for the cumulative sum of the
    input series data

    :param series: Input Time Series data
    :type series: :mod:`pandas.Series`

    :param window: A tuple indicating a time window
    :type window: tuple

    :param filter_gaps: If set, a process being switched out
        for :mod:`bart.sched.functions.WINDOW_SIZE` is
        ignored. This is helpful when small interruptions need
        to be ignored to compare overall correlation
    :type filter_gaps: bool
    """

    if filter_gaps:
        series = filter_small_gaps(series)

    series = series.cumsum()
    return select_window(series, window)

def filter_small_gaps(series):
    """A helper function that does filtering of gaps
    in residency series < :mod:`bart.sched.functions.WINDOW_SIZE`

    :param series: Input Time Series data
    :type series: :mod:`pandas.Series`
    """

    start = None
    for index, value in series.items():

        if value == SCHED_SWITCH_IN:
            if start == None:
                continue

            if index - start < WINDOW_SIZE:
                series[start] = NO_EVENT
                series[index] = NO_EVENT
            start = None

        if value == SCHED_SWITCH_OUT:
            start = index

    return series

def first_cpu(series, window=None):
    """:func:`aggfunc` to calculate the time of
    the first switch in event in the series
    This is returned as a vector of unit length
    so that it can be aggregated and reduced across
    nodes to find the first cpu of a task

    :param series: Input Time Series data
    :type series: :mod:`pandas.Series`

    :param window: A tuple indicating a time window
    :type window: tuple
    """
    series = select_window(series, window)
    series = series[series == SCHED_SWITCH_IN]
    if len(series):
        return [series.index.values[0]]
    else:
        return [float("inf")]

def last_cpu(series, window=None):
    """:func:`aggfunc` to calculate the time of
    the last switch out event in the series
    This is returned as a vector of unit length
    so that it can be aggregated and reduced across
    nodes to find the last cpu of a task

    :param series: Input Time Series data
    :type series: :mod:`pandas.Series`

    :param window: A tuple indicating a time window
    :type window: tuple
    """
    series = select_window(series, window)
    series = series[series == SCHED_SWITCH_OUT]

    if len(series):
        return [series.index.values[-1]]
    else:
        return [0]

def select_window(series, window):
    """Helper Function to select a portion of
    pandas time series

    :param series: Input Time Series data
    :type series: :mod:`pandas.Series`

    :param window: A tuple indicating a time window
    :type window: tuple
    """

    if not window:
        return series

    start, stop = window
    ix = series.index
    selector = ((ix >= start) & (ix <= stop))
    window_series = series[selector]
    return window_series

def residency_sum(series, window=None):
    """:func:`aggfunc` to calculate the total
    residency

    The input series is processed for
    intervals between a :mod:`bart.sched.functions.SCHED_SWITCH_OUT`
    and :mod:`bart.sched.functions.SCHED_SWITCH_IN` to track
    additive residency of a task

    .. math::

        S_{in} = i_{1}, i_{2}...i_{N} \\\\
        S_{out} = o_{1}, o_{2}...o_{N} \\\\
        R_{total} = \sum_{k}^{N}\Delta_k = \sum_{k}^{N}(o_{k} - i_{k})

    :param series: Input Time Series data
    :type series: :mod:`pandas.Series`

    :param window: A tuple indicating a time window
    :type window: tuple

    :return: A scalar float value
    """

    if not len(series):
        return 0.0

    org_series = series
    series = select_window(series, window)
    series = sanitize_asymmetry(series, window)

    s_in = series[series == SCHED_SWITCH_IN]
    s_out = series[series == SCHED_SWITCH_OUT]

    if not (len(s_in) and len(s_out)):
        try:
            org_series = sanitize_asymmetry(org_series)
            running = select_window(org_series.cumsum(), window)
            if running.values[0] == TASK_RUNNING and running.values[-1] == TASK_RUNNING:
                return window[1] - window[0]
        except Exception as e:
            pass

    if len(s_in) != len(s_out):
        raise RuntimeError(
            "Unexpected Lengths: s_in={}, s_out={}".format(
                len(s_in),
                len(s_out)))
    else:
        return np.sum(s_out.index.values - s_in.index.values)

def first_time(series, value, window=None):
    """:func:`aggfunc` to:

    - Return the first index where the
      series == value

    - If no such index is found
      +inf is returned

    :param series: Input Time Series data
    :type series: :mod:`pandas.Series`

    :param window: A tuple indicating a time window
    :type window: tuple

    :return: A vector of Unit Length
    """

    series = select_window(series, window)
    series = series[series == value]

    if not len(series):
        return [float("inf")]

    return [series.index.values[0]]

def period(series, align="start", window=None):
    """This :func:`aggfunc` returns a tuple
    of the average duration between two triggers:

        - When :code:`align=start` the :code:`SCHED_IN`
          trigger is used

        - When :code:`align=end` the :code:`SCHED_OUT`
          trigger is used

    .. math::

        E = e_{1}, e_{2}...e_{N} \\\\
        T_p = \\frac{\sum_{j}^{\lfloor N/2 \\rfloor}(e_{2j + 1} - e_{2j})}{N}

    :param series: Input Time Series data
    :type series: :mod:`pandas.Series`

    :param window: A tuple indicating a time window
    :type window: tuple

    :return:
        A list of deltas of successive starts/stops
        of a task

    """

    series = select_window(series, window)
    series = sanitize_asymmetry(series, window)

    if align == "start":
        series = series[series == SCHED_SWITCH_IN]
    elif align == "end":
        series = series[series == SCHED_SWITCH_OUT]

    if len(series) % 2 == 0:
        series = series[:1]

    if not len(series):
        return []

    return list(np.diff(series.index.values))

def last_time(series, value, window=None):
    """:func:`aggfunc` to:

    - The first index where the
      series == value

    - If no such index is found
      :mod:`bart.sched.functions.TIME_INVAL`
      is returned

    :param series: Input Time Series data
    :type series: :mod:`pandas.Series`

    :param window: A tuple indicating a time window
    :type window: tuple

    :return: A vector of Unit Length
    """

    series = select_window(series, window)
    series = series[series == value]
    if not len(series):
        return [TIME_INVAL]

    return [series.index.values[-1]]

def binary_correlate(series_x, series_y):
    """Helper function to Correlate binary Data

    Both the series should have same indices

    For binary time series data:

    .. math::

        \\alpha_{corr} = \\frac{N_{agree} - N_{disagree}}{N}

    :param series_x: First time Series data
    :type series_x: :mod:`pandas.Series`

    :param series_y: Second time Series data
    :type series_y: :mod:`pandas.Series`
    """

    if len(series_x) != len(series_y):
        raise ValueError("Cannot compute binary correlation for \
                          unequal vectors")

    agree = len(series_x[series_x == series_y])
    disagree = len(series_x[series_x != series_y])

    return (agree - disagree) / len(series_x)

def get_pids_for_process(ftrace, execname, cls=None):
    """Get the PIDs for a given process

    :param ftrace: A ftrace object with a sched_switch
        event
    :type ftrace: :mod:`trappy.ftrace.FTrace`

    :param execname: The name of the process
    :type execname: str

    :param cls: The SchedSwitch event class (required if
        a different event is to be used)
    :type cls: :mod:`trappy.base.Base`

    :return: The set of PIDs for the execname
    """

    if not cls:
        try:
            df = ftrace.sched_switch.data_frame
        except AttributeError:
            raise ValueError("SchedSwitch event not found in ftrace")

        if len(df) == 0:
            raise ValueError("SchedSwitch event not found in ftrace")
    else:
        event = getattr(ftrace, cls.name)
        df = event.data_frame

    mask = df["next_comm"].apply(lambda x : True if x == execname else False)
    return list(np.unique(df[mask]["next_pid"].values))

def get_task_name(ftrace, pid, cls=None):
    """Returns the execname for pid

    :param ftrace: A ftrace object with a sched_switch
        event
    :type ftrace: :mod:`trappy.ftrace.FTrace`

    :param pid: The PID of the process
    :type pid: int

    :param cls: The SchedSwitch event class (required if
        a different event is to be used)
    :type cls: :mod:`trappy.base.Base`

    :return: The execname for the PID
    """

    if not cls:
        try:
            df = ftrace.sched_switch.data_frame
        except AttributeError:
           raise ValueError("SchedSwitch event not found in ftrace")
    else:
        event = getattr(ftrace, cls.name)
        df = event.data_frame

    df = df[df["next_pid"] == pid]
    if not len(df):
        return ""
    else:
        return df["next_comm"].values[0]

def sched_triggers(ftrace, pid, sched_switch_class):
    """Returns the list of sched_switch triggers

    :param ftrace: A ftrace object with a sched_switch
        event
    :type ftrace: :mod:`trappy.ftrace.FTrace`

    :param pid: The PID of the associated process
    :type pid: int

    :param sched_switch_class: The SchedSwitch event class
    :type sched_switch_class: :mod:`trappy.base.Base`

    :return: List of triggers, such that
        ::

            triggers[0] = switch_in_trigger
            triggers[1] = switch_out_trigger
    """

    if not hasattr(ftrace, "sched_switch"):
        raise ValueError("SchedSwitch event not found in ftrace")

    triggers = []
    triggers.append(sched_switch_in_trigger(ftrace, pid, sched_switch_class))
    triggers.append(sched_switch_out_trigger(ftrace, pid, sched_switch_class))
    return triggers

def sched_switch_in_trigger(ftrace, pid, sched_switch_class):
    """
    :param ftrace: A ftrace object with a sched_switch
        event
    :type ftrace: :mod:`trappy.ftrace.FTrace`

    :param pid: The PID of the associated process
    :type pid: int

    :param sched_switch_class: The SchedSwitch event class
    :type sched_switch_class: :mod:`trappy.base.Base`

    :return: :mod:`trappy.stats.Trigger.Trigger` on
        the SchedSwitch: IN for the given PID
    """

    task_in = {}
    task_in[NEXT_PID_FIELD] = pid

    return Trigger(ftrace,
                   sched_switch_class,              # trappy Event Class
                   task_in,                         # Filter Dictionary
                   SCHED_SWITCH_IN,                 # Trigger Value
                   CPU_FIELD)                       # Primary Pivot

def sched_switch_out_trigger(ftrace, pid, sched_switch_class):
    """
    :param ftrace: A ftrace object with a sched_switch
        event
    :type ftrace: :mod:`trappy.ftrace.FTrace`

    :param pid: The PID of the associated process
    :type pid: int

    :param sched_switch_class: The SchedSwitch event class
    :type sched_switch_class: :mod:`trappy.base.Base`

    :return: :mod:`trappy.stats.Trigger.Trigger` on
        the SchedSwitch: OUT for the given PID
    """

    task_out = {}
    task_out[PREV_PID_FIELD] = pid

    return Trigger(ftrace,
                   sched_switch_class,              # trappy Event Class
                   task_out,                        # Filter Dictionary
                   SCHED_SWITCH_OUT,                # Trigger Value
                   CPU_FIELD)                       # Primary Pivot

def trace_event(series, window=None):
    """
    :func:`aggfunc` to be used for plotting
    the process residency data using
    :mod:`trappy.plotter.EventPlot`

    :param series: Input Time Series data
    :type series: :mod:`pandas.Series`

    :param window: A tuple indicating a time window
    :type window: tuple

    :return: A list of events
        of the type:
        ::

            [
                [start_time_1, stop_time_1],
                [start_time_2, stop_time_2],
                #
                #
                [start_time_N, stop_time_N],
            ]
    """
    rects = []
    series = select_window(series, window)
    series = sanitize_asymmetry(series, window)

    s_in = series[series == SCHED_SWITCH_IN]
    s_out = series[series == SCHED_SWITCH_OUT]

    if not len(s_in):
        return rects

    if len(s_in) != len(s_out):
        raise RuntimeError(
            "Unexpected Lengths: s_in={}, s_out={}".format(
                len(s_in),
                len(s_out)))

    return np.column_stack((s_in.index.values, s_out.index.values))

max_stars_repo_name: JaimeVHArm/lisa

max_issues_repo_name: JaimeVHArm/lisa

max_forks_repo_name: JaimeVHArm/lisa


text:
import numpy as np

def function(x):
    value = 0.
    for i in range(0,ninputs):
        value+=x[i]**2
    return value

nsamples = 4
ninputs = 2
noutputs = 1

fh_in = open('TIN', 'w')
fh_out = open('TOUT', 'w')
for d in range(0,nsamples):
    x = np.random.uniform(low=-1., high=1., size=ninputs)
    #print(x)
    value = function(x)
    #print('value: %f' % (value))
    for i in range(0,ninputs):
        fh_in.write('%f ' % (x[i]))
    fh_in.write('\n')
    fh_out.write('%f\n' % (value))

fh_in.close()
fh_out.close()

max_stars_repo_name: rohskopf/modulant

max_issues_repo_name: rohskopf/modulant

max_forks_repo_name: rohskopf/modulant


text:
from logging import getLogger
import os

import tqdm
import numpy as np
from sklearn.cluster import KMeans
import joblib
import minerl

from IPython import embed
logger = getLogger(__name__)

class _KMeansCacheNotFound(FileNotFoundError):
    pass

def cached_kmeans(cache_dir, env_id, n_clusters, random_state):
    if cache_dir is None:  # ignore cache
        logger.info('Load dataset & do kmeans')
        kmeans = _do_kmeans(env_id=env_id, n_clusters=n_clusters, random_state=random_state)
    else:
        filepath = os.path.join(cache_dir, env_id, f'n_clusters_{n_clusters}', f'random_state_{random_state}', 'kmeans.joblib')
        try:
            kmeans = _load_kmeans_result_cache(filepath)
            logger.info('found kmeans cache')
        except _KMeansCacheNotFound:
            logger.info('kmeans cache not found. Load dataset & do kmeans & save result as cache')
            kmeans = _do_kmeans(env_id=env_id, n_clusters=n_clusters, random_state=random_state)
            _save_kmeans_result_cache(kmeans, filepath)
    return kmeans

def _do_kmeans(env_id, n_clusters, random_state):
    logger.debug(f'loading data...')

    MINERL_DATA_ROOT = os.getenv('MINERL_DATA_ROOT', '/home/usuaris/imatge/juan.jose.nieto/mineRL/data/')
    dat = minerl.data.make(env_id, data_dir=MINERL_DATA_ROOT)
    act_vectors = []
    for _, act, _, _, _ in tqdm.tqdm(dat.batch_iter(batch_size=16, seq_len=32, num_epochs=1, preload_buffer_size=32, seed=random_state)):
        act_vectors.append(act['vector'])
    acts = np.concatenate(act_vectors).reshape(-1, 64)
    logger.debug(f'loading data... done.')
    logger.debug(f'executing keamns...')
    kmeans = KMeans(n_clusters=n_clusters, random_state=random_state).fit(acts)
    logger.debug(f'executing keamns... done.')
    return kmeans

# def _describe_kmeans_result(kmeans):
#     result = [(obf_a, minerl.herobraine.envs.MINERL_TREECHOP_OBF_V0.unwrap_action({'vector': obf_a})) for obf_a in kmeans.cluster_centers_]
#     logger.debug(result)
#     return result

def _save_kmeans_result_cache(kmeans, filepath):
    os.makedirs(os.path.dirname(filepath), exist_ok=True)
    joblib.dump(kmeans, filepath)
    logger.info(f'saved kmeans {filepath}')

def _load_kmeans_result_cache(filepath):
    if not os.path.exists(filepath):
        raise _KMeansCacheNotFound
    logger.debug(f'loading kmeans {filepath}')
    return joblib.load(filepath)

max_stars_repo_name: yuyecreus/Habitat-PixelEDL

max_issues_repo_name: imatge-upc/pixelcoordEDL

max_forks_repo_name: imatge-upc/pixelcoordEDL