text:
import numba

from NN.Errors import *
from NN.Basic.Optimizers import *

# TODO: Support 'SAME' padding

@numba.jit([
    "void(int64, int64, int64, int64, float32[:,:,:,:],"
    "int64, int64, int64, float32[:,:,:,:], float32[:,:,:,:])"
], nopython=True)
def conv_bp(n, n_filters, out_h, out_w, dx_padded,
            filter_height, filter_width, sd, inner_weight, delta):
    for i in range(n):
        for f in range(n_filters):
            for j in range(out_h):
                for k in range(out_w):
                    for h in range(dx_padded.shape[1]):
                        jsd, ksd = j * sd, k * sd
                        for p in range(filter_height):
                            for q in range(filter_width):
                                dx_padded[i, h, jsd+p, ksd+q] += (
                                    inner_weight[f][h][p][q] * delta[i, f, j, k]
                                )

@numba.jit([
    "void(int64, int64, int64, int64, float32[:,:,:,:], float32[:,:,:,:],"
    "int64, int64, int64, int32[:,:,:,:,:])"
], nopython=True)
def max_pool(n, n_channels, out_h, out_w, x, out,
             pool_height, pool_width, sd, pos_cache):
    for i in range(n):
        for j in range(n_channels):
            for k in range(out_h):
                for l in range(out_w):
                    ksd, lsd = k * sd, l * sd
                    _max = x[i, j, ksd, lsd]
                    pos = (0, 0)
                    for p in range(pool_height):
                        for q in range(pool_width):
                            if x[i, j, ksd+p, lsd+q] > _max:
                                _max = x[i, j, ksd+p, lsd+q]
                                pos = (p, q)
                    pos_cache[i, j, k, l] = pos
                    out[i, j, k, l] = _max

@numba.jit([
    "void(int64, int64, int64, int64, int64, float32[:,:,:,:], float32[:,:,:,:], int32[:,:,:,:,:])"
], nopython=True)
def max_pool_bp(n, n_channels, out_h, out_w, sd, dx, delta, pos_cache):
    for i in range(n):
        for j in range(n_channels):
            for k in range(out_h):
                for l in range(out_w):
                    ksd, lsd = k * sd, l * sd
                    pos = pos_cache[i, j, k, l]
                    dx[i, j, ksd+pos[0], lsd+pos[1]] = delta[i, j, k, l]

# Abstract Layers

class Layer:
    LayerTiming = Timing()

    def __init__(self, shape):
        """
        :param shape: shape[0] = units of previous layer
                      shape[1] = units of current layer (self)
        """
        self._shape = shape
        self.parent = None
        self.child = None
        self.is_fc = False
        self.is_fc_base = False
        self.is_sub_layer = False
        self._last_sub_layer = None

    def __str__(self):
        return self.__class__.__name__

    def __repr__(self):
        return str(self)

    @property
    def name(self):
        return str(self)

    @property
    def shape(self):
        return self._shape

    @shape.setter
    def shape(self, value):
        self._shape = value

    @property
    def params(self):
        return self._shape,

    @property
    def special_params(self):
        return

    def set_special_params(self, dic):
        for key, value in dic.items():
            setattr(self, key, value)

    @property
    def root(self):
        return self

    @root.setter
    def root(self, value):
        raise BuildLayerError("Setting Layer's root is not permitted")

    @property
    def last_sub_layer(self):
        child = self.child
        if not child:
            return None
        while child.child:
            child = child.child
        return child

    @last_sub_layer.setter
    def last_sub_layer(self, value):
            self._last_sub_layer = value

    # Core

    def derivative(self, y, delta=None):
        return self._derivative(y, delta)

    @LayerTiming.timeit(level=1, prefix="[Core] ")
    def activate(self, x, w, bias=None, predict=False):
        if self.is_fc:
            x = x.reshape(x.shape[0], -1)
        if self.is_sub_layer:
            if bias is None:
                return self._activate(x, predict)
            return self._activate(x + bias, predict)
        if bias is None:
            return self._activate(x.dot(w), predict)
        return self._activate(x.dot(w) + bias, predict)

    @LayerTiming.timeit(level=1, prefix="[Core] ")
    def bp(self, y, w, prev_delta):
        if self.child is not None and isinstance(self.child, SubLayer):
            if not isinstance(self, SubLayer):
                return prev_delta
            return self._derivative(y, prev_delta)
        if isinstance(self, SubLayer):
            return self._derivative(y, prev_delta.dot(w.T) * self._root.derivative(y))
        return prev_delta.dot(w.T) * self._derivative(y)

    def _activate(self, x, predict):
        pass

    def _derivative(self, y, delta=None):
        pass

    # Util

    @staticmethod
    @LayerTiming.timeit(level=2, prefix="[Core Util] ")
    def safe_exp(y):
        return np.exp(y - np.max(y, axis=1, keepdims=True))

class SubLayer(Layer):
    def __init__(self, parent, shape):
        super(SubLayer, self).__init__(shape)
        self.parent = parent
        parent.child = self
        self._root = None
        self.description = ""

    @property
    def root(self):
        parent = self.parent
        while parent.parent:
            parent = parent.parent
        return parent

    @root.setter
    def root(self, value):
        self._root = value

    def get_params(self):
        pass

    @property
    def params(self):
        return self.get_params()

    def _activate(self, x, predict):
        raise NotImplementedError("Please implement activation function for " + self.name)

    def _derivative(self, y, delta=None):
        raise NotImplementedError("Please implement derivative function for " + self.name)

class ConvLayer(Layer):
    LayerTiming = Timing()

    def __init__(self, shape, stride=1, padding=0, parent=None):
        """
        :param shape:    shape[0] = shape of previous layer           c x h x w
                         shape[1] = shape of current layer's weight   f x c x h x w
        :param stride:   stride
        :param padding:  zero-padding
        """
        if parent is not None:
            parent = parent.root if parent.is_sub_layer else parent
            shape = parent.shape
        Layer.__init__(self, shape)
        self._stride, self._padding = stride, padding
        if len(shape) == 1:
            self.n_channels = self.n_filters = self.out_h = self.out_w = None
        else:
            self.feed_shape(shape)
        self.x_cache = self.x_col_cache = None
        self.inner_weight = None

    def feed_shape(self, shape):
        self._shape = shape
        self.n_channels, height, width = shape[0]
        self.n_filters, filter_height, filter_width = shape[1]
        full_height, full_width = width + 2 * self._padding, height + 2 * self._padding
        if (
            (full_height - filter_height) % self._stride != 0 or
            (full_width - filter_width) % self._stride != 0
        ):
            raise BuildLayerError(
                "Weight shape does not work, "
                "shape: {} - stride: {} - padding: {} not compatible with {}".format(
                    self._shape[1][1:], self._stride, self._padding, (height, width)
                ))
        self.out_h = int((height + 2 * self._padding - filter_height) / self._stride) + 1
        self.out_w = int((width + 2 * self._padding - filter_width) / self._stride) + 1

    @property
    def params(self):
        return self._shape, self._stride, self._padding

    @property
    def stride(self):
        return self._stride

    @property
    def padding(self):
        return self._padding

    def _activate(self, x, predict):
        raise NotImplementedError("Please implement activation function for " + self.name)

    def _derivative(self, y, delta=None):
        raise NotImplementedError("Please implement derivative function for " + self.name)

class ConvPoolLayer(ConvLayer):
    LayerTiming = Timing()

    def __init__(self, shape, stride=1, padding=0):
        """
        :param shape:    shape[0] = shape of previous layer           c x h x w
                         shape[1] = shape of pool window              c x ph x pw
        :param stride:   stride
        :param padding:  zero-padding
        """
        ConvLayer.__init__(self, shape, stride, padding)
        self._pool_cache = {}

    @property
    def params(self):
        return (self._shape[0], self._shape[1][1:]), self._stride, self._padding

    def feed_shape(self, shape):
        if len(shape[1]) == 2:
            shape = (shape[0], (shape[0][0], *shape[1]))
        ConvLayer.feed_shape(self, shape)

    @LayerTiming.timeit(level=1, prefix="[Core] ")
    def activate(self, x, w, bias=None, predict=False):
        return self._activate(x, w, bias, predict)

    def _activate(self, x, *args):
        raise NotImplementedError("Please implement activation function for " + self.name)

    def _derivative(self, y, *args):
        raise NotImplementedError("Please implement derivative function for " + self.name)

    @LayerTiming.timeit(level=1, prefix="[Core] ")
    def bp(self, y, w, prev_delta):
        return self._derivative(y, w, prev_delta)

# noinspection PyProtectedMember
class ConvMeta(type):
    def __new__(mcs, *args, **kwargs):
        name, bases, attr = args[:3]
        conv_layer, layer = bases

        def __init__(self, shape, stride=1, padding=0):
            conv_layer.__init__(self, shape, stride, padding)

        def _activate(self, x, w, bias, predict):
            self.x_cache, self.inner_weight = x, w
            n, n_channels, height, width = x.shape
            n_filters, _, filter_height, filter_width = w.shape

            p, sd = self._padding, self._stride
            x_padded = np.pad(x, ((0, 0), (0, 0), (p, p), (p, p)), mode='constant')

            height += 2 * p
            width += 2 * p

            shape = (n_channels, filter_height, filter_width, n, self.out_h, self.out_w)
            strides = (height * width, width, 1, n_channels * height * width, sd * width, sd)
            strides = x.itemsize * np.asarray(strides)
            x_cols = np.lib.stride_tricks.as_strided(x_padded, shape=shape, strides=strides).reshape(
                n_channels * filter_height * filter_width, n * self.out_h * self.out_w)
            self.x_col_cache = x_cols

            if bias is None:
                res = w.reshape(n_filters, -1).dot(x_cols)
            else:
                res = w.reshape(n_filters, -1).dot(x_cols) + bias.reshape(-1, 1)
            res.shape = (n_filters, n, self.out_h, self.out_w)
            return layer._activate(self, res.transpose(1, 0, 2, 3), predict)

        def _derivative(self, y, w, prev_delta):
            n = len(y)
            n_channels, height, width = self._shape[0]
            n_filters, filter_height, filter_width = self._shape[1]

            p, sd = self._padding, self._stride
            if isinstance(prev_delta, tuple):
                prev_delta = prev_delta[0]

            __derivative = self.LayerTiming.timeit(level=1, func_name="bp", cls_name=name, prefix="[Core] ")(
                layer._derivative)
            if self.is_fc_base:
                delta = __derivative(self, y) * prev_delta.dot(w.T).reshape(y.shape)
            else:
                delta = __derivative(self, y) * prev_delta

            dw = delta.transpose(1, 0, 2, 3).reshape(n_filters, -1).dot(
                self.x_col_cache.T).reshape(self.inner_weight.shape)
            db = np.sum(delta, axis=(0, 2, 3))

            n_filters, _, filter_height, filter_width = self.inner_weight.shape
            *_, out_h, out_w = delta.shape

            dx_padded = np.zeros((n, n_channels, height + 2 * p, width + 2 * p), dtype=np.float32)
            conv_bp(
                n, n_filters, out_h, out_w, dx_padded,
                filter_height, filter_width, sd, self.inner_weight, delta
            )
            dx = dx_padded[..., p:-p, p:-p] if p > 0 else dx_padded
            return dx, dw, db

        def activate(self, x, w, bias=None, predict=False):
            return self.LayerTiming.timeit(level=1, func_name="activate", cls_name=name, prefix="[Core] ")(
                _activate)(self, x, w, bias, predict)

        def bp(self, y, w, prev_delta):
            return self.LayerTiming.timeit(level=1, func_name="bp", cls_name=name, prefix="[Core] ")(
                _derivative)(self, y, w, prev_delta)

        for key, value in locals().items():
            if str(value).find("function") >= 0 or str(value).find("property"):
                attr[key] = value

        return type(name, bases, attr)

