text:
from __future__ import absolute_import
import torch
import numpy as np
import copy
import multiprocessing
from pysurvival import HAS_GPU
from pysurvival import utils
from pysurvival.utils import neural_networks as nn
from pysurvival.utils import optimization as opt
from pysurvival.models import BaseModel
import sys

# %matplotlib inline

class BaseMultiTaskModel(BaseModel):
    """ Base class for all  Multi-Task estimators:
        * Multi-Task Logistic Regression model (MTLR)
        * Neural Multi-Task Logistic Regression model (N-MTLR)
    BaseMultiTaskModel shouldn't be used as is.
    The underlying model is written in PyTorch.

    The original Multi-Task model, a.k.a the Multi-Task Logistic Regression
    model (MTLR), was first introduced by  Chun-Nam Yu et al. in
    *Learning Patient-Specific Cancer Survival Distributions as a Sequence of
    Dependent Regressors*

    The Neural Multi-Task Logistic Regression model (N-MTLR) was developed
    by S. Fotso in the paper *Deep Neural Networks for Survival
    Analysis Based on a Multi-Task Framework*, allowing the use of
    Neural Networks within the original design.

    Parameters
    ----------
    * `structure`:  **list of dictionaries** --
        Provides the structure of the MLP built within the N-MTLR.

        ex: `structure = [ {'activation': 'ReLU', 'num_units': 128}, ]`.

        Each dictionary corresponds to a fully connected hidden layer:

        * `num_units` is the number of hidden units in this layer
        * `activation` is the activation function that will be used.
        The list of all available activation functions can be found :
            * Atan
            * BentIdentity
            * BipolarSigmoid
            * CosReLU
            * ELU
            * Gaussian
            * Hardtanh
            * Identity
            * InverseSqrt
            * LeakyReLU
            * LeCunTanh
            * LogLog
            * LogSigmoid
            * ReLU
            * SELU
            * Sigmoid
            * Sinc
            * SinReLU
            * Softmax
            * Softplus
            * Softsign
            * Swish
            * Tanh

        In case there are more than one dictionary,
        each hidden layer will be applied in the resulting MLP,
        using the order it is provided in the structure:
        ex: structure = [ {'activation': 'relu', 'num_units': 128},
                          {'activation': 'tanh', 'num_units': 128}, ]

    * `bins`: **int** *(default=100)* --
         Number of subdivisions of the time axis

    * `auto_scaler`: **boolean** *(default=True)* --
        Determines whether a sklearn scaler should be automatically applied
    """

    def __init__(self, structure, bins=100, auto_scaler=True):

        # Saving the attributes
        self.loss_values = []
        self.bins = bins
        self.structure = structure

        # Initializing the elements from BaseModel
        super(BaseMultiTaskModel, self).__init__(auto_scaler)

    def get_times(self, T, is_min_time_zero=True, extra_pct_time=0.1):
        """ Building the time axis (self.times) as well as the time intervals
            ( all the [ t(k-1), t(k) ) in the time axis.
        """

        # Setting the min_time and max_time
        max_time = max(T)
        if is_min_time_zero:
            min_time = 0.
        else:
            min_time = min(T)

        # Setting optional extra percentage time
        if 0. <= extra_pct_time <= 1.:
            p = extra_pct_time
        else:
            raise Exception("extra_pct_time has to be between [0, 1].")

            # Building time points and time buckets
        self.times = np.linspace(min_time, max_time * (1. + p), self.bins)
        self.get_time_buckets()
        self.num_times = len(self.time_buckets)

    def compute_XY(self, X, T, E, is_min_time_zero, extra_pct_time):
        """ Given the survival_times, events and time_points vectors,
            it returns a ndarray of the encodings for all units
            such that:
                Y = [[0, 0, 1, 0, 0],  # unit experienced an event at t = 3
                     [0, 1, 0, 0, 0],  # unit experienced an event at t = 2
                     [0, 1, 1, 1, 1],] # unit was censored at t = 2
        """

        # building times axis
        self.get_times(T, is_min_time_zero, extra_pct_time)
        n_units = T.shape[0]

        # Initializing the output variable
        Y_cens, Y_uncens = [], []
        X_cens, X_uncens = [], []

        if isinstance(X, list):
            for input_ in X:
                X_cens.append([])
                X_uncens.append([])

        # Building the output variable
        for i, (t, e) in enumerate(zip(T, E)):
            y = np.zeros(self.num_times + 1)
            min_abs_value = [abs(a_j_1 - t) for (a_j_1, a_j) in self.time_buckets]
            index = np.argmin(min_abs_value)

            if e == 1:
                y[index] = 1.
                if isinstance(X, list):
                    for j, input_ in enumerate(X):
                        X_uncens[j].append(input_[i, :].tolist())
                else:
                    X_uncens.append(X[i, :].tolist())
                Y_uncens.append(y.tolist())

            else:
                y[(index):] = 1.

                if isinstance(X, list):
                    for j, input_ in enumerate(X):
                        X_cens[j].append(input_[i, :].tolist())
                else:
                    X_cens.append(X[i, :].tolist())
                Y_cens.append(y.tolist())

        # Transform into torch.Tensor
        if isinstance(X, list):
            for j, input_ in enumerate(X_cens):
                X_cens[j] = torch.FloatTensor(input_)
                if torch.cuda.is_available():
                    X_cens[j] = X_cens[j].cuda()
            for j, input_ in enumerate(X_uncens):
                X_uncens[j] = torch.FloatTensor(input_)
                if torch.cuda.is_available():
                    X_uncens[j] = X_uncens[j].cuda()
        else:
            X_cens = torch.FloatTensor(X_cens)
            X_uncens = torch.FloatTensor(X_uncens)
            if torch.cuda.is_available():
                X_cens = X_cens.cuda()
                X_uncens = X_uncens.cuda()

        Y_cens = torch.FloatTensor(Y_cens)
        Y_uncens = torch.FloatTensor(Y_uncens)
        if torch.cuda.is_available():
            Y_cens = Y_cens.cuda()
            Y_uncens = Y_uncens.cuda()

        return X_cens, X_uncens, Y_cens, Y_uncens

    def loss_function(self, model, X_cens, X_uncens, Y_cens, Y_uncens,
                      Triangle, l2_reg, l2_smooth, min_clamp_value=1e-8,
                      max_clamp_value=torch.finfo(torch.float32).max-1):
        """ Computes the loss function of the any MTLR model.
            All the operations have been vectorized to ensure optimal speed
        """

        score_cens = model(X_cens)
        score_uncens = model(X_uncens)
        # Likelihood Calculations -- Uncensored
        temp = torch.clamp(torch.mm(score_uncens, Triangle), min=1e-8, max=88.5)
        phi_uncens = torch.clamp(torch.exp(temp), max=max_clamp_value)
        reduc_phi_uncens = torch.sum(phi_uncens * Y_uncens, dim=1)

        # Likelihood Calculations -- Censored
        temp = torch.clamp(torch.mm(score_cens, Triangle), min=1e-8, max=88.5)
        phi_cens = torch.clamp(torch.exp(temp), max=max_clamp_value)
        reduc_phi_cens = torch.sum(phi_cens * Y_cens, dim=1)

        # Likelihood Calculations -- Normalization
        temp = torch.clamp(torch.mm(score_uncens, Triangle), min=1e-8, max=88.5)
        z_uncens = torch.clamp(torch.exp(temp), max=max_clamp_value)
        reduc_z_uncens = torch.sum(z_uncens, dim=1)

        temp = torch.clamp(torch.mm(score_cens, Triangle), min=1e-8, max=88.5)
        z_cens = torch.clamp(torch.exp(temp), max=max_clamp_value)
        reduc_z_cens = torch.sum(z_cens, dim=1)

        reduc_phi_uncens = torch.clamp(reduc_phi_uncens, min=min_clamp_value, max=max_clamp_value)
        reduc_phi_cens = torch.clamp(reduc_phi_cens, min=min_clamp_value, max=max_clamp_value)
        reduc_z_uncens = torch.clamp(reduc_z_uncens, min=min_clamp_value, max=max_clamp_value)
        reduc_z_cens = torch.clamp(reduc_z_cens, min=min_clamp_value, max=max_clamp_value)

        # MTLR cost function
        loss = - (
                torch.sum(torch.log(reduc_phi_uncens))
                + torch.sum(torch.log(reduc_phi_cens))

                - torch.sum(torch.log(reduc_z_uncens))
                - torch.sum(torch.log(reduc_z_cens))
        )

        # print(f"loss_loss: {loss}")
        # Adding the regularized loss
        nb_set_parameters = len(list(model.parameters()))
        for i, w in enumerate(model.parameters()):
            loss += l2_reg * torch.sum(w * w) / 2.

            if i >= nb_set_parameters - 2:
                loss += l2_smooth * norm_diff(w)

        return loss

    def fit(self, X, T, E, init_method='glorot_uniform', optimizer='adam',
            lr=1e-4, num_epochs=1000, dropout=0.2, l2_reg=1e-2,
            l2_smooth=1e-2, batch_normalization=False, bn_and_dropout=False,
            verbose=True, extra_pct_time=0.1, is_min_time_zero=True, max_norm=1.0,
            min_clamp_value=1e-8, max_clamp_value=torch.finfo(torch.float32).max -1):
        """ Fit the estimator based on the given parameters.

        Parameters:
        -----------
        * `X` : **array-like**, *shape=(n_samples, n_features)* --
            The input samples.

        * `T` : **array-like** --
            The target values describing when the event of interest or censoring
            occurred.

        * `E` : **array-like** --
            The values that indicate if the event of interest occurred i.e.:
            E[i]=1 corresponds to an event, and E[i] = 0 means censoring,
            for all i.

        * `init_method` : **str** *(default = 'glorot_uniform')* --
            Initialization method to use. Here are the possible options:

            * `glorot_uniform`: Glorot/Xavier uniform initializer
            * `he_uniform`: He uniform variance scaling initializer
            * `uniform`: Initializing tensors with uniform (-1, 1) distribution
            * `glorot_normal`: Glorot normal initializer,
            * `he_normal`: He normal initializer.
            * `normal`: Initializing tensors with standard normal distribution
            * `ones`: Initializing tensors to 1
            * `zeros`: Initializing tensors to 0
            * `orthogonal`: Initializing tensors with a orthogonal matrix,

        * `optimizer`:  **str** *(default = 'adam')* --
            iterative method for optimizing a differentiable objective function.
            Here are the possible options:

            - `adadelta`
            - `adagrad`
            - `adam`
            - `adamax`
            - `rmsprop`
            - `sparseadam`
            - `sgd`

        * `lr`: **float** *(default=1e-4)* --
            learning rate used in the optimization

        * `num_epochs`: **int** *(default=1000)* --
            The number of iterations in the optimization

        * `dropout`: **float** *(default=0.5)* --
            Randomly sets a fraction rate of input units to 0
            at each update during training time, which helps prevent overfitting.

        * `l2_reg`: **float** *(default=1e-4)* --
            L2 regularization parameter for the model coefficients

        * `l2_smooth`: **float** *(default=1e-4)* --
            Second L2 regularizer that ensures the parameters vary smoothly
            across consecutive time points.

        * `batch_normalization`: **bool** *(default=True)* --
            Applying Batch Normalization or not

        * `bn_and_dropout`: **bool** *(default=False)* --
            Applying Batch Normalization and Dropout at the same time

        * `display_loss`: **bool** *(default=True)* --
            Whether or not showing the loss function values at each update

        * `verbose`: **bool** *(default=True)* --
            Whether or not producing detailed logging about the modeling

        * `extra_pct_time`: **float** *(default=0.1)* --
            Providing an extra fraction of time in the time axis

        * `is_min_time_zero`: **bool** *(default=True)* --
            Whether the the time axis starts at 0

        * `max_norm`: **float** *(default=1.0)* --
            Max l2 norm for gradient clipping

        **Returns:**

        * self : object

        Example:
        --------

        #### 1 - Importing packages
        import numpy as np
        import pandas as pd
        from matplotlib import pyplot as plt
        from sklearn.model_selection import train_test_split
        from pysurvival.models.simulations import SimulationModel
        from pysurvival.models.multi_task import LinearMultiTaskModel
        from pysurvival.utils.metrics import concordance_index
        #%matplotlib inline  # To use with Jupyter notebooks

        #### 2 - Generating the dataset from a Weibull parametric model
        # Initializing the simulation model
        sim = SimulationModel( survival_distribution = 'Weibull',
                               risk_type = 'linear',
                               censored_parameter = 10.0,
                               alpha = .01, beta = 3.0 )

        # Generating N random samples
        N = 1000
        dataset = sim.generate_data(num_samples = N, num_features = 3)

        # Showing a few data-points
        time_column = 'time'
        event_column = 'event'
        dataset.head(2)

        #### 3 - Creating the modeling dataset
        # Defining the features
        features = sim.features

        # Building training and testing sets #
        index_train, index_test = train_test_split( range(N), test_size = 0.2)
        data_train = dataset.loc[index_train].reset_index( drop = True )
        data_test  = dataset.loc[index_test].reset_index( drop = True )

        # Creating the X, T and E input
        X_train, X_test = data_train[features], data_test[features]
        T_train, T_test = data_train['time'].values, data_test['time'].values
        E_train, E_test = data_train['event'].values, data_test['event'].values

        #### 4 - Initializing a MTLR model and fitting the data.
        # Building a Linear model
        mtlr = LinearMultiTaskModel(bins=50)
        mtlr.fit(X_train, T_train, E_train, lr=5e-3, init_method='orthogonal')

        # Building a Neural MTLR
        # structure = [ {'activation': 'Swish', 'num_units': 150},  ]
        # mtlr = NeuralMultiTaskModel(structure=structure, bins=150)
        # mtlr.fit(X_train, T_train, E_train, lr=5e-3, init_method='adam')

        #### 5 - Cross Validation / Model Performances
        c_index = concordance_index(mtlr, X_test, T_test, E_test) #0.95
        print('C-index: {:.2f}'.format(c_index))

        """

        # Checking data format (i.e.: transforming into numpy array)
        X, T, E = utils.check_data(X, T, E)

        input_shape = []
        # Extracting data parameters
        if isinstance(X, list):
            nb_inputs = len(X)
            for data in X:
                nb_units, num_vars = data.shape
                input_shape.append(num_vars)
            # Scaling data
            if self.auto_scaler:
                for index, data in enumerate(X):
                    X[index] = self.scaler.fit_transform(data)
        else:
            nb_inputs = 1
            nb_units, self.num_vars = X.shape
            input_shape.append(self.num_vars)
            # Scaling data
            if self.auto_scaler:
                X = self.scaler.fit_transform(X)

        # Building the time axis, time buckets and output Y
        X_cens, X_uncens, Y_cens, Y_uncens \
            = self.compute_XY(X, T, E, is_min_time_zero, extra_pct_time)

        # Initializing the model
        model = nn.NeuralNet(input_shape, self.num_times, self.structure,
                             init_method, dropout, batch_normalization,
                             bn_and_dropout)

        # Creating the Triangular matrix
        Triangle = np.tri(self.num_times, self.num_times + 1, dtype=np.float32)
        Triangle = torch.FloatTensor(Triangle)

        if torch.cuda.is_available():
            model = model.cuda()
            Triangle = Triangle.cuda()

        # Performing order 1 optimization
        model, loss_values = opt.optimize(self.loss_function, model, optimizer,
                                          lr, num_epochs, verbose, X_cens=X_cens, X_uncens=X_uncens,
                                          Y_cens=Y_cens, Y_uncens=Y_uncens, Triangle=Triangle,
                                          l2_reg=l2_reg, l2_smooth=l2_smooth, max_norm=max_norm,
                                          min_clamp_value=min_clamp_value, max_clamp_value=max_clamp_value)

        # Saving attributes
        self.model = model.eval()
        self.loss_values = loss_values

        return self

    def predict(self, x, t=None):
        """ Predicting the hazard, density and survival functions

        Parameters:
        ----------
        * `x` : **array-like** *shape=(n_samples, n_features)* --
            array-like representing the datapoints.
            x should not be standardized before, the model
            will take care of it

        * `t`: **double** *(default=None)* --
             time at which the prediction should be performed.
             If None, then return the function for all available t.
        """

        # Convert x into the right format
        x = utils.check_data(x)

        # Scaling the data
        if self.auto_scaler:
            if x.ndim == 1:
                x = self.scaler.transform(x.reshape(1, -1))
            elif x.ndim == 2:
                x = self.scaler.transform(x)
            x = torch.FloatTensor(x)
            if torch.cuda.is_available():
                x = x.cuda()
        else:
            # Ensuring x has 2 dimensions
            if isinstance(x, list):
                for index in range(len(x)):
                    if x[index].ndim == 1:
                        x[index] = np.reshape(x[index], (1, -1))
                    # Transforming into pytorch objects
                    x[index] = torch.FloatTensor(x[index])
                    if torch.cuda.is_available():
                        x[index] = x[index].cuda()
            else:
                if x.ndim == 1:
                    x = np.reshape(x, (1, -1))
                # Transforming into pytorch objects
                x = torch.FloatTensor(x)
                if torch.cuda.is_available():
                    x = x.cuda()

        # Predicting using linear/nonlinear function
        score_torch = self.model(x)
        score = score_torch.data.cpu().numpy()

        # Cretaing the time triangles
        Triangle1 = np.tri(self.num_times, self.num_times + 1)
        Triangle2 = np.tri(self.num_times + 1, self.num_times + 1)

        # Calculating the score, density, hazard and Survival
        phi = np.exp(np.dot(score, Triangle1))
        div = np.repeat(np.sum(phi, 1).reshape(-1, 1), phi.shape[1], axis=1)
        density = (phi / div)
        Survival = np.dot(density, Triangle2)
        hazard = density[:, :-1] / Survival[:, 1:]

        # Returning the full functions of just one time point
        if t is None:
            return hazard, density, Survival
        else:
            min_abs_value = [abs(a_j_1 - t) for (a_j_1, a_j) in self.time_buckets]
            index = np.argmin(min_abs_value)
            return hazard[:, index], density[:, index], Survival[:, index]

    def predict_risk(self, x, use_log=False):
        """ Computing the risk score

        Parameters:
        -----------
        * `x` : **array-like** *shape=(n_samples, n_features)* --
            array-like representing the datapoints.
            x should not be standardized before, the model
            will take care of it

        * `use_log`: **bool** *(default=True)* --
            Applies the log function to the risk values

        """

        risk = super(BaseMultiTaskModel, self).predict_risk(x)
        if use_log:
            return np.log(risk)
        else:
            return risk

