text:
import numpy as np

def layer_sizes_test_case():
    np.random.seed(1)
    X_assess = np.random.randn(5, 3)
    Y_assess = np.random.randn(2, 3)
    return X_assess, Y_assess

def initialize_parameters_test_case():
    n_x, n_h, n_y = 2, 4, 1
    return n_x, n_h, n_y

def forward_propagation_test_case():
    np.random.seed(1)
    X_assess = np.random.randn(2, 3)

    parameters = {'W1': np.array([[-0.00416758, -0.00056267],
        [-0.02136196,  0.01640271],
        [-0.01793436, -0.00841747],
        [ 0.00502881, -0.01245288]]),
     'W2': np.array([[-0.01057952, -0.00909008,  0.00551454,  0.02292208]]),
     'b1': np.array([[ 0.],
        [ 0.],
        [ 0.],
        [ 0.]]),
     'b2': np.array([[ 0.]])}

    return X_assess, parameters

def compute_cost_test_case():
    np.random.seed(1)
    Y_assess = np.random.randn(1, 3)
    parameters = {'W1': np.array([[-0.00416758, -0.00056267],
        [-0.02136196,  0.01640271],
        [-0.01793436, -0.00841747],
        [ 0.00502881, -0.01245288]]),
     'W2': np.array([[-0.01057952, -0.00909008,  0.00551454,  0.02292208]]),
     'b1': np.array([[ 0.],
        [ 0.],
        [ 0.],
        [ 0.]]),
     'b2': np.array([[ 0.]])}

    a2 = (np.array([[ 0.5002307 ,  0.49985831,  0.50023963]]))

    return a2, Y_assess, parameters

def backward_propagation_test_case():
    np.random.seed(1)
    X_assess = np.random.randn(2, 3)
    Y_assess = np.random.randn(1, 3)
    parameters = {'W1': np.array([[-0.00416758, -0.00056267],
        [-0.02136196,  0.01640271],
        [-0.01793436, -0.00841747],
        [ 0.00502881, -0.01245288]]),
     'W2': np.array([[-0.01057952, -0.00909008,  0.00551454,  0.02292208]]),
     'b1': np.array([[ 0.],
        [ 0.],
        [ 0.],
        [ 0.]]),
     'b2': np.array([[ 0.]])}

    cache = {'A1': np.array([[-0.00616578,  0.0020626 ,  0.00349619],
         [-0.05225116,  0.02725659, -0.02646251],
         [-0.02009721,  0.0036869 ,  0.02883756],
         [ 0.02152675, -0.01385234,  0.02599885]]),
     'A2': np.array([[ 0.5002307 ,  0.49985831,  0.50023963]]),
     'Z1': np.array([[-0.00616586,  0.0020626 ,  0.0034962 ],
         [-0.05229879,  0.02726335, -0.02646869],
         [-0.02009991,  0.00368692,  0.02884556],
         [ 0.02153007, -0.01385322,  0.02600471]]),
     'Z2': np.array([[ 0.00092281, -0.00056678,  0.00095853]])}

    return parameters, cache, X_assess, Y_assess

def update_parameters_test_case():
    parameters = {'W1': np.array([[-0.00615039,  0.0169021 ],
        [-0.02311792,  0.03137121],
        [-0.0169217 , -0.01752545],
        [ 0.00935436, -0.05018221]]),
     'W2': np.array([[-0.0104319 , -0.04019007,  0.01607211,  0.04440255]]),
     'b1': np.array([[ -8.97523455e-07],
        [  8.15562092e-06],
        [  6.04810633e-07],
        [ -2.54560700e-06]]),
     'b2': np.array([[  9.14954378e-05]])}

    grads = {'dW1': np.array([[ 0.00023322, -0.00205423],
        [ 0.00082222, -0.00700776],
        [-0.00031831,  0.0028636 ],
        [-0.00092857,  0.00809933]]),
     'dW2': np.array([[ -1.75740039e-05,   3.70231337e-03,  -1.25683095e-03,
          -2.55715317e-03]]),
     'db1': np.array([[  1.05570087e-07],
        [ -3.81814487e-06],
        [ -1.90155145e-07],
        [  5.46467802e-07]]),
     'db2': np.array([[ -1.08923140e-05]])}

    return parameters, grads

def nn_model_test_case():
    np.random.seed(1)
    X_assess = np.random.randn(2, 3)
    Y_assess = np.random.randn(1, 3)
    return X_assess, Y_assess

def predict_test_case():
    np.random.seed(1)
    X_assess = np.random.randn(2, 3)
    parameters = {'W1': np.array([[-0.00615039,  0.0169021 ],
        [-0.02311792,  0.03137121],
        [-0.0169217 , -0.01752545],
        [ 0.00935436, -0.05018221]]),
     'W2': np.array([[-0.0104319 , -0.04019007,  0.01607211,  0.04440255]]),
     'b1': np.array([[ -8.97523455e-07],
        [  8.15562092e-06],
        [  6.04810633e-07],
        [ -2.54560700e-06]]),
     'b2': np.array([[  9.14954378e-05]])}

    return parameters, X_assess

max_stars_repo_name: pedro-abundio-wang/deep-learning

max_issues_repo_name: pedro-abundio-wang/deep-learning


text:
import numpy as np

from abc import ABC
from abc import abstractmethod

from lagom.core.plotter import GridImage

try:  # workaround on server without fake screen but still running other things well
    from lagom.core.plotter import ImageViewer
except ImportError:
    import warnings
    warnings.warn('ImageViewer failed to import due to pyglet. ')

class VecEnv(ABC):
    r"""Base class for all asynchronous, vectorized environments.

    A vectorized environment handles batched data with its sub-environments. Each observation
    returned from vectorized environment is a batch of observations for each sub-environment.
    And :meth:`step` is expected to receive a batch of actions for each sub-environment.

    .. note::

        All sub-environments should share the identical observation and action spaces.
        In other words, a vector of multiple different environments is not supported.

    See docstring in :class:`Env` for more details about common APIs.

    The subclass should implement at least the following:

    - :meth:`step_async`
    - :meth:`step_wait`
    - :meth:`reset`
    - :meth:`get_images`
    - :meth:`close_extras`
    - :meth:`T`
    - :meth:`max_episode_reward`
    - :meth:`reward_range`

    """
    def __init__(self, list_make_env, observation_space, action_space):
        r"""Initialize the vectorized environment.

        .. note::

            The random seeds for all environments should be handled within each make_env function.
            And it should not handled here in general, because of APIs for parallelized environments.

        Args:
            list_make_env (list): a list of functions each returns an instantiated enviroment.
            observation_space (Space): observation space of the environment
            action_space (Space): action space of the environment
        """
        self.list_make_env = list_make_env
        self.num_env = len(self.list_make_env)

        self._observation_space = observation_space
        self._action_space = action_space

        # Some settings
        self.closed = False
        self.viewer = None  # rendering

    @abstractmethod
    def step_async(self, actions):
        r"""Ask all the environments to take a step with a list of actions, each for one environment.

        .. note::

            Call :meth:`step_wait` to get the step execution results.

        .. note::

            Do not call this function when it is already pending.

        Args:
            actions (list): a list of actions, each for one environment.
        """
        pass

    @abstractmethod
    def step_wait(self):
        r"""Wait for the jobs in :meth:`step_async` to finish and return all results.

        Returns
        -------
        observations : list
            a list of observations, each returned from one environment after executing the given action.
        rewards : list
            a list of scalar rewards, each returned from one environment.
        dones : list
            a list of booleans indicating whether the episode terminates, each returned from one environment.
        infos : list
            a list of dictionaries of additional informations, each returned from one environment.
        """
        pass

    def step(self, actions):
        r"""Take an action in each environment for one time step through the environments' dynamics.

        It firstly calls :meth:`step_async` to distribute actions to all environments and then
        calls :meth:`step_wait` to get all the results.

        See docstrings in :meth:`step_async` and :meth:`step_wait` for more details.
        """
        assert len(actions) == self.num_env, f'expected length {self.num_env}, got {len(actions)}'
        # Send actions to all environments asynchronously
        self.step_async(actions)

        # wait to receive all results and return them
        return self.step_wait()

    @abstractmethod
    def reset(self):
        r"""Reset all the environments and return a list of initial observations from each environment.

        .. warning::

            If :meth:`step_async` is still working, then it will be aborted.

        Returns
        -------
        observations : list
            a list of initial observations from all environments.
        """
        pass

    def render(self, mode='human'):
        r"""Render all the environments.

        It firstly retrieve RGB images from all environments and use :class:`GridImage`
        to make a grid of them as a single image. Then it either returns the image array
        or display the image to the screen by using :class:`ImageViewer`.

        See docstring in :class:`Env` for more detais about rendering.
        """
        # Get images from all environments with shape [N, H, W, C]
        imgs = self.get_images()
        imgs = np.stack(imgs)
        # Make a grid of images
        grid = GridImage(ncol=5, padding=5, pad_value=0)
        imgs = imgs.transpose(0, 3, 1, 2)  # to shape [N, C, H, W]
        grid.add(imgs)
        gridimg = np.asarray(grid())
        gridimg = gridimg.transpose(0, 2, 3, 1)  # back to shape [N, H, W, C]

        # render the grid of image
        if mode == 'human':
            self.get_viewer()(gridimg)
        elif mode == 'rgb_array':
            return gridimg
        else:
            raise ValueError(f'expected human or rgb_array, got {mode}')

    @abstractmethod
    def get_images(self):
        r"""Returns a batched RGB array with shape [N, H, W, C] from all environments.

        Returns
        -------
        imgs : ndarray
            a batched RGB array with shape [N, H, W, C]
        """
        pass

    def get_viewer(self):
        r"""Returns an instantiated :class:`ImageViewer`.

        Returns
        -------
        viewer : ImageViewer
            an image viewer
        """
        if self.viewer is None:  # create viewer is not existed
            self.viewer = ImageViewer(max_width=500)  # set a max width here

        return self.viewer

    @abstractmethod
    def close_extras(self):
        r"""Clean up the extra resources e.g. beyond what's in this base class. """
        pass

    def close(self):
        r"""Close all environments.

        It closes all the existing image viewers, then calls :meth:`close_extras` and set
        :attr:`closed` as ``True``.

        .. warning::

            This function itself does not close the environments, it should be handled
            in :meth:`close_extras`. This is useful for parallelized environments.

        .. note::

            This will be automatically called when garbage collected or program exited.

        """
        if self.closed:  # already closed
            return

        # Close viewer if existed
        if self.viewer is not None:
            self.viewer.close()

        # Close extra
        self.close_extras()

        # Set closed flag
        self.closed = True

    @property
    def unwrapped(self):
        r"""Unwrap this vectorized environment.

        Useful for sequential wrappers applied, it can access information from the original
        vectorized environment.
        """
        if isinstance(self, VecEnvWrapper):
            return self.venv.unwrapped
        else:
            return self

    @property
    def observation_space(self):
        r"""Returns a :class:`Space` object to define the observation space. """
        return self._observation_space

    @property
    def action_space(self):
        r"""Returns a :class:`Space` object to define the action space. """
        return self._action_space

    @property
    @abstractmethod
    def T(self):
        r"""Maximum horizon of the environment, if available. """
        pass

    @property
    @abstractmethod
    def max_episode_reward(self):
        r"""Maximum episodic rewards of the environment, if available. """
        pass

    @property
    @abstractmethod
    def reward_range(self):
        r"""Returns a tuple of min and max possible rewards.

        .. note::

            By default, it could be infinite range i.e. ``(-float('inf'), float('inf'))``.
            One can also set a narrower range.

        """
        pass

    def __repr__(self):
        return f'<{self.__class__.__name__}: {self.list_make_env[0].keywords["env_id"]}, n: {self.num_env}>'

