text:
"""
Robust linear models with support for the M-estimators  listed under
:ref:`norms <norms>`.

References
----------
PJ Huber.  'Robust Statistics' John Wiley and Sons, Inc., New York.  1981.

PJ Huber.  1973,  'The 1972 Wald Memorial Lectures: Robust Regression:
    Asymptotics, Conjectures, and Monte Carlo.'  The Annals of Statistics,
    1.5, 799-821.

R Venables, B Ripley. 'Modern Applied Statistics in S'  Springer, New York,
    2002.
"""
import numpy as np
import scipy.stats as stats

from statsmodels.tools.decorators import cache_readonly
from statsmodels.tools.sm_exceptions import ConvergenceWarning
import statsmodels.regression.linear_model as lm
import statsmodels.regression._tools as reg_tools
import statsmodels.robust.norms as norms
import statsmodels.robust.scale as scale
import statsmodels.base.model as base
import statsmodels.base.wrapper as wrap

__all__ = ['RLM']

def _check_convergence(criterion, iteration, tol, maxiter):
    cond = np.abs(criterion[iteration] - criterion[iteration - 1])
    return not (np.any(cond > tol) and iteration < maxiter)

class RLM(base.LikelihoodModel):
    __doc__ = """
    Robust Linear Model

    Estimate a robust linear model via iteratively reweighted least squares
    given a robust criterion estimator.

    %(params)s
    M : statsmodels.robust.norms.RobustNorm, optional
        The robust criterion function for downweighting outliers.
        The current options are LeastSquares, HuberT, RamsayE, AndrewWave,
        TrimmedMean, Hampel, and TukeyBiweight.  The default is HuberT().
        See statsmodels.robust.norms for more information.
    %(extra_params)s

    Attributes
    ----------

    df_model : float
        The degrees of freedom of the model.  The number of regressors p less
        one for the intercept.  Note that the reported model degrees
        of freedom does not count the intercept as a regressor, though
        the model is assumed to have an intercept.
    df_resid : float
        The residual degrees of freedom.  The number of observations n
        less the number of regressors p.  Note that here p does include
        the intercept as using a degree of freedom.
    endog : ndarray
        See above.  Note that endog is a reference to the data so that if
        data is already an array and it is changed, then `endog` changes
        as well.
    exog : ndarray
        See above.  Note that endog is a reference to the data so that if
        data is already an array and it is changed, then `endog` changes
        as well.
    M : statsmodels.robust.norms.RobustNorm
         See above.  Robust estimator instance instantiated.
    nobs : float
        The number of observations n
    pinv_wexog : ndarray
        The pseudoinverse of the design / exogenous data array.  Note that
        RLM has no whiten method, so this is just the pseudo inverse of the
        design.
    normalized_cov_params : ndarray
        The p x p normalized covariance of the design / exogenous data.
        This is approximately equal to (X.T X)^(-1)

    Examples
    --------
    >>> import statsmodels.api as sm
    >>> data = sm.datasets.stackloss.load(as_pandas=False)
    >>> data.exog = sm.add_constant(data.exog)
    >>> rlm_model = sm.RLM(data.endog, data.exog, \
                           M=sm.robust.norms.HuberT())

    >>> rlm_results = rlm_model.fit()
    >>> rlm_results.params
    array([  0.82938433,   0.92606597,  -0.12784672, -41.02649835])
    >>> rlm_results.bse
    array([ 0.11100521,  0.30293016,  0.12864961,  9.79189854])
    >>> rlm_results_HC2 = rlm_model.fit(cov="H2")
    >>> rlm_results_HC2.params
    array([  0.82938433,   0.92606597,  -0.12784672, -41.02649835])
    >>> rlm_results_HC2.bse
    array([ 0.11945975,  0.32235497,  0.11796313,  9.08950419])
    >>> mod = sm.RLM(data.endog, data.exog, M=sm.robust.norms.Hampel())
    >>> rlm_hamp_hub = mod.fit(scale_est=sm.robust.scale.HuberScale())
    >>> rlm_hamp_hub.params
    array([  0.73175452,   1.25082038,  -0.14794399, -40.27122257])
    """ % {'params': base._model_params_doc,
           'extra_params': base._missing_param_doc}

    def __init__(self, endog, exog, M=None, missing='none',
                 **kwargs):
        self.M = M if M is not None else norms.HuberT()
        super(base.LikelihoodModel, self).__init__(endog, exog,
                                                   missing=missing, **kwargs)
        self._initialize()
        # things to remove_data
        self._data_attr.extend(['weights', 'pinv_wexog'])

    def _initialize(self):
        """
        Initializes the model for the IRLS fit.

        Resets the history and number of iterations.
        """
        self.pinv_wexog = np.linalg.pinv(self.exog)
        self.normalized_cov_params = np.dot(self.pinv_wexog,
                                            np.transpose(self.pinv_wexog))
        self.df_resid = (float(self.exog.shape[0] -
                               np.linalg.matrix_rank(self.exog)))
        self.df_model = float(np.linalg.matrix_rank(self.exog) - 1)
        self.nobs = float(self.endog.shape[0])

    def score(self, params):
        raise NotImplementedError

    def information(self, params):
        raise NotImplementedError

    def predict(self, params, exog=None):
        """
        Return linear predicted values from a design matrix.

        Parameters
        ----------
        params : array_like
            Parameters of a linear model
        exog : array_like, optional.
            Design / exogenous data. Model exog is used if None.

        Returns
        -------
        An array of fitted values
        """
        # copied from linear_model  # TODO: then is it needed?
        if exog is None:
            exog = self.exog
        return np.dot(exog, params)

    def loglike(self, params):
        raise NotImplementedError

    def deviance(self, tmp_results):
        """
        Returns the (unnormalized) log-likelihood from the M estimator.
        """
        tmp_resid = self.endog - tmp_results.fittedvalues
        return self.M(tmp_resid / tmp_results.scale).sum()

    def _update_history(self, tmp_results, history, conv):
        history['params'].append(tmp_results.params)
        history['scale'].append(tmp_results.scale)
        if conv == 'dev':
            history['deviance'].append(self.deviance(tmp_results))
        elif conv == 'sresid':
            history['sresid'].append(tmp_results.resid / tmp_results.scale)
        elif conv == 'weights':
            history['weights'].append(tmp_results.model.weights)
        return history

    def _estimate_scale(self, resid):
        """
        Estimates the scale based on the option provided to the fit method.
        """
        if isinstance(self.scale_est, str):
            if self.scale_est.lower() == 'mad':
                return scale.mad(resid, center=0)
            else:
                raise ValueError("Option %s for scale_est not understood" %
                                 self.scale_est)
        elif isinstance(self.scale_est, scale.HuberScale):
            return self.scale_est(self.df_resid, self.nobs, resid)
        else:
            return scale.scale_est(self, resid) ** 2

    def fit(self, maxiter=50, tol=1e-8, scale_est='mad', init=None, cov='H1',
            update_scale=True, conv='dev', start_params=None):
        """
        Fits the model using iteratively reweighted least squares.

        The IRLS routine runs until the specified objective converges to `tol`
        or `maxiter` has been reached.

        Parameters
        ----------
        conv : str
            Indicates the convergence criteria.
            Available options are "coefs" (the coefficients), "weights" (the
            weights in the iteration), "sresid" (the standardized residuals),
            and "dev" (the un-normalized log-likelihood for the M
            estimator).  The default is "dev".
        cov : str, optional
            'H1', 'H2', or 'H3'
            Indicates how the covariance matrix is estimated.  Default is 'H1'.
            See rlm.RLMResults for more information.
        init : str
            Specifies method for the initial estimates of the parameters.
            Default is None, which means that the least squares estimate
            is used.  Currently it is the only available choice.
        maxiter : int
            The maximum number of iterations to try. Default is 50.
        scale_est : str or HuberScale()
            'mad' or HuberScale()
            Indicates the estimate to use for scaling the weights in the IRLS.
            The default is 'mad' (median absolute deviation.  Other options are
            'HuberScale' for Huber's proposal 2. Huber's proposal 2 has
            optional keyword arguments d, tol, and maxiter for specifying the
            tuning constant, the convergence tolerance, and the maximum number
            of iterations. See statsmodels.robust.scale for more information.
        tol : float
            The convergence tolerance of the estimate.  Default is 1e-8.
        update_scale : Bool
            If `update_scale` is False then the scale estimate for the
            weights is held constant over the iteration.  Otherwise, it
            is updated for each fit in the iteration.  Default is True.
        start_params : array_like, optional
            Initial guess of the solution of the optimizer. If not provided,
            the initial parameters are computed using OLS.

        Returns
        -------
        results : statsmodels.rlm.RLMresults
            Results instance
        """
        if cov.upper() not in ["H1", "H2", "H3"]:
            raise ValueError("Covariance matrix %s not understood" % cov)
        else:
            self.cov = cov.upper()
        conv = conv.lower()
        if conv not in ["weights", "coefs", "dev", "sresid"]:
            raise ValueError("Convergence argument %s not understood" % conv)
        self.scale_est = scale_est

        if start_params is None:
            wls_results = lm.WLS(self.endog, self.exog).fit()
        else:
            start_params = np.asarray(start_params, dtype=np.double).squeeze()
            if (start_params.shape[0] != self.exog.shape[1] or
                    start_params.ndim != 1):
                raise ValueError('start_params must by a 1-d array with {0} '
                                 'values'.format(self.exog.shape[1]))
            fake_wls = reg_tools._MinimalWLS(self.endog, self.exog,
                                             weights=np.ones_like(self.endog),
                                             check_weights=False)
            wls_results = fake_wls.results(start_params)

        if not init:
            self.scale = self._estimate_scale(wls_results.resid)

        history = dict(params=[np.inf], scale=[])
        if conv == 'coefs':
            criterion = history['params']
        elif conv == 'dev':
            history.update(dict(deviance=[np.inf]))
            criterion = history['deviance']
        elif conv == 'sresid':
            history.update(dict(sresid=[np.inf]))
            criterion = history['sresid']
        elif conv == 'weights':
            history.update(dict(weights=[np.inf]))
            criterion = history['weights']

        # done one iteration so update
        history = self._update_history(wls_results, history, conv)
        iteration = 1
        converged = 0
        while not converged:
            if self.scale == 0.0:
                import warnings
                warnings.warn('Estimated scale is 0.0 indicating that the most'
                              ' last iteration produced a perfect fit of the '
                              'weighted data.', ConvergenceWarning)
                break
            self.weights = self.M.weights(wls_results.resid / self.scale)
            wls_results = reg_tools._MinimalWLS(self.endog, self.exog,
                                                weights=self.weights,
                                                check_weights=True).fit()
            if update_scale is True:
                self.scale = self._estimate_scale(wls_results.resid)
            history = self._update_history(wls_results, history, conv)
            iteration += 1
            converged = _check_convergence(criterion, iteration, tol, maxiter)
        results = RLMResults(self, wls_results.params,
                             self.normalized_cov_params, self.scale)

        history['iteration'] = iteration
        results.fit_history = history
        results.fit_options = dict(cov=cov.upper(), scale_est=scale_est,
                                   norm=self.M.__class__.__name__, conv=conv)
        # norm is not changed in fit, no old state

        # doing the next causes exception
        # self.cov = self.scale_est = None #reset for additional fits
        # iteration and history could contain wrong state with repeated fit
        return RLMResultsWrapper(results)