# noinspection PyProtectedMember
class ConvSubMeta(type):
    def __new__(mcs, *args, **kwargs):
        name, bases, attr = args[:3]
        conv_layer, sub_layer = bases

        def __init__(self, parent, shape, *_args, **_kwargs):
            conv_layer.__init__(self, None, parent=parent, **_kwargs)
            self.out_h, self.out_w = parent.out_h, parent.out_w
            sub_layer.__init__(self, parent, shape, *_args, **_kwargs)
            self._shape = ((shape[0][0], self.out_h, self.out_w), shape[0])
            if name == "ConvNorm":
                self.gamma = np.ones(self.n_filters, dtype=np.float32)
                self.beta = np.ones(self.n_filters, dtype=np.float32)
                self.init_optimizers()

        def _activate(self, x, predict):
            n, n_channels, height, width = x.shape
            out = sub_layer._activate(self, x.transpose(0, 2, 3, 1).reshape(-1, n_channels), predict)
            return out.reshape(n, height, width, n_channels).transpose(0, 3, 1, 2)

        def _derivative(self, y, w, delta=None):
            if self.is_fc_base:
                delta = delta.dot(w.T).reshape(y.shape)
            n, n_channels, height, width = delta.shape
            dx = sub_layer._derivative(self, y, delta.transpose(0, 2, 3, 1).reshape(-1, n_channels))
            return dx.reshape(n, height, width, n_channels).transpose(0, 3, 1, 2)

        # noinspection PyUnusedLocal
        def activate(self, x, w, bias=None, predict=False):
            return self.LayerTiming.timeit(level=1, func_name="activate", cls_name=name, prefix="[Core] ")(
                _activate)(self, x, predict)

        def bp(self, y, w, prev_delta):
            if isinstance(prev_delta, tuple):
                prev_delta = prev_delta[0]
            return self.LayerTiming.timeit(level=1, func_name="bp", cls_name=name, prefix="[Core] ")(
                _derivative)(self, y, w, prev_delta)

        @property
        def params(self):
            return sub_layer.get_params(self)

        for key, value in locals().items():
            if str(value).find("function") >= 0 or str(value).find("property"):
                attr[key] = value

        return type(name, bases, attr)

# Activation Layers

class Tanh(Layer):
    def _activate(self, x, predict):
        return np.tanh(x)

    def _derivative(self, y, delta=None):
        return 1 - y ** 2

class Sigmoid(Layer):
    def _activate(self, x, predict):
        return 1 / (1 + np.exp(-x))

    def _derivative(self, y, delta=None):
        return y * (1 - y)

class ELU(Layer):
    def _activate(self, x, predict):
        rs, mask = x.copy(), x < 0
        rs[mask] = np.exp(rs[mask]) - 1
        return rs

    def _derivative(self, y, delta=None):
        _rs, _indices = np.ones(y.shape), y < 0
        _rs[_indices] = y[_indices] + 1
        return _rs

class ReLU(Layer):
    def _activate(self, x, predict):
        return np.maximum(0, x)

    def _derivative(self, y, delta=None):
        return y > 0

class Softplus(Layer):
    def _activate(self, x, predict):
        return np.log(1 + np.exp(x))

    def _derivative(self, y, delta=None):
        return 1 - 1 / np.exp(y)

class Identical(Layer):
    def _activate(self, x, predict):
        return x

    def _derivative(self, y, delta=None):
        return 1

# Convolution Layers

class ConvTanh(ConvLayer, Tanh, metaclass=ConvMeta):
    pass

class ConvSigmoid(ConvLayer, Sigmoid, metaclass=ConvMeta):
    pass

class ConvELU(ConvLayer, ELU, metaclass=ConvMeta):
    pass

class ConvReLU(ConvLayer, ReLU, metaclass=ConvMeta):
    pass

class ConvSoftplus(ConvLayer, Softplus, metaclass=ConvMeta):
    pass

class ConvIdentical(ConvLayer, Identical, metaclass=ConvMeta):
    pass

# Pooling Layers

class MaxPool(ConvPoolLayer):
    def _activate(self, x, *args):
        self.x_cache = x
        sd = self._stride
        n, n_channels, height, width = x.shape
        # noinspection PyTupleAssignmentBalance
        _, pool_height, pool_width = self._shape[1]
        same_size = pool_height == pool_width == sd
        tiles = height % pool_height == 0 and width % pool_width == 0
        if same_size and tiles:
            x_reshaped = x.reshape(n, n_channels, int(height / pool_height), pool_height,
                                   int(width / pool_width), pool_width)
            self._pool_cache["x_reshaped"] = x_reshaped
            out = x_reshaped.max(axis=3).max(axis=4)
            self._pool_cache["method"] = "reshape"
        else:
            out = np.zeros((n, n_channels, self.out_h, self.out_w), dtype=np.float32)
            pos_cache = np.zeros((n, n_channels, self.out_h, self.out_w, 2), dtype=np.int32)
            max_pool(
                n, n_channels, self.out_h, self.out_w, x, out,
                pool_height, pool_width, sd, pos_cache
            )
            self._pool_cache["method"] = "original"
            self._pool_cache["pos_cache"] = pos_cache
        return out

    def _derivative(self, y, *args):
        w, prev_delta = args
        if isinstance(prev_delta, tuple):
            prev_delta = prev_delta[0]
        if self.is_fc_base:
            delta = prev_delta.dot(w.T).reshape(y.shape)
        else:
            delta = prev_delta
        method = self._pool_cache["method"]
        if method == "reshape":
            x_reshaped_cache = self._pool_cache["x_reshaped"]
            dx_reshaped = np.zeros_like(x_reshaped_cache)
            out_newaxis = y[..., None, :, None]
            mask = (x_reshaped_cache == out_newaxis)  # type: np.ndarray
            dout_newaxis = delta[..., None, :, None]
            dout_broadcast, _ = np.broadcast_arrays(dout_newaxis, dx_reshaped)
            dx_reshaped[mask] = dout_broadcast[mask]
            dx_reshaped /= np.sum(mask, axis=(3, 5), keepdims=True)
            dx = dx_reshaped.reshape(self.x_cache.shape)
        elif method == "original":
            sd = self._stride
            dx = np.zeros_like(self.x_cache)
            n, n_channels, *_ = self.x_cache.shape
            max_pool_bp(
                n, n_channels, self.out_h, self.out_w, sd, dx, delta, self._pool_cache["pos_cache"]
            )
        else:
            raise LayerError("Undefined pooling method '{}' found".format(method))
        return dx, None, None

# Special Layer

class Dropout(SubLayer):
    def __init__(self, parent, shape, keep_prob=0.5):
        if keep_prob < 0 or keep_prob >= 1:
            raise BuildLayerError("Probability of Dropout should be a positive float smaller than 1")
        SubLayer.__init__(self, parent, shape)
        self._mask = None
        self._prob = keep_prob
        self._prob_inv = 1 / keep_prob
        self.description = "(Keep prob: {})".format(keep_prob)

    def get_params(self):
        return self._prob,

    def _activate(self, x, predict):
        if not predict:
            # noinspection PyTypeChecker
            self._mask = np.random.binomial(
                [np.ones(x.shape)], self._prob
            )[0].astype(np.float32) * self._prob_inv
            return x * self._mask
        return x

    def _derivative(self, y, delta=None):
        return delta * self._mask

class Normalize(SubLayer):
    def __init__(self, parent, shape, lr=0.001, eps=1e-8, momentum=0.9, optimizers=None):
        SubLayer.__init__(self, parent, shape)
        self.sample_mean, self.sample_var = None, None
        self.running_mean, self.running_var = None, None
        self.x_cache, self.x_normalized_cache = None, None
        self._lr, self._eps = lr, eps
        if optimizers is None:
            self._g_optimizer, self._b_optimizer = Adam(self._lr), Adam(self._lr)
        else:
            self._g_optimizer, self._b_optimizer = optimizers
        self.gamma = np.ones(self.shape[1], dtype=np.float32)
        self.beta = np.ones(self.shape[1], dtype=np.float32)
        self._momentum = momentum
        self.init_optimizers()
        self.description = "(lr: {}, eps: {}, momentum: {}, optimizer: ({}, {}))".format(
            lr, eps, momentum, self._g_optimizer.name, self._b_optimizer.name
        )

    def get_params(self):
        return self._lr, self._eps, self._momentum, (self._g_optimizer.name, self._b_optimizer.name)

    @property
    def special_params(self):
        return {
            "gamma": self.gamma, "beta": self.beta,
            "running_mean": self.running_mean, "running_var": self.running_var,
            "_g_optimizer": self._g_optimizer, "_b_optimizer": self._b_optimizer
        }

    def init_optimizers(self):
        _opt_fac = OptFactory()
        if not isinstance(self._g_optimizer, Optimizer):
            self._g_optimizer = _opt_fac.get_optimizer_by_name(
                self._g_optimizer, None, self._lr, None
            )
        if not isinstance(self._b_optimizer, Optimizer):
            self._b_optimizer = _opt_fac.get_optimizer_by_name(
                self._b_optimizer, None, self._lr, None
            )
        self._g_optimizer.feed_variables([self.gamma])
        self._b_optimizer.feed_variables([self.beta])

    # noinspection PyTypeChecker
    def _activate(self, x, predict):
        if self.running_mean is None or self.running_var is None:
            self.running_mean = np.zeros(x.shape[1], dtype=np.float32)
            self.running_var = np.zeros(x.shape[1], dtype=np.float32)
        if not predict:
            self.sample_mean = np.mean(x, axis=0, keepdims=True)
            self.sample_var = np.var(x, axis=0, keepdims=True)
            x_normalized = (x - self.sample_mean) / np.sqrt(self.sample_var + self._eps)
            self.x_cache, self.x_normalized_cache = x, x_normalized
            out = self.gamma * x_normalized + self.beta
            self.running_mean = self._momentum * self.running_mean + (1 - self._momentum) * self.sample_mean
            self.running_var = self._momentum * self.running_var + (1 - self._momentum) * self.sample_var
        else:
            x_normalized = (x - self.running_mean) / np.sqrt(self.running_var + self._eps)
            out = self.gamma * x_normalized + self.beta
        return out

    def _derivative(self, y, delta=None):
        n, d = self.x_cache.shape
        dx_normalized = delta * self.gamma
        x_mu = self.x_cache - self.sample_mean
        sample_std_inv = 1.0 / np.sqrt(self.sample_var + self._eps)
        ds_var = -0.5 * np.sum(dx_normalized * x_mu, axis=0, keepdims=True) * sample_std_inv ** 3
        ds_mean = (-1.0 * np.sum(dx_normalized * sample_std_inv, axis=0, keepdims=True) - 2.0 *
                   ds_var * np.mean(x_mu, axis=0, keepdims=True))
        dx1 = dx_normalized * sample_std_inv
        dx2 = 2.0 / n * ds_var * x_mu
        dx = dx1 + dx2 + 1.0 / n * ds_mean
        dg = np.sum(delta * self.x_normalized_cache, axis=0)
        db = np.sum(delta, axis=0)
        self.gamma += self._g_optimizer.run(0, dg)
        self.beta += self._b_optimizer.run(0, db)
        self._g_optimizer.update()
        self._b_optimizer.update()
        return dx

class ConvDrop(ConvLayer, Dropout, metaclass=ConvSubMeta):
    pass

class ConvNorm(ConvLayer, Normalize, metaclass=ConvSubMeta):
    pass

# Cost Layer

