text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you 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.

"""
Tests for individual operators
This module contains operator tests which currently do not exist on
ONNX backend test framework. Once we have PRs on the ONNX repo and get
those PRs merged, this file will get EOL'ed.
"""
# pylint: disable=too-many-locals,wrong-import-position,import-error
from __future__ import absolute_import
import sys
import os
import unittest
import logging
import tarfile
from collections import namedtuple
import numpy as np
import numpy.testing as npt
from onnx import numpy_helper, helper, load_model
from onnx import TensorProto
from mxnet.test_utils import download
from mxnet.contrib import onnx as onnx_mxnet
import mxnet as mx
import backend

CURR_PATH = os.path.dirname(os.path.abspath(os.path.expanduser(__file__)))
sys.path.insert(0, os.path.join(CURR_PATH, '../../python/unittest'))

logger = logging.getLogger()
logger.setLevel(logging.DEBUG)

def get_rnd(shape, low=-1.0, high=1.0, dtype=np.float32):
    if dtype == np.float32:
        return (np.random.uniform(low, high,
                                  np.prod(shape)).reshape(shape).astype(np.float32))
    elif dtype == np.int32:
        return (np.random.randint(low, high,
                                  np.prod(shape)).reshape(shape).astype(np.float32))
    elif dtype == np.bool_:
        return np.random.choice(a=[False, True], size=shape).astype(np.float32)

def forward_pass(sym, arg, aux, data_names, input_data):
    """ Perform forward pass on given data
    :param sym: Symbol
    :param arg: Arg params
    :param aux: Aux params
    :param data_names: Input names (list)
    :param input_data: Input data (list). If there is only one input,
                        pass it as a list. For example, if input is [1, 2],
                        pass input_data=[[1, 2]]
    :return: result of forward pass
    """
    data_shapes = []
    data_forward = []
    for idx in range(len(data_names)):
        val = input_data[idx]
        data_shapes.append((data_names[idx], np.shape(val)))
        data_forward.append(mx.nd.array(val))
    # create module
    mod = mx.mod.Module(symbol=sym, data_names=data_names, context=mx.cpu(), label_names=None)
    mod.bind(for_training=False, data_shapes=data_shapes, label_shapes=None)
    if not arg and not aux:
        mod.init_params()
    else:
        mod.set_params(arg_params=arg, aux_params=aux,
                       allow_missing=True, allow_extra=True)
    # run inference
    batch = namedtuple('Batch', ['data'])
    mod.forward(batch(data_forward), is_train=False)

    return mod.get_outputs()[0].asnumpy()

class TestNode(unittest.TestCase):
    """ Tests for models.
    Tests are dynamically added.
    Therefore edit test_models to add more tests.
    """

    def test_import_export(self):
        def get_input_tensors(input_data):
            input_tensor = []
            input_names = []
            input_sym = []
            for idx, ip in enumerate(input_data):
                name = "input" + str(idx + 1)
                input_sym.append(mx.sym.Variable(name))
                input_names.append(name)
                input_tensor.append(helper.make_tensor_value_info(name,
                                                                  TensorProto.FLOAT, shape=np.shape(ip)))
            return input_names, input_tensor, input_sym

        def get_onnx_graph(testname, input_names, inputs, output_name, output_shape, attr):
            outputs = [helper.make_tensor_value_info("output", TensorProto.FLOAT, shape=output_shape)]

            nodes = [helper.make_node(output_name, input_names, ["output"], **attr)]

            graph = helper.make_graph(nodes, testname, inputs, outputs)

            model = helper.make_model(graph)
            return model

        for test in test_cases:
            test_name, mxnet_op, onnx_name, inputs, attrs, mxnet_specific = test
            with self.subTest(test_name):
                names, input_tensors, inputsym = get_input_tensors(inputs)
                test_op = mxnet_op(*inputsym, **attrs)
                mxnet_output = forward_pass(test_op, None, None, names, inputs)
                outputshape = np.shape(mxnet_output)

                if mxnet_specific:
                    onnxmodelfile = onnx_mxnet.export_model(test_op, {}, [np.shape(ip) for ip in inputs],
                                                            np.float32,
                                                            onnx_name + ".onnx")
                    onnxmodel = load_model(onnxmodelfile)
                else:
                    onnxmodel = get_onnx_graph(test_name, names, input_tensors, onnx_name, outputshape, attrs)

                bkd_rep = backend.prepare(onnxmodel, operation='export')
                output = bkd_rep.run(inputs)

                npt.assert_almost_equal(output[0], mxnet_output)

# test_case = ("test_case_name", mxnet op, "ONNX_op_name", [input_list], attribute map, MXNet_specific=True/False)
test_cases = [
    ("test_equal", mx.sym.broadcast_equal, "Equal", [get_rnd((1, 3, 4, 5)), get_rnd((1, 5))], {}, False),
    ("test_greater", mx.sym.broadcast_greater, "Greater", [get_rnd((1, 3, 4, 5)), get_rnd((1, 5))], {}, False),
    ("test_less", mx.sym.broadcast_lesser, "Less", [get_rnd((1, 3, 4, 5)), get_rnd((1, 5))], {}, False),
    ("test_and", mx.sym.broadcast_logical_and, "And",
     [get_rnd((3, 4, 5), dtype=np.bool_), get_rnd((3, 4, 5), dtype=np.bool_)], {}, False),
    ("test_xor", mx.sym.broadcast_logical_xor, "Xor",
     [get_rnd((3, 4, 5), dtype=np.bool_), get_rnd((3, 4, 5), dtype=np.bool_)], {}, False),
    ("test_or", mx.sym.broadcast_logical_or, "Or",
     [get_rnd((3, 4, 5), dtype=np.bool_), get_rnd((3, 4, 5), dtype=np.bool_)], {}, False),
    ("test_not", mx.sym.logical_not, "Not", [get_rnd((3, 4, 5), dtype=np.bool_)], {}, False),
    ("test_square", mx.sym.square, "Pow", [get_rnd((2, 3), dtype=np.int32)], {}, True),
    ("test_spacetodepth", mx.sym.space_to_depth, "SpaceToDepth", [get_rnd((1, 1, 4, 6))],
     {'block_size': 2}, False),
    ("test_softmax", mx.sym.SoftmaxOutput, "Softmax", [get_rnd((1000, 1000)), get_rnd(1000)],
     {'ignore_label': 0, 'use_ignore': False}, True),
    ("test_fullyconnected", mx.sym.FullyConnected, "Gemm", [get_rnd((4,3)), get_rnd((4, 3)), get_rnd(4)],
     {'num_hidden': 4, 'name': 'FC'}, True)
]

if __name__ == '__main__':
    unittest.main()

max_stars_repo_name: ymjiang/incubator-mxnet

max_issues_repo_name: ymjiang/incubator-mxnet

max_forks_repo_name: ymjiang/incubator-mxnet


text:
import itertools
import copy

import numpy as np
import scipy
import scipy.linalg
from functools import reduce

__all__ = ['ImplicitSurface',
           'ImplicitCollection',
           'ImplicitPlane',
           'ImplicitSphere',
           'ImplicitXAlignedCylinder',
           'ImplicitEllipse',
           'ImplicitIntersection',
           'ImplicitUnion',
           'ImplicitDifference',
           'ImplicitComplement',
           'GridMapBase',
           'GridMap',
           'GridSlip'
           ]

class ImplicitSurface(object):

    def __init__(self):
        pass

    def __call__(self):
        raise NotImplementedError('Must be implemented by subclass.')

    def interior(self, grid, asarray=False):
        val = self.__call__(grid)
        if asarray:
            retval = np.zeros_like(val)
            retval[np.where(val < 0.0)] = 1.0
            return retval
        else:
            return np.where(val < 0.0)

    def exterior(self, grid, asarray=False):
        val = self.__call__(grid)
        if asarray:
            retval = np.zeros_like(val)
            retval[np.where(val >= 0.0)] = 1.0
            return retval
        else:
            return np.where(val >= 0.0)

class ImplicitCollection(ImplicitSurface):

    def __init__(self, *items):
        if np.iterable(items[0]):
            self.items = list(items[0])
        else:
            self.items = list(items)

    def __call__(self):
        raise NotImplementedError('Must be implemented by subclass.')

class ImplicitPlane(ImplicitSurface):
    def __init__(self, p, n):
        self.p = np.array(p)
        self.n = np.array(n)
        self.n = self.n/np.linalg.norm(self.n)
        self.d = -np.dot(self.p,self.n)

    def __call__(self, grid):
        return self.d + reduce(lambda x,y:x+y, list(map(lambda x,y:x*y,self.n,grid)))

class ImplicitSphere(ImplicitSurface):
    def __init__(self, c=None, r=1.0):
        if c is None:
            self.c = None
        else:
            self.c = np.array(c)
        self.r = r

    def __call__(self, grid):
        if self.c is None:
            c = np.zeros_like(grid[0].shape)
        else:
            c = self.c
        return reduce(lambda x,y:x+y, list(map(lambda x,y:(y-x)**2,c,grid))) - self.r**2

class ImplicitXAlignedCylinder(ImplicitSurface):
    def __init__(self, c=None, length=1.0, r=1.0):
        if c is None:
            self.c = None
        else:
            self.c = np.array(c)
        self.len = length
        self.r = r

    def __call__(self, grid):
        if self.c is None:
            c = np.zeros_like(grid[0].shape)
        else:
            c = self.c

        g = grid[1:]
        cc = c[1:]
#       longways = (grid[1]-c[1])**2 - self.r**2
        longways =  reduce(lambda x,y:x+y, list(map(lambda x,y:(y-x)**2,cc,g))) - self.r**2
        cutoff   = np.abs(grid[0] - c[0]) -  self.len/2
        return np.maximum(longways, cutoff)

class ImplicitEllipse(ImplicitSurface):
    def __init__(self, c=None, a=None, r=1.0):
        if c is None:
            self.c = None
        else:
            self.c = np.array(c)
        if a is None:
            self.a = None
        else:
            self.a = np.array(a)
        self.r = r

    def __call__(self, grid):
        if self.c is None:
            c = np.zeros(len(grid))
        else:
            c = self.c
        if self.a is None:
            a = np.ones(len(grid))
        else:
            a = self.a
        return reduce(lambda x,y:x+y, list(map(lambda x,y,z:(((y-x)**2)/z), c,grid,a)) ) - self.r**2

class ImplicitIntersection(ImplicitCollection):
    def __init__(self, *items):
        ImplicitCollection.__init__(self, *items)

    def __call__(self, grid):
        return reduce(lambda x,y: np.maximum(x,y), [x(grid) for x in self.items])

class ImplicitUnion(ImplicitCollection):
    def __init__(self, *items):
        ImplicitCollection.__init__(self, *items)

    def __call__(self, grid):
        return reduce(lambda x,y: np.minimum(x,y), [x(grid) for x in self.items])

class ImplicitDifference(ImplicitCollection):

    # Maybe sometime, this should just take *items, and pop the first one for
    # base.  This way, we can allow a single list to be passed.  For now, whatever.
    def __init__(self, base, *items):
        ImplicitCollection.__init__(self, *items)
        self.base = base

    def __call__(self, grid):
        items = [self.base] + self.items
        return reduce(lambda x,y: np.maximum(x,-y), [x(grid) for x in items])