class RLMResults(base.LikelihoodModelResults):
    """
    Class to contain RLM results

    Attributes
    ----------

    bcov_scaled : ndarray
        p x p scaled covariance matrix specified in the model fit method.
        The default is H1. H1 is defined as
        ``k**2 * (1/df_resid*sum(M.psi(sresid)**2)*scale**2)/
        ((1/nobs*sum(M.psi_deriv(sresid)))**2) * (X.T X)^(-1)``

        where ``k = 1 + (df_model +1)/nobs * var_psiprime/m**2``
        where ``m = mean(M.psi_deriv(sresid))`` and
        ``var_psiprime = var(M.psi_deriv(sresid))``

        H2 is defined as
        ``k * (1/df_resid) * sum(M.psi(sresid)**2) *scale**2/
        ((1/nobs)*sum(M.psi_deriv(sresid)))*W_inv``

        H3 is defined as
        ``1/k * (1/df_resid * sum(M.psi(sresid)**2)*scale**2 *
        (W_inv X.T X W_inv))``

        where `k` is defined as above and
        ``W_inv = (M.psi_deriv(sresid) exog.T exog)^(-1)``

        See the technical documentation for cleaner formulae.
    bcov_unscaled : ndarray
        The usual p x p covariance matrix with scale set equal to 1.  It
        is then just equivalent to normalized_cov_params.
    bse : ndarray
        An array of the standard errors of the parameters.  The standard
        errors are taken from the robust covariance matrix specified in the
        argument to fit.
    chisq : ndarray
        An array of the chi-squared values of the parameter estimates.
    df_model
        See RLM.df_model
    df_resid
        See RLM.df_resid
    fit_history : dict
        Contains information about the iterations. Its keys are `deviance`,
        `params`, `iteration` and the convergence criteria specified in
        `RLM.fit`, if different from `deviance` or `params`.
    fit_options : dict
        Contains the options given to fit.
    fittedvalues : ndarray
        The linear predicted values.  dot(exog, params)
    model : statsmodels.rlm.RLM
        A reference to the model instance
    nobs : float
        The number of observations n
    normalized_cov_params : ndarray
        See RLM.normalized_cov_params
    params : ndarray
        The coefficients of the fitted model
    pinv_wexog : ndarray
        See RLM.pinv_wexog
    pvalues : ndarray
        The p values associated with `tvalues`. Note that `tvalues` are assumed
        to be distributed standard normal rather than Student's t.
    resid : ndarray
        The residuals of the fitted model.  endog - fittedvalues
    scale : float
        The type of scale is determined in the arguments to the fit method in
        RLM.  The reported scale is taken from the residuals of the weighted
        least squares in the last IRLS iteration if update_scale is True.  If
        update_scale is False, then it is the scale given by the first OLS
        fit before the IRLS iterations.
    sresid : ndarray
        The scaled residuals.
    tvalues : ndarray
        The "t-statistics" of params. These are defined as params/bse where
        bse are taken from the robust covariance matrix specified in the
        argument to fit.
    weights : ndarray
        The reported weights are determined by passing the scaled residuals
        from the last weighted least squares fit in the IRLS algorithm.

    See Also
    --------
    statsmodels.base.model.LikelihoodModelResults
    """

    def __init__(self, model, params, normalized_cov_params, scale):
        super(RLMResults, self).__init__(model, params,
                                         normalized_cov_params, scale)
        self.model = model
        self.df_model = model.df_model
        self.df_resid = model.df_resid
        self.nobs = model.nobs
        self._cache = {}
        # for remove_data
        self._data_in_cache.extend(['sresid'])

        self.cov_params_default = self.bcov_scaled
        # TODO: "pvals" should come from chisq on bse?

    @cache_readonly
    def fittedvalues(self):
        return np.dot(self.model.exog, self.params)

    @cache_readonly
    def resid(self):
        return self.model.endog - self.fittedvalues  # before bcov

    @cache_readonly
    def sresid(self):
        if self.scale == 0.0:
            sresid = self.resid.copy()
            sresid[:] = 0.0
            return sresid
        return self.resid / self.scale

    @cache_readonly
    def bcov_unscaled(self):
        return self.normalized_cov_params

    @cache_readonly
    def weights(self):
        return self.model.weights

    @cache_readonly
    def bcov_scaled(self):
        model = self.model
        m = np.mean(model.M.psi_deriv(self.sresid))
        var_psiprime = np.var(model.M.psi_deriv(self.sresid))
        k = 1 + (self.df_model + 1) / self.nobs * var_psiprime / m ** 2

        if model.cov == "H1":
            ss_psi = np.sum(model.M.psi(self.sresid) ** 2)
            s_psi_deriv = np.sum(model.M.psi_deriv(self.sresid))
            return k ** 2 * (1 / self.df_resid * ss_psi * self.scale ** 2) /\
                ((1 / self.nobs * s_psi_deriv) ** 2) *\
                model.normalized_cov_params
        else:
            W = np.dot(model.M.psi_deriv(self.sresid) * model.exog.T,
                       model.exog)
            W_inv = np.linalg.inv(W)
            # [W_jk]^-1 = [SUM(psi_deriv(Sr_i)*x_ij*x_jk)]^-1
            # where Sr are the standardized residuals
            if model.cov == "H2":
                # These are correct, based on Huber (1973) 8.13
                return k * (1 / self.df_resid) * np.sum(
                    model.M.psi(self.sresid) ** 2) * self.scale ** 2 \
                       / ((1 / self.nobs) *
                          np.sum(model.M.psi_deriv(self.sresid))) * W_inv
            elif model.cov == "H3":
                return k ** -1 * 1 / self.df_resid * np.sum(
                    model.M.psi(self.sresid) ** 2) * self.scale ** 2 \
                       * np.dot(
                    np.dot(W_inv, np.dot(model.exog.T, model.exog)),
                    W_inv)

    @cache_readonly
    def pvalues(self):
        return stats.norm.sf(np.abs(self.tvalues)) * 2

    @cache_readonly
    def bse(self):
        return np.sqrt(np.diag(self.bcov_scaled))

    @cache_readonly
    def chisq(self):
        return (self.params / self.bse) ** 2

    def summary(self, yname=None, xname=None, title=0, alpha=.05,
                return_fmt='text'):
        """
        This is for testing the new summary setup
        """
        top_left = [('Dep. Variable:', None),
                    ('Model:', None),
                    ('Method:', ['IRLS']),
                    ('Norm:', [self.fit_options['norm']]),
                    ('Scale Est.:', [self.fit_options['scale_est']]),
                    ('Cov Type:', [self.fit_options['cov']]),
                    ('Date:', None),
                    ('Time:', None),
                    ('No. Iterations:', ["%d" % self.fit_history['iteration']])
                    ]
        top_right = [('No. Observations:', None),
                     ('Df Residuals:', None),
                     ('Df Model:', None)
                     ]

        if title is not None:
            title = "Robust linear Model Regression Results"

        # boiler plate
        from statsmodels.iolib.summary import Summary
        smry = Summary()
        smry.add_table_2cols(self, gleft=top_left, gright=top_right,
                             yname=yname, xname=xname, title=title)
        smry.add_table_params(self, yname=yname, xname=xname, alpha=alpha,
                              use_t=self.use_t)

        # add warnings/notes, added to text format only
        etext = []
        wstr = ("If the model instance has been used for another fit with "
                "different fit parameters, then the fit options might not be "
                "the correct ones anymore .")
        etext.append(wstr)

        if etext:
            smry.add_extra_txt(etext)

        return smry

    def summary2(self, xname=None, yname=None, title=None, alpha=.05,
                 float_format="%.4f"):
        """Experimental summary function for regression results

        Parameters
        ----------
        yname : str
            Name of the dependent variable (optional)
        xname : list[str], optional
            Names for the exogenous variables. Default is `var_##` for ## in
            the number of regressors. Must match the number of parameters
            in the model
        title : str, optional
            Title for the top table. If not None, then this replaces the
            default title
        alpha : float
            significance level for the confidence intervals
        float_format : str
            print format for floats in parameters summary

        Returns
        -------
        smry : Summary instance
            this holds the summary tables and text, which can be printed or
            converted to various output formats.

        See Also
        --------
        statsmodels.iolib.summary2.Summary : class to hold summary results
        """
        from statsmodels.iolib import summary2
        smry = summary2.Summary()
        smry.add_base(results=self, alpha=alpha, float_format=float_format,
                      xname=xname, yname=yname, title=title)

        return smry

class RLMResultsWrapper(lm.RegressionResultsWrapper):
    pass

wrap.populate_wrapper(RLMResultsWrapper, RLMResults)  # noqa:E305

max_stars_repo_name: EkremBayar/bayar

max_issues_repo_name: EkremBayar/bayar

max_forks_repo_name: EkremBayar/bayar

reason: import numpy,import scipy,import statsmodels,from statsmodels


text:
"""
Variables
---------
Module description
"""
from abc import ABC, abstractmethod
import operator
import numbers
import collections
from collections.abc import Iterable

import chainer
import chainer.links as L
import chainer.functions as F
import numpy as np

from brancher.utilities import join_dicts_list, join_sets_list
from brancher.utilities import flatten_list
from brancher.utilities import partial_broadcast
from brancher.utilities import coerce_to_dtype
from brancher.utilities import broadcast_parent_values
from brancher.utilities import split_dict
from brancher.utilities import reformat_sampler_input
from brancher.utilities import tile_parameter
from brancher.utilities import get_model_mapping
from brancher.utilities import reassign_samples

from brancher.pandas_interface import reformat_sample_to_pandas
from brancher.pandas_interface import reformat_model_summary
from brancher.pandas_interface import pandas_frame2dict
from brancher.pandas_interface import pandas_frame2value

class BrancherClass(ABC):
    """
    BrancherClass is the abstract superclass of all Brancher variables and models.
    """
    @abstractmethod
    def _flatten(self):
        """
        Abstract method. It returs a list of all the variables contained in the model.
        """
        pass

    def flatten(self):
        return set(self._flatten())

    def get_variable(self, var_name):
        """
        It returns the variable in the model with the requested name.

        Args:
            var_name: String. Name  of the requested variable.

        Returns:
            brancher.Variable.

        """
        flat_list = self._flatten()
        try:
            return {var.name: var for var in flat_list}[var_name]
        except ValueError:
            raise ValueError("The variable {} is not present in the model".format(var_name))