class CostLayer(Layer):
    def __init__(self, shape, cost_function="MSE", transform=None):
        super(CostLayer, self).__init__(shape)
        self._available_cost_functions = {
            "MSE": CostLayer._mse,
            "SVM": CostLayer._svm,
            "CrossEntropy": CostLayer._cross_entropy
        }
        self._available_transform_functions = {
            "Softmax": CostLayer._softmax,
            "Sigmoid": CostLayer._sigmoid
        }
        if cost_function not in self._available_cost_functions:
            raise LayerError("Cost function '{}' not implemented".format(cost_function))
        self._cost_function_name = cost_function
        self._cost_function = self._available_cost_functions[cost_function]
        if transform is None and cost_function == "CrossEntropy":
            self._transform = "Softmax"
            self._transform_function = CostLayer._softmax
        else:
            self._transform = transform
            self._transform_function = self._available_transform_functions.get(transform, None)

    def __str__(self):
        return self._cost_function_name

    def _activate(self, x, predict):
        if self._transform_function is None:
            return x
        return self._transform_function(x)

    def _derivative(self, y, delta=None):
        raise LayerError("derivative function should not be called in CostLayer")

    def bp_first(self, y, y_pred):
        if self._cost_function_name == "CrossEntropy" and (
                self._transform == "Softmax" or self._transform == "Sigmoid"):
            return y - y_pred
        dy = -self._cost_function(y, y_pred)
        if self._transform_function is None:
            return dy
        return dy * self._transform_function(y_pred, diff=True)

    @property
    def calculate(self):
        return lambda y, y_pred: self._cost_function(y, y_pred, False)

    @property
    def cost_function(self):
        return self._cost_function_name

    @cost_function.setter
    def cost_function(self, value):
        if value not in self._available_cost_functions:
            raise LayerError("'{}' is not implemented".format(value))
        self._cost_function_name = value
        self._cost_function = self._available_cost_functions[value]

    def set_cost_function_derivative(self, func, name=None):
        name = "Custom Cost Function" if name is None else name
        self._cost_function_name = name
        self._cost_function = func

    # Transform Functions

    @staticmethod
    def _softmax(y, diff=False):
        if diff:
            return y * (1 - y)
        exp_y = CostLayer.safe_exp(y)
        return exp_y / np.sum(exp_y, axis=1, keepdims=True)

    @staticmethod
    def _sigmoid(y, diff=False):
        if diff:
            return y * (1 - y)
        return 1 / (1 + np.exp(-y))

    # Cost Functions

    @staticmethod
    def _mse(y, y_pred, diff=True):
        if diff:
            return -y + y_pred
        return 0.5 * np.average((y - y_pred) ** 2)

    @staticmethod
    def _svm(y, y_pred, diff=True):
        n, y = y_pred.shape[0], np.argmax(y, axis=1)
        correct_class_scores = y_pred[np.arange(n), y]
        margins = np.maximum(0, y_pred - correct_class_scores[..., None] + 1.0)
        margins[np.arange(n), y] = 0
        loss = np.sum(margins) / n
        num_pos = np.sum(margins > 0, axis=1)
        if not diff:
            return loss
        dx = np.zeros_like(y_pred)
        dx[margins > 0] = 1
        dx[np.arange(n), y] -= num_pos
        dx /= n
        return dx

    @staticmethod
    def _cross_entropy(y, y_pred, diff=True):
        if diff:
            return -y / y_pred + (1 - y) / (1 - y_pred)
        # noinspection PyTypeChecker
        return np.average(-y * np.log(np.maximum(y_pred, 1e-12)) - (1 - y) * np.log(np.maximum(1 - y_pred, 1e-12)))

# Factory

class LayerFactory:
    available_root_layers = {
        "Tanh": Tanh, "Sigmoid": Sigmoid,
        "ELU": ELU, "ReLU": ReLU, "Softplus": Softplus,
        "Identical": Identical,
        "ConvTanh": ConvTanh, "ConvSigmoid": ConvSigmoid,
        "ConvELU": ConvELU, "ConvReLU": ConvReLU, "ConvSoftplus": ConvSoftplus,
        "ConvIdentical": ConvIdentical,
        "MaxPool": MaxPool,
        "MSE": CostLayer, "SVM": CostLayer, "CrossEntropy": CostLayer
    }
    available_sub_layers = {
        "Dropout", "Normalize", "ConvNorm", "ConvDrop"
    }
    available_cost_functions = {
        "MSE", "SVM", "CrossEntropy"
    }
    available_special_layers = {
        "Dropout": Dropout,
        "Normalize": Normalize,
        "ConvDrop": ConvDrop,
        "ConvNorm": ConvNorm
    }
    special_layer_default_params = {
        "Dropout": (0.5, ),
        "Normalize": (0.001, 1e-8, 0.9),
        "ConvDrop": (0.5, ),
        "ConvNorm": (0.001, 1e-8, 0.9)
    }

    def get_root_layer_by_name(self, name, *args, **kwargs):
        if name not in self.available_sub_layers:
            if name in self.available_root_layers:
                if name in self.available_cost_functions:
                    kwargs["cost_function"] = name
                name = self.available_root_layers[name]
            else:
                raise BuildNetworkError("Undefined layer '{}' found".format(name))
            return name(*args, **kwargs)
        return None

    def get_layer_by_name(self, name, parent, current_dimension, *args, **kwargs):
        _layer = self.get_root_layer_by_name(name, *args, **kwargs)
        if _layer:
            return _layer, None
        _current, _next = parent.shape[1], current_dimension
        layer_param = self.special_layer_default_params[name]
        _layer = self.available_special_layers[name]
        if args or kwargs:
            _layer = _layer(parent, (_current, _next), *args, **kwargs)
        else:
            _layer = _layer(parent, (_current, _next), *layer_param)
        return _layer, (_current, _next)

max_stars_repo_name: wangjiangtao-NJPI/MachineLearning

max_issues_repo_name: wangjiangtao-NJPI/MachineLearning

max_forks_repo_name: wangjiangtao-NJPI/MachineLearning

reason: import numba


text:
from sympy import *

x, y, z = symbols("x,y,z")

func = raw_input("Ingrese la Funcion: \n")
pa = raw_input("Ingrese el punto: \n")
tol = 0.0001

func = sympify(func)

fdiff = diff(func, x)
ev = 1
resu = 0
ite = 0

if fdiff.subs(x, 0.4) < 1:

  while ev > 0.1:

    ite += 1
    fevala = func.subs(x, pa)
    eu = abs(pa - fevala)

    if eu < tol:
      ev = 0
      print(str(resu))

    elif ite == 100:
      print("No hay Punto Fijo")
      ev = 0

    else:
      pa = fevala
    resu = pa

else:
  print ("No converge en el punto dado")

ext: py

max_stars_repo_path: Python/Codigo Fuente/Sin Vista/Punto Fijo.py

max_stars_repo_name: Broomva/Numerical-Methods

max_issues_repo_name: Broomva/Numerical-Methods

max_forks_repo_name: Broomva/Numerical-Methods

reason: from sympy


text:
#!/usr/bin/env python

import numpy as np
import matplotlib.pyplot as plt
min = -1.0
max = 1.0
npts = 1000
inc = (max-min)/npts
t = np.arange(min, max, inc)
nplots = 3
fig, ax = plt.subplots(ncols=nplots,tight_layout=True,subplot_kw={"projection": "3d"})
for i in range(nplots):
    ax[i].set_xlabel("x")
    ax[i].set_ylabel("y")
    ax[i].set_zlabel("z")

a=1.0
b=1.0
c=2.0
x = a*np.cos(2*np.pi*t*c)
y = b*np.sin(2*np.pi*t*c)
ax[0].set_title("circle")
ax[0].plot3D(x, y, t, 'black')

a=1.0
b=3.0
c=2.0
x = a*np.cos(2*np.pi*t*c)
y = b*np.sin(2*np.pi*t*c)
ax[1].set_title("ellipse")
#ax[1].set_aspect("equal") not supported
ax[1].plot3D(x, y, t, 'blue')

x = t*a*np.cos(2*np.pi*t*c)
y = t*b*np.sin(2*np.pi*t*c)
ax[2].set_title("spiral")
ax[2].plot3D(x, y, t, 'green')
plt.show()

max_stars_repo_name: mrkraimer/testPvaPy

max_issues_repo_name: mrkraimer/testPvaPy

max_forks_repo_name: mrkraimer/testPvaPy

reason: import numpy


text:
import numpy as np
import matplotlib.pyplot as plt
from packages.mdp.mdp import MDP
from scipy.signal import convolve2d
from scipy.interpolate import interp1d
from itertools import groupby, cycle
from collections.abc import Iterable