class VecEnvWrapper(VecEnv):
    r"""Wraps the vectorized environment to allow a modular transformation.

    This class is the base class for all wrappers for vectorized environments. The subclass
    could override some methods to change the behavior of the original vectorized environment
    without touching the original code.

    .. note::

        Don't forget to call ``super().__init__(venv)`` if the subclass overrides :meth:`__init__`.

    The subclass should implement at least the following:

    - :meth:`step_wait`
    - :meth:`reset`
    - :meth:`close_extras`

    """
    def __init__(self, venv):
        assert isinstance(venv, VecEnv), f'expected VecEnv type, got {type(venv)}'

        self.venv = venv

        # Call VecEnv constructor to make some settings as direct attributes in this wrapper class
        super().__init__(list_make_env=venv.list_make_env,
                         observation_space=venv.observation_space,
                         action_space=venv.action_space)

    def step_async(self, actions):
        self.venv.step_async(actions)

    @abstractmethod
    def step_wait(self):
        pass

    @abstractmethod
    def reset(self):
        pass

    def get_images(self):
        self.venv.get_images()

    @abstractmethod
    def close_extras(self):
        pass

    @property
    def T(self):
        return self.venv.T

    @property
    def max_episode_reward(self):
        return self.venv.max_episode_reward

    @property
    def reward_range(self):
        return self.venv.reward_range

ext: py

max_stars_repo_name: dkorduban/lagom

max_issues_repo_name: dkorduban/lagom

max_forks_repo_name: dkorduban/lagom

reason: import numpy


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Apr 26 08:48:00 2018