class Variable(BrancherClass):
    """
    Variable is the abstract superclass of deterministic and random variables. Variables are the building blocks of
    all probabilistic models in Brancher.
    """
    @abstractmethod
    def calculate_log_probability(self, values, reevaluate):
        """
        Abstract method. It returns the log probability of the values given the model.

        Args:
            values: Dictionary(brancher.Variable: chainer.Variable). A dictionary having the brancher.variables of the
            model as keys and chainer.Variables as values. This dictionary has to provide values for all variables of
            the model except for the deterministic variables.

            reevaluate: Bool. If false it returns the output of the latest call. It avoid unnecessary computations when
            multiple children variables ask for the log probability of the same paternt variable.

        Returns:
            chainer.Variable. the log probability of the input values given the model.

        """
        pass

    @abstractmethod
    def _get_sample(self, number_samples, resample, observed, input_values):
        """
        Abstract method. It returns samples from the joint distribution specified by the model. If an input is provided
        it only samples the variables that are not contained in the input.

        Args:
            number_samples: Int.

            resample: Bool. If true it returns its previously stored sampled. It is used when multiple children variables
            ask for a sample to the same parent. In this case the resample variable is False since the children should be
            fed with the same value of the parent.

            observed: Bool. It specifies whether the samples should be interpreted frequentistically as samples from the
            observations of as Bayesian samples from the prior model. The first batch dimension is reserved to Bayesian
            samples while the second batch dimension is reserved to observation samples. Theerefore, this boolean
            changes the shape of the resulting sampled array.

            input_values: Dictionary(brancher.Variable, chainer.Variable). A dictionary having the brancher.variables of the
            model as keys and chainer.Variables as values. This dictionary has to provide values for all variables of
            the model that do not need to be sampled. Using an input allows to use a probabilistic model as a random
            function.

        Returns:
            Dictionary(brancher.Variable: chainer.Variable). A dictionary of samples from all the variables of the model

        """
        pass

    def get_sample(self, number_samples, input_values={}):
        reformatted_input_values = reformat_sampler_input(pandas_frame2dict(input_values),
                                                          number_samples=number_samples)
        raw_sample = {self: self._get_sample(number_samples, resample=False,
                                             observed=self.is_observed, input_values=reformatted_input_values)[self]}
        sample = reformat_sample_to_pandas(raw_sample, number_samples)
        self.reset()
        return sample

    @abstractmethod
    def reset(self):
        """
        Abstract method. It recursively reset the self._evalueted and self._current_value attributes of the variable and
        all downstream variables. It is used after sampling and evaluating the log probability of a model.

        Args: None.

        Returns: None.
        """
        pass

    @property
    @abstractmethod
    def is_observed(self):
        """
        Abstract property method. It returns True if the variable is observed and False otherwise.

        Args: None.

        Returns: Bool.
        """
        pass

    def __str__(self):
        """
        Method.

        Args: None

        Returns: String
        """
        return self.name

    def _apply_operator(self, other, op):
        """
        Method. It is used for using operations between variables symbolically. It always returns a partialLink object
        that define a mathematical operation between variables. The vars attribute of the link is the set of variables
        that are used in the operation. The fn attribute is a lambda that specify the operation as a functions between the
        values of the variables in vars and a numeric output. This is required for defining the forward pass of the model.

        Args:
            other: PartialLink, RandomVariable, numeric or np.array.

            op: Binary operator.

        Returns: PartialLink
        """
        if isinstance(other, PartialLink):
            vars = other.vars
            vars.add(self)
            fn = lambda values: op(values[self], other.fn(values))
            links = other.links
        elif isinstance(other, Variable):
            vars = {self, other}
            fn = lambda values: op(values[self], values[other])
            links = set()
        elif isinstance(other, (numbers.Number, np.ndarray)):
            vars = {self}
            fn = lambda values: op(values[self], other)
            links = set()
        else:
            raise TypeError('') #TODO

        return PartialLink(vars=vars, fn=fn, links=links)

    def __add__(self, other):
        return self._apply_operator(other, operator.add)

    def __radd__(self, other):
        return self.__add__(other)

    def __sub__(self, other):
        return self._apply_operator(other, operator.sub)

    def __rsub__(self, other):
        return -1*self.__sub__(other)

    def __mul__(self, other):
        return self._apply_operator(other, operator.mul)

    def __rmul__(self, other):
        return self.__mul__(other)

    def __truediv__(self, other):
        return self._apply_operator(other, operator.truediv)

    def __rtruediv__(self, other):
        raise NotImplementedError

    def __pow__(self, other):
        return self._apply_operator(other, operator.pow)

    def __rpow__(self, other):
        raise NotImplementedError

    def __getitem__(self, key):
        if isinstance(key, collections.Iterable):
            variable_slice = (slice(None, None, None), *key)
        else:
            variable_slice = (slice(None, None, None), key)
        vars = {self}
        fn = lambda values: values[self][variable_slice]
        links = set()
        return PartialLink(vars=vars, fn=fn, links=links)

    def shape(self):
        vars = {self}
        fn = lambda values: values[self].shape
        links = set()
        return PartialLink(vars=vars, fn=fn, links=links)

class DeterministicVariable(Variable):
    """
    Deterministic variables are a subclass of random variables that always return the same value. The hyper-parameters of
    a probabilistic model are usually encoded as DeterministicVariables. When the user input a parameter as a Numeric value or
    an array, Brancher created a DeterministicVariable that store its value.

    Parameters
    ----------
    data : chainer.Variable, numeric, or np.ndarray. The value of the variable. It gets stored in the self.current value
    attribute.

    name : String. The name of the variable.

    learnable : Bool. This boolean value specify if the value of the DeterministicVariable can be updated during treaning.

    """
    def __init__(self, data, name, learnable=False, is_observed=False):
        self._current_value = coerce_to_dtype(data, is_observed)
        self.name = name
        self._observed = is_observed
        self.parents = ()
        self._type = "Deterministic"
        self.learnable = learnable
        if learnable:
            if not isinstance(self._current_value, collections.abc.Iterable):
                self.link = L.Bias(axis=1, shape=self._current_value.shape[1:]) #TODO: For Julia: this can be implemented as parameter
            else:
                self.learnable = False #TODO: Warning?

    def calculate_log_probability(self, values, reevaluate=True, for_gradient=False, normalized=True):
        """
        Method. It returns the log probability of the values given the model. This value is always 0 since the probability
        of a deterministic variable having its value is always 1.

        Args:
            values: Dictionary(brancher.Variable: chainer.Variable). A dictionary having the brancher.variables of the
            model as keys and chainer.Variables as values. This dictionary has to provide values for all variables of
            the model except for the deterministic variables.

            reevaluate: Bool. If false it returns the output of the latest call. It avoid unnecessary computations when
            multiple children variables ask for the log probability of the same paternt variable.

        Returns:
            chainer.Variable. the log probability of the input values given the model.
        """
        return 0.

    @property
    def value(self):
        assert self._current_value is not None
        if self.learnable:
            return self.link(self._current_value) #TODO: For Julia: this can be implemented as parameter
        return self._current_value

    @value.setter
    def value(self, val):
        self._current_value = coerce_to_dtype(val, is_observed=self.is_observed)

    @property
    def is_observed(self):
        return self._observed

    def _get_sample(self, number_samples, resample=False, observed=False, input_values={}):
        if self in input_values:
            value = input_values[self]
        else:
            value = self.value
        if isinstance(value, chainer.Variable):
            return {self: tile_parameter(value, number_samples=number_samples)}
        else:
            return {self: value} #TODO: This is for allowing discrete data, temporary? (for Julia)

    def reset(self):
        pass

    def _flatten(self):
        return [self]

class RandomVariable(Variable):
    """
    Random variables are the main building blocks of probabilistic models.

    Parameters
    ----------
    distribution : brancher.Distribution
        Summary
    name : str
        Summary
    parents : tuple of brancher variables
        Summary
    link : callable
        Summary
    """
    def __init__(self, distribution, name, parents, link):
        self.name = name
        self.distribution = distribution
        self.link = link
        self.parents = parents
        self._type = "Random"
        self.samples = []

        self._evaluated = False
        self._observed = False
        self._observed_value = None
        self._current_value = None
        self.dataset = None
        self.has_random_dataset = False
        self.has_observed_value = False

    @property
    def value(self):
        if self._observed:
            return self._observed_value
        return self._current_value

    @value.setter
    def value(self, val):
        self._current_value = coerce_to_dtype(val)

    @property
    def is_observed(self):
        return self._observed

    def _apply_link(self, parents_values):  #TODO: This is for allowing discrete data, temporary? (for julia) #For Julia: Very important method
        cont_values, discrete_values = split_dict(parents_values,
                                                  condition=lambda key, val: isinstance(val, chainer.Variable))
        if cont_values:
            reshaped_dict, number_samples, number_datapoints = broadcast_parent_values(cont_values)
            reshaped_dict.update(discrete_values)
        else:
            reshaped_dict = discrete_values
        reshaped_output = self.link(reshaped_dict)
        output = {key: F.reshape(val, (number_samples, number_datapoints) + val.shape[1:])
                  if isinstance(val, chainer.Variable) else val
                  for key, val in reshaped_output.items()}
        return output

    def calculate_log_probability(self, input_values, reevaluate=True, for_gradient=False,
                                  include_parents=True, normalized=True):
        """
        Method. It returns the log probability of the values given the model. This value is always 0 since the probability
        of a deterministic variable having its value is always 1.

        Args:
            values: Dictionary(brancher.Variable: chainer.Variable). A dictionary having the brancher.variables of the
            model as keys and chainer.Variables as values. This dictionary has to provide values for all variables of
            the model except for the deterministic variables.

            reevaluate: Bool. If false it returns the output of the latest call. It avoid unnecessary computations when
            multiple children variables ask for the log probability of the same paternt variable.

        Returns:
            chainer.Variable. the log probability of the input values given the model.

        """
        if self._evaluated and not reevaluate:
            return 0.
        if self in input_values:
            value = input_values[self]
        else:
            value = self.value

        self._evaluated = True
        deterministic_parents_values = {parent: parent.value for parent in self.parents
                                        if (type(parent) is DeterministicVariable)}
        parents_input_values = {parent: parent_input for parent, parent_input in input_values.items() if parent in self.parents}
        parents_values = {**parents_input_values, **deterministic_parents_values}
        parameters_dict = self._apply_link(parents_values)
        log_probability = self.distribution.calculate_log_probability(value, **parameters_dict)
        parents_log_probability = sum([parent.calculate_log_probability(input_values, reevaluate, for_gradient,
                                                                        normalized=normalized)
                                       for parent in self.parents])
        if self.is_observed:
            log_probability = F.sum(log_probability, axis=1, keepdims=True)
        if type(log_probability) is chainer.Variable and type(parents_log_probability) is chainer.Variable:
            log_probability, parents_log_probability = partial_broadcast(log_probability, parents_log_probability)
        if include_parents:
            return log_probability + parents_log_probability
        else:
            return log_probability

    def _get_sample(self, number_samples=1, resample=True, observed=False, input_values={}):
        """
        Summary
        """
        if self.samples and not resample:
            return {self: self.samples[-1]}
        if not observed:
            if self in input_values:
                return {self: input_values[self]}
            else:
                var_to_sample = self
        else:
            if self.has_observed_value:
                return {self: self._observed_value}
            elif self.has_random_dataset:
                var_to_sample = self.dataset
            else:
                var_to_sample = self
        parents_samples_dict = join_dicts_list([parent._get_sample(number_samples, resample, observed, input_values)
                                                for parent in var_to_sample.parents])
        input_dict = {parent: parents_samples_dict[parent] for parent in var_to_sample.parents}
        parameters_dict = var_to_sample._apply_link(input_dict)
        sample = var_to_sample.distribution.get_sample(**parameters_dict, number_samples=number_samples)
        self.samples = [sample] #TODO: to fix
        return {**parents_samples_dict, self: sample}

    def observe(self, data, random_indices=()):
        """
        Summary
        """
        data = pandas_frame2value(data, self.name)
        if isinstance(data, RandomVariable):
            self.dataset = data
            self.has_random_dataset = True
        else:
            self._observed_value = coerce_to_dtype(data, is_observed=True)
            self.has_observed_value = True
        self._observed = True

    def unobserve(self):
        self._observed = False
        self.has_observed_value = False
        self.has_random_dataset = False
        self._observed_value = None
        self.dataset = None

    def reset(self):
        """
        Summary
        """
        self.samples = []
        self._evaluated = False
        self._current_value = None
        for parent in self.parents:
            parent.reset()

    def _flatten(self):
        return flatten_list([parent._flatten() for parent in self.parents]) + [self]