class ImplicitComplement(ImplicitSurface):

    def __init__(self, base):
        self.base = base

    def __call__(self, grid):
        return -1.0*self.base(grid)

    # These must be defined so that the equality is switched if the complement
    # is the calling surface
    def interior(self, grid, asarray=False):
        val = self.__call__(grid)
        if asarray:
            retval = np.zeros_like(val)
            retval[np.where(val <= 0.0)] = 1.0
            return retval
        else:
            return np.where(val <= 0.0)

    def exterior(self, grid, asarray=False):
        val = self.__call__(grid)
        if asarray:
            retval = np.zeros_like(val)
            retval[np.where(val > 0.0)] = 1.0
            return retval
        else:
            return np.where(val > 0.0)

class GridMapBase(object):
    def __init__(self):
        pass

    def __call__(self):
        raise NotImplementedError('Must be implemented by subclass.')

class GridMap(GridMapBase):
    def __init__(self, funcs):
        self.funcs = funcs

    def __call__(self, grid):
        new_grid = []
        for f,g in zip(self.funcs,grid):
            if f is not None:
                new_grid.append(f(g))
            else:
                new_grid.append(g)

        return tuple(new_grid)

class GridSlip(GridMapBase):
    # Creates a slip or a fault along the specifid plane
    def __init__(self, p, n):
        self.p = np.array(p)
        self.n = np.array(n)
        self.n = self.n/np.linalg.norm(self.n)
        self.d = -np.dot(self.p,self.n)

        self.basis = None

    def __call__(self, grid, direction, amount):
        # amount is a scalar
        # direction is a vector that will be normalized
        # the slip occurs along that vector's projection onto the plane, i.e., it is orthogonal to the normal
        # the exterior of the slip plane (the part the normal points to) is the part that is modified.
        d = np.array(direction)
        d = d/np.linalg.norm(d)

        dim = len(d)

        if self.basis is None:
            basis = [self.n]

            # Gramm schmidt
            while len(basis) != dim:
                c = np.random.rand(3)
                for b in basis:
                    c -= np.dot(b,c)*b
                cn = np.linalg.norm(c)
                if cn > 1e-4:
                    basis.append(c/cn)

            self.basis = basis[1:]

        # Project the slip direction onto the basis
        proj=0
        for b in self.basis:
            proj += np.dot(d,b)*b

        # Scale the projection
        proj = amount*proj

        # Evaluate the plane to figure out what components of the old grid are shifted.
        val = self.d + reduce(lambda x,y:x+y, list(map(lambda x,y:x*y,self.n,grid)))
        loc = np.where(val >= 0.0)

        # Perform the shift.
        new_grid = [copy.deepcopy(g) for g in grid]
        for ng,p in zip(new_grid,proj):
            ng[loc] -=p

        return tuple(new_grid)

class Weird(ImplicitSurface):
    def __init__(self, p,n, c=None, r=1.0):
        if c is None:
            self.c = None
        else:
            self.c = np.array(c)
        self.r = r

        self.p = np.array(p)
        self.n = np.array(n)
        self.n = self.n/np.linalg.norm(self.n)
        self.d = -np.dot(self.p,self.n)

    def __call__(self, grid):
        if self.c is None:
            c = np.zeros_like(grid[0].shape)
        else:
            c = self.c

# Creates flat eye thingies
#       val1 = reduce(lambda x,y:x+y, map(lambda x,y:(y-x)**2,c,grid)) - self.r**2
#
#       val2 = self.d + reduce(lambda x,y:x+y, map(lambda x,y:x*y,self.n,grid))
#
#       return 4*val1/np.max(val1)+val2

        val1 = reduce(lambda x,y:x+y, list(map(lambda x,y:(y-x)**2,c,grid))) - self.r**2

        val2 = self.d + reduce(lambda x,y:x+y, list(map(lambda x,y:x*y,self.n,grid)))

        return 4*val1/np.max(val1)+val2

class Hyperbola(ImplicitSurface):
    def __init__(self, c=None, r=1.0):
        if c is None:
            self.c = None
        else:
            self.c = np.array(c)
        self.r = r

    def __call__(self, grid):
        if self.c is None:
            c = np.zeros_like(grid[0].shape)
        else:
            c = self.c

        return  reduce(lambda x,y:-x+y, list(map(lambda x,y:(y-x)**2,c,grid))) - self.r**2

class Weird2(ImplicitSurface):
    def __init__(self, c=None, r=1.0, s=None):
        if c is None:
            self.c = None
        else:
            self.c = np.array(c)

        if s is None:
            self.s = None
        else:
            self.s = np.array(s)

        self.r = r

    def __call__(self, grid):
        c = np.zeros_like(grid[0].shape) if self.c is None else self.c
        s = np.ones_like(grid[0].shape) if self.s is None else self.s

        return ((grid[0]-c[0])**2)/s[0] + ((grid[1]-c[1])**1)/s[1] - self.r**2
        #reduce(lambda x,y:-x+y, map(lambda x,y:(y-x)**2,c,grid)) - self.r**2

# FIXME
#if __name__ == "__main__":
#
#   import numpy as np
#   import matplotlib.pyplot as mpl
#
#   from pysit import *
#
#   x_lbc = PML(0.0, 100.0)
#   x_rbc = PML(0.0, 100.0)
#   z_lbc = PML(0.0, 100.0)
#   z_rbc = PML(0.0, 100.0)
#
#   xmin, xmax = 0.0, 20
#   zmin, zmax = 0.0, 10
#   nx, nz = 500,250
#
#   x_config = (xmin, xmax, nx, x_lbc, x_rbc)
#   z_config = (zmin, zmax, nz, z_lbc, z_rbc)
#
#   d = Domain((x_config, z_config))
#
#   grid = d.generate_grid()
#   grid = grid[0].reshape(d.shape).T, grid[1].reshape(d.shape).T
#
#   p1 = ImplicitPlane((0.0,4.0),(0.0,-1.0))
#   p2 = ImplicitPlane((0.0,6.0),(0.0,1.0))
#
#   w1 = Weird((0.0,5.0),(0.0,-1.0), (10,5))
#   w2 = Weird((0.0,6.0),(0.0,1.0), (10,5))
#
#   w3 = Weird2((10,4),1.5,(-100.0,1.))
#   w4 = Weird2((10,4),0.5,( 100.0,1.))
#
#   band = ImplicitIntersection(w1,w2)
#   band = ImplicitIntersection(p1,p2)
#   band = ImplicitDifference(w3,w4)
#
#   g = GridSlip((10.0,5.0), (1,0.5))
#   d = (0.0,-1.0)
#   a = 0.5
#   new_grid = g(grid,d,a)
#
#
#   plt.figure()
#   plt.imshow(band(grid))
#   plt.colorbar()
#
#   plt.figure()
#   plt.imshow(band(new_grid))
#   plt.colorbar()
#
#   plt.figure()
#   plt.imshow(band.interior(grid))
#   plt.colorbar()
#
#   plt.figure()
#   plt.imshow(band.interior(new_grid))
#   plt.colorbar()
#
#   plt.show()

ext: py

max_stars_repo_name: zfang-slim/pysit

max_issues_repo_name: zfang-slim/pysit

max_forks_repo_name: zfang-slim/pysit

reason: import numpy,import scipy


text:
import os

import numpy
import moderngl
import moderngl_window

from moderngl_window.conf import settings
from moderngl_window.resources import programs
from moderngl_window.meta import ProgramDescription

settings.PROGRAM_DIRS.append(os.path.join(os.path.dirname(__file__), 'programs'))

class MeshProgram:
    """
    Describes how a mesh is rendered using a specific shader program
    """

    def __init__(self, program: moderngl.Program = None, **kwargs):
        """Initialize.

        Args:
            program: The moderngl program
        """
        self.program = program

    @property
    def ctx(self):
        """moderngl.Context: The current context"""
        return moderngl_window.ctx()

    def draw(self, mesh, projection_matrix: numpy.ndarray = None,
             model_matrix: numpy.ndarray = None, camera_matrix: numpy.ndarray = None, time=0.0):
        """Draw code for the mesh

        Args:
            mesh (Mesh): The mesh to render
        Keyword Args:
            projection_matrix (numpy.ndarray): projection_matrix (bytes)
            model_matrix (numpy.ndarray): view_matrix (bytes)
            camera_matrix (numpy.ndarray): camera_matrix (bytes)
            time (float): The current time
        """
        self.program["m_proj"].write(projection_matrix)
        self.program["m_mv"].write(model_matrix)
        mesh.vao.render(self.program)

    def apply(self, mesh):
        """
        Determine if this ``MeshProgram`` should be applied to the mesh.
        Can return self or some ``MeshProgram`` instance to support dynamic ``MeshProgram`` creation

        Args:
            mesh: The mesh to inspect
        """
        raise NotImplementedError("apply is not implemented. Please override the MeshProgram method")

class ColorProgram(MeshProgram):
    """
    Simple color program
    """
    def __init__(self, program=None, **kwargs):
        super().__init__(program=None)
        self.program = programs.load(ProgramDescription(path="scene_default/color.glsl"))

    def draw(self, mesh, projection_matrix=None, model_matrix=None, camera_matrix=None, time=0):

        if mesh.material:
            # if mesh.material.double_sided:
            #     self.ctx.disable(moderngl.CULL_FACE)
            # else:
            #     self.ctx.enable(moderngl.CULL_FACE)

            if mesh.material.color:
                self.program["color"].value = tuple(mesh.material.color)
            else:
                self.program["color"].value = (1.0, 1.0, 1.0, 1.0)

        self.program["m_proj"].write(projection_matrix)
        self.program["m_model"].write(model_matrix)
        self.program["m_cam"].write(camera_matrix)
        mesh.vao.render(self.program)

    def apply(self, mesh):
        if not mesh.material:
            return None

        if not mesh.attributes.get("NORMAL"):
            return None

        if mesh.material.mat_texture is None:
            return self

        return None

class TextureProgram(MeshProgram):
    """
    Simple texture program
    """
    def __init__(self, program=None, **kwargs):
        super().__init__(program=None)
        self.program = programs.load(ProgramDescription(path="scene_default/texture.glsl"))

    def draw(self, mesh, projection_matrix=None, model_matrix=None, camera_matrix=None, time=0):
        # if mesh.material.double_sided:
        #     self.ctx.disable(moderngl.CULL_FACE)
        # else:
        #     self.ctx.enable(moderngl.CULL_FACE)

        mesh.material.mat_texture.texture.use()
        self.program["texture0"].value = 0
        self.program["m_proj"].write(projection_matrix)
        self.program["m_model"].write(model_matrix)
        self.program["m_cam"].write(camera_matrix)
        mesh.vao.render(self.program)

    def apply(self, mesh):
        if not mesh.material:
            return None

        if not mesh.attributes.get("NORMAL"):
            return None

        if mesh.material.mat_texture is not None:
            return self

        return None

class FallbackProgram(MeshProgram):
    """
    Fallback program only rendering positions in white
    """
    def __init__(self, program=None, **kwargs):
        super().__init__(program=None)
        self.program = programs.load(ProgramDescription(path="scene_default/fallback.glsl"))

    def draw(self, mesh, projection_matrix=None, model_matrix=None, camera_matrix=None, time=0):

        self.program["m_proj"].write(projection_matrix)
        self.program["m_model"].write(model_matrix)
        self.program["m_cam"].write(camera_matrix)

        if mesh.material:
            self.program["color"].value = tuple(mesh.material.color[0:3])
        else:
            self.program["color"].value = (1.0, 1.0, 1.0)

        mesh.vao.render(self.program)

    def apply(self, mesh):
        return self

max_stars_repo_name: minuJeong/moderngl-window

max_issues_repo_name: minuJeong/moderngl-window

max_forks_repo_name: minuJeong/moderngl-window

reason: import numpy


text:
from pso.GPSO import GPSO
import numpy as np
import time
import pandas as pd