@author: prmiles
"""

from pymcmcstat.structures.ParameterSet import ParameterSet
from pymcmcstat.samplers.Metropolis import Metropolis
import test.general_functions as gf

import unittest
from mock import patch
import numpy as np

# --------------------------
# Evaluation
# --------------------------
class UnpackSet(unittest.TestCase):

    def test_unpack_set(self):
        CL = {'theta':1.0, 'ss': 1.0, 'prior':0.0, 'sigma2': 0.0}
        parset = ParameterSet(theta = CL['theta'], ss = CL['ss'], prior = CL['prior'], sigma2 = CL['sigma2'])

        MA = Metropolis()
        oldpar, ss, oldprior, sigma2 = MA.unpack_set(parset)
        NL = {'theta':oldpar, 'ss': ss, 'prior':oldprior, 'sigma2': sigma2}
        self.assertDictEqual(CL,NL)

# --------------------------
class CalculatePosteriorRatio(unittest.TestCase):
    @classmethod
    def setup_size_2(cls):
        MA = Metropolis()
        ss1 = np.array([1., 2.])
        ss2 = np.array([1.1, 2.4])
        newprior = np.array([0.,0.])
        oldprior = np.array([0.,0.])
        sigma2 = np.array([1.,1.])
        return MA, ss1, ss2, newprior, oldprior, sigma2

    def test_size_of_alpha_for_1d_nsos(self):
        MA = Metropolis()
        ss1 = np.array([1.])
        ss2 = np.array([1.1])
        newprior = np.array([0.])
        oldprior = np.array([0.])
        sigma2 = np.array([1.])
        alpha = MA.calculate_posterior_ratio(ss1, ss2, sigma2, newprior, oldprior)
        self.assertEqual(alpha.size, 1)

    def test_size_of_alpha_for_2d_nsos(self):
        MA, ss1, ss2, newprior, oldprior, sigma2 = self.setup_size_2()
        alpha = MA.calculate_posterior_ratio(ss1, ss2, sigma2, newprior, oldprior)
        self.assertEqual(alpha.size, 1)

    def test_likelihood_goes_up(self):
        MA, ss1, ss2, newprior, oldprior, sigma2 = self.setup_size_2()
        alpha1 = MA.calculate_posterior_ratio(ss1, ss2, sigma2, newprior, oldprior)

        ss3 = np.array([0.4, 0.5])
        alpha2 = MA.calculate_posterior_ratio(ss3, ss2, sigma2, newprior, oldprior)
        self.assertTrue(alpha1 < alpha2)

    def test_likelihood_goes_down(self):
        MA, ss1, ss2, newprior, oldprior, sigma2 = self.setup_size_2()
        alpha1 = MA.calculate_posterior_ratio(ss1, ss2, sigma2, newprior, oldprior)

        ss3 = np.array([1.4, 2.5])
        alpha2 = MA.calculate_posterior_ratio(ss3, ss2, sigma2, newprior, oldprior)
        self.assertTrue(alpha2 < alpha1)

    def test_alpha_value(self):
        MA, ss1, ss2, newprior, oldprior, sigma2 = self.setup_size_2()
        alpha1 = MA.calculate_posterior_ratio(ss1, ss2, sigma2, newprior, oldprior)
        self.assertTrue(np.allclose(alpha1, np.array([2.568050833375483])), msg = str('alpha = {}'.format(alpha1)))

# --------------------------
def setup_CL(theta=1.0, ss=1.0, prior=0.0, sigma2=1.0):
    return {'theta':theta, 'ss': np.array([ss]), 'prior': prior, 'sigma2': sigma2}

class RunMetropolisStep(unittest.TestCase):
    @classmethod
    def setup_rms(cls, CL):
        sos_object, prior_object, parameters = gf.setup_mcmc_case_mh()
        MS = Metropolis()
        R = np.array([[0.4, 0.2],[0, 0.3]])
        parset = ParameterSet(theta=CL['theta'], ss=CL['ss'],
                              prior=CL['prior'], sigma2=CL['sigma2'])
        accept, _, outbound, npar_sample_from_normal = MS.run_metropolis_step(
                old_set=parset, parameters=parameters, R=R,
                prior_object=prior_object, sos_object=sos_object)
        return accept, outbound

    @patch('pymcmcstat.samplers.Metropolis.is_sample_outside_bounds', return_value = True)
    def test_run_step_outside_bounds(self, mock_1):
        accept, outbound = self.setup_rms(setup_CL(sigma2 = 1.0))
        self.assertEqual(outbound, 1, msg='outbound set to 1')
        self.assertEqual(accept, 0, msg='Not accepted because outside bounds')

    @patch('pymcmcstat.samplers.Metropolis.is_sample_outside_bounds', return_value = False)
    @patch('pymcmcstat.samplers.Metropolis.calculate_log_posterior_ratio',
           return_value=np.log(1.1))
    def test_run_step_inside_bounds(self, mock_1, mock_2):
        accept, outbound = self.setup_rms(setup_CL())
        self.assertEqual(outbound, 0, msg='outbound set to 0')
        self.assertTrue(accept, msg='Accepted because loglikelihood > 0')

    @patch('pymcmcstat.samplers.Metropolis.is_sample_outside_bounds',
           return_value=False)
    @patch('pymcmcstat.samplers.Metropolis.calculate_log_posterior_ratio',
           return_value=np.log(0.5))
    @patch('numpy.random.rand',
           return_value=0.4)
    def test_run_step_inside_bounds_test_accept(self, mock_1, mock_2, mock_3):
        accept, outbound = self.setup_rms(setup_CL())
        self.assertEqual(outbound, 0, msg='outbound set to 0')
        self.assertTrue(accept, msg='Accepted because 0.5 > 0.4')

    @patch('pymcmcstat.samplers.Metropolis.is_sample_outside_bounds', return_value = False)
    @patch('pymcmcstat.samplers.Metropolis.Metropolis.calculate_posterior_ratio', return_value = 0.3)
    @patch('numpy.random.rand', return_value = 0.4)
    def test_run_step_inside_bounds_test_accept_fail(self, mock_1, mock_2, mock_3):
        accept, outbound = self.setup_rms(setup_CL())
        self.assertEqual(outbound, 0, msg='outbound set to 0')
        self.assertEqual(accept, 0, msg='Accepted because 0.3 < 0.4')

max_stars_repo_name: vishalbelsare/pymcmcstat

max_issues_repo_name: vishalbelsare/pymcmcstat

max_forks_repo_name: vishalbelsare/pymcmcstat


text:
# sympy/galgebra/vector.py

"""
vector.py is a helper class for the MV class that defines the basis
vectors and metric and calulates derivatives of the basis vectors for
the MV class.
"""

import itertools
import copy

from sympy import Symbol, S, Matrix, trigsimp, diff, expand
from sympy.core.compatibility import range

from sympy.galgebra.printing import GA_Printer
from sympy.galgebra.stringarrays import str_array
from sympy.galgebra.ncutil import linear_derivation, bilinear_product
from sympy.galgebra.debug import oprint

def flatten(lst):
    return list(itertools.chain(*lst))

def TrigSimp(x):
    return trigsimp(x, recursive=True)

class Vector(object):
    """
    Vector class.

    Setup is done by defining a set of basis vectors in static function 'Bases'.
    The linear combination of scalar (commutative) sympy quatities and the
    basis vectors form the vector space.  If the number of basis vectors
    is 'n' the metric tensor is formed as an n by n sympy matrix of scalar
    symbols and represents the dot products of pairs of basis vectors.
    """

    is_orthogonal = False

    @staticmethod
    def setup(base, n=None, metric=None, coords=None, curv=(None, None), debug=False):
        """
        Generate basis of vector space as tuple of vectors and
        associated metric tensor as Matrix.  See str_array(base,n) for
        usage of base and n and str_array(metric) for usage of metric.

        To overide elements in the default metric use the character '#'
        in the metric string.  For example if one wishes the diagonal
        elements of the metric tensor to be zero enter metric = '0 #,# 0'.

        If the basis vectors are e1 and e2 then the default metric -

            Vector.metric = ((dot(e1,e1),dot(e1,e2)),dot(e2,e1),dot(e2,e2))

        becomes -

            Vector.metric = ((0,dot(e1,e2)),(dot(e2,e1),0)).

        The function dot returns a Symbol and is symmetric.

        The functions 'Bases' calculates the global quantities: -

            Vector.basis
                tuple of basis vectors
            Vector.base_to_index
                dictionary to convert base to base inded
            Vector.metric
                metric tensor represented as a matrix of symbols and numbers

        """
        Vector.is_orthogonal = False
        Vector.coords = coords
        Vector.subscripts = []
        base_name_lst = base.split(' ')

        # Define basis vectors

        if '*' in base:
            base_lst = base.split('*')
            base = base_lst[0]
            Vector.subscripts = base_lst[1].split('|')
            base_name_lst = []
            for subscript in Vector.subscripts:
                base_name_lst.append(base + '_' + subscript)
        else:
            if len(base_name_lst) > 1:
                Vector.subscripts = []
                for base_name in base_name_lst:
                    tmp = base_name.split('_')
                    Vector.subscripts.append(tmp[-1])
            elif len(base_name_lst) == 1 and Vector.coords is not None:
                base_name_lst = []
                for coord in Vector.coords:
                    Vector.subscripts.append(str(coord))
                    base_name_lst.append(base + '_' + str(coord))
            else:
                raise TypeError("'%s' does not define basis vectors" % base)

        basis = []
        base_to_index = {}
        index = 0
        for base_name in base_name_lst:
            basis_vec = Vector(base_name)
            basis.append(basis_vec)
            base_to_index[basis_vec.obj] = index
            index += 1

        Vector.base_to_index = base_to_index
        Vector.basis = tuple(basis)

        # define metric tensor

        default_metric = []
        for bv1 in Vector.basis:
            row = []
            for bv2 in Vector.basis:
                row.append(Vector.basic_dot(bv1, bv2))
            default_metric.append(row)
        Vector.metric = Matrix(default_metric)
        if metric is not None:
            if metric[0] == '[' and metric[-1] == ']':
                Vector.is_orthogonal = True
                metric_str_lst = metric[1:-1].split(',')
                Vector.metric = []
                for g_ii in metric_str_lst:
                    Vector.metric.append(S(g_ii))
                Vector.metric = Matrix(Vector.metric)
            else:
                metric_str_lst = flatten(str_array(metric))
                for index in range(len(metric_str_lst)):
                    if metric_str_lst[index] != '#':
                        Vector.metric[index] = S(metric_str_lst[index])

        Vector.metric_dict = {}  # Used to calculate dot product
        N = range(len(Vector.basis))
        if Vector.is_orthogonal:
            for ii in N:
                    Vector.metric_dict[Vector.basis[ii].obj] = Vector.metric[ii]
        else:
            for irow in N:
                for icol in N:
                    Vector.metric_dict[(Vector.basis[irow].obj, Vector.basis[icol].obj)] = Vector.metric[irow, icol]

        # calculate tangent vectors and metric for curvilinear basis

        if curv != (None, None):
            X = S.Zero
            for (coef, base) in zip(curv[0], Vector.basis):
                X += coef * base.obj
            Vector.tangents = []
            for (coord, norm) in zip(Vector.coords, curv[1]):
                tau = diff(X, coord)
                tau = trigsimp(tau)
                tau /= norm
                tau = expand(tau)
                Vtau = Vector()
                Vtau.obj = tau
                Vector.tangents.append(Vtau)
            metric = []
            for tv1 in Vector.tangents:
                row = []
                for tv2 in Vector.tangents:
                    row.append(tv1 * tv2)
                metric.append(row)
            metric = Matrix(metric)
            metric = metric.applyfunc(TrigSimp)
            Vector.metric_dict = {}
            if metric.is_diagonal:
                Vector.is_orthogonal = True
                tmp_metric = []
                for ii in N:
                    tmp_metric.append(metric[ii, ii])
                    Vector.metric_dict[Vector.basis[ii].obj] = metric[ii, ii]
                Vector.metric = Matrix(tmp_metric)
            else:
                Vector.is_orthogonal = False
                Vector.metric = metric
                for irow in N:
                    for icol in N:
                        Vector.metric_dict[(Vector.basis[irow].obj, Vector.basis[icol].obj)] = Vector.metric[irow, icol]
            Vector.norm = curv[1]

            if debug:
                oprint('Tangent Vectors', Vector.tangents,
                       'Metric', Vector.metric,
                       'Metric Dictionary', Vector.metric_dict,
                       'Normalization', Vector.norm, dict_mode=True)

            # calculate derivatives of tangent vectors

            Vector.dtau_dict = None
            dtau_dict = {}

            for x in Vector.coords:
                for (tau, base) in zip(Vector.tangents, Vector.basis):
                    dtau = tau.diff(x).applyfunc(TrigSimp)
                    result = S.Zero
                    for (t, b) in zip(Vector.tangents, Vector.basis):
                        t_dtau = TrigSimp(t * dtau)
                        result += t_dtau * b.obj
                    dtau_dict[(base.obj, x)] = result

            Vector.dtau_dict = dtau_dict

            if debug:
                oprint('Basis Derivatives', Vector.dtau_dict, dict_mode=True)

        return tuple(Vector.basis)

    def __init__(self, basis_str=None):
        if isinstance(basis_str, Vector):
            self.obj = basis_str
        else:
            if basis_str is None or basis_str == '0':
                self.obj = S(0)
            else:
                self.obj = Symbol(basis_str, commutative=False)

    """
    def diff(self, x):
        (coefs, bases) = linear_expand(self.obj)
        result = S.Zero
        for (coef, base) in zip(coefs, bases):
            result += diff(coef, x) * base
        return result
    """

    def diff(self, x):
        Dself = Vector()
        if isinstance(Vector.dtau_dict, dict):
            Dself.obj = linear_derivation(self.obj, Vector.Diff, x)
        else:
            Dself.obj = diff(self.obj, x)
        return Dself

    @staticmethod
    def basic_dot(v1, v2):
        """
        Dot product of two basis vectors returns a Symbol
        """
        i1 = list(Vector.basis).index(v1)  # Python 2.5
        i2 = list(Vector.basis).index(v2)  # Python 2.5
        if i1 < i2:
            dot_str = '(' + str(Vector.basis[i1]) + '.' + str(Vector.basis[i2]) + ')'
        else:
            dot_str = '(' + str(Vector.basis[i2]) + '.' + str(Vector.basis[i1]) + ')'
        return Symbol(dot_str)

    @staticmethod
    def dot(b1, b2):
        if Vector.is_orthogonal:
            if b1 != b2:
                return S.Zero
            else:
                return Vector.metric_dict[b1]
        else:
            return Vector.metric_dict[(b1, b2)]

    @staticmethod
    def Diff(b, x):
        return Vector.dtau_dict[(b, x)]

    ######################## Operator Definitions#######################

    def __str__(self):
        return GA_Printer().doprint(self)

    def __mul__(self, v):
        if not isinstance(v, Vector):
            self_x_v = Vector()
            self_x_v.obj = self.obj * v
            return self_x_v
        else:
            result = expand(self.obj * v.obj)
            result = bilinear_product(result, Vector.dot)
            return result

    def __rmul__(self, s):
        s_x_self = Vector()
        s_x_self.obj = s * self.obj
        return s_x_self

    def __add__(self, v):
        self_p_v = Vector()
        self_p_v.obj = self.obj + v.obj
        return self_p_v

    def __add_ab__(self, v):
        self.obj += v.obj
        return

    def __sub__(self, v):
        self_m_v = Vector()
        self_m_v.obj = self.obj - v.obj
        return self_m_v

    def __sub_ab__(self, v):
        self.obj -= v.obj
        return

    def __pos__(self):
        return self

    def __neg__(self):
        n_self = copy.deepcopy(self)
        n_self.obj = -self.obj
        return n_self

    def applyfunc(self, fct):
        fct_self = Vector()
        fct_self.obj = fct(self.obj)
        return fct_self

max_stars_repo_name: toolforger/sympy

max_issues_repo_name: toolforger/sympy

max_forks_repo_name: toolforger/sympy

reason: from sympy


text:
"""Transform module tests"""

import unittest

import numpy as np

from vibromaf.signal.perception import PerceptualThreshold
from vibromaf.signal.transform import (
    BlockBuilder,
    PerceptualSpectrumBuilder,
    compute_block_dct,
    compute_block_dft,
    cut_off_strategy,
    preprocess_input_signal,
    zero_padding_strategy,
)

class TransformTest(unittest.TestCase):
    """Transform Test"""

    def test_block_dft__ones_block__constant_block(self):
        input_block = np.array([1, 0, 0, 0, 0, 0, 0, 0, 0, 0])
        result = compute_block_dft(input_block)
        self.assertListEqual([-10] * 10, list(result))

    def test_block_dft__block_with_legth_10__block_with_same_length(self):
        input_block = np.linspace(1, 1, 10)
        result = compute_block_dft(input_block)
        self.assertEqual(10, np.size(result))

    def test_block_dft__parseval_theorem(self):
        input_block = np.linspace(1, 1, 16)
        spectrum = compute_block_dft(input_block)
        spectrum = np.power(10, spectrum / 20)
        result = np.sum(np.power(spectrum, 2))
        self.assertAlmostEqual(24, float(result), delta=0.05)

    def test_block_dct__block_of_length_10__block_of_length_10(self):
        input_block = np.linspace(1, 1, 10)
        result = compute_block_dct(input_block)
        self.assertEqual(10, np.size(result))

    def test_block_dct__ones_block_of_length_10__block_with_2_and_zeros(self):
        input_block = np.ones(10)
        result = compute_block_dct(input_block)
        self.assertListEqual([-np.inf] * 9, list(result[1:]))
        self.assertAlmostEqual(10, result[0])

    def test_cut_off_strategy__signal_more_than_block_length__cut_off(self):
        input_signal = np.array([1, 2, 3, 4, 5])
        result = cut_off_strategy(input_signal, 3)
        self.assertListEqual([1, 2, 3], list(result))

    def test_cut_off_strategy__signal_more_than_multiple_block_length__cut_off(self):
        input_signal = np.array([1, 2, 3, 4, 5])
        result = cut_off_strategy(input_signal, 2)
        self.assertListEqual([1, 2, 3, 4], list(result))

    def test_cut_off_strategy__signal_smaller_than_block_length__empty(self):
        input_signal = np.array([1, 2])
        result = cut_off_strategy(input_signal, 3)
        self.assertListEqual([], list(result))

    def test_cut_off_strategy__signal_equal_to_block_length__unchanged(self):
        input_signal = np.array([1, 2, 3])
        result = cut_off_strategy(input_signal, 3)
        self.assertListEqual([1, 2, 3], list(result))

    def test_zero_padding_strategy__signal_more_than_block_length__fill_with_zeros(
        self,
    ):
        input_signal = np.array([1, 2, 3, 4])
        result = zero_padding_strategy(input_signal, 3)
        self.assertListEqual([1, 2, 3, 4, 0, 0], list(result))

    def test_zero_padding_strategy__signal_smaller_than_block_length__zeros_only(
        self,
    ):
        input_signal = np.array([1, 2])
        result = zero_padding_strategy(input_signal, 3)
        self.assertListEqual([1, 2, 0], list(result))

    def test_zero_padding_strategy__signal_equal_to_block_length__unchanged(
        self,
    ):
        input_signal = np.array([1, 2, 3])
        result = zero_padding_strategy(input_signal, 3)
        self.assertListEqual([1, 2, 3], list(result))

    def test_preprocess_input_signal__dist_larger_than_ref__dist_should_be_truncated_and_warn(
        self,
    ):
        signal = np.array([0, 1])
        dist = np.array([0, 1, 2])
        with self.assertWarnsRegex(
            RuntimeWarning,
            r"Truncating distorted signal .* since longer than reference",
        ):
            result = preprocess_input_signal(dist, signal)
        self.assertListEqual(list(signal), list(result))

    def test_preprocess_input_signal__dist_shorter_than_ref__should_throw(self):
        signal = np.array([0, 1, 2])
        dist = np.array([0, 1])
        with self.assertRaisesRegex(ValueError, r"Distorted .* must not be shorter"):
            preprocess_input_signal(dist, signal)

class BlockBuilderTest(unittest.TestCase):
    """Block Builder Test"""

    def test_divide__empty_array__should_throw(self):
        unit = BlockBuilder(10)
        with self.assertRaises(ValueError):
            unit.divide(np.array([]))

    def test_divide__array_shorter_than_block_length__should_throw(self):
        unit = BlockBuilder(10)
        with self.assertRaises(ValueError):
            unit.divide(np.array([1, 2, 3]))

    def test_divide__perfect_matching_array__perfect_division(self):
        input_signal = np.array([1, 2, 3, 4, 5, 6])
        unit = BlockBuilder(2)
        result = unit.divide(input_signal)
        self.assertListEqual([1, 2], list(result[0]))
        self.assertListEqual([3, 4], list(result[1]))
        self.assertListEqual([5, 6], list(result[2]))

    def test_divide__perfect_not_matching_array__division_truncated(self):
        input_signal = np.array([1, 2, 3, 4, 5])
        unit = BlockBuilder(2)
        result = unit.divide(input_signal)
        self.assertListEqual([1, 2], list(result[0]))
        self.assertListEqual([3, 4], list(result[1]))
        self.assertEqual((2, 2), result.shape)

    def test_divide_and_normalize__periodic_array__array_divided_by_two(self):
        input_signal = np.array([-2, 2, -2, 2])
        unit = BlockBuilder(4)
        result = unit.divide_and_normalize(input_signal)
        self.assertListEqual([-1, 1, -1, 1], list(result[0]))

    def test_divide_and_normalize__periodic_array__array_divided_by_two_and_shifted(
        self,
    ):
        input_signal = np.array([0, 4, 0, 4])
        unit = BlockBuilder(4)
        result = unit.divide_and_normalize(input_signal)
        self.assertListEqual([-1, 1, -1, 1], list(result[0]))

    def test_divide_and_normalize__multiple_blocks__correct_reshaped(self):
        input_signal = np.array([-2, 2, -2, 2, -2, 2, -2, 2])
        unit = BlockBuilder(2)
        result = unit.divide_and_normalize(input_signal)
        self.assertTrue(
            np.array_equal(np.array([[-1, 1], [-1, 1], [-1, 1], [-1, 1]]), result)
        )

    def test_divide_and_normalize__input_signals_identical_zero__should_be_not_nan(
        self,
    ):
        input_signal = np.zeros(10)
        unit = BlockBuilder(5)
        result = unit.divide_and_normalize(input_signal)
        self.assertFalse(np.isnan(result).any())

class PerceptualSpectrumBuilderTest(unittest.TestCase):
    """Perceptual Spectrum Builder Test"""

    def test_compute_perceptual_spectrum__signal_length_equals_block_length__signal_shape_unchanged(
        self,
    ):
        input_signal = np.array([1, 2, 3, 4])
        unit = PerceptualSpectrumBuilder(BlockBuilder(4))
        result = unit.compute_perceptual_spectrum(input_signal)
        self.assertEqual(input_signal.size, result.size)
        self.assertEqual(result.shape, (1, 4))

    def test_compute_perceptual_spectrum__signal_length_equals_twice_block_length__signal_shape_changed(
        self,
    ):
        input_signal = np.array([1, 2, 3, 4])
        unit = PerceptualSpectrumBuilder(BlockBuilder(2))
        result = unit.compute_perceptual_spectrum(input_signal)
        self.assertEqual(input_signal.size, result.size)
        self.assertEqual(result.shape, (2, 2))

    def test_compute_perceptual_spectrum__signal_with_ones__block_wise_dct(self):
        input_signal = np.array([1, 1])
        threshold = PerceptualThreshold(0)

        unit = PerceptualSpectrumBuilder(
            block_builder=BlockBuilder(2),
            block_transform_strategy=compute_block_dct,
            perceptual_threshold=threshold,
        )
        result = unit.compute_perceptual_spectrum(input_signal)

        result += threshold.calculate(2)
        self.assertAlmostEqual(3.01, result[0, 0], delta=0.001)  # sqrt(2) in dB
        self.assertEqual(-np.inf, result[0, 1])

max_stars_repo_name: hofbi/vibromaf

max_issues_repo_name: hofbi/vibromaf

max_forks_repo_name: hofbi/vibromaf


text:
# -*- coding: utf-8 -*-

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

__author__          = "Vanhee Laurent (and Marée Raphaël) <laurent.vanhee@student.uliege.ac.be>"
__copyright__       = "Copyright 2010-2017 University of Liège, Belgium, http://www.cytomine.be/"

import csv
import datetime
import math
import os
import socket
import time
import numpy as np
from cytomine.models import *
from shapely.wkt import loads
from shutil import copyfile
import config
from cytomine import Cytomine
from PIL import Image
import sys
import inspect
##TODO : add ConnectionHistory to python client and use it

def get_data(project_dir, id_project, users_metadata_file, id_ref_user, im_subset=None, us_subset=None,
             cytomine_host=config.cytomine_host, cytomine_public_key=config.cytomine_public_key,
             cytomine_private_key=config.cytomine_private_key, modules=None):
    """

    :param project_dir: gold, silver
    :param id_project: 2338 -> silver, 1197608-> gold
    :param users_metadata_file: stats/students_gold.csv, stats/students_silver.csv
    :param id_ref_user : 1590 -> gold, 1611 -> silver
    :param im_subset : optional, makes script only fetch info on a subset of images
    :param us_subset : optional, makes script only fetch info on a subset of users
    :param cytomine_host: optional, Cytomine host address
    :param cytomine_private_key: optional, Cytomine user private key
    :param cytomine_public_key: optional, Cytomine user public key
    :param modules: optional, file directory containing modules
    :return:
    """

    #Connection to Cytomine Core
    conn = Cytomine(cytomine_host, cytomine_public_key, cytomine_private_key, base_path = '/api/', working_path = '/tmp/', verbose= False)

    rescaled_size = 1024
    working_path = config.WORKING_DIRECTORY  # directory should exist

    maxperpage = 500  # unamur
    opening_delay = 10000  # 10s

    timestep = 86400  # 1 day (in seconds)

    id_user = 27389949  # DEMO-LANDMARK-ZEBRAFISH
    start_time = config.start_time
    end_time = config.end_time

    if '/' not in project_dir:
        project_dir = project_dir + "/"
    # --------------------------------------------------------------------------------------
    # --------------------------------------------------------------------------------------
    if not os.path.exists(working_path + project_dir):
        os.makedirs(working_path + project_dir)

    # We convert start/end string time to python timestamps (multiply by 1000 as python expects seconds)
    start_timestamp = long(1000 * time.mktime(datetime.datetime.strptime(start_time, "%Y-%m-%d %H:%M:%S").timetuple()))
    end_timestamp = long(1000 * time.mktime(datetime.datetime.strptime(end_time, "%Y-%m-%d %H:%M:%S").timetuple()))

    # Create project stats file
    csv_stats_filename = 'stats.csv'
    stats_file = os.path.join(working_path + project_dir, csv_stats_filename)
    fstats = open(stats_file, "wb")
    csvoutstats = csv.writer(fstats)
    csvoutstats.writerow([start_timestamp, end_timestamp])
    csvoutstats.writerow(
        ['id_project', 'id_image', 'id_user', 'username', 'email', 'nb_user_annotations', 'central_positions_nb_opens',
         'nb_positions', 'zoom_frequencies', 'nb_reference_points'])

    # Get all project users:
    id_users = conn.get_project_users(id_project)  # of from arglist

    # Get all image instances from project
    # Here we should check arg, if null then get all images from project
    image_instances = ImageInstanceCollection()
    image_instances.project = id_project
    image_instances = conn.fetch(image_instances)
    images = image_instances.data()
    print "Nb images in project: %d" % len(images)

    if not os.path.exists(working_path + project_dir + "images/"):
        os.makedirs(working_path + project_dir + "images/")

    # Create user id file
    if modules is not None:
        csv_user_filename = 'timeline.csv'
        modules_file = os.path.join(working_path + project_dir, csv_user_filename)
        copyfile(modules, modules_file)

    # Create user id file
    csv_user_filename = 'users.csv'
    users_file = os.path.join(working_path + project_dir, csv_user_filename)
    copyfile(users_metadata_file, users_file)

    fusers = open(users_metadata_file, "rb")
    csvoutusers = csv.reader(fusers)
    data_users = list(csvoutusers)
    data_users.pop(0)
    header = data_users.pop(0)
    user_idx = 0
    for i in range(len(header)):
        if header[i] == "CYTOMINE ID":
            user_idx = i
    userlist = []
    print user_idx
    for u_tmp in data_users:
        userlist.append(int(u_tmp[user_idx]))
    fusers.close()

    if im_subset is None:
        im_subset = []
        for image in images:
            im_subset.append(image.id)

    #Go through all images
    for image in images:
        id_image=image.id
        print "Downloading data associated to image %d" %id_image
        if id_image in im_subset:
            image_dir = "images/image_" + str(id_image)
            if not os.path.exists(working_path + project_dir + image_dir):
                os.makedirs(working_path + project_dir + image_dir)

            #get abstractimage thumb, compute rescaling factor (original image size / thumb size)
            image_instance = conn.get_image_instance(id_image)
            max_dim = max(image_instance.height,image_instance.width)
            image_depth = image_instance.depth
            if max_dim>rescaled_size:
                rescale_factor = max_dim/rescaled_size
            else:
                rescale_factor=1

            url = image_instance.preview[0:image_instance.preview.index('?')] + "?maxSize=" + str(rescaled_size)
            filename = working_path + project_dir + image_dir + "/image.png"
            if not os.path.exists(filename):
                conn.fetch_url_into_file(url, filename, override=True)

            rescaled_im = Image.open(filename)
            rescaled_width, rescaled_height = rescaled_im.size

            # Get reference annotations in this image (to detect if reference regions were visualized by user)
            nb_ref_annotations = 0
            if id_ref_user:
                ref_annotations = conn.get_annotations(id_image=id_image,
                                                       id_user=id_ref_user,
                                                       id_project=id_project,
                                                       showWKT=True
                                                       )
                nb_ref_annotations = len(ref_annotations.data())
            else:
                ref_annotations = None
            #save the center of the reference annotations in a csv file
            csv_filename = 'reference_cytomine_annotations.csv'
            if nb_ref_annotations > 0:

                output_annotation_file = os.path.join(working_path + project_dir + image_dir + "/", csv_filename)
                f = open(output_annotation_file, "wb")
                csv_annotations = csv.writer(f)
                csv_annotations.writerow(['type', 'x_center', 'y_center', 'annotationIdent', 'localIdent'])
                for a in ref_annotations.data():

                    descr = conn.get_annotation_properties(a.id)
                    l_id = 0
                    for prop in descr.data():
                        if prop.key == 'n':
                            l_id = prop.value

                    geom = loads(a.location)

                    if geom.type == 'Point':
                        csv_annotations.writerow([geom.type, round(geom.x / rescale_factor), round(geom.y / rescale_factor), a.id, l_id])
                    else:
                        csv_annotations.writerow([geom.type, round(geom.centroid.x / rescale_factor), round(geom.centroid.y / rescale_factor), a.id, l_id])
                f.close()

            # for this image, go through project's users (except those not in provided userlist)
            for u in id_users.data():
                if u.id in userlist and (us_subset is None or u.id in us_subset ):

                    zooms = np.zeros(image_instance.depth)
                    id_user = u.id
                    if not os.path.exists(working_path + project_dir + image_dir + "/user_positions"):
                        os.makedirs(working_path + project_dir +  image_dir + "/user_positions")
                    # Get_positions for this user in this image: using paging (maxperpage) and using start/end timestamps
                    pos_success = False
                    while (not pos_success):
                        # Retry if we got error
                        try:
                            positions = conn.get_positions(id_image=id_image,
                                                           id_user=id_user,
                                                           maxperpage=maxperpage,
                                                           afterthan=start_timestamp,
                                                           beforethan=end_timestamp,
                                                           showDetails=True)
                            pos_success = True
                        except socket.error:
                            print socket.error
                            time.sleep(1)
                            continue
                        except socket.timeout:
                            print socket.timeout
                            time.sleep(1)
                            continue
                        except ValueError:
                            print socket.timeout
                            time.sleep(1)
                            continue
                        except Exception:
                            time.sleep(1)
                            continue

                    nb_positions = len(positions.data())
                    nb_opens = 0
                    if nb_positions == 0:
                        fstats.flush()
                    elif nb_positions > 0:
                        csv_filename = str(id_user) + "_" + str(u.username) + '_cytomine_positions.csv'
                        #create output csv file to store positions
                        output_position_file = os.path.join(working_path + project_dir + image_dir + "/user_positions" , csv_filename)
                        f = open(output_position_file, "wb")
                        csvout = csv.writer(f)

                        #Create vector for distribution of zoom levels
                        #Create vector for distribution of days (end_time - start_time)
                        nb_time_intervals = int(math.ceil((end_timestamp-start_timestamp)/timestep))/1000
                        time_intervals = np.zeros(nb_time_intervals+1)
                        central_positions = np.zeros(nb_time_intervals+1)

                        #Filter obtained positions based on start/end timestamp (only write in csv positions included in the given time interval)
                        #Save every position in a csv file
                        previous_central=float(start_timestamp)
                        csvout.writerow(['corners', 'center', 'zoom', 'created'])
                        for p in positions.data():
                            if float(p.created) > float(start_timestamp) and float(p.created) < float(end_timestamp):

                                in_timeinterval = int(math.floor((float(p.created) - float(start_timestamp)) / timestep)) / 1000
                                geom = loads(p.location)
                                if (p.x != image_instance.width / 2) and (p.y != image_instance.height / 2):
                                    rescaled_corners = list(geom.exterior.coords)
                                    rescaled_corners.pop()
                                    rescaled_corners = [tuple(map(lambda divide : round(divide / rescale_factor), corner)) for corner in rescaled_corners]
                                    #csvout.writerow([list(geom.exterior.coords),(int(round(p.x / rescale_factor)), int(round(p.y / rescale_factor))), int(p.zoom), float(p.created)])
                                    csvout.writerow([rescaled_corners, (round(geom.centroid.x / rescale_factor), round(geom.centroid.y / rescale_factor)), int(p.zoom), float(p.created)])
                                    zooms[p.zoom - 1] += 1
                                    time_intervals[in_timeinterval] += 1
                                else:
                                    central_positions[in_timeinterval] += 1
                                    if float(p.created) - previous_central > opening_delay:  # if > 10s we assume user opens the image again
                                        nb_opens += 1
                                    previous_central = float(p.created)
                            else:
                                print "Point removed because not in the timeframe"  # should not happen

                        f.close()

                    #Get Annotations
                    #Get user annotations in this image (to generate statistics about annotation creation).
                    nb_filtered_annotations=0
                    annotations = conn.get_annotations(id_image=id_image,
                                                       id_user=id_user,
                                                       id_project=id_project)
                    nb_annotations = len(annotations.data())
                    csv_filename = str(id_user) + "_" + str(u.username) + '_cytomine_annotations.csv'
                    if not os.path.exists(working_path + project_dir + image_dir + "/user_annotations"):
                        os.makedirs(working_path + project_dir + image_dir + "/user_annotations")
                    if nb_annotations > 0:
                        print "We actually have at least 1 annotation"
                        output_annotation_file = os.path.join(working_path + project_dir + image_dir + "/user_annotations", csv_filename)
                        f = open(output_annotation_file, "wb")
                        csv_annotations = csv.writer(f)
                        csv_annotations.writerow(['type', 'x_center', 'y_center', 'annotationIdent'])
                        for a in annotations.data():
                            if float(a.created) > float(start_timestamp) and float(a.created) < float(end_timestamp):
                                geom = loads(a.location)
                                nb_filtered_annotations += 1
                                if geom.type == 'Point':
                                    csv_annotations.writerow([geom.type,geom.x,geom.y,a.id])
                                else:
                                    csv_annotations.writerow([geom.type, geom.centroid.x, geom.centroid.y,a.id])
                        f.close()

                    #Get AnnotationActions
                    if not os.path.exists(working_path + project_dir + image_dir + "/user_actions"):
                        os.makedirs(working_path + project_dir + image_dir + "/user_actions")
                    pos_success = False
                    while (not pos_success):
                        # Retry if we got error
                        try:
                            ann_actions = conn.get_annoationactions(id_image=id_image,
                                                           id_user=id_user,
                                                           maxperpage=maxperpage,
                                                           afterthan=start_timestamp,
                                                           beforethan=end_timestamp,
                                                           showDetails=True)
                            pos_success = True
                        except socket.error:
                            print socket.error
                            time.sleep(1)
                            continue
                        except socket.timeout:
                            print socket.timeout
                            time.sleep(1)
                            continue
                        except ValueError:
                            print socket.timeout
                            time.sleep(1)
                            continue
                        except Exception:
                            time.sleep(1)
                            continue
                    nb_actions = len(ann_actions.data())
                    if nb_actions == 0:
                        fstats.flush()
                    elif nb_actions > 0:
                        csv_filename = str(id_user) + "_" + str(u.username) + '_cytomine_actions.csv'
                        #create output csv file to store actions
                        output_action_file = os.path.join(working_path + project_dir + image_dir + "/user_actions" , csv_filename)
                        f = open(output_action_file, "wb")
                        csvout = csv.writer(f)

                        #Create vector for distribution of zoom levels
                        #Create vector for distribution of days (end_time - start_time)
                        nb_time_intervals = int(math.ceil((end_timestamp-start_timestamp)/timestep))/1000
                        time_intervals = np.zeros(nb_time_intervals+1)
                        central_positions = np.zeros(nb_time_intervals+1)

                        #Filter obtained positions based on start/end timestamp (only write in csv positions included in the given time interval)
                        #Save every position in a csv file
                        previous_central=float(start_timestamp)
                        csvout.writerow(['annotationIdent', 'created', 'action'])
                        for a in ann_actions.data():
                            if float(a.created) > float(start_timestamp) and float(a.created) < float(end_timestamp):
                                csvout.writerow([a.annotationIdent, a.created, a.action])
                        f.close()

                    csvoutstats.writerow(
                        [id_project, id_image, id_user, u.username, u.email, nb_annotations, nb_opens, sum(zooms), zooms,
                         nb_ref_annotations])

def handle_args(args):

    users = False
    user_info = None
    pub_key = config.cytomine_public_key
    priv_key = config.cytomine_private_key
    host = config.cytomine_host
    images = False
    image_info = None
    modules = None
    try:
        name = str(args[1])
        id_proj = int(args[2])
        users_files = str(args[3])
        ref_user = int(args[4])

        i = 5
        while i < len(args):
            if args[i] == "-U" and (i + 1) < len(args):
                users = True
                user_info = str(args[i + 1])
                i += 2
            elif args[i] == "-I" and (i + 1) < len(args):
                images = True
                image_info = str(args[i + 1])
                i += 2
            elif args[i] == "-PR" and (i + 1) < len(args):
                priv_key = str(args[i + 1])
                i += 2
            elif args[i] == "-PU" and (i + 1) < len(args):
                pub_key = str(args[i + 1])
                i += 2
            elif args[i] == "-H" and (i + 1) < len(args):
                host = str(args[i + 1])
                i += 2
            elif args[i] == "-M" and (i + 1) < len(args):
                modules = str(args[i + 1])
                i += 2

    except:
        error_msg()
        return

    image_list = None
    if images is True:
        f_image_list = open(image_info, "rb")
        csv_in = csv.reader(f_image_list)
        data_im = list(csv_in)
        image_list = []
        for i_tmp in data_im:
            image_list.append(str(i_tmp[0]))
        f_image_list.close()

    user_list = None
    if users is True:
        f_user_list = open(user_info, 'rb')
        csv_in = csv.reader(f_user_list)
        data_user = list(csv_in)
        user_list = []
        for u_tmp in data_user:
            user_list.append(str(u_tmp[0]))
        f_user_list.close()

    get_data(name, id_proj, users_files, ref_user, im_subset=image_list, us_subset=user_list, cytomine_host=host, cytomine_private_key=priv_key,
             cytomine_public_key=pub_key, modules=modules)

def error_msg():
    """
    Output error msg
    :return:
    """

    print "Format : download_data.py <project_name> <project_id> <users_file> <ref_user_id>"
    print "Options (-I and -U uses default settings if not specified):"
    print "  -H <host>:\n    The Cytomine host address\n"
    print "  -PR <cytomine_pr_key>\n    The user's Cytomine private key"
    print "  -PU <cytomine_pu_key>\n    The user's Cytomine public key"
    print "  -I <image_id_file_dir> :\n    CSV file with image IDs, Default takes all images in the project\n    Gets data on the subset of images\n"
    print "  -U <user_id file_dir> :\n    CSV file with user IDs, Default takes all the users in the users metadata file\n    Gets data on the subset of users\n"
    print "  -/m <module file_dir> :\n    CSV file with moduless, Default no modules\n    copies this file\n"

if __name__ == '__main__':

    handle_args(sys.argv)


text:
import numpy as np
from scipy.linalg import block_diag
import cmath;

def homogeneous_module(Kx, Ky, e_r, m_r = 1):
    '''
    homogeneous layer is much simpler to do, so we will create an isolated module to deal with it
    :return:
    '''
    assert type(Kx) == np.ndarray, 'not np.array'
    assert type(Ky) == np.ndarray, 'not np.array'
    j = cmath.sqrt(-1);
    N = len(Kx);
    I = np.identity(N);
    P = (e_r**-1)*np.block([[Kx*Ky, e_r*m_r*I-Kx**2], [Ky**2-m_r*e_r*I, -Ky*Kx]])
    Q = (e_r/m_r)*P;
    W = np.identity(2*N)
    arg = (m_r*e_r*I-Kx**2-Ky**2); #arg is +kz^2
    arg = arg.astype('complex');
    Kz = np.conj(np.sqrt(arg)); #conjugate enforces the negative sign convention (we also have to conjugate er and mur if they are complex)
    eigenvalues = block_diag(j*Kz, j*Kz) #determining the modes of ex, ey... so it appears eigenvalue order MATTERS...
    #W is just identity matrix
    V = Q@np.linalg.inv(eigenvalues); #eigenvalue order is arbitrary (hard to compare with matlab
    #alternative V with no inverse
    #V = np.matmul(np.linalg.inv(P),np.matmul(Q,W)); apparently, this fails because P is singular
    return W,V,Kz

def homogeneous_1D(Kx, e_r, m_r = 1):
    '''
    efficient homogeneous 1D module
    :param Kx:
    :param e_r:
    :param m_r:
    :return:
    '''
    j = cmath.sqrt(-1);

    I = np.identity(len(Kx));
    P = e_r*I - Kx**2;
    Q = I;
    arg = (m_r*e_r*I-Kx**2); #arg is +kz^2
    arg = arg.astype('complex');

    Kz = np.conj(np.sqrt(arg)); #conjugate enforces the negative sign convention (we also have to conjugate er and mur if they are complex)
    eigenvalues = j*Kz #determining the modes of ex, ey... so it appears eigenvalue order MATTERS...
    V = np.matmul(Q,eigenvalues); #eigenvalue order is arbitrary (hard to compare with matlab
    return I,V, Kz

max_stars_repo_name: zhaonat/RCWA

reason: import numpy,from scipy


text:
"""Basic module in NvTK.
This module provides

