text:
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import scipy.spatial

# Gradient descent optimization
# The learning rate is specified by eta
class GDOptimizer(object):
    def __init__(self, eta):
        self.eta = eta

    def initialize(self, layers):
        pass

    # This function performs one gradient descent step
    # layers is a list of dense layers in the network
    # g is a list of gradients going into each layer before the nonlinear activation
    # a is a list of of the activations of each node in the previous layer going
    def update(self, layers, g, a):
        m = a[0].shape[1]
        for layer, curGrad, curA in zip(layers, g, a):
            layer.updateWeights(-self.eta / m * np.dot(curGrad, curA.T))
            layer.updateBias(-self.eta / m * np.sum(curGrad, 1).reshape(layer.b.shape))

# Cost function used to compute prediction errors
class QuadraticCost(object):

    # Compute the squared error between the prediction yp and the observation y
    # This method should compute the cost per element such that the output is the
    # same shape as y and yp
    @staticmethod
    def fx(y, yp):
        return 0.5 * np.square(yp - y)

    # Derivative of the cost function with respect to yp
    @staticmethod
    def dx(y, yp):
        return y - yp

# Sigmoid function fully implemented as an example
class SigmoidActivation(object):
    @staticmethod
    def fx(z):
        return 1 / (1 + np.exp(-z))

    @staticmethod
    def dx(z):
        return SigmoidActivation.fx(z) * (1 - SigmoidActivation.fx(z))

# Hyperbolic tangent function
class TanhActivation(object):

    # Compute tanh for each element in the input z
    @staticmethod
    def fx(z):
        return np.tanh(z)

    # Compute the derivative of the tanh function with respect to z
    @staticmethod
    def dx(z):
        return 1 - np.square(np.tanh(z))

# Rectified linear unit
class ReLUActivation(object):
    @staticmethod
    def fx(z):
        return np.maximum(0, z)

    @staticmethod
    def dx(z):
        return (z > 0).astype('float')

# Linear activation
class LinearActivation(object):
    @staticmethod
    def fx(z):
        return z

    @staticmethod
    def dx(z):
        return np.ones(z.shape)

# This class represents a single hidden or output layer in the neural network
class DenseLayer(object):

    # numNodes: number of hidden units in the layer
    # activation: the activation function to use in this layer
    def __init__(self, numNodes, activation):
        self.numNodes = numNodes
        self.activation = activation

    def getNumNodes(self):
        return self.numNodes

    # Initialize the weight matrix of this layer based on the size of the matrix W
    def initialize(self, fanIn, scale=1.0):
        s = scale * np.sqrt(6.0 / (self.numNodes + fanIn))
        self.W = np.random.normal(0, s, (self.numNodes, fanIn))
        self.b = np.random.uniform(-1, 1, (self.numNodes, 1))

    # Apply the activation function of the layer on the input z
    def a(self, z):
        return self.activation.fx(z)

    # Compute the linear part of the layer
    # The input a is an n x k matrix where n is the number of samples
    # and k is the dimension of the previous layer (or the input to the network)
    def z(self, a):
        return self.W.dot(a) + self.b  # Note, this is implemented where we assume a is k x n

    # Compute the derivative of the layer's activation function with respect to z
    # where z is the output of the above function.
    # This derivative does not contain the derivative of the matrix multiplication
    # in the layer.  That part is computed below in the model class.
    def dx(self, z):
        return self.activation.dx(z)

    # Update the weights of the layer by adding dW to the weights
    def updateWeights(self, dW):
        self.W = self.W + dW

    # Update the bias of the layer by adding db to the bias
    def updateBias(self, db):
        self.b = self.b + db

# This class handles stacking layers together to form the completed neural network
class Model(object):

    # inputSize: the dimension of the inputs that go into the network
    def __init__(self, inputSize):
        self.layers = []
        self.inputSize = inputSize

    # Add a layer to the end of the network
    def addLayer(self, layer):
        self.layers.append(layer)

    # Get the output size of the layer at the given index
    def getLayerSize(self, index):
        if index >= len(self.layers):
            return self.layers[-1].getNumNodes()
        elif index < 0:
            return self.inputSize
        else:
            return self.layers[index].getNumNodes()

    # Initialize the weights of all of the layers in the network and set the cost
    # function to use for optimization
    def initialize(self, cost, initializeLayers=True):
        self.cost = cost
        if initializeLayers:
            for i in range(0, len(self.layers)):
                if i == len(self.layers) - 1:
                    self.layers[i].initialize(self.getLayerSize(i - 1))
                else:
                    self.layers[i].initialize(self.getLayerSize(i - 1))

    # Compute the output of the network given some input a
    # The matrix a has shape n x k where n is the number of samples and
    # k is the dimension
    # This function returns
    # yp - the output of the network
    # a - a list of inputs for each layer of the newtork where
    #     a[i] is the input to layer i
    # z - a list of values for each layer after evaluating layer.z(a) but
    #     before evaluating the nonlinear function for the layer
    def evaluate(self, x):
        curA = x.T
        a = [curA]
        z = []
        for layer in self.layers:
            z.append(layer.z(curA))
            curA = layer.a(z[-1])
            a.append(curA)
        yp = a.pop()
        return yp, a, z

    # Compute the output of the network given some input a
    # The matrix a has shape n x k where n is the number of samples and
    # k is the dimension
    def predict(self, a):
        a, _, _ = self.evaluate(a)
        return a.T

    # Train the network given the inputs x and the corresponding observations y
    # The network should be trained for numEpochs iterations using the supplied
    # optimizer
    def train(self, x, y, numEpochs, optimizer):

        # Initialize some stuff
        hist = []
        optimizer.initialize(self.layers)

        # Run for the specified number of epochs
        for epoch in range(0, numEpochs):

            # Feed forward
            # Save the output of each layer in the list a
            # After the network has been evaluated, a should contain the
            # input x and the output of each layer except for the last layer
            yp, a, z = self.evaluate(x)

            # Compute the error
            C = self.cost.fx(yp, y.T)
            d = self.cost.dx(yp, y.T)
            grad = []

            # Backpropogate the error
            idx = len(self.layers)
            for layer, curZ in zip(reversed(self.layers), reversed(z)):
                idx = idx - 1
                # Here, we compute dMSE/dz_i because in the update
                # function for the optimizer, we do not give it
                # the z values we compute from evaluating the network
                grad.insert(0, np.multiply(d, layer.dx(curZ)))
                d = np.dot(layer.W.T, grad[0])

            # Update the errors
            optimizer.update(self.layers, grad, a)

            # Compute the error at the end of the epoch
            yh = self.predict(x)
            C = self.cost.fx(yh, y)
            C = np.mean(C)
            hist.append(C)
        return hist

    def trainBatch(self, x, y, batchSize, numEpochs, optimizer):

        # Copy the data so that we don't affect the original one when shuffling
        x = x.copy()
        y = y.copy()
        hist = []
        n = x.shape[0]

        for epoch in np.arange(0, numEpochs):

            # Shuffle the data
            r = np.arange(0, x.shape[0])
            x = x[r, :]
            y = y[r, :]
            e = []

            # Split the data in chunks and run SGD
            for i in range(0, n, batchSize):
                end = min(i + batchSize, n)
                batchX = x[i:end, :]
                batchY = y[i:end, :]
                e += self.train(batchX, batchY, 1, optimizer)
            hist.append(np.mean(e))

        return hist

########################################################################
######### Part b #######################################################
########################################################################

########################################################################
#########  Gradient Computing and MLE ##################################
########################################################################
def compute_gradient_of_likelihood(single_obj_loc, sensor_loc, single_distance, noise=1):
    """
    Compute the gradient of the loglikelihood function for part a.

    Input:
    single_obj_loc: 1 * d numpy array.
    Location of the single object.

    sensor_loc: k * d numpy array.
    Location of sensor.

    single_distance: k dimensional numpy array.
    Observed distance of the object.

    Output:
    grad: d-dimensional numpy array.

    """
    loc_difference = single_obj_loc - sensor_loc  # k * d.
    phi = np.linalg.norm(loc_difference, axis=1)  # k.
    weight = (phi - single_distance) / phi  # k.

    grad = -np.sum(np.expand_dims(weight, 1) * loc_difference, axis=0) / noise**2  # d
    return grad

########################################################################
######### Part c #################################################
########################################################################
def log_likelihood(obj_loc, sensor_loc, distance, noise=1):
    """
    This function computes the log likelihood (as expressed in Part a).
    Input:
    obj_loc: shape [1,2]
    sensor_loc: shape [7,2]
    distance: shape [7]
    Output:
    The log likelihood function value.
    """
    diff_distance = np.sqrt(np.sum((sensor_loc - obj_loc)**2, axis=1)) - distance
    func_value = -sum((diff_distance)**2) / (2 * noise**2)
    return func_value

########################################################################
######### Part e, f, g #################################################
########################################################################

########################################################################
#########  Gradient Computing and MLE ##################################
########################################################################
def compute_grad_likelihood_part_e(sensor_loc, obj_loc, distance, noise=1):
    """
    Compute the gradient of the loglikelihood function for part d.

    Input:
    sensor_loc: k * d numpy array.
    Location of sensors.

    obj_loc: n * d numpy array.
    Location of the objects.

    distance: n * k dimensional numpy array.
    Observed distance of the object.

    Output:
    grad: k * d numpy array.
    """
    grad = np.zeros(sensor_loc.shape)
    for i, single_sensor_loc in enumerate(sensor_loc):
        single_distance = distance[:, i]
        grad[i] = compute_gradient_of_likelihood(single_sensor_loc, obj_loc, single_distance,
                                                 noise)

    return grad

def find_mle_by_grad_descent_part_e(initial_sensor_loc,
                                    obj_loc,
                                    distance,
                                    noise=1,
                                    lr=0.001,
                                    num_iters=1000):
    """
    Compute the gradient of the loglikelihood function for part a.

    Input:
    initial_sensor_loc: k * d numpy array.
    Initialized Location of the sensors.

    obj_loc: n * d numpy array. Location of the n objects.

    distance: n * k dimensional numpy array.
    Observed distance of the n object.

    Output:
    sensor_loc: k * d numpy array. The mle for the location of the object.

    """
    sensor_loc = initial_sensor_loc
    for t in range(num_iters):
        sensor_loc += lr * compute_grad_likelihood_part_e(sensor_loc, obj_loc, distance, noise)

    return sensor_loc

########################################################################

#########  Estimate distance given estimated sensor locations. #########
########################################################################

def compute_distance_with_sensor_and_obj_loc(sensor_loc, obj_loc):
    """
    stimate distance given estimated sensor locations.

    Input:
    sensor_loc: k * d numpy array.
    Location of the sensors.

    obj_loc: n * d numpy array. Location of the n objects.

    Output:
    distance: n * k dimensional numpy array.
    """
    estimated_distance = scipy.spatial.distance.cdist(obj_loc, sensor_loc, metric='euclidean')
    return estimated_distance

########################################################################
#########  Data Generating Functions ###################################
########################################################################

def generate_sensors(num_sensors=7, spatial_dim=2):
    """
    Generate sensor locations.
    Input:
    num_sensors: The number of sensors.
    spatial_dim: The spatial dimension.
    Output:
    sensor_loc: num_sensors * spatial_dim numpy array.
    """
    sensor_loc = 100 * np.random.randn(num_sensors, spatial_dim)
    return sensor_loc

def generate_dataset(sensor_loc,
                     num_sensors=7,
                     spatial_dim=2,
                     num_data=1,
                     original_dist=True,
                     noise=1):
    """
    Generate the locations of n points.

    Input:
    sensor_loc: num_sensors * spatial_dim numpy array. Location of sensor.
    num_sensors: The number of sensors.
    spatial_dim: The spatial dimension.
    num_data: The number of points.
    original_dist: Whether the data are generated from the original
    distribution.

    Output:
    obj_loc: num_data * spatial_dim numpy array. The location of the num_data objects.
    distance: num_data * num_sensors numpy array. The distance between object and
    the num_sensors sensors.
    """
    assert num_sensors, spatial_dim == sensor_loc.shape

    obj_loc = 100 * np.random.randn(num_data, spatial_dim)
    if not original_dist:
        obj_loc += 1000

    distance = scipy.spatial.distance.cdist(obj_loc, sensor_loc, metric='euclidean')
    distance += np.random.randn(num_data, num_sensors) * noise
    return distance, obj_loc

max_issues_repo_name: QinganZhao/Machine-Learning

max_forks_repo_name: QinganZhao/Machine-Learning