np.random.seed(42)

# f1 完成
def Sphere(p):
    # Sphere函数
    out_put = 0
    for i in p:
        out_put += i ** 2
    return out_put
# f2 完成
def Sch222(x):
    out_put = 0
    out_put01 = 1
    for i in x:
        out_put += abs(i)
        out_put01 = abs(i)*out_put01
    out_put = out_put01+out_put
    return out_put

# f3 完成
def Quadric(x):
    output = 0
    # print(x.shape[0])
    for i in range(x.shape[0]):
        output += (np.sum(x[0:i+1])) ** 2
        # print(np.square(np.sum(x[0:i+1])))
    return output

# f4 完成
def Schl(x):
    # print(np.max(np.abs(x)))
    return np.max(np.abs(x))

# f5 完成
def Step(x):
    output = 0
    for i in x:
        output += (np.floor(i+0.5))**2
    return output

# f6 完成
def Noise(x):
    output = 0
    cnt = 1
    for i in x:
        output = cnt * (i**4) + output
        cnt += 1
    output += np.random.rand()
    return output
# f7 完成
def Rosenbrock(p):
    '''
    -2.048<=xi<=2.048
    函数全局最优点在一个平滑、狭长的抛物线山谷内，使算法很难辨别搜索方向，查找最优也变得十分困难
    在(1,...,1)处可以找到极小值0
    :param p:
    :return:
    '''
    n_dim = len(p)
    res = 0
    for i in range(n_dim - 1):
        res += 100 * np.square(np.square(p[i]) - p[i + 1]) + np.square(p[i] - 1)
    return res

# f8 有问题，忽略,这个是APSO的f8
def Schewel(x):
    out_put = 0
    for i in x:
        out_put += -i*np.sin(np.sqrt(abs(i)))
    return out_put

# f9 完成
def Rastrigin(p):
    '''
    多峰值函数，也是典型的非线性多模态函数
    -5.12<=xi<=5.12
    在范围内有10n个局部最小值，峰形高低起伏不定跳跃。很难找到全局最优
    has a global minimum at x = 0  where f(x) = 0
    '''
    return np.sum([np.square(x) - 10 * np.cos(2 * np.pi * x) + 10 for x in p])

# f10
def Ackley(x):
    part1 = 0
    part2 = 0
    for i in x:
        part1 += (i**2)
        part2 += np.cos(2 * np.pi * i)
    left = 20 * np.exp(-0.2 * ((part1 / x.shape[0]) ** .5))
    right = np.exp(part2 / x.shape[0])
    return  -left - right + 20 + np.e

# f11  ok
def Griewank(p):
    '''
    存在多个局部最小值点，数目与问题的维度有关。
    此函数是典型的非线性多模态函数，具有广泛的搜索空间，是优化算法很难处理的复杂多模态问题。
    在(0,...,0)处取的全局最小值0
    -600<=xi<=600
    '''
    part1 = [np.square(x) / 4000 for x in p]
    part2 = [np.cos(x / np.sqrt(i + 1)) for i, x in enumerate(p)]
    return np.sum(part1) - np.prod(part2) + 1

g = 10000
times = 30
table = np.zeros((2, 10))
gBest = np.zeros((10, 30))  # 1010个函数的30次的最优值
for i in range(times):
    optimizer = GPSO(func=Sphere, dim=30, pop=20, max_iter=g, lb=np.ones(30) * (-100), ub=np.ones(30) * 100,
                     w=0.9, c1=2, c2=2, acceptance=0.01)
    start = time.time()
    optimizer.run()
    end = time.time()
    print('Sphere:', optimizer.gbest_y)
    table[0, 0] += optimizer.gbest_y
    table[1, 0] += end - start
    gBest[0, i] = optimizer.gbest_y

    optimizer = GPSO(func=Sch222, dim=30, pop=20, max_iter=g, lb=np.ones(30) * (-10), ub=np.ones(30) * 10,
                     w=0.9, c1=2, c2=2, acceptance=0.01)
    start = time.time()
    optimizer.run()
    end = time.time()
    print('Sch222:', optimizer.gbest_y)
    table[0, 1] += optimizer.gbest_y
    table[1, 1] += end - start
    gBest[1, i] = optimizer.gbest_y

    optimizer = GPSO(func=Quadric, dim=30, pop=20, max_iter=g, lb=np.ones(30) * (-100), ub=np.ones(30) * 100,
                     w=0.9, c1=2, c2=2, acceptance=100)
    start = time.time()
    optimizer.run()
    end = time.time()
    print('Quadric:', optimizer.gbest_y)
    table[0, 2] += optimizer.gbest_y
    table[1, 2] += end - start
    gBest[2, i] = optimizer.gbest_y

    optimizer = GPSO(func=Rosenbrock, dim=30, pop=20, max_iter=g, lb=np.ones(30) * (-10), ub=np.ones(30) * 10,
                     w=0.9, c1=2, c2=2, acceptance=100)
    start = time.time()
    optimizer.run()
    end = time.time()
    print('Rosenbrock:', optimizer.gbest_y)
    table[0, 3] += optimizer.gbest_y
    table[1, 3] += end - start
    gBest[3, i] = optimizer.gbest_y

    optimizer = GPSO(func=Step, dim=30, pop=20, max_iter=g, lb=np.ones(30) * (-100), ub=np.ones(30) * 100,
                     w=0.9, c1=2, c2=2, acceptance=0)
    start = time.time()
    optimizer.run()
    end = time.time()
    print('Step:', optimizer.gbest_y)
    table[0, 4] += optimizer.gbest_y
    table[1, 4] += end - start
    gBest[4, i] = optimizer.gbest_y

    optimizer = GPSO(func=Noise, dim=30, pop=20, max_iter=g, lb=np.ones(30) * (-1.28), ub=np.ones(30) * 1.28,
                     w=0.9, c1=2, c2=2, acceptance=0.01)
    start = time.time()
    optimizer.run()
    end = time.time()
    print('Noise:', optimizer.gbest_y)
    table[0, 5] += optimizer.gbest_y
    table[1, 5] += end - start
    gBest[5, i] = optimizer.gbest_y

    optimizer = GPSO(func=Schewel, dim=30, pop=20, max_iter=g, lb=np.ones(30) * (-500), ub=np.ones(30) * 500,
                     w=0.9, c1=2, c2=2, acceptance=-10000)
    start = time.time()
    optimizer.run()
    end = time.time()
    print('Schewel:', optimizer.gbest_y)
    table[0, 6] += optimizer.gbest_y
    table[1, 6] += end - start
    gBest[6, i] = optimizer.gbest_y

    optimizer = GPSO(func=Rastrigin, dim=30, pop=20, max_iter=g, lb=np.ones(30) * (-5.12), ub=np.ones(30) * 5.12,
                     w=0.9, c1=2, c2=2, acceptance=50)
    start = time.time()
    optimizer.run()
    end = time.time()
    print('Rastrigin:', optimizer.gbest_y)
    table[0, 7] += optimizer.gbest_y
    table[1, 7] += end - start
    gBest[7, i] = optimizer.gbest_y

    optimizer = GPSO(func=Ackley, dim=30, pop=20, max_iter=g, lb=np.ones(30) * (-32), ub=np.ones(30) * 32,
                     w=0.9, c1=2, c2=2, acceptance=0.01)
    start = time.time()
    optimizer.run()
    end = time.time()
    print('Ackley:', optimizer.gbest_y)
    table[0, 8] += optimizer.gbest_y
    table[1, 8] += end - start
    gBest[8, i] = optimizer.gbest_y

    optimizer = GPSO(func=Griewank, dim=30, pop=20, max_iter=g, lb=np.ones(30) * (-600), ub=np.ones(30) * 600,
                     w=0.9, c1=2, c2=2, acceptance=0.01)
    start = time.time()
    optimizer.run()
    end = time.time()
    print('Griewank:', optimizer.gbest_y)
    table[0, 9] += optimizer.gbest_y
    table[1, 9] += end - start
    gBest[9, i] = optimizer.gbest_y

table = table / times
table = pd.DataFrame(table)
table.columns = ['Sphere', 'Schwefel_P222', 'Quadric', 'Rosenbrock', 'Step', 'Quadric_Noise', 'Schwefel',
                 'Rastrigin', 'Ackley', 'Griewank']
table.index = ['mean score', 'mean time']
print(table)
print('10个测试函数的30次std:', np.std(gBest, axis=1))
print('10个测试函数的30次best:', np.min(gBest, axis=1))

max_stars_repo_name: gmjustforfun/code

max_issues_repo_name: gmjustforfun/code

max_forks_repo_name: gmjustforfun/code


text:
import math
import sys
import random
import numpy as np
import collections

# sys.stdin = open("testdata.in", "r")
sys.stdout = open("testdata.in", "w+")

N = 100000
print(N)
for i in range(1, N):
    j = random.randint(1, 100000)
    print(j, end=' ')


text:
#OS LIBRARIES
import csv
import sys
import os
# Numerical libraries
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.cluster import KMeans
# to handle data in form of rows and columns
import pandas as pd
# importing ploting libraries
#from matplotlib import pyplot as plt
import matplotlib.pyplot as plt
import folium
#%matplotlib inline
#importing seaborn for statistical plots
import seaborn as sns
from sklearn import metrics
from matplotlib import style
import matplotlib.ticker as ticker
import seaborn as sns
# Configuración matplotlib
# ==============================================================================
plt.rcParams['image.cmap'] = "bwr"
#plt.rcParams['figure.dpi'] = "100"
plt.rcParams['savefig.bbox'] = "tight"
style.use('ggplot') or plt.style.use('ggplot')

# Configuración warnings
# ==============================================================================
import warnings
warnings.filterwarnings('ignore')

#FILES
strFile = "C:/Users/Helen Jurado/Documents/GitHub/SIAMS/public/img/img_arules/codo_jambu.png"
strFile1 = "C:/Users/Helen Jurado/Documents/GitHub/SIAMS/public/img/img_arules/correlacion.png"
strFile2 = "C:/Users/Helen Jurado/Documents/GitHub/SIAMS/public/img/img_arules/distribucion.png"
strFile3 = "C:/Users/Helen Jurado/Documents/GitHub/SIAMS/public/img/img_arules/density.png"
#ENTRADAS
arch_inicial= sys.argv[1]
cluster=sys.argv[2]
#file_url = 'C:/Users/Helen Jurado/Documents/GitHub/SIAMS/storage/archivos_apriori/coordenadas-120921_21_09_38.csv'
#data = pd.read_csv(file_url)
data = pd.read_csv(arch_inicial)
n_cltrs=int(cluster)
#n_cltrs=4
data.info()
data.describe()
features = data[['latitud', 'longitud']]
print(features)
distribution = data[['orden', 'id_trayectoria', 'velocidad']]
wcss = []
for i in range(1, n_cltrs+2):
    kmeans = KMeans(n_clusters = i, max_iter = 300)
    kmeans.fit(features)
    wcss.append(kmeans.inertia_)
plt.plot(range(1, n_cltrs+2,1), wcss,marker='*')
plt.title("Codo de Jambú")
plt.xlabel('Número de clusters')
plt.ylabel('WCSS')
#strFile = "C:/Users/Helen Jurado/Documents/GitHub/SIAMS/public/img/img_arules/codo_jambu.png"
if os.path.isfile(strFile):
   os.remove(strFile)