1.  `BasicModule` class - the general abstract class

2.  `BasicConv1d` class - Basic Convolutional Module (1d)

3.  `BasicRNNModule` class - Basic RNN(LSTM) Module in batch-first style

4.  `BasicLinearModule`

5.  `BasicPredictor` Module

6.  `BasicLoss` Module

and supporting methods.
"""

import random, time, logging
import numpy as np

import torch
from torch import nn

__all__ = ["BasicModule", "BasicConv1d", "BasicRNNModule", "BasicLinearModule", "BasicPredictor", "BasicLoss", "Flatten"]

class BasicModule(nn.Module):
    """Basic module class in NvTK."""

    def __init__(self):
        super(BasicModule,self).__init__()

    def initialize_weights(self):
        """initialize module parameters.

        Conv module weight will be initialize in xavier_normal_,
        bias will be initialize in zero_

        Linear module weight will be initialize in xavier_normal_,
        bias will be initialize in zero_

        BatchNorm module weight will be initialize in constant = 1,
        bias will be initialize in constant = 0

        LSTM module weight will be initialize in orthogonal_
        """
        for m in self.modules():
            if isinstance(m, (nn.Conv1d, nn.Conv2d)):
                nn.init.xavier_normal_(m.weight)
                if m.bias is not None:
                    m.bias.data.zero_()
                logging.debug("init Conv param...")

            elif isinstance(m, nn.Linear):
                nn.init.xavier_normal_(m.weight)
                if m.bias is not None:
                    m.bias.data.zero_()
                logging.debug("init Linear param...")

            elif isinstance(m, (nn.BatchNorm1d, nn.BatchNorm2d)):
                nn.init.constant_(m.weight, 1)
                nn.init.constant_(m.bias, 0)
                logging.debug("init BatchNorm param...")

            elif isinstance(m, nn.LSTM):
                nn.init.orthogonal_(m.all_weights[0][0])
                nn.init.orthogonal_(m.all_weights[0][1])
                nn.init.orthogonal_(m.all_weights[1][0])
                nn.init.orthogonal_(m.all_weights[1][1])
                logging.debug("init LSTM param...")

    def initialize_weights_from_pretrained(self, pretrained_net_fname):
        """initialize module weights from pretrained model

        Parameters
        ----------
        pretrained_net_fname : str
            the pretrained model file path (e.g. `checkpoint.pth`).
        """
        pretrained_dict = torch.load(pretrained_net_fname)
        net_state_dict = self.state_dict()
        pretrained_dict = {k: v for k, v in pretrained_dict.items() if k in net_state_dict}
        net_state_dict.update(pretrained_dict)
        self.load_state_dict(net_state_dict)
        logging.info("params loaded from: %s" % pretrained_net_fname)

    def load(self, path):
        """load module weights from saved model

        Parameters
        ----------
        path : str
            the saved model file path (e.g. `checkpoint.pth`).
        """
        self.load_state_dict(torch.load(path, map_location=torch.device('cpu')))

    def save(self, fname=None):
        """save module weights to file

        Parameters
        ----------
        fname : str, optional
            Specify the saved model file path.
            Default is "None". Saved file will be formatted as "model.time.pth".
        """
        if fname is None:
            fname = time.strftime("model" + '%m%d_%H:%M:%S.pth')
        torch.save(self.state_dict(), fname)
        return fname

    # TODO (jiaqili, 20220417): test the tensor flow in module forward
    def test(self, input_size):
        device = list(self.parameters)[0].device
        x = torch.zeros(input_size).to(device)
        self.forward(x)
        logging.info("Test: all the tensor flow shape reported")

class Flatten(nn.Module):
    """Flatten Module: flatten the tensor as (batch_size, -1)."""
    def forward(self, x):
        return x.view(x.size(0), -1)

class EXP(nn.Module):
    """Exp Module: calculate the exp of tensor as `x.exp()`"""
    def forward(self, x):
        return x.exp()

class BasicConv1d(BasicModule):
    """Basic Convolutional Module (1d) in NvTK.

    Parameters
    ----------
    in_planes : int
        Number of input channels
    out_planes : int
        Number of output channels produced by the convolution
    kernel_size : int, optional
        Size of the convolving kernel
    conv_args : dict, optional
        Other convolutional args, Default is dict().
        Will be pass to `torch.nn.Conv1d(**conv_args)`
        (e.g. `conv_args={'dilation':1}`)
    bn : bool, optional
        Whether to use BatchNorm1d, Default is True.
    activation : nn.Module, optional
        Activation Module, Default is nn.ReLU.
    activation_args : dict, optional
        Other activation args, Default is dict().
        Will be pass to `activation(**activation_args)`
        (e.g. `activation=nn.LeakyReLU, activation_args={'p':0.2}`)
    dropout : bool, optional
        Whether to use Dropout, Default is True.
    dropout_args : dict, optional
        Dropout args, Default is {'p':0.5}.
        Will be pass to `nn.Dropout(**dropout_args)` if dropout
        (e.g. `dropout=True, dropout_args={'p':0.5}`)
    pool : nn.Module, optional
        Pool Module (1d), Default is nn.AvgPool1d.
    pool_args : dict, optional
        Other pool args, Default is {'kernel_size': 3}.
        Will be pass to `pool(**pool_args)`
        (e.g. `pool=nn.AvgPool1d, pool_args={'kernel_size': 3}`)

    Attributes
    ----------
    in_channels : int

    out_channels : int

    conv : nn.Conv1d
        The convolutional neural network component of the model.
    bn : nn.BatchNorm1d
        The Batch Normalization
    activation : nn.Module
        The activation Module
    dropout : nn.Dropout
        The Dropout Module
    pool : nn.Module
        The pool Module

    Tensor flows
    ----------
    -> conv(x)

    -> bn(x) if bn

    -> activation(x) if activation

    -> dropout(x) if dropout

    -> pool(x) if pool

    """
    def __init__(self, in_planes, out_planes, kernel_size=3, conv_args={},
                    bn=True,
                    activation=nn.ReLU, activation_args={},
                    dropout=True, dropout_args={'p':0.5},
                    pool=nn.AvgPool1d, pool_args={'kernel_size': 3}):
        super().__init__()
        self.in_channels = in_planes
        self.out_channels = out_planes
        self.conv = nn.Conv1d(in_planes, out_planes, kernel_size=kernel_size, **conv_args)
        self.bn = nn.BatchNorm1d(out_planes, eps=1e-5, momentum=0.01, affine=True) if bn else None
        self.activation = activation(**activation_args) if activation is not None else None
        self.dropout = nn.Dropout(**dropout_args) if dropout else None
        self.pool = pool(**pool_args) if pool else None

    def forward(self, x):
        x = self.conv(x)
        if self.bn is not None:
            x = self.bn(x)
        if self.activation is not None:
            x = self.activation(x)
        if self.dropout is not None:
            x = self.dropout(x)
        if self.pool is not None:
            x = self.pool(x)
        logging.debug(x.shape)
        return x

class BasicRNNModule(BasicModule):
    """
    Basic RNN(LSTM) Module in batch-first style
    """
    def __init__(self, LSTM_input_size=512, LSTM_hidden_size=512, LSTM_hidden_layes=2):
        super().__init__()
        self.rnn_hidden_state = None
        self.rnn = nn.LSTM(
            input_size=LSTM_input_size,
            hidden_size=LSTM_hidden_size,
            num_layers=LSTM_hidden_layes,
            batch_first=True, # batch, seq, feature
            bidirectional=True,
        )
    def forward(self, input):
        output, self.rnn_hidden_state = self.rnn(input, None)
        logging.debug(output.shape)
        return output

class BasicLinearModule(BasicModule):
    """
    Basic Linear Module in NvTK.

    Parameters
    ----------
    input_size : int
        Number of input size
    output_size : int
        Number of output size produced by the Linear
    bias : bool, optional
        Bias of the Linear, Default is True.
        It could be False when use BatchNorm.
    bn : bool, optional
        Whether to use BatchNorm1d, Default is True.
    activation : nn.Module, optional
        Activation Module, Default is nn.ReLU.
    activation_args : dict, optional
        Other activation args, Default is dict().
        Will be pass to `activation(**activation_args)`
        (e.g. `activation=nn.LeakyReLU, activation_args={'p':0.2}`)
    dropout : bool, optional
        Whether to use Dropout, Default is True.
    dropout_args : dict, optional
        Dropout args, Default is {'p':0.5}.
        Will be pass to `nn.Dropout(**dropout_args)` if dropout
        (e.g. `dropout=True, dropout_args={'p':0.5}`)

    Attributes
    ----------
    input_size : int

    output_size : int

    linear : nn.Linear
        The Linear neural network component
    bn : nn.BatchNorm1d
        The Batch Normalization
    activation : nn.Module
        The activation Module
    dropout : nn.Dropout
        The Dropout Module

    Tensor flows
    ----------
    -> linear(x)

    -> bn(x) if bn

    -> activation(x) if activation

    -> dropout(x) if dropout

    """
    def __init__(self, input_size, output_size, bias=True, bn=True,
                    activation=nn.ReLU, activation_args={},
                    dropout=True, dropout_args={'p':0.5}):
        super().__init__()
        self.input_size = input_size
        self.output_size = output_size
        self.linear = nn.Linear(input_size, output_size, bias=bias)
        self.bn = nn.BatchNorm1d(output_size, eps=1e-5, momentum=0.01, affine=True) if bn else None
        self.activation = activation(**activation_args) if activation is not None else None
        self.dropout = nn.Dropout(**dropout_args) if dropout else None

    def forward(self, x):
        x = self.linear(x)
        if self.bn is not None:
            x = self.bn(x)
        if self.activation is not None:
            x = self.activation(x)
        if self.dropout is not None:
            x = self.dropout(x)
        logging.debug(x.shape)
        return x

class BasicPredictor(BasicModule):
    """BasicPredictor Module in NvTK.

    BasicPredictor support task types of 'none', 'binary_classification', 'classification', 'regression';
    1. 'none' : nullify the whole BaiscPredictor with identity
    2. 'binary_classification' : activate with Sigmoid
    3. 'classification' : activate with Softmax(dim=1)
    4. 'regression' : Identity

    Parameters
    ----------
    input_size : int
        Number of input size
    output_size : int
        Number of output size (task numbers)
    tasktype : str, optional
        Specify the task type, Default is "binary_classification".
        (e.g. `tasktype="regression"`)

    Attributes
    ----------
    supported_tasks : currently supported task types
    tasktype : task type of Predictor
    input_size : int
    Map : nn.Linear
        The Linear Module Mapping input to output.
    Pred: nn.Module
        The Activation Module in specified task type.

    Tensor flow
    ----------
    -> Map(x)

    -> Pred(x)

    """
    def __init__(self, input_size, output_size, tasktype='binary_classification'):
        super().__init__()
        self.supported_tasks = ['none', 'binary_classification', 'classification', 'regression']

        self.input_size = input_size
        self.tasktype = tasktype

        self.Map = nn.Linear(input_size, output_size, bias=True)
        self.switch_task(tasktype) # init self.Map

    def forward(self, x):
        return self.Pred(self.Map(x))

    def switch_task(self, tasktype):
        """switch to specified task type

        Parameters
        ----------
        tasktype : str
            Specify the task type (e.g. `tasktype="regression"`)
        """

        msg = 'tasktype: %s not supported, check the document' % tasktype
        assert tasktype in self.supported_tasks, msg

        if tasktype == 'none':
            self.Map = nn.Sequential()
            self.Pred = nn.Sequential()
        elif tasktype == 'classification':
            self.Pred = nn.Softmax(dim=1)
        elif tasktype == 'binary_classification':
            self.Pred = nn.Sigmoid()
        elif tasktype == 'regression':
            self.Pred = nn.Sequential()

        self.tasktype = tasktype

    # TODO
    # def finetune(self, output_size, tasktype='binary_classification'):
    #     self.Map = nn.Linear(input_size, output_size, bias=True)
    #     self.Pred = nn.Sequential()

    #     self.tasktype = tasktype
    #     self.switch_task(tasktype)

    def current_task(self):
        """return current task type"""
        return self.tasktype

    def remove(self):
        """Predictor.remove: replace predictor with null Sequential,
        same as switch_task('none')
        """
        self.switch_task('none')

class BasicLoss(nn.Module):
    """BasicLoss Module in NvTK.

    BasicLoss support task types of 'binary_classification', 'classification', 'regression';
    1. 'binary_classification' : BCELoss function
    2. 'classification' : CrossEntropyLoss function
    3. 'regression' : MSELoss function

    Parameters
    ----------
    tasktype : str, optional
        Specify the task type, Default is "binary_classification".
        (e.g. `tasktype="regression"`)
    reduction : str, optional
        Specifies the reduction to apply to the output: `'none'` | `'mean'` | `'sum'`.

    Attributes
    ----------
    supported_tasks : currently supported task types
    tasktype : task type of Predictor
    loss : loss function

    """
    def __init__(self, tasktype='binary_classification', reduction='mean'):
        super().__init__()
        self.supported_tasks = ['binary_classification', 'classification', 'regression']

        self.tasktype = tasktype
        self.reduction = reduction

        self.switch_task(tasktype) # init self.loss

    def forward(self, pred, target):
        return self.loss(pred, target)

    def switch_task(self, tasktype):
        """switch to specified task type

        Parameters
        ----------
        tasktype : str
            Specify the task type (e.g. `tasktype="regression"`)
        """

        msg = 'tasktype: %s not supported, check the document' % tasktype
        assert tasktype in self.supported_tasks, msg

        if tasktype == 'classification':
            self.loss = nn.CrossEntropyLoss(reduction=self.reduction)
        elif tasktype == 'binary_classification':
            self.loss = nn.BCELoss(reduction=self.reduction)
        elif tasktype == 'regression':
            self.loss = nn.MSELoss(reduction=self.reduction)

        self.tasktype = tasktype

max_stars_repo_name: JiaqiLiZju/NvTK

max_issues_repo_name: JiaqiLiZju/NvTK

max_forks_repo_name: JiaqiLiZju/NvTK


text:
import os, glob
import h5py
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def main():
    num = extract_num_trajectory_from_out_files()
    an = eval_an_trajectory_at_num_time_points( num )
    plot_trajectories( num , an )

def extract_num_trajectory_from_out_files():
    out_files = find_necessary_out_files()

    num_trajectory = []
    for f in out_files:
        num_trajectory.append( extract_time_pos_mom( f ) )

    num_trajectory = remove_empty_and_sort_by_time( num_trajectory )
    num_trajectory = np.array( num_trajectory,
                               dtype=[('t','float'),
                                      ('x','float'), ('y','float'), ('z','float'),
                                      ('px','float'), ('py','float'), ('pz','float') ] )
    return( num_trajectory )

def remove_empty_and_sort_by_time( num_trajectory ):
    removed_empty = [ x for x in num_trajectory if x ]
    sorted_by_time = sorted( removed_empty, key = lambda x: x[0] )
    return ( sorted_by_time )

def find_necessary_out_files():
    os.chdir("./")
    h5files = []
    for file in glob.glob("single_particle_electric_field_[0-9]*.h5"):
        h5files.append( file )
    return h5files

def extract_time_pos_mom( h5file ):
    h5 = h5py.File( h5file, mode="r")
    t = h5["/Time_grid"].attrs["current_time"][0]
    t_pos_mom = ()
    if ( len(h5["/Particle_sources/emit_single_particle"]) > 0 ):
        x = h5["/Particle_sources/emit_single_particle/position_x"][0]
        y = h5["/Particle_sources/emit_single_particle/position_y"][0]
        z = h5["/Particle_sources/emit_single_particle/position_z"][0]
        px = h5["/Particle_sources/emit_single_particle/momentum_x"][0]
        py = h5["/Particle_sources/emit_single_particle/momentum_y"][0]
        pz = h5["/Particle_sources/emit_single_particle/momentum_z"][0]
        t_pos_mom = (t, x, y, z, px, py, pz)
    h5.close()
    return( t_pos_mom )

def eval_an_trajectory_at_num_time_points( num_trajectory ):
    global mass, charge, E0
    mass, charge, x0, y0, z0, px0, py0, pz0 = get_mass_charge_and_initial_pos_and_mom()
    E0 = eval_field_amplitude()

    an_trajectory = np.empty_like( num_trajectory )
    for i, t in enumerate( num_trajectory['t'] ):
        x, y, z = coords( t, x0, y0, z0, px0, py0, pz0 )
        px, py, pz = momenta( t, px0, py0, pz0 )
        an_trajectory[i] = ( t, x, y ,z, px, py, pz )

    return( an_trajectory )

def get_mass_charge_and_initial_pos_and_mom():
    initial_out_file = "single_particle_electric_field_0000000.h5"
    h5 = h5py.File( initial_out_file, mode="r")
    m = h5["/Particle_sources/emit_single_particle"].attrs["mass"][0]
    q = h5["/Particle_sources/emit_single_particle"].attrs["charge"][0]
    x0 = h5["/Particle_sources/emit_single_particle/position_x"][0]
    y0 = h5["/Particle_sources/emit_single_particle/position_y"][0]
    z0 = h5["/Particle_sources/emit_single_particle/position_z"][0]
    px0 = h5["/Particle_sources/emit_single_particle/momentum_x"][0]
    py0 = h5["/Particle_sources/emit_single_particle/momentum_y"][0]
    pz0 = h5["/Particle_sources/emit_single_particle/momentum_z"][0]
    h5.close()
    return( m, q, x0, y0, z0, px0, py0, pz0 )

def eval_field_amplitude():
    initial_out_file = "single_particle_electric_field_0000000.h5"
    h5 = h5py.File( initial_out_file, mode="r")
    # todo: boundary conditions are not saved explicitly in out file.
    # Attempt to extract them results in cryptic manipulations
    # with number of nodes in different directions of mesh.
    # Need to fix it somehow.
    x_n_nodes = h5["/Spatial_mesh"].attrs["x_n_nodes"][0]
    y_n_nodes = h5["/Spatial_mesh"].attrs["y_n_nodes"][0]
    z_n_nodes = h5["/Spatial_mesh"].attrs["z_n_nodes"][0]
    z_volume_size = h5["/Spatial_mesh"].attrs["z_volume_size"][0]
    # first point with X != 0, Y != 0, Z == 0
    phi_near  = h5["/Spatial_mesh/potential"][ z_n_nodes * y_n_nodes + z_n_nodes ]
    # first point with X != 0, Y != 0, Z == Z_size
    phi_far  = h5["/Spatial_mesh/potential"][ z_n_nodes * y_n_nodes + z_n_nodes + z_n_nodes - 1 ]
    E0 = - ( phi_far - phi_near ) / z_volume_size
    return E0

def momenta( t, px0, py0, pz0 ):
    global mass, charge, E0
    px = px0
    py = py0
    pz = pz0 + charge / mass * E0 * t
    return ( px, py, pz )

def coords( t, x0, y0, z0, px0, py0, pz0 ):
    global mass, charge, E0
    x = x0 + px0 / mass * t
    y = y0 + py0 / mass * t
    z = z0 + pz0 / mass * t + 1 / 2 * charge / mass * E0 * t * t
    return ( x, y, z )

def plot_trajectories( num , an ):
    plot_3d( num, an )
    plot_2d( num, an )
    plot_kin_en( num , an )

def plot_3d( num, an ):
    fig = plt.figure()
    ax = fig.gca( projection='3d' )
    ax.plot( num['x'], num['y'], num['z'], '.r', label = "Num" )
    ax.plot( an['x'], an['y'], an['z'], label = "An" )
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z')
    plt.legend( loc = 'upper left', title="3d" )
    #plt.show()
    print( 'Saving 3d trajectory plot to "3d.png"' )
    plt.savefig('3d.png')

def plot_2d( num, an ):
    plt.figure(1)
    #XY
    plt.subplot(131)
    plt.plot( num['x'], num['y'],
              linestyle='', marker='o',
              label = "Num" )
    plt.plot( an['x'], an['y'],
              linestyle='-', marker='', lw = 2,
              label = "An" )
    plt.legend( loc = 'upper right', title="XY" )
    #XZ
    plt.subplot(132)
    plt.plot( num['x'], num['z'],
        linestyle='', marker='o',
        label = "Num" )
    plt.plot( an['x'], an['z'],
              linestyle='-', marker='', lw = 2,
              label = "An" )
    plt.legend( loc = 'upper right', title="XZ" )
    #YZ
    plt.subplot(133)
    plt.plot( num['y'], num['z'],
        linestyle='', marker='o',
        label = "Num" )
    plt.plot( an['y'], an['z'],
              linestyle='-', marker='', lw = 2,
              label = "An" )
    plt.legend( loc = 'upper right', title="YZ" )
    print( 'Saving 2d trajectory projection plots to "2d.png"' )
    plt.savefig('2d.png')

def plot_kin_en( num , an ):
    global mass
    E_num = ( num['px']**2 + num['py']**2 + num['pz']**2 ) / ( 2 * mass )
    E_an = ( an['px']**2 + an['py']**2 + an['pz']**2 ) / ( 2 * mass )
    t = num['t']
    plt.figure()
    axes = plt.gca()
    axes.set_xlabel( "t [s]" )
    axes.set_ylabel( "E [erg]" )
    # axes.set_ylim( [min( E_an.min(), E_num.min() ),
    #                 max( E_an.max(), E_num.max() ) ] )
    line, = plt.plot( t, E_num, 'o' )
    line.set_label( "Num" )
    line, = plt.plot( t, E_an, ls = 'solid', lw = 3 )
    line.set_label( "An" )
    plt.legend( loc = 'upper right' )
    print( 'Saving kinetic energy comparison plot to "kin_en.png"' )
    plt.savefig( 'kin_en.png' )

main()


text:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy.stats as scistats
from matplotlib.backends.backend_pdf import PdfPages
from matplotlib.colors import LogNorm
import mne.time_frequency as mtf

#FFT that will zero pad a signal with NANs.
#This allows for greater resolution in the frequency domain
def nanFFT(x,n):
    fft = np.fft.fft(np.nan_to_num(x),n=n) #fft of non-nan components
    return fft

#Simple calculation for the power of an fft bin
def powerSpec(x,n=None):
    if n is None:
        n = len(x)
    x = nanFFT(x,n)
    return np.square(abs(x))

#Will save the short time fourier transform
def save_stft(powMap,ws,delta,grid,fs=None,channels=None,name=None):

    if fs is None:
        fs = grid.fs()
    if channels is None:
        channels = grid.wc()
    if name is None:
        name = 'STFT.pdf'
    elif name.find('.pdf') == -1:
        name += '.pdf'

    freq = np.fft.fftfreq(int(ws*fs),1.0/fs)
    freq = freq[:len(freq)/2]
    (x,y,z) = powMap.shape
    #pp = PdfPages(grid.wc()[0]+'-'+grid.wc()[-1]+'_'+name)
    fig = plt.figure()

    for cnum,chan in enumerate(channels):
        #
        plt.pcolor(powMap[cnum,:,:10].T,norm=LogNorm(vmin=powMap[cnum,:,:10].min(), vmax=powMap[cnum,:,:10].max()),cmap='blues')
        #
        xTickMultiple=200
        xlabels = [t for (i,t) in enumerate(np.arange(y)*delta/fs) if not i%xTickMultiple]
        plt.xticks(np.arange(len(xlabels)),xlabels,rotation=45)
        plt.axes().xaxis.set_major_locator(plt.MultipleLocator(xTickMultiple))
        #
        yTickMultiple = 1
        ylabels = [hz for (i,hz) in enumerate(freq[:11]) if not i%yTickMultiple]
        plt.yticks(np.arange(len(ylabels)),ylabels)
        plt.axes().yaxis.set_major_locator(plt.MultipleLocator(yTickMultiple))

        plt.ylabel('Hz')
        #
        plt.xlabel('Seconds')
        plt.title('STFT Power '+chan)
        plt.colorbar()
        #
        print 'Saving channel '+chan
        plt.savefig('STFT_'+chan)
        #pp.savefig()
        fig.clf()
        fig = plt.figure(name)

    #pp.close()

def morlet(grid,dec=None,frequencies=None,fs=None,events=None,nc=None,normFreqs=None):
    #Performs a morlet wavelet analysis across events

    if events is None:
        events = grid.event_matrix()
    if fs is None:
        fs = grid.fs()
    if dec is None:
        dec = 1
    if frequencies is None:
        frequencies= np.arange(1,101,dec)
    if nc is None:
        nc = 0.12*np.arange(1,len(frequencies)+1)
        print nc

    sigs = grid.wd()
    maxtimediff = (events['pulse.off']-events['pulse.on']).max()
    maxtimepoints =  round(maxtimediff*fs)+1
    data = np.zeros([len(events),len(grid.wc()),maxtimepoints])

    for k,chan in enumerate(grid.wc()):
        for j,(on,off) in enumerate(events.values):
            print 'Pre-processing Channel:'+chan + ' | Event '+str(j)+'/'+str(len(events.values))
            d=grid.splice(sigs[chan], times=[on,off])
            data[j,k,:len(d)] = d
    pows,phase=mtf.induced_power(data, Fs=float(fs), frequencies=frequencies, use_fft=False, n_cycles=nc, decim=dec, n_jobs=1,normFreqs=normFreqs)
    return pows,phase

def mne_stft(grid,ws,delta,sigs=None,fs=None,events=None):
    #Performs a stft transform across the events????????????

    if events is None:
        events = grid.event_matrix()
    if fs is None:
        fs = grid.fs()
    if sigs is None:
        sigs = grid.wd()

    ws*=fs
    delta*=fs
    maxtimediff = (events['pulse.off']-events['pulse.on']).max()
    maxtimepoints =  round(maxtimediff*fs/delta) - np.floor(ws/delta)
    powMap = np.zeros([len(events),len(grid.wc()),ws])
    data = np.zeros([len(events),len(grid.wc()),maxtimepoints])

    for j,(on,off) in enumerate(events.values):
        for k,chan in enumerate(grid.wc()):
            d=grid.splice(sigs[chan], times=[on,off])
            data[j,k,:len(d)] = mtf.stft(x=d,wsize=ws,tstep=delta)

    return powMap.mean(axis=2)

def plot_morlet(morlet,grid,ws,frequencies=None,fs=None,dec=None,name=None):

    if fs is None:
        fs = grid.fs()
    if dec is None:
        dec = 1
    if frequencies is None:
        frequencies = np.arange(1,101,1)
    if name is None:
        name = 'Morlet'

    (x,y,z) = morlet.shape
    times = np.arange(0,z/float(fs),ws)

    for k,chan in enumerate(grid.wc()):

        x = morlet[k,:,:]
        fig, (ax1, ax2) = plt.subplots(nrows=2)
        data, freqs, bins, im = ax1.specgram(x)
        ax1.axis('tight')
        # We need to explictly set the linear threshold in this case...
        # Ideally you should calculate this from your bin size...
        ax2.set_yscale('symlog', linthreshy=0.01)
        ax2.pcolormesh(bins, freqs, 10 * np.log10(data))
        ax2.axis('tight')
        plt.xlabel('Time (s)')
        plt.ylabel('Frequency (Hz)')
        plt.title('Induced power (%s)' % chan)
        plt.savefig(name+chan)
        plt.clf()
        plt.figure()

def stft(grid,ws,delta,sigs=None,fs=None,events=None):
    #Performs a short time sliding window fourier transform
    #Averages a moving window across signals
    #Note that the algorithm processes an entire signal
    #and moves to the next instead of loading all signals
    #into memory and averaging down the signals

    if events is None:
        events = grid.event_matrix()
    if fs is None:
        fs = grid.fs()
    if sigs is None:
        sigs = grid.wd()

    ws*=fs
    delta*=fs
    maxtimediff = (events['pulse.off']-events['pulse.on']).max()
    maxtimepoints =   round(maxtimediff*fs/delta) - np.floor(ws/delta)

    powMap = np.zeros([len(grid.wc()),maxtimepoints,ws])

    for k,chan in enumerate(grid.wc()):
        #
        #Iterate events
        counter = np.zeros([maxtimepoints,ws])#Keep track of the number of events that contribute to the total time band
        for j,(on,off) in enumerate(events.values):
            print 'Processing Channel:'+chan + ' | Event '+str(j)+'/'+str(len(events.values))
            #
            #Move along each event, shifting the window by delta
            d=grid.splice(sigs[chan], times=[on,off])
            starts = np.linspace(0, int(len(d))-ws, maxtimepoints)#range(0,int(len(d)),int(delta))
            #
            for i,start in enumerate(starts):
                #
                #Get the current window splice
                s = d[start:start+ws]
                if not len(s)%2: #Ensure d is evenly sized for efficiency
                    s = s[:-1]
                #
                pwr = powerSpec(s)
                powMap[k,i,:len(pwr)] += pwr
                counter[i,:len(pwr)]+= np.ones(len(pwr))
        powMap[k,:,:]= (powMap[k,:,:]/counter) #Get the mean across the channel
    powMap[np.isnan(powMap)]=0 #Remove nan columns
    return powMap

def calcFFT(grid,events,ws=None):
    from signalUtils import shortest_event

    if ws is None:
        ws = shortest_event(grid,events)
        print '\tUsing minimum event size of '+str(ws)+' as size of FFT window'

    if ws%2:
        ws -= 1
        print '\tReducing FFT window size to be even (efficiency)'

    #FFT across all signals and events
    blockFFTs = np.zeros([len(grid.wc()),len(events),ws])

    for i,chan in enumerate(grid.wc()):
        sig  = grid.wd(channels=[chan])[chan]

        print "\tCalulcating fft for Channel: "+chan
        for row,(on,off) in enumerate(events[['pulse.on','pulse.off']].values):
            block = grid.splice(data=sig,times=[on,off])[:ws] #Slice out the event and truncate to the length of the shortest
            blockFFTs[i,row,:] = block
            blockFFTs[i,row,:] = np.power(np.abs(np.fft.fft(blockFFTs[i,row,:])),2)
    return blockFFTs

def save_fft(grid,stims,isis,name='',maxFreq=None):

    pp = PdfPages(name+'.pdf')
    (x,y,z) = stims.shape
    (x2,y2,z2) = isis.shape
    if x!=x2:
        print 'Number of channels does not match'
        print x,x2
        return
    elif z != z2:
        print 'FFTs are uneven lengths'
        print z,z2
        return
    elif y!= y2:
        print 'Warning - Number of stimuli events does not match the number of baseline periods'

    chans = grid.wc()
    freq = np.fft.fftfreq(z,1.0/grid.fs())
    freq = freq[:z/2]
    maxFreq = freq[-1] if maxFreq is None else maxFreq
    freq = freq[np.where(freq<maxFreq)]
    z = len(freq)

    for i in range(x):
        print 'Saving channel '+chans[i]
        plt.figure()
        plt.hold(True)
        plt.title(chans[i])

        stimprc50 = np.percentile(stims[i,:,:z],50,axis=0)
        stimprc25 = np.percentile(stims[i,:,:z],25,axis=0)
        stimprc75 = np.percentile(stims[i,:,:z],75,axis=0)

        isiprc50 = np.percentile(isis[i,:,:z],50,axis=0)
        isiprc25 = np.percentile(isis[i,:,:z],25,axis=0)
        isiprc75 = np.percentile(isis[i,:,:z],75,axis=0)

        plt.loglog(freq,stimprc50,'b',label='Mean stimulus')
        plt.loglog(freq,stimprc25,'b.')
        plt.loglog(freq,stimprc75,'b--')

        plt.loglog(freq,isiprc50,'g',label='Mean ISI')
        plt.loglog(freq,isiprc25,'g.')
        plt.loglog(freq,isiprc75,'g--')

        plt.legend(loc=1)
        plt.grid(True,which='both')
        pp.savefig()
        plt.clf()
        plt.close()

    pp.close()

reason: import numpy,import scipy


text:
# -*- coding: utf-8 -*-
"""
Created on Mon Nov 12 13:04:40 2018