text:
"""
  This software is governed by the CeCILL-B license under French law and
  abiding by the rules of distribution of free software.  You can  use,
  modify and/ or redistribute the software under the terms of the CeCILL-B
  license as circulated by CEA, CNRS and INRIA at the following URL
  "http://www.cecill.info".
  As a counterpart to the access to the source code and  rights to copy,
  modify and redistribute granted by the license, users are provided only
  with a limited warranty  and the software's author,  the holder of the
  economic rights,  and the successive licensors  have only  limited
  liability.
  In this respect, the user's attention is drawn to the risks associated
  with loading,  using,  modifying and/or developing or reproducing the
  software by the user in light of its specific status of free software,
  that may mean  that it is complicated to manipulate,  and  that  also
  therefore means  that it is reserved for developers  and  experienced
  professionals having in-depth computer knowledge. Users are therefore
  encouraged to load and test the software's suitability as regards their
  requirements in conditions enabling the security of their systems and/or
  data to be ensured and,  more generally, to use and operate it in the
  same conditions as regards security.
  The fact that you are presently reading this means that you have had
  knowledge of the CeCILL-B license and that you accept its terms.
"""

import numpy as np
from functools import partial

from keras.models import Model
from keras.layers import Input, LeakyReLU, Reshape
from keras.layers import Conv3D, Add, UpSampling3D, Activation
from keras.optimizers import Adam
from keras.initializers import lecun_normal
gen_initializer = lecun_normal()

import sys
sys.path.insert(0, './utils')
from layers import wasserstein_loss, ReflectPadding3D, gradient_penalty_loss, InstanceNormalization3D, activation_SegSRGAN, charbonnier_loss
from Adam_lr_mult import LR_Adam

def resnet_blocks(input_res, kernel, name):

    in_res_1 = ReflectPadding3D(padding=1)(input_res)
    out_res_1 = Conv3D(kernel, 3, strides=1, kernel_initializer=gen_initializer,
                       use_bias=False,
                       name=name+'_conv_a',
                       data_format='channels_first')(in_res_1)
    out_res_1 = InstanceNormalization3D(name=name+'_isnorm_a')(out_res_1)
    out_res_1 = Activation('relu')(out_res_1)

    in_res_2 = ReflectPadding3D(padding=1)(out_res_1)
    out_res_2 = Conv3D(kernel, 3, strides=1, kernel_initializer=gen_initializer,
                       use_bias=False,
                       name=name+'_conv_b',
                       data_format='channels_first')(in_res_2)
    out_res_2 = InstanceNormalization3D(name=name+'_isnorm_b')(out_res_2)

    out_res = Add()([out_res_2,input_res])
    return out_res

class SegSRGAN(object):
    def __init__(self, ImageRow=64, ImageColumn=64, ImageDepth=64,
                 FirstDiscriminatorKernel = 32, FirstGeneratorKernel = 16,
                 lamb_rec = 1, lamb_adv = 0.001, lamb_gp = 10,
                 lr_DisModel = 0.0001, lr_GenModel = 0.0001):
        self.ImageRow = ImageRow
        self.ImageColumn = ImageColumn
        self.ImageDepth = ImageDepth
        self.D = None                   # discriminator
        self.G = None                   # generator
        self.DisModel = None            # discriminator model
        self.GenModel = None            # generator model
        self.DiscriminatorKernel = FirstDiscriminatorKernel
        self.GeneratorKernel = FirstGeneratorKernel
        self.lamb_adv = lamb_adv
        self.lamb_rec = lamb_rec
        self.lamb_gp = lamb_gp
        self.lr_DisModel = lr_DisModel
        self.lr_GenModel = lr_GenModel

    def discriminator_block(self, name):
        """Creates a discriminator model that takes an image as input and outputs a single value, representing whether
        the input is real or generated. Unlike normal GANs, the output is not sigmoid and does not represent a probability!
        Instead, the output should be as large and negative as possible for generated inputs and as large and positive
        as possible for real inputs.
        Note that the improved WGAN paper suggests that BatchNormalization should not be used in the discriminator."""

        # In:
        inputs = Input(shape=(2, self.ImageRow, self.ImageColumn, self.ImageDepth), name='dis_input')

        # Input 64
        disnet = Conv3D(self.DiscriminatorKernel*1, 4, strides=2,
                        padding = 'same',
                        kernel_initializer='he_normal',
                        data_format='channels_first',
                        name=name+'_conv_dis_1')(inputs)
        disnet = LeakyReLU(0.01)(disnet)

        # Hidden 1 : 32
        disnet = Conv3D(self.DiscriminatorKernel*2, 4, strides=2,
                        padding = 'same',
                        kernel_initializer='he_normal',
                        data_format='channels_first',
                        name=name+'_conv_dis_2')(disnet)
        disnet = LeakyReLU(0.01)(disnet)

        # Hidden 2 : 16
        disnet = Conv3D(self.DiscriminatorKernel*4, 4, strides=2,
                        padding = 'same',
                        kernel_initializer='he_normal',
                        data_format='channels_first',
                        name=name+'_conv_dis_3')(disnet)
        disnet = LeakyReLU(0.01)(disnet)

        # Hidden 3 : 8
        disnet = Conv3D(self.DiscriminatorKernel*8, 4, strides=2,
                        padding = 'same',
                        kernel_initializer='he_normal',
                        data_format='channels_first',
                        name=name+'_conv_dis_4')(disnet)
        disnet = LeakyReLU(0.01)(disnet)

        # Hidden 4 : 4
        disnet = Conv3D(self.DiscriminatorKernel*16, 4, strides=2,
                        padding = 'same',
                        kernel_initializer='he_normal',
                        data_format='channels_first',
                        name=name+'_conv_dis_5')(disnet)
        disnet = LeakyReLU(0.01)(disnet)

        # Decision : 2
        decision = Conv3D(1, 2, strides=1,
                          use_bias=False,
                          kernel_initializer='he_normal',
                          data_format='channels_first',
                          name='dis_decision')(disnet)
        decision = Reshape((1,))(decision)

        model = Model(inputs=[inputs],outputs=[decision],name=name)
        return model

    def generator_block(self,name):
        #
        inputs = Input(shape=(1, self.ImageRow, self.ImageColumn, self.ImageDepth))

        # Representation
        gennet = ReflectPadding3D(padding=3)(inputs)
        gennet = Conv3D(self.GeneratorKernel, 7, strides=1, kernel_initializer=gen_initializer,
                        use_bias=False,
                        name=name+'_gen_conv1',
                        data_format='channels_first')(gennet)
        gennet = InstanceNormalization3D(name=name+'_gen_isnorm_conv1')(gennet)
        gennet = Activation('relu')(gennet)

        # Downsampling 1
        gennet = ReflectPadding3D(padding=1)(gennet)
        gennet = Conv3D(self.GeneratorKernel*2, 3, strides=2, kernel_initializer=gen_initializer,
                        use_bias=False,
                        name=name+'_gen_conv2',
                        data_format='channels_first')(gennet)
        gennet = InstanceNormalization3D(name=name+'_gen_isnorm_conv2')(gennet)
        gennet = Activation('relu')(gennet)

        # Downsampling 2
        gennet = ReflectPadding3D(padding=1)(gennet)
        gennet = Conv3D(self.GeneratorKernel*4, 3, strides=2, kernel_initializer=gen_initializer,
                        use_bias=False,
                        name=name+'_gen_conv3',
                        data_format='channels_first')(gennet)
        gennet = InstanceNormalization3D(name=name+'_gen_isnorm_conv3')(gennet)
        gennet = Activation('relu')(gennet)

        # Resnet blocks : 6, 8*4 = 32
        gennet = resnet_blocks(gennet, self.GeneratorKernel*4, name=name+'_gen_block1')
        gennet = resnet_blocks(gennet, self.GeneratorKernel*4, name=name+'_gen_block2')
        gennet = resnet_blocks(gennet, self.GeneratorKernel*4, name=name+'_gen_block3')
        gennet = resnet_blocks(gennet, self.GeneratorKernel*4, name=name+'_gen_block4')
        gennet = resnet_blocks(gennet, self.GeneratorKernel*4, name=name+'_gen_block5')
        gennet = resnet_blocks(gennet, self.GeneratorKernel*4, name=name+'_gen_block6')

        # Upsampling 1
        gennet = UpSampling3D(size=(2, 2, 2),
                              data_format='channels_first')(gennet)
        gennet = ReflectPadding3D(padding=1)(gennet)
        gennet = Conv3D(self.GeneratorKernel*2, 3, strides=1, kernel_initializer=gen_initializer,
                        use_bias=False,
                        name=name+'_gen_deconv1',
                        data_format='channels_first')(gennet)
        gennet = InstanceNormalization3D(name=name+'_gen_isnorm_deconv1')(gennet)
        gennet = Activation('relu')(gennet)

        # Upsampling 2
        gennet = UpSampling3D(size=(2, 2, 2),
                              data_format='channels_first')(gennet)
        gennet = ReflectPadding3D(padding=1)(gennet)
        gennet = Conv3D(self.GeneratorKernel, 3, strides=1, kernel_initializer=gen_initializer,
                        use_bias=False,
                        name=name+'_gen_deconv2',
                        data_format='channels_first')(gennet)
        gennet = InstanceNormalization3D(name=name+'_gen_isnorm_deconv2')(gennet)
        gennet = Activation('relu')(gennet)

        # Reconstruction
        gennet = ReflectPadding3D(padding=3)(gennet)
        gennet = Conv3D(2, 7, strides=1, kernel_initializer=gen_initializer,
                        use_bias=False,
                        name=name+'_gen_1conv',
                        data_format='channels_first')(gennet)

        predictions = gennet
        predictions = activation_SegSRGAN()([predictions,inputs])

        model = Model(inputs=inputs,outputs=predictions,name=name)
        return model

    def generator(self):
        if self.G:
            return self.G

        self.G = self.generator_block('G')
        return self.G

    def discriminator(self):
        if self.D:
            return self.D

        self.D = self.discriminator_block('DX')
        return self.D

    def generator_model(self):
        if self.GenModel:
            return self.GenModel

        print "We freeze the weights of Discriminator by setting their learning rate as 0 when updating Generator !"
        # We freeze the weights of Discriminator by setting their learning rate as 0 when updating Generator !
        AllParameters = 63
        GeneratorParameters = 52
        multipliers = np.ones(AllParameters)
        for idx in range(GeneratorParameters,AllParameters):
            multipliers[idx]= 0.0

        # Input
        input_gen = Input(shape=(1, self.ImageRow, self.ImageColumn, self.ImageDepth), name='input_gen')

        # Network
        Gx_gen = self.generator()(input_gen)               # Fake X
        fool_decision = self.discriminator()(Gx_gen)      # Fooling D

        # Model
        self.GenModel = Model(input_gen,[fool_decision,Gx_gen])
        self.GenModel.compile(LR_Adam(lr=self.lr_GenModel, beta_1=0.5, beta_2=0.999, multipliers=multipliers),
                            loss=[wasserstein_loss, charbonnier_loss],
                            loss_weights=[self.lamb_adv, self.lamb_rec])

        return self.GenModel

    def discriminator_model(self):
        if self.DisModel:
            return self.DisModel

        # Input
        real_dis = Input(shape=(2, self.ImageRow, self.ImageColumn, self.ImageDepth), name='real_dis')
        fake_dis = Input(shape=(2, self.ImageRow, self.ImageColumn, self.ImageDepth), name='fake_dis')
        interp_dis = Input(shape=(2, self.ImageRow, self.ImageColumn, self.ImageDepth), name='interp_dis')

        # Discriminator
        real_decision = self.discriminator()(real_dis)            # Real X
        fake_decision = self.discriminator()(fake_dis)            # Fake X
        interp_decision = self.discriminator()(interp_dis)        # interpolation X

        # GP loss
        partial_gp_loss = partial(gradient_penalty_loss,
                          averaged_samples=interp_dis,
                          gradient_penalty_weight=self.lamb_gp)
        partial_gp_loss.__name__ = 'gradient_penalty'  # Functions need names or Keras will throw an error*

        # Model
        self.DisModel = Model([real_dis,fake_dis,interp_dis], [real_decision, fake_decision,interp_decision])
        self.DisModel.compile(Adam(lr=self.lr_DisModel, beta_1=0.5, beta_2=0.999),
                            loss=[wasserstein_loss,wasserstein_loss,partial_gp_loss],
                            loss_weights=[1,1,self.lamb_gp])

        return self.DisModel

ext: py

max_stars_repo_name: rousseau/deepBrain

max_issues_repo_name: rousseau/deepBrain

max_forks_repo_name: rousseau/deepBrain

reason: import numpy


text:
"""Depth Measures Module.

This module includes different methods to order functional data,
from the center (larger values) outwards(smaller ones)."""

from functools import reduce
import math

import scipy.integrate
from scipy.stats import rankdata

import numpy as np

__author__ = "Amanda Hernando Bernabé"
__email__ = "amanda.hernando@estudiante.uam.es"