class ProbabilisticModel(BrancherClass):
    """
    Summary

    Parameters
    ----------
    variables : tuple of brancher variables
        Summary
    """
    def __init__(self, variables):
        self.variables = self._validate_variables(variables)
        self._set_summary()
        self.posterior_model = None
        self.posterior_sampler = None
        self.observed_submodel = None
        self.diagnostics = {}
        if not all([var.is_observed for var in self.variables]): #TODO: this is not elegant
            self.update_observed_submodel()
        else:
            self.observed_submodel = self

    def __str__(self): #TODO: Work in progress
        """
        Method.

        Args: None

        Returns: String
        """
        return self.model_summary.__str__()

    @staticmethod
    def _validate_variables(variables):
        for var in variables:
            if not isinstance(var, (DeterministicVariable, RandomVariable, ProbabilisticModel)):
                raise ValueError("Invalid input type: {}".format(type(var)))
        return variables

    def _set_summary(self): #TODO: Work in progress
        feature_list = ["Distribution", "Parents", "Observed"]
        var_list = self.flatten()
        var_names = [var.name for var in var_list]
        summary_data = [[var._type, var.parents, var.is_observed]
                         for var in var_list]
        self._model_summary = reformat_model_summary(summary_data, var_names, feature_list)

    @property
    def model_summary(self):
        self._set_summary()
        return self._model_summary

    @property
    def value(self):
        return tuple(var.value for var in self.variables)

    @value.setter
    def value(self, val):
        raise AttributeError("The value of a probabilistic model cannot be explicitly set.")

    @property
    def is_observed(self):
        return all([var.is_observed for var in self._flatten()])

    def update_observed_submodel(self):
        """
        Summary

        Parameters
        ---------
        """
        flattened_model = self._flatten()
        observed_variables = [var for var in flattened_model if var.is_observed]
        self.observed_submodel = ProbabilisticModel(observed_variables)

    def set_posterior_model(self, model, sampler=None): #TODO: Clean up code duplication
        self.posterior_model = PosteriorModel(posterior_model=model, joint_model=self)
        if sampler:
            if isinstance(sampler, ProbabilisticModel):
                self.posterior_sampler = PosteriorModel(sampler, joint_model=self)
            elif isinstance(sampler, Variable):
                self.posterior_sampler = PosteriorModel(ProbabilisticModel([sampler]), joint_model=self)
            elif isinstance(sampler, Iterable) and all([isinstance(subsampler, (ProbabilisticModel, Variable))
                                                        for subsampler in sampler]):
                self.posterior_sampler = [PosteriorModel(ProbabilisticModel([var]), joint_model=self)
                                          if isinstance(var, Variable) else PosteriorModel(var, joint_model=self)
                                          for var in sampler]
            else:
                raise ValueError("The sampler should be ither a probabilistic model, a brancher variable or an iterable of variables and/or models")

    def calculate_log_probability(self, rv_values, for_gradient=False, normalized=True):
        """
        Summary
        """
        log_probability = sum([var.calculate_log_probability(rv_values, reevaluate=False,
                                                             for_gradient=for_gradient,
                                                             normalized=normalized)
                               for var in self.variables])
        self.reset()
        return log_probability

    def _get_sample(self, number_samples, observed=False, input_values={}):
        """
        Summary
        """
        joint_sample = join_dicts_list([var._get_sample(number_samples=number_samples, resample=False,
                                                        observed=observed, input_values=input_values)
                                        for var in self.variables])
        joint_sample.update(input_values)
        self.reset()
        return joint_sample

    def get_sample(self, number_samples, input_values={}):
        reformatted_input_values = reformat_sampler_input(pandas_frame2dict(input_values),
                                                                            number_samples=number_samples)
        raw_sample = self._get_sample(number_samples, observed=False, input_values=reformatted_input_values)
        sample = reformat_sample_to_pandas(raw_sample, number_samples=number_samples)
        return sample

    def check_posterior_model(self):
        """
        Summary
        """
        if not self.posterior_model:
            raise AttributeError("The posterior model has not been initialized.")
#        elif self.posterior_model._is_trained is False:
#            raise AttributeError("The posterior model needs to be trained before sampling.")

    def _get_posterior_sample(self, number_samples, input_values={}):
        """
        Summary
        """
        self.check_posterior_model()
        posterior_sample = self.posterior_model._get_posterior_sample(number_samples=number_samples,
                                                                      input_values=input_values)
        sample = self._get_sample(number_samples, input_values=posterior_sample)
        return sample

    def get_posterior_sample(self, number_samples, input_values={}):
        reformatted_input_values = reformat_sampler_input(pandas_frame2dict(input_values),
                                                                            number_samples=number_samples)
        raw_sample = self._get_posterior_sample(number_samples, input_values=reformatted_input_values)
        sample = reformat_sample_to_pandas(raw_sample, number_samples=number_samples)
        return sample

    def get_p_and_q_log_probabilities(self, q_samples, q_model, empirical_samples={},
                                      for_gradient=False, normalized=True):  #TODO: Work in progress
        q_log_prob = q_model.calculate_log_probability(q_samples,
                                                       for_gradient=for_gradient, normalized=normalized)
        p_samples = reassign_samples(q_samples, source_model=q_model, target_model=self)
        p_samples.update(empirical_samples)
        p_log_prob = self.calculate_log_probability(p_samples, for_gradient=for_gradient, normalized=normalized)
        return q_log_prob, p_log_prob

    def get_importance_weights(self, q_samples, q_model, empirical_samples={},
                               for_gradient=False, give_normalization=False):
        if not empirical_samples:
            empirical_samples = self.observed_submodel._get_sample(1, observed=True)
        q_log_prob, p_log_prob = self.get_p_and_q_log_probabilities(q_samples=q_samples,
                                                                    q_model=q_model,
                                                                    empirical_samples=empirical_samples,
                                                                    for_gradient=for_gradient,
                                                                    normalized=False)
        log_weights = (p_log_prob - q_log_prob).data
        alpha = np.max(log_weights)
        norm_log_weights = log_weights - alpha
        weights = np.exp(norm_log_weights) #TODO: for Julia: this should either be numpy or cupy
        norm = np.mean(weights)
        weights /= norm
        if not give_normalization:
            return weights
        else:
            return weights, norm*np.exp(alpha)

    def estimate_log_model_evidence(self, number_samples, method="ELBO", input_values={}, for_gradient=False, posterior_model=()):
        if not posterior_model:
            self.check_posterior_model()
            posterior_model = self.posterior_model
        if method is "ELBO":
            empirical_samples = self.observed_submodel._get_sample(1, observed=True) #TODO: You need to correct for subsampling
            posterior_samples = posterior_model._get_sample(number_samples=number_samples,
                                                            observed=False, input_values=input_values)
            posterior_log_prob, joint_log_prob = self.get_p_and_q_log_probabilities(q_samples=posterior_samples,
                                                                                empirical_samples=empirical_samples,
                                                                                for_gradient=for_gradient,
                                                                                q_model=posterior_model)
            log_model_evidence = F.mean(joint_log_prob - posterior_log_prob)
            return log_model_evidence
        else:
            raise NotImplementedError("The requested estimation method is currently not implemented.")

    def reset(self):
        """
        Summary
        """
        for variable in self.variables:
            variable.reset()

    def _flatten(self):
        return flatten_list([var._flatten() for var in self.variables])

class PosteriorModel(ProbabilisticModel):
    """
    Summary

    Parameters
    ----------
    variables : tuple of brancher variables
        Summary
    """
    def __init__(self, posterior_model, joint_model):
        super().__init__(posterior_model.variables)
        self.posterior_model = None
        self.model_mapping = get_model_mapping(self, joint_model)

        self._is_trained = False

    def posterior_sample2joint_sample(self, posterior_sample):
        return reassign_samples(posterior_sample, self.model_mapping)

    def _get_posterior_sample(self, number_samples, observed=False, input_values={}):
        sample = self.posterior_sample2joint_sample(self._get_sample(number_samples, observed, input_values))
        sample.update(input_values)
        return sample

def var2link(var):
    if isinstance(var, Variable):
        vars = {var}
        fn = lambda values: values[var]
    elif isinstance(var, (numbers.Number, np.ndarray)):
        vars = {}
        fn = lambda values: var
    elif isinstance(var, (tuple, list)) and all([isinstance(v, (Variable, PartialLink)) for v in var]):
        vars = join_sets_list([{v} if isinstance(v, Variable) else v.vars for v in var])
        fn = lambda values: tuple([values[v] if isinstance(v, Variable) else v.fn(values) for v in var])
    else:
        return var
    return PartialLink(vars=vars, fn=fn, links=set())

class PartialLink(BrancherClass): #TODO: This should become "ProbabilisticProgram?"

    def __init__(self, vars, fn, links):
        self.vars = vars
        self.fn = fn
        self.links = links

    def _apply_operator(self, other, op):
        other = var2link(other)
        return PartialLink(vars=self.vars.union(other.vars),
                           fn=lambda values: op(self.fn(values), other.fn(values)),
                           links=self.links.union(other.links))

    def __add__(self, other):
        return self._apply_operator(other, operator.add)

    def __radd__(self, other):
        return self.__add__(other)

    def __sub__(self, other):
        return self._apply_operator(other, operator.sub)

    def __rsub__(self, other):
        return -1*self.__sub__(other)

    def __mul__(self, other):
        return self._apply_operator(other, operator.mul)

    def __rmul__(self, other):
        return self.__mul__(other)

    def __truediv__(self, other):
        return self._apply_operator(other, operator.truediv)

    def __rtruediv__(self, other):
        raise NotImplementedError

    def __pow__(self, other):
        return self._apply_operator(other, operator.pow)

    def __rpow__(self, other):
        raise NotImplementedError

    def __getitem__(self, key):
        if isinstance(key, collections.Iterable) and all([isinstance(k, int) for k in key]):
            variable_slice = (slice(None, None, None), *key)
        elif isinstance(key, int):
            variable_slice = (slice(None, None, None), key)
        elif isinstance(key, collections.Hashable):
            variable_slice = key
        else:
            raise ValueError("The input to __getitem__ is neither numeric nor a hashabble key")

        vars = self.vars
        fn = lambda values: self.fn(values)[variable_slice]
        links = set()
        return PartialLink(vars=vars,
                           fn=fn,
                           links=self.links)

    def shape(self):
        vars = self.vars
        fn = lambda values: self.fn(values).shape
        links = set()
        return PartialLink(vars=vars,
                           fn=fn,
                           links=self.links)

    def _flatten(self):
        return flatten_list([var._flatten() for var in self.vars]) + [self]