class Gridworld(MDP):

    def __init__(self, R=None, discount=None, motionPatterns=None, motionMethod='stay', walls=None, shape=None,
                 circular=False):

        # initialize reward and discount through base class constructor
        MDP.__init__(self, R=R, discount=discount)

        # store motion specifications
        if motionMethod not in ('stay', 'renormalize'):
            raise ValueError(f"motion method '{motionMethod}' not available")
        self.motionMethod = motionMethod
        self.motionPatterns = motionPatterns
        self.circular = circular

        # construct map
        self.walls_2d = self.constructMap(walls, shape)

    @MDP.T.getter
    def T(self):
        T = MDP.T.fget(self)
        if T is None:
            T, self._Pmaps = self.constructTransitionModel(self.motionPatterns, self.motionMethod, self.circular)
            self._T = T
        return T

    @property
    def Pmaps(self):
        if self._Pmaps is None:
            MDP._T, self._Pmaps = self.constructTransitionModel(self.motionPatterns, self.motionMethod, self.circular)
            return self._Pmaps
        else:
            return None

    @property
    def shape(self):
        return self.walls_2d.shape if self.walls_2d is not None else None

    @property
    def size(self):
        return np.prod(self.shape)

    @property
    def nWalls(self):
        return self.walls_2d.sum() if self.walls_2d is not None else None

    @property
    def nStates(self):
        try:
            return self.size - self.nWalls
        except AttributeError:
            return None

    @property
    def nActions(self):
        return self.motionPatterns.shape[0] if self.motionPatterns is not None else None

    @property
    def statePositions(self):
        if self.walls_2d is None:
            return None
        else:
            # get 2d positions of the valid states
            gridIndices = np.c_[np.unravel_index(range(self.size), self.shape)]
            return gridIndices[~self.walls_2d.ravel(), :]

    @property
    def map(self):
        if self.walls_2d is None:
            return None
        else:
            # create map from wall specifications
            m = np.zeros_like(self.walls_2d, dtype=float)
            m[self.walls_2d] = np.nan
            m[~self.walls_2d] = range(self.nStates)

    @property
    def walls(self):
        # set walls as property without setter to forbid manual changes without calling the constructor
        return self.walls

    @property
    def motionPatterns(self):
        return self._motionPatterns

    @motionPatterns.setter
    def motionPatterns(self, x):
        if x is None:
            self._motionPatterns = None
        else:
            # convert to numpy array
            x = np.array(x)

            # assert that the pattern shapes are odd (current state assumed in center)
            try:
                assert np.allclose(x.sum(axis=(1, 2)), 1)
                assert np.all(x >= 0)
                assert np.all(np.mod(x.shape[1:3], 2))
            except AssertionError:
                raise ValueError('invalid motion patterns')

            # store patterns
            self._motionPatterns = x

    def constructMap(self, walls, shape):

        # assert that size information is provided
        if (walls is None) and (shape is None):
            raise ValueError('either shape or walls must be specified')

        # assert that size information is consistent
        if (walls is not None) and (shape is not None) and (np.asarray(walls).shape != shape):
            raise ValueError('inconsistent shape and wall specifications')

        # always convert wall specifications to 2d boolean map
        if walls is None:
            walls_2d = np.full(shape, False)
        else:
            walls = np.array(walls)
            if walls.ndim == 1:
                walls_2d = np.full(self.shape, False)
                walls_2d[np.unravel_index(walls, self.shape)] = True
            else:
                walls_2d = walls

        return walls_2d

    def constructTransitionModel(self, motionPatterns, invalidMovesMethod, circularFlag):

        # convert to numpy array
        motionPatterns = np.array(motionPatterns)

        # containers for the transition model / probability maps
        Pmaps = np.zeros((*self.shape, self.nStates, self.nActions))
        T = np.zeros((self.nStates, self.nStates, self.nActions))

        # boundary conditions for convolution of motion kernel
        boundary = 'wrap' if circularFlag else 'fill'

        # iterate over actions
        for a, pattern in enumerate(motionPatterns):

            # iterate over all possible current states
            for s in range(self.nStates):
                # create a map indicating the current state
                Pmap_before = np.zeros(self.shape)
                Pmap_before[tuple(self.statePositions[s, :])] = 1

                # apply motion pattern and add walls
                Pmap = convolve2d(Pmap_before, pattern, mode='same', boundary=boundary)
                Pmap[self.walls_2d] = np.nan

                # compute leftover probability mass
                validMass = np.nansum(Pmap)
                if invalidMovesMethod == 'renormalize':
                    if validMass == 0:
                        # if no motion possible for the given pattern, stay at current state
                        Pmap_after = Pmap_before
                        Pmap_after[self.walls_2d] = np.nan
                    else:
                        # otherwise, renormalize the leftover probability mass
                        Pmap_after = Pmap / validMass
                elif invalidMovesMethod == 'stay':
                    # shift the invalid probability mass to the current state
                    Pmap_after = Pmap
                    Pmap_after[tuple(self.statePositions[s, :])] += 1 - validMass

                # store current map and extract corresponding transition vector
                Pmaps[:, :, s, a] = Pmap_after
                T[:, s, a] = Pmap_after[~self.walls_2d]

        return T, Pmaps

    def plot(self, values=None, cmap=None):
        fig, ax = plt.subplots()

        if values is None:
            if cmap is None:
                cmap = plt.get_cmap('gray').reversed()
            plt.imshow(self.walls_2d, cmap=cmap)
        elif np.ndim(values) == 2:
            im = np.reshape(values, (*self.shape, -1))
            plt.imshow(im, cmap=cmap)
        else:
            if cmap is None:
                cmap = plt.get_cmap('viridis')
                cmap.set_bad(color='k')
            map = np.full(self.shape, np.nan)
            map[tuple(self.statePositions.T)] = values
            plt.imshow(map, cmap=cmap)
        ax.set_xticks(np.arange(self.shape[1])-0.5, minor=True)
        ax.set_xticks(np.arange(self.shape[1]), minor=False)
        ax.xaxis.grid(True, which='minor', linewidth=2)
        ax.xaxis.grid(False, which='major')
        ax.set_yticks(np.arange(self.shape[0])-0.5, minor=True)
        ax.set_yticks(np.arange(self.shape[0]), minor=False)
        ax.yaxis.grid(True, which='minor', linewidth=2)
        ax.yaxis.grid(False, which='major')

    def plotLabels(self):
        self.plot()
        for s in range(self.nStates):
            pos = self.statePositions[s, :]
            plt.text(pos[1], pos[0], s, horizontalalignment='center', verticalalignment='center')

    def plotTrajectories(self, trajectories, **kwargs):

        # always represent trajectories as nested lists
        if isinstance(trajectories, np.ndarray):
            trajectories = trajectories.tolist()
        if not isinstance(trajectories, Iterable):
            trajectories = [[trajectories]]
        elif not isinstance(trajectories[0], Iterable):
            trajectories = [trajectories]

        # plot the grid world
        self.plot()

        # iterate over all trajectories
        for trajectory in trajectories:
            # remove duplicate successive states
            trajectory = [s for s,_ in groupby(trajectory)]

            # get 2d positions of trajectory states
            positions = self.statePositions[trajectory, :]

            # depending on the trajectory length ...
            if len(trajectory) == 1:
                # plot a single point
                plt.plot(positions[:, 1], positions[:, 0], marker='o', **kwargs)
            elif len(trajectory) == 2:
                # plot a straight line
                plt.plot(positions[:, 1], positions[:, 0], **kwargs)
            else:
                # interpolate along the state positions

                # linear length along the points
                distance = np.cumsum(np.sqrt(np.sum(np.diff(positions, axis=0) ** 2, axis=1)))
                distance = np.insert(distance, 0, 0) / distance[-1]

                # interpolator object
                interpolator = interp1d(distance, positions, kind='quadratic', axis=0)

                # interpolation grid and interpolated values
                grid = np.linspace(0, 1, 10 * len(trajectory)-1)
                interpolation = interpolator(grid)

                # plot the interpolated line
                plt.plot(interpolation[:, 1], interpolation[:, 0], **kwargs)

    def plotArrows(self, angles, mags=None, length=0.5, head_width=0.2, head_length=0.2, facecolor='k', **kwargs):

        # if no magnitudes specified, use constant magnitude of 1
        if mags is None:
            mags = cycle([1])

        # convert degrees to radians
        angles_rad = np.deg2rad(angles)

        # plot all arrows
        for pos, ang, mag in zip(self.statePositions, angles_rad, mags):
            if mag == 0:
                continue
            dx = -length * mag * np.sin(ang)
            dy = -length * mag * np.cos(ang)
            plt.arrow(pos[1] - dx/2, pos[0] - dy/2, dx, dy,
                      head_width=mag*head_width, head_length=mag*head_length,
                      length_includes_head=True, facecolor=facecolor, **kwargs)

    def plotPolicy(self, pi, localpolicy2arrow, localpolicy2value=None, relMag=True, values=None, **kwargs):

        # convert local action distributions to arrows / values
        angles, mags = zip(*[localpolicy2arrow(a) for a in pi])
        if localpolicy2value:
            values = [localpolicy2value(a) for a in pi]

        # plot the grid world
        self.plot(values=values, cmap=kwargs.pop('cmap', None))

        # normalize magnitudes
        if relMag:
            mags = np.array(mags, dtype=float)
            mags /= mags.max()

        # plot arrows
        self.plotArrows(angles, mags, **kwargs)

reason: import numpy,from scipy


text:
import numpy

def affine_transform_inverse(np_transform):
    """Returns a new matrix that is the inverse of the given np_transform. Rather than simply using
        numpy.linalg.inv(), use this for homegeneous matrices to get the proper inverse.

    Args:
        np_transform(numpy.nparray): A 4x4 matrix representing an affine transform.

    Returns:
        (numpy.nparray): The inverse of the given homogeneous matrix.
    """
    rotation = np_transform[:3, :3]
    translation = np_transform[:3, 3]
    rotation_inv = numpy.linalg.inv(rotation)
    translation_inv = -1 * numpy.dot(rotation_inv, translation)
    result = numpy.identity(4)
    result[:3, :3] = rotation_inv
    result[:3, 3] = translation_inv.flatten()
    return result

max_stars_repo_name: SRI-IPS/general-tools

max_issues_repo_name: SRI-IPS/general-tools

max_forks_repo_name: SRI-IPS/general-tools


text:
import numpy as np
import matplotlib as mp

import matplotlib.pyplot as plt
import tensorflow as tf
import tensorflow.contrib.slim as slim
from tensorflow.examples.tutorials.mnist import input_data
import math
import os
import datetime

mnist = input_data.read_data_sets("MNIST_data/", one_hot=True)

tf.reset_default_graph()
batchSize = 128
learning_rate = 0.001
OUTPUT_DIR = "tmp"
logs_path = os.path.join(
            OUTPUT_DIR, datetime.datetime.now().strftime('%Y_%m_%d_%H_%M'))
if not os.path.exists(logs_path):
    os.makedirs(logs_path)

x = tf.placeholder(tf.float32, [None, 784],name="x-in")
true_y = tf.placeholder(tf.float32, [None, 10],name="y-in")
keep_prob = tf.placeholder("float")

x_image = tf.reshape(x,[-1,28,28,1])
tf.summary.image('input', x_image, 3)

hidden_1 = slim.conv2d(x_image,32,[3,3],scope="hidden_1")
pool_1 = slim.max_pool2d(hidden_1,[2,2],scope="pool_1")

hidden_2 = slim.conv2d(pool_1,64,[3,3],scope="hidden_2")
pool_2 = slim.max_pool2d(hidden_2,[2,2],scope="pool_2")

out_y = slim.fully_connected(slim.flatten(pool_2),10,activation_fn=tf.nn.softmax,scope="out_y")

with tf.name_scope('Accuracy'):
	correct_prediction = tf.equal(tf.argmax(out_y,1), tf.argmax(true_y,1))
	accuracy = tf.reduce_mean(tf.cast(correct_prediction, "float"))

# specify cost function
with tf.name_scope('cross_entropy'):
	cross_entropy = -tf.reduce_sum(true_y*tf.log(out_y))

train_step = tf.train.AdamOptimizer(learning_rate).minimize(cross_entropy,name="Adam")

# create a summary for our cost and accuracy
tf.summary.scalar("cost", cross_entropy)
tf.summary.scalar("accuracy", accuracy)

# merge all summaries into a single "operation" which we can execute in a session
summary_op = tf.summary.merge_all()
saver = tf.train.Saver()

with tf.Session() as sess:
	# variables need to be initialized before we can use them
	sess.run(tf.global_variables_initializer())
	# create log writer object
	writer = tf.summary.FileWriter(logs_path, graph=tf.get_default_graph())

	for i in range(1001):
	    batch = mnist.train.next_batch(batchSize)
	    _,summary = sess.run([train_step,summary_op], feed_dict={x:batch[0],true_y:batch[1], keep_prob:0.5})
	    writer.add_summary(summary,i)
	    if i % 100 == 0 and i != 0:
	        trainAccuracy = sess.run(accuracy, feed_dict={x:batch[0],true_y:batch[1], keep_prob:1.0})
	        saver.save(sess,"tmp/my-model",i)
	        print("step %d, training accuracy %g"%(i, trainAccuracy))

	testAccuracy = sess.run(accuracy, feed_dict={x:mnist.test.images[:256],true_y:mnist.test.labels[:256], keep_prob:1.0})
	print("test accuracy %g"%(testAccuracy))

	Use = mnist.test.images[10]
	plt.imshow(np.reshape(Use,[28,28]), interpolation="nearest", cmap="gray")
	r = 10
	print("Prediction: ", sess.run(
	    tf.argmax(out_y, 1), feed_dict={x: mnist.test.images[r:r + 1]}))

	print("Probs: ", sess.run(
	    out_y, feed_dict={x: mnist.test.images[r:r + 1]}))

max_stars_repo_name: PeterChenYijie/DeepLearningZeroToALL

max_issues_repo_name: PeterChenYijie/DeepLearning

max_forks_repo_name: PeterChenYijie/DeepLearning


text:
import click
from subprocess import call
from sys import stdout
from pathlib import Path
import tempfile
import shutil
import os

project_dir = Path(__file__).resolve().parents[2]

# git submodule init
# git submodule update
# cd dependencies/snap/examples/rolx && make

@click.group()
def cli():
    pass

@cli.command()
@click.option(
    "--rolx-bin",
    default=str(project_dir / "dependencies/snap/examples/rolx/testrolx"),
    type=click.Path(exists=True),
)
@click.option("--input", type=click.Path(exists=True), required=True)
@click.option("--output", type=click.Path(exists=False), required=True)
def role_mine(rolx_bin, input, output):
    folder = tempfile.mkdtemp()
    call([rolx_bin, "-i:" + str(Path(input).resolve())], stdout=stdout, cwd=folder)

    shutil.copyfile(folder + "/mappings.txt", output + "-mappings")
    shutil.copyfile(folder + "/v.txt", output + "-v")

import pandas as pd
from sklearn.decomposition import NMF
import numpy as np

@cli.command()
@click.option("--path", required=True)
@click.option("--n-components", default=16)
def role_factorize(path, n_components):
    mappings_file = path + "-mappings"
    v_file = path + "-v"
    assert os.path.exists(mappings_file) and os.path.exists(v_file)

    rolx_vec = pd.read_csv(v_file, header=None, sep=" ")
    rolx_mapping = pd.read_csv(mappings_file, header=None, skiprows=1, sep=" ")
    df = pd.concat([rolx_mapping, rolx_vec], axis=1).iloc[:, 1:]

    print("factorizing {}".format(v_file))
    nmf = NMF(n_components=16, solver="mu", beta_loss="kullback-leibler", max_iter=1000)

    X = df.values[:, 1:]
    W = nmf.fit_transform(X)
    H = nmf.components_

    print("factorized matrix G and F")
    pd.concat([df.iloc[:, 0], pd.DataFrame(W)], axis=1).to_csv(
        path + "-nmf-G.csv", index=False
    )
    pd.DataFrame(H).to_csv(path + "-nmf-F.csv", index=False)