class LinearMultiTaskModel(BaseMultiTaskModel):
    """ LinearMultiTaskModel is the original Multi-Task model,
        a.k.a the Multi-Task Logistic Regression model (MTLR).
        It was first introduced by  Chun-Nam Yu et al. in
        Learning Patient-Specific Cancer Survival Distributions
        as a Sequence of Dependent Regressors

        Reference:
        ----------
            * http://www.cs.cornell.edu/~cnyu/papers/nips11_survival.pdf

        Parameters:
        ----------
            * bins: int
                Number of subdivisions of the time axis

            * auto_scaler: boolean (default=True)
                Determines whether a sklearn scaler should be automatically
                applied

    """

    def __init__(self, bins=100, auto_scaler=True):
        super(LinearMultiTaskModel, self).__init__(
            structure=None, bins=bins, auto_scaler=auto_scaler)

    def fit(self, X, T, E, init_method='glorot_uniform', optimizer='adam',
            lr=1e-4, num_epochs=1000, l2_reg=1e-2, l2_smooth=1e-2,
            verbose=True, extra_pct_time=0.1, is_min_time_zero=True, max_norm=1.0,
            min_clamp_value=1e-8, max_clamp_value=torch.finfo(torch.float32).max-1):
        super(LinearMultiTaskModel, self).fit(X=X, T=T, E=E,
                                              init_method=init_method, optimizer=optimizer,
                                              lr=lr, num_epochs=num_epochs, dropout=None, l2_reg=l2_reg,
                                              l2_smooth=l2_smooth, batch_normalization=False,
                                              bn_and_dropout=False, verbose=verbose,
                                              extra_pct_time=extra_pct_time, is_min_time_zero=is_min_time_zero,
                                              max_norm=max_norm,
                                              min_clamp_value=min_clamp_value, max_clamp_value=max_clamp_value)

        return self

class NeuralMultiTaskModel(BaseMultiTaskModel):
    """ NeuralMultiTaskModel is the Neural Multi-Task Logistic Regression
        model (N-MTLR) was developed by Fotso S. in
        Deep Neural Networks for Survival Analysis Based on a
        Multi-Task Framework,
        allowing the use of Neural Networks within the original design.

    Reference:
    ----------
    * https://arxiv.org/pdf/1801.05512

    Parameters:
    ----------
    * `structure`:  **list of dictionaries** --
        Provides the structure of the MLP built within the N-MTLR.

        ex: `structure = [ {'activation': 'ReLU', 'num_units': 128}, ]`.

        Each dictionary corresponds to a fully connected hidden layer:

        * `units` is the number of hidden units in this layer
        * `activation` is the activation function that will be used.
        The list of all available activation functions can be found :
            * Atan
            * BentIdentity
            * BipolarSigmoid
            * CosReLU
            * ELU
            * Gaussian
            * Hardtanh
            * Identity
            * InverseSqrt
            * LeakyReLU
            * LeCunTanh
            * LogLog
            * LogSigmoid
            * ReLU
            * SELU
            * Sigmoid
            * Sinc
            * SinReLU
            * Softmax
            * Softplus
            * Softsign
            * Swish
            * Tanh

        In case there are more than one dictionary,
        each hidden layer will be applied in the resulting MLP,
        using the order it is provided in the structure:
        ex: structure = [ {'activation': 'relu', 'num_units': 128},
                          {'activation': 'tanh', 'num_units': 128}, ]

    * `bins`: **int** *(default=100)* --
         Number of subdivisions of the time axis

    * `auto_scaler`: **boolean** *(default=True)* --
        Determines whether a sklearn scaler should be automatically applied
    """

    def __init__(self, structure, bins=100, auto_scaler=True):

        # Checking the validity of structure
        structure = nn.check_mlp_structure(structure)

        # print(structure)

        # Initializing the instance
        super(NeuralMultiTaskModel, self).__init__(
            structure=structure, bins=bins, auto_scaler=auto_scaler)

    def __repr__(self):
        """ Representing the class object """

        if self.structure is None:
            super(NeuralMultiTaskModel, self).__repr__()
            return self.name

        else:
            S = len(self.structure)
            self.name = self.__class__.__name__
            empty = len(self.name)
            self.name += '( '
            for i, s in enumerate(self.structure):
                if isinstance(s, list):
                    for s_ in s:
                        n = 'Layer({}): '.format(i + 1)
                        activation = nn.activation_function(s_['activation'],
                                                            return_text=True)
                        n += 'activation = {}, '.format(s_['activation'])
                        if 'num_units' in s_.keys():
                            n += 'units = {} '.format(s_['num_units'])

                        if i != S - 1:
                            self.name += n + '; \n'
                            self.name += empty * ' ' + '  '
                        else:
                            self.name += n
                    self.name = self.name + ')'
                else:
                    n = 'Layer({}): '.format(i + 1)
                    activation = nn.activation_function(s['activation'],
                                                        return_text=True)
                    n += 'activation = {}, '.format(s['activation'])
                    if 'num_units' in s.keys():
                        n += 'units = {} '.format(s['num_units'])

                    if i != S - 1:
                        self.name += n + '; \n'
                        self.name += empty * ' ' + '  '
                    else:
                        self.name += n
                self.name = self.name + ')'

            return self.name