ext: py

reason: import numpy


text:
import argparse
import os

import numpy as np
from scipy.interpolate import interp1d
from collections import defaultdict
from shutil import copy

def post_process_track_and_save(txt_file, post_processed_txt_file,  min_track_length_filter=2):
    track_dict = defaultdict(list)
    with open(txt_file, 'r+') as f1:
        for line in f1:
            infos = line.rstrip()
            infos = infos.split(',')
            frame = int(infos[0])
            tid = int(round(float(infos[1])))
            bbx = [float(infos[2]), float(infos[3]), float(infos[4]), float(infos[5])]

            track_dict[tid].append(
                [frame, bbx]
            )

    smoothed_tracklet = post_processing_trajectory_smoothness(track_dict, min_track_length_filter=min_track_length_filter)

    track_infos = []
    for frame, bbxs in smoothed_tracklet.items():
        bbxs = sorted(bbxs, key=lambda x: x[1])
        for bbx in bbxs:
            track_id = bbx[1]
            bb_left = bbx[0][0]
            bb_top = bbx[0][1]
            bb_width = bbx[0][2]
            bb_height = bbx[0][3]
            track_infos.append(
                f"{frame},{track_id},{bb_left},{bb_top},{bb_width},{bb_height},-1,-1,-1,-1\n"
            )

    with open(post_processed_txt_file, 'w') as f1:
        f1.writelines(track_infos)

def post_processing_trajectory_smoothness(track_res, min_track_length_filter=2):
    def interpolate_boxes(frames, boxes, target_frames):
        uniq_idx = [0] + [i for i in range(1, len(frames)) if frames[i] > frames[i - 1]]
        uniq_frames = [frames[i] for i in uniq_idx]
        if len(uniq_frames) == 1:
            return target_frames, boxes[uniq_idx]

        ret = []
        for i in range(boxes.shape[1]):
            f = interp1d(uniq_frames, boxes[uniq_idx, i], assume_sorted=True)
            ret.append(f(target_frames))
        return target_frames, np.array(ret).T

    track_res_stats = {}
    for track_id, item in track_res.items():
        frames = []
        bboxes = []

        for frame_bbox in item:
            frame, bbox = frame_bbox

            frames.append(frame)
            bboxes.append(bbox)

        track_res_stats.setdefault(track_id, {})
        track_res_stats[track_id]["frames"] = frames
        track_res_stats[track_id]["boxes"] = bboxes

    track_res_smoothness = {}
    for track_id in track_res_stats:
        frames = track_res_stats[track_id]['frames']
        assert sorted(frames) == frames
        if len(frames) < min_track_length_filter:
            continue
        boxes = track_res_stats[track_id]['boxes']
        sampled_frames = [ii for ii in range(min(frames), max(frames) + 1)]
        tframes, tboxes = interpolate_boxes(np.array(frames), np.array(boxes), sampled_frames)
        valid = [i for i, box in enumerate(tboxes) if not np.any(np.isnan(box))]
        tframes = [tframes[i] for i in valid]
        tboxes = [tboxes[i].tolist() for i in valid]
        if len(tframes) < min_track_length_filter:
            continue
        for ii, frame in enumerate(tframes):
            track_res_smoothness.setdefault(frame, []).append([tboxes[ii], track_id])
    return track_res_smoothness

if __name__ == '__main__':
    parser = argparse.ArgumentParser()
    parser.add_argument("--predict_dir", type=str, required=True, help="predict directory")
    parser.add_argument('--out_path', type=str, required=True, help="output directory")
    parser.add_argument("--min_track_length_filter", type=int, default=2, help="maximum tracklet frame gap can be tolerate")
    parser.add_argument("--not_include_move_camera", type=int)
    args = parser.parse_args()

    if not os.path.exists(args.out_path):
        os.makedirs(args.out_path, exist_ok=True)

    detectors = ['SDP', 'DPM', 'FRCNN']
    static_camera_ids = ['01', '02', '03', '04', '08', '09']
    move_camera_ids = ['05', '06', '07', '10', '11', '12', '13', '14']

    static_cameras = ['MOT17-{}-{}'.format(static_camera_id, detector) for detector in detectors for static_camera_id in
                      static_camera_ids]
    move_cameras = ['MOT17-{}-{}'.format(move_camera_id, detector) for detector in detectors for move_camera_id in
                    move_camera_ids]

    for result_txt in os.listdir(args.predict_dir):
        if result_txt.endswith(".txt"):
            video_name = os.path.basename(result_txt).split('.')[0]

            output_txt_file = os.path.join(args.out_path, result_txt)
            predict_txt_file = os.path.join(args.predict_dir, result_txt)

            if video_name in move_cameras and args.not_include_move_camera == 1:
                copy(predict_txt_file, output_txt_file)
            elif "MOT20" in video_name:
                post_process_track_and_save(predict_txt_file, output_txt_file, args.min_track_length_filter)
            else:
                post_process_track_and_save(predict_txt_file, output_txt_file, args.min_track_length_filter)

reason: import numpy,from scipy


text:
'''

 Copyright (c) 2019 Fractus IT d.o.o. <http://fractus.io>

'''

# importing the libraries
from keras.models import model_from_json
import numpy as np

def main():

    # define dataset
    X = np.array([[0, 0], [0, 1], [1, 0], [1, 1]])
    Y = np.array([[0], [1], [1], [0]])

    print("Loaingd model from disk...")

    # load json and create model
    json_file = open('model.json', 'r')
    loaded_model_json = json_file.read()

    json_file.close()
    model = model_from_json(loaded_model_json)
    # load weights into new model
    model.load_weights("model.h5")
    print("Loaded model from disk")

    model.compile(optimizer = 'adam', loss = 'binary_crossentropy', metrics = ['accuracy'])

        # evaluate the model
    scores = model.evaluate(X, Y)

    # print accuracy
    print("\n%s: %.2f%%" % (model.metrics_names[1], scores[1]*100))

    print("predicting output:")
    print(model.predict_proba(X))

    print("model summary:")
    print(model.summary())

if __name__ == '__main__':
    main()

max_stars_repo_name: fractus-io/keras-tutorial

max_issues_repo_name: fractus-io/keras-tutorial

max_forks_repo_name: fractus-io/keras-tutorial


text:
"""
Calculations provided by aiida_sirius.

Register calculations via the "aiida.calculations" entry point in setup.json.
"""
from __future__ import absolute_import

import json
import tempfile
from copy import deepcopy

import numpy as np
import six

from aiida.common import datastructures
# from aiida.common.constants import bohr_to_ang
from aiida.engine import CalcJob
from aiida.orm import StructureData, UpfData
from aiida.plugins import DataFactory

from ..upf_to_json import upf_to_json

SiriusParameters = DataFactory('sirius.scf')
SinglefileData = DataFactory('singlefile')
KpointsData = DataFactory('array.kpoints')
Dict = DataFactory('dict')

# CODATA 2018
bohr_to_ang = 0.529177210903

def read_structure(structure, magnetization):
    """
    Convert obj of type `StructureData` to sirius input format (e.g. a dictionary)
    Return:
    - cell
    - atomic coordinates
    """
    cell = structure.cell
    # cell in atomic units
    cell = np.array(cell) / bohr_to_ang
    # convert back to list of list
    cell = [list(x) for x in cell]
    sites = structure.sites
    atom_types = set([x.kind_name for x in sites])
    atomic_coordinates = {}
    for atom_type in atom_types:
        ltypes = list(filter(lambda x: x.kind_name == atom_type, sites))
        lpos = np.array([x.position for x in ltypes])
        try:
            lmag = np.array(magnetization[atom_type])
        except KeyError:
            lmag = np.zeros_like(lpos)
        lpos_with_mag = np.concatenate((lpos, lmag), axis=1)
        atomic_coordinates[atom_type] = [list(x) for x in lpos_with_mag]
    return cell, atomic_coordinates

def make_sirius_json(parameters, structure, kpoints, magnetization):
    """
    DEPRECATED
    Keyword Arguments:
    structure -- structure
    pseudos   -- dictionary of UpfData
    """
    raise Exception('deprecated')
    sirius_cell, sirius_pos = read_structure(structure, magnetization.get_dict())
    sirius_json = {'parameters': deepcopy(parameters)}
    sirius_json['unit_cell']['lattice_vectors'] = sirius_cell
    sirius_json['unit_cell']['atom_types'] = list(sirius_pos.keys())
    sirius_json['unit_cell']['atom_coordinate_units'] = 'A'
    sirius_json['unit_cell']['atoms'] = sirius_pos
    if 'mesh' in kpoints.attributes:
        sirius_json['parameters']['ngridk'] = kpoints.attributes['mesh']
        sirius_json['parameters']['shiftk'] = kpoints.attributes['offset']
    else:
        print('setting kpoints as list(list)')
        print(kpoints.get_kpoints())
        sirius_json['parameters']['vk'] = [list(x) for x in kpoints.get_kpoints()]

    return sirius_json

def add_cell_kpoints_mag_to_sirius(sirius_params, structure, magnetization, kpoints):
    """
    Add atoms, lattice vectors, kpoints and magnetization to sirius_params.

    :param sirius_params: sirius json
    :param structure: aiida structure
    :param magnetization: {ATOM: [x, y, z]},
        note magnetization is given as a 3d vector
    :param kpoints: aiida kpoints
    """
    sout = deepcopy(sirius_params)

    sout['unit_cell'] = {}
    # add lattice vectors (unit is bohr)
    sout['unit_cell']['lattice_vectors'] = [
        list(x) for x in np.array(structure.attributes['cell']) / bohr_to_ang
    ]
    sout['unit_cell']['lattice_vectors_scale'] = 1

    elems = set([site['kind_name'] for site in structure.attributes['sites']])

    sout['unit_cell']['atom_types'] = list(elems)
    sout['unit_cell']['atom_files'] = {elem: elem + '.json' for elem in elems}
    sout['unit_cell']['atom_coordinate_units'] = 'A'

    if 'mesh' in kpoints.attributes:
        sout['parameters']['ngridk'] = kpoints.attributes['mesh']
        sout['parameters']['shiftk'] = kpoints.attributes['offset']
    else:
        sout['parameters']['vk'] = [list(x) for x in kpoints.get_array('kpoints')]

    sout['unit_cell']['atoms'] = {}
    for atom_type in elems:
        angstrom_coords = []
        for site in filter(lambda x: x['kind_name'] == atom_type,
                           structure.attributes['sites']):
            angstrom_coords.append(site['position'])

        if atom_type in dict(magnetization):
            coords_magnetization = np.hstack((np.array(angstrom_coords),
                                              np.array(magnetization[atom_type])))
        else:
            coords_magnetization = angstrom_coords

        sout['unit_cell']['atoms'][atom_type] = [
            list(x) for x in coords_magnetization
        ]

    return sout