from src.data.snap_import_user_projection import UnimodalUserProjection
from pyspark.sql import SparkSession, functions as F
import pandas as pd

@cli.command()
@click.option(
    "--enwiki-meta",
    type=click.Path(exists=True),
    default=str(project_dir / "data/processed/enwiki-meta-compact"),
)
@click.option(
    "--rolx-roles",
    type=click.Path(exists=True),
    default=str(project_dir / "data/processed/rolx-roles"),
)
@click.option(
    "--output",
    type=click.Path(exists=False),
    default=str(project_dir / "data/processed/role-features"),
)
def average_roles(enwiki_meta, rolx_roles, output):
    spark = SparkSession.builder.getOrCreate()
    model = UnimodalUserProjection(spark).extract(enwiki_meta).transform()

    role_df = spark.read.csv(
        "data/processed/rolx-roles",
        schema="user_id INT, role_id INT",
        sep="\t",
        comment="-",
        ignoreLeadingWhiteSpace=True,
    )

    bipartite = (
        spark.table("bipartite")
        .join(role_df, on="user_id", how="left")
        .na.fill({"role_id": -1})
    )

    article_roles = (
        bipartite.groupby("article_id", "edit_date")
        .pivot("role_id")
        .agg(F.count("user_id").alias("n_users"))
        .fillna(0)
    )

    totals = bipartite.groupby("article_id", "edit_date").agg(
        F.count("user_id").alias("deg")
    )

    normalized = article_roles.join(totals, on=["article_id", "edit_date"]).select(
        "article_id",
        "edit_date",
        *[
            (F.col(x) / F.col("deg")).alias("role_{}".format(x))
            for x in article_roles.columns[2:]
        ],
    )

    user_roles = (
        bipartite.join(normalized, on=["article_id", "edit_date"])
        .groupby("user_id")
        .agg(*[F.sum(x).alias(x) for x in normalized.columns[2:]])
    )

    # run row-wise normalization
    avg_roles = user_roles.toPandas()
    df = avg_roles.iloc[:, 1:]
    # https://stackoverflow.com/a/35679163
    df = pd.concat([avg_roles.iloc[:, 1], df.div(df.sum(axis=0), axis=1)], axis=1)
    df.to_csv(output)

if __name__ == "__main__":
    cli()


text:
# Unlabeled data generation
#
# Licensed under the BSD 3-Clause License
# Copyright (c) 2020, Yuriy Sverchkov

from typing import Collection, Iterable, Optional, Protocol
from abc import abstractmethod
from logging import getLogger

import numpy as np
from scipy.stats import ks_2samp, chi2_contingency

from generalizedtrees.constraints import Constraint, test
from generalizedtrees.features import FeatureSpec

logger = getLogger()

# Utility functions

def same_distribution(data_1, data_2, /, feature_spec, alpha):
    """
    Performs statistical test to determine if data are from the
    same distribution.

    data_1 and data_2 need to be numpy matrices
    """
    n_tests = 0
    min_p = 1.0
    for i in range(len(feature_spec)):

        if feature_spec[i] & FeatureSpec.DISCRETE:
            # Get frequencies.
            # Note: we're assuming that the union of values present in the samples is
            # the set of possible values. This is not all possible values that the
            # variable could originally take.
            v1, c1 = np.unique(data_1[:, i], return_counts=True)
            map1 = {v1[j]: c1[j] for j in range(len(v1))}

            v2, c2 = np.unique(data_2[:, i], return_counts=True)
            map2 = {v2[j]: c2[j] for j in range(len(v2))}

            values = np.union1d(v1, v2)
            k = len(values)

            # If only one value is present skip this test
            if k > 1:

                freq = [[map1[v] if v in map1 else 0 for v in values],
                        [map2[v] if v in map2 else 0 for v in values]]

                _, p, _, _ = chi2_contingency(observed=freq)

                if p < min_p:
                    min_p = p
                n_tests += 1

        else:
            # KS-test
            _, p = ks_2samp(data_1[:, i], data_2[:, i])
            if p < min_p:
                min_p = p
            n_tests += 1

    return min_p < alpha/n_tests

##################################
# Data Factory Learner Component #
##################################

# Interface definition
class DataFactoryLC(Protocol):

    feature_spec: Collection[FeatureSpec]

    def refit(self, data_matrix: np.ndarray) -> 'DataFactoryLC':
        return self

    @abstractmethod
    def generate(self, n: int, constraints: Iterable[Constraint] = ()) -> np.ndarray:
        raise NotImplementedError

# Base classes
class ConstraintFreeSamplingFactoryLC(DataFactoryLC):

    data_matrix: Optional[np.ndarray] = None

    max_attempts: int
    max_sample: int
    on_timeout: str = 'partial'

    def generate(self, n: int, constraints: Iterable[Constraint] = ()) -> np.ndarray:

        n = max(0, n)

        logger.debug(f'Drawing sample of size {n}')

        result = np.empty((n, self.data_matrix.shape[1]))

        if n < 1:
            return result

        else:

            i = 0
            oversample_prop = 1

            for _ in range(self.max_attempts):
                needed = n - i
                n_sampled = round(min(needed * oversample_prop, self.max_sample))
                sample = self._generate(n_sampled)
                accepted = test(constraints, sample)

                #if not all(accepted):
                #    logger.debug(f'Constraints: {constraints}')
                #    for j in np.nonzero(accepted):
                #        logger.debug(f'Rejected sample: {sample[j,:]}')

                n_acc = sum(accepted)
                n_fit = min(n_acc, needed)
                logger.debug(
                    f'Sampling loop: i={i}, n={n}. '
                    f'Needed {needed} samples, sampled {n_sampled}, accepted {n_acc}, fit {n_fit}')
                result[i:(i+n_fit),:] = sample[accepted,:][0:n_fit,:]
                i += n_fit
                if i >= n:
                    return result
                oversample_prop = max(oversample_prop, n_sampled / n_acc) if n_acc > 0 else np.inf

            logger.critical('Could not generate an acceptable sample within a reasonable time.')
            logger.debug(f'Failed to generate a sample for constraints: {constraints}')
            if self.on_timeout == 'partial':
                logger.warning('Returning partial sample')
                return result[1:i,:]
            elif self.on_timeout == 'dirty':
                logger.warning('Matrix with uninitialized elements')
                return result
            else: # on_timeout == 'raise':
                raise RuntimeError('Could not generate an acceptable sample within a reasonable time.')

    @abstractmethod
    def _generate(self, n: int) -> np.ndarray:
        raise NotImplementedError

# Implementations
class TrepanDataFactoryLC(ConstraintFreeSamplingFactoryLC):

    alpha: float

    rng: np.random.Generator

    def __init__(
        self,
        alpha: float = 0.05,
        max_attempts: int = 1000,
        max_sample: int = 100000,
        on_timeout: str = 'partial',
        rng: np.random.Generator = np.random.default_rng()) -> None:

        self.alpha = alpha
        self.max_attempts = max_attempts
        self.max_sample = max_sample
        self.on_timeout = on_timeout
        self.rng = rng

    def copy(self) -> 'TrepanDataFactoryLC':
        """
        Make a deep copy
        """
        clone = TrepanDataFactoryLC(
            alpha=self.alpha,
            max_attempts=self.max_attempts,
            max_sample=self.max_sample,
            on_timeout=self.on_timeout,
            rng=self.rng)
        clone.feature_spec = self.feature_spec

        return clone

    def refit(self, data_matrix: np.ndarray) -> 'DataFactoryLC':

        if len(data_matrix) < 1 or (
            self.data_matrix is not None and
            same_distribution(data_matrix, self.data_matrix, feature_spec=self.feature_spec, alpha = self.alpha)
        ):
            return self

        clone = self.copy()
        clone.data_matrix = data_matrix

        return clone

    def _generate(self, n: int) -> np.ndarray:

        # The Trepan generator independently generates the individual feature values.
        return np.column_stack([self._generate_feature(i, n) for i in range(self.data_matrix.shape[1])])

    def _generate_feature(self, i: int, n: int) -> np.ndarray:

        if self.feature_spec[i] is FeatureSpec.CONTINUOUS:
            # Sample from a KDE.
            # We use Generator and not RandomState but KDE implementations use RandomState
            # so it's more reliable to just implement the sampling here.
            return self.rng.normal(
                loc = self.rng.choice(self.data_matrix[:, i], size=n),
                scale = 1/np.sqrt(self.data_matrix.shape[0]),
                size = n)

        elif self.feature_spec[i] & FeatureSpec.DISCRETE:
            # Sample from the empirical distribution
            values, counts = np.unique(self.data_matrix[:,i], return_counts=True)
            return self.rng.choice(values, p=counts/self.data_matrix.shape[0], size=n)

        else:
            raise ValueError(f"I don't know how to handle feature spec {self.feature_spec[i]}")

class SmearingDataFactoryLC(ConstraintFreeSamplingFactoryLC):

    p_alt: float
    rng: np.random.Generator

    def __init__(
        self,
        p_alt: float = 0.5,
        max_attempts: int = 1000,
        max_sample: int = 100000,
        on_timeout: str = 'partial',
        rng: np.random.Generator = np.random.default_rng()
    ) -> None:

        self.p_alt = p_alt
        self.max_attempts = max_attempts
        self.max_sample = max_sample
        self.on_timeout = on_timeout
        self.rng = rng

    def refit(self, data_matrix: np.ndarray) -> DataFactoryLC:
        if self.data_matrix is None:
            self.data_matrix = data_matrix

        return self

    def _generate(self, n) -> np.ndarray:

        m, d = self.data_matrix.shape

        matrix_a = self.data_matrix[self.rng.integers(m, size=n, dtype=np.intp),:]
        matrix_b = self.data_matrix[self.rng.integers(m, size=n, dtype=np.intp),:]
        replacemask = self.rng.choice([False, True], size=n*d, p=[1-self.p_alt, self.p_alt]).reshape(n,d)
        matrix_a[replacemask] = matrix_b[replacemask]

        return matrix_a


text:
# TODO: needs to be standardized with BaseTask

import pandas as pd
import numpy as np
import os
from tophat.constants import SEED
from config.common import *

from gensim.models.doc2vec import TaggedDocument, FAST_VERSION
from gensim.models import Doc2Vec

import fastavro as avro
import itertools as it
from tophat.schemas import factors_avro
from typing import Iterator, Dict, Any, Callable, Tuple

import multiprocessing
n_cores = multiprocessing.cpu_count()

assert FAST_VERSION > -1, "This will be painfully slow otherwise"

def fit_interactions(interactions_df: pd.DataFrame,
                     user_col: str='ops_user_id',
                     item_col: str='ops_product_id',
                     emb_dim: int=16,
                     ):
    # Note: gensim requires the tag(s) and words to be str for silly reasons
    # For now, the only doctag will be the user_id,
    # but you can imagine user features being tags
    interactions_df[item_col+'str'] = interactions_df[item_col].astype(str)

    docs = interactions_df\
        .groupby(user_col)[item_col+'str'].apply(list).reset_index()\
        .apply(lambda row: TaggedDocument(row[item_col+'str'],
                                          tags=[str(row[user_col])]),
               axis=1)\
        .tolist()

    interactions_df.drop(item_col+'str', axis=1, inplace=True)

    model = Doc2Vec(docs,
                    size=emb_dim,
                    dm=1,
                    window=4, min_count=5, negative=5,
                    iter=10,
                    workers=1,  # should be `n_cores`, but see issue gensim#336
                    seed=SEED,
                    )

    return model