def outlyingness_to_depth(outlyingness, *, supreme=None):
    r"""Convert outlyingness function to depth function.

    An outlyingness function :math:`O(x)` can be converted to a depth
    function as

    .. math::
        D(x) = \frac{1}{1 + O(x)}

    if :math:`O(x)` is unbounded or as

    .. math::
        D(x) = 1 - \frac{O(x)}{\sup O(x)}

    if :math:`O(x)` is bounded ([Se06]_).

    Args:
        outlyingness (Callable): Outlyingness function.
        supreme (float, optional): Supreme value of the outlyingness function.

    Returns:
        Callable: The corresponding depth function.

    References:
        .. [Se06] Serfling, R. (2006). Depth functions in nonparametric
           multivariate inference. DIMACS Series in Discrete Mathematics and
           Theoretical Computer Science, 72, 1.
    """

    if supreme is None or math.isinf(supreme):
        def depth(*args, **kwargs):
            return 1 / (1 + outlyingness(*args, **kwargs))
    else:
        def depth(*args, **kwargs):
            return 1 - outlyingness(*args, **kwargs) / supreme

    return depth

def _rank_samples(fdatagrid):
    """Ranks the he samples in the FDataGrid at each point of discretisation.

    Args:
        fdatagrid (FDataGrid): Object whose samples are ranked.

    Returns:
        numpy.darray: Array containing the ranks of the sample points.

    Examples:
        Univariate setting:

        >>> import skfda
        >>>
        >>> data_matrix = [[1, 1, 2, 3, 2.5, 2],
        ...                [0.5, 0.5, 1, 2, 1.5, 1],
        ...                [-1, -1, -0.5, 1, 1, 0.5],
        ...                [-0.5, -0.5, -0.5, -1, -1, -1]]
        >>> sample_points = [0, 2, 4, 6, 8, 10]
        >>> fd = skfda.FDataGrid(data_matrix, sample_points)
        >>> _rank_samples(fd)
        array([[ 4.,  4.,  4.,  4.,  4.,  4.],
               [ 3.,  3.,  3.,  3.,  3.,  3.],
               [ 1.,  1.,  2.,  2.,  2.,  2.],
               [ 2.,  2.,  2.,  1.,  1.,  1.]])

        Several input dimensions:

        >>> data_matrix = [[[[1], [0.7], [1]],
        ...                 [[4], [0.4], [5]]],
        ...                [[[2], [0.5], [2]],
        ...                 [[3], [0.6], [3]]]]
        >>> sample_points = [[2, 4], [3, 6, 8]]
        >>> fd = skfda.FDataGrid(data_matrix, sample_points)
        >>> _rank_samples(fd)
        array([[[ 1.,  2.,  1.],
                [ 2.,  1.,  2.]],
               [[ 2.,  1.,  2.],
                [ 1.,  2.,  1.]]])

    """
    if fdatagrid.dim_codomain > 1:
        raise ValueError("Currently multivariate data is not allowed")

    ranks = np.zeros(fdatagrid.data_matrix.shape[:-1])

    for index, _ in np.ndenumerate(ranks[0]):
        ranks[(slice(None),) + index] = rankdata(
            fdatagrid.data_matrix[(slice(None),) + index + (0,)], method='max')
    return ranks

def band_depth(fdatagrid, *, pointwise=False):
    """Implementation of Band Depth for functional data.

    The band depth of each sample is obtained by computing the fraction of the
    bands determined by two sample curves containing the whole graph of the
    first one. In the case the fdatagrid domain dimension is 2, instead of
    curves, surfaces determine the bands. In larger dimensions, the hyperplanes
    determine the bands.

    Args:
        fdatagrid (FDataGrid): Object over whose samples the band depth is
            going to be calculated.
        pointwise (boolean, optional): Indicates if the pointwise depth is
            returned instead. Defaults to False.

    Returns:
        depth (numpy.darray): Array containing the band depth of the samples,
            or the band depth of the samples at each point of discretization
            if pointwise equals to True.

    Examples:

        >>> import skfda
        >>>
        >>> data_matrix = [[1, 1, 2, 3, 2.5, 2],
        ...                [0.5, 0.5, 1, 2, 1.5, 1],
        ...                [-1, -1, -0.5, 1, 1, 0.5],
        ...                [-0.5, -0.5, -0.5, -1, -1, -1]]
        >>> sample_points = [0, 2, 4, 6, 8, 10]
        >>> fd = skfda.FDataGrid(data_matrix, sample_points)
        >>> band_depth(fd)
        array([ 0.5       ,  0.83333333,  0.5       ,  0.5       ])

    """
    if pointwise:
        return modified_band_depth(fdatagrid, pointwise)
    else:
        n = fdatagrid.n_samples
        nchoose2 = n * (n - 1) / 2

        ranks = _rank_samples(fdatagrid)
        axis = tuple(range(1, fdatagrid.dim_domain + 1))
        n_samples_above = fdatagrid.n_samples - np.amax(ranks, axis=axis)
        n_samples_below = np.amin(ranks, axis=axis) - 1
        depth = ((n_samples_below * n_samples_above + fdatagrid.n_samples - 1)
                 / nchoose2)

        return depth

def modified_band_depth(fdatagrid, *, pointwise=False):
    """Implementation of Modified Band Depth for functional data.

    The band depth of each sample is obtained by computing the fraction of time
    its graph is contained in the bands determined by two sample curves.
    In the case the fdatagrid domain dimension is 2, instead of curves,
    surfaces determine the bands. In larger dimensions, the hyperplanes
    determine the bands.

    Args:
        fdatagrid (FDataGrid): Object over whose samples the modified band
            depth is going to be calculated.
        pointwise (boolean, optional): Indicates if the pointwise depth is
            returned instead. Defaults to False.

    Returns:
        depth (numpy.darray): Array containing the modified band depth of the
            samples, or the modified band depth of the samples at each point
            of discretization if pointwise equals to True.

    Examples:

        >>> import skfda
        >>>
        >>> data_matrix = [[1, 1, 2, 3, 2.5, 2],
        ...                [0.5, 0.5, 1, 2, 1.5, 1],
        ...                [-1, -1, -0.5, 1, 1, 0.5],
        ...                [-0.5, -0.5, -0.5, -1, -1, -1]]
        >>> sample_points = [0, 2, 4, 6, 8, 10]
        >>> fd = skfda.FDataGrid(data_matrix, sample_points)
        >>> depth = modified_band_depth(fd)
        >>> depth.round(2)
        array([ 0.5 ,  0.83,  0.72,  0.67])
        >>> pointwise = modified_band_depth(fd, pointwise = True)
        >>> pointwise.round(2)
        array([[ 0.5 ,  0.5 ,  0.5 ,  0.5 ,  0.5 ,  0.5 ],
               [ 0.83,  0.83,  0.83,  0.83,  0.83,  0.83],
               [ 0.5 ,  0.5 ,  0.83,  0.83,  0.83,  0.83],
               [ 0.83,  0.83,  0.83,  0.5 ,  0.5 ,  0.5 ]])

    """
    n = fdatagrid.n_samples
    nchoose2 = n * (n - 1) / 2

    ranks = _rank_samples(fdatagrid)
    n_samples_above = fdatagrid.n_samples - ranks
    n_samples_below = ranks - 1
    match = n_samples_above * n_samples_below
    axis = tuple(range(1, fdatagrid.dim_domain + 1))

    if pointwise:
        depth_pointwise = (match + fdatagrid.n_samples - 1) / nchoose2

        return depth_pointwise
    else:
        npoints_sample = reduce(lambda x, y: x * len(y),
                                fdatagrid.sample_points, 1)
        proportion = match.sum(axis=axis) / npoints_sample
        depth = (proportion + fdatagrid.n_samples - 1) / nchoose2

        return depth

def _cumulative_distribution(column):
    """Calculates the cumulative distribution function of the values passed to
    the function and evaluates it at each point.

    Args:
        column (numpy.darray): Array containing the values over which the
            distribution function is calculated.

    Returns:
        numpy.darray: Array containing the evaluation at each point of the
            distribution function.

    Examples:
        >>> _cumulative_distribution(np.array([1, 4, 5, 1, 2, 2, 4, 1, 1, 3]))
        array([ 0.4,  0.9,  1. ,  0.4,  0.6,  0.6,  0.9,  0.4,  0.4,  0.7])

    """
    if len(column.shape) != 1:
        raise ValueError("Only supported 1 dimensional arrays.")
    _, indexes, counts = np.unique(column, return_inverse=True,
                                   return_counts=True)
    count_cumulative = np.cumsum(counts) / len(column)
    return count_cumulative[indexes].reshape(column.shape)

def fraiman_muniz_depth(fdatagrid, *, pointwise=False):
    r"""Implementation of Fraiman and Muniz (FM) Depth for functional data.

    Each column is considered as the samples of an aleatory variable.
    The univariate depth of each of the samples of each column is calculated
    as follows:

    .. math::
        D(x) = 1 - \left\lvert \frac{1}{2}- F(x)\right\rvert

    Where :math:`F` stands for the marginal univariate distribution function of
    each column.

    The depth of a sample is the result of integrating the previously computed
    depth for each of its points and normalizing dividing by the length of
    the interval.

    Args:
        fdatagrid (FDataGrid): Object over whose samples the FM depth is going
            to be calculated.
        pointwise (boolean, optional): Indicates if the pointwise depth is
             returned instead. Defaults to False.

    Returns:
        depth (numpy.darray): Array containing the Fraiman-Muniz depth of the
            samples, or the Fraiman-Muniz of the samples at each point
            of discretization if pointwise equals to True.

    Examples:
        Currently, this depth function can only be used
        for univariate functional data:

        >>> import skfda
        >>>
        >>> data_matrix = [[1, 1, 2, 3, 2.5, 2],
        ...                [0.5, 0.5, 1, 2, 1.5, 1],
        ...                [-1, -1, -0.5, 1, 1, 0.5],
        ...                [-0.5, -0.5, -0.5, -1, -1, -1]]
        >>> sample_points = [0, 2, 4, 6, 8, 10]
        >>> fd = skfda.FDataGrid(data_matrix, sample_points)
        >>> fraiman_muniz_depth(fd)
        array([ 0.5  ,  0.75 ,  0.925,  0.875])

        You can use ``pointwise`` to obtain the pointwise depth,
        before the integral is applied.

        >>> pointwise = fraiman_muniz_depth(fd, pointwise = True)
        >>> pointwise
        array([[ 0.5 ,  0.5 ,  0.5 ,  0.5 ,  0.5 ,  0.5 ],
               [ 0.75,  0.75,  0.75,  0.75,  0.75,  0.75],
               [ 0.75,  0.75,  1.  ,  1.  ,  1.  ,  1.  ],
               [ 1.  ,  1.  ,  1.  ,  0.75,  0.75,  0.75]])

    """
    if fdatagrid.dim_domain > 1 or fdatagrid.dim_codomain > 1:
        raise ValueError("Currently multivariate data is not allowed")

    pointwise_depth = np.array([
        1 - abs(0.5 - _cumulative_distribution(
            fdatagrid.data_matrix[:, i, 0])
        ) for i in range(len(fdatagrid.sample_points[0]))]).T

    if pointwise:
        return pointwise_depth
    else:

        interval_len = (fdatagrid.domain_range[0][1]
                        - fdatagrid.domain_range[0][0])

        depth = (scipy.integrate.simps(pointwise_depth,
                                       fdatagrid.sample_points[0])
                 / interval_len)

        return depth

max_issues_repo_name: KonstantinKlepikov/scikit-fda

max_forks_repo_name: KonstantinKlepikov/scikit-fda

reason: import numpy,import scipy,from scipy


text:
from dolfin import *
from FIAT import *
import numpy

def tangential_edge_integral(mesh, f, quad_degree):
    """
    This function uses quadrature to compute the average of
    of the dot product  f.t  on each edge of the given mesh.
    Here f is a given vector valued function and t is the vector
    pointing from the first to the second vertex of the edge.
    The returned value is a numpy array edge_int where
        edge_int[edge.index()] = avg of f.t on that edge
    """
    # Reference interval
    L = reference_element.DefaultLine()
    # Coordinates for the reference interval
    x = numpy.array(L.get_vertices())
    # Make quadrature rule on the line
    Q = make_quadrature(L, quad_degree)
    # Quadrature points
    qp = Q.get_points()
    # Quadrature weights
    qw = Q.get_weights()
    # Get mesh vertex coordinates
    coor = mesh.coordinates()
    # Create an array to save the result
    edge_int = numpy.zeros(mesh.num_edges())
    # Integrate over all the edges
    for e in edges(mesh):
        # Extract end points of the edge
        end_pts = e.entities(0)
        y = numpy.vstack([coor[end_pts[0]], coor[end_pts[1]]])
        # Compute the tangent.  Following the FEniCS convention this
        # is the vector pointing from the first to the end point
        # of the edge.
        t = y[1] - y[0]
        # Make the affine map from the reference line
        (A, b) = reference_element.make_affine_mapping(x, y)
        # Quadrature points on the physical edge
        phys_qp = map(lambda x: A.dot(x)+b, qp)
        # Evaluate f.t at quadrature points
        vals = map(lambda x: t.dot(f(x)), phys_qp)
        # Approximate average of f.t on the edge using the quadrature rule
        edge_int[e.index()] = qw.dot(vals)/2.0
    return edge_int