plt.savefig(strFile)#'C:/Users/Helen Jurado/Documents/GitHub/SIAMS/public/img/img_arules/codo_jambu.png')
#plt.show()
# — — — — — — — — — — — — — — — -Heat map to identify highly correlated variables — — — — — — — — — — — — -
#-------------------------------Heat map to identify highly correlated variables-------------------------
plt.figure(figsize=(10,8))
sns.heatmap(features.corr(),
            annot=True,
            linewidths=.5,
            center=0,
            cbar=False,
            cmap="YlGnBu")
#strFile1 = "C:/Users/Helen Jurado/Documents/GitHub/SIAMS/public/img/img_arules/correlacion.png"
if os.path.isfile(strFile1):
   os.remove(strFile1)
plt.savefig(strFile1)#'C:/Users/Helen Jurado/Documents/GitHub/SIAMS/public/img/img_arules/correlacion.png')
#plt.show()
mydata = data
mydata.drop(columns = {'id_det_tra','id_trayectoria','orden','longitud','latitud','fecha','coordenadas','tipo_coordenada'}, inplace=True)
data = features
#--Checking Outliers
#plt.figure(figsize=(15,10))
#pos = 1
#for i in mydata.columns:
#    plt.subplot(3, 3, pos)
#    sns.boxplot(data=mydata[i])
#    pos += 1
fig, axes = plt.subplots(nrows=3, ncols=1, figsize=(6, 6))
sns.distplot(
    distribution.id_trayectoria,
    hist    = False,
    rug     = True,
    color   = "blue",
    kde_kws = {'shade': True, 'linewidth': 1},
    ax      = axes[0]
)
axes[0].set_title("Distribución original", fontsize = 'medium')
axes[0].set_xlabel('Trayectorias', fontsize='small')
axes[0].tick_params(labelsize = 6)

sns.distplot(
    np.sqrt(distribution.id_trayectoria),
    hist    = False,
    rug     = True,
    color   = "blue",
    kde_kws = {'shade': True, 'linewidth': 1},
    ax      = axes[1]
)
axes[1].set_title("Transformación raíz cuadrada", fontsize = 'medium')
axes[1].set_xlabel('sqrt(Trayectorias)', fontsize='small')
axes[1].tick_params(labelsize = 6)

sns.distplot(
    np.log(distribution.id_trayectoria),
    hist    = False,
    rug     = True,
    color   = "blue",
    kde_kws = {'shade': True, 'linewidth': 1},
    ax      = axes[2]
)
axes[2].set_title("Transformación logarítmica", fontsize = 'medium')
axes[2].set_xlabel('log(Trayectorias)', fontsize='small')
axes[2].tick_params(labelsize = 6)
fig.tight_layout()
#strFile3 = "C:/Users/Helen Jurado/Documents/GitHub/SIAMS/public/img/img_arules/density.png"
if os.path.isfile(strFile3):
   os.remove(strFile3)
plt.savefig(strFile3)#'C:/Users/Helen Jurado/Documents/GitHub/SIAMS/public/img/img_arules/density.png')
##Scale the data
from scipy.stats import zscore
mydata_z = data.apply(zscore)
mydata_z.head()
from sklearn.cluster import KMeans
# create kmeans model/object
kmeans = KMeans(
    init="random",
    n_clusters=n_cltrs,
    n_init=10,
    max_iter=300,
    random_state=42
)
# do clustering
kmeans.fit(features)
# save results
labels = kmeans.labels_
# send back into dataframe and display it

data['cluster'] = labels
# display the number of mamber each clustering
_clusters = data.groupby('cluster')['latitud'].count()
print(_clusters)
data
# List to store cluster and intra cluster distance
clusters = []
inertia_vals = []
# Since creating one cluster is similar to observing the data as a whole, multiple values of K are utilized to come up with the optimum cluster value
#Note: Cluster number and intra cluster distance is appended for plotting the elbow curve
for k in range(1, n_cltrs, 1):

    # train clustering with the specified K
    model = KMeans(n_clusters=k, random_state=7, n_jobs=10)
    model.fit(mydata_z)
# append model to cluster list
    clusters.append(model)
    inertia_vals.append(model.inertia_)
p_data = data
print(p_data)
df = pd.DataFrame(data = p_data, columns = ['latitud', 'longitud', 'cluster'])
df.head()
colors = ['red', 'blue', 'green', 'purple', 'orange', 'darkred', \
     'beige', 'darkblue', 'darkgreen', 'cadetblue', \
     'pink', 'lightblue', 'lightgreen', 'gray', \
     'black', 'lightgray', 'red', 'blue', 'green', 'purple', \
     'orange', 'darkred', 'lightred', 'beige', 'darkblue', \
     'darkgreen', 'cadetblue', 'darkpurple','pink', 'lightblue', \
     'lightgreen', 'gray', 'black', 'lightgray' ]

fig = plt.figure(figsize=(7,5))
ax = fig.add_subplot(1,1,1)
ax.set_xlabel('Latitud', fontsize = 15)
ax.set_ylabel('Longitud', fontsize = 15)
ax.set_title('SIAMS UG', fontsize = 20)
color_theme = np.array(colors)
ax.scatter(df['latitud'],df['longitud'], c=color_theme[df.cluster], s =50)
#strFile2 = "C:/Users/Helen Jurado/Documents/GitHub/SIAMS/public/img/img_arules/distribucion.png"
if os.path.isfile(strFile2):
   os.remove(strFile2)
plt.savefig(strFile2)#'C:/Users/Helen Jurado/Documents/GitHub/SIAMS/public/img/img_arules/distribucion.png')
plt.show()

wcss = []
for i in range(1, n_cltrs+2):
    kmeans = KMeans(n_clusters = i, max_iter = 300)
    kmeans.fit(features)
    wcss.append(kmeans.inertia_)

lat = data.iloc[0]['latitud']
lng = data.iloc[0]['longitud']
map = folium.Map(location=[lng, lat], zoom_start=12)
for _, row in data.iterrows():
    folium.CircleMarker(
        location=[row["latitud"], row["longitud"]],
        radius=12,
        weight=2,
        fill=True,
        fill_color=colors[int(row["cluster"])],
        color=colors[int(row["cluster"])]
    ).add_to(map)
map

max_stars_repo_name: XavierJuradoM/SIAMS

max_issues_repo_name: XavierJuradoM/SIAMS

max_forks_repo_name: XavierJuradoM/SIAMS

reason: import numpy,from scipy


text:
from numpy.linalg import norm
from scipy.stats import norm

def SmoothStep(x,h):
    # This function computes the step function smoothed via Gaussian cdf
    # INPUTS
    #   x         : [matrix] (n_ x n_) grid points
    #   h         : [scalar] smoothing parameter
    # OUTPUTS
    #   y         : [matrix] (n_ x n_) smoothed values ad grid points

    ## Code
    y = norm.cdf(x,0,h)
    return y

max_stars_repo_name: dpopadic/arpmRes

max_issues_repo_name: dpopadic/arpmRes

max_forks_repo_name: dpopadic/arpmRes

reason: from numpy,from scipy


text:
# Arda Mavi

import os
import numpy as np
from os import listdir
from random import shuffle
from get_dataset import read_npy_dataset, split_npy_dataset
from keras.callbacks import ModelCheckpoint, TensorBoard
from get_models import get_segment_model, get_Discriminator, get_GAN, get_Generator, save_model

epochs = 100
batch_size = 4
test_size = 0.2

# Dice Coefficient Loss Function
def dice_coefficient(y_true, y_pred):
    smoothing_factor = 1
    flat_y_true = y_true.flatten()
    flat_y_pred = y_pred.flatten()
    return (2. * np.sum(flat_y_true * flat_y_pred) + smoothing_factor) / (np.sum(flat_y_true) + np.sum(flat_y_pred) + smoothing_factor)

# Training Segment Model:
def train_seg_model(model, splitted_npy_dataset_path, test_path, epochs):
    test_XY = np.load(test_path+'/test.npy')
    X_test, Y_test = test_XY[0], test_XY[1]

    batch_dirs = listdir(splitted_npy_dataset_path)
    len_batch_dirs = len(batch_dirs)

    if not os.path.exists('Data/Checkpoints/'):
        os.makedirs('Data/Checkpoints/')
    checkpoints = []
    checkpoints.append(ModelCheckpoint('Data/Checkpoints/best_weights.h5', monitor='val_loss', verbose=0, save_best_only=False, save_weights_only=True, mode='auto', period=1))
    checkpoints.append(TensorBoard(log_dir='Data/Checkpoints/./logs', histogram_freq=0, write_graph=True, write_images=False, embeddings_freq=0, embeddings_layer_names=None, embeddings_metadata=None))

    for epoch in range(epochs):
        print('Epoch: {0}/{1}'.format(epoch+1, epochs))
        model.fit_generator(data_gen(splitted_npy_dataset_path), steps_per_epoch=batch_size, epochs=int(len_batch_dirs/batch_size), callbacks=checkpoints)

        scores = model.evaluate(X_test, Y_test)
        dice_score = dice_coefficient(model.predict(X_test), Y_test)
        print('Test loss:', scores[0], '\nTest accuracy:', scores[1], '\nDice Coefficient Accuracy:', dice_score)
    return model

# Training GAN:
def train_gan(Generator, Encoder, Discriminator, GAN, splitted_npy_dataset_path, test_path, epochs):
    # Dataset sample sharing:
    batch_dirs = listdir(splitted_npy_dataset_path)
    batch_num = len(batch_dirs)

    test_XY = np.load(test_path+'/test.npy')
    X_test, Y_test = test_XY[0], test_XY[1]

    data_generator = data_gen(splitted_npy_dataset_path)

    for epoch in range(epochs):
        print('Epoch: {0}/{1}'.format(epoch+1, epochs))

        print('Discriminator Training:')
        # Discriminator Training:
        for i in range(int(batch_num/4)):
            # Getting generated data:
            getted_generator_data = next(data_generator)
            generator_pred = Generator.predict(getted_generator_data[0])

            # Getting real data:
            getted_real_data = next(data_generator)

            # Getting batch:
            # Xs:
            dis_batch_X1 = np.concatenate((getted_generator_data[0], getted_real_data[0]), axis=0)
            dis_batch_X2 = np.concatenate((generator_pred, getted_real_data[1]), axis=0)
            # Ys:
            dis_batch_Y = np.concatenate((np.zeros((getted_generator_data[0].shape[0], 1)), np.ones((getted_real_data[0].shape[0], 1))), axis=0)

            Discriminator.trainable = True
            Discriminator.fit([dis_batch_X1, dis_batch_X2], dis_batch_Y, batch_size=1, epochs=1, shuffle=True)

        print('Generator Training:')
        # Generator Training:
        for i in range(int(batch_num/2)):
            # Getting generated data:
            getted_generator_data = next(data_generator)

            # Getting batch:
            # Xs:
            gan_batch_X = getted_generator_data[0]
            # Ys:
            gan_batch_Y = np.ones((gan_batch_X.shape[0], 1))

            Discriminator.trainable = False
            Encoder.trainable = True
            GAN.fit(gan_batch_X, gan_batch_Y, batch_size=1, epochs=1, shuffle=True)

        scores = Generator.evaluate(X_test, Y_test)
        dice_score = dice_coefficient(Generator.predict(X_test), Y_test)
        print('Test Loss:', scores[0], '\nTest Accuracy:', scores[1], '\nDice Coefficient Accuracy:', dice_score)
        save_model(Generator, path='Data/Checkpoints/GAN-Models/Generator/', model_name = 'model', weights_name = 'weights')
        print('Segmentation model checkpoints saved to "Data/Chackpoints/GAN-Models/Generator/"')

    return Generator, Encoder, Discriminator