class SiriusBaseCalculation(CalcJob):
    """
    AiiDA calculation plugin wrapping the diff executable.

    Simple AiiDA plugin wrapper for 'diffing' two files.
    """

    @classmethod
    def define(cls, spec):
        """Define inputs and outputs of the calculation."""
        # yapf: disable
        super(SiriusBaseCalculation, cls).define(spec)
        spec.input('metadata.options.resources', valid_type=dict, default={'num_machines': 1, 'num_mpiprocs_per_machine': 1})
        spec.input('structure', valid_type=StructureData, help='The input structure')
        spec.input('kpoints', valid_type=KpointsData, help='kpoints')
        spec.input('sirius_config', valid_type=SiriusParameters, help='sirius parameters')
        spec.input('magnetization', valid_type=Dict, default=Dict(dict={}))
        spec.input_namespace('pseudos', valid_type=UpfData, dynamic=True,
                             help='A mapping of `UpfData` nodes onto the kind name to which they should apply.')

    def _read_pseudos(self, sirius_json):
        # parse pseudos
        for atom_label in self.inputs.pseudos:
            upf = self.inputs.pseudos[atom_label]
            with upf.open() as fh:
                upf_json = upf_to_json(fh.read(), fname=upf.filename)
            sirius_json['unit_cell']['atom_files'][atom_label] = json.dumps(upf_json)
        return sirius_json

class SiriusSCFCalculation(SiriusBaseCalculation):
    """
    AiiDA calculation plugin wrapping the diff executable.

    Simple AiiDA plugin wrapper for 'diffing' two files.
    """

    @classmethod
    def define(cls, spec):
        """Define inputs and outputs of the calculation."""
        # yapf: disable
        super(SiriusSCFCalculation, cls).define(spec)
        spec.input('metadata.options.parser_name', valid_type=six.string_types, default='sirius.scf')
        spec.input('metadata.options.output_filename', valid_type=six.string_types, default='sirius.scf.out')
        spec.output('sirius', valid_type=SinglefileData, help='standard output')
        spec.output('output', valid_type=Dict, help='sirius.scf json output')
        spec.exit_code(100, 'ERROR_MISSING_OUTPUT_FILES', message='Calculation did not produce all expected output files.')
        spec.exit_code(201, 'ERROR_NOT_CONVERGED', message='Calculation did not converge.')

    def prepare_for_submission(self, folder):
        """
        Create input files.
        sirius.json,
        nlcg.yaml

        :param folder: an `aiida.common.folders.Folder` where the plugin should temporarily place all files needed by
            the calculation.
        :return: `aiida.common.datastructures.CalcInfo` instance
        """
        codeinfo = datastructures.CodeInfo()
        output_filename = self.metadata.options.output_filename
        codeinfo.cmdline_params = ['--output=output.json']
        codeinfo.code_uuid = self.inputs.code.uuid
        codeinfo.stdout_name = self.metadata.options.output_filename
        codeinfo.withmpi = self.inputs.metadata.options.withmpi

        # with config from input
        structure = self.inputs.structure
        kpoints = self.inputs.kpoints
        magnetization = self.inputs.magnetization
        sirius_json = make_sirius_json(self.inputs.sirius_config.get_dict()['parameters'],
                                       structure, kpoints, magnetization)
        with tempfile.NamedTemporaryFile(mode='w', suffix='.json', delete=False) as sirius_tmpfile:
            sirius_json = self._read_pseudos(sirius_json)
            sirius_tmpfile_name = sirius_tmpfile.name
            # merge with settings given from outside
            sirius_json = {**sirius_json, **self.inputs.sirius_config.get_dict()}
            # dump to file
            json.dump(sirius_json, sirius_tmpfile)
        sirius_config = SinglefileData(file=sirius_tmpfile_name)
        sirius_config.store()

        # Prepare a `CalcInfo` to be returned to the engine
        calcinfo = datastructures.CalcInfo()
        calcinfo.codes_info = [codeinfo]
        calcinfo.local_copy_list = [
            (sirius_config.uuid, sirius_config.filename, 'sirius.json')
        ]
        calcinfo.retrieve_list = [self.metadata.options.output_filename, 'output.json']

        return calcinfo

max_stars_repo_name: simonpintarelli/aiida-sirius

max_issues_repo_name: simonpintarelli/aiida-sirius

max_forks_repo_name: simonpintarelli/aiida-sirius


text:
import numpy as np
import graph_utils

def read_data_file(filename):
    with open(filename, "r") as f:
        s = f.read()
    s = s.strip()
    vals = []
    for v in s.split():
        try:
            vals.append(int(v))
        except ValueError:
            continue
    return vals

def printstats(filename):
    vals = read_data_file(filename)
    print("Mean:\t%.3f" % np.mean(vals))
    print("Stddev:\t%.3f" % np.std(vals))

def countbig(filename, big=150):
    vals = read_data_file(filename)
    big_vals = 0
    for v in vals:
        if v >= big:
            big_vals +=1

    return big_vals

def plot_cdf_set_all(exps,
                     out_fname,
                     title="Graph title",
                     xlabel="X Label", ylabel="Y Label",
                     xlim=None, ylim=None,
                     legend_title="Config",
                     cmap_n="tab10"):
    cdfs = {}
    for expname in exps:
        exp = exps[expname]
        if "filename" in exp:
            fname = exp["filename"]
            vals = read_data_file(fname)
        if "data" in exp:
            vals = exp["data"]

        ex, ey = graph_utils.get_cdf(vals)
        cdfs[expname] = {"label": expname,
                         "line": (ex, ey)}

    graph_utils.plot_results(cdfs, out_fname,
                             title=title,
                             xlabel=xlabel,
                             ylabel=ylabel,
                             xlim=xlim,
                             ylim=ylim,
                             cmap_n=cmap_n,
                             markers=None,
                             markersize=0.0,
                             legend_title=legend_title,
                             drawstyle='steps-post')

import yaml
def plot_ethtool_delta(yaml_before, yaml_after):
    with open(yaml_before, "r") as f:
        before = yaml.load(f)["NIC statistics"]
    with open(yaml_after, "r") as f:
        after = yaml.load(f)["NIC statistics"]

    for k in before:
        try:
            b = int(before[k])
            a = int(after[k])
            if a - b != 0:
                print("%s: %d" % (k, a-b))
        except:
            continue

max_stars_repo_name: TritonNetworking/ondemand-netstack

max_issues_repo_name: TritonNetworking/ondemand-netstack

max_forks_repo_name: TritonNetworking/ondemand-netstack


text:
import numpy as np
from matplotlib import pyplot as plt
from matplotlib import animation

# First set up the figure, the axis, and the plot element we want to animate
fig, ax = plt.subplots(subplot_kw=dict(projection='polar', frameon=True), figsize=(15, 15))

# initialization function: plot the background of each frame
def init():
    ax.set_theta_zero_location("N")
    ax.set_theta_direction(-1)
    ax.get_xaxis().set_visible(False)
    ax.get_yaxis().set_visible(False)
    return ax

# animation function.  This is called sequentially
def animate(i):
    fig1, ax1, cax1 = vis.plot_polar_contour(cop_proc[i], cmap='spring', levels=200)
    cax = cax1
    return cax

# call the animator.  blit=True means only re-draw the parts that have changed.
anim = animation.FuncAnimation(fig, animate, init_func=init,
                               frames=30, interval=2000, blit=True)
anim.save('test.html')


text:
import numpy as np
from simplification.cutil import simplify_coords
import matplotlib.pyplot as plt

def give_ratio_index2(seq):
    val = np.zeros(len(seq))*np.nan
    val[np.array(seq) == "T"] = 0
    val[np.array(seq) == "B"] = 1
    index = [ind for ind, v in enumerate(val) if v in [0, 1]]
    return index, val[~np.isnan(val)]

def Tsites(seq) :
    val = np.zeros(len(seq))*np.nan
    val[np.array(seq) == "T"] = 1
    val[np.array(seq) == "B"] = 1
    val[np.array(seq) == "X"] = 1
    index = [ind for ind, v in enumerate(val) if v in [1]]
    return index

def runningMean(x, N):
    filt = np.convolve(x, np.ones((N,))/N)[(N-1):]
    filt1 = np.convolve(x, np.ones((N,))/N)[:-N+1]
    filt = 0.5*(filt+filt1)
    return filt

def Simplify(x,y, stdev,smoothing, MinAmplitude,Sparam):
    Xs = []
    Ys =[]
    Yr=[]
    if np.std(y)> stdev :
        Y = runningMean(y,smoothing)
        if max(Y)-min(Y) > MinAmplitude :
            Yr = (np.array(Y)-min(Y))/(max(Y)-min(Y)) # rescale between 0 and 1
            coords = list(zip(x,Yr))
            simplified = np.array(simplify_coords(coords, Sparam))
            simplified = simplified[simplified[:,0].argsort()]
            Xs = simplified[:,0]
            Ys = simplified[:,1]
            for i in range(0,len(Xs)-2) :
                if abs(Xs[i+2]-Xs[i]) < 1000 and abs(Ys[i+2]-Ys[i]) < 0.3 and (abs(Ys[i+2]-Ys[i+1])>0.6 or abs(Ys[i]-Ys[i+1])>0.6):
                    Ys[i+1] = (Ys[i+2]+Ys[i])/2
    return Xs,Ys,Yr