def norm_diff(W):
    """ Special norm function for the last layer of the MTLR """
    dims = len(W.shape)
    if dims == 1:
        diff = W[1:] - W[:-1]
    elif dims == 2:
        diff = W[1:, :] - W[:-1, :]
    return torch.sum(diff * diff)


text:
import sympy as sp

x, y, psi, v, dpsi, T = sp.symbols('x y psi v dpsi T')

state = sp.Matrix([x,
                   y,
                   psi,
                   v,
                   dpsi])

F = sp.Matrix([[x + (v/dpsi) * (sp.sin(T * dpsi + psi) - sp.sin(psi))],
               [y + (v/dpsi) * (sp.cos(psi) - sp.cos(T * dpsi + psi))],
               [T * dpsi + psi],
               [v],
               [dpsi]])

jF = F.jacobian(state)

# print(F[0], '\n', F[1], '\n', F[2], '\n', F[3], '\n',F[4])

# for i in range(0, 5):
#     print(jF[i, 0], ' ', jF[i, 1], ' ', jF[i, 2],
#           ' ', jF[i, 3], ' ', jF[i, 4], ' ')

h = sp.Matrix([sp.sqrt(x), sp.sqrt(y)])
jH = h.jacobian(state)

print(jF)

ext: py

max_issues_repo_name: raghuramshankar/kalman-filter

max_forks_repo_name: raghuramshankar/kalman-filter

reason: import sympy


text:
#------------------------------------------------------------------------------
#
#   DQN Demo
#
#   Author : Igor Janos
#
#------------------------------------------------------------------------------

import numpy as np
import torch
import random
import math

from skimage.transform import resize

#------------------------------------------------------------------------------
#
#   RandomPolicy class
#
#------------------------------------------------------------------------------

class RandomPolicy:
    def __init__(self, nActions):
        self.nActions = nActions

    def selectAction(self, state):
        return torch.tensor([[random.randrange(self.nActions)]], dtype=torch.long)

    def reset(self):
        pass

#------------------------------------------------------------------------------
#
#   GreedyPolicy class
#
#------------------------------------------------------------------------------

class GreedyPolicy:

    def __init__(self, nActions, netQ):
        self.nActions = nActions
        self.netQ = netQ

    def selectAction(self, state):
        with torch.no_grad():
            actionDistribution = self.netQ(state)
            act = actionDistribution.max(1)     # Highest value
            act = act[1].view(1, 1)             # Index of highest value
            return act

    def reset(self):
        pass

#------------------------------------------------------------------------------
#
#   EpsilonGreedyPolicy class
#
#------------------------------------------------------------------------------

class EpsilonGreedyPolicy:

    def __init__(self, nActions, netQ, epsilon=0.1, epsilonUpper=0.1, decay=0.0):
        self.nActions = nActions
        self.netQ = netQ
        self.epsilon = epsilon
        self.epsilonUpper = epsilonUpper
        self.decay = decay
        self.steps = 0
        self.currentEpsilon = epsilon

    def selectAction(self, state):

        # Casom klesa ?
        if (self.decay > 0.0):
            epsilon = self.epsilon + (self.epsilonUpper - self.epsilon)*math.exp(-1.0 * self.steps / self.decay)
        else:
            epsilon = self.epsilon

        self.steps += 1
        self.currentEpsilon = epsilon

        sample = random.random()
        if (sample > epsilon):
            with torch.no_grad():
                actionDistribution = self.netQ(state)
                act = actionDistribution.max(1)     # Highest value
                act = act[1].view(1, 1)             # Index of highest value
            return act
        else:
            return torch.tensor([[random.randrange(self.nActions)]], dtype=torch.long)

#------------------------------------------------------------------------------
#
#   EnvState class
#
#------------------------------------------------------------------------------

class EnvState:
    def __init__(self, nFrames, shape):
        self.nFrames = nFrames
        self.shape = shape
        self.frames = []

    def reset(self):
        self.frames = []

    def push(self, frame, lastFrame):

        # vypadne (h,w)
        frame = preprocessFrame(frame, shape=self.shape)
        lastFrame = preprocessFrame(lastFrame, shape=self.shape)
        frame = np.maximum(frame, lastFrame)

        # Vyhodime prvy
        if (len(self.frames) > 0): self.frames.pop(0)

        # Dokladame na koniec - uchovavame ako uint8 v pamati
        while (len(self.frames) < self.nFrames):
            self.frames.append(torch.from_numpy(frame).unsqueeze(0))

        # Vratime von - 4 posledne frejmy
        frames = [
            self.frames[3],
            self._difFrame(self.frames[3], self.frames[2]),
            self._difFrame(self.frames[2], self.frames[1]),
            self._difFrame(self.frames[1], self.frames[0])
        ]

        result = torch.cat(frames, dim=0)
        result = result.unsqueeze(0)
        result = result.to(torch.uint8)
        return result

    def _difFrame(self, f1, f2):
        return (128 + 0.5*(f1.float()-f2.float())).to(torch.uint8)

def preprocessFrame(frame, shape=(84,84)):

    '''
        Frame - (h, w, ch)
        Skonvertujeme do GrayScale, a potom orezeme na
        192 x 160 rozmer
    '''

    h, w, ch = frame.shape

    # Odrezeme nadbytocnosti
    nLines = 192
    nFirst = (h - nLines) // 2
    frame = frame[nFirst:(nFirst+nLines),:,:]

    # Skonvertujeme na ciernobiely obrazok <0;1>
    frame = np.max(frame, axis=2)
    frame = frame.astype(np.float32) / 255.0
    frame = resize(frame, shape, anti_aliasing=False)
    frame = ((frame > 0) * 255.0).astype(np.uint8)

    return frame

def visualize(state):

    # Skonvertujeme stav do podoby vhodnej na vizualizaciu
    bs, ch, h, w = state.shape

    frames = [ state[0,0], state[0,1], state[0,2], state[0,3] ]
    singleFrame = torch.cat(frames, dim=1)
    return singleFrame.cpu().detach().numpy()

reason: import numpy


text:
#!/usr/bin/env python
#
# Copyright (C) 2018 Intel Corporation
#
# SPDX-License-Identifier: MIT

from __future__ import absolute_import, division, print_function

import argparse
import glog as log
import numpy as np
import os.path as osp
import json
import cv2
import sys
from lxml import etree
from tqdm import tqdm
from skimage import measure
from pycocotools import mask as mask_util
from pycocotools import coco as coco_loader

def parse_args():
    """Parse arguments of command line"""
    parser = argparse.ArgumentParser(
        description='Convert CVAT annotation with instance segmentation ''to COCO representation'
    )
    parser.add_argument(
        '--cvat-xml', required=True,
        help='input file with CVAT annotation in *.xml format'
    )
    parser.add_argument(
        '--output', default=None,
        help='output annotation file. If not defined, the output file name will be created from input file name'
    )
    parser.add_argument(
        '--image-dir', required=True,
        help='directory with images from annotation'
    )
    parser.add_argument(
        '--labels', default=None,
        help='path to file with labels'
    )
    parser.add_argument(
        '--draw', default=None,
        help='directory to save images with its segments. By default is disabled'
    )
    parser.add_argument(
        '--draw_labels', action='store_true',
        help='insert in output images labels of objects. By default is false'
    )
    parser.add_argument(
        '--use_background_label', action='store_true',
        help='insert in output annotation objects with label \'background\'. By default is false'
    )
    parser.add_argument(
        '--polygon-area-threshold', type=int, default=1,
        help='polygons with area less than this value will be ignored. By default set to 1'
    )
    return parser.parse_args()

def mask_to_polygon(mask, tolerance=1.0, area_threshold=1):
    """Convert object's mask to polygon [[x1,y1, x2,y2 ...], [...]]
    Args:
        mask: object's mask presented as 2D array of 0 and 1
        tolerance: maximum distance from original points of polygon to approximated
        area_threshold: if area of a polygon is less than this value, remove this small object
    """
    polygons = []
    # pad mask with 0 around borders
    padded_mask = np.pad(mask, pad_width=1, mode='constant', constant_values=0)
    contours = measure.find_contours(padded_mask, 0.5)
    # Fix coordinates after padding
    contours = np.subtract(contours, 1)
    for contour in contours:
        if not np.array_equal(contour[0], contour[-1]):
            contour = np.vstack((contour, contour[0]))
        contour = measure.approximate_polygon(contour, tolerance)
        if len(contour) > 2:
            contour = np.flip(contour, axis=1)
            reshaped_contour = []
            for xy in contour:
                reshaped_contour.append(xy[0])
                reshaped_contour.append(xy[1])
            for i in range(0, len(reshaped_contour)):
                if reshaped_contour[i] < 0:
                    reshaped_contour[i] = 0
            # Check if area of a polygon is enough
            rle = mask_util.frPyObjects([reshaped_contour], mask.shape[0], mask.shape[1])
            area = mask_util.area(rle)
            if sum(area) > area_threshold:
                polygons.append(reshaped_contour)
    return polygons

def draw_polygons(polygons, img_name, input_dir, output_dir, draw_labels):
    """Draw on image contours of its objects and save
    Args:
        polygons: all objects on image represented as 2D array of objects' contours
        img_name: name of image file
        input_dir: path to directory with images from annotation
        output_dir: directory to save images
    """
    name = osp.basename(img_name)
    input_file = osp.join(input_dir, name)
    output_file = osp.join(output_dir, name)
    img = cv2.imread(input_file)
    yellow = (0, 255, 255)
    red = (0, 0, 255)
    for poly in polygons:
        label = poly['label']
        _, bbox = polygon_area_and_bbox(poly['points'], img.shape[0], img.shape[1])
        for j in range(0, len(poly['points'])):
            i = 0
            points = []
            while i < len(poly['points'][j]):
                x = int(poly['points'][j][i])
                y = int(poly['points'][j][i + 1])
                points.append([x, y])
                i += 2
            bbox = [int(value) for value in bbox]
            img = cv2.polylines(img, np.int32([points]), True, yellow, 1)
            img = cv2.rectangle(img, (bbox[0], bbox[1]), (bbox[0] + bbox[2], bbox[1] + bbox[3]), red, 2)
            if draw_labels:
                x = bbox[0] + bbox[2] // 4
                y = bbox[1] + bbox[3] // 2
                cv2.putText(img, label, (x, y), cv2.FONT_HERSHEY_COMPLEX_SMALL, 1, red, 1)
    cv2.imwrite(output_file, img)