def data_gen(splitted_npy_dataset_path):
    batch_dirs = listdir(splitted_npy_dataset_path)
    while True:
        shuffle(batch_dirs)
        for batch_path in batch_dirs:
            batch_XY = np.load(splitted_npy_dataset_path+'/'+batch_path)
            X_batch, Y_batch = batch_XY[0], batch_XY[1]
            yield X_batch, Y_batch

def main(train_gan_model = True):
    if train_gan_model:
        # Getting Generator:
        Generator, Encoder = get_segment_model(data_shape = (128, 128, 128, 1))
        Discriminator = get_Discriminator(input_shape_1 = (128,128,128,1), input_shape_2 = (128, 128, 128, 1), Encoder = Encoder)
        GAN = get_GAN((128, 128, 128, 1), Generator, Discriminator)

        # Saving non-trained models:
        save_model(Generator, path='Data/GAN-Models/Generator/', model_name = 'model', weights_name = 'weights')
        save_model(Encoder, path='Data/GAN-Models/Encoder/', model_name = 'model', weights_name = 'weights')
        save_model(Discriminator, path='Data/GAN-Models/Discriminator/', model_name = 'model', weights_name = 'weights')
        print('Non-Trained model saved to "Data/GAN-Models"!')

        # Train:
        Generator, Encoder, Discriminator = train_gan(Generator, Encoder, Discriminator, GAN, splitted_npy_dataset_path='Data/npy_dataset/splitted_npy_dataset', test_path = 'Data/npy_dataset/test_npy', epochs = epochs)

        # Saving trained models:
        save_model(Generator, path='Data/GAN-Models/Generator/', model_name = 'model', weights_name = 'weights')
        save_model(Encoder, path='Data/GAN-Models/Encoder/', model_name = 'model', weights_name = 'weights')
        save_model(Discriminator, path='Data/GAN-Models/Discriminator/', model_name = 'model', weights_name = 'weights')
        print('Trained model saved to "Data/GAN-Models"!')
        return Generator
    else:
        segment_model, _ = get_segment_model(data_shape = (128, 128, 128, 1))
        save_model(segment_model, path='Data/Model/', model_name = 'model', weights_name = 'weights')
        print('Non-Trained model saved to "Data/Model"!')
        model = train_seg_model(segment_model, splitted_npy_dataset_path='Data/npy_dataset/splitted_npy_dataset', test_path = 'Data/npy_dataset/test_npy', epochs = epochs)
        save_model(segment_model, path='Data/Model/', model_name = 'model', weights_name = 'weights')
        print('Trained model saved to "Data/Model"!')
        return segment_model

if __name__ == '__main__':
    main(train_gan_model = False)


text:
#ce fichier regroupe des fonctions absolument pas pensees pour etre generales
#mais utilisees au cours du post-doc pour exploiter les calculs CFX
try: from paraview import vtk
except: import vtk

import numpy

from fonctions_basiques import *
from calculs import *
from ecriture import *

#__________________________________________________________________________________________
def lire_fichier_out_cfx(acces_fichier):
    """Fonction extremement speciale.
    Lecture d'un fichier .out CFX et va ensuite rajouter omega aux blocs
    identitifes par leurs noms de domaine

    """
    #Lecture des omega par noms de domaine dans le fichier .out
    f = file(acces_fichier)
    info = f.read()
    info = info[:info.rfind('END')]

    #lecture des parametres
    parametres = info[info.find('EXPRESSIONS:'):]
    parametres = parametres[:parametres.find('END')]
    parametres = parametres.split('\n')[1:]
    for param in parametres:
        if '[' in param and ']' in param:
            param = param[:param.find('[')]
            param = param.strip()
            exec param

    #creation du dictionnaire des informations, pour chacun des domaines et chacunes des frontieres
    domaines = info.split('DOMAIN: ')[1:]
    for k in domaines:
        k = k[:k.rfind('END')]

    dict_domaines = {}
    for dom in domaines:
        nom = dom[:dom.find('\n')]

    return omega
#__________________________________________________________________________________________

#__________________________________________________________________________________________
def reconstruire_serie_fourier(input, variable, time, periode, dephasage=0, suppression_coeff_fourier=True, liste_rangs=None):
    """fonction qui reconstruit une serie de Fourier
    Les coefficient doivent etre stockes au points.
    Le nom de la variable doit contenir <variable>
    Le coefficient est determine par la fin du nom de la variable
        A0, A1, ..., A15, ... , B0, ..., B17, ...

    <input> est un objet VTK, qui peut etre un multibloc

    la phase est calculee avec
        phase = k * 2 * numpy.pi * (time + dephasage) / periode
    ou k est le range de l'harmonique qui est en train d'etre calculee
    """
    ####################################
    # cas multibloc
    if isinstance(input, vtk.vtkMultiBlockDataSet):
        output = vtk_new_instance(input)
        for numbloc in get_numeros_blocs_non_vides(input):
            output.SetBlock(numbloc, reconstruire_serie_fourier(input.GetBlock(numbloc), variable, time, periode, dephasage, suppression_coeff_fourier, liste_rangs))

    ####################################
    # cas monobloc - on etablit la liste des coefficients de Fourier presents
    else:
        # on se base sur le coefficient A0
        base_variable_name = None
        for var in get_noms_arrays_presents(input, 'points'):
            if var[-2:] == 'A0' and variable in var:
                base_variable_name = var[:-2]
                break
        if base_variable_name is None:
            raise Exception, "n'a pas reussi a determiner le nom generique de la variable..."

        # on cherche tous les coefficients disponibles
        liste_coeff = []
        for var in get_noms_arrays_presents(input, 'points'):
            if base_variable_name in var:
                liste_coeff.append(var)

        # on ecrit la formule
        formule = None
        for coeff in liste_coeff:
            rang = coeff.split(base_variable_name)[1]
            if rang[0] == 'A':
                op = 'cos'
            elif rang[0] == 'B':
                op = 'sin'
            else:
                raise Exception, 'Pas reussi a determiner a quoi correspond le coeff {0}'.format(coeff)

            try:
                rang = int(rang[1:])
            except:
                raise Exception, 'Pas reussi a determiner le rang du coefficient {0}'.format(coeff)

            phase = rang * 2 * numpy.pi * (time + dephasage) / periode

            if liste_rangs is None or rang in liste_rangs:
                if formule is None:
                    formule = '{0}*{1}({2})'.format(coeff, op, phase)
                else:
                    formule += '+{0}*{1}({2})'.format(coeff, op, phase)

        print 'reconstruction de la serie de Fourier \n {0} ...'.format(formule[:150])

        # si les coefficients sont de types scalaires, alors c'est un scalaire que l'on calcule
        if input.GetPointData().GetArray(liste_coeff[0]).GetNumberOfComponents() == 1:
            output = CalculetteGenerique(input = input,
                formule = formule,
                nom_du_resultat = variable,
                variables_scalaires = liste_coeff
                ).get_output()
        # sinon c'est un vecteur
        elif input.GetPointData().GetArray(liste_coeff[0]).GetNumberOfComponents() == 3:
            output = CalculetteGenerique(input = input,
                formule = formule,
                nom_du_resultat = variable,
                variables_vectorielles = liste_coeff
                ).get_output()
        else:
            raise Exception, 'pas compris si la variable etait un scalaire ou un vecteur'

        if suppression_coeff_fourier is True:
            print 'suppression des arrays {0} ...'.format(liste_coeff[:5])
            for coeff_fourier in liste_coeff:
                output.GetPointData().RemoveArray(coeff_fourier)

    return output
#__________________________________________________________________________________________

#__________________________________________________________________________________________
def applatir_surface(input, axe=2, etheta_min=None, etheta_max=None):
    """Fonction qui ramene une surface dans le plan (iHat, jHat) en supprimant la composante tangentielle.
    Traite les vecteurs correctement de facon a conserver la topologie parietale.

    Indiquer l'axe de rotation:
    - 0 pour iHat
    - 1 pour jHat
    - 2 pour kHat

    si etheta_min etheta_max sont indiques (ou seulement l'un des deux) alors
    threshold sur normals.etheta pour selectionner l'une des faces de l'aubage
    """
    # cas multibloc - recursif
    if isinstance(input, vtk.vtkMultiBlockDataSet):
        output = vtk_new_instance(input)
        for numbloc in get_numeros_blocs_non_vides(input):
            output.SetBlock(numbloc, applatir_surface(input.GetBlock(numbloc), axe, etheta_min, etheta_max))

        return output

    # cas monobloc
    # calcul des vecteurs tangents a la surface qui forment un repere
    surface = convertir_en_polydata(input)
    surface = calculer_vecteur_normal(surface)
    normals = get_vtk_array_as_numpy_array(surface, 'Normals')
    surface = CalculettePyturbo(surface, ['er', 'coordr'], axe = axe).get_output()
    er = get_vtk_array_as_numpy_array(surface, 'er')
    v1 = numpy.cross(normals, er)

    v1 = v1 / numpy.sqrt(numpy.sum(v1**2, axis = 1))[:, None]
    v2 = numpy.cross(normals, v1)

    surface = ajouter_numpy_array_as_vtk_array(surface, v1, 'v1')
    surface = ajouter_numpy_array_as_vtk_array(surface, v2, 'v2')

    if etheta_min is not None or etheta_max is not None:
        surface = CalculettePyturbo(surface, 'etheta.Normals').get_output()
        surface = VTKThreshold(surface, 'etheta.Normals', etheta_min, etheta_max)
        surface = convertir_en_polydata(surface)

    # changement de repere
    coords = get_vtk_array_as_numpy_array(surface, 'coords')
    coordr = get_vtk_array_as_numpy_array(surface, 'coordr')
    coords = numpy.concatenate(
        (coords[:, axe], coordr, numpy.zeros(coordr.size)),
        axis = None).transpose()

    vtkArray = numpy_support.numpy_to_vtk(numpy.ascontiguousarray(coords), deep = 1)
    points = vtk.vtkPoints()
    points.SetData(vtkArray)

    # initialisation de output
    output = vtk_new_shallowcopy(surface)
    output.SetPoints(points)
    # on doit recalculer le vecteur normal pour le rendu graphique
    output = calculer_vecteur_normal(output)

    # on traite les vecteurs aux points
    for array in get_noms_arrays_presents(output):
        if output.GetPointData().GetArray(array).GetNumberOfComponents() == 3:
            output = CalculettePyturbo(output,
                a_calculer = "({0}.v1)*((v1.{1})*iHat+(v1.er)*jHat)+({0}.v2)*((v2.{1})*iHat+(v2.er)*jHat)"
                    .format(array, ['iHat', 'jHat', 'kHat'][axe]),
                nom_resultat = array).get_output()
    output = calculer_vecteur_normal(output)

    return output
#__________________________________________________________________________________________

max_stars_repo_name: aurmarsan/pyturbo

max_issues_repo_name: aurmarsan/pyturbo

max_forks_repo_name: aurmarsan/pyturbo


text:
"""

Usage:
    generate_pairs_predictions.py --tmf=<TestMentionsFile> --tef=<TestEmbedFile> --mf=<ModelFile> --out=<OurPredFile>
            [--cuda=<y>] [--topic=<type>] [--em=<ExtractMethod>]

Options:
    -h --help                   Show this screen.
    --cuda=<y>                  True/False - Whether to use cuda device or not [default: True]
    --topic=<type>              subtopic/topic/corpus - relevant only to ECB+, take pairs only from the same sub-topic, topic or corpus wide [default: subtopic]
    --em=<ExtractMethod>        pairwize/head_lemma/exact_string - model type to run [default: pairwize]
"""

import logging
import ntpath
import pickle
from itertools import product