def Detection(Xs,Ys,MinJump):
        slope = {}
        jump = {}
        Tracts = {}
        ID = 0
        Counted = []
        for a in range(0, len(Xs)-1) :
            slope[a] = (Ys[a+1]-Ys[a])/(Xs[a+1]-Xs[a])*100000
            jump[a] = (Ys[a+1]-Ys[a])
        for i in range(0, len(Xs)-2) :
            if jump[i] > MinJump or (jump[i]+jump[i+1]>MinJump and jump[i] <= MinJump and jump[i+1] <= MinJump):
                startTract = Xs[i]
                rest = 0
                for k in range(i+1, len(Xs)-1) :
                    if jump[k]+rest < -MinJump*0.6:
                        try :
                            if jump[k+1]<0 and slope[k+1] < -50 :
                                k+=1 # problem when decrease is splitted into 2 tracts
                        except : pass
                        endTract = Xs[k+1]
                        Counted.append(k)
                        Tracts[ID] = {}
                        Tracts[ID]['start'] = startTract
                        Tracts[ID]['end'] = endTract
                        Tracts[ID]['JumpScore'] = max((jump[i]-MinJump)*10, (jump[i]+jump[i+1]-MinJump)*10, (-jump[k]-rest-MinJump)*10)
                        Oriented = False
                        try :
                            if endTract == Tracts[ID-1]['end'] and Tracts[ID-1]['dir'] == '?':
                                ID -=1
                                Tracts[ID]['start'] = startTract
                                del Tracts[ID+1]
                        except : pass
                        maxi = Ys[i:k+1].argmax()
                        slope1 = (Ys[i+maxi]-Ys[i])/(Xs[i+maxi]-Xs[i])*100000
                        slope2 = (Ys[k+1]-Ys[i+maxi])/(Xs[k+1]-Xs[i+maxi])*100000

                        if -slope2*1.5 < slope1:
                            Tracts[ID]['dir'] = 'R'
                            S = slope1/-slope2
                            Y = Ys[k+1]-Ys[i]
                            Tracts[ID]['score'] = max(S,20*Y)
                            Oriented = True

                        if Ys[k+1]-Ys[i] > 0.075 and k+1 != len(Xs)-1:
                            Tracts[ID]['dir'] = 'R'
                            S = slope1/-slope2
                            Y = Ys[k+1]-Ys[i]
                            Tracts[ID]['score'] = max(S,20*Y)
                            Oriented = True

                        if Ys[i]-Ys[k+1] > 0.075 and i !=0 :
                            if Oriented == False :
                                Tracts[ID]['dir'] = 'L'
                                S = slope2/slope1
                                Y = Ys[k+1]-Ys[i]
                                Tracts[ID]['score'] = min(S,20*Y)
                            else  :
                                Tracts[ID]['dir'] = '?'
                                Tracts[ID]['score'] = 0

                        if slope1*1.5< -slope2:
                            if Oriented == False :
                                Tracts[ID]['dir'] = 'L'
                                S = slope2/slope1
                                Y = Ys[k+1]-Ys[i]
                                Tracts[ID]['score'] = min(S,20*Y)
                            else  :
                                Tracts[ID]['dir'] = '?'
                                Tracts[ID]['score'] = 0

                        if 'dir' not in Tracts[ID].keys() :
                            Tracts[ID]['dir'] = '?'
                            Tracts[ID]['score'] = 0
                        ID+=1
                        break
                    rest+=jump[k]

        for k in range(0, len(Xs)-2) :
            if k not in Counted :
                if jump[k] < -MinJump or (jump[k]+jump[k+1]<-MinJump and jump[k]>=-MinJump and jump[k+1]>=-MinJump ):
                    endTract = Xs[k+1]
                    rest = 0
                    for i in range(0, k)[::-1] :
                        if jump[i]+rest > MinJump*0.6:
                            startTract = Xs[i]
                            conflict = False
                            maxi = Ys[i:k+1].argmax()
                            slope1 = (Ys[i+maxi]-Ys[i])/(Xs[i+maxi]-Xs[i])*100000
                            slope2 = (Ys[k+1]-Ys[i+maxi])/(Xs[k+1]-Xs[i+maxi])*100000
                            S = slope2/slope1
                            if k+1 == len(Xs)-1 or i == 0 :
                                Y = 0
                            else : Y = Ys[k+1]-Ys[i]
                            score = min(S,20*Y)

                            for item in Tracts :
                                if startTract == Tracts[item]['start'] :
                                    ID = item
                                    if  Tracts[item]['dir'] =='R' :
                                        if -score*2 > Tracts[item]['score'] :
                                            if -score < 2*Tracts[item]['score'] :
                                                conflict = True
                                            else : conflict = False
                                        if -score*2 <= Tracts[item]['score'] : conflict = 'Solved'
                                    if  Tracts[item]['dir'] =='?' :
                                        conflict = True
                                    if  Tracts[item]['dir'] =='L' :
                                        if Tracts[item]['score'] < score : conflict = 'Solved'
                                        else :
                                            del Tracts[item]
                                    break

                                if startTract < Tracts[item]['start']:
                                    ID = item
                                    if endTract <= Tracts[item]['start']:
                                        for it2 in range(item, max(Tracts.keys())+1)[::-1] :
                                            Tracts[it2+1] = Tracts[it2]
                                        break
                                    else :
                                        if  Tracts[item]['dir'] =='R' :
                                            if -score*2 > Tracts[item]['score'] :
                                                if -score < 2*Tracts[item]['score'] :conflict = True
                                                else : conflict = False
                                            if -score*2 <= Tracts[item]['score'] : conflict = 'Solved'
                            if Tracts!={} and startTract > Tracts[max(Tracts.keys())]['start']:
                                ID=max(Tracts.keys())+1

                            if conflict == False:
                                Tracts[ID] = {}
                                Tracts[ID]['start'] = startTract
                                Tracts[ID]['end'] = endTract
                                if Ys[i]-Ys[k+1] > 0.075 or slope1*1.5< -slope2:
                                    Tracts[ID]['dir'] = 'L'
                                    Tracts[ID]['score'] = min(S,20*Y)
                                else:
                                    Tracts[ID]['dir'] = '?'
                                    Tracts[ID]['score'] = 0
                                Tracts[ID]['JumpScore'] = max((-jump[k]-MinJump)*10, (-jump[k]-jump[k+1]-MinJump)*10, (jump[i]+rest-MinJump)*10)

                            if conflict == True:
                                Tracts[ID] = {}
                                Tracts[ID]['start'] = startTract
                                Tracts[ID]['end'] = endTract
                                Tracts[ID]['dir'] = '?'
                                Tracts[ID]['score'] = 0
                                Tracts[ID]['JumpScore'] = max((-jump[k]-MinJump)*10, (-jump[k]-jump[k+1]-MinJump)*10, (jump[i]+rest-MinJump)*10)

                            ID+=1
                            break
                        rest+=jump[i]

        Finish = False
        while Finish == False: # removing duplicates
            TractCopy = Tracts.copy()
            Finish = True
            for item in TractCopy:
                for item2 in TractCopy:
                    if item != item2 and item in Tracts:
                        if  TractCopy[item]['start'] >= TractCopy[item2]['start'] and TractCopy[item]['start'] < TractCopy[item2]['end'] :
                            Finish = False
                            if abs(TractCopy[item]['score']) > abs(TractCopy[item2]['score']) :
                                del Tracts[item2]
                                for it in range(item2,len(Tracts)):
                                    Tracts[it] = Tracts[it+1]
                            else :
                                del Tracts[item]
                                for it in range(item,len(Tracts)):
                                    Tracts[it] = Tracts[it+1]
        return Tracts

def DetectInits(Tracts, Xs, X, Yr, LowPlateau,MinDist, score, jumpscore):
    Initiations = []
    if Tracts !={} :
        for ID in range(min(Tracts.keys()), max(Tracts.keys())):
            if abs(Tracts[ID+1]['JumpScore']) > jumpscore/3.0 and abs(Tracts[ID]['JumpScore']) > jumpscore/3.0 :
                if (Tracts[ID]['dir'] == 'L' and Tracts[ID]['score']<-score  and (Tracts[ID+1]['dir'] == 'R' or (Tracts[ID+1]['dir'] != 'R'and abs(Tracts[ID+1]['score'])<score and abs(Tracts[ID+1]['JumpScore']) > jumpscore and abs(Tracts[ID]['JumpScore']) > jumpscore))) or (Tracts[ID+1]['dir'] == 'R' and Tracts[ID+1]['score']>score and (Tracts[ID]['dir'] == 'L' or (Tracts[ID]['dir'] != 'L'and abs(Tracts[ID]['score'])>-score and abs(Tracts[ID+1]['JumpScore']) > jumpscore and abs(Tracts[ID]['JumpScore']) > jumpscore))):
                        INIT = (Tracts[ID+1]['start'] + Tracts[ID]['end'])/2
                        I1 = (np.abs(X-Tracts[ID]['end'])).argmin()
                        I2 = (np.abs(X-Tracts[ID+1]['start'])).argmin()
                        M = []
                        for x in range(min(I1, I2), max(I1,I2)) :
                            M.append(Yr[x])
                        if  M != [] and np.array(M).mean()< LowPlateau and Tracts[ID+1]['start'] - Tracts[ID]['end'] > MinDist:
                            Initiations.append([INIT,Tracts[ID]['start'], Tracts[ID]['end'],Tracts[ID+1]['start'], Tracts[ID+1]['end']])
                            if Tracts[ID]['dir'] == '?' or Tracts[ID]['dir']=='R':  # correct tract direction when initiation is detected
                                Tracts[ID]['dir'] = 'L'
                                Tracts[ID]['score'] = 5
                                Tracts[ID]['JumpScore'] = 5
                            if Tracts[ID+1]['dir']== '?' or Tracts[ID+1]['dir']=='L':
                                Tracts[ID+1]['dir'] = 'R'
                                Tracts[ID+1]['score'] = 5
                                Tracts[ID+1]['JumpScore'] = 5
    return Initiations, Tracts

def DetectTermsFilter3(Tracts, Xs, X, Yr, LowPlateau, jumpscore, score):
    Terminaisons = []
    if Tracts !={} :
        for ID in range(min(Tracts.keys()), max(Tracts.keys())):
            if abs(Tracts[ID+1]['JumpScore']) > jumpscore/3.0 and abs(Tracts[ID]['JumpScore']) > jumpscore/3.0 :
                if (Tracts[ID+1]['dir'] == 'L' and Tracts[ID+1]['score']<-score  and (Tracts[ID]['dir'] == 'R' or (Tracts[ID]['dir'] != 'R'and abs(Tracts[ID]['score'])<score and abs(Tracts[ID]['JumpScore']) > jumpscore and abs(Tracts[ID+1]['JumpScore']) > jumpscore))) or (Tracts[ID]['dir'] == 'R' and Tracts[ID]['score']>score and (Tracts[ID+1]['dir'] == 'L' or (Tracts[ID+1]['dir'] != 'L'and abs(Tracts[ID+1]['score'])>-score and abs(Tracts[ID]['JumpScore']) > jumpscore and abs(Tracts[ID+1]['JumpScore']) > jumpscore))):
                    TERM = (Tracts[ID]['start'] + Tracts[ID+1]['end'])/2
                    B1 = (np.abs(X-Tracts[ID]['start'])).argmin()
                    B2 = (np.abs(X-(Tracts[ID]['start']-500))).argmin()
                    A1 = (np.abs(X-Tracts[ID+1]['end'])).argmin()
                    A2 = (np.abs(X-(Tracts[ID+1]['end']+500))).argmin()
                    if (Yr[B1]< LowPlateau or Yr[B2]< LowPlateau) and (Yr[A1]< LowPlateau or Yr[A2]< LowPlateau) :
                        Terminaisons.append([TERM, Tracts[ID]['start'],Tracts[ID]['end'],Tracts[ID+1]['start'], Tracts[ID+1]['end']])
    return Terminaisons

def ExportBedForks(Tracts, output,chrom, Fasta, readname,strand, jumpscore, score):
    for ID in Tracts:
        if Tracts[ID]['dir']!='?' and Tracts[ID]['JumpScore'] >= jumpscore and Tracts[ID]['score']>=score :
            output.write( "%s\t%d\t%d\t%s\t%s\t%s\t%s\n"%(chrom,Tracts[ID]['start'], Tracts[ID]['end'], Tracts[ID]['dir'], Fasta, readname,strand))