def fix_segments_intersections(polygons, height, width, img_name, use_background_label,
                               threshold=0.0, ratio_tolerance=0.001, area_threshold=1):
    """Find all intersected regions and crop contour for back object by objects which
        are in front of the first one. It is related to a specialty of segmentation
        in CVAT annotation. Intersection is calculated via function 'iou' from cocoapi
    Args:
        polygons: all objects on image represented as 2D array of objects' contours
        height: height of image
        width: width of image
        img_name: name of image file
        threshold: threshold of intersection over union of two objects.
            By default is set to 0 and processes any two intersected objects
        ratio_tolerance: used for situation when one object is fully or almost fully
            inside another one and we don't want make "hole" in one of objects
    """
    converted_polygons = []
    empty_polygon = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    # Convert points of polygons from string to coco's array.
    # All polygons must be sorted in order from bottom to top
    for polygon in polygons:
        label = polygon['label']
        points = polygon['points'].split(';')
        new_polygon = []
        for xy in points:
            x = float(xy.split(',')[0])
            y = float(xy.split(',')[1])
            new_polygon.append(x)
            new_polygon.append(y)
        converted_polygons.append({'label': label, 'points': new_polygon})

    for i in range(0, len(converted_polygons)):
        rle_bottom = mask_util.frPyObjects([converted_polygons[i]['points']], height, width)
        segment_overlapped = False
        for j in range(i + 1, len(converted_polygons)):
            rle_top = mask_util.frPyObjects([converted_polygons[j]['points']], height, width)
            iou = mask_util.iou(rle_bottom, rle_top, [0, 0])
            area_top = sum(mask_util.area(rle_top))
            area_bottom = sum(mask_util.area(rle_bottom))
            if area_bottom == 0:
                continue
            area_ratio = area_top / area_bottom
            sum_iou = sum(iou)

            # If segment is fully inside another one, save this segment as is
            if area_ratio - ratio_tolerance < sum_iou[0] < area_ratio + ratio_tolerance:
                continue
            # Check situation when bottom segment is fully inside top.
            # It means that in annotation is mistake. Save this segment as is
            if 1 / area_ratio - ratio_tolerance < sum_iou[0] < 1 / area_ratio + ratio_tolerance:
                continue

            if sum_iou[0] > threshold:
                segment_overlapped = True
                bottom_mask = np.array(mask_util.decode(rle_bottom), dtype=np.uint8)
                top_mask = np.array(mask_util.decode(rle_top), dtype=np.uint8)

                bottom_mask = np.subtract(bottom_mask, top_mask)
                bottom_mask[bottom_mask > 1] = 0

                bottom_mask = np.sum(bottom_mask, axis=2)
                bottom_mask = np.array(bottom_mask > 0, dtype=np.uint8)
                converted_polygons[i]['points'] = mask_to_polygon(bottom_mask, area_threshold=area_threshold)
                # If some segment is empty, do small fix to avoid error in cocoapi function
                if len(converted_polygons[i]['points']) == 0:
                    converted_polygons[i]['points'] = [empty_polygon]
                rle_bottom = mask_util.frPyObjects(converted_polygons[i]['points'], height, width)
        if not segment_overlapped:
            converted_polygons[i]['points'] = [converted_polygons[i]['points']]

    output_polygons = []
    for i in range(0, len(converted_polygons)):
        if not use_background_label and converted_polygons[i]['label'] == 'background':
            continue
        poly_len = len(converted_polygons[i]['points'])
        if poly_len == 0 or converted_polygons[i]['points'] == [empty_polygon]:
            log.warning('Image <{}> has an empty polygon with label <{}>. '
                        'Perhaps there is a mistake in annotation'.
                        format(img_name, converted_polygons[i]['label']))
        else:
            output_polygons.append(converted_polygons[i])

    return output_polygons

def polygon_area_and_bbox(polygon, height, width):
    """Calculate area of object's polygon and bounding box around it
    Args:
        polygon: objects contour represented as 2D array
        height: height of object's region (use full image)
        width: width of object's region (use full image)
    """
    rle = mask_util.frPyObjects(polygon, height, width)
    area = mask_util.area(rle)
    bbox = mask_util.toBbox(rle)
    bbox = [min(bbox[:, 0]),
            min(bbox[:, 1]),
            max(bbox[:, 0] + bbox[:, 2]) - min(bbox[:, 0]),
            max(bbox[:, 1] + bbox[:, 3]) - min(bbox[:, 1])]
    return area, bbox

def insert_license_data(result_annotation):
    """Fill license fields in annotation by blank data
    Args:
        result_annotation: output annotation in COCO representation
    """
    result_annotation['licenses'].append({
        'name': '',
        'id': 0,
        'url': ''
    })

def insert_info_data(xml_root, result_annotation):
    """Fill available information of annotation
    Args:
        xml_root: root for xml parser
        result_annotation: output annotation in COCO representation
    """
    log.info('Reading information data...')
    version = ''
    date = ''
    description = ''
    year = ''
    for child in xml_root:
        if child.tag == 'version':
            version = child.text
        if child.tag == 'meta':
            for task in child:
                for entry in task:
                    if entry.tag == 'name':
                        description = entry.text
                    if entry.tag == 'created':
                        date = entry.text
    date = date.split(' ')[0]
    year = date.split('-')[0]
    result_annotation['info'] = {
        'contributor': '',
        'date_created': date,
        'description': description,
        'url': '',
        'version': version,
        'year': year
    }
    log.info('Found the next information data: {}'.format(result_annotation['info']))

def insert_categories_data(xml_root, use_background_label, result_annotation, labels_file=None):
    """Get labels from input annotation and fill categories field in output annotation
    Args:
        xml_root: root for xml parser
        use_background_label: key to enable using label background
        result_annotation: output annotation in COCO representation
        labels_file: path to file with labels names.
                     If not defined, parse annotation to get labels names
    """
    def get_categories(names, bg_found, use_background_label, sort=False):
        bg_used = False
        category_map = {}
        categories = []
        # Sort labels by its names to make the same order of ids for different annotations
        if sort:
            names.sort()
        # Always use id = 0 for background
        if bg_found and use_background_label:
            category_map['background'] = 0
            bg_used = True
        cat_id = 1
        # Define id for all labels beginning from 1
        for name in names:
            if name == 'background':
                continue
            category_map[name] = cat_id
            categories.append({'id': cat_id, 'name': name, 'supercategory': ''})
            cat_id += 1
        return category_map, categories, bg_used

    categories = []
    category_map = {}
    label_names = []
    bg_found = False
    bg_used = False

    if labels_file is None:
        log.info('Reading labels from annotation...')
        for label in xml_root.iter('label'):
            for name in label.findall("./name"):
                if name.text == 'background':
                    bg_found = True
                    continue
                label_names.append(name.text)
        if len(label_names) == 0:
            log.info('Labels in annotation were not found. Please use \'--labels\' argument to define file with labels.')
        else:
            category_map, categories, bg_used = get_categories(label_names, bg_found, use_background_label, sort=True)
    else:
        log.info('Parsing labels from file <{}>...'.format(labels_file))
        with open(labels_file, 'r') as file:
            string = '  '
            while string != '' and string != '\n':
                string = file.readline()
                labels = string.split(' ')
                for label in labels:
                    if label == '\n' or label == '':
                        continue
                    label = label.replace('\n', '')
                    if label == 'background':
                        bg_found = True
                        continue
                    label_names.append(label)
            category_map, categories, bg_used = get_categories(label_names, bg_found, use_background_label)

    if len(categories) == 0:
        raise ValueError('Categories list is empty. Something wrong.')

    result_annotation['categories'] = categories
    log.info('Found the next labels: {}'.format(category_map))
    if bg_found and not bg_used:
        log.warning('Label <background> was found but not used. '
                    'To enable it should use command line argument [--use_background_label]')
    return category_map

def insert_image_data(image, result_annotation):
    """Get data from input annotation for image and fill fields for this image in output annotation
    Args:
        image: dictionary with data for image from original annotation
        path_to_images: path to directory with images from annotation
        result_annotation: output annotation in COCO representation
    """
    new_img = {}
    new_img['coco_url'] = ''
    new_img['date_captured'] = ''
    new_img['flickr_url'] = ''
    new_img['license'] = 0
    new_img['id'] = image['id']
    new_img['file_name'] = osp.basename(image['name'])
    new_img['height'] = int(image['height'])
    new_img['width'] = int(image['width'])
    result_annotation['images'].append(new_img)

def insert_annotation_data(image, category_map, segm_id, object, img_dims, result_annotation):
    """Get data from input annotation for object and fill fields for this object in output annotation
    Args:
        image: dictionary with data for image from input CVAT annotation
        category_map: map for categories represented in the annotation {name: id}
        segm_id: identificator of current object
        object: includes data for the object [label, polygon]
        img_dims: dimensions of image [height, width]
        result_annotation: output annotation in COCO representation
    """
    new_anno = {}
    new_anno['category_id'] = category_map[object['label']]
    new_anno['id'] = segm_id
    new_anno['image_id'] = image['id']
    new_anno['iscrowd'] = 0
    new_anno['segmentation'] = object['points']
    area, bbox = polygon_area_and_bbox(object['points'], img_dims[0], img_dims[1])
    new_anno['area'] = float(np.sum(area))
    new_anno['bbox'] = bbox
    result_annotation['annotations'].append(new_anno)

def main():
    args = parse_args()
    xml_file_name = args.cvat_xml
    if args.output is not None:
        output_file_name = args.output
    else:
        output_file_name = args.cvat_xml.split('.xml')[0] + '.json'
        log.info('Output file name set to: {}'.format(output_file_name))
    root = etree.parse(xml_file_name).getroot()

    if args.draw is not None:
        log.info('Draw key was enabled. Images will be saved in directory <{}>'.format(args.draw))

    result_annotation = {
        'licenses': [],
        'info': {},
        'categories': [],
        'images': [],
        'annotations': []
    }

    insert_license_data(result_annotation)
    insert_info_data(root, result_annotation)
    category_map = insert_categories_data(root, args.use_background_label, result_annotation, labels_file=args.labels)

    if len(category_map) == 0:
        sys.exit('Labels were not found. Be sure that annotation <{}> includes field <labels> or '
                 'annotation directory includes file <labels.txt>'.format(xml_file_name))

    segm_id = 0
    z_order_off_counter = 0
    # Parse original annotation
    for img in tqdm(root.iter('image'), desc='Processing images from ' + xml_file_name):
        image = {}
        for key, value in img.items():
            image[key] = value
        img_name = osp.join(args.image_dir, osp.basename(image['name']))
        if not osp.isfile(img_name):
            log.warning('Image <{}> is not available'.format(img_name))
        image['polygon'] = []
        z_order_on_counter = 0
        polygon_counter = 0
        for poly in img.iter('polygon'):
            polygon = {}
            for key, value in poly.items():
                polygon[key] = value
                if key == 'z_order':
                    z_order_on_counter += 1
            polygon_counter += 1
            image['polygon'].append(polygon)
        # If at least one of polygons on image does not have field 'z_order' do not sort them
        if z_order_on_counter == polygon_counter:
            image['polygon'].sort(key=lambda x: int(x['z_order']))
        else:
            z_order_off_counter += 1

        # Create new image
        image['id'] = int(image['id'])
        insert_image_data(image, result_annotation)
        height = result_annotation['images'][-1]['height']
        width = result_annotation['images'][-1]['width']
        image['polygon'] = fix_segments_intersections(image['polygon'], height, width,
                                                      image['name'], args.use_background_label,
                                                      area_threshold=args.polygon_area_threshold)

        # Create new annotation for this image
        for poly in image['polygon']:
            insert_annotation_data(image, category_map, segm_id, poly, [height, width], result_annotation)
            segm_id += 1

        # Draw contours of objects on image
        if args.draw is not None:
            draw_polygons(image['polygon'], image['name'], args.image_dir, args.draw, args.draw_labels)

    log.info('Processed images: {}'.format(len(result_annotation['images'])))
    log.info('Processed objects: {}'.format(len(result_annotation['annotations'])))
    if z_order_off_counter > 0:
        log.warning('Annotation does not have a field \'z_order\' for {} image(s). '
                    'Overlapped objects may be cropped incorrectly!'. format(z_order_off_counter))

    # Save created annotation
    log.info('Saving annotation...')
    with open(output_file_name, 'w') as outfile:
        json.dump(result_annotation, outfile)
    log.info('Annotation was saved in <{}> successfully'.format(output_file_name))

    # Try to load created annotation via cocoapi
    try:
        log.info('Trying to load annotation <{}> via cocoapi...'.format(output_file_name))
        coco_loader.COCO(output_file_name)
    except:
        raise
    else:
        log.info('Conversion <{}> --> <{}> has finished successfully!'.format(xml_file_name, output_file_name))

if __name__ == "__main__":
    main()

max_stars_repo_name: raunilillemets/cvat

max_issues_repo_name: raunilillemets/cvat

max_forks_repo_name: raunilillemets/cvat


text:
# -*- coding: utf-8 -*-
"""
Created on Mon Apr 22 09:33:04 2019

@author: if715029
"""

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from sklearn.preprocessing import PolynomialFeatures
from sklearn import linear_model
from sklearn.metrics import (accuracy_score,precision_score,recall_score)

#%% Importar datos
data = pd.read_csv('../data/ex2data2.txt',header=None)

X = data.iloc[:,0:2]
Y = data.iloc[:,2]

plt.scatter(X[0],X[1],c=Y)
plt.show()

#%% Preparar datos (crear polinomio 'X*')
ngrado = 10 # grado del polinomio
poly = PolynomialFeatures(ngrado)
Xa = poly.fit_transform(X) # X*