@author: Arvinder Shinh
"""

import pandas as pd
import numpy as np

def transformer(frame, catagoryColumn='tag_category',
                tagColumn='matched_tag',
                text='text',
                shortDescription='preview_text',
                shortDescriptionLen='200',
                columnList=['heading','preview_text','timestamp','source_link','html_file','tag_category','matched_tag'],
                noTagColumn='other_tags',
                noTagValue='other'):

    catagories = [x for x in frame.columns if x.endswith("tags")]
    dataColumns = [x for x in frame.columns if not x.endswith("tags")]

    frame[noTagColumn]=np.NaN

    condition = True
    for catagory in catagories:
           condition = condition & pd.isna(frame[catagory]).values

    frame.loc[condition,[noTagColumn]]=noTagValue

    catagories.append(noTagColumn)

    X, Y = len(frame.index), len(catagories)
    data = {dataColumn : np.tile(np.asarray(frame[dataColumn]), Y) for dataColumn in dataColumns}

    data[catagoryColumn] = np.asarray(catagories).repeat(X)
    data[tagColumn] = frame.loc[:, catagories].values.ravel('F')

    dataFrame = pd.DataFrame(data)
    dataFrame.dropna(subset=[tagColumn], inplace=True)

    dataDict = {}
    j = 0
    for i in range(dataFrame.shape[0]):
             series=dataFrame.iloc[i]

             for tag in series[tagColumn].split('#'):
                 if not (tag.strip() == ""):
                    s = series.copy()
                    s[tagColumn] = tag.strip()
                    dataDict[j] = s
                    j = j+1

    pattern='(^.{'+shortDescriptionLen+'})*'

    dataFrameTransformed = pd.DataFrame(dataDict).transpose()
    dataFrameTransformed[shortDescription]=dataFrameTransformed[text].str.extract(pattern, expand=False)

    return dataFrameTransformed.loc[:,columnList]

dfOriginal = pd.read_csv('C:/Workspace/PythonProject/Django/UrlTagMap/subset_compiled_data.csv')
dfTransformed = transformer(dfOriginal)

dfTransformed.to_csv('C:/Workspace/PythonProject/Django/UrlTagMap/subset_compiled_data_Out.csv',  index=False)

max_stars_repo_name: arvindershinh/WebSearchEngine

max_issues_repo_name: arvindershinh/WebSearchEngine

max_forks_repo_name: arvindershinh/WebSearchEngine


text:
# -*- coding: utf-8 -*-
#
# Licensed under the terms of the MIT License

"""."""

import numpy as np
from numpy import ndarray, arange, angle, rad2deg, unwrap
from numpy.fft import fft

def db2(data, db_range=60.0, db_cut=None, b_normalize=True):
    """Convert 20log10(data).

    Parameters
    ----------
    data : ndarray
        dB 변환 할 데이터
    db_range : float
        dB Dynamic Range (the default is 60.0)
    db_cut : float, optional
        dB 변환 시 최소값 (the default is None, None이면 미 수행), (b_normalize가 True 이면 Normalize 이후 cut 값)
    b_normalize : bool, optional
        최대값을 0으로 변환 여부 (the default is True)

    Returns
    -------
    ndarray
        20log10(data)
    """

    data_abs = np.abs(data)
    max_v = data_abs.max()
    db_max = 20 * np.log10(max_v)

    with np.errstate(divide='ignore'):
        db_data = 20 * np.log10(data_abs)
    db_data[db_data < (db_max - db_range)] = db_max - db_range

    if b_normalize:
        db_data -= db_max

    if db_cut:
        db_data[db_data < db_cut] = db_cut
    return db_data

def frequency_response(data, fs, uprate):
    """

    Parameters
    ----------
    data: ndarray
    fs: float
        sampling frequency
    uprate: int

    Returns
    -------
    (ndarray, ndarray, ndarray)
        mag_fre_db, phase_fre_rad, frequnecy
    """
    response = fft(data, len(data) * uprate)
    mag_db = db2(response, 90)[:len(response)//2]
    phase_rad = angle(response)[:len(response) // 2]
    phase_rad = unwrap(phase_rad)

    n_sample = len(mag_db)
    frequency = arange(n_sample) / n_sample * fs / 2

    return mag_db, phase_rad, frequency

reason: import numpy,from numpy