def build_edge2dof_map(V):
    """
    This function takes a N1Curl(1) space and return an integer valued array edge2dof.
    This array has the number of edges as its length. In particular
        edge2dof[i] = j
    means that dof #i, that is u.vector()[i], is associated to edge #j.
    """
    # Extract the cell to edge map (given an cell index, it returns the indices of its edges)
    cell2edges = V.mesh().topology()(3, 1)
    # Extract the cell dofmap (given a cell index, it returns the dof numbers)
    cell2dofs = V.dofmap().cell_dofs
    # Array to save the result
    edge2dof = numpy.zeros(mesh.num_edges(), dtype="int")
    # Iterate over cells, associating the edges to the dofs for that cell
    for c in range(mesh.num_cells()):
        # get the global edge numbers for this cell
        c_edges = cell2edges(c)
        # get the global dof numbers for this cell
        c_dofs = cell2dofs(c)
        # associate the edge numbers to the corresponding dof numbers
        edge2dof[c_dofs] = c_edges
    # This algorithm might not look fast as it does quite some redundant work. In actual
    # runs, for most meshes, this is not the most time consuming step and does not take
    # more than a milisecond.
    return edge2dof

def n1curl_1_canonical_projection(f, mesh, quad_degree):
    # Initialize the mesh
    mesh.init()
    # Compute the average edge integrals
    edge_int = tangential_edge_integral(mesh, f, quad_degree)
    # Create the return value
    V = FunctionSpace(mesh, "N1curl", 1)
    u = Function(V)
    # Permute the average edge integrals to match the order of the dofs.
    u.vector()[:] = edge_int.take(build_edge2dof_map(V))
    return u

f = Expression(("1.0 - x[1] + 2*x[2]", "3.0 + x[0] + 3*x[2]", "2.0 - 2*x[0] - 3*x[1]"))
domain = Box(0., 0., 0., 1., 1., 1.)
n = 8; mesh = BoxMesh(0., 0., 0., 1., 1., 1., n,n,n)
quad_degree = 2
V = FunctionSpace(mesh, "N1curl", 1)
u = n1curl_1_canonical_projection(f, mesh, quad_degree)
v = interpolate(f, V)

print(assemble(dot(u-v,u-v)*dx))

max_stars_repo_name: wathen/PhD

max_issues_repo_name: wathen/PhD

max_forks_repo_name: wathen/PhD


text:
from mu import *
from scipy import optimize,misc
import time
import sys

r_pixel = [0.04,0.07,0.11]
r_tib=[0.26,0.32,0.43,0.62]
r_tob=[0.71,0.79,0.88,0.97,1.07,1.14]
r_tracker=np.append(np.append(r_pixel,r_tib), r_tob)

r_mb1= [4.20,4.22,4.24,4.26,4.42,4.44,4.46,4.48]
r_mb2= [5.00,5.02,5.04,5.06,5.24,5.26,5.28,5.30]
r_mb3= [6.08,6.10,6.12,6.14,6.32,6.34,6.36,6.38]
r_mb4= [7.10,7.12,7.14,7.16,7.34,7.36,7.38,7.40]

r_drift=np.append(np.append(r_mb1,r_mb2),np.append(r_mb3,r_mb4))

# Errors

err_1_t=10.*10**(-6) # Pixel detector transverse error, units: meters
err_1_z=10.*10**(-6) # Pixel detector z error, units: meters
err_1 = [err_1_t, err_1_z]

err_2_t=20.*10**(-6) # Internal Silicon Strip Tracker transverse error, units: meters
err_2_z=20.*10**(-6) # Internal Silicon Strip Tracker z error, units: meters
err_2 = [err_2_t, err_2_z]

err_3_t=30.*10**(-6) # External Silicon Strip Tracker transverse error, units: meters
err_3_z=30.*10**(-6) # External Silicon Strip Tracker z error, units: meters
err_3 = [err_3_t, err_3_z]

# Incorrect values used for testing
err_4_t=100.*10**(-6) # MB transverse error, units: meters
err_4_z=150.*10**(-6) # MB z error, units: meters
err_4 = [err_4_t, err_4_z]

#Generates data (tracker hit locations) from particle with specified initial conditions.
def gen_data(x0=0., y0=0., z0=0., phi0=0., pt=3., pz=0.5, q=-1):
    xhit, yhit, zhit = np.asarray(muon_path(x0, y0, z0, phi0, pt, pz, q, plot = False))
    rhit = np.sqrt(xhit*xhit + yhit*yhit)
    phihit = np.arctan2(yhit, xhit)
    return rhit, phihit, zhit

#Smears data by assigning each data point a value based on a normal distribution
#created using the tracking device's standard measuring error.
def smear_func(initial, final, err, rhit, phihit, zhit):
    i = initial
    f = i + final
    r_smear = rhit[i:f]
    phi_smear = np.random.normal(phihit[i:f], err[0]/rhit[i:f])
    x_smear = r_smear * np.cos(phi_smear)
    y_smear = r_smear * np.sin(phi_smear)
    z_smear = np.random.normal(zhit[i:f], err[1])
    return f, x_smear, y_smear, z_smear

#Calculates the actual hits of the particle in case its initial position is
#beyond the innermost tracker.
def determine_start(rhit):
    r_pixel_new = []
    r_tib_new = []
    r_tob_new = []
    r_drift_new = []
    for j in range(len(r_pixel)):
        if rhit[0]<(r_pixel[j]+0.01):
            r_pixel_new = np.append(r_pixel_new, r_pixel[j])
    if len(r_pixel_new)>0:
        return r_pixel_new, r_tib, r_tob, r_drift
    else:
        for j in range(len(r_tib)):
            if rhit[0]<(r_tib[j]+0.01):
                r_tib_new = np.append(r_tib_new, r_tib[j])
        if len(r_tib_new)>0:
            return r_pixel_new, r_tib_new, r_tob, r_drift
        else:
            for j in range(len(r_tob)):
                if rhit[0]<(r_tob[j]+0.01):
                    r_tob_new = np.append(r_tob_new, r_tob[j])
            if len(r_tob_new)>0:
                return r_pixel_new, r_tib_new, r_tob_new, r_drift
            else:
                for j in range(len(r_drift)):
                    if rhit[0]<(r_drift[j]+0.01):
                        r_drift_new = np.append(r_drift_new, r_drift[j])
                return r_pixel_new, r_tib_new, r_tob_new, r_drift_new

#Smears each data point in the appropriate way (data from multiple trackers).
def smear_data(rhit, phihit, zhit, tube):
    x_smear, y_smear, z_smear = [],[],[]
    err_layers = [err_1, err_2, err_3, err_4]
    det_layers = [[], r_pixel, r_tib, r_tob, r_drift]
    if rhit[0]>r_pixel[0]:
        temporary = determine_start(rhit)
        det_layers = [[], temporary[0], temporary[1], temporary[2], temporary[3]]
        r_tracker2 = np.concatenate((det_layers[0], det_layers[1], det_layers[2], det_layers[3]))
    temp = [0]
    if tube == 1:
        for j in range(len(det_layers)-2):
            temp = smear_func(temp[0], len(det_layers[j+1]), err_layers[j], rhit, phihit, zhit)
            x_smear = np.append(x_smear, temp[1])
            y_smear = np.append(y_smear, temp[2])
            z_smear = np.append(z_smear, temp[3])
    if tube == 2:
        if rhit[0] > r_pixel[0]: temp = smear_func(len(r_tracker2), len(det_layers[4]), err_layers[3], rhit, phihit, zhit)
        else: temp = smear_func(len(r_tracker), len(det_layers[4]), err_layers[3], rhit, phihit, zhit)
        x_smear = np.append(x_smear, temp[1])
        y_smear = np.append(y_smear, temp[2])
        z_smear = np.append(z_smear, temp[3])
    return x_smear, y_smear, z_smear

#The following functions are used to create a least squares circle fit.
#======================================================================#
def calc_R(x,y, xc, yc):
    # Calculate the distance of each 2D points from the centre (xc, yc).
    return np.sqrt((x-xc)**2 + (y-yc)**2)

def f(c, x, y):
    # Calculate the algebraic distance between the data points and the mean circle centred at c=(xc, yc).
    Ri = calc_R(x, y, *c)
    return Ri - Ri.mean()

def leastsq_circle(x,y):
    # Coordinates of the barycenter
    x_m = np.mean(x)
    y_m = np.mean(y)
    center_estimate = x_m, y_m
    center, ier = optimize.leastsq(f, center_estimate, args=(x,y))
    xc, yc = center
    Ri       = calc_R(x, y, *center)
    R        = Ri.mean()
    #R_err    = Ri.std()/np.sqrt(len(Ri))
    #residu   = np.sum((Ri - R)**2)
    return xc, yc, R

def plot_data_circle(x,y, xc, yc, R):
    f = plt.figure( facecolor='white')  #figsize=(7, 5.4), dpi=72,
    plt.axis('equal')

    theta_fit = np.linspace(-np.pi, np.pi, 180)

    x_fit = xc + R*np.cos(theta_fit)
    y_fit = yc + R*np.sin(theta_fit)
    plt.plot(x_fit, y_fit, 'b-' , label="fitted circle", lw=2)
    plt.plot([xc], [yc], 'bD', mec='y', mew=1)
    plt.xlabel('x')
    plt.ylabel('y')
    # plot data
    plt.plot(x, y, 'r.', label='data', mew=1)

    plt.legend(loc='best',labelspacing=0.1 )
    plt.grid()
    plt.title('Least Squares Circle')
    plt.show()

#===================================================================#

#Calculates transverse momentum by fitting circular path in smeared data.
def pt_calc(rhit, phihit, zhit,tube):
    x_data, y_data, z_data = smear_data(rhit, phihit, zhit, tube)
    if len(x_data)>0: xc, yc, Rc = leastsq_circle(x_data, y_data)
    if len(x_data)==0: Rc=0
    if tube == 1:
        pt = 0.3*B*Rc
    if tube == 2:
        pt = 0.3*B2*Rc
    if tube == 1:
        k = np.polyfit(r_tracker, z_data,1)[0]
    if tube == 2:
        k = np.polyfit(r_drift, z_data,1)[0]
    p = pt*np.sqrt(k*k+1)
   #print Rc
    return pt, p

#Iterates calculation of transverse momentum over a given amount of times,
#and calculates a mean and a standard error for the estimated transverse
#momentum for this particle.
def pt_datapoint(pt,pz, iter_num):
    #rhit, phihit, zhit = gen_data(pt=pt, pz=pz)
    pt_data,pt_data2,p_data,p_data2 = [],[],[],[]
    for i in range(iter_num):
        rhit, phihit, zhit = gen_data(pt=pt, pz=pz)
        pt_data = np.append(pt_data, pt_calc(rhit, phihit, zhit,1)[0])
        pt_data2 = np.append(pt_data2, pt_calc(rhit, phihit, zhit,2)[0])
        p_data = np.append(p_data, pt_calc(rhit, phihit, zhit,1)[1])
        p_data2 = np.append(p_data2, pt_calc(rhit, phihit, zhit,2)[1])
    p_mean = np.mean(p_data)
    p_err = (np.std(p_data))/(np.sqrt(iter_num))
    p_mean2 = np.mean(p_data2) - E_loss #Adjusting for calorimeter energy loss.
    p_err2 = (np.std(p_data2))/(np.sqrt(iter_num))
    pt_mean = np.mean(pt_data)
    pt_err = (np.std(pt_data))/(np.sqrt(iter_num))
    pt_mean2 = np.mean(pt_data2) -(pt_mean/p_mean)*E_loss #Adjusting for calorimeter energy loss.
    pt_err2 = (np.std(pt_data2))/(np.sqrt(iter_num))
    return pt_mean, pt_err, pt_mean2, pt_err2, p_mean, p_err, p_mean2, p_err2