def ExportBedForksNoFilter(Tracts, output,chrom, Fasta, readname,strand):
    for ID in Tracts:
        if Tracts[ID]['dir']!='?' :
            output.write( "%s\t%d\t%d\t%s\t%s\t%s\t%s\t%.2f\t%.2f\n"%(chrom,Tracts[ID]['start'], Tracts[ID]['end'], Tracts[ID]['dir'], Fasta, readname,strand,  Tracts[ID]['JumpScore'], Tracts[ID]['score']))

def ExportInits(InitsTable, output,chrom, Fasta, readname,strand):
    for I in InitsTable :
        output.write("%s\t%d\t%d\t%d\t%d\t%d\t%s\t%s\t%s\n" % (chrom, I[0], I[1], I[2], I[3], I[4], Fasta, readname,strand))

def Attributes(l):
    attribute = l.split('{')[1].split(',')
    Attri = {}
    for i in range(0,len(attribute)) :
        Attri[attribute[i].split(':')[0].split('\'')[1]] = attribute[i].split(':')[1].split('}')[0]
    strand = Attri['mapped_strand'].split('\'')[1]
    start = float(Attri['mapped_start'])
    end = float(Attri['mapped_end'])
    chrom = Attri['mapped_chrom'].split('\'')[1]
    return start, end, strand, chrom

def PlotTracts(Tracts, color1,color2,color3,y, head=500):
    if Tracts !={} :
        for ID in Tracts  :
            if abs(Tracts[ID]['score']) >= 2 :
                if Tracts[ID]['dir'] =='R' :
                    plt.arrow(Tracts[ID]['start'],y, Tracts[ID]['end']-Tracts[ID]['start']-head, 0, fc=color1, ec=color1,head_width=0.1, head_length=head, lw=1.5)
                if  Tracts[ID]['dir'] =='L' :
                    #print('LL1', Tracts[ID]['start'],Tracts[ID]['end'])
                    plt.arrow(Tracts[ID]['end'],y, -Tracts[ID]['end']+Tracts[ID]['start']+head, 0, fc=color1, ec=color1,head_width=0.1, head_length=head, lw=1.5)
            else :
                #print 'score', Tracts[ID]['score']
                if Tracts[ID]['dir'] =='R' :
                    plt.arrow(Tracts[ID]['start'],y, Tracts[ID]['end']-Tracts[ID]['start']-head, 0, fc=color2, ec=color2,head_width=0.1, head_length=head, lw=1.5) #, linestyle = 'dashed')
                    #print('RR', Tracts[ID]['start'],Tracts[ID]['end'])
                if  Tracts[ID]['dir'] =='L' :
                    plt.arrow(Tracts[ID]['end'],y, -Tracts[ID]['end']+Tracts[ID]['start']+head, 0, fc=color2, ec=color2,head_width=0.1, head_length=head, lw=1.5) #, linestyle = 'dashed')
                if Tracts[ID]['dir'] =='?' :
                    plt.arrow(Tracts[ID]['end'],y, -Tracts[ID]['end']+Tracts[ID]['start'], 0, fc=color2, ec=color2,head_width=0.0, head_length=0, lw=1.5) #, linestyle = 'dashed')

def PlotInits(Inits, color,y):
    if Inits!= None:
        for i in Inits:
            plt.plot(i[0], y, color,marker = '^', markersize=8, alpha = 1)

def PlotTerms(Inits, color,y):
    if Inits!= None:
        for i in Inits:
            plt.plot(i[0], y, color,marker = 'v', markersize=8, alpha = 1)

max_stars_repo_name: organic-chemistry/RepNano

max_issues_repo_name: organic-chemistry/RepNano

max_forks_repo_name: organic-chemistry/RepNano


text:
"""Define all the core classes for an RL problem."""
from abc import ABCMeta, abstractmethod, abstractproperty
from operator import itemgetter

import logging
import numpy as np

from collections import namedtuple, MutableSequence

from torch.autograd import Variable
from torch.nn import Module

from gtd.utils import cached_property, softmax

class State(object):
    pass

class HiddenState(object):
    pass

class Action(object):
    __metaclass__ = ABCMeta

    @property
    def justification(self):
        """Return a Justification object."""
        try:
            return self._justification
        except AttributeError:
            return None

    @justification.setter
    def justification(self, j):
        self._justification = j

    @abstractmethod
    def __eq__(self, other):
        raise NotImplementedError()

    @abstractmethod
    def __hash__(self):
        raise NotImplementedError()

    def __ne__(self, other):
        return not self.__eq__(other)

class Trace(object):
    __metaclass__ = ABCMeta

    @abstractmethod
    def to_json_dict(self):
        raise NotImplementedError

    @abstractmethod
    def dumps(self):
        """Return a unicode representation of the object."""
        raise NotImplementedError()

class Justification(Trace):
    """An object containing debug information explaining why an Action was chosen by the Policy."""
    __metaclass__ = ABCMeta

    @abstractmethod
    def dumps(self):
        """Return a unicode representation of the object."""
        raise NotImplementedError()

class BaseExperience(object):
    pass  # marker class

class Experience(namedtuple('Experience', ['state', 'action', 'undiscounted_reward', 'metadata']),
                 BaseExperience):
    pass

class ScoredExperience(namedtuple('ScoredExperience',
    ['state', 'action', 'undiscounted_reward', 'log_prob', 'state_value', 'metadata']),
    BaseExperience):
    """Like an Experience, but with an extra log_prob attribute.

    Attributes:
        state (State)
        action (Action)
        undiscounted_reward (float)
        log_prob (Variable): a scalar Variable
        state_value (Variable): a scalar Variable
        metadata (dict[str->any])
    """
    pass

class Episode(MutableSequence):
    """A list of Experiences."""
    __slots__ = ['_experiences']
    discount_negative_reward = False

    def __init__(self, iterable=None):
        self._experiences = []
        if iterable:
            for item in iterable:
                self.append(item)

    @classmethod
    def configure(cls, discount_negative_reward):
        cls.discount_negative_reward = discount_negative_reward

    def __getitem__(self, item):
        return self._experiences[item]

    def __setitem__(self, key, value):
        assert isinstance(value, BaseExperience)
        self._experiences[key] = value

    def __delitem__(self, key):
        del self._experiences[key]

    def __len__(self):
        return len(self._experiences)

    def insert(self, index, value):
        assert isinstance(value, BaseExperience)
        return self._experiences.insert(index, value)

    def append(self, experience):
        assert isinstance(experience, BaseExperience)
        self._experiences.append(experience)

    def discounted_return(self, t, gamma):
        """Returns G_t, the discounted return.

        Args:
            t (int): index of the episode (supports negative indexing from
                back)
            gamma (float): the discount factor

        Returns:
            float
        """
        def discounted_reward(undiscounted, index):
            if undiscounted < 0 and not Episode.discount_negative_reward:
                return undiscounted
            else:
                return undiscounted * np.power(gamma, index)

        if t < -len(self._experiences) or t > len(self._experiences):
            raise ValueError("Index t = {} is out of bounds".format(t))

        return sum(discounted_reward(experience.undiscounted_reward, i)
                   for i, experience in enumerate(self._experiences[t:]))

    def __str__(self):
        experiences = "{}".format(self._experiences)[:50]
        return "Episode({}..., undiscounted return: {})".format(
                experiences, self.discounted_return(0, 1.0))
    __repr__ = __str__

class ScoredEpisode(Episode):
    """An Episode holding ScoredExperiences, rather than plain Experiences.

    Mostly just a marker class.
    """
    pass

class ActionScores(object):
    """A map from Actions to scores (log probabilities)."""
    def __init__(self, d, state_value):
        """Construct ActionScores.

        Args:
            d (dict[Action, Variable[FloatTensor]]): a map from Actions to scalar PyTorch Variables.
            state_value (Variable[FloatTensor]): a scalar value for the state we are at
        """
        for v in d.values():
            assert isinstance(v, Variable)
        assert isinstance(state_value, Variable)

        self._vars = d
        # self._floats = {action: v.data.cpu()[0] for action, v in self._vars.items()}
        # for action, v in self._vars.items():
        #     print(v.item())
        self._floats = {action: v.item() for action, v in self._vars.items()}
        self._state_value = state_value

    @property
    def state_value(self):
        return self._state_value

    @property
    def as_floats(self):
        return self._floats

    @property
    def as_variables(self):
        return self._vars

    @property
    def best_action(self):
        if len(self.as_floats) == 0:
            return None
        return max(self.as_floats.items(), key=itemgetter(1))[0]

    def sample_action(self):
        actions, log_probs = zip(*self.as_floats.items())
        probs = softmax(log_probs)
        # TODO: Bring this back when the action score hack is fixed
        #if not np.isclose(total, 1.0):
        #    logging.warn('Action probs do not sum to 1: {} (will normalize).'.format(total))

        return np.random.choice(actions, p=probs)

    @property
    def justification(self):
        """Return a Justification object."""
        try:
            return self._justification
        except AttributeError:
            return None

    @justification.setter
    def justification(self, j):
        self._justification = j

class Policy(Module):
    """A parameterized RL policy mapping states to actions."""
    __metaclass__ = ABCMeta

    @abstractmethod
    def act(self, states, test=False):
        """Given a batch of states, return a corresponding batch of actions.

        Args:
            states (list[State])
            test (bool): True if test time (False for train time)

        Returns:
            actions (list[Action])
        """
        raise NotImplementedError()

    @abstractmethod
    def score_actions(self, states, test=False):
        """Given a batch of states, return a corresponding batch of action scores (log probs).

        Args:
            states (list[State])
            test (bool): True if test time (False for train time)

        Returns:
            action_scores_batch (list[ActionScores]): a batch of ActionScores
        """
        raise NotImplementedError()

    @abstractmethod
    def update_from_episodes(self, episodes, gamma, take_grad_step):
        """Updates the Policy based on a batch of episodes.

        Args:
            episodes (list[Episode])
            gamma (float): discount factor
            take_grad_step (Callable): takes a loss Variable and takes a
                gradient step on the loss
        """
        raise NotImplementedError()

    @abstractmethod
    def update_from_replay_buffer(self, replay, gamma):
        """Updates the Policy from sampling from the replay buffer.

        Args:
            replay (ReplayBuffer):
            gamma (float): discount factor
            take_grad_step (Callable): takes a loss Variable and takes a
                gradient step on the loss

        Returns:
            trace (Trace)
        """
        raise NotImplementedError()

    @abstractmethod
    def end_episodes(self):
        """Signals the end of episodes to reset any episode state."""
        raise NotImplementedError()

    @abstractproperty
    def has_attention(self):
        """Returns True if the property action_attention is implemented."""
        raise NotImplementedError()

    @property
    def action_attention(self):
        raise NotImplementedError()

max_stars_repo_name: SAGNIKMJR/autoweb

max_issues_repo_name: SAGNIKMJR/autoweb

max_forks_repo_name: SAGNIKMJR/autoweb