#%% Crear y entrenar el modelo de clasificación
#modelo = linear_model.LogisticRegression(C=1e20) # Con sobre-ajuste.
modelo = linear_model.LogisticRegression(C=10) # Sin sobre-ajuste.
modelo.fit(Xa,Y)

Yhat = modelo.predict(Xa)

#%% desempeño
accuracy_score(Y,Yhat)

#%% precisión
precision_score(Y,Yhat)

#%% recall
recall_score(Y,Yhat)

#%% ver modelo de prediccion
h = 0.01
x0min,x0max,x1min,x1max = X[0].min(),X[0].max(),X[1].min(),X[1].max()
xx,yy = np.meshgrid(np.arange(x0min,x0max,h),np.arange(x1min,x1max,h))

Xnew = pd.DataFrame(np.c_[xx.ravel(),yy.ravel()])

Xa_new = poly.fit_transform(Xnew)
Z = modelo.predict(Xa_new)
Z = Z.reshape(xx.shape)

plt.contour(xx,yy,Z)
plt.scatter(X[0],X[1],c=Y)
plt.show()

#%%
W = modelo.coef_
plt.bar(np.arange(len(W[0])),W[0])
plt.show()


text:
# Licensed under a 3-clause BSD style license - see LICENSE.rst
import pytest
import numpy as np
from numpy.testing import assert_allclose
import astropy.units as u
from astropy.coordinates import SkyCoord
from astropy.units import Unit
from gammapy.data import DataStore
from gammapy.irf import PSF3D, EffectiveAreaTable2D, PSFMap
from gammapy.makers.utils import make_map_exposure_true_energy, make_psf_map
from gammapy.maps import MapAxis, MapCoord, WcsGeom, RegionGeom, Map
from gammapy.utils.testing import requires_data, requires_dependency, mpl_plot_check

@pytest.fixture(scope="session")
def data_store():
    return DataStore.from_dir("$GAMMAPY_DATA/hess-dl3-dr1/")

def fake_psf3d(sigma=0.15 * u.deg, shape="gauss"):
    offset_axis = MapAxis.from_nodes([0, 1, 2, 3] * u.deg, name="offset")

    energy_axis_true = MapAxis.from_energy_bounds(
        "0.1 TeV", "10 TeV", nbin=4, name="energy_true"
    )

    rad = np.linspace(0, 1.0, 101) * u.deg
    rad_axis = MapAxis.from_edges(rad, name="rad")

    O, R, E = np.meshgrid(offset_axis.center, rad_axis.edges, energy_axis_true.center)

    Rmid = 0.5 * (R[:-1] + R[1:])
    if shape == "gauss":
        val = np.exp(-0.5 * Rmid ** 2 / sigma ** 2)
    else:
        val = Rmid < sigma

    drad = 2 * np.pi * (np.cos(R[:-1]) - np.cos(R[1:])) * u.Unit("sr")
    psf_value = val / ((val * drad).sum(0)[0])

    return PSF3D(
        axes=[energy_axis_true, offset_axis, rad_axis],
        data=psf_value.T.value,
        unit=psf_value.unit
    )

def fake_aeff2d(area=1e6 * u.m ** 2):
    energy_axis_true = MapAxis.from_energy_bounds(
        "0.1 TeV", "10 TeV", nbin=4, name="energy_true"
    )

    offset_axis = MapAxis.from_edges(
        [0.0, 1.0, 2.0, 3.0] * u.deg, name="offset"
    )

    return EffectiveAreaTable2D(
        axes=[energy_axis_true, offset_axis], data=area.value, unit=area.unit
     )

def test_make_psf_map():
    psf = fake_psf3d(0.3 * u.deg)

    pointing = SkyCoord(0, 0, unit="deg")
    energy_axis = MapAxis(
        nodes=[0.2, 0.7, 1.5, 2.0, 10.0], unit="TeV", name="energy_true"
    )
    rad_axis = MapAxis(nodes=np.linspace(0.0, 1.0, 51), unit="deg", name="rad")

    geom = WcsGeom.create(
        skydir=pointing, binsz=0.2, width=5, axes=[rad_axis, energy_axis]
    )

    psfmap = make_psf_map(psf, pointing, geom)

    assert psfmap.psf_map.geom.axes[0] == rad_axis
    assert psfmap.psf_map.geom.axes[1] == energy_axis
    assert psfmap.psf_map.unit == "deg-2"
    assert psfmap.psf_map.data.shape == (4, 50, 25, 25)

def make_test_psfmap(size, shape="gauss"):
    psf = fake_psf3d(size, shape)
    aeff2d = fake_aeff2d()

    pointing = SkyCoord(0, 0, unit="deg")
    energy_axis = MapAxis(
        nodes=[0.2, 0.7, 1.5, 2.0, 10.0], unit="TeV", name="energy_true"
    )
    rad_axis = MapAxis.from_edges(
        edges=np.linspace(0.0, 1, 101), unit="deg", name="rad"
    )

    geom = WcsGeom.create(
        skydir=pointing, binsz=0.2, width=5, axes=[rad_axis, energy_axis]
    )

    exposure_geom = geom.squash(axis_name="rad")

    exposure_map = make_map_exposure_true_energy(pointing, "1 h", aeff2d, exposure_geom)

    return make_psf_map(psf, pointing, geom, exposure_map)

def test_psf_map_containment_radius():
    psf_map = make_test_psfmap(0.15 * u.deg)
    psf = fake_psf3d(0.15 * u.deg)

    position = SkyCoord(0, 0, unit="deg")

    # Check that containment radius is consistent between psf_table and psf3d
    assert_allclose(
        psf_map.containment_radius(
            energy_true=1 * u.TeV, position=position, fraction=0.9
        ),
        psf.containment_radius(
            energy_true=1 * u.TeV, offset=0 * u.deg, fraction=0.9
        ),
        rtol=1e-2,
    )
    assert_allclose(
        psf_map.containment_radius(
            energy_true=1 * u.TeV, position=position, fraction=0.5
        ),
        psf.containment_radius(
            energy_true=1 * u.TeV, offset=0 * u.deg, fraction=0.5
        ),
        rtol=1e-2,
    )

def test_psf_map_containment():
    psf_map = make_test_psfmap(0.15 * u.deg)
    assert_allclose(
        psf_map.containment(rad=10 * u.deg, energy_true=[10] * u.TeV), 1
    )

def test_psfmap_to_psf_kernel():
    psfmap = make_test_psfmap(0.15 * u.deg)

    energy_axis = psfmap.psf_map.geom.axes[1]
    # create PSFKernel
    kern_geom = WcsGeom.create(binsz=0.02, width=5.0, axes=[energy_axis])
    psfkernel = psfmap.get_psf_kernel(
        position=SkyCoord(1, 1, unit="deg"), geom=kern_geom, max_radius=1 * u.deg
    )
    assert_allclose(psfkernel.psf_kernel_map.data.sum(axis=(1, 2)), 1.0, atol=1e-7)

def test_psfmap_to_from_hdulist():
    psfmap = make_test_psfmap(0.15 * u.deg)
    hdulist = psfmap.to_hdulist()
    assert "PSF" in hdulist
    assert "PSF_BANDS" in hdulist
    assert "PSF_EXPOSURE" in hdulist
    assert "PSF_EXPOSURE_BANDS" in hdulist

    new_psfmap = PSFMap.from_hdulist(hdulist)
    assert_allclose(psfmap.psf_map.data, new_psfmap.psf_map.data)
    assert new_psfmap.psf_map.geom == psfmap.psf_map.geom
    assert new_psfmap.exposure_map.geom == psfmap.exposure_map.geom

def test_psfmap_read_write(tmp_path):
    psfmap = make_test_psfmap(0.15 * u.deg)

    psfmap.write(tmp_path / "tmp.fits")
    new_psfmap = PSFMap.read(tmp_path / "tmp.fits")

    assert_allclose(psfmap.psf_map.quantity, new_psfmap.psf_map.quantity)

def test_containment_radius_map():
    psf = fake_psf3d(0.15 * u.deg)
    pointing = SkyCoord(0, 0, unit="deg")
    energy_axis = MapAxis(nodes=[0.2, 1, 2], unit="TeV", name="energy_true")
    psf_theta_axis = MapAxis(nodes=np.linspace(0.0, 0.6, 30), unit="deg", name="rad")
    geom = WcsGeom.create(
        skydir=pointing, binsz=0.5, width=(4, 3), axes=[psf_theta_axis, energy_axis]
    )

    psfmap = make_psf_map(psf=psf, pointing=pointing, geom=geom)
    m = psfmap.containment_radius_map(energy_true=1 * u.TeV)
    coord = SkyCoord(0.3, 0, unit="deg")
    val = m.interp_by_coord(coord)
    assert_allclose(val, 0.226477, rtol=1e-2)

def test_psfmap_stacking():
    psfmap1 = make_test_psfmap(0.1 * u.deg, shape="flat")
    psfmap2 = make_test_psfmap(0.1 * u.deg, shape="flat")
    psfmap2.exposure_map.quantity *= 2

    psfmap_stack = psfmap1.copy()
    psfmap_stack.stack(psfmap2)
    assert_allclose(psfmap_stack.psf_map.data, psfmap1.psf_map.data)
    assert_allclose(psfmap_stack.exposure_map.data, psfmap1.exposure_map.data * 3)

    psfmap3 = make_test_psfmap(0.3 * u.deg, shape="flat")

    psfmap_stack = psfmap1.copy()
    psfmap_stack.stack(psfmap3)

    assert_allclose(psfmap_stack.psf_map.data[0, 40, 20, 20], 0.0)
    assert_allclose(psfmap_stack.psf_map.data[0, 20, 20, 20], 1.768388, rtol=1e-6)
    assert_allclose(psfmap_stack.psf_map.data[0, 0, 20, 20], 17.683883, rtol=1e-6)

# TODO: add a test comparing make_mean_psf and PSFMap.stack for a set of observations in an Observations

def test_sample_coord():
    psf_map = make_test_psfmap(0.1 * u.deg, shape="gauss")

    coords_in = MapCoord(
        {"lon": [0, 0] * u.deg, "lat": [0, 0.5] * u.deg, "energy_true": [1, 3] * u.TeV},
        frame="icrs",
    )

    coords = psf_map.sample_coord(map_coord=coords_in)
    assert coords.frame == "icrs"
    assert len(coords.lon) == 2
    assert_allclose(coords.lon, [0.074855, 0.042655], rtol=1e-3)
    assert_allclose(coords.lat, [-0.101561, 0.347365], rtol=1e-3)

def test_sample_coord_gauss():
    psf_map = make_test_psfmap(0.1 * u.deg, shape="gauss")

    lon, lat = np.zeros(10000) * u.deg, np.zeros(10000) * u.deg
    energy = np.ones(10000) * u.TeV
    coords_in = MapCoord.create(
        {"lon": lon, "lat": lat, "energy_true": energy}, frame="icrs"
    )
    coords = psf_map.sample_coord(coords_in)

    assert_allclose(np.mean(coords.skycoord.data.lon.wrap_at("180d").deg), 0, atol=2e-3)
    assert_allclose(np.mean(coords.lat), 0, atol=2e-3)

def make_psf_map_obs(geom, obs):
    exposure_map = make_map_exposure_true_energy(
        geom=geom.squash(axis_name="rad"),
        pointing=obs.pointing_radec,
        aeff=obs.aeff,
        livetime=obs.observation_live_time_duration,
    )

    psf_map = make_psf_map(
        geom=geom, psf=obs.psf, pointing=obs.pointing_radec, exposure_map=exposure_map
    )
    return psf_map