#Iterates previous function over a given range of transverse momenta.
def pt_data(pt_i, pt_f, pz, point_num, iter_num):
    pt, pt_actual, pt_err, pt2, pt_err2, p, p_err, p2, p_err2 = [], [], [], [], [], [], [], [], []
    stepsize = float((pt_f - pt_i))/point_num
    for i in range(point_num):
        temp1, temp2, temp3, temp4, temp5, temp6, temp7, temp8 = pt_datapoint(iter_num = iter_num, pt = pt_i + i*stepsize, pz=pz)
        pt = np.append(pt, temp1)
        pt_err = np.append(pt_err, temp2)
        pt2 = np.append(pt2, temp3)
        pt_err2 = np.append(pt_err2, temp4)
        pt_actual = np.append(pt_actual, pt_i +i*stepsize)
        p = np.append(p,temp5)
        p_err = np.append(p_err,temp6)
        p2 = np.append(p2,temp7)
        p_err2 = np.append(p_err2,temp8)
        #Progress bar.
        sys.stdout.write("\r%d%%" % i)
        sys.stdout.flush()
    return pt, pt_err, pt2, pt_err2, p, p_err, p2, p_err2, pt_actual

#Plots the standard error in the transverse momentum and its deviation from
#the actual value against the transverse momentum.
def plot_data(pt_i=5, pt_f=20, pz=0.5, point_num=100, iter_num=50): #Does not work if muon trapped in tracker.
    pt, pt_err, pt2, pt_err2, p, p_err, p2, p_err2, pt_actual = pt_data(pt_i, pt_f, pz, point_num, iter_num)
    p_actual = np.sqrt(pt_actual*pt_actual +pz*pz)
    if all(pt) == 0: pt_dev = [0]*len(pt)
    else: pt_dev = (pt-pt_actual)
    pt2_dev = (pt2-pt_actual)
    p_dev = (p-p_actual)
    p2_dev = (p2-p_actual)
    ##Plotting transverse momentum.
    #plt.figure(1)
    #plt.subplot(211)
    #plt.title('$p_t$ error vs $p_t$, Tracker')
    #plt.ylabel('$\Delta p_t$')
    #plt.grid()
    #plt.plot(pt_actual, pt_err/pt, 'b.' , label="$p_t$ Error")
    #plt.plot(pt_actual, pt_dev/pt, 'r.' , label="$p_t$ Deviation")
    #plt.legend(loc='best',labelspacing=0.1,numpoints=1)
    #plt.subplot(212)
    #plt.title('$p_t$ error vs $p_t$, DT')
    #plt.xlabel('$p_t$')
    #plt.ylabel('$\Delta p_t$')
    #plt.grid()
    #plt.plot(pt_actual, pt_err2/pt2, 'b.' , label="$p_t$ Error")
    #plt.plot(pt_actual, pt2_dev/pt2, 'r.' , label="$p_t$ Deviation")

    #Plotting fractional momentum.
    plt.figure(1)
    plt.subplot(211)
    plt.ylabel('$\Delta p$')
    plt.grid()
    plt.title('$p$ error vs $p$, Tracker')
    plt.plot(p_actual, p_err/p, 'b.' , label="$p$ Error")
    plt.plot(p_actual, p_dev/p, 'r.' , label="$p$ Deviation")
    plt.legend(loc='best',labelspacing=0.1,numpoints=1)
    plt.subplot(212)
    plt.xlabel('$p$')
    plt.ylabel('$\Delta p$')
    plt.grid()
    plt.title('$p$ error vs $p$, DT')
    plt.plot(p_actual, p_err2/p2, 'b.' , label="$p$ Error")
    plt.plot(p_actual, p2_dev/p2, 'r.' , label="$p$ Deviation")
    plt.show()

    ##Plotting momentum.
    #plt.figure(2)
    #plt.subplot(211)
    #plt.ylabel('$\Delta p$')
    #plt.grid()
    #plt.title('$p$ error vs $p$, Tracker')
    #plt.plot(p_actual, p_err, 'b.' , label="$p$ Error")
    #plt.plot(p_actual, p_dev, 'r.' , label="$p$ Deviation")
    #plt.legend(loc='best',labelspacing=0.1,numpoints=1)
    #plt.subplot(212)
    #plt.xlabel('$p$')
    #plt.ylabel('$\Delta p$')
    #plt.grid()
    #plt.title('$p$ error vs $p$, DT')
    #plt.plot(p_actual, p_err2, 'b.' , label="$p$ Error")
    #plt.plot(p_actual, p2_dev, 'r.' , label="$p$ Deviation")
    #plt.show()

    #Plotting actual deviations.
    plt.figure(2)
    plt.ylabel('$\Delta p$')
    plt.grid()
    plt.title('$p$ deviation from actual value')
    plt.plot(p_actual, p_dev/p_actual, 'b.' , label="$p$ Deviation, Tracker")
    plt.plot(p_actual, p2_dev/p_actual, 'r.' , label="$p$ Deviation, DT")
    plt.legend(loc='best',labelspacing=0.1,numpoints=1)
    plt.show()

#Plots histograms of momentum measurements for specific actual momentum, useful to deduce distributions.
def poutsoktonos(pt,pz, iter_num):
    pt_tracker_list, pt_dt_list = [],[]
  #  rhit, phihit, zhit = gen_data(pt=pt, pz=pz)
    for i in range(iter_num):
        rhit, phihit, zhit = gen_data(pt=pt, pz=pz)
        temp = pt_calc(rhit,phihit,zhit,1)
        pt_dt = pt_calc(rhit,phihit,zhit,2)[0] - (temp[0]/temp[1])*E_loss

        pt_tracker_list = np.append(pt_tracker_list,pt_calc(rhit,phihit,zhit,1)[0])
        pt_dt_list = np.append(pt_dt_list, pt_dt)

        #Progress bar.
        sys.stdout.write("\r%d%%" % i)
        sys.stdout.flush()
    bins=np.histogram(np.hstack((pt_tracker_list,pt_dt_list)), bins=100)[1]
    plt.hist(pt_tracker_list,bins=bins, color='b', alpha=1, label='Tracker Data')
    plt.hist(pt_dt_list, bins=bins, color='r', alpha=0.5, label='Drift Tube Data')
    plt.title("$p_t$ Histogram")
    plt.xlabel("$p_t$ / $GeVc^{-1}$")
    plt.ylabel("Frequency")
    plt.axvline(x=pt, color='k', linewidth=2,ls='dashed', label='Actual Transverse Momentum')
    plt.legend()
    print 'pt tracker mean' , np.mean(pt_tracker_list)
    print 'pt drift tube mean' , np.mean(pt_dt_list)
    plt.show()

#=============================================#
#IMPORTING DATA#

def import_data():
    raw = np.fromfile('C:\Users\Nikos\Documents\GitHub\Muon-Tracking\muon_list.txt', sep=' ')
    raw_split = np.split(raw,40000)
    return raw_split

def pt_calc2(rhit, phihit, zhit,tube):
    x_data, y_data, z_data = smear_data(rhit, phihit, zhit, tube)
    if len(x_data)>0: xc, yc, Rc = leastsq_circle(x_data, y_data)
    if len(x_data)==0: Rc=0
    if tube == 1:
        pt = 0.3*B*Rc
    if tube == 2:
        pt = 0.3*B2*Rc
    if tube == 1:
        k = np.polyfit(r_tracker, z_data,1)[0]
    if tube == 2:
        k = np.polyfit(r_drift, z_data,1)[0]
    pz= pt*k
    return pt, pz

def pt_datapoint2(pt,pz,phi,charge, iter_num):
    rhit, phihit, zhit = gen_data(pt=pt, pz=pz, phi0=phi, q=charge)
    pt_data,pt_data2,pz_data,pz_data2 = [],[],[],[]
    for i in range(iter_num):
        pt_data = np.append(pt_data, pt_calc2(rhit, phihit, zhit,1)[0])
        pt_data2 = np.append(pt_data2, pt_calc2(rhit, phihit, zhit,2)[0])
        pz_data = np.append(pz_data, pt_calc2(rhit, phihit, zhit,1)[1])
        pz_data2 = np.append(pz_data2, pt_calc2(rhit, phihit, zhit,2)[1])
    phi0 = phihit[0]
    if phihit[1]>phihit[0] : q = -1
    if phihit[1]<phihit[0] : q = 1
    p2 = np.sqrt(pt_data2*pt_data2 + pz_data2*pz_data2)
    p2_mean = np.mean(p2)
    pz_mean = np.mean(pz_data)
    pz_mean2 = np.mean(pz_data2) - (pz_mean/p2_mean)*E_loss #Adjusting for calorimeter energy loss.
    pt_mean = np.mean(pt_data)
    pt_mean2 = np.mean(pt_data2) -(pt_mean/p2_mean)*E_loss #Adjusting for calorimeter energy loss.
    return pt_mean, pt_mean2, pz_mean, pz_mean2, phi0, q

def find_fourvector(pt_d,pz_d,phi,charge,iter_num):
    pt, pt2, pz, pz2, phi0, q = pt_datapoint2(pt_d,pz_d,phi,charge,iter_num)
    p_d= np.sqrt(pt_d*pt_d + pz_d*pz_d)
    if p_d <200:
        p= np.sqrt(pt*pt + pz*pz)
        px = pt*np.cos(phi0)
        py = pt*np.sin(phi0)
        vector = np.array([p, px, py, pz])
    if p_d >= 200:
        p= np.sqrt(pt2*pt2 + pz2*pz2)
        px = pt2*np.cos(phi0)
        py = pt2*np.sin(phi0)
        vector = np.array([p, px, py, pz2])
    return vector

def pt_d_calc(data, i):
    pt = np.sqrt((data[i])[1]**2 +(data[i])[2]**2)
    return pt

def phi0_d_calc(data,i):
    phi0 = np.arctan2((data[i])[2],(data[i])[1])
    return phi0

def find_higgs(points, iter_num):
    mass_list, mass_listz1, mass_listz2 = [],[],[]
    data = import_data()
    for i in range(points):
        muon1 = find_fourvector(pt_d_calc(data,4*i),(data[4*i])[3],phi0_d_calc(data,4*i), (data[4*i])[0], iter_num)
        muon2 = find_fourvector(pt_d_calc(data,4*i+1),(data[4*i+1])[3],phi0_d_calc(data,4*i+1), (data[4*i+1])[0], iter_num)
        muon3 = find_fourvector(pt_d_calc(data,4*i+2),(data[4*i+2])[3],phi0_d_calc(data,4*i+2), (data[4*i+2])[0], iter_num)
        muon4 = find_fourvector(pt_d_calc(data,4*i+3),(data[4*i+3])[3],phi0_d_calc(data,4*i+3), (data[4*i+3])[0], iter_num)

        z_boson1 = muon1 + muon2
        z_boson2 = muon3 + muon4
        mass_z1 = np.sqrt((z_boson1[0])**2 - ((z_boson1[1])**2 +(z_boson1[2])**2 +(z_boson1[3])**2))
        #print mass_z1
        #print z_boson1
        #print z_boson2
        mass_listz1 = np.append(mass_listz1, mass_z1)
        mass_z2 = np.sqrt((z_boson2[0])**2 - ((z_boson2[1])**2 +(z_boson2[2])**2 +(z_boson2[3])**2))
        #print mass_z2
        mass_listz2 = np.append(mass_listz2, mass_z2)

        higgs = z_boson1 + z_boson2
        #print higgs
        mass_h = np.sqrt((higgs[0])**2 - ((higgs[1])**2 +(higgs[2])**2 +(higgs[3])**2))
        #print mass_h
        mass_list = np.append(mass_list, mass_h)

        #Progress bar.
        sys.stdout.write("\r%d%%" % i)
        sys.stdout.flush()
    mass_mean = np.mean(mass_list)
    mass_mean_z1 = np.mean(mass_listz1)
    mass_mean_z2 = np.mean(mass_listz2)
    return mass_list, mass_listz1, mass_listz2

def plot_hist_higgs(points, iter_no):
    h_list, z1_list, z2_list = find_higgs(points, iter_no)
    plt.figure(1)
    bins=np.histogram(z1_list, bins=50)[1]
    plt.hist(z1_list,bins=bins, color='b',normed = False, alpha=1, label='Z Boson')
    #plt.hist(z2_list, color='r',normed = False, alpha=0.5, label='Z Boson 2')
    plt.title("Z Mass Histogram")
    plt.xlabel("$m$ / $GeVc^{-2}$")
    plt.ylabel("Frequency")
    plt.legend()
    plt.show()

    plt.figure(2)
    plt.hist(h_list, color='g',normed = False, alpha=1, label='Higgs')
    plt.title("Higgs Mass Histogram")
    plt.xlabel("$m$ / $GeVc^{-2}$")
    plt.ylabel("Frequency")
    plt.axvline(x=126, color='k', linewidth=2,ls='dashed', label='Actual Higgs Mass')
    plt.legend()
    print 'higgs mean' , np.mean(h_list)
    print 'z1 and z2 mean' , np.mean(z1_list), np.mean(z2_list)
    plt.show()
    return h_list, z1_list