import numpy as np
import random
import torch
from docopt import docopt

from src.coref_system.relation_extraction import HeadLemmaRelationExtractor, RelationTypeEnum
from src.coref_system.relation_extraction import RelationExtraction
from src.dataobjs.dataset import TopicConfig
from src.dataobjs.topics import Topics
from src.utils.embed_utils import EmbedFromFile

logger = logging.getLogger(__name__)
MAX_ALLOWED_BATCH_SIZE = 20000

def generate_pred_matrix(inference_model, topic):
    all_pairs = list(product(topic.mentions, repeat=2))
    pairs_chunks = [all_pairs]
    if len(all_pairs) > MAX_ALLOWED_BATCH_SIZE:
        pairs_chunks = [all_pairs[i:i + MAX_ALLOWED_BATCH_SIZE] for i in
                        range(0, len(all_pairs), MAX_ALLOWED_BATCH_SIZE)]
    predictions = np.empty(0)
    with torch.no_grad():
        for chunk in pairs_chunks:
            chunk_predictions, _ = inference_model.predict(chunk, bs=len(chunk))
            predictions = np.append(predictions, chunk_predictions.detach().cpu().numpy())
    predictions = 1 - predictions
    pred_matrix = predictions.reshape(len(topic.mentions), len(topic.mentions))
    return pred_matrix

def predict_and_save():
    all_predictions = list()
    for topic in _event_topics.topics_dict.values():
        logger.info("Evaluating Topic No:" + str(topic.topic_id))
        all_predictions.append((topic, generate_pred_matrix(_model, topic)))
    logger.info("Generating prediction file-" + _outfile)
    pickle.dump(all_predictions, open(_outfile, "wb"))

def get_pairwise_model():
    pairwize_model = torch.load(_model_file)
    pairwize_model.set_embed_utils(EmbedFromFile([_embed_file]))

    if _use_cuda:
        pairwize_model.cuda()

    pairwize_model.eval()
    return pairwize_model

if __name__ == '__main__':
    logging.basicConfig(level=logging.INFO)
    _arguments = docopt(__doc__, argv=None, help=True, version=None, options_first=False)
    print(_arguments)
    _mentions_file = _arguments.get("--tmf")
    _embed_file = _arguments.get("--tef")
    _model_file = _arguments.get("--mf")
    _outfile = _arguments.get("--out")
    _use_cuda = True if _arguments.get("--cuda").lower() == "true" else False
    _topic_arg = _arguments.get("--topic")
    _extract_method_str = _arguments.get("--em")

    _topic_config = Topics.get_topic_config(_topic_arg)
    _extract_method = RelationExtraction.get_extract_method(_extract_method_str)

    torch.manual_seed(1)
    random.seed(1)
    np.random.seed(1)

    logger.info("loading model from-" + ntpath.basename(_model_file))
    _event_topics = Topics()
    _event_topics.create_from_file(_mentions_file, True)

    if _topic_config == TopicConfig.Corpus and len(_event_topics.topics_dict) > 1:
        _event_topics.to_single_topic()

    _cluster_algo = None
    _model = None
    if _extract_method == RelationTypeEnum.PAIRWISE:
        _model = get_pairwise_model()
    elif _extract_method == RelationTypeEnum.SAME_HEAD_LEMMA:
        _model = HeadLemmaRelationExtractor()

    logger.info("Running agglomerative clustering with model:" + type(_model).__name__)
    predict_and_save()

max_stars_repo_name: AlonEirew/cross-doc-event-coref

max_issues_repo_name: AlonEirew/cross-doc-event-coref

max_forks_repo_name: AlonEirew/cross-doc-event-coref


text:
#test spatial neighbors
from src import spatial
from src.generate import run
from src import main
from src import data
from src.models import Hang2020
import tempfile
import glob
import geopandas as gpd
import numpy as np
import os
import pytest
import pandas as pd

ROOT = os.path.dirname(os.path.dirname(spatial.__file__))

@pytest.fixture()
def config():
    #Turn of CHM filtering for the moment
    config = data.read_config(config_path="{}/config.yml".format(ROOT))
    config["min_CHM_height"] = None
    config["iterations"] = 1
    config["rgb_sensor_pool"] = "{}/tests/data/*.tif".format(ROOT)
    config["HSI_sensor_pool"] = "{}/tests/data/*.tif".format(ROOT)
    config["min_samples"] = 1
    config["crop_dir"] = tempfile.gettempdir()
    config["bands"] = 3
    config["classes"] = 5
    config["top_k"] = 1
    config["convert_h5"] = False
    config["megaplot_dir"] = None
    config["plot_n_individuals"] = None

    return config

#Make sure deepforect crown boxes exist
@pytest.fixture()
def deepforest_boxes():
    df = gpd.read_file("{}/tests/data/sample.shp".format(ROOT))
    rgb_pool = glob.glob("{}/tests/data/*.tif".format(ROOT))
    for plot in df.plotID.unique():
        run(plot=plot, df=df, savedir=None, raw_box_savedir="{}/tests/data/interim/".format(ROOT), rgb_pool=rgb_pool)

def test_spatial_neighbors(deepforest_boxes, config):
    gdf = gpd.read_file("{}/tests/data/sample.shp".format(ROOT))
    m =  Hang2020.vanilla_CNN(bands=3, classes=5)
    #Taxon to ID dict and the reverse
    species_label_dict = {}
    for index, taxonID in enumerate(gdf.taxonID.unique()):
        species_label_dict[taxonID] = index
    m = main.TreeModel(model=m, classes=5, config=config, label_dict=species_label_dict)
    m.ROOT = "{}/tests/".format(ROOT)
    rgb_pool = glob.glob("{}/tests/data/*.tif".format(ROOT))
    neighbors = spatial.spatial_neighbors(gdf, buffer=4, HSI_pool=rgb_pool, data_dir="{}/tests/data/".format(ROOT), model=m, image_size=11)
    assert len(neighbors) == gdf.shape[0]

def test_spatial_smooth(deepforest_boxes, config):
    gdf = gpd.read_file("{}/tests/data/sample.shp".format(ROOT))
    m =  Hang2020.vanilla_CNN(bands=3, classes=5)
    #Taxon to ID dict and the reverse
    species_label_dict = {}
    for index, taxonID in enumerate(gdf.taxonID.unique()):
        species_label_dict[taxonID] = index
    m = main.TreeModel(model=m, classes=5, config=config, label_dict=species_label_dict)
    m.ROOT = "{}/tests/".format(ROOT)
    rgb_pool = glob.glob("{}/tests/data/*.tif".format(ROOT))
    neighbors = spatial.spatial_neighbors(gdf, buffer=4, HSI_pool=rgb_pool, data_dir="{}/tests/data/".format(ROOT), model=m, image_size=11)

    #Create a set of features
    features = np.zeros((gdf.shape[0],len(gdf.taxonID.unique())))
    features[0,1] = 0.95
    features[1,1] = 0.75
    labels, scores = spatial.spatial_smooth(neighbors, features)
    assert len(labels) == gdf.shape[0]
    assert len(scores) == gdf.shape[0]

max_stars_repo_name: henrykironde/DeepTreeAttention

max_issues_repo_name: henrykironde/DeepTreeAttention

max_forks_repo_name: henrykironde/DeepTreeAttention


text:
#import the necessary packages
from sklearn.neighbors import KNeighborsClassifier
from sklearn.svm import SVC
from skimage import exposure
from skimage import feature
from imutils import paths

import argparse
import imutils
import numpy as np
import cv2

# Red Color Range
HSV_RED_LOWER = np.array([0, 100, 100])
HSV_RED_UPPER = np.array([10, 255, 255])
HSV_RED_LOWER1 = np.array([160, 100, 100])
HSV_RED_UPPER1 = np.array([179, 255, 255])

# Yellow Color Range
HSV_YELLOW_LOWER = np.array([10, 80, 120])
HSV_YELLOW_UPPER = np.array([40, 255, 255])

# Blue Color Range
HSV_BLUE_LOWER = np.array([80, 160, 65])
HSV_BLUE_UPPER = np.array([140, 255, 255])

#construct the argument parse and parse command line arguments
ap = argparse.ArgumentParser()
ap.add_argument("-d", "--training", required=True, help="Path to the logos training dataset")
args = vars(ap.parse_args())

# initialize the data matrix and labels
print ("[INFO] extracting features...")

data = []
labels = []

# loop over the image paths in the training set
for imagePath in paths.list_images(args["training"]):
	# extract the make of the car
	make = imagePath.split("/")[-2]

	# load the image, convert it to grayscale, and detect edges
	# grayscale
	image = cv2.imread(imagePath)
	gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)

	# resize
	sign = cv2.resize(gray, (128, 128))

	# extract Histogram of Oriented Gradients from the sign
	H = feature.hog(sign, orientations=8, pixels_per_cell=(12, 12),
		cells_per_block=(2, 2), transform_sqrt=True, block_norm="L2")

	# update the data and labels
	data.append(H)
	labels.append(make)

# "train" the nearest neighbors classifier
print("[INFO] training classifier...")
model = KNeighborsClassifier(n_neighbors=1)

model.fit(data, labels)
print(model)
print("[INFO] evaluating...")

cap = cv2.VideoCapture(0);

while True:
	ret, img = cap.read();
	cv2.imshow("frame", img)

	#gray and hsv transform
	gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
	hsv = cv2.cvtColor(img, cv2.COLOR_BGR2HSV)

	# binary
	# 1. red
	redBinary = cv2.inRange(hsv, HSV_RED_LOWER, HSV_RED_UPPER)
	redBinary1 = cv2.inRange(hsv, HSV_RED_LOWER1, HSV_RED_UPPER1)
	redBinary = cv2.bitwise_or(redBinary, redBinary1)

	# 2. yellow
	yellowBinary = cv2.inRange(hsv, HSV_YELLOW_LOWER, HSV_YELLOW_UPPER)

	# 3. blue
	blueBinary = cv2.inRange(hsv, HSV_BLUE_LOWER, HSV_BLUE_UPPER)

	# red || blow
	binary = cv2.bitwise_or( blueBinary, (redBinary))

	# find contours
	image, contours, hierachy = cv2.findContours(binary, cv2.RETR_TREE,cv2.CHAIN_APPROX_SIMPLE)

	# 찾은 contours안을 흰색으로 채운다.(중복된 contours를 흰색으로 채워 다시 contours를 검출 한다.)
	# contours안에 contours가 잡힌 경우와 같은 영역이지만 중복해서 잡힌 경우를 예외처리한다.
	for cnt in contours:
		area = cv2.contourArea(cnt)
		binary = cv2.drawContours(binary, [cnt], -1, (255,255,255), -1)

	# 그 후 채워진 binary를 다시 contours검출
	image, goodContours, hierachy = cv2.findContours(binary, cv2.RETR_TREE,cv2.CHAIN_APPROX_SIMPLE)

	kernel = np.ones((3, 3), np.uint8)

	# 원본 grayscale와 binary를 and연산해여 mask를 씌워준다.
	gray = cv2.bitwise_and(binary, gray)

	# 추가 전처리 - binary mask를 팽창하여 mask의 크기를 넓혀준다.
	# gray = cv2.bitwise_and(cv2.dilate(binary, kernel, iterations = 1), gray)
	# cv2.imshow("hi", gray)

	for cnt in goodContours:
		area = cv2.contourArea(cnt)
		print(area)
		# 면적이 2000이상이고 가로/세로 비율이 0.6~1.4인 후보군을 bounding해준다.
		if area > 2000.0 :
			x, y, w, h = cv2.boundingRect(cnt)
			rate = w / h
			if rate > 0.8 and rate < 1.2 :
				cv2.rectangle(img, (x, y), (x+w, y+h), (200, 152, 50), 2)
				inputImage = gray[y:y+h, x:x+w]
				# kernel = np.ones((1, 1), np.uint8)
				# erosion = cv2.erode(inputImage, kernel, iterations=1)

				# 후보군을 resize해준 후 feature vector(Hog)를 추출해준다.
				sign = cv2.resize(inputImage, (128, 128))
				cv2.imshow("sign", sign)
				(H, hogImage) = feature.hog(sign, orientations=8, pixels_per_cell=(12, 12), \
					cells_per_block=(2, 2), transform_sqrt=True, block_norm="L2", visualise=True)

				cv2.imshow("hog", hogImage)

				# 학습된 모델에 넣는다.
				pred = model.predict(H.reshape(1, -1))[0]

				cv2.putText(img, pred.title(), (x, y), cv2.FONT_HERSHEY_SIMPLEX, 1, \
					(70, 255, 0), 2)

	cv2.imshow("candidates", img)
	if cv2.waitKey(1) == 27:
		break;