@requires_data()
@pytest.mark.parametrize(
    "pars",
    [
        {
            "energy": None,
            "rad": None,
            "energy_shape": 32,
            "psf_energy": 0.8659643,
            "rad_shape": 144,
            "psf_rad": 0.0015362848,
            "psf_exposure": 3.14711e12 * u.Unit("cm2 s"),
            "psf_value_shape": (32, 144),
            "psf_value": 4369.96391 * u.Unit("sr-1"),
        },
        {
            "energy": MapAxis.from_energy_bounds(1, 10, 100, "TeV", name="energy_true"),
            "rad": None,
            "energy_shape": 100,
            "psf_energy": 1.428893959,
            "rad_shape": 144,
            "psf_rad": 0.0015362848,
            "psf_exposure": 4.723409e12 * u.Unit("cm2 s"),
            "psf_value_shape": (100, 144),
            "psf_value": 3714.303683 * u.Unit("sr-1"),
        },
        {
            "energy": None,
            "rad": MapAxis.from_nodes(np.arange(0, 2, 0.002), unit="deg", name="rad"),
            "energy_shape": 32,
            "psf_energy": 0.8659643,
            "rad_shape": 1000,
            "psf_rad": 0.000524,
            "psf_exposure": 3.14711e12 * u.Unit("cm2 s"),
            "psf_value_shape": (32, 1000),
            "psf_value": 2.619861 * u.Unit("deg-2"),
        },
        {
            "energy": MapAxis.from_energy_bounds(1, 10, 100, "TeV", name="energy_true"),
            "rad": MapAxis.from_nodes(np.arange(0, 2, 0.002), unit="deg", name="rad"),
            "energy_shape": 100,
            "psf_energy": 1.428893959,
            "rad_shape": 1000,
            "psf_rad": 0.000524,
            "psf_exposure": 4.723409e12 * u.Unit("cm2 s"),
            "psf_value_shape": (100, 1000),
            "psf_value": 2.238232 * u.Unit("deg-2"),
        },
    ],
)
def test_make_psf(pars, data_store):
    obs = data_store.obs(23523)
    psf = obs.psf

    if pars["energy"] is None:
        energy_axis = psf.axes["energy_true"]
    else:
        energy_axis = pars["energy"]

    if pars["rad"] is None:
        rad_axis = psf.axes["rad"]
    else:
        rad_axis = pars["rad"]

    position = SkyCoord(83.63, 22.01, unit="deg")

    geom = WcsGeom.create(
        skydir=position, npix=(3, 3), axes=[rad_axis, energy_axis], binsz=0.2
    )

    psf_map = make_psf_map_obs(geom, obs)
    psf = psf_map.to_region_nd_map(position)

    axis = psf.psf_map.geom.axes["energy_true"]
    assert axis.unit == "TeV"
    assert axis.nbin == pars["energy_shape"]
    assert_allclose(axis.center.value[15], pars["psf_energy"], rtol=1e-3)

    rad_axis = psf.psf_map.geom.axes["rad"]
    assert rad_axis.unit == "deg"
    assert rad_axis.nbin == pars["rad_shape"]
    assert_allclose(rad_axis.center.to_value("rad")[15], pars["psf_rad"], rtol=1e-3)

    exposure = psf.exposure_map.quantity.squeeze()
    assert exposure.unit == "m2 s"
    assert exposure.shape == (pars["energy_shape"],)
    assert_allclose(exposure[15], pars["psf_exposure"], rtol=1e-3)

    data = psf.psf_map.quantity.squeeze()
    assert data.unit == "deg-2"
    assert data.shape == pars["psf_value_shape"]
    assert_allclose(data[15, 50], pars["psf_value"], rtol=1e-3)

@requires_data()
def test_make_mean_psf(data_store):
    observations = data_store.get_observations([23523, 23526])
    position = SkyCoord(83.63, 22.01, unit="deg")

    psf = observations[0].psf

    geom = WcsGeom.create(
        skydir=position,
        npix=(3, 3),
        axes=psf.axes[["rad", "energy_true"]],
        binsz=0.2,
    )

    psf_map_1 = make_psf_map_obs(geom, observations[0])
    psf_map_2 = make_psf_map_obs(geom, observations[1])

    stacked_psf = psf_map_1.copy()
    stacked_psf.stack(psf_map_2)

    psf = stacked_psf.to_region_nd_map(position).psf_map

    assert not np.isnan(psf.quantity.squeeze()).any()
    assert_allclose(psf.quantity.squeeze()[22, 22], 12206.1665 / u.sr, rtol=1e-3)

@requires_data()
@pytest.mark.parametrize("position", ["0d 0d", "180d 0d", "0d 90d", "180d -90d"])
def test_psf_map_read(position):
    position = SkyCoord(position)
    filename = "$GAMMAPY_DATA/fermi_3fhl/fermi_3fhl_psf_gc.fits.gz"
    psf = PSFMap.read(filename, format="gtpsf")

    value = psf.containment(
        position=position, energy_true=100 * u.GeV, rad=0.1 * u.deg
    )

    assert_allclose(value, 0.682022, rtol=1e-5)
    assert psf.psf_map.unit == "sr-1"

def test_psf_map_write_gtpsf(tmpdir):
    energy_axis_true = MapAxis.from_energy_bounds(
        "1 TeV", "10 TeV", nbin=3, name="energy_true"
    )
    geom = RegionGeom.create("icrs;circle(0, 0, 0.1)")
    psf = PSFMap.from_gauss(
        energy_axis_true=energy_axis_true, sigma=[0.1, 0.2, 0.3] * u.deg, geom=geom
    )
    psf.exposure_map = Map.from_geom(geom.to_cube([energy_axis_true]), unit="cm2 s")

    filename = tmpdir / "test_psf.fits"
    psf.write(filename, format="gtpsf")

    psf = PSFMap.read(filename, format="gtpsf")

    value = psf.containment_radius(
        energy_true=energy_axis_true.center, fraction=0.394
    )

    assert_allclose(value, [0.1, 0.2, 0.3] * u.deg, rtol=1e-5)
    assert psf.psf_map.unit == "sr-1"

def test_to_image():
    psfmap = make_test_psfmap(0.15 * u.deg)

    psf2D = psfmap.to_image()
    assert_allclose(psf2D.psf_map.geom.data_shape, (1, 100, 25, 25))
    assert_allclose(psf2D.exposure_map.geom.data_shape, (1, 1, 25, 25))
    assert_allclose(psf2D.psf_map.data[0][0][12][12], 7.068315, rtol=1e-2)

def test_psf_map_from_gauss():
    energy_axis = MapAxis.from_nodes(
        [1, 3, 10], name="energy_true", interp="log", unit="TeV"
    )
    rad = np.linspace(0, 1.5, 100) * u.deg
    rad_axis = MapAxis.from_nodes(rad, name="rad", unit="deg")

    # define sigmas starting at 0.1 in steps of 0.1 deg
    sigma = [0.1, 0.2, 0.4] * u.deg

    # with energy-dependent sigma
    psfmap = PSFMap.from_gauss(energy_axis, rad_axis, sigma)

    assert psfmap.psf_map.geom.axes[0] == rad_axis
    assert psfmap.psf_map.geom.axes[1] == energy_axis
    assert psfmap.psf_map.unit == "sr-1"
    assert psfmap.psf_map.data.shape == (3, 100, 1, 2)

    radius = psfmap.containment_radius(
        fraction=0.394, energy_true=[1, 3, 10] * u.TeV
    )
    assert_allclose(radius, sigma, rtol=0.01)

    # test that it won't work with different number of sigmas and energies
    with pytest.raises(ValueError):
        PSFMap.from_gauss(energy_axis, rad_axis, sigma=[1, 2] * u.deg)

def test_psf_map_from_gauss_const_sigma():
    energy_axis = MapAxis.from_nodes(
        [1, 3, 10], name="energy_true", interp="log", unit="TeV"
    )
    rad = np.linspace(0, 1.5, 100) * u.deg
    rad_axis = MapAxis.from_nodes(rad, name="rad", unit="deg")

    # with constant sigma
    psfmap = PSFMap.from_gauss(energy_axis, rad_axis, sigma=0.1 * u.deg)
    assert psfmap.psf_map.geom.axes[0] == rad_axis
    assert psfmap.psf_map.geom.axes[1] == energy_axis
    assert psfmap.psf_map.unit == Unit("sr-1")
    assert psfmap.psf_map.data.shape == (3, 100, 1, 2)

    radius = psfmap.containment_radius(
        energy_true=[1, 3, 10] * u.TeV, fraction=0.394
    )
    assert_allclose(radius, 0.1 * u.deg, rtol=0.01)

@requires_data()
@requires_dependency("matplotlib")
def test_psf_map_plot_containment_radius():
    filename = "$GAMMAPY_DATA/fermi_3fhl/fermi_3fhl_psf_gc.fits.gz"
    psf = PSFMap.read(filename, format="gtpsf")

    with mpl_plot_check():
        psf.plot_containment_radius_vs_energy()

@requires_data()
@requires_dependency("matplotlib")
def test_psf_map_plot_psf_vs_rad():
    filename = "$GAMMAPY_DATA/fermi_3fhl/fermi_3fhl_psf_gc.fits.gz"
    psf = PSFMap.read(filename, format="gtpsf")

    with mpl_plot_check():
        psf.plot_psf_vs_rad()

max_stars_repo_name: mcerruti/gammapy

max_issues_repo_name: mcerruti/gammapy

max_forks_repo_name: mcerruti/gammapy

reason: import numpy,from numpy,import astropy,from astropy


text:
from flatland.evaluators.client import FlatlandRemoteClient
from flatland.core.env_observation_builder import DummyObservationBuilder
# from my_observation_builder import CustomObservationBuilder
from src.graph_observations import GraphObsForRailEnv
from src.predictions import ShortestPathPredictorForRailEnv
import numpy as np
import time
from src.util.tree_builder import Agent_Tree
from run_local_3 import naive_solver, reroute_solver
from flatland.utils.rendertools import RenderTool

max_prediction_depth = 80
#####################################################################
# Instantiate a Remote Client
#####################################################################
remote_client = FlatlandRemoteClient()

#####################################################################
# Define your custom controller
#
# which can take an observation, and the number of agents and
# compute the necessary action for this step for all (or even some)
# of the agents
#####################################################################

def my_controller(env, obs, number_of_agents, agent_idx, cell_sequence):
    _action = reroute_solver(cell_sequence, naive_solver(env, obs), env, agent_idx)
    # _action = naive_solver(env, obs)
    for k in _action.keys():

        if env.agents[k].status == 1:
            # env.dev_pred_dict[k] = cell_sequence[k][agent_idx[k]:]
            agent_idx[k] += env.agents[k].speed_data['speed']
    # _action = {}
    # for _idx in range(number_of_agents):
    #     _action[_idx] = np.random.randint(0, 5)
    return _action

#####################################################################
# Instantiate your custom Observation Builder
#
# You can build your own Observation Builder by following
# the example here :
# https://gitlab.aicrowd.com/flatland/flatland/blob/master/flatland/envs/observations.py#L14
#####################################################################
my_observation_builder = observation_builder = GraphObsForRailEnv(predictor=ShortestPathPredictorForRailEnv(max_depth=max_prediction_depth),
                                             bfs_depth=200)

# Or if you want to use your own approach to build the observation from the env_step,
# please feel free to pass a DummyObservationBuilder() object as mentioned below,
# and that will just return a placeholder True for all observation, and you
# can build your own Observation for all the agents as your please.
# my_observation_builder = DummyObservationBuilder()