reason: from scipy


text:
import os
from argparse import ArgumentParser

import cv2
import numpy

from capture.capture import VideoCapture
from face_detector import FaceDetection
from face_embedding import FaceEmbedding
from util import cosine_similarity, draw_square, show_image

RED = (255, 0, 0)
GREEN = (0, 255, 0)

def main():
    args = ArgumentParser()
    args.add_argument('-c', '--camera_url', default=0, type=str, help='0 - local camera')
    args.add_argument('-dt', '--detect_threshold', default=0.975, type=float, help="Threshold of face detection")
    args.add_argument('-rf', '--recognized_threshold', default=0.8, type=float, help="Threshold of face recognition")
    args.add_argument('--device', default='cuda:0', type=str, help="Device run model. `cuda:<id>` or `cpu`")
    args.add_argument('--detect_face_model', default='data/pretrained/mobilenet_header.pth',
                      type=str, help="Face detector model path")
    args.add_argument('--detect_face_backbone', default='data/pretrained/mobile_backbone.tar',
                      type=str, help="Face detector backbone path")
    args.add_argument('--recognized_model', default='data/pretrained/embedder_resnet50_asia.pth'
                      , type=str, help="Face embedding model path")
    args.add_argument('--model_registered', default='model_faces.npy', type=str, help="Model contain face's vectors")
    args.add_argument('--model_ids', default='model_face_ids.npy', type=str, help="Model contain face's ids")
    args = args.parse_args()

    try:
        args.camera_url = int(args.camera_url)
    except:
        pass

    if not (os.path.isfile(args.model_registered) and os.path.isfile(args.model_ids)):
        face_model = numpy.zeros((0, 512), dtype=numpy.float32)
        ids_model = []
    else:
        face_model = numpy.load(args.model_registered, allow_pickle=True)
        ids_model = numpy.load(args.model_ids, allow_pickle=True).tolist()

    detector = FaceDetection(args.detect_face_model, args.detect_face_backbone, scale_size=480, device=args.device)
    embedder = FaceEmbedding(args.recognized_model, device=args.device)

    # recognize
    video = VideoCapture(args.camera_url)

    for frame in video:
        faces = detector(frame)
        faces = embedder(faces)

        for face in faces:
            txt = "None"
            color = RED

            scores = cosine_similarity(face.embedding.reshape(1, 512), face_model, skip_normalize=True).ravel()
            args_idx = numpy.argmax(scores)

            if scores[args_idx] >= args.recognized_threshold:
                txt = ids_model[args_idx]
                color = GREEN

            frame = draw_square(frame, face.box.astype(numpy.int), color=color)
            frame = cv2.putText(frame, f"EID: {txt}",
                                (int(face.box[0]), int(face.box[1] - 20)), cv2.FONT_HERSHEY_PLAIN, 1,
                                GREEN)

        if not show_image(frame, 'Face Recognition', windows_size=(1920, 1080)):
            break

    video.stop()

if __name__ == '__main__':
    main()


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

# __authors__ = "Adrien Guille, Antoine Gourru"
# __email__ = "adrien.guille@univ-lyon2.fr, antoine.gourru@univ-lyon2.fr"

import pandas as pd
import numpy as np
import scipy as sp

from utils import read_data,evaluate, find_optimal_C

from sklearn.decomposition import TruncatedSVD
from sklearn.preprocessing import normalize

from gensim.models import Word2Vec

from models import train_rle,train_w2v,train_geld,prepare_data

ratios = 0.5
methods = ["RLE","GELD"]#"RLE","GELD"
dataset= "cora2"
d = 160

print("\nLoading %s..." % dataset)

tf_idf, groups, A, graph,voc,raw,tf = read_data(dataset)
n = A.shape[0]
density = 100 * 2 * len(list(graph.edges())) / (n * (n - 1))
degrees = np.array([degree for _, degree in graph.degree()])
mean_degree = degrees.mean()
print("%d nodes, %d features, %d edges (density: %.3f%%, mean degree: %.1f)" % (tf_idf.shape[0],
                                                                                tf_idf.shape[1],
                                                                                len(list(graph.edges())),
                                                                                density,
                                                                                mean_degree))
embeddings = {}

print("\nTraining embeddings in dimension %d..." % d)

# -------------------------
# Learn the representations
# -------------------------
data_graph,data_text,sigma_init,D_init,U = prepare_data(d,tf,A,voc,raw)

if "GELD" in methods:
    optimal_alpha = {"cora2": 0.99, "dblp": 0.8,"nyt":0.95}
    D,D_norm,sigma = train_geld(d,data_graph,data_text,U,D_init,sigma_init,n_epoch=40,lamb=None,alpha=optimal_alpha[dataset])
    embeddings["GELD"] = D_norm

if "RLE" in methods:
    optimal_window = {"cora2": 15, "dblp": 5,"nyt":10}
    embeddings["RLE"], _ = train_rle(A,tf,U,d,0.7,verbose=True)

# ------------------------------------------------------------------------------
# Train and evaluate the SVM for test/train ratios in [0.5, 0.9]
# ------------------------------------------------------------------------------

print("\nEvaluating embeddings...")
ratios = [0.9,0.5]
results = pd.DataFrame(ratios, columns=["ratio"])
for method in methods:
    optimal_C = find_optimal_C(embeddings[method], groups)
    results[method] = results["ratio"].apply(lambda ratio: evaluate(embeddings[method],
                                                                    groups,
                                                                    ratio,
                                                                    C=optimal_C,
                                                                    verbose=False))

results2 = pd.DataFrame(ratios, columns=["ratio"])

for method in methods:
    #print(method+"%.1f",results[method].tolist()[0])
    results2[[method+"_accuracy_mean", method+"_accuracy_std"]] = pd.DataFrame(results[method].tolist(), index=results.index)
print(results2.T)

results2.to_csv(dataset+"_classification.csv")

max_stars_repo_name: AntoineGourru/GELD

max_issues_repo_name: AntoineGourru/GELD

max_forks_repo_name: AntoineGourru/GELD

reason: import numpy,import scipy


text:
import os
os.environ['TF_CPP_MIN_LOG_LEVEL'] = '3' #supress tensorflow messages
os.environ['FOR_DISABLE_CONSOLE_CTRL_HANDLER'] = '1' #prevent error 200 in console

import numpy as np
import time

from tensorflow.keras.utils import to_categorical
from tensorflow.keras.optimizers import Adam
from sklearn.model_selection import train_test_split

from random import seed as set_py_seed
from numpy.random import seed as set_np_seed
from tensorflow.random import set_seed as set_tf_seed

from utils import play_sound, calc_accuracy
from wifinet import WIFINet
from ann import ANN
from cnn import CNN
from mlmodels import KNN, SVM, Tree, Forest, LogReg, GNB, AdaBoost

def main():
    '''Loads data and train various models.'''
    train, test, train_labels, test_labels = load_data(normalize = True, norm_type = 1) #load data
    #train, test, train_labels, test_labels = np.random.randint(0, 90, size=(1000, 20000, 1)), np.random.randint(0, 90, size=(100, 20000, 1)), np.random.randint(0, 7, size=(1000,)), np.random.randint(0, 7, size=(100,)) #fake data for faster testing

    #train 6 wifinet models using various sample lengths
    slen = [5000, 10000, 15000, 20000, 25000, 30000]
    for i in range(len(slen)):
        wnet = WIFINet((slen[i], 1), (7, ), dilation_depth = 4, num_filters = 32, O1 = 32, O2 = 10, D1 = 1024, D2 = 128) #create model
        opt = Adam(learning_rate = 1e-4, beta_1 = 0.9, beta_2 = 0.999, epsilon = 1e-07, amsgrad = False) #setup optimizer
        wnet.fit(train[:, 0:slen[i], :], to_categorical(train_labels), validation_data = (test[:, 0:slen[i], :], to_categorical(test_labels)),
                 epochs = 10, batch_size = 32, optimizer = opt, save = True, model_name = 'wifinet_L{}'.format(slen[i]), hist_name = 'wifinet_history_L{}'.format(slen[i])) #fit model

    #train 3 additional wifinet models of various depths
    slen = 20000
    depth = [2, 8, 10]
    for i in range(len(depth)):
        wnet = WIFINet((slen, 1), (7, ), dilation_depth = depth[i], num_filters = 32, O1 = 32, O2 = 10, D1 = 1024, D2 = 128) #create model
        opt = Adam(learning_rate = 1e-4, beta_1 = 0.9, beta_2 = 0.999, epsilon = 1e-07, amsgrad = False) #setup optimizer
        wnet.fit(train[:, 0:slen, :], to_categorical(train_labels), validation_data = (test[:, 0:slen, :], to_categorical(test_labels)),
                 epochs = 10, batch_size = 32, optimizer = opt, save = True, model_name = 'wifinet_D{}'.format(depth[i]), hist_name = 'wifinet_history_D{}'.format(depth[i])) #fit model

    #train 6 ANNs of various depths
    depth = [[1024, 128], [2048, 128], [4096, 128], [1024, 512, 128], [2048, 512, 128], [4096, 512, 128]]
    for i in range(len(depth)):
        ann = ANN((slen, ), (7, ), depth = depth[i]) #create model
        opt = Adam(learning_rate = 1e-4, beta_1 = 0.9, beta_2 = 0.999, epsilon = 1e-07, amsgrad = False) #setup optimizer
        ann.fit(train[:, 0:slen, :], to_categorical(train_labels), validation_data = (test[:, 0:slen, :], to_categorical(test_labels)),
                epochs = 10, batch_size = 32, optimizer = opt, save = True, model_name = 'ann_D{}'.format(depth[i]), hist_name = 'ann_history_D{}'.format(depth[i])) #fit model

    #train 3 CNNs of various depths
    depth = [[1, 2, 1], [1,3,2], [2, 5, 3]]
    for i in range(len(depth)):
        cnn = CNN((slen, 1), (7, ), depth = depth[i], filters = [16, 32, 64], kernels = [3, 5, 3], D1 = 1024, D2 = 128) #create model
        opt = Adam(learning_rate = 1e-4, beta_1 = 0.9, beta_2 = 0.999, epsilon = 1e-07, amsgrad = False) #setup optimizer
        cnn.fit(train[:, 0:slen, :], to_categorical(train_labels), validation_data = (test[:, 0:slen, :], to_categorical(test_labels)),
                epochs = 10, batch_size = 32, optimizer = opt, save = True, model_name = 'cnn_D{}'.format(depth[i]), hist_name = 'cnn_history_D{}'.format(depth[i])) #fit model

    #train 12 KNN models of various number of neighbors
    neighbors = [1, 2, 3, 5, 7, 10, 15, 20, 25, 30, 40, 50]
    for i in range(len(neighbors)):
        knn = KNN(num_neighbors = neighbors[i])
        knn.fit(np.squeeze(train)[:, 0:slen], train_labels)
        knn.save(path = './model/knn_{}.pkl'.format(neighbors[i]))

    #run prediction for all KNN models on 1000 samples only, too long otherwise even on 24 cores for this many models
    num_samples = 1000
    for i in range(len(neighbors)):
        knn = KNN.load('./model/knn_{}.pkl'.format(neighbors[i]))
        pred, time = knn.predict(np.squeeze(test)[0:num_samples, 0:slen])
        print('The prediction took {} seconds to predict {} samples.'.format(time, num_samples))
        print('Accuracy: {}%'.format(calc_accuracy(np.argmax(pred, axis = 1), test_labels[0:num_samples])))

    #train SVM model and run prediction, wrapped with CalibratedClassifierCV.
    svm = SVM()
    svm.fit(np.squeeze(train)[:, 0:slen], train_labels)
    svm.save()
    pred, time = svm.predict(np.squeeze(test)[:, 0:slen])
    print('The prediction took {} seconds.'.format(time))
    print('Accuracy: {}%'.format(calc_accuracy(np.argmax(pred, axis = 1), test_labels)))

    #train tree model and run prediction
    leaves = [1, 2, 3, 4, 5, 10, 15, 20, 25, 30, 40, 50, 75, 100]
    for i in range(len(leaves)):
        tree = Tree(leaves[i])
        tree.fit(np.squeeze(train)[:, 0:slen], train_labels)
        tree.save(path = './model/tree_{}.pkl'.format(leaves[i]))
        pred, time = tree.predict(np.squeeze(test)[:, 0:slen])
        print('The prediction took {} seconds.'.format(time))
        print('Accuracy: {}%'.format(calc_accuracy(np.argmax(pred, axis = 1), test_labels)))

    #train 15 forest models with cvarious numbers of trees and run prediction
    trees = [1, 5, 10, 15, 20, 25, 30, 40, 50, 75, 100, 150, 200, 250, 300]
    for i in range(len(trees)):
        forest = Forest(n_trees = trees[i])
        forest.fit(np.squeeze(train)[:, 0:slen], train_labels)
        forest.save(path = './model/forest_{}.pkl'.format(trees[i]))
        pred, time = forest.predict(np.squeeze(test)[:, 0:slen])
        print('The prediction took {} seconds.'.format(time))
        print('Accuracy: {}%'.format(calc_accuracy(np.argmax(pred, axis = 1), test_labels)))

    #train logistic regression model and run prediction
    logreg = LogReg()
    logreg.fit(np.squeeze(train)[:, 0:slen], train_labels)
    logreg.save()
    pred, time = logreg.predict(np.squeeze(test)[:, 0:slen])
    print('The prediction took {} seconds.'.format(time))
    print('Accuracy: {}%'.format(calc_accuracy(np.argmax(pred, axis = 1), test_labels)))

    #train guassian naive bayes model and run prediction
    gnb = GNB()
    gnb.fit(np.squeeze(train)[:, 0:slen], train_labels)
    gnb.save()
    pred, time = gnb.predict(np.squeeze(test)[:, 0:slen])
    print('The prediction took {} seconds.'.format(time))
    print('Accuracy: {}%'.format(calc_accuracy(np.argmax(pred, axis = 1), test_labels)))

    #train 13 adaboost ensemble model and run prediction
    ests = [1, 5, 10, 20, 30, 40, 50, 75, 100, 200, 300]
    for i in range(len(ests)):
        ada = AdaBoost(n_estimators = ests[i])
        ada.fit(np.squeeze(train)[:, 0:slen], train_labels)
        ada.save(path = './model/ada_{}.pkl'.format(ests[i]))
        pred, time = ada.predict(np.squeeze(test)[:, 0:slen])
        print('The prediction took {} seconds.'.format(time))
        print('Accuracy: {}%'.format(calc_accuracy(np.argmax(pred, axis = 1), test_labels)))