def model_to_dfs(d2v_model) -> Tuple[pd.DataFrame, pd.DataFrame]:
    dv = d2v_model.docvecs
    wv = d2v_model.wv

    user_keys = list(dv.doctags.keys())
    user_df = pd.DataFrame(index=pd.Index(user_keys, name='id', dtype=str))
    user_df.reset_index(inplace=True)  # downstream expects `id` column
    user_df['factors'] = dv[user_keys].astype(np.float32).tolist()
    user_df['bias'] = 0.

    item_keys = list(wv.vocab.keys())
    item_df = pd.DataFrame(index=pd.Index(item_keys, name='id', dtype=str))
    item_df.reset_index(inplace=True)
    item_df['factors'] = [wv.word_vec(w) for w in item_keys]
    item_df['bias'] = 0.

    return user_df, item_df

def rec_generator(keys: Iterator,
                  get_vec_fn: Callable,
                  partition: int = 0,
                  n_partitions: int = 1) -> Iterator[Dict[str, Any]]:
    sliced_keys = it.islice(keys, partition, None, n_partitions)
    for k in sliced_keys:
        record = {
            'id': str(k),
            'factors': list(get_vec_fn(k)),
            'bias': 0.,
        }

        yield record

dv_rec_generator = lambda dv, **kwargs: rec_generator(
    keys=dv.doctags.keys(), get_vec_fn=lambda doctag: dv[doctag], **kwargs
)

wv_rec_generator = lambda wv, **kwargs: rec_generator(
    keys=wv.vocab.keys(), get_vec_fn=wv.word_vec, **kwargs
)

def export(d2v_model: Doc2Vec, dir_export: str):
    # Record generators
    user_rec_gen = dv_rec_generator(d2v_model.docvecs)
    item_rec_gen = wv_rec_generator(d2v_model.wv)

    with open(os.path.join(dir_export, 'user_docvecs.avro'), 'wb') as f_out:
        avro.writer(f_out, factors_avro, user_rec_gen, codec='snappy')

    with open(os.path.join(dir_export, 'item_wordvecs.avro'), 'wb') as f_out:
        avro.writer(f_out, factors_avro, item_rec_gen, codec='snappy')

max_stars_repo_name: JasonTam/tophat

max_issues_repo_name: JasonTam/tophat

max_forks_repo_name: JasonTam/tophat


text:
"""
-------------------------------------------------------------------------
Example script for the Fourier transform usage (BCs)
-------------------------------------------------------------------------
This is an example script for differences between symmetric BCs and smoothed
or finite dots in artificial generated images.

.. moduleauthor:: Michael Barbier <michael@unam.bilkent.edu.tr>
"""

# -------------------------------------------------------------------------
# Import libraries
# -------------------------------------------------------------------------
import cv2 as cv
import numpy as np
from matplotlib import pyplot as plt
from skimage import io, filters

# -------------------------------------------------------------------------
# Functions
# -------------------------------------------------------------------------
def generate_lattice_square(a, n, dot_radius, shift_x, shift_y):
    """
    Generates an image with dots on a period square lattice with lattice parameter a
    and dot size given by the dot_radius.

    :param a: The lattice length parameter
    :param n: Number of the repeating the basic lattice in both x and y direction
    :param dot_radius: radius of the plotted dots (should be non-negative < a/2), if zero, single pixels are used
    :param shift_y: Shift the dot in x direction (should be non-negative < a)
    :param shift_x: Shift the dot in y direction (should be non-negative < a)
    :return: The generated img
    """
    sz = a * n
    shape = (sz, sz)

    img = np.zeros(shape, np.dtype('float'))
    for ix in range(0, int(sz/a)):
        for iy in range(0, int(sz/a)):
            if dot_radius == 0:
                img[ix*a+shift_x, iy*a+shift_y] = 1.0
            else:
                cv.circle(img, (ix*a+shift_x, iy*a+shift_y), radius=dot_radius, color=1, thickness=cv.FILLED)
    return img

def show_fft(img):
    """
    Shows the Fourier spectrum together with the original image.
    The phase, magnitude and the magnitude in decibel is shown.
    Remark that the log of the magnitude (in db) is not defined if it is zero (resulting in wrong white pixels)

    :param img: input image
    :return: Nothing
    """
    plt.subplot(221), plt.imshow(img, cmap='gray')
    plt.title('Input Image'), plt.xticks([]), plt.yticks([])

    f = np.fft.fft2(img)
    fshift = np.fft.fftshift(f)
    magnitude_spectrum = np.abs(fshift)
    magnitude_spectrum_db = 20 * np.log(np.abs(fshift))
    phase_spectrum = np.angle(fshift)

    plt.subplot(222), plt.imshow(phase_spectrum, cmap='gray')
    plt.title('Phase Spectrum'), plt.xticks([]), plt.yticks([])

    plt.subplot(223), plt.imshow(magnitude_spectrum, cmap='gray')
    plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])

    plt.subplot(224), plt.imshow(magnitude_spectrum_db, cmap='gray')
    plt.title('Magnitude Spectrum (in dB)'), plt.xticks([]), plt.yticks([])
    plt.show()

# -------------------------------------------------------------------------
# Script
# -------------------------------------------------------------------------

# At a region of square lattice
x0 = 0
y0 = 0
a = 16
n = 4
sz = a*n
size_crop = sz-1
shift_x = round(a/2)
shift_y = round(a/2)

# show the original image with Fourier spectra
dot_radius = 0
img_ori = generate_lattice_square(a, n, dot_radius, shift_x, shift_y)
show_fft(img_ori)

# show the image with asymmetric BCs due to crop: the Fourier spectrum is not perfect
img_asymmetric_1 = img_ori[y0:y0+sz-1, x0:x0+sz-1]
img_asymmetric_2 = img_ori[y0:y0+sz-2, x0:x0+sz-2]
img_asymmetric_3 = img_ori[y0:y0+sz-3, x0:x0+sz-3]
show_fft(img_asymmetric_1)
show_fft(img_asymmetric_2)
show_fft(img_asymmetric_3)

# Use of finite size dots instead of single pixel dots
dot_radius = 3
img_ori_4 = generate_lattice_square(a, n, dot_radius, shift_x, shift_y)
show_fft(img_ori_4)

# Smooth the image with a Gaussian
sigma = 2
img_smooth = filters.gaussian(img_ori_4, sigma=sigma)
show_fft(img_smooth)


text:
#!/usr/bin/python -O

import sys
import math
import numpy
import radermacher
from pyami import mem, imagefun
from scipy import ndimage
from appionlib.apImage import imagefile, imagenorm

def radonTransform(image, stepsize=2, maskrad=None):
	"""
	computes Radon transform of image
	"""
	radonlist = []
	if maskrad is None:
		maskrad = image.shape[0]/2
	blackcircle = imagefun.filled_circle(image.shape, maskrad)
	mask = 1 - blackcircle
	maskline = mask.sum(axis=0) + 1

	nsteps = int(math.ceil(180/stepsize))
	for i in range(nsteps):
		rotated = ndimage.rotate(image, -i*stepsize, reshape=False, order=1)
		rotated = mask*rotated
		line = rotated.sum(axis=0)
		radonlist.append(line/maskline)
	radon = numpy.array(radonlist)
	#radonlr = numpy.fliplr(radon)
	#radon = numpy.vstack((radon, radonlr))
	return radon

def createDots(angle=0.0, noiselevel=3.0, shape=(128,128), shift=None):
	a = numpy.zeros(shape, dtype=numpy.float64)
	a30 = int(min(shape)*.30) #25%
	a40 = int(2*min(shape)/5.0) #40%
	a50 = int(min(shape)/2.0) #50%
	a60 = int(3*min(shape)/5.0) #60%
	a70 = int(min(shape)*.70) #75%
	a[a30:a40,a30:a40] = 1
	a[a50:a60,a30:a40] = 1
	a[a60:a70,a40:a50] = 1
	a[a40:a50,a50:a60] = 1
	a[a30:a40,a60:a70] = 1
	a[a60:a70,a60:a70] = 1
	b = a
	b = ndimage.rotate(b, angle, reshape=False, order=1)
	if shift is not None:
		b = ndimage.shift(b, shift=shift, mode='wrap', order=1)
	bnoise = b + noiselevel*numpy.random.random(shape)
	bnoise = ndimage.median_filter(bnoise, size=2)
	bnoise = imagenorm.normStdev(bnoise)

	return bnoise

if __name__ == "__main__":
	u = mem.used()
	stepsize = 1.0
	shape = (256,256)
	noiselevel = 2.0

	a = createDots(angle=0., noiselevel=noiselevel, shape=shape)
	b = createDots(angle=-30., noiselevel=noiselevel, shape=shape, shift=(3,4))
	imagefile.arrayToJpeg(a, "imagea.jpg")
	imagefile.arrayToJpeg(b, "imageb.jpg")

	ar = radonTransform(a, stepsize)
	br = radonTransform(b, stepsize)
	imagefile.arrayToJpeg(ar, "radona.jpg")
	imagefile.arrayToJpeg(br, "radonb.jpg")

	#print "%.3f"%(a[1,0])
	#print "%.3f"%(b[1,0])
	for radius in range(0,11,2):
		e = radermacher.radonshift(ar,br,radius)
		e = e / (ar.std()*br.std()) # result is not normalized properly
		val = e.argmax()
		y = val%e.shape[1]
		x = int(val/e.shape[1])
		anglelist = radermacher.getAngles(radius)
		shiftangle = anglelist[y]
		xshift = radius*math.sin(shiftangle)
		yshift = radius*math.cos(shiftangle)
		rotangle = stepsize*x
		print "Radius %02d, Max value %.16f at rot=%d, shift=%.1f,%.1f"%(radius, e.max(), rotangle, xshift, yshift)
		imagefile.arrayToJpeg(e, "radoncc%02d.jpg"%(radius), msg=False)
	print "mem used = %.1f MB"%((mem.used()-u)/1024.)
	print "done"
	#sys.exit(1)


text:
#!/usr/bin/python
# -*- coding: utf-8 -*-
import torch
import torch.nn as nn
import torch.nn.functional as F
import math
import numpy as np
import torch.backends.cudnn as cudnn

def point_form(boxes):
    """ Convert prior_boxes to (xmin, ymin, xmax, ymax)
    representation for comparison to point form ground truth data.
    Args:
        boxes: (tensor) center-size default boxes from priorbox layers.
    Return:
        boxes: (tensor) Converted xmin, ymin, xmax, ymax form of boxes.
    """

    return torch.cat((boxes[:, :2] - boxes[:, 2:] / 2, boxes[:, :2]
                     + boxes[:, 2:] / 2), 1)  # xmin, ymin
                                              # xmax, ymax

def center_size(boxes):
    """ Convert prior_boxes to (cx, cy, w, h)
    representation for comparison to center-size form ground truth data.
    Args:
        boxes: (tensor) point_form boxes
    Return:
        boxes: (tensor) Converted xmin, ymin, xmax, ymax form of boxes.
    """

    return torch.cat(((boxes[:, 2:] + boxes[:, :2]) / 2, boxes[:, 2:]
                     - boxes[:, :2]), 1)  # cx, cy
                                          # w, h

def intersect(box_a, box_b):
    """ We resize both tensors to [A,B,2] without new malloc:
    [A,2] -> [A,1,2] -> [A,B,2]
    [B,2] -> [1,B,2] -> [A,B,2]
    Then we compute the area of intersect between box_a and box_b.
    Args:
      box_a: (tensor) bounding boxes, Shape: [A,4].
      box_b: (tensor) bounding boxes, Shape: [B,4].
    Return:
      (tensor) intersection area, Shape: [A,B].
    """

    A = box_a.size(0)
    B = box_b.size(0)
    max_xy = torch.min(box_a[:, 2:].unsqueeze(1).expand(A, B, 2),
                       box_b[:, 2:].unsqueeze(0).expand(A, B, 2))
    min_xy = torch.max(box_a[:, :2].unsqueeze(1).expand(A, B, 2),
                       box_b[:, :2].unsqueeze(0).expand(A, B, 2))
    inter = torch.clamp(max_xy - min_xy, min=0)
    return inter[:, :, 0] * inter[:, :, 1]