#####################################################################
# Main evaluation loop
#
# This iterates over an arbitrary number of env evaluations
#####################################################################
evaluation_number = 0
while True:

    evaluation_number += 1
    # Switch to a new evaluation environemnt
    #
    # a remote_client.env_create is similar to instantiating a
    # RailEnv and then doing a env.reset()
    # hence it returns the first observation from the
    # env.reset()
    #
    # You can also pass your custom observation_builder object
    # to allow you to have as much control as you wish
    # over the observation of your choice.
    time_start = time.time()
    observation, info = remote_client.env_create(
        obs_builder_object=my_observation_builder
    )
    env_creation_time = time.time() - time_start
    if not observation:
        #
        # If the remote_client returns False on a `env_create` call,
        # then it basically means that your agent has already been
        # evaluated on all the required evaluation environments,
        # and hence its safe to break out of the main evaluation loop
        break

    print("Evaluation Number : {}".format(evaluation_number))

    #####################################################################
    # Access to a local copy of the environment
    #
    #####################################################################
    # Note: You can access a local copy of the environment
    # by using :
    #       remote_client.env
    #
    # But please ensure to not make any changes (or perform any action) on
    # the local copy of the env, as then it will diverge from
    # the state of the remote copy of the env, and the observations and
    # rewards, etc will behave unexpectedly
    #
    # You can however probe the local_env instance to get any information
    # you need from the environment. It is a valid RailEnv instance.
    local_env = remote_client.env
    number_of_agents = len(local_env.agents)

    # Now we enter into another infinite loop where we
    # compute the actions for all the individual steps in this episode
    # until the episode is `done`
    #
    # An episode is considered done when either all the agents have
    # reached their target destination
    # or when the number of time steps has exceed max_time_steps, which
    # is defined by :
    #
    # max_time_steps = int(4 * 2 * (env.width + env.height + 20))
    #
    time_taken_by_controller = []
    time_taken_per_step = []
    steps = 0

    env_renderer = RenderTool(local_env)

    img = env_renderer.render_env(show=True,
                                  show_inactive_agents=False,
                                  show_predictions=True,
                                  show_observations=True,
                                  frames=True,
                                  return_image=True)

    tree_dict = {}
    agent_idx = {}
    for idx, agent in enumerate(local_env.agents):
        tree = Agent_Tree(idx, agent.initial_position)
        tree.build_tree(observation, local_env)
        tree_dict[idx] = tree
        agent_idx[idx] = 0

    success, rout_list = observation_builder.rerouting(observation, tree_dict)  # returns list of re-routes
    cell_sequence = observation_builder.cells_sequence.copy()
    # for id, route in rout_list:  #viualising the re-routes
    #     env.dev_pred_dict[id] = route
    if success:
        print("Success!")
    else:
        print("Failed!")

    while True:
        #####################################################################
        # Evaluation of a single episode
        #
        #####################################################################
        # Compute the action for this step by using the previously
        # defined controller
        time_start = time.time()
        action = my_controller(local_env, observation, number_of_agents, agent_idx, cell_sequence)
        time_taken = time.time() - time_start
        time_taken_by_controller.append(time_taken)

        img = env_renderer.render_env(show=True,
                                      show_inactive_agents=False,
                                      show_predictions=True,
                                      show_observations=True,
                                      frames=True,
                                      return_image=True)

        # Perform the chosen action on the environment.
        # The action gets applied to both the local and the remote copy
        # of the environment instance, and the observation is what is
        # returned by the local copy of the env, and the rewards, and done and info
        # are returned by the remote copy of the env
        time_start = time.time()
        observation, all_rewards, done, info = remote_client.env_step(action)
        steps += 1
        time_taken = time.time() - time_start
        time_taken_per_step.append(time_taken)

        if done['__all__']:
            print("Reward : ", sum(list(all_rewards.values())))
            #
            # When done['__all__'] == True, then the evaluation of this
            # particular Env instantiation is complete, and we can break out
            # of this loop, and move onto the next Env evaluation
            break

    np_time_taken_by_controller = np.array(time_taken_by_controller)
    np_time_taken_per_step = np.array(time_taken_per_step)
    print("=" * 100)
    print("=" * 100)
    print("Evaluation Number : ", evaluation_number)
    print("Current Env Path : ", remote_client.current_env_path)
    print("Env Creation Time : ", env_creation_time)
    print("Number of Steps : ", steps)
    print("Mean/Std of Time taken by Controller : ", np_time_taken_by_controller.mean(),
          np_time_taken_by_controller.std())
    print("Mean/Std of Time per Step : ", np_time_taken_per_step.mean(), np_time_taken_per_step.std())
    print("=" * 100)

print("Evaluation of all environments complete...")
########################################################################
# Submit your Results
#
# Please do not forget to include this call, as this triggers the
# final computation of the score statistics, video generation, etc
# and is necesaary to have your submission marked as successfully evaluated
########################################################################
print(remote_client.submit())

max_stars_repo_name: RaiAnant/Flatland-rl

max_issues_repo_name: RaiAnant/Flatland-rl

max_forks_repo_name: RaiAnant/Flatland-rl


text:
from nose.tools import *
import analysis.core as wdt
import numpy as np

class TestClass:

    @classmethod
    def setup_class(cls):
        filename = './tests/wdt_runs/S0100/W0100/runs/run1_res.m'
        cls.data = wdt.DataFile(filename)

    def test_DataFile_value(self):
        """ Verify DataFile is getting values correctly """
        ana_keff = np.array([[ 1.77253  ,  1.76093  ,  0.0114943]])
        assert_true(np.allclose(ana_keff, self.data.get_data('ANA_KEFF')))

    def test_DataFile_error(self):
        """ Verify DataFile is getting errors correctly """
        ana_keff = np.array([[ 0.00024,  0.0002 ,  0.00431]])
        assert_true(np.allclose(ana_keff, self.data.get_data('ANA_KEFF', err=True)))

    def test_DataFile_reshape_size(self):
        """ Verify matrix data, when reshaped, forms a square matrix of
        the correct size """
        inf_s0 = self.data.get_data('INF_S0',reshape=True)
        assert_true(np.shape(inf_s0) == (11,11))

    def test_DataFile_reshape_vals(self):
        """ Verify values after reshape """
        inf_s0_0 = np.array([ 9.71843000e-02, 5.00055000e-02,
                              8.21536000e-05, 2.94249000e-08, 0.00000000e+00,
                              0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
                              0.00000000e+00, 0.00000000e+00, 0.00000000e+00])
        assert_true(np.allclose(inf_s0_0,self.data.get_data('INF_S0',reshape=True)[0]))

    def test_DataFile_reshape_err(self):
        """ Verify errors after reshape """
        inf_s0_0 = np.array([ 2.70000000e-04, 5.60000000e-04,
                              1.54100000e-02, 1.00000000e+00, 0.00000000e+00,
                              0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
                              0.00000000e+00, 0.00000000e+00, 0.00000000e+00])
        assert_true(np.allclose(inf_s0_0,self.data.get_data('INF_S0',reshape=True, err=True)[0]))

    def test_DataFile_CPU(self):
        """ DataFile gets cpu value correctly """
        assert_true(np.allclose(49.391399999999997, self.data.get_cpu()))

    @raises(KeyError)
    def test_DataFile_bad_label(self):
        """ Sending a bad Serpent parameter should return a key error """
        ana_keff = self.data.get_data('WRONG_LABEL')

    def test_DataFile_FOM_of_constant(self):
        """ Requesting the FOM of a constant should return 0 """
        eq_(0, self.data.get_fom('TOT_CPU_TIME'))

    def test_DataFile_FOM_values(self):
        """ DataFile should return correct value when calculating FOM """
        fom = np.array([[ 15578.97788826, 133112.68680858,
                          140210.80099436, 133112.68680858, 147892.18161859,
                          114775.73505434, 140210.80099436, 258245.40387226,
                          1198014.18127723, 277728.93914383, 56240.11017663]])
        ok_(np.allclose(fom, self.data.get_fom('INF_FLX')))

    def test_DataFile_FOM_values_cycles(self):
        """ DataFile should return correct value when calculating FOM using cycles """
        fom = np.array([[ 7694.6752847 , 65746.21959237,
                          69252.07756233, 65746.21959237, 73046.01899196,
                          56689.34240363, 69252.07756233, 127551.02040816,
                          591715.97633136, 137174.21124829, 27777.77777778]])
        ok_(np.allclose(fom, self.data.get_fom('INF_FLX', cpu=False)))

    def test_DataFile_FOM_reshape_shape(self):
        """ DataFile should return the correct shape when calculating
        FOM using reshape """

        eq_((11,11), np.shape(self.data.get_fom('INF_SP0',reshape=True)))
    def test_DataFile_FOM_reshape_values(self):
        """ DataFile should return the correct values when calculating FOM using reshape """
        fom = np.array([ 2.77728939e+05, 6.23159116e+04,
                         8.52596187e+01, 2.02464397e-02, 0.00000000e+00,
                         0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
                         0.00000000e+00, 0.00000000e+00, 0.00000000e+00])
        ok_(np.allclose(fom, self.data.get_fom('INF_SP0', reshape=True)[0]))


text:
import apache_beam as beam
import logging
from joblib import load
import numpy as np
import pandas as pd

from google.cloud import storage
from apache_beam.options.pipeline_options import StandardOptions, GoogleCloudOptions, SetupOptions, PipelineOptions

from sklearn.ensemble import RandomForestClassifier

#setup global
dummy_dict = {'Yes': 1, 'No': 0}
internet_dict = {'No': 0, 'No internet service': 1, 'Yes': 2}
yesno_cols = ['Partner', 'Dependents', 'PhoneService', 'PaperlessBilling']
internet_cols = ['OnlineSecurity', 'OnlineBackup', 'DeviceProtection', 'TechSupport', 'StreamingTV', 'StreamingMovies']
# =============================================================================
# Build and run the pipeline
# =============================================================================
def run(argv=None):
    pipeline_options = PipelineOptions(flags=argv)

    google_cloud_options = pipeline_options.view_as(GoogleCloudOptions)
    google_cloud_options.project = 'your-project' #change this
    google_cloud_options.job_name = 'telco-churn-prediction'
    google_cloud_options.staging_location = 'gs://path/to/your/staging' #change this
    google_cloud_options.temp_location = 'gs://path/to/your/temp' #change this
    pipeline_options.view_as(StandardOptions).runner = 'DataflowRunner'
    pipeline_options.view_as(SetupOptions).save_main_session = True
    pipeline_options.view_as(SetupOptions).setup_file = "./setup.py"
    logging.info("Pipeline arguments: {}".format(pipeline_options))

    # table_schema = 'customerID: STRING, prediction: FLOAT'
    query = ('select * from `your-project.Telco_Churn.input`')
    bq_source = beam.io.BigQuerySource(query=query, use_standard_sql=True)
    p = beam.Pipeline(options=pipeline_options)
    (p
     | "Read data from BQ" >> beam.io.Read(bq_source)
     | "Preprocess data" >> beam.ParDo(FormatInput())
     | "predicting" >> beam.ParDo(
                PredictSklearn(project='your-project', bucket_name='your-bucket-name', model_path='/path/to/model_rf.joblib',
                               destination_name='model_rf.joblib'))
     | "Write data to BQ" >> beam.io.WriteToBigQuery(table='prediction', dataset='Telco_Churn', project='your-project',
                                                     # schema=table_schema,
                                                     # create_disposition=beam.io.BigQueryDisposition.CREATE_IF_NEEDED,
                                                     write_disposition=beam.io.BigQueryDisposition.WRITE_TRUNCATE)
     )

    result = p.run()
    result.wait_until_finish()

# =============================================================================
# Function to download model from bucket
# =============================================================================
def download_blob(bucket_name=None, source_blob_name=None, project=None, destination_file_name=None):
    storage_client = storage.Client(project)
    bucket = storage_client.get_bucket(bucket_name)
    blob = bucket.blob(source_blob_name)
    blob.download_to_filename(destination_file_name)

class FormatInput(beam.DoFn):
    def process(self, element):
        """ Format the input to the desired shape"""
        df = pd.DataFrame([element], columns=element.keys())
        df[yesno_cols] = df[yesno_cols].apply(lambda x: x.map(dummy_dict))
        df[internet_cols] = df[internet_cols].apply(lambda x: x.map(internet_dict))
        df['gender'] = df['gender'].map({'Female': 0, 'Male': 1})
        df['MultipleLines'] = df['MultipleLines'].map({'No': 0, 'No phone service': 1, 'Yes': 2})
        df['InternetService'] = df['InternetService'].map({'DSL': 0, 'Fiber optic': 1, 'No': 2})
        df['Contract'] = df['Contract'].map({'Month-to-month': 0, 'One year': 1, 'Two year': 2})
        df['PaymentMethod'] = df['PaymentMethod'].map(
            {'Bank transfer (automatic)': 0, 'Credit card (automatic)': 1, 'Electronic check': 2, 'Mailed check': 3})
        output = df.to_dict('records')
        return output #return a dict for easier comprehension