def load_data(normalize = True, norm_type = 1, for_training = True):
    '''Load data. Normalize if true using specified method. If data is for training, expand dims and split into train and test sets.'''
    path = ['data/data.npy', 'data/labels.npy']
    print('Loading data from {} and {} ...'.format(path[0], path[1]))
    data = np.load(path[0])[:, 0, :] #load data and keep only I component
    info = np.load(path[1]) #load labels

    #normalize data
    if normalize:
        if norm_type == 1: #offset by threshold to move to zero and scale to [-1, 1] while preserving zero location
            threshold, abs_max = 59.0, np.max(np.abs(data))
            print('Normalizing data using threshold offset {} and absolute max value {} ...'.format(threshold, abs_max))
            data = (data + threshold) / abs_max
        elif norm_type == 2: #offset by threshold and an scale to [0, 1]
            x_min, x_max = np.min(data), np.max(data)
            print('Normalizing data using min value {} and max value {} ...'.format(x_min, x_max))
            data = (data - x_min) / (x_max - x_min)
        elif norm_type == 3: #offset by threshold and scale to [-1, 1]
            x_min, x_max = np.min(data), np.max(data)
            print('Normalizing data using min value {} and max value {} ...'.format(x_min, x_max))
            data = 2 * (data - x_min) / (x_max - x_min) - 1

    if for_training:
        #expand dims
        print('Expanding dims and splitting into train and test sets ...')
        data = np.expand_dims(data, axis = 2)
        #split dataset into train and test
        labels = info[:, 0]
        train, test, train_labels, test_labels = train_test_split(data, labels, test_size = 0.20, random_state = 42, shuffle = True, stratify = labels)
        print('Train Shape:', train.shape, 'Train Labels Shape:', train_labels.shape, 'Test Shape:', test.shape, 'Test Labels Shape:', test_labels.shape)
        return train, test, train_labels, test_labels
    else:
        print('Data Shape:', data.shape, 'Labels Shape:', info.shape)
        return data, info

def set_seeds(seed=42):
    '''Set all seeds to make training deterministic.'''
    os.environ['PYTHONHASHSEED'] = str(seed) #set python hash seed
    os.environ['TF_CUDNN_DETERMINISM'] = '1'
    os.environ['TF_DETERMINISTIC_OPS'] = '1'
    set_py_seed(seed) #set python seed
    set_np_seed(seed) #set numpy seed
    set_tf_seed(seed) #set tensorflow seed

if __name__ == '__main__':
    try:
        set_seeds()
        main()
        play_sound()
    except Exception as err:
        print('exception occurred ...')
        print(err)
        play_sound()

max_stars_repo_name: rjalnev/WIFINET

max_issues_repo_name: rjalnev/WIFINET

max_forks_repo_name: rjalnev/WIFINET

reason: import numpy,from numpy


text:
# encoding:utf-8
# from dataset.data import read_image  # 图片读取方法，可以自己写，我是用的baseline里自带的
from data.datasets.dataset_loader import read_image
import os
import torch
import numpy as np
import json
from data.datasets.eval_reid import eval_func
import torchvision.transforms as T
from modeling.aligned_baseline import AlignedBaseline
from modeling.pcb_baseline import PCBBaseline
from utils.re_ranking import re_ranking
from utils.aligned_reranking import aligned_re_ranking
from torchvision import transforms
# from sklearn.decomposition import PCA
import torch.nn as nn

os.environ['CUDA_VISIBLE_DEVICES'] = '2,6'  # 指定gpu
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
from config import cfg
from data import make_data_loader
import re

def normalize(nparray, order=2, axis=0):
    """Normalize a N-D numpy array along the specified axis."""
    norm = np.linalg.norm(nparray, ord=order, axis=axis, keepdims=True)
    return nparray / (norm + np.finfo(np.float32).eps)

def compute_dist(array1, array2, type='euclidean'):
    """Compute the euclidean or cosine distance of all pairs.
    Args:
      array1: numpy array with shape [m1, n]
      array2: numpy array with shape [m2, n]
      type: one of ['cosine', 'euclidean']
    Returns:
      numpy array with shape [m1, m2]
    """
    # array1 = array1.cpu().numpy()
    # array2 = array2.cpu().numpy()
    assert type in ['cosine', 'euclidean']
    if type == 'cosine':
        array1 = normalize(array1, axis=1)
        array2 = normalize(array2, axis=1)
        dist = np.matmul(array1, array2.T)
        return dist
    else:
        # shape [m1, 1]
        square1 = np.sum(np.square(array1), axis=1)[..., np.newaxis]
        # shape [1, m2]
        square2 = np.sum(np.square(array2), axis=1)[np.newaxis, ...]
        squared_dist = - 2 * np.matmul(array1, array2.T) + square1 + square2
        squared_dist[squared_dist < 0] = 0
        dist = np.sqrt(squared_dist)
        return dist

def euclidean_dist(x, y):
    m, n = x.size(0), y.size(0)
    xx = torch.pow(x, 2).sum(1, keepdim=True).expand(m, n)
    yy = torch.pow(y, 2).sum(1, keepdim=True).expand(n, m).t()
    dist = xx + yy
    dist.addmm_(1, -2, x, y.t())
    dist = dist.clamp(min=1e-12).sqrt()  # for numerical stability
    return dist

def re_rank(qf, gf):
    return re_ranking(qf, gf, k1=6, k2=2, lambda_value=0.4)

def transform(img):
    normalize_transform = T.Normalize(mean=[0.485, 0.456, 0.406], std=[0.229, 0.224, 0.225])
    result = T.Compose([
        T.Resize([384, 128]),
        T.ToTensor(),
        normalize_transform
    ])
    result = result(img)
    return result

def compute_qg_dist(all_feature, query_num):
    gf = all_feature[query_num:]
    qf = all_feature[:query_num]

    q_g_dist = compute_dist(
        qf.cpu().detach().numpy(), gf.cpu().detach().numpy(), type='euclidean')
    g_g_dist = compute_dist(
        gf.cpu().detach().numpy(), gf.cpu().detach().numpy(), type='euclidean')
    q_q_dist = compute_dist(
        qf.cpu().detach().numpy(), qf.cpu().detach().numpy(), type='euclidean')
    return q_g_dist, q_q_dist, g_g_dist

def get_model(traindata_num_classes, modelname):
    pretrain_path = r"/home/zxh/reid-baseline/new_experiment/r101_pcb_new/resnext101_ibn_a_model_70.pth"
    model = PCBBaseline(traindata_num_classes, 1, pretrain_path, 'bnneck', 'after', 'resnext101_ibn_a', 'self', 6, False,
                        False, False, False, False, False)
    model.load_param(modelname)
    return model

def write_json(dist, query_name, gallery_name, save_dir='', topk=10):
    if not osp.exists(save_dir):
        os.makedirs(save_dir)
    num_q, num_g = dist.shape

    print(dist.shape)

    print('# query: {}\n# gallery {}'.format(num_q, num_g))
    print('Writing top-{} ranks ...'.format(topk))

    num_query, num_gallery = dist.shape
    print(dist.shape)
    dist_line = dist.reshape(-1)
    print(dist_line.shape)
    dist_sorted_indices = np.argsort(dist_line)
    flags = np.zeros(len(gallery_name))
    result = {}
    for i in range(len(dist)):
        qimg_path, qpid, qcamid = query_name[i]
        qimg_name = self.parse_filename(osp.basename(qimg_path))
        result[qimg_name] = []
    dist_cp = dist.copy()
    dist_cp.sort(1)
    dist_r1 = dist_cp[:, 0]
    rank1 = np.argsort(dist_r1)
    dist_r1.sort()
    thr = dist_r1[int(len(rank1) * 0.72)]  # TODO
    print(thr)
    num = 0
    num_jump = 0
    for i in tqdm(range(len(dist_sorted_indices))):
        index_query = dist_sorted_indices[i] // num_gallery
        index_gallery = dist_sorted_indices[i] % num_gallery
        d = dist_line[dist_sorted_indices[i]]
        qimg_path, qpid, qcamid = query_name[index_query]
        qimg_name = self.parse_filename(osp.basename(qimg_path))
        gimg_path, gpid, gcamid = gallery_name[index_gallery]
        gimg_name = self.parse_filename(osp.basename(gimg_path))
        if d < thr:
            if flags[index_gallery] == 0:
                flags[index_gallery] = 1
                if len(result[qimg_name]) < topk:
                    num += 1
                    result[qimg_name].append(gimg_name)
            else:
                num_jump += 1
        else:
            if flags[index_gallery] == 0:
                if len(result[qimg_name]) < topk:
                    num += 1
                    result[qimg_name].append(gimg_name)
            else:
                num_jump += 1
        if num >= num_query * topk:
            break
    print(num)
    print(num_jump)
    with open(os.path.join(save_dir, r'submission_second_a.json'), 'w', encoding='utf-8') as f:
        json.dump(result, f)
    print("Json Save Done")

def adjust_rerank_function(dist_list, query_pid_list, gallery_pid_list, query_camid_list, gallery_camid_list):
    q_pids = np.asarray(query_pid_list)
    g_pids = np.asarray(gallery_pid_list)
    q_camids = np.asarray(query_camid_list)
    g_camids = np.asarray(gallery_camid_list)

    max = 0
    plist = []
    global_q_g_dist, global_q_q_dist, global_g_g_dist = dist_list[0:3]
    local_q_g_dist, local_q_q_dist, local_g_g_dist = dist_list[3:6]
    flip_global_q_g_dist, flip_global_q_q_dist, flip_global_g_g_dist = dist_list[6:9]
    flip_local_q_g_dist, flip_local_q_q_dist, flip_local_g_g_dist = dist_list[9:12]

    for k1 in range(6, 8, 1):
        for k2 in range(3, 5, 1):
            for l in [0.77, 0.78, 0.79, 0.80, 0.81, 0.82, 0.83, 0.84, 0.85, 0.86, 0.87, 0.88]:  #
                for l_w in [0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97]:

                    distmat_global = aligned_re_ranking(
                        global_q_g_dist, global_q_q_dist, global_g_g_dist, k1=k1, k2=k2,
                        lambda_value=l)
                    print("Global dismat is computed done")

                    distmat_local = aligned_re_ranking(
                        local_q_g_dist, local_q_q_dist, local_g_g_dist, k1=k1, k2=k2,
                        lambda_value=l)

                    print("Local dismat is computed done")

                    flip_distmat_global = aligned_re_ranking(
                        flip_global_q_g_dist, flip_global_q_q_dist, flip_global_g_g_dist, k1=k1, k2=k2,
                        lambda_value=l)

                    print("Global_flip dismat is computed done")

                    flip_distmat_local = aligned_re_ranking(
                        flip_local_q_g_dist, flip_local_q_q_dist, flip_local_g_g_dist, k1=k1, k2=k2,
                        lambda_value=l)

                    print("Local_flip dismat is computed done")

                    distmat = l_w * distmat_global + (1 - l_w) * distmat_local

                    flip_distmat = l_w * flip_distmat_global + (1 - l_w) * flip_distmat_local

                    distmat = (flip_distmat + distmat) / 2

                    cmc, mAP = eval_func(distmat, q_pids, g_pids, q_camids, g_camids)
                    for r in [1]:
                        if max < (mAP + cmc[r - 1]) / 2:
                            max = (mAP + cmc[r - 1]) / 2
                            plist = [k1, k2, l, mAP, cmc[r - 1]]
                        print("====k1=%d=====k2=%d=====l=%f=====l_w=%f" % (k1, k2, l, l_w))
                        print("CMC curve, Rank-%d:%.4f, map:%.4f, final: %.4f" % (
                            r, cmc[r - 1], mAP, (mAP + cmc[r - 1]) / 2))

    print(plist)