def jaccard(box_a, box_b):
    '''Compute the jaccard overlap of two sets of boxes.  The jaccard overlap\n    is simply the intersection over union of two boxes.  Here we operate on\n    ground truth boxes and default boxes.\n    E.g.:\n        A \xe2\x88\xa9 B / A \xe2\x88\xaa B = A \xe2\x88\xa9 B / (area(A) + area(B) - A \xe2\x88\xa9 B)\n    Args:\n        box_a: (tensor) Ground truth bounding boxes, Shape: [num_objects,4]\n        box_b: (tensor) Prior boxes from priorbox layers, Shape: [num_priors,4]\n    Return:\n        jaccard overlap: (tensor) Shape: [box_a.size(0), box_b.size(0)]\n    '''

    inter = intersect(box_a, box_b)
    area_a = ((box_a[:, 2] - box_a[:, 0]) * (box_a[:, 3] - box_a[:,
              1])).unsqueeze(1).expand_as(inter)  # [A,B]
    area_b = ((box_b[:, 2] - box_b[:, 0]) * (box_b[:, 3] - box_b[:,
              1])).unsqueeze(0).expand_as(inter)  # [A,B]
    union = area_a + area_b - inter
    return inter / union  # [A,B]

def match(
    truths,
    priors,
    labels,
    loc_t,
    conf_t,
    overlap_t,
    idx,
    ):
    """Match each prior box with the ground truth box of the highest jaccard
    overlap, encode the bounding boxes, then return the matched indices
    corresponding to both confidence and location preds.
    Args:
        truths: (tensor) Ground truth boxes, Shape: [num_obj, 4].
        priors: (tensor) Prior boxes from priorbox layers, Shape: [n_priors, 4].
        labels: (tensor) All the class labels for the image, Shape: [num_obj].
        loc_t: (tensor) Tensor to be filled w/ endcoded location targets.
        conf_t: (tensor) Tensor to be filled w/ matched indices for conf preds.
        overlap_t: (tensor) Tensor to be filled w/ iou score for each priors.
        overlap_t: (tensor) Tensor to be filled w/ match object idx for each priors.
        idx: (int) current batch index
    """

    # jaccard index

    overlaps = jaccard(truths, point_form(priors))

    # [1,num_priors] best ground truth for each prior

    (best_truth_overlap, best_truth_idx) = overlaps.max(0)

    # (Bipartite Matching)
    # [1,num_objects] best prior for each ground truth

    (best_prior_overlap, best_prior_idx) = overlaps.max(1)
    best_truth_overlap.index_fill_(0, best_prior_idx, 1)  # ensure best prior
    for j in range(best_prior_idx.size(0)):
        best_truth_idx[best_prior_idx[j]] = j

    overlap_t[idx] = best_truth_overlap  # [num_priors] jaccord for each prior
    conf_t[idx] = labels[best_truth_idx]  # [num_priors] top class label for each prior
    loc_t[idx] = truths[best_truth_idx]  # Shape: [num_priors,4]

def mutual_match(
    truths,
    priors,
    regress,
    classif,
    labels,
    loc_t,
    conf_t,
    overlap_t,
    pred_t,
    idx,
    ):
    """Classify to regress and regress to classify, Mutual Match for label assignement.
    Args:
        truths: (tensor) Ground truth boxes, Shape: [num_obj, 4].
        priors: (tensor) Prior boxes from priorbox layers, Shape: [num_priors, 4].
        regress: (tensor) Regression prediction, Shape: [num_priors, 4].
        classif: (tensor) Classification prediction, Shape: [num_priors, num_classes].
        labels: (tensor) All the class labels for the image, Shape: [num_obj].
        loc_t: (tensor) Tensor to be filled w/ endcoded location targets.
        conf_t: (tensor) Tensor to be filled w/ matched indices for conf preds.
        overlap_t: (tensor) Tensor to be filled w/ iou score for each priors.
        overlap_t: (tensor) Tensor to be filled w/ pred score for each priors.
        idx: (int) current batch index
    """

    num_obj = truths.size()[0]
    acr_overlaps = jaccard(truths, point_form(priors))
    reg_overlaps = jaccard(truths, decode(regress, priors))
    pred_classifs = classif.sigmoid().t()[labels - 1, :]
    sigma = 2.0
    pred_classifs = acr_overlaps ** ((sigma - pred_classifs) / sigma)

    # pred_classifs = sigma * acr_overlaps / (sigma-pred_classifs)

    # ## at least 1 anchor per object ###

    acr_overlaps[torch.arange(num_obj), acr_overlaps.max(1)[1]] = 1.0
    reg_overlaps[torch.arange(num_obj), reg_overlaps.max(1)[1]] = 1.0
    pred_classifs[torch.arange(num_obj), pred_classifs.max(1)[1]] = 1.0
    acr_overlaps[acr_overlaps != acr_overlaps.max(dim=0,
                 keepdim=True)[0]] = 0.0
    reg_overlaps[reg_overlaps != reg_overlaps.max(dim=0,
                 keepdim=True)[0]] = 0.0
    pred_classifs[pred_classifs != pred_classifs.max(dim=0,
                  keepdim=True)[0]] = 0.0

    # ## assign pos and ign nums according to acr_overlaps ###

    for (reg_overlap, pred_classif, acr_overlap) in zip(reg_overlaps,
            pred_classifs, acr_overlaps):
        num_ign = (acr_overlap >= 0.4).sum()
        num_pos = (acr_overlap >= 0.5).sum()

        ign_mask = torch.topk(reg_overlap, num_ign, largest=True)[1]
        pos_mask = torch.topk(reg_overlap, num_pos, largest=True)[1]
        reg_overlap[ign_mask] = 2.0
        reg_overlap[pos_mask] = 3.0

        pos_mask = torch.topk(pred_classif, num_pos, largest=True)[1]
        pred_classif[pos_mask] = 3.0

    # ## for classification ###

    (best_truth_overlap, best_truth_idx) = reg_overlaps.max(dim=0)
    overlap_t[idx] = best_truth_overlap  # [num_priors] jaccord for each prior
    conf_t[idx] = labels[best_truth_idx]  # [num_priors] top class label for each prior

    # ## for regression ###

    (best_truth_overlap, best_truth_idx) = pred_classifs.max(dim=0)
    pred_t[idx] = best_truth_overlap  # [num_priors] jaccord for each prior
    loc_t[idx] = truths[best_truth_idx]  # Shape: [num_priors,4]

def encode(matched, priors, variances=[0.1, 0.2]):
    """Encode the variances from the priorbox layers into the ground truth boxes
    we have matched (based on jaccard overlap) with the prior boxes.
    Args:
        matched: (tensor) Coords of ground truth for each prior in point-form
            Shape: [num_priors, 4].
        priors: (tensor) Prior boxes in center-offset form
            Shape: [num_priors,4].
        variances: (list[float]) Variances of priorboxes
    Return:
        encoded boxes (tensor), Shape: [num_priors, 4]
    """

    # dist b/t match center and prior's center

    g_cxcy = (matched[:, :2] + matched[:, 2:]) / 2 - priors[:, :2]

    # encode variance

    g_cxcy /= variances[0] * priors[:, 2:]

    # match wh / prior wh

    g_wh = (matched[:, 2:] - matched[:, :2]) / priors[:, 2:]
    g_wh = torch.log(g_wh) / variances[1]

    # return target for smooth_l1_loss

    targets = torch.cat([g_cxcy, g_wh], 1)  # [num_priors,4]
    return targets

def decode(loc, priors, variances=[0.1, 0.2]):
    """Decode locations from predictions using priors to undo
    the encoding we did for offset regression at train time.
    Args:
        loc (tensor): location predictions for loc layers,
            Shape: [num_priors,4]
        priors (tensor): Prior boxes in center-offset form.
            Shape: [num_priors,4].
        variances: (list[float]) Variances of priorboxes
    Return:
        decoded bounding box predictions
    """

    boxes = torch.cat((priors[:, :2] + loc[:, :2] * variances[0]
                      * priors[:, 2:], priors[:, 2:] * torch.exp(loc[:,
                      2:] * variances[1])), 1)
    boxes[:, :2] -= boxes[:, 2:] / 2
    boxes[:, 2:] += boxes[:, :2]
    return boxes

def matrix_iou(a, b):
    """
    return iou of a and b, numpy version for data augenmentation
    """

    lt = np.maximum(a[:, np.newaxis, :2], b[:, :2])
    rb = np.minimum(a[:, np.newaxis, 2:], b[:, 2:])

    area_i = np.prod(rb - lt, axis=2) * (lt < rb).all(axis=2)
    area_a = np.prod(a[:, 2:] - a[:, :2], axis=1)
    area_b = np.prod(b[:, 2:] - b[:, :2], axis=1)
    return area_i / (area_a[:, np.newaxis] + area_b - area_i)

def nms(dets, thresh=0.5):
    """Python version Non maximun suppression.
        See: https://github.com/rbgirshick/fast-rcnn/blob/master/lib/utils/nms.py
    Args:
        dets (numpy arrays): detected bounding boxes
        thresh (float): iou threshold
        mode (string): iou or ciou
    """

    x1 = dets[:, 0]
    y1 = dets[:, 1]
    x2 = dets[:, 2]
    y2 = dets[:, 3]
    scores = dets[:, 4]

    areas = (x2 - x1 + 1) * (y2 - y1 + 1)
    order = scores.argsort()[::-1]

    keep = []
    while order.size > 0:
        i = order[0]
        keep.append(i)
        xx1 = np.maximum(x1[i], x1[order[1:]])
        yy1 = np.maximum(y1[i], y1[order[1:]])
        xx2 = np.minimum(x2[i], x2[order[1:]])
        yy2 = np.minimum(y2[i], y2[order[1:]])

        w = np.maximum(0.0, xx2 - xx1 + 1)
        h = np.maximum(0.0, yy2 - yy1 + 1)
        inter = w * h
        ovr = inter / (areas[i] + areas[order[1:]] - inter)
        inds = np.where(ovr <= thresh)[0]
        order = order[inds + 1]

    return keep

max_stars_repo_name: XrosLiang/MutualGuide

max_issues_repo_name: XrosLiang/MutualGuide

max_forks_repo_name: XrosLiang/MutualGuide


text:
import torch.nn as nn
import math
import torch
import torchvision.models as torchmodels
import re
from torch.autograd import Variable
import torch.nn.functional as F
import numpy as np
import time
import random
import torch.nn.init as torchinit
import math
from torch.nn import init, Parameter
import copy
import torchvision
from models import base_single
import utils
from utils import *
import tqdm
import gc

def arch_is_picked(n, policy):
    n_split = n.split('.')

    if len(n_split) >= 5:

        num_layer = int(n_split[1])
        action = int(n_split[2])

        if policy.shape[0] == 1:
            return policy[0][num_layer][action].item() > 0.
        else:
            tmp = torch.mean(policy, 0)
            return tmp[num_layer][action].item() > 0.
    else:
        return True