class PredictSklearn(beam.DoFn):

    def __init__(self, project=None, bucket_name=None, model_path=None, destination_name=None):
        self._model = None
        self._project = project
        self._bucket_name = bucket_name
        self._model_path = model_path
        self._destination_name = destination_name

    def setup(self):
        """Download sklearn model from GCS"""
        logging.info(
            "Sklearn model initialisation {}".format(self._model_path))
        download_blob(bucket_name=self._bucket_name, source_blob_name=self._model_path,
                      project=self._project, destination_file_name=self._destination_name)
        # unpickle sklearn model
        self._model = load(self._destination_name)

    def process(self, element):
        """Predicting using developed model"""
        input_dat = {k: element[k] for k in element.keys() if k not in ['customerID']}
        tmp = np.array(list(i for i in input_dat.values()))
        tmp = tmp.reshape(1, -1)
        element["prediction"] = self._model.predict_proba(tmp)[:,1].item()
        output = {k: element[k] for k in element.keys() if k in ['customerID', 'prediction']}
        output['customerID'] = str(output['customerID'])
        return [output]

# log the output
if __name__ == '__main__':
    logging.getLogger().setLevel(logging.INFO)
    run()

max_stars_repo_name: nghiamab/dataflow-demo

max_issues_repo_name: damioke/dataflow-demo

max_forks_repo_name: damioke/dataflow-demo


text:
from argparse import ArgumentParser
from data_acquirers import RealDataAcquirer, SimulatedDataAcquirer, RecordedDataAcquirer
from polarimeter import Polarimeter
import time
import os
import numpy as np
import json

class Runner():
    def __init__(self, config):
        self.config = config

    def experiment(self):
        capture_time = 1
        acquirer = RealDataAcquirer(capture_time, self.config)
        p = Polarimeter()
        phase_diff, stdev = p.measure(*acquirer.acquire())

        print('{} +/- {} at {}'.format(phase_diff, stdev, acquirer.lastDataTimestamp()))

    def exper_time(self):
        capture_time = 1
        acquirer = RealDataAcquirer(capture_time, self.config)
        p = Polarimeter()
        x = 5 # how many flushes to file
        y = 10 # experiment count per one flush
        timestamp = int(time.time())
        results_file = os.path.join(self.config["experiment_results_folder"] + str(timestamp) + ".txt")

        with open(results_file, "w") as f:
            for _ in range(x):
                for _ in range(y):
                    data = acquirer.acquire()
                    phase_diff, stdev = p.measure(*data)
                    f.write("{},{},{}\n".format(acquirer.lastDataTimestamp(), phase_diff, stdev))
                #time.sleep(300)
                f.flush()

    def exper_time2(self):
        phase_difference = 60
        acquirer = SimulatedDataAcquirer(phase_difference=phase_difference)
        p = Polarimeter()

        x = 6 # how many flushes to file, no of data points
        y = 10 # experiments taken into mean
        timestamp = int(time.time())
        results_file = os.path.join(self.config["experiment_results_folder"] + str(timestamp) + ".txt")

        with open(results_file, "w") as f:
            for _ in range(x):
                t = int(time.time())
                data = []
                for _ in range(y):
                    phase_diff, stdev = p.measure(*acquirer.acquire())
                    data.append([acquirer.lastDataTimestamp(), phase_diff, stdev])
                phase = []
                stdev = []
                for i in range(len(data)):
                    phase.append(data[i][1])
                    stdev.append(data[i][2])
                f.write("{},{},{}\n".format(t, np.mean(phase), np.mean(stdev)))
                f.flush()

    def rerunning(self):
        results_file = self.config["rerunning_results_file"]
        with open(results_file, 'w') as f:
            data_folder_path = self.config["rerunning_data_folder"]
            acquirer = RecordedDataAcquirer(data_folder_path)
            p = Polarimeter()
            while acquirer.isDataAvailable():
                data = acquirer.acquire()
                phase_diff, stdev = p.measure(*data)
                f.write("{},{},{}\n".format(acquirer.lastDataTimestamp(), phase_diff, stdev))

    def laser_stability(self):
        x = 7000 #how many flushes to file
        y = 10 #no of experiments in a single flush
        timestamp = int(time.time())
        results_file = os.path.join(self.config["laser_folder"] + str(timestamp) + "_stab.txt")

        with open(results_file, "w") as f:
            for _ in range(x):
                for _ in range(y):
                    times = int(time.time())
                    acquirer = RealDataAcquirer(1, self.config)
                    one_acq = acquirer.acquire()
                    chA = one_acq[:,1]
                    chB = one_acq[:,2]
                    final_chA = np.mean(chA)
                    final_chB = np.mean(chB)
                    f.write("{},{},{}\n".format(times, final_chA, final_chB))
            f.flush()

def main():
    with open('config.json') as json_data_file:
        config = json.load(json_data_file)

    parser = ArgumentParser()
    parser.add_argument("scenario", nargs=1)
    args = parser.parse_args()
    scenario = args.scenario[0]
    if scenario == 'experiment':
        Runner(config).experiment()
    elif scenario == 'exper_time':
        Runner(config).exper_time()
    elif scenario == 'exper_time2':
        Runner(config).exper_time2()
    elif scenario == 'laser_stability':
        Runner(config).laser_stability()
    elif scenario == 'rerunning':
        Runner(config).rerunning()
    else:
        raise Exception("Uknown scenario: {}".format(args.scenario))

if __name__ == "__main__":
    main()

max_stars_repo_name: malgorzatakim/polarimeter

max_issues_repo_name: malgorzatakim/polarimeter

max_forks_repo_name: malgorzatakim/polarimeter


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""Create template graphs.

Usage:
  graficop2.py [options]

Options:
  -h --help           Show this screen.
  -u                  Print also unics.

"""

import csv, re, pdb, ast, time, os, sys

from docopt import docopt

import datetime

import matplotlib.pyplot as plot

import numpy as np

import matplotlib.dates as mdates

init = datetime.date(2001, 1, 1)

def allgraph(unic):

    csv.field_size_limit(sys.maxsize)

    alldict = csv.reader(open("alldict.csv", "r"))

    alltemplates = {}

    print("Creating output dictionary...")

    for line in alldict:

        page = line[0]

        pagedict = ast.literal_eval(line[1])

        sort = sorted(pagedict)

        oldcount = {}

        for s in sort:

            new = 0

            for tmp in pagedict[s]:

                if unic:

                    new = 1

                else:

                    d = pagedict[s]

                    new = d[tmp]

                try:

                    old = oldcount[tmp]

                except KeyError:

                    old = 0

                add = new - old

                if (add != 0):

                    try:

                        singletemplate = alltemplates[tmp]

                    except KeyError:

                        singletemplate = {}

                    try:

                        singletemplate[s] = singletemplate[s] + add

                    except KeyError:

                        singletemplate[s] = add

                    alltemplates[tmp] = singletemplate

                oldcount[tmp] = new

            for o in oldcount:

                if o not in pagedict[s] and oldcount[o] != 0:

                    singletemplate = alltemplates[o]

                    try:

                        singletemplate[s] = singletemplate[s] - oldcount[o]

                    except KeyError:

                        singletemplate[s] = - oldcount[o]

                    oldcount[o] = 0

                    alltemplates[o] = singletemplate

    print("Creating .csv files...")

    if not os.path.exists("grafici"):

        os.makedirs("grafici")

    if unic:

        if not os.path.exists("grafici/unic"):

            os.makedirs("grafici/unic")

    if not os.path.exists("grafici/total"):

        os.makedirs("grafici/total")

    alldict = {}

    for singletemplate in alltemplates:

        name = singletemplate

        if len(name) > 130:

            continue

        singletemplate = alltemplates[singletemplate]

        if unic:

            outu = csv.writer(open("grafici/unic/" + name + ".csv", "w"))

        out = csv.writer(open("grafici/total/" + name + ".csv", "w"))

        sort = sorted(singletemplate)

        value = 0

        valueu = 0

        x = []

        y = []

        yu = []

        for s in sort:

            value = value + singletemplate[s]

            valueu = valueu + 1

            x.append(s)

            y.append(value)

            y.append(valueu)

        l = len(x)

        n = 0

        while n < l:

            if unic:

                outu.writerow([x[n], yu[n]])

            out.writerow([x[n], y[n]])

            n = n + 1

if __name__ == "__main__":

    arguments = docopt(__doc__)

    allgraph(arguments['-u'])

max_stars_repo_name: Tihacker/Wikipedia-Templates-Analysis

max_issues_repo_name: Tihacker/Wikipedia-Templates-Analysis

max_forks_repo_name: Tihacker/Wikipedia-Templates-Analysis


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Feb  5 14:47:50 2020

@author: valery
"""
from pathlib import Path

def set_paths(root, case, run, sens='grad'):
    '''
    Prepare folders for Spyking Circus

    Parameters
    ----------
    root : str
        Full path to the case
    case : str
        Patient ID
    run : str
        Recording designation which will be analyzed
    sens : str, optional
        'grad' or 'mag'. The default is 'grad'.

    '''
    from pkg_resources import resource_filename
    ### Main
    paths = {}
    paths['code_source'] = Path.cwd()
    paths['case_meg']    = Path(root)
    paths['case']        = case
    paths['SPC_root']    = paths['case_meg'] / 'CIRCUS'
    if not paths['SPC_root'].is_dir():
        paths['SPC_root'].mkdir(exist_ok=True)

    ### Package paths
    paths['circus_pkg'] = Path(resource_filename('circus','clustering.py')).parent
    paths['circus_updates'] =  paths['code_source'] / 'spyking_circus_updates'

    ### Files
    paths['art_cor']         = paths['case_meg']/'MEG_data'/'tsss_mc_artefact_correction'
    if sorted(paths['art_cor'].glob('*{}*'.format(run))) != []:
        paths['fif_file']    = sorted(paths['art_cor'].glob('*{}*'.format(run)))[0]
        paths['SPC']         = paths['SPC_root'] / paths['fif_file'].stem #folder
        paths['npy_file']    = paths['SPC']/'{}.npy'.format(paths['fif_file'].stem)
        paths['SPC_output']  = paths['SPC']/'{}_CIRCUS_{}'.format(paths['case'], sens)
        paths['SPC_params']  = paths['npy_file'].with_suffix('.params')
        paths['SPC_results'] = paths['SPC_root'] / 'Results'
        if not paths['SPC'].is_dir():
             paths['SPC'].mkdir(exist_ok=True)
        if not paths['SPC_output'].is_dir():
            paths['SPC_output'].mkdir(exist_ok=True)
        if not paths['SPC_results'].is_dir():
            paths['SPC_results'].mkdir(exist_ok=True)
    else: print("Error. No fif file")

    ### Folders for plots
    paths['Templates']     = paths['SPC'] / 'Templates_{}'.format(sens)
    paths['temporal_svg']  = paths['SPC'] / 'temporal_svg'
    if not paths['Templates'].is_dir():
         paths['Templates'].mkdir(exist_ok=True)
    if not paths['temporal_svg'].is_dir():
         paths['temporal_svg'].mkdir(exist_ok=True)

    return paths

def fif_to_npy(fif, npy):
    '''
    Convert and save fif-file to numpy file

    Parameters
    ----------
    fif : pathlib.PosixPath
        The path to fif-file
    npy : pathlib.PosixPath
        The path to numpy file

    '''
    import mne
    import numpy as np
    data = mne.io.read_raw_fif(str(fif), preload=True, verbose=False)
    data.pick_types(meg=True, eeg=False, stim=False, eog=False, ecg=False)
    np.save(npy, (data.get_data(reject_by_annotation='omit')).astype(float))
    del data