text:
"""
Class for converting between sphere to image coordinates using map projections.

Based on Calabretta & Greisen 2002, A&A, 357, 1077-1122
http://adsabs.harvard.edu/abs/2002A%26A...395.1077C
"""

import numpy as np

from ugali.utils.logger import logger
from ugali.utils.mlab import isstring

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

class SphericalRotator:
    """
    Base class for rotating points on a sphere.

    The input is a fiducial point (deg) which becomes (0, 0) in rotated coordinates.
    """

    def __init__(self, lon_ref, lat_ref, zenithal=False):
        self.setReference(lon_ref, lat_ref, zenithal)

    def setReference(self, lon_ref, lat_ref, zenithal=False):

        if zenithal:
            phi = (np.pi / 2.) + np.radians(lon_ref)
            theta = (np.pi / 2.) - np.radians(lat_ref)
            psi = 0.
        if not zenithal:
            phi = (-np.pi / 2.) + np.radians(lon_ref)
            theta = np.radians(lat_ref)
            psi = np.radians(90.) # psi = 90 corresponds to (0, 0) psi = -90 corresponds to (180, 0)

        cos_psi,sin_psi = np.cos(psi),np.sin(psi)
        cos_phi,sin_phi = np.cos(phi),np.sin(phi)
        cos_theta,sin_theta = np.cos(theta),np.sin(theta)

        self.rotation_matrix = np.matrix([
            [cos_psi * cos_phi - cos_theta * sin_phi * sin_psi,
             cos_psi * sin_phi + cos_theta * cos_phi * sin_psi,
             sin_psi * sin_theta],
            [-sin_psi * cos_phi - cos_theta * sin_phi * cos_psi,
             -sin_psi * sin_phi + cos_theta * cos_phi * cos_psi,
             cos_psi * sin_theta],
            [sin_theta * sin_phi,
             -sin_theta * cos_phi,
             cos_theta]
        ])

        self.inverted_rotation_matrix = np.linalg.inv(self.rotation_matrix)

    def cartesian(self,lon,lat):
        lon = np.radians(lon)
        lat = np.radians(lat)

        x = np.cos(lat) * np.cos(lon)
        y = np.cos(lat) * np.sin(lon)
        z =  np.sin(lat)
        return np.array([x,y,z])

    def rotate(self, lon, lat, invert=False):
        vec = self.cartesian(lon,lat)

        if invert:
            vec_prime = np.dot(np.array(self.inverted_rotation_matrix), vec)
        else:
            vec_prime = np.dot(np.array(self.rotation_matrix), vec)

        lon_prime = np.arctan2(vec_prime[1], vec_prime[0])
        lat_prime = np.arcsin(vec_prime[2])

        return (np.degrees(lon_prime) % 360.), np.degrees(lat_prime)

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

class Projector:
    """
    Class for performing two-dimensional map projections from the celestial sphere.
    """

    def __init__(self, lon_ref, lat_ref, proj_type = 'ait'):
        self.lon_ref = lon_ref
        self.lat_ref = lat_ref
        self.proj_type = proj_type

        if proj_type.lower() == 'ait':
            self.rotator = SphericalRotator(lon_ref, lat_ref, zenithal=False)
            self.sphere_to_image_func = aitoffSphereToImage
            self.image_to_sphere_func = aitoffImageToSphere
        elif proj_type.lower() == 'tan':
            self.rotator = SphericalRotator(lon_ref, lat_ref, zenithal=True)
            self.sphere_to_image_func = gnomonicSphereToImage
            self.image_to_sphere_func = gnomonicImageToSphere
        elif proj_type.lower() == 'car':
            def rotate(lon,lat,invert=False):
                if invert:
                    return lon + np.array([lon_ref]), lat + np.array([lat_ref])
                else:
                    return lon - np.array([lon_ref]), lat - np.array([lat_ref])
            self.rotator = SphericalRotator(lon_ref, lat_ref, zenithal=False)
            # Monkey patch the rotate function
            self.rotator.rotate = rotate
            self.sphere_to_image_func = cartesianSphereToImage
            self.image_to_sphere_func = cartesianImageToSphere
        else:
            logger.warn('%s not recognized'%(proj_type))

    def sphereToImage(self, lon, lat):
        scalar = np.isscalar(lon)

        lon, lat = np.asarray(lon), np.asarray(lat)
        lon_rotated, lat_rotated = self.rotator.rotate(lon.flat, lat.flat)
        x, y = self.sphere_to_image_func(lon_rotated, lat_rotated)

        if scalar: return np.asscalar(x), np.asscalar(y)
        else:      return x.reshape(lon.shape), y.reshape(lat.shape)

    sphere2image = sphereToImage

    def imageToSphere(self, x, y):
        scalar = np.isscalar(x)

        x, y = np.asarray(x), np.asarray(y)
        lon_rotated, lat_rotated = self.image_to_sphere_func(x.flat, y.flat)
        lon, lat = self.rotator.rotate(lon_rotated, lat_rotated, invert = True)

        if scalar: return np.asscalar(lon), np.asscalar(lat)
        else:      return lon.reshape(x.shape), lat.reshape(y.shape)

    image2sphere = imageToSphere

def sphere2image(lon_ref,lat_ref,lon,lat):
    proj = Projector(lon_ref,lat_ref)
    return proj.sphere2image(lon,lat)

def image2sphere(lon_ref,lat_ref,x,y):
    proj = Projector(lon_ref,lat_ref)
    return proj.image2sphere(x,y)

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

# ADW: Unteseted dummy projection.
def cartesianSphereToImage(lon, lat):
    lon = lon - 360.*(lon>180)
    x,y = lon,lat
    return x,y

def cartesianImageToSphere(x,y):
    x = x - 360.*(x>180)
    lon,lat = x,y
    return lon,lat

### ADW: Speed up and fixed some issues for conversion of
### lon = 180 and lon = 360 (returned by rotator)
def aitoffSphereToImage(lon, lat):
    """
    Hammer-Aitoff projection (deg).
    """
    lon = lon - 360.*(lon>180)
    lon = np.radians(lon)
    lat = np.radians(lat)

    half_lon = lon/2.
    cos_lat = np.cos(lat)

    gamma = (180. / np.pi) * np.sqrt(2. / (1. + (cos_lat * np.cos(half_lon))))
    x = 2. * gamma * cos_lat * np.sin(half_lon)
    y = gamma * np.sin(lat)
    return x, y

def aitoffImageToSphere(x, y):
    """
    Inverse Hammer-Aitoff projection (deg).
    """
    x = x - 360.*(x>180)
    x = np.asarray(np.radians(x))
    y = np.asarray(np.radians(y))
    z = np.sqrt(1. - (x / 4.)**2 - (y / 2.)**2) # rad
    lon = 2. * np.arctan2((2. * z**2) - 1, (z / 2.) * x)
    lat = np.arcsin( y * z)
    return ((180. - np.degrees(lon)) % 360.), np.degrees(lat)

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

def gnomonicSphereToImage(lon, lat):
    """
    Gnomonic projection (deg).
    """
    # Convert angle to [-180, 180] interval
    lon = lon - 360.*(lon>180)
    lon = np.radians(lon)
    lat = np.radians(lat)
    r_theta = (180. / np.pi) / np.tan(lat)
    x = r_theta * np.cos(lon)
    y = r_theta * np.sin(lon)
    return x, y

def gnomonicImageToSphere(x, y):
    """
    Inverse gnomonic projection (deg).
    """
    # Convert angle to [-180, 180] interval
    x = x - 360.*(x>180)
    x = np.asarray(x)
    y = np.asarray(y)
    lon = np.degrees(np.arctan2(y, x))
    r_theta = np.sqrt(x**2 + y**2)
    lat = np.degrees(np.arctan(180. / (np.pi * r_theta)))
    return lon, lat

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

def angsep2(lon_1, lat_1, lon_2, lat_2):
    """
    Angular separation (deg) between two sky coordinates.
    """
    import healpy

    v10, v11, v12 = healpy.ang2vec(np.radians(90. - lat_1), np.radians(lon_1)).transpose()
    v20, v21, v22 = healpy.ang2vec(np.radians(90. - lat_2), np.radians(lon_2)).transpose()
    val = (v10 * v20) + (v11 * v21) + (v12 * v22)
    val = np.clip(val, -1., 1.)
    return np.degrees(np.arccos(val))

def angsep(lon1,lat1,lon2,lat2):
    """
    Angular separation (deg) between two sky coordinates.
    Borrowed from astropy (www.astropy.org)

    Notes
    -----
    The angular separation is calculated using the Vincenty formula [1],
    which is slighly more complex and computationally expensive than
    some alternatives, but is stable at at all distances, including the
    poles and antipodes.

    [1] http://en.wikipedia.org/wiki/Great-circle_distance
    """
    lon1,lat1 = np.radians([lon1,lat1])
    lon2,lat2 = np.radians([lon2,lat2])

    sdlon = np.sin(lon2 - lon1)
    cdlon = np.cos(lon2 - lon1)
    slat1 = np.sin(lat1)
    slat2 = np.sin(lat2)
    clat1 = np.cos(lat1)
    clat2 = np.cos(lat2)

    num1 = clat2 * sdlon
    num2 = clat1 * slat2 - slat1 * clat2 * cdlon
    denominator = slat1 * slat2 + clat1 * clat2 * cdlon

    return np.degrees(np.arctan2(np.hypot(num1,num2), denominator))

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

# ADW: Reduce numpy array operations for speed
def galToCel(ll, bb):
    """
    Converts Galactic (deg) to Celestial J2000 (deg) coordinates
    """
    bb = np.radians(bb)
    sin_bb = np.sin(bb)
    cos_bb = np.cos(bb)

    ll = np.radians(ll)
    ra_gp = np.radians(192.85948)
    de_gp = np.radians(27.12825)
    lcp = np.radians(122.932)

    sin_lcp_ll = np.sin(lcp - ll)
    cos_lcp_ll = np.cos(lcp - ll)

    sin_d = (np.sin(de_gp) * sin_bb) \
            + (np.cos(de_gp) * cos_bb * cos_lcp_ll)
    ramragp = np.arctan2(cos_bb * sin_lcp_ll,
                            (np.cos(de_gp) * sin_bb) \
                            - (np.sin(de_gp) * cos_bb * cos_lcp_ll))
    dec = np.arcsin(sin_d)
    ra = (ramragp + ra_gp + (2. * np.pi)) % (2. * np.pi)
    return np.degrees(ra), np.degrees(dec)