class MetaGraph(nn.Module):
    '''
    For EWC related code, we modified the code from https://github.com/GMvandeVen/continual-learning
    '''

    def __init__(self):
        super().__init__()
        # -SI:
        self.si_c = 0           #-> hyperparam: how strong to weigh SI-loss ("regularisation strength")
        self.epsilon = 0.1      #-> dampening parameter: bounds 'omega' when squared parameter-change goes to 0

        # -EWC:
        self.gamma = 1.         #-> hyperparam (online EWC): decay-term for old tasks' contribution to quadratic term
        self.online = True      #-> "online" (=single quadratic term) or "offline" (=quadratic term per task) EWC
        self.emp_FI = False     #-> if True, use provided labels to calculate FI ("empirical FI"); else predicted labels
        self.EWC_task_count = 0 #-> keeps track of number of quadratic loss terms (for "offline EWC")

    #----------------- EWC-specifc functions -----------------#

    def estimate_fisher_update_by_policy(self, inputs, targets, policy, allowed_classes=None, collate_fn=None):
        '''After completing training on a task, estimate diagonal of Fisher Information matrix.
        [dataset]:          <DataSet> to be used to estimate FI-matrix
        [allowed_classes]:  <list> with class-indeces of 'allowed' or 'active' classes'''

        # Prepare <dict> to store estimated Fisher Information matrix
        est_fisher_info = {}
        for n, p in self.named_parameters():
            if p.requires_grad:
                n = n.replace('.', '__')
                est_fisher_info[n] = p.detach().clone().zero_()

        # Set model to evaluation mode
        mode = self.training
        self.eval()

        # Create data-loader to give batches of size 1
        #data_loader = torch.utils.data.DataLoader(dataset, batch_size=policies[0].shape[0], shuffle=False)

        mask = [0.] * self.num_classes
        if allowed_classes is not None:
            for c in allowed_classes:
                mask[c] = 1.0
        else:
            mask = [1.] * self.num_classes

        #for index,(x,y) in tqdm.tqdm(enumerate(data_loader)):
        if True:
            x = inputs
            y = targets
            # run forward pass of model
            if allowed_classes is None:
                output, _ = self.forward(x, policy)
            else:
                out, _ = self.forward(x, policy)
                #output =  out[:, allowed_classes]
                output = out * Variable(torch.tensor(mask)).cuda(non_blocking=True)

            if self.emp_FI:
                # -use provided label to calculate loglikelihood --> "empirical Fisher":
                label = torch.LongTensor([y]) if type(y)==int else y
                if allowed_classes is not None:
                    label = [int(np.where(i == allowed_classes)[0][0]) for i in label.numpy()]
                    label = torch.LongTensor(label)
                label = label.to(self._device())
            else:
                # -use predicted label to calculate loglikelihood:
                label = output.max(1)[1]
            # calculate negative log-likelihood
            negloglikelihood = F.nll_loss(F.log_softmax(output, dim=1), label)

            # Calculate gradient of negative loglikelihood
            self.zero_grad()
            negloglikelihood.backward()

            # Square gradients and keep running sum
            for n, p in self.named_parameters():
                if p.requires_grad: #and arch_is_picked(n, policies):
                    n = n.replace('.', '__')
                    if p.grad is not None:
                        est_fisher_info[n] += p.grad.detach() ** 2

        # Normalize by sample size used for estimation
        est_fisher_info = {n: p for n, p in est_fisher_info.items()}

        # Store new values in the network
        for n, p in self.named_parameters():
            if p.requires_grad:

                n = n.replace('.', '__')
                # -mode (=MAP parameter estimate)
                self.register_buffer('{}_EWC_prev_task{}'.format(n, "" if self.online else self.EWC_task_count+1),
                                     p.detach().clone())
                # -precision (approximated by diagonal Fisher Information matrix)
                if self.online and self.EWC_task_count==1:
                    existing_values = getattr(self, '{}_EWC_estimated_fisher'.format(n))
                    est_fisher_info[n] += self.gamma * existing_values
                self.register_buffer('{}_EWC_estimated_fisher{}'.format(n, "" if self.online else self.EWC_task_count+1),
                                     est_fisher_info[n])

        # If "offline EWC", increase task-count (for "online EWC", set it to 1 to indicate EWC-loss can be calculated)
        self.EWC_task_count = 1 if self.online else self.EWC_task_count + 1

        # Set model back to its initial mode
        self.train(mode=mode)

        del est_fisher_info
        gc.collect()

    def ewc_loss(self, policies):
        '''Calculate EWC-loss.'''
        if self.EWC_task_count>0:
            losses = []

            # If "offline EWC", loop over all previous tasks (if "online EWC", [EWC_task_count]=1 so only 1 iteration)
            for task in range(1, self.EWC_task_count+1):
                #print('---------')
                #print(policies)
                for n, p in self.named_parameters():
                    if p.requires_grad :#and arch_is_picked(n, policies):
                        # Retrieve stored mode (MAP estimate) and precision (Fisher Information matrix)
                        n = n.replace('.', '__')
                        mean = getattr(self, '{}_EWC_prev_task{}'.format(n, "" if self.online else task))
                        fisher = getattr(self, '{}_EWC_estimated_fisher{}'.format(n, "" if self.online else task))
                        # If "online EWC", apply decay-term to the running sum of the Fisher Information matrices
                        fisher = self.gamma*fisher if self.online else fisher
                        # Calculate EWC-loss
                        losses.append((fisher * (p-mean)**2).sum())

            # Sum EWC-loss from all parameters (and from all tasks, if "offline EWC")
            gc.collect()
            return (1./2)*sum(losses)#, losses_dict
        else:
            # EWC-loss is 0 if there are no stored mode and precision yet
            #return torch.tensor(0., device=self._device())
            return Variable(torch.tensor(0.)).cuda()

    #------------------MetaGraph-----------------------
    def forward(self, x, policy, drop_path_prob=0):

        x = F.relu(self.norm1(self.conv1(x)))
        t = 0
        lat = Variable(torch.zeros(x.size(0)), requires_grad=False).cuda().float()

        for out_planes, num_blocks, stride in self.cfg:
            for idx in range(num_blocks):
                action = policy[:, t, :].contiguous()

                # early termination if all actions in the batch are zero
                if action[:, :].data.sum() == 0:
                    if idx != 0:
                        t += 1
                        continue
                    else:
                        feature_map = [self.layers[t].forward(x, 0)]
                        lat_in_this_block = [1. * (action[:, 0]+1)]
                        x = sum(feature_map)
                        lat += sum(lat_in_this_block).float()
                else:
                    prev_x = x
                    action_mask     = [action[:, i].contiguous().float().view(-1, 1, 1, 1) for i in range(action.size(1))]
                    feature_map_raw = [self.layers[t].forward(x, i) for i in range(action.size(1))]
                    feature_map     = [feature_map_raw[i] * action_mask[i] for i in range(action.size(1))]
                    lat_in_this_block = [1. * action[:, i].float() for i in range(action.size(1))]
                    x = sum(feature_map)
                    lat_delta = sum(lat_in_this_block).float()
                    is_no_action = (action.data==0).all(1)
                    if is_no_action.any():
                        if idx!=0:
                            # Exactly use previous feature map
                            x[is_no_action] = prev_x[is_no_action]
                            lat_delta[is_no_action] = 0
                        else:
                            # Feature map reshaped, use default layer
                            x[is_no_action] = feature_map_raw[2][is_no_action]
                            lat_delta[is_no_action] = 1.
                    lat += lat_delta

                t += 1

        x = F.relu(x)
        if self.dset_name == 'Tiny':
            x = F.avg_pool2d(x, 7)
        else:
            x = F.avg_pool2d(x, 4)
        x = x.view(x.size(0), -1)

        x = self.linear(x)

        return x, lat

    def store_back(self):
        for t in range(self.num_of_blocks):
            self.layers[t].store_back()

    def assign_mask(self):
        for t in range(self.num_of_blocks):
            self.layers[t].assign_mask()

    def count_real_params(self):
        num_param = 0
        for t in range(self.num_of_blocks):
            num_param += self.layers[t].count_parameters()

        print(' [*] Total num of real parameter: %.2f M' % (num_param/1e6))

class ResNet18(MetaGraph):
    # out_planes, num_blocks, stride
    cfg = [
           (64, 2, 1),
           (128, 2, 2),
           (256, 2, 2),
           (512, 2, 2)]
    def __init__(self, config=None, num_classes=10, detailed=True, norm_type='GroupNorm'):
        super(ResNet18, self).__init__()
        self.detailed = detailed
        self.dset_name = 'cifar'
        self.num_of_actions = 4
        self.num_of_blocks = sum([num_blocks for out_planes, num_blocks, stride in self.cfg])
        self.num_classes = num_classes
        # stem
        self.conv1 = nn.Conv2d(3, 64, kernel_size=3, stride=1, padding=1, bias=False)
        self.norm_type = norm_type
        if norm_type == 'GroupNorm':
            self.norm1 = nn.GroupNorm(2, 64)
        else:
            self.norm1 = nn.BatchNorm2d(64)

        self.layers = self._make_layers(in_planes=64)
        self.linear = nn.Linear(512, num_classes)

        if self.detailed:
            pytorch_total_params = sum(p.numel() for p in self.parameters())
            '''
            For easier implementation, we duplicate the weight from the superkernel to create different options during forwarding and backwarding.
            One may implemented it by reusing the same weight without creating dummy weights.
            '''
            print(' [*] Total num of parameters with dummy weights: %.2f M' % (pytorch_total_params/1e6))
            self.count_real_params()

    def _make_layers(self, in_planes):
        layers = []
        for out_planes, num_blocks, stride in self.cfg:
            strides = [stride] + [1]*(num_blocks-1)
            for stride in strides:
                block = self._make_action(in_planes, out_planes, stride)
                layers.append(block)
                in_planes = out_planes

        self.num_of_layers = len(layers)
        if self.detailed:
            print(" [*] Total num of layers : ", self.num_of_layers)
        return nn.Sequential(*layers)

    def _make_action(self, inp, oup, stride):
        action = base_single.GroupBasicBlock(inp, oup, stride, 7, self.num_of_actions, self.num_of_blocks, groups=4, norm_type=self.norm_type)
        return action

class ResNet18_64(MetaGraph):
    # out_planes, num_blocks, stride
    cfg = [
           (64, 2, 1),
           (128, 2, 2),
           (256, 2, 2),
           (512, 2, 2)]
    def __init__(self, config=None, num_classes=200, detailed=True, norm_type='GroupNorm'):
        super(ResNet18_64, self).__init__()
        self.num_classes = num_classes
        self.dset_name = 'Tiny'
        self.num_of_actions = 4
        self.num_of_blocks = sum([num_blocks for out_planes, num_blocks, stride in self.cfg])
        self.detailed = detailed
        # stem
        self.conv1 = nn.Conv2d(3, 64, kernel_size=3, stride=1, padding=1, bias=False)
        if norm_type == 'GroupNorm':
            self.norm1 = nn.GroupNorm(2, 64)
        else:
            self.norm1 = nn.BatchNorm2d(64)
        self.layers = self._make_layers(in_planes=64)
        self.linear = nn.Linear(512, num_classes)

        if self.detailed:
            '''
            For easier implementation, we duplicate the weight from the superkernel to create different options during forwarding and backwarding.
            One may implemented it by reusing the same weight without creating dummy weights.
            '''
            pytorch_total_params = sum(p.numel() for p in self.parameters())
            print(' [*] Total num of parameters with dummy weights: %.2f M' % (pytorch_total_params/1e6))
            self.count_real_params()

    def _make_layers(self, in_planes):
        layers = []
        for out_planes, num_blocks, stride in self.cfg:
            strides = [stride] + [1]*(num_blocks-1)
            for stride in strides:
                block = self._make_action(in_planes, out_planes, stride)
                layers.append(block)
                in_planes = out_planes
        self.num_of_layers = len(layers)
        if self.detailed:
            print(" [*] Total num of layers: ", self.num_of_layers)
        return nn.Sequential(*layers)

    def _make_action(self, inp, oup, stride):
        action = base_single.GroupBasicBlock(inp, oup, stride, 7, self.num_of_actions, self.num_of_blocks, groups=4, norm_type=self.norm_type)

        return action

max_stars_repo_name: andyqmongo/InstAParam

max_issues_repo_name: andyqmongo/InstAParam

max_forks_repo_name: andyqmongo/InstAParam