def inference_samples(model, batch_size, query_list, gallery_list, query_pid_list, gallery_pid_list, query_camid_list, \
                          gallery_camid_list, adjust_rerank, npy_save_path):  # 传入模型，数据预处理方法，batch_size
    query_num = len(query_list)
    model = nn.DataParallel(model)
    model = model.to(device)
    model.eval()
    iter_n_query = (len(query_list)) // batch_size
    iter_n_gallery = (len(gallery_list)) // batch_size
    if (len(query_list)) % batch_size != 0:
        iter_n_query += 1
    if (len(gallery_list)) % batch_size != 0:
        iter_n_gallery += 1
    all_feature = list()
    local_feature = list()
    flip_all_feature = list()
    flip_local_feature = list()
    print("ALL_QUERY_BACTH:", str(iter_n_query))
    for i in range(iter_n_query):
        img_list = list()
        flip_img_list = list()
        for q_img in query_list[i * batch_size:(i + 1) * batch_size]:
            q_img = read_image(q_img)
            flip_q_img = transforms.RandomHorizontalFlip(p=1.0)(q_img)
            q_img = transform(q_img)
            flip_q_img = transform(flip_q_img)
            img_list.append(q_img)
            flip_img_list.append(flip_q_img)
        img_data = torch.Tensor([t.numpy() for t in img_list]).to(device)
        flip_img_data = torch.Tensor([t.numpy() for t in flip_img_list]).to(device)
        print("batch ----%d----" % (i))
        batch_data = img_data
        flip_batch_data = flip_img_data
        with torch.no_grad():
            batch_feature = model(batch_data, None)[0].detach().cuda()
            local_batch_feature = model(batch_data, None)[1].detach().cuda()
            all_feature.append(batch_feature)
            local_feature.append(local_batch_feature)

            flip_batch_feature = model(flip_batch_data, None)[0].detach().cuda()
            local_flip_batch_feature = model(flip_batch_data, None)[1].detach().cuda()
            flip_all_feature.append(flip_batch_feature)
            flip_local_feature.append(local_flip_batch_feature)
    print("ALL_GALLERY_BACTH:", str(iter_n_gallery))
    for i in range(iter_n_gallery):
        img_list = list()
        flip_img_list = list()
        for g_img in gallery_list[i * batch_size:(i + 1) * batch_size]:
            g_img = read_image(g_img)
            flip_g_img = transforms.RandomHorizontalFlip(p=1.0)(g_img)
            g_img = transform(g_img)
            flip_g_img = transform(flip_g_img)
            img_list.append(g_img)
            flip_img_list.append(flip_g_img)
        img_data = torch.Tensor([t.numpy() for t in img_list]).to(device)
        flip_img_data = torch.Tensor([t.numpy() for t in flip_img_list]).to(device)
        batch_data = img_data
        flip_batch_data = flip_img_data
        with torch.no_grad():
            print("batch ----%d----" % (i))
            batch_feature = model(batch_data, None)[0].detach().cuda()
            local_batch_feature = model(batch_data, None)[1].detach().cuda()
            all_feature.append(batch_feature)
            local_feature.append(local_batch_feature)

            flip_batch_feature = model(flip_batch_data, None)[0].detach().cuda()
            local_flip_batch_feature = model(flip_batch_data, None)[1].detach().cuda()
            flip_all_feature.append(flip_batch_feature)
            flip_local_feature.append(local_flip_batch_feature)

    all_feature = torch.cat(all_feature, dim=0)
    local_feature = torch.cat(local_feature, dim=0)
    flip_all_feature = torch.cat(flip_all_feature, dim=0)
    flip_local_feature = torch.cat(flip_local_feature, dim=0)

    #Flip Cat
    # all_feature = torch.cat((all_feature, flip_all_feature), dim=1)
    # local_feature = torch.cat((local_feature, flip_local_feature), dim=1)
    ######
    all_feature = torch.nn.functional.normalize(all_feature, dim=1, p=2)
    local_feature = torch.nn.functional.normalize(local_feature, dim=1, p=2)
    flip_all_feature = torch.nn.functional.normalize(flip_all_feature, dim=1, p=2)
    flip_local_feature = torch.nn.functional.normalize(flip_local_feature, dim=1, p=2)

    global_q_g_dist, global_g_g_dist, global_q_q_dist = compute_qg_dist(all_feature, query_num)
    local_q_g_dist, local_g_g_dist, local_q_q_dist = compute_qg_dist(local_feature, query_num)
    flip_global_q_g_dist, flip_global_g_g_dist, flip_global_q_q_dist = compute_qg_dist(flip_all_feature, query_num)
    flip_local_q_g_dist, flip_local_g_g_dist, flip_local_q_q_dist = compute_qg_dist(flip_local_feature, query_num)

    if adjust_rerank:
        dist_list = [global_q_g_dist, global_g_g_dist, global_q_q_dist, local_q_g_dist, local_g_g_dist, local_q_q_dist, \
                      flip_global_q_g_dist, flip_global_g_g_dist, flip_global_q_q_dist, flip_local_q_g_dist,
                      flip_local_g_g_dist, flip_local_q_q_dist]
        adjust_rerank_function(dist_list, query_pid_list, gallery_pid_list, query_camid_list, gallery_camid_list)
    else:
        l_w = 0.95
        k1, k2, l = 6, 3, 0.80
        distmat_global = aligned_re_ranking(
            global_q_g_dist, global_q_q_dist, global_g_g_dist, k1=k1, k2=k2,
            lambda_value=l)
        print("Global dismat is computed done")

        distmat_local = aligned_re_ranking(
            local_q_g_dist, local_q_q_dist, local_g_g_dist, k1=k1, k2=k2,
            lambda_value=l)

        print("Local dismat is computed done")

        flip_distmat_global = aligned_re_ranking(
            flip_global_q_g_dist, flip_global_q_q_dist, flip_global_g_g_dist, k1=k1, k2=k2,
            lambda_value=l)

        print("Global_flip dismat is computed done")

        flip_distmat_local = aligned_re_ranking(
            flip_local_q_g_dist, flip_local_q_q_dist, flip_local_g_g_dist, k1=k1, k2=k2,
            lambda_value=l)

        print("Local_flip dismat is computed done")

        distmat = l_w * distmat_global + (1 - l_w) * distmat_local

        flip_distmat = l_w * flip_distmat_global + (1 - l_w) * flip_distmat_local

        distmat = (flip_distmat + distmat) / 2
        np.save(npy_save_path, distmat)

def merge_npy(npy_paths):
    '''
    需要自动化调整这个函数的写法
    :param npy_paths: 包含所有npypath的list
    :return:
    '''
    npy_dist_1 = np.load(npy_paths[0])
    npy_dist_2 = np.load(npy_paths[1])
    distmat = npy_dist_1+npy_dist_2
    write_json(dist=distmat, query_name=query_list, gallery_name=gallery_list, save_dir='./', topk=200)

def main(traindata_num_classes, modelname, batch_size, merge_npy_paths, adjust_rerank , npy_save_path):
    if len(merge_npy_paths) > 0:
        merge_npy(merge_npy_paths)
    else:
        query_list = [os.path.join(r'/home/zxh/datasets/NAICReID/clean_query', x) for x in
                      # 测试集gallery文件夹
                      os.listdir(r'/home/zxh/datasets/NAICReID/clean_query')]

        gallery_list = [os.path.join(r'/home/zxh/datasets/NAICReID/clean_bounding_box_test', x) for x in
                        # 测试集gallery文件夹
                        os.listdir(r'/home/zxh/datasets/NAICReID/clean_bounding_box_test')]
        pattern = re.compile(r'([-\d]+)_c(\d)')

        query_pid_list = []
        query_camid_list = []
        for img_path in query_list:
            pid, camid = map(int, pattern.search(img_path).groups())
            query_pid_list.append(pid)
            query_camid_list.append(camid)

        gallery_camid_list = []
        gallery_pid_list = []
        for img_path in gallery_list:
            pid, camid = map(int, pattern.search(img_path).groups())
            gallery_pid_list.append(pid)
            gallery_camid_list.append(camid)

        model = get_model(traindata_num_classes, modelname)
        inference_samples(model, batch_size, query_list, gallery_list, query_pid_list, gallery_pid_list, query_camid_list, \
                          gallery_camid_list, adjust_rerank, npy_save_path)

if __name__ == "__main__":
    """
    需要更换的有 modelname, adjust_rerank 如果不是调整rerank 则需要给出需要融合的 npy_save_path
    """
    traindata_num_classes = 9968
    modelname = "/home/zxh/reid-baseline/new_experiment/r101_pcb_new/resnext101_ibn_a_model_70.pth"
    batch_size = 128*2
    merge_npy_paths = []
    adjust_rerank = True
    npy_save_path = ''
    main(traindata_num_classes, modelname, batch_size, merge_npy_paths, adjust_rerank, npy_save_path)

max_stars_repo_name: ZHUXUHAN/reid-baseline

max_issues_repo_name: ZHUXUHAN/reid-baseline

max_forks_repo_name: ZHUXUHAN/reid-baseline


text:
# Extend Gabor features version

import numpy as np
from skimage.filters import gabor_kernel
from scipy import ndimage

def GrayScaleNormalization(imgArray, imgMax,imgMin):

    imgRange = imgMax - imgMin
    imgArray = (imgArray - imgMin) * (255.0 / imgRange)

    # transfer to closest int
    imgArray = np.rint(imgArray).astype(np.int16)

    return imgArray

def genKernelBank(sigma_range, freq_range, out_kernel_bank):
    for i in range(4):
        theta = i / 4. * np.pi
        kernel_bank_per_theta = []

        for sigma in sigma_range:
            for freq in freq_range:
                kernel = np.real(gabor_kernel(freq, theta=theta, sigma_x=sigma, sigma_y=sigma))
                kernel_bank_per_theta.append(kernel)

        out_kernel_bank.append(kernel_bank_per_theta)
    return out_kernel_bank

def calcFeatures(dicomimg, xcoord, ycoord, width, height, kernel_bank,imgMax,imgMin):
    # xcoord and ycoord start from leftmost coordinate

    n_kernel_per_theta = len(kernel_bank[0])

    resultMean = np.zeros((4, n_kernel_per_theta))
    resultStd = np.zeros((4, n_kernel_per_theta))

    for theta in range(4):
        for i, kernel in enumerate(kernel_bank[theta]):

            x_ext_radius = (kernel.shape[0] + 1) / 2
            y_ext_radius = (kernel.shape[1] + 1) / 2

            subimageGabor = dicomimg[ycoord - y_ext_radius:(ycoord + height) + y_ext_radius, xcoord - x_ext_radius:(xcoord + width) + x_ext_radius]

            # print xcoord,ycoord
            # print y_ext_radius,x_ext_radius
            # gaborsize = np.shape(subimageGabor)
            # print 'gabor size:',gaborsize

            roi_in = GrayScaleNormalization(subimageGabor,imgMax,imgMin)

            roi_out = np.zeros(roi_in.shape)

            ndimage.filters.convolve(roi_in, kernel, output=roi_out, mode='constant', cval=0.0)

            zoom_roi_out = roi_out[y_ext_radius:y_ext_radius + height, x_ext_radius:x_ext_radius+ width]

            # zoominsize = np.shape(zoom_roi_out)
            # print 'zoom in size:',zoominsize

            resultMean[theta][i] = zoom_roi_out.mean()
            resultStd[theta][i] = zoom_roi_out.std()

    out_mean_vec = np.mean(resultMean, axis=0)
    out_std_vec = np.mean(resultStd, axis=0)
    gaborfeatures = np.column_stack((out_mean_vec,out_std_vec))

    return gaborfeatures

max_stars_repo_name: joshlyman/TextureAnalysis

reason: import numpy,from scipy