gal2cel = galToCel

def celToGal(ra, dec):
    """
    Converts Celestial J2000 (deg) to Calactic (deg) coordinates
    """
    dec = np.radians(dec)
    sin_dec = np.sin(dec)
    cos_dec = np.cos(dec)

    ra = np.radians(ra)
    ra_gp = np.radians(192.85948)
    de_gp = np.radians(27.12825)

    sin_ra_gp = np.sin(ra - ra_gp)
    cos_ra_gp = np.cos(ra - ra_gp)

    lcp = np.radians(122.932)
    sin_b = (np.sin(de_gp) * sin_dec) \
            + (np.cos(de_gp) * cos_dec * cos_ra_gp)
    lcpml = np.arctan2(cos_dec * sin_ra_gp,
                          (np.cos(de_gp) * sin_dec) \
                          - (np.sin(de_gp) * cos_dec * cos_ra_gp))
    bb = np.arcsin(sin_b)
    ll = (lcp - lcpml + (2. * np.pi)) % (2. * np.pi)
    return np.degrees(ll), np.degrees(bb)

cel2gal = celToGal

def estimate_angle(angle, origin, new_frame, offset=1e-7):
    """
    https://github.com/astropy/astropy/issues/3093
    """
    from astropy.coordinates import SkyCoord
    import astropy.units as u
    angle_deg = angle*np.pi/180
    newlat = offset * np.cos(angle_deg) + origin.data.lat.degree
    newlon = (offset * np.sin(angle_deg) / np.cos(newlat * np.pi/180) + origin.data.lon.degree)
    sc = SkyCoord(newlon, newlat, unit='degree', frame=origin.frame.name)
    new_origin = origin.transform_to(new_frame)
    new_sc = sc.transform_to(new_frame)
    return new_origin.position_angle(new_sc).deg

def gal2cel_angle(glon,glat,angle,offset=1e-7):
    from astropy.coordinates import SkyCoord
    import astropy.units as u
    origin = SkyCoord(glon,glat,unit=u.deg,frame='galactic')
    return estimate_angle(angle,origin,'fk5',offset)

def cel2gal_angle(ra,dec,angle,offset=1e-7):
    from astropy.coordinates import SkyCoord
    import astropy.units as u
    origin = SkyCoord(ra,dec,unit=u.deg,frame='fk5')
    return estimate_angle(angle,origin,'galactic',offset)

### ADW: Probably works, remember 90-pa for kernel convention...
### def gal2cel_angle(glon,glat,angle):
###     """
###     WARNING: Not vectorized
###     """
###     gal_angle = np.radians(angle)
###     galproj = Projector(glon,glat)
###     x,y = np.sin(gal_angle),np.cos(gal_angle)
###     alon,alat = galproj.imageToSphere(0.1*x,0.1*y)
###
###     ra,dec = gal2cel(glon,glat)
###     ara,adec = gal2cel(alon,alat)
###     celproj = Projector(ra,dec)
###     cel_x,cel_y = celproj.sphereToImage(ara,adec)
###     cel_angle = np.degrees(np.arctan2(cel_x,cel_y))
###     return cel_angle + 360*(cel_angle<0)

### ADW: WARNING DOESN'T WORK YET
### def cel2gal_angle(ra,dec,angle):
###     """
###     WARNING: Not vectorized
###     """
###     cel_angle = np.radians(angle)
###     celproj = Projector(ra,dec)
###     x,y = np.sin(cel_angle),np.cos(cel_angle)
###     angle_ra,angle_dec = celproj.imageToSphere(1e-2*x,1e-2*y)
###
###     glon,glat = gal2cel(ra,dec)
###     angle_glon,angle_glat = cel2gal(angle_ra,angle_dec)
###     galproj = Projector(glon,glat)
###     gal_x,gal_y = galproj.sphereToImage(angle_glon,angle_glat)
###     gal_angle = np.degrees(np.arctan2(gal_x,gal_y))
###     return gal_angle + 360*(gal_angle<0)

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

def dec2hms(dec):
    """
    ADW: This should really be replaced by astropy
    """
    DEGREE = 360.
    HOUR = 24.
    MINUTE = 60.
    SECOND = 3600.

    dec = float(dec)
    fhour = dec*(HOUR/DEGREE)
    hour = int(fhour)

    fminute = (fhour - hour)*MINUTE
    minute = int(fminute)

    second = (fminute - minute)*MINUTE
    return (hour, minute, second)

def dec2dms(dec):
    """
    ADW: This should really be replaced by astropy
    """
    DEGREE = 360.
    HOUR = 24.
    MINUTE = 60.
    SECOND = 3600.

    dec = float(dec)
    sign = np.copysign(1.0,dec)

    fdeg = np.abs(dec)
    deg = int(fdeg)

    fminute = (fdeg - deg)*MINUTE
    minute = int(fminute)

    second = (fminute - minute)*MINUTE

    deg = int(deg * sign)
    return (deg, minute, second)

def hms2dec(hms):
    """
    Convert longitude from hours,minutes,seconds in string or 3-array
    format to decimal degrees.

    ADW: This really should be replaced by astropy
    """
    DEGREE = 360.
    HOUR = 24.
    MINUTE = 60.
    SECOND = 3600.

    if isstring(hms):
        hour,minute,second = np.array(re.split('[hms]',hms))[:3].astype(float)
    else:
        hour,minute,second = hms.T

    decimal = (hour + minute * 1./MINUTE + second * 1./SECOND)*(DEGREE/HOUR)
    return decimal

def dms2dec(dms):
    """
    Convert latitude from degrees,minutes,seconds in string or 3-array
   format to decimal degrees.
    """
    DEGREE = 360.
    HOUR = 24.
    MINUTE = 60.
    SECOND = 3600.

    # Be careful here, degree needs to be a float so that negative zero
    # can have its signbit set:
    # http://docs.scipy.org/doc/numpy-1.7.0/reference/c-api.coremath.html#NPY_NZERO

    if isstring(dms):
        degree,minute,second = np.array(re.split('[dms]',hms))[:3].astype(float)
    else:
        degree,minute,second = dms.T

    sign = np.copysign(1.0,degree)
    decimal = np.abs(degree) + minute * 1./MINUTE + second * 1./SECOND
    decimal *= sign
    return decimal

def sr2deg(solid_angle):
    return np.degrees(np.degrees(solid_angle))

def deg2sr(solid_angle):
    return np.radians(np.radians(solid_angle))

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

def distanceToDistanceModulus(distance):
    """
    Return distance modulus for a given distance (kpc).
    """
    return 5. * (np.log10(np.array(distance) * 1.e3) - 1.)

dist2mod = distanceToDistanceModulus

def distanceModulusToDistance(distance_modulus):
    """
    Return distance (kpc) for a given distance modulus.
    """
    return 10**((0.2 * np.array(distance_modulus)) - 2.)

mod2dist = distanceModulusToDistance

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

def ang2const(lon,lat,coord='gal'):
    import ephem

    scalar = np.isscalar(lon)
    lon = np.array(lon,copy=False,ndmin=1)
    lat = np.array(lat,copy=False,ndmin=1)

    if coord.lower() == 'cel':
        ra,dec = lon,lat
    elif coord.lower() == 'gal':
        ra,dec = gal2cel(lon,lat)
    else:
        msg = "Unrecognized coordinate"
        raise Exception(msg)

    x,y = np.radians([ra,dec])
    const = [ephem.constellation(coord) for coord in zip(x,y)]
    if scalar: return const[0]
    return const

def ang2iau(lon,lat,coord='gal'):
    # Default name formatting
    # http://cdsarc.u-strasbg.fr/ftp/pub/iau/
    # http://cds.u-strasbg.fr/vizier/Dic/iau-spec.htx
    fmt = "J%(hour)02i%(hmin)02i%(deg)+03i%(dmin)02i"

    scalar = np.isscalar(lon)
    lon = np.array(lon,copy=False,ndmin=1)
    lat = np.array(lat,copy=False,ndmin=1)

    if coord.lower() == 'cel':
        ra,dec = lon,lat
    elif coord.lower() == 'gal':
        ra,dec = gal2cel(lon,lat)
    else:
        msg = "Unrecognized coordinate"
        raise Exception(msg)

    x,y = np.radians([ra,dec])
    iau = []
    for _ra,_dec in zip(ra,dec):
        hms = dec2hms(_ra); dms = dec2dms(_dec)
        params = dict(hour=hms[0],hmin=hms[1],
                      deg=dms[0],dmin=dms[1])
        iau.append(fmt%params)
    if scalar: return iau[0]
    return iau

def match(lon1, lat1, lon2, lat2, tol=None, nnearest=1):
    """
    Adapted from Eric Tollerud.
    Finds matches in one catalog to another.

    Parameters
    lon1 : array-like
        Longitude of the first catalog (degrees)
    lat1 : array-like
        Latitude of the first catalog (shape of array must match `lon1`)
    lon2 : array-like
        Longitude of the second catalog
    lat2 : array-like
        Latitude of the second catalog (shape of array must match `lon2`)
    tol : float or None, optional
        Proximity (degrees) of a match to count as a match.  If None,
        all nearest neighbors for the first catalog will be returned.
    nnearest : int, optional
        The nth neighbor to find.  E.g., 1 for the nearest nearby, 2 for the
        second nearest neighbor, etc.  Particularly useful if you want to get
        the nearest *non-self* neighbor of a catalog.  To do this, use:
        ``spherematch(lon, lat, lon, lat, nnearest=2)``

    Returns
    -------
    idx1 : int array
        Indices into the first catalog of the matches. Will never be
        larger than `lon1`/`lat1`.
    idx2 : int array
        Indices into the second catalog of the matches. Will never be
        larger than `lon2`/`lat2`.
    ds : float array
        Distance (in degrees) between the matches
    """
    from scipy.spatial import cKDTree

    lon1 = np.asarray(lon1)
    lat1 = np.asarray(lat1)
    lon2 = np.asarray(lon2)
    lat2 = np.asarray(lat2)

    if lon1.shape != lat1.shape:
        raise ValueError('lon1 and lat1 do not match!')
    if lon2.shape != lat2.shape:
        raise ValueError('lon2 and lat2 do not match!')

    rotator = SphericalRotator(0,0)

    # This is equivalent, but faster than just doing np.array([x1, y1, z1]).T
    x1, y1, z1 = rotator.cartesian(lon1.ravel(),lat1.ravel())
    coords1 = np.empty((x1.size, 3))
    coords1[:, 0] = x1
    coords1[:, 1] = y1
    coords1[:, 2] = z1

    x2, y2, z2 = rotator.cartesian(lon2.ravel(),lat2.ravel())
    coords2 = np.empty((x2.size, 3))
    coords2[:, 0] = x2
    coords2[:, 1] = y2
    coords2[:, 2] = z2

    tree = cKDTree(coords2)
    if nnearest == 1:
        idxs2 = tree.query(coords1)[1]
    elif nnearest > 1:
        idxs2 = tree.query(coords1, nnearest)[1][:, -1]
    else:
        raise ValueError('invalid nnearest ' + str(nnearest))

    ds = angsep(lon1, lat1, lon2[idxs2], lat2[idxs2])

    idxs1 = np.arange(lon1.size)

    if tol is not None:
        msk = ds < tol
        idxs1 = idxs1[msk]
        idxs2 = idxs2[msk]
        ds = ds[msk]

    return idxs1, idxs2, ds

max_stars_repo_name: mcnanna/ugali

max_issues_repo_name: mcnanna/ugali

max_forks_repo_name: mcnanna/ugali

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