text:
""" Imaging is based on used of the FFT to perform Fourier transforms efficiently. Since the observed visibility data_models
do not arrive naturally on grid points, the sampled points are resampled on the FFT grid using a convolution function to
smear out the sample points. The resulting grid points are then FFT'ed. The result can be corrected for the griddata
convolution function by division in the image plane of the transform.

This approach may be extended to include image plane effect such as the w term and the antenna/station primary beam.

This module contains functions for performing the griddata process and the inverse degridding process.
"""

import logging

import numpy

log = logging.getLogger(__name__)

def coordinateBounds(npixel):
    r""" Returns lowest and highest coordinates of an image/grid given:

    1. Step size is :math:`1/npixel`:

       .. math:: \frac{high-low}{npixel-1} = \frac{1}{npixel}

    2. The coordinate :math:`\lfloor npixel/2\rfloor` falls exactly on zero:

       .. math:: low + \left\lfloor\frac{npixel}{2}\right\rfloor * (high-low) = 0

    This is the coordinate system for shifted FFTs.
    """
    if npixel % 2 == 0:
        return -0.5, 0.5 * (npixel - 2) / npixel
    else:
        return -0.5 * (npixel - 1) / npixel, 0.5 * (npixel - 1) / npixel

def coordinates(npixel: int):
    """ 1D array which spans [-.5,.5[ with 0 at position npixel/2

    """
    return (numpy.arange(npixel) - npixel // 2) / npixel

def coordinates2(npixel: int):
    """Two dimensional grids of coordinates spanning -1 to 1 in each dimension

    1. a step size of 2/npixel and
    2. (0,0) at pixel (floor(n/2),floor(n/2))
    """
    return (numpy.mgrid[0:npixel, 0:npixel] - npixel // 2) / npixel

def coordinates2Offset(npixel: int, cx: int, cy: int, quadrant=False):
    """Two dimensional grids of coordinates centred on an arbitrary point.

    This is used for A and w beams.

    1. a step size of 2/npixel and
    2. (0,0) at pixel (cx, cy,floor(n/2))
    """
    if cx is None:
        cx = npixel // 2
    if cy is None:
        cy = npixel // 2
    if quadrant == False:
        mg = numpy.mgrid[0:npixel, 0:npixel]
    else:
        # If npixel is even, we should create a grid with npixel//2+1
        mg = numpy.mgrid[0:npixel//2+1, 0:npixel//2+1]
    return (mg[0] - cy) / npixel, (mg[1] - cx) / npixel

def grdsf(nu):
    """Calculate PSWF using an old SDE routine re-written in Python

    Find Spheroidal function with M = 6, alpha = 1 using the rational
    approximations discussed by Fred Schwab in 'Indirect Imaging'.
    This routine was checked against Fred's SPHFN routine, and agreed
    to about the 7th significant digit.
    The griddata function is (1-NU**2)*GRDSF(NU) where NU is the distance
    to the edge. The grid correction function is just 1/GRDSF(NU) where NU
    is now the distance to the edge of the image.
    """
    p = numpy.array([[8.203343e-2, -3.644705e-1, 6.278660e-1, -5.335581e-1, 2.312756e-1],
                     [4.028559e-3, -3.697768e-2, 1.021332e-1, -1.201436e-1, 6.412774e-2]])
    q = numpy.array([[1.0000000e0, 8.212018e-1, 2.078043e-1],
                     [1.0000000e0, 9.599102e-1, 2.918724e-1]])

    _, np = p.shape
    _, nq = q.shape

    nu = numpy.abs(nu)

    nuend = numpy.zeros_like(nu)
    part = numpy.zeros(len(nu), dtype='int')
    part[(nu >= 0.0) & (nu < 0.75)] = 0
    part[(nu >= 0.75) & (nu <= 1.0)] = 1
    nuend[(nu >= 0.0) & (nu < 0.75)] = 0.75
    nuend[(nu >= 0.75) & (nu <= 1.0)] = 1.0

    delnusq = nu ** 2 - nuend ** 2

    top = p[part, 0]
    for k in range(1, np):
        top += p[part, k] * numpy.power(delnusq, k)

    bot = q[part, 0]
    for k in range(1, nq):
        bot += q[part, k] * numpy.power(delnusq, k)

    grdsf = numpy.zeros_like(nu)
    ok = (bot > 0.0)
    grdsf[ok] = top[ok] / bot[ok]
    ok = numpy.abs(nu > 1.0)
    grdsf[ok] = 0.0

    # Return the griddata function and the grid correction function
    return grdsf, (1 - nu ** 2) * grdsf

def w_beam(npixel, field_of_view, w, cx=None, cy=None, remove_shift=False):
    """ W beam, the fresnel diffraction pattern arising from non-coplanar baselines

    :param npixel: Size of the grid in pixels
    :param field_w_beamof_view: Field of view
    :param w: Baseline distance to the projection plane
    :param cx: location of delay centre def :npixel//2
    :param cy: location of delay centre def :npixel//2
    :param remove_shift: Remove overall phase shift at the centre of the image
    :return: npixel x npixel array with the far field
    """
    if cx is None:
        cx = npixel // 2
    if cy is None:
        cy = npixel // 2

    # Original codes
    # ly, mx = coordinates2Offset(npixel, cx, cy)
    # r2 = field_of_view**2*(ly ** 2 + mx ** 2)
    # ph = numpy.zeros_like(r2)
    # ph[r2 < 1.0] = -2 * numpy.pi * w * (1 - numpy.sqrt(1.0 - r2[r2 < 1.0]))
    # cp = numpy.zeros_like(r2, dtype='complex')
    # cp[r2 < 1.0] = numpy.exp(1j * ph[r2 < 1.0])
    # cp[r2 == 0] = 1.0 + 0j
    # if remove_shift:
    #     cp /= cp[npixel // 2, npixel // 2]

    # numpy.putmask
    # ly, mx = coordinates2Offset(npixel, cx, cy)
    # r2 = field_of_view**2*(ly ** 2 + mx ** 2)
    # ph = numpy.zeros_like(r2)
    # m = r2 < 1.0
    # cp = numpy.zeros_like(r2, dtype='complex')
    # numpy.putmask(ph, m, -2 * numpy.pi * w * (1 - numpy.sqrt(1.0 - r2)))
    # numpy.putmask(cp, m, numpy.exp(1j * ph))
    # numpy.putmask(cp, r2 == 0, 1.0 + 0j)
    # if remove_shift:
    #     cp /= cp[npixel // 2, npixel // 2]

    # numpy.putmask - 2
    # ly, mx = coordinates2Offset(npixel, cx, cy)
    # r2 = field_of_view ** 2 * (ly ** 2 + mx ** 2)
    # ph = -2 * numpy.pi * w * (1 - numpy.sqrt(1.0 - r2))
    # numpy.putmask(ph, r2 >= 1.0, 0)
    # cp = numpy.zeros_like(r2, dtype='complex')
    # cp = numpy.exp(1j * ph)
    # numpy.putmask(cp, r2 >= 1.0, 0 + 0j)
    # numpy.putmask(cp, r2 == 0, 1.0 + 0j)
    # if remove_shift:
    #     cp /= cp[npixel // 2, npixel // 2]

    # SubArray Copy Symmetrically
    ly, mx = coordinates2Offset(npixel, cx, cy, quadrant=True)
    r2 = field_of_view ** 2 * (ly ** 2 + mx ** 2)
    ph = -2 * numpy.pi * w * (1 - numpy.sqrt(1.0 - r2))
    numpy.putmask(ph, r2 >= 1.0, 0)
    cp = numpy.zeros_like(r2, dtype='complex')
    cp = numpy.exp(1j * ph)
    numpy.putmask(cp, r2 >= 1.0, 0 + 0j)
    numpy.putmask(cp, r2 == 0, 1.0 + 0j)
    # Correct for linear phase shift in faceting
    if remove_shift:
        cp /= cp[-1, -1]

    cp = numpy.pad(cp, ((0, int(cx) + npixel % 2 - 1), (0, int(cy) + npixel % 2 - 1)), 'reflect')

    # assert((cp==cp1).all())

    return cp

def frac_coord(npixel, kernel_oversampling, p):
    """ Compute whole and fractional parts of coordinates, rounded to
    kernel_oversampling-th fraction of pixel size

    The fractional values are rounded to nearest 1/kernel_oversampling pixel value. At
    fractional values greater than (kernel_oversampling-0.5)/kernel_oversampling coordinates are
    rounded to next integer index.

    :param npixel: Number of pixels in total
    :param kernel_oversampling: Fractional values to round to
    :param p: Coordinate in range [-.5,.5[
    """
    assert numpy.array(p >= -0.5).all() and numpy.array(
        p < 0.5).all(), "Cellsize is too large: uv overflows grid uv= %s" % str(p)
    x = npixel // 2 + p * npixel
    flx = numpy.floor(x + 0.5 / kernel_oversampling)
    fracx = numpy.around((x - flx) * kernel_oversampling)
    return flx.astype(int), fracx.astype(int)

def weight_gridding(shape, visweights, vuvwmap, vfrequencymap, vpolarisationmap=None, weighting='uniform'):
    """Reweight data using one of a number of algorithms

    :param shape:
    :param visweights: Visibility weights
    :param vuvwmap: map uvw to grid fractions
    :param vfrequencymap: map frequency to image channels
    :param vpolarisationmap: map polarisation to image polarisation
    :param weighting: '' | 'uniform'
    :return: visweights, density, densitygrid
    """
    densitygrid = numpy.zeros(shape, dtype='float')
    if weighting == 'uniform':
        log.debug("weight_gridding: Performing uniform weighting")
        inchan, inpol, ny, nx = shape

        wts = visweights[...]
        # uvw -> fraction of grid mapping
        for flip in [-1.0, 1.0]:
            y, yf = frac_coord(ny, 1.0, flip * vuvwmap[:, 1])
            x, xf = frac_coord(nx, 1.0, flip * vuvwmap[:, 0])
            coords = list(vfrequencymap), x, y
            for pol in range(inpol):
                for vwt, chan, x, y in zip(wts, *coords):
                    densitygrid[chan, pol, y, x] += vwt[..., pol]

        # Find the total weight per sample counting redundancies with other samples
        newvisweights = numpy.zeros_like(visweights)
        density = numpy.zeros_like(visweights)
        y, _ = frac_coord(ny, 1.0, vuvwmap[:, 1])
        x, _ = frac_coord(nx, 1.0, vuvwmap[:, 0])
        coords = list(vfrequencymap), x, y
        for pol in range(inpol):
            density[..., pol] += [densitygrid[chan, pol, y, x] for chan, x, y in zip(*coords)]

        # Normalise each visibility weight to sum to one in a grid cell
        if numpy.sum(density[:, 0] > 0.0) < visweights.shape[0]:
            log.warning("weight_gridding: Losing samples in weighting")

        newvisweights[density > 0.0] = visweights[density > 0.0] / density[density > 0.0]
        return newvisweights, density, densitygrid
    else:
        return visweights, None, None

def visibility_recentre(uvw, dl, dm):
    """ Compensate for kernel re-centering - see `w_kernel_function`.

    :param uvw: Visibility coordinates
    :param dl: Horizontal shift to compensate for
    :param dm: Vertical shift to compensate for
    :returns: Visibility coordinates re-centrered on the peak of their w-kernel
    """

    u, v, w = numpy.hsplit(uvw, 3)  # pylint: disable=unbalanced-tuple-unpacking
    return numpy.hstack([u - w * dl, v - w * dm, w])


text:
# Custom transform class for pytorch dataloader

import cv2
import torch
import time
import numpy as np
from PIL import Image
from functools import partial
from scipy.ndimage import map_coordinates

np_pad_name = {'zero':'constant',
               'circular': 'wrap',
               'edge': 'edge'}

class CustomPad(object):
    """
    Pad image with given mode
    input:
        pad: int or 4-elements tuple
        img: PIL image or ndarray
    """

    def __init__(self, pad, mode, **kwargs):
        assert mode in np_pad_name.keys(), \
        'Padding mode zero and circular are implemented for now'
        self.mode = np_pad_name[mode]

        assert isinstance(pad, (int, tuple))
        if isinstance(pad, int):
            self.pad = (pad, pad, pad, pad)
        else:
            assert len(pad) == 4, \
            'Padding length must be an integer or a 4-element tuple'
            self.pad = pad

        self.kwargs = kwargs
        self.numpy_mode = False

    def __call__(self, img):
        '''
        img: ndarray or PIL image
        '''
        if not isinstance(img, (Image.Image, np.ndarray)):
            raise TypeError('img should be either PIL Image or ndarray. Got {}'.format(type(img)))
        if isinstance(img, np.ndarray):
            self.numpy_mode = True

        pad_left = self.pad[0]
        pad_top = self.pad[1]
        pad_right = self.pad[2]
        pad_bottom = self.pad[3]

        if not self.numpy_mode and img.mode == 'P':
            palette = img.getpalette()
            img = np.asarray(img)
            img = np.pad(img, ((pad_top, pad_bottom), (pad_left, pad_right)), self.mode, **self.kwargs)
            img = Image.fromarray(img)
            img.putpalette(palette)
            return img

        # Turn to RGB if is PIL image
        if not self.numpy_mode:
            img = np.asarray(img)

        if len(img.shape) == 3:
            img = np.pad(img, ((pad_top, pad_bottom), (pad_left, pad_right), (0, 0)), self.mode, **self.kwargs)

        # Grayscale image
        if len(img.shape) == 2:
            img = np.pad(img, ((pad_top, pad_bottom), (pad_left, pad_right)), self.mode, **self.kwargs)

        if not self.numpy_mode:
            return Image.fromarray(img)

        return img

class Resize(object):
    '''
    Resize ndarray to target size
    Take ndarray as input
    '''
    def __init__(self, size, interpolation=cv2.INTER_LINEAR):
        '''
        size: HxW
        '''
        if isinstance(size, (tuple, list)):
            self.size = (size[1], size[0]) # Cv2 uses WxH
        if isinstance(size, int):
            self.size = (size, size)
        self.interpolation = interpolation

    def __call__(self, img):
        if not isinstance(img, np.ndarray):
            raise TypeError("Custom ToTensor takes numpy array as input")

        return cv2.resize(img, dsize=self.size, interpolation=self.interpolation)

class LogPolartoPolar(object):
    '''
    Transfer x from log to linear space
    Take ndarray as input
    '''
    require_grid = True

    def __call__(self, log_grid, img):
        '''
        Arround .7 sec per call on i7-8750h
        '''
        nx, ny, nchan = img.shape

        # Scale log_grid numbers to actual image pixel index
        log_grid_expanded = (log_grid - log_grid[0]) / (log_grid[-1] - log_grid[0]) * ny
        digitized_x = np.digitize(np.linspace(0, ny, ny), log_grid_expanded)-1
        rp1, rp2 = np.meshgrid(np.linspace(0, nx, nx), digitized_x)
        one_chan_imgs = list(map(lambda chan: np.expand_dims(map_coordinates(img[:,:,chan], (rp2, rp1)), 2), range(nchan)))
        polar_data = np.concatenate(one_chan_imgs, axis=2)

        return polar_data

class Normalize(object):
    '''
    Normalize img values channel-wise
    '''
    def __init__(self, mean, std):
        self.mean = mean
        self.std = std

    def __call__(self, img):
        '''
        img: ndarray in HWC, or tensor in CHW
        '''
        ret = (img - self.mean) / self.std

        return ret

class ToTensor(object):
    '''
    Convet numpy array of HWC to tensor CHW
    Take ndarray as input
    '''
    def __call__(self, img):
        if not isinstance(img, np.ndarray):
            raise TypeError("Custom ToTensor takes numpy array as input")

        return torch.from_numpy(np.transpose(img, (2, 0, 1)))

class Crop(object):
    '''
    Crop image
    '''
    def __init__(self, top_left, size):
        '''
        tl: top left conner
        size: output img size of HxW
        '''
        self.tl = top_left
        self.br = (top_left[0] + size[0], top_left[1] + size[1])
        self.size = size

    def __call__(self, img):
        '''
        img: ndarray
        '''
        h, w = img.shape[:-1]
        if self.br[0] > h or self.br[1] > w:
            raise ValueError("Impossible to crop form input img of size {h}x{w}".format())

        return img[self.tl[0]:self.br[0], self.tl[1]:self.br[1]]

class chan():
    '''
    Return a channel as output
    '''
    def __init__(self, chan):
        self.chan = chan

    def __call__(self, img):
        if isinstance(img, np.ndarray):
            return img[:,:,self.chan]
        if isinstance(img, torch.Tensor):
            return img[self.chan]

        raise ValueError(
            "To pick a channel, input img must be either an ndarray or torch tensor")

class GroupRandomCrop(object):
    '''
    Return a random crop at the same position of three input imgs
    input: 2 input images and 1 label image
    '''
    group_tran = True

    def __init__(self, size, label_size = None):
        '''
        size: height x width
        label_size: label image size
        '''
        if isinstance(size, int):
            self.size = (size, size)
        else:
            self.size = size
        if isinstance(label_size, int):
            self.label_size = (label_size, label_size)
        else:
            self.label_size = label_size

        self.inner_pad_size = (int(0.5 * (self.size[0] - self.label_size[0])), int(0.5 * (self.size[1] - self.label_size[1])))
        if self.inner_pad_size[0] * self.inner_pad_size[1] < 0:
            raise ValueError(
                "Label size {} must be smaller then img size {}".format(self.label_size, self.size))

        self.numpy_mode = False

    def get_params(self, h, w):
        '''
        Return i, j
        i, j: coordinates of upper left conner
        '''
        if h == self.size[0]:
            i = 0
        else:
            i = np.random.randint(0, h - self.size[0])
        if w == self.size[1]:
            j = 0
        else:
            j = np.random.randint(0, w - self.size[1])

        return i, j

    def __call__(self, low, high, mid):
        '''
        l, h, m: img_t0, img_t1, img_t0.5
        '''

        if not isinstance(low, (Image.Image, np.ndarray)):
            raise TypeError(
                'img should be either PIL Image or ndarray. Got {}'.format(type(low)))
        if isinstance(low, np.ndarray):
            self.numpy_mode = True

        if self.numpy_mode:
            h, w = low.shape[:2]
        else:
            w, h = low.size

        if self.size[0] > h or self.size[1] > w:
            raise ValueError(
                "Crop output size {}x{} must be <= {}x{}, or pad more first".format(self.size[0], self.size[1], h, w))

        i, j = self.get_params(h, w)
        br = (i + self.size[0], j + self.size[1])

        if not self.numpy_mode:
            return low.crop((j, i, br[1], br[0])), \
                   high.crop((j, i, br[1], br[0])), \
                   mid.crop((j+self.inner_pad_size[1], i+self.inner_pad_size[0], j+self.label_size[1], i+self.label_size[0]))

        return low[i:br[0], j:br[1]], \
               high[i:br[0], j:br[1]], \
               mid[i+self.inner_pad_size[0]:i+self.inner_pad_size[0]+self.label_size[0], j+self.inner_pad_size[1]:j+self.inner_pad_size[1]+self.label_size[1]]

    def __repr__(self):
        return self.__class__.__name__ + '(size={0}, label_size={})'.format(self.size, self.label_size)

# class GroupRandomCrop2(object):
#     '''
#     Return a random crop at the same position of three input imgs
#     input: 2 input images and 1 label image
#     '''
#     def __init__(self, sizes):
#         '''
#         sizes: ((height, width), (height, width), ...) or (height, width)
#         '''
#         assert isinstance(sizes, (tuple, list)), \
#         "Sizes must be in format" + \
#         "((height, width), (height, width), ...) or (height, width)"

#         self.sizes = []
#         if len(sizes) == 2 and isinstance(sizes[0], int) and isinstance(sizes[1], int):
#             self.sizes.append(list(sizes))
#         else:
#             for pair in sizes:
#                 if not isinstance(pair, (tuple, list)) or \
#                 len(pair) != 2 or \
#                 not isinstance(pair[0], int) or \
#                 not isinstance(pair[1], int):
#                     raise ValueError(
#                         "Sizes must be in format" +
#                         "((height, width), (height, width), ...) or (height, width)")
#                 self.sizes.append(list(pair))

#         self.numpy_mode = False

    # def get_params(self, h, w):
    #     '''
    #     Return i, j
    #     i, j: coordinates of upper left conner
    #     '''
    #     out_h_max = np.array(self.sizes)
    #     i = np.random.randint(0, h - self.size[0])
    #     j = np.random.randint(0, h - self.size[1])

    #     return i, j

    # def __call__(self, imgs):
    #     '''
    #     l, h, m: img_t0, img_t1, img_t0.5
    #     '''
    #     if not isinstance(low, (Image.Image, np.ndarray)):
    #         raise TypeError(
    #             'img should be either PIL Image or ndarray. Got {}'.format(type(l)))
    #     if isinstance(low, np.ndarray):
    #         self.numpy_mode = True

    #     if self.numpy_mode:
    #         h, w = imgs[0].shape[:2]
    #     else:
    #         w, h = imgs[0].size

    #     for pair in self.sizes:
    #         if self.pair[0] > h or self.pair[1] > w:
    #             raise ValueError(
    #                 "Crop output sizes {} must all <= {}, or pad more first".format(self.sizes, (h, w)))

    #     i, j = self.get_params(h, w)
    #     if not self.numpy_mode:
    #         return low.crop((j, i, j+self.size[1], i+self.size[0])), high.crop((j, i, j+self.size[1], i+self.size[0])), mid.crop((j, i, j+self.size[1], i+self.size[0]))

    #     return low[i:i+self.size[0], j:j+self.size[1]], high[i:i+self.size[0], j:j+self.size[1]], mid[i:i+self.size[0], j:j+self.size[1]]

    # def __repr__(self):
    #     return self.__class__.__name__ + '(size={0})'.format(self.size)

ext: py

max_stars_repo_name: superpenshine/skylark

max_issues_repo_name: superpenshine/skylark

max_forks_repo_name: superpenshine/skylark

reason: import numpy,from scipy


text:
import numpy as np
import torch
import logging
import losses
import json
from tqdm import tqdm
import torch.nn.functional as F
import math
import time
from sklearn.metrics.pairwise import pairwise_distances

def re_rank(q_g_dist, q_q_dist, g_g_dist, k1=50, k2=20, lambda_value=0.94):

    original_dist = np.concatenate(
      [np.concatenate([q_q_dist, q_g_dist], axis=1),
       np.concatenate([q_g_dist.T, g_g_dist], axis=1)],
      axis=0)
    original_dist = np.power(original_dist, 2).astype(np.float32)
    original_dist = np.transpose(1. * original_dist/np.max(original_dist,axis = 0))
    V = np.zeros_like(original_dist).astype(np.float32)
    initial_rank = np.argsort(original_dist).astype(np.int32)

    query_num = q_g_dist.shape[0]
    gallery_num = q_g_dist.shape[0] + q_g_dist.shape[1]
    all_num = gallery_num

    for i in range(all_num):
        # k-reciprocal neighbors
        forward_k_neigh_index = initial_rank[i,:k1+1]
        backward_k_neigh_index = initial_rank[forward_k_neigh_index,:k1+1]
        fi = np.where(backward_k_neigh_index==i)[0]
        k_reciprocal_index = forward_k_neigh_index[fi]
        k_reciprocal_expansion_index = k_reciprocal_index
        for j in range(len(k_reciprocal_index)):
            candidate = k_reciprocal_index[j]
            candidate_forward_k_neigh_index = initial_rank[candidate,:int(np.around(k1/2.))+1]
            candidate_backward_k_neigh_index = initial_rank[candidate_forward_k_neigh_index,:int(np.around(k1/2.))+1]
            fi_candidate = np.where(candidate_backward_k_neigh_index == candidate)[0]
            candidate_k_reciprocal_index = candidate_forward_k_neigh_index[fi_candidate]
            if len(np.intersect1d(candidate_k_reciprocal_index,k_reciprocal_index))> 2./3*len(candidate_k_reciprocal_index):
                k_reciprocal_expansion_index = np.append(k_reciprocal_expansion_index,candidate_k_reciprocal_index)

        k_reciprocal_expansion_index = np.unique(k_reciprocal_expansion_index)
        weight = np.exp(-original_dist[i,k_reciprocal_expansion_index])
        V[i,k_reciprocal_expansion_index] = 1.*weight/np.sum(weight)
    original_dist = original_dist[:query_num,]
    if k2 != 1:
        V_qe = np.zeros_like(V,dtype=np.float32)
        for i in range(all_num):
            V_qe[i,:] = np.mean(V[initial_rank[i,:k2],:],axis=0)
        V = V_qe
        del V_qe
    del initial_rank
    invIndex = []
    for i in range(gallery_num):
        invIndex.append(np.where(V[:,i] != 0)[0])

    jaccard_dist = np.zeros_like(original_dist,dtype = np.float32)

    for i in range(query_num):
        temp_min = np.zeros(shape=[1,gallery_num],dtype=np.float32)
        indNonZero = np.where(V[i,:] != 0)[0]
        indImages = []
        indImages = [invIndex[ind] for ind in indNonZero]
        for j in range(len(indNonZero)):
            temp_min[0,indImages[j]] = temp_min[0,indImages[j]]+ np.minimum(V[i,indNonZero[j]],V[indImages[j],indNonZero[j]])
        jaccard_dist[i] = 1-temp_min/(2.-temp_min)

    final_dist = jaccard_dist*(1-lambda_value) + original_dist*lambda_value
    del original_dist
    del V
    del jaccard_dist
    final_dist = final_dist[:query_num,query_num:]
    return final_dist

def l2_norm(input):
    input_size = input.size()
    buffer = torch.pow(input, 2)
    normp = torch.sum(buffer, 1).add_(1e-12)
    norm = torch.sqrt(normp)
    _output = torch.div(input, norm.view(-1, 1).expand_as(input))
    output = _output.view(input_size)

    return output

def calc_recall_at_k(T, Y, k):
    """
    T : [nb_samples] (target labels)
    Y : [nb_samples x k] (k predicted labels/neighbours)
    """

    s = 0
    for t,y in zip(T,Y):
        if t in torch.Tensor(y).long()[:k]:
            s += 1
    return s / (1. * len(T))

def predict_batchwise(model, dataloader):
    device = "cuda"
    model_is_training = model.training
    model.eval()

    ds = dataloader.dataset
    A = [[] for i in range(len(ds[0]))]
    with torch.no_grad():
        # extract batches (A becomes list of samples)
        for batch in tqdm(dataloader):
            for i, J in enumerate(batch):
                # i = 0: sz_batch * images
                # i = 1: sz_batch * labels
                # i = 2: sz_batch * indices
                if i == 0:
                    # move images to device of model (approximate device)
                    J = model(J.cuda())

                for j in J:
                    A[i].append(j)
    model.train()
    model.train(model_is_training) # revert to previous training state

    return [torch.stack(A[i]) for i in range(len(A))]

def proxy_init_calc(model, dataloader):
    nb_classes = dataloader.dataset.nb_classes()
    X, T, *_ = predict_batchwise(model, dataloader)

    proxy_mean = torch.stack([X[T==class_idx].mean(0) for class_idx in range(nb_classes)])

    return proxy_mean

def evaluate_cos(model, dataloader,rerank=True):
    nb_classes = dataloader.dataset.nb_classes()

    # calculate embeddings with model and get targets
    X, T = predict_batchwise(model, dataloader)
    X = l2_norm(X)
    X=X.cpu()
    if(rerank):
         distances_q_g=pairwise_distances(X,X)
         distances_g_g= distances_q_g
         distances_q_q=distances_q_g
         distances=re_rank(distances_q_g, distances_q_q, distances_g_g, k1=50, k2=20, lambda_value=.94)
    else:
        distances=pairwise_distances(X,X)
    # get predictions by assigning nearest 8 neighbors with cosine
    K = 32
    Y = []
    xs = []
    similarities= torch.from_numpy(-distances)
    Y = T[similarities.topk(1 + K)[1][:,1:]]
    Y = Y.float().cpu()

    recall = []
    for k in [1, 2, 4, 8, 16, 32]:
        r_at_k = calc_recall_at_k(T, Y, k)
        recall.append(r_at_k)
        print("R@{} : {:.3f}".format(k, 100 * r_at_k))

    return recall

def evaluate_Inshop(model, query_dataloader, gallery_dataloader,rerank=False):
    # nb_classes = query_dataloader.dataset.nb_classes()

    # # calculate embeddings with model and get targets
    query_X, query_T = predict_batchwise_inshop(model, query_dataloader)
    gallery_X, gallery_T = predict_batchwise_inshop(model, gallery_dataloader)

    query_X=query_X.cpu()
    gallery_X=gallery_X.cpu()

    if(rerank):
         distances_q_g = pairwise_distances(query_X, gallery_X)
         distances_g_g = pairwise_distances(gallery_X, gallery_X)
         distances_q_q = pairwise_distances(query_X, query_X)
         distances=re_rank(distances_q_g, distances_q_q, distances_g_g, k1=6, k2=2, lambda_value=0)
    else:
         distances=pairwise_distances(query_X,gallery_X)

    distances=torch.from_numpy(distances)
    # calculate recall @ 1, 2, 4, 8
    recall = []

    def calc_recall_at_k_inshop(distances, query_T, gallery_T, k):
        m = len(distances)
        match_counter = 0

        for i in range(m):
            pos_dist = distances[i][gallery_T == query_T[i]]
            neg_dist = distances[i][gallery_T != query_T[i]]
            thresh = torch.min(pos_dist).item()

            if torch.sum(neg_dist < thresh) < k:
                match_counter += 1

        return match_counter / m

    for k in [1, 10, 20, 30, 40, 50]:
        r_at_k = calc_recall_at_k_inshop(distances, query_T, gallery_T, k)
        recall.append(r_at_k)
        print("R@{} : {:.3f}".format(k, 100 * r_at_k))

    return recall

def evaluate_cos_SOP(model, dataloader):
    nb_classes = dataloader.dataset.nb_classes()

    # calculate embeddings with model and get targets
    X, T = predict_batchwise(model, dataloader)
    X = l2_norm(X)

    # get predictions by assigning nearest 8 neighbors with cosine
    K = 1000
    Y = []
    xs = []
    for x in X:
        if len(xs)<10000:
            xs.append(x)
        else:
            xs.append(x)
            xs = torch.stack(xs,dim=0)
            cos_sim = F.linear(xs,X)
            y = T[cos_sim.topk(1 + K)[1][:,1:]]
            Y.append(y.float().cpu())
            xs = []

    # Last Loop
    xs = torch.stack(xs,dim=0)
    cos_sim = F.linear(xs,X)
    y = T[cos_sim.topk(1 + K)[1][:,1:]]
    Y.append(y.float().cpu())
    Y = torch.cat(Y, dim=0)

    # calculate recall @ 1, 2, 4, 8
    recall = []
    for k in [1, 10, 100, 1000]:
        r_at_k = calc_recall_at_k(T, Y, k)
        recall.append(r_at_k)
        print("R@{} : {:.3f}".format(k, 100 * r_at_k))
    return recall

reason: import numpy


text:
#! /usr/bin/env python

##############################################################################
import numpy as np
import emcee
import os, sys
import scipy.interpolate as si

import lmfit as LM

import SpectralFitting
import SpectralFitting_functs as SF

from emcee.utils import MPIPool

#Likelihood function here: if we put it in the Spectral Fitting class
#we'd need to pickle the fit_settings dictionary, which saves pickling the parameters dictionary
#And speeds things up by a factor of ~30
def lnprob(T, theta, var_names, bounds, ret_specs=False):

    #Log prob function. T is an array of values

    assert len(T)==len(var_names), 'Error! The number of variables and walker position shapes are different'

    #Prior information comes from the parameter bounds now
    if np.any(T > bounds[:, 1]) or np.any(T < bounds[:, 0]):
        return -np.inf

    #make theta from the emcee walker positions
    for name, val in zip(var_names, T):
        theta[name].value = val

    if ret_specs==False:
        ll=SF.lnlike(theta, fit.fit_settings)
        return ll
    else:
        return SF.lnlike(theta, fit.fit_settings, ret_specs=True)

#Can select either Kroupa or Salpeter to use with the SSP models
element_imf='kroupa'

####################################################
#Read in the data
datafile='data/example_spectrum.txt'

lamdas, flux, errors, instrumental_resolution=np.genfromtxt(datafile, unpack=True)

# The instrumental resolution can be included if it's known. We need a value of sigma_inst in km/s for every pixel
# Otherwise leave it as None
#instrumental_resolution=None

# Sky Spectra
# Give a list of 1D sky spectra to be scaled and subtracted during the fit
# Otherwise leave sky as None
skyspecs=None
# ######################################################

# Mask out regions that we don't want to fit, e.g. around telluric residuals, particularly nasty skylines, etc
# THESE SHOULD BE OBSERVED WAVELENGTHS
telluric_lam_1=np.array([[6862, 6952]])
telluric_lam_2=np.array([[7586, 7694]])
skylines=np.array([[8819, 8834], [8878.0, 8893], [8911, 8925], [8948, 8961]])

masked_wavelengths=np.vstack([telluric_lam_1, telluric_lam_2, skylines]).reshape(-1, 1, 2)
string_masked_wavelengths=["{} to {}".format(pair[0][0], pair[0][1]) for pair in masked_wavelengths]

#Mask pixels we don't want
pixel_mask=np.ones_like(flux, dtype=bool)
for array in masked_wavelengths:
    m=SF.make_mask(lamdas, array)
    pixel_mask= m & pixel_mask

#Now switch the weights of these pixels to 0
pixel_weights=np.ones_like(flux)
pixel_weights[~pixel_mask]=0.0

#Wavelengths we'll fit between.
#Split into 4 to make the multiplicative polynomials faster
fit_wavelengths=np.array([[4600, 5600], [5600, 6800], [6800, 8000], [8000,  9000], [9700, 10500]])
string_fit_wavelengths=["{} to {}".format(pair[0], pair[1]) for pair in fit_wavelengths]

#FWHM. Should make a way to measure this!
FWHM_gal=2.5

#Now set up the spectral fitting class
print 'Setting up the fit'
fit=SpectralFitting.SpectralFit(lamdas, flux, errors, pixel_weights, fit_wavelengths, FWHM_gal, instrumental_resolution=instrumental_resolution, skyspecs=skyspecs, element_imf=element_imf)
fit.set_up_fit()

# MPI usage taken from the emcee docs
with MPIPool() as pool:
    if not pool.is_master():
        pool.wait()
        sys.exit(0)

    #Here are the available fit parameters
    #They can easily be switched off by changing vary to False
    #The min and max values act as flat priors
    theta=LM.Parameters()
    #LOSVD parameters
    theta.add('Vel', value=1800, min=-1000.0, max=10000.0)
    theta.add('sigma', value=330.0, min=10.0, max=500.0)

    #Abundance of Na. Treat this separately, since it can vary up to +1.0 dex
    theta.add('Na', value=0.5, min=-0.45, max=1.0, vary=True)

    #Abundance of elements which can vary positively and negatively
    theta.add('Ca', value=0.0,  min=-0.45, max=0.45, vary=True)
    theta.add('Fe', value=0.0, min=-0.45, max=0.45, vary=True)
    theta.add('C', value=0.0, min=-0.45, max=0.45, vary=True)
    theta.add('N', value=0.0, min=-0.45, max=0.45, vary=True)
    theta.add('Ti', value=0.0, min=-0.45, max=0.45, vary=True)
    theta.add('Mg', value=0.0, min=-0.45, max=0.45, vary=True)
    theta.add('Si', value=0.0, min=-0.45, max=0.45, vary=True)
    theta.add('Ba', value=0.0, min=-0.45, max=0.45, vary=True)

    #Abundance of elements which can only vary above 0.0
    theta.add('as_Fe', value=0.0, min=0.0, max=0.45, vary=True)
    theta.add('Cr', value=0.0, min=0.0, max=0.45, vary=False)
    theta.add('Mn', value=0.0, min=0.0, max=0.45, vary=False)
    theta.add('Ni', value=0.0, min=0.0, max=0.45, vary=False)
    theta.add('Co', value=0.0, min=0.0, max=0.45, vary=False)
    theta.add('Eu', value=0.0, min=0.0, max=0.45, vary=False)
    theta.add('Sr', value=0.0, min=0.0, max=0.45, vary=False)
    theta.add('K', value=0.0, min=0.0, max=0.45, vary=False)
    theta.add('V', value=0.0, min=0.0, max=0.45, vary=False)
    theta.add('Cu', value=0.0, min=0.0, max=0.45, vary=False)

    #Emission line kinematics
    #Each line is fixed to the same velocity and sigma
    theta.add('Vel_em', value=1800, min=0.0, max=10000)
    theta.add('sig_em', value=200.0, min=10.0, max=500.0)

    #Emission line strengths
    #These are log flux- they get exponentiated in the likelihood function
    theta.add('Ha', value=-1.5, min=-10.0, max=10.0)
    theta.add('Hb', value=-2.0, min=-10.0, max=10.0)
    theta.add('NII', value=-1.0, min=-10.0, max=10.0)
    theta.add('SII_6716', value=-2.0, min=-10.0, max=10.0)
    theta.add('SII_6731', value=-2.0, min=-10.0, max=10.0)
    theta.add('OIII', value=-2.0, min=-10.0, max=10.0)
    theta.add('OI', value=-2.0, min=-10.0, max=10.0)

    #Base population parameters
    #Age, Metallicity, and the two IMF slopes
    theta.add('age', value=13.0, min=1.0, max=14.0)
    theta.add('Z', value=0.0, min=-1.0, max=0.2)
    theta.add('imf_x1', value=2.35, min=0.5, max=3.5)
    theta.add('imf_x2', value=2.35, min=0.5, max=3.5)

    #Strengths of skylines
    theta.add('O2_Scale', value=0.0, min=-100000000, max=100000000, vary=False)
    theta.add('sky_Scale', value=0.0, min=-100000000, max=100000000, vary=False)
    theta.add('OH_Scale', value=0.0, min=-100000000, max=100000000, vary=False)
    theta.add('NaD_sky_scale', value=0.0, min=-100000000, max=100000000, vary=False)

    #Option to rescale the error bars up or down
    theta.add('ln_f', value=0.0, min=-5.0, max=5.0, vary=True)

    #Select the parameters we're varying, ignore the fixed ones
    variables=[thing for thing in theta if theta[thing].vary]
    ndim=len(variables)
    #Vice versa, plus add in the fixed value
    fixed=[ "{}={},".format(thing, theta[thing].value) for thing in theta if not theta[thing].vary]

    ###################################################################################################
    #Set up the initial positions of the walkers as a ball with a different standard deviation in each dimension
    nwalkers=200
    nsteps=30000

    #Get the spread of the starting positions
    stds=[]
    n_general=9
    n_positive=1
    n_emission_lines=7

    #Add in all these standard deviations
    #Kinematic parameters
    stds.extend([100.0, 50.0])
    #General parameters
    stds.extend([0.1]*n_general)
    #Positive parameters
    stds.extend([0.1]*n_positive)
    #Emission lines
    stds.extend([100.0, 50.0])
    stds.extend([1.0]*n_emission_lines)
    #Age
    stds.extend([1.0])
    #Z, imf1, imf2
    stds.extend([0.1, 0.1, 0.1])
    #Sky
    #stds.extend([100.0,  100.0,  100.0, 100.0])
    #ln_f
    stds.extend([0.5])

    stds=np.array(stds)

    assert len(stds)==len(variables), "You must have the same number of dimensions for the Gaussian ball as variables!"

    #Not get the starting values for each parameter, as well as the prior bounds
    start_values, bounds=SF.get_start_vals_and_bounds(theta)
    p0=SF.get_starting_poitions_for_walkers(start_values, stds, nwalkers)
    #Check everything is within the bounds
    #Make sure the positive parameters stay positive
    p0[2+n_general:2+n_general+n_positive, :]=np.abs(p0[2+n_general:2+n_general+n_positive, :])
    #This checks to see if any rows of the array have values which are too high, and replaces them with the upper bound value
    #Add the machine epsilon to deal with cases where we end up with, for example, one walker set to be -0.20000000000000001 instead of -0.2
    p0[p0<bounds[:, 0, None]]=bounds[np.any(p0<bounds[:, 0, None], axis=1), 0]+10*np.finfo(np.float64).eps
    #And the same for any values which are too low
    p0[p0>bounds[:, 1, None]]=bounds[np.any(p0>bounds[:, 1, None], axis=1), 1]-10*np.finfo(np.float64).eps

    assert np.all((p0>bounds[:, 0, None])&(p0<bounds[:, 1, None])), 'Some walkers are starting in bad places of parameter space!'
    p0=p0.T
    ###################################################################################################

    ###################################################################################################
    #Do the sampling
    #This may take a while!

    #Notice the pool argument here
    sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=[theta, variables, bounds], pool=pool)
    result = sampler.run_mcmc(p0, nsteps)

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

    #get rid of the burn-in
    burnin=nsteps-5000
    samples = sampler.chain[:, burnin:, :].reshape((-1, ndim))
    print("\tDone")

    np.savetxt('samples.txt', samples)

max_stars_repo_name: WillCollier/PyStaff

max_issues_repo_name: WillCollier/PyStaff

max_forks_repo_name: WillCollier/PyStaff

reason: import numpy,import scipy


text:
# Common imports
import numpy as np
import pandas as pd
from math import *
import matplotlib.pyplot as plt
import os

def freefall_acceleration (g, D, v):
    return -g + D*v*v

def freefall_yanalytic (g, D, y_0, t, v_0=0):
    vT = sqrt(g/D)
    return y_0-(vT*vT/g)*log(cosh(g*t/vT))

def freefall_no_drag_acceleration  (g, D, v): #Arguments must be in the same order
    return -g

def freefall_no_drag_yanalytic (g, D, y_0, t, v_0):
    return -0.5*g*t**2 + v_0*t + y_0

class ClassicalMechanicsSolvers ():
    def __init__(self, g, D, DeltaT, tfinal, v_0, y_0):
        self.g = g
        self.D = D
        self.set_up_time (DeltaT, tfinal)
        self.set_up_initial_arrays (v_0, y_0)

    ###############################################
    # SETTING UP TIME
    def set_up_time (self, DeltaT, tfinal):
        self.DeltaT = DeltaT
        self.n = ceil(tfinal/DeltaT)
        self.t = np.zeros(self.n)
    ###############################################

    ###############################################
    # SETTING UP THE INITIAL ARRAYS
    def set_up_initial_arrays (self, v_0, y_0): ## No longer need n as an argument
        # set up arrays for t, a, v, and y and we can compare our results with analytical ones
        self.a = np.zeros(self.n)
        self.v = np.zeros(self.n)
        self.y = np.zeros(self.n)
        self.yanalytic = np.zeros(self.n)
        # Initial conditions
        self.v[0] = v_0  #m/s
        self.y[0] = y_0 #m
        self.yanalytic[0] = self.y[0]
    ###############################################

    ###############################################
    # EULER'S METHOD
    def euler (self, min_position, acceleration_eq, yanalytic_eq):
        i_stop = self.n
        # Start integrating using Euler's method
        for i in range(self.n-1):
            # expression for acceleration
            self.a[i] = acceleration_eq (self.g, self.D, self.v[i])
            # update velocity and position
            self.y[i+1] = self.y[i] + self.DeltaT*self.v[i]
            self.v[i+1] = self.v[i] + self.DeltaT*self.a[i]
            # update time to next time step and compute analytical answer
            self.t[i+1] = self.t[i] + self.DeltaT
            self.yanalytic[i+1] = yanalytic_eq(self.g, self.D, self.y[0], self.t[i+1], self.v[0])
            if ( self.y[i+1] < min_position):
                stop_i = i+2
                break
        if stop_i != self.n:
            self.t = self.t[0:stop_i]
            self.a = self.a[0:stop_i]
            self.v = self.v[0:stop_i]
            self.y = self.y[0:stop_i]
            self.yanalytic = self.yanalytic[0:stop_i]
        self.a[-1] = acceleration_eq (self.g, self.D, self.v[-1])
    ###############################################

    ###############################################
    # DISPSLAY THE DATA
    def display_data (self):
        data = {'t[s]': self.t,
                'y[m]': self.y-self.yanalytic,
                'v[m/s]': self.v,
                'a[m/s^2]': self.a
                }
        NewData = pd.DataFrame(data)
        display(NewData.head())
    ###############################################

    ###############################################
    # GRAPH THE DATA
    #finally we plot the data
    def graph_data (self):
        fig, axs = plt.subplots(3, 1)
        axs[0].plot(self.t, self.y, self.t, self.yanalytic)
        axs[0].set_ylabel('y and exact')
        axs[1].plot(self.t, self.v)
        axs[1].set_ylabel('v[m/s]')
        axs[2].plot(self.t, self.a)
        axs[2].set_xlabel('time[s]')
        axs[2].set_ylabel('a[m/s^2]')
        fig.tight_layout()
        plt.show()
    ###############################################


text:
import numpy as np

from asym_rlpo.data import Episode, EpisodeBuffer

def make_episode(num_timesteps: int):
    return Episode(
        states=np.zeros(num_timesteps),
        observations=np.zeros(num_timesteps),
        actions=np.zeros(num_timesteps, dtype=int),
        rewards=np.zeros(num_timesteps),
    )

def test_episode_buffer_max_timesteps():
    episode_buffer = EpisodeBuffer(100)
    assert episode_buffer.num_interactions() == 0
    assert episode_buffer.num_episodes() == 0

    episode_buffer.append_episode(make_episode(101))
    assert episode_buffer.num_interactions() == 0
    assert episode_buffer.num_episodes() == 0

    episode_buffer.append_episode(make_episode(10))
    assert episode_buffer.num_interactions() == 10
    assert episode_buffer.num_episodes() == 1

    episode_buffer.append_episode(make_episode(20))
    assert episode_buffer.num_interactions() == 30
    assert episode_buffer.num_episodes() == 2

    episode_buffer.append_episode(make_episode(80))
    assert episode_buffer.num_interactions() == 100
    assert episode_buffer.num_episodes() == 2

    episode_buffer.append_episode(make_episode(1))
    assert episode_buffer.num_interactions() == 81
    assert episode_buffer.num_episodes() == 2

max_stars_repo_name: abaisero/asym-porl

max_issues_repo_name: abaisero/asym-porl

max_forks_repo_name: abaisero/asym-porl


text:
#!/usr/bin/env python3

#numpyをインポート（行列を扱う各種関数を含む）
import numpy as np

#可視化のためにmatplotlibモジュールをインポート
import matplotlib.pyplot as plt

#アニメーション用のモジュールをインポート
import matplotlib.animation as animation

#確率密度関数の描画用
from scipy.stats import norm

#混合ガウス分布に基づく乱数を生成する
np.random.seed(0)

#各ガウス分布のサンプル数
n_samples=[200,400,400]

#各ガウス分布の平均
means=[-2,3,5]

#各ガウス分布の分散
sigmas=[3,2,0.5]

#GMMに従う乱数を生成
x=None
for n,mean,sigma in zip(n_samples,means,sigmas):
    samples_for_each_dist=np.random.normal(mean,sigma, int(n))
    if x is None:
        x=samples_for_each_dist
    else:
        x=np.concatenate((x,samples_for_each_dist))

#モデルパラメータを初期化する
alpha=np.array([1./3.,1./3.,1./3.])
var=np.array([1.,1.,1.])
mu=np.array([-1,0,1])

#GMMを構成するガウス分布の数
n_clusters=len(alpha)

#繰り返し計算でパラメータを最適化する(ここでは100回繰り返す）
n_iterations=101
log_likelihood=np.zeros(n_iterations)
ims=[]
for t in range(n_iterations):
    print("t{}".format(t))
    if t==0:
        alpha_buf=alpha[None,:]
        var_buf=var[None,:]
        mu_buf=mu[None,:]
    else:
        alpha_buf=np.concatenate((alpha_buf,alpha[None,:]),axis=0)
        var_buf=np.concatenate((var_buf,var[None,:]),axis=0)
        mu_buf=np.concatenate((mu_buf,mu[None,:]),axis=0)
    #Eステップ

    #係数部
    coef=alpha/np.sqrt(2.*np.pi*var)

    #exponent: n_sample, n_clusters
    exponent=-1.*np.power(x[:,None]-mu[None,:],2.)/(2*var[None,:])

    #βを求める
    beta=coef[None,:]*np.exp(exponent)
    likelihood_each_sample=np.maximum(np.sum(beta,axis=1,keepdims=True),1.e-18)
    beta=beta/likelihood_each_sample

    #対数尤度を求める
    current_log_likelihood=np.average(np.log(likelihood_each_sample))
    log_likelihood[t]=current_log_likelihood

    #Mステップ(パラメータを更新する)
    N=np.maximum(np.sum(beta,axis=0),1.e-18)

    #事前確率を更新
    alpha=N/np.sum(N)

    #平均値を更新
    mu=np.einsum("ij,i->j",beta,x)/N
    #分散を計算
    var=np.einsum("ij,ij->j",beta,np.power(x[:,None]-mu[None,:],2.))/N
    var=np.maximum(var,1.e-18)

#対数ゆう度をグラフ化する
plt.figure(figsize=(10,4))
plt.plot(np.arange(0,n_iterations,1),log_likelihood,color="black",linewidth=1,label="log likelihood")
plt.xlabel("Number of iterations")
plt.legend()
plt.savefig("./log_likelihood_gmm.png")
plt.show()

#パラメータ更新の様子をアニメーションで表示
def animation_update(t):
    plt.cla()
    plt.hist(x,bins=50,normed=True,label="observed samples")

    xmin=-20
    xmax=20
    plt.xlim([xmin,xmax])
    plt.ylim([0,0.40])
    x_range=np.arange(-20,20,0.01)

    #GMMを描画する
    pdf=None
    for alpha,var, mu in zip(alpha_buf[t],var_buf[t],mu_buf[t]):
        pdf_each=alpha*norm.pdf(x_range,mu,np.sqrt(var))
        if pdf is None:
            pdf=pdf_each
        else:
            pdf=pdf+pdf_each

    plt.plot(x_range,pdf,color="black",linewidth=1,label=r"$p(x|\theta^{(t="+str(t)+")})$")
    plt.legend()

    if t==0:
        plt.savefig("./initialized_gmm.png")
    if t==n_iterations-1:
        plt.savefig("./learned_gmm.png")

#音声データをプロットする
fig=plt.figure(figsize=(10,4))

ani=animation.FuncAnimation(fig,animation_update,interval=200,frames=n_iterations)
plt.show()


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Tests for nussl/utils.py
"""

import unittest
import nussl
import numpy as np
import os
import librosa

class TestUtils(unittest.TestCase):
    """

    """

    def test_find_peak_indices(self):
        array = np.arange(0, 100)
        peak = nussl.utils.find_peak_indices(array, 1)[0]
        assert peak == 99

        array = np.arange(0, 100).reshape(10, 10)
        peak = nussl.utils.find_peak_indices(array, 3, min_dist=0)
        assert peak == [[9, 9], [9, 8], [9, 7]]

    def test_find_peak_values(self):
        array = np.arange(0, 100)
        peak = nussl.utils.find_peak_values(array, 1)[0]
        assert peak == 99

        array = np.arange(0, 100).reshape(10, 10)
        peak = nussl.utils.find_peak_values(array, 3, min_dist=0)
        assert peak == [99, 98, 97]

    def test_add_mismatched_arrays(self):
        long_array = np.ones((20,))
        short_array = np.arange(10)
        expected_result = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
                                    1, 1, 1, 1, 1, 1, 1, 1, 1, 1], dtype=float)

        # Test basic cases
        result = nussl.utils.add_mismatched_arrays(long_array, short_array)
        assert all(np.equal(result, expected_result))

        result = nussl.utils.add_mismatched_arrays(short_array, long_array)
        assert all(np.equal(result, expected_result))

        expected_result = expected_result[:len(short_array)]

        result = nussl.utils.add_mismatched_arrays(long_array, short_array, truncate=True)
        assert all(np.equal(result, expected_result))

        result = nussl.utils.add_mismatched_arrays(short_array, long_array, truncate=True)
        assert all(np.equal(result, expected_result))

        # Test complex casting
        short_array = np.arange(10, dtype=complex)
        expected_result = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
                                    1, 1, 1, 1, 1, 1, 1, 1, 1, 1], dtype=complex)

        result = nussl.utils.add_mismatched_arrays(long_array, short_array)
        assert all(np.equal(result, expected_result))

        result = nussl.utils.add_mismatched_arrays(short_array, long_array)
        assert all(np.equal(result, expected_result))

        expected_result = expected_result[:len(short_array)]

        result = nussl.utils.add_mismatched_arrays(long_array, short_array, truncate=True)
        assert all(np.equal(result, expected_result))

        result = nussl.utils.add_mismatched_arrays(short_array, long_array, truncate=True)
        assert all(np.equal(result, expected_result))

        # Test case where arrays are equal length
        short_array = np.ones((15,))
        expected_result = short_array * 2

        result = nussl.utils.add_mismatched_arrays(short_array, short_array)
        assert all(np.equal(result, expected_result))

        result = nussl.utils.add_mismatched_arrays(short_array, short_array, truncate=True)
        assert all(np.equal(result, expected_result))

    def test_musdb_track(self):
        signal = nussl.AudioSignal(librosa.util.example_audio_file())

        repet = nussl.Repet(signal)
        repet.run()

        bg, fg = repet.make_audio_signals()

        src_dict = {'vocals': fg, 'accompaniment': bg}
        target = nussl.core.constants.STEM_TARGET_DICT
        track = nussl.utils.audio_signals_to_musdb_track(signal, src_dict, target)

        i = 0

    def test_all_algorithms(self):
        l = len(nussl.separation.all_separation_algorithms)
        i = 0

max_stars_repo_name: KingStorm/nussl

max_issues_repo_name: KingStorm/nussl

max_forks_repo_name: KingStorm/nussl


text:
"""Test code for vision package"""
import numpy as np
import tvm
import topi
import math

from topi.vision import ssd, nms

def verify_roi_align(batch_size, in_channel, size, num_roi, pooled_size, spatial_scale):
    data_shape = (batch_size, in_channel, size, size)
    rois_shape = (num_roi, 5)
    data=tvm.placeholder(data_shape)
    rois=tvm.placeholder(rois_shape)
    np_data = np.random.uniform(size=data_shape).reshape(data_shape).astype('float32') * size
    np_rois = np.random.uniform(size=rois_shape).astype('float32') * size
    np_rois[:, 0] = 0

    def check_device(device):
        ctx = tvm.context(device, 0)
        if not ctx.exist:
            print("Skip because %s is not enabled" % device)
            return
        print("Running on target: %s" % device)
        with tvm.target.create(device):
            out = topi.vision.rcnn.roi_align_v2(data, rois, pooled_size=pooled_size, spatial_scale=spatial_scale)
            s = topi.generic.schedule_roi_align_v2(out)

        tvm_data = tvm.nd.array(np_data, ctx)
        tvm_rois = tvm.nd.array(np_rois, ctx)
        tvm_out = tvm.nd.array(np.zeros((num_roi, in_channel, pooled_size, pooled_size)).astype(out.dtype), ctx=ctx)
        f = tvm.build(s, [data, rois, out], device)
        f(tvm_data, tvm_rois, tvm_out)

        import mxnet
        mx_ctx = mxnet.gpu(0)
        mx_data = mxnet.nd.array(np_data, mx_ctx)
        mx_rois = mxnet.nd.array(np_rois, mx_ctx)
        mx_out = mxnet.nd.contrib.ROIAlign_v2(mx_data, mx_rois, pooled_size=(pooled_size, pooled_size), spatial_scale=spatial_scale)
        mx_out = mx_out.asnumpy()

        tvm_out = tvm_out.asnumpy()

        np.testing.assert_allclose(tvm_out, mx_out, rtol=1e-3)

    for device in ['cuda', 'llvm']:
        check_device(device)

def test_roi_align_v2():
    verify_roi_align(1, 1, 14, 64, 7, 1.)
    verify_roi_align(1, 1, 14, 64, 7, 0.5)

if __name__ == "__main__":
    test_roi_align_v2()


text:
from builtins import zip
import numpy as np
import unittest
import lsst.utils.tests
from sn_fitter.fit_sn_cosmo import Fit_LC
from sn_tools.sn_telescope import Telescope
from sn_fit.process_fit import Fitting
from sn_fit.mbcov import MbCov
import os
import h5py
from astropy.table import Table, vstack

main_repo = 'https://me.lsst.eu/gris/DESC_SN_pipeline'

def getFile(refdir, fname):
    fullname = '{}/{}/{}'.format(main_repo, refdir, fname)

    # check whether the file is available; if not-> get it!
    if not os.path.isfile(fname):
        print('wget path:', fullname)
        cmd = 'wget --no-clobber --no-verbose {}'.format(fullname)
        os.system(cmd)

def getRefDir(dirname):
    fullname = '{}/{}'.format(main_repo, dirname)

    if not os.path.exists(dirname):
        print('wget path:', fullname)
        cmd = 'wget --no-verbose --recursive {} --directory-prefix={} --no-clobber --no-parent -nH --cut-dirs=3 -R \'index.html*\''.format(
            fullname+'/', dirname)
        os.system(cmd)

def getconfig(prodid, simu_name, lc_name, saveData=False):

    conf = {}
    conf['ProductionID'] = prodid

    conf['Instrument'] = {}
    conf['Instrument']['name'] = 'LSST'
    conf['Instrument']['throughput_dir'] = 'LSST_THROUGHPUTS_BASELINE'
    conf['Instrument']['atmos_dir'] = 'THROUGHPUTS_DIR'
    conf['Instrument']['airmass'] = 1.1
    conf['Instrument']['atmos'] = True
    conf['Instrument']['aerosol'] = False

    conf['Simulations'] = {}
    conf['Simulations']['filename'] = simu_name
    conf['Simulations']['dirname'] = '.'

    conf['Fitter'] = {}
    conf['Fitter']['name'] = 'sn_fitter.fit_sn_cosmo'
    conf['Fitter']['model'] = 'salt2-extended'
    conf['Fitter']['version'] = 1.0
    conf['Fitter']['covmb'] = 1

    conf['Display'] = 0

    conf['Output'] = {}
    conf['Output']['directory'] = '.'
    conf['Output']['save'] = saveData

    conf['Multiprocessing'] = {}
    conf['Multiprocessing']['nproc'] = 1

    return conf

class TestSNFitcosmo(unittest.TestCase):

    def testFit(self):

        # fit instance
        telescope = Telescope(airmass=1.2)
        _fit = Fit_LC(telescope=telescope)
        print('test')

        # grab LCs
        simuLC = 'LC_sncosmo_Fake_Fake_DESC_seas_-1_-2.0_0.2.hdf5'

        # get this file from server if necessary
        getFile('unittests', simuLC)

        # open LC file from simulation
        fFile = h5py.File(simuLC, 'r')
        keys = list(fFile.keys())

        # fit LC and store results in data - an astropy Table
        data = Table()
        for key in keys:
            # get LC
            lcsimu = Table.read(fFile, path=key)
            # print(lcsimu.meta)
            # fit LC
            resfit = _fit(lcsimu)
            data = vstack([data, resfit])

        # select a set of results to make the test
        idx = data['z'] < 0.9
        sel = data[idx]
        sel['z'] = np.round(sel['z'], 2)
        ll = np.round(list(np.arange(0.1, 0.9, 0.02)), 2)
        selb = sel[np.in1d(sel['z'], ll)]

        """
        for col in ['z', 'Cov_colorcolor', 'mbfit']:
            print('dictRef[\'{}\']='.format(col), selb[col].tolist())
        """
        dictRef = {}
        dictRef['z'] = [0.1, 0.12, 0.14, 0.16, 0.18, 0.2, 0.22, 0.24, 0.26, 0.28, 0.3, 0.32, 0.34, 0.36, 0.38, 0.4, 0.42, 0.44, 0.46,
                        0.48, 0.5, 0.52, 0.54, 0.56, 0.58, 0.6, 0.62, 0.64, 0.66, 0.68, 0.7, 0.72, 0.74, 0.76, 0.78, 0.8, 0.82, 0.84, 0.86, 0.88]

        dictRef['Cov_colorcolor'] = [7.847411468307449e-07, 1.2926429633578325e-06, 2.182689313770883e-06, 3.4064682540206656e-06, 5.094755132641675e-06, 7.418516965058399e-06, 1.0449342672175223e-05, 1.4530838906902605e-05, 1.9838442223282924e-05, 5.488615265833384e-05, 6.79687434393572e-05, 8.154211056708287e-05, 9.585755607011092e-05, 0.00011249971780490519, 0.0001332468082708582, 0.0001583244598719656, 0.00018974040864400726, 0.00023307069810571821, 0.0002929759524889111,
                                     0.00036750366827063667, 0.0004505526493535665, 0.0005418476331991809, 0.0006556090670946579, 0.0007993634180891086, 0.0009706867250971499, 0.0011684978183072229, 0.0013961226214366633, 0.0054156739068237205, 0.006306687006515156, 0.007213958623759652, 0.00819136400938683, 0.00939627974924506, 0.010844715095289384, 0.01246028538030565, 0.01424202236460317, 0.016338703851064878, 0.018937841331963138, 0.022182703639333848, 0.025822864238779813, 0.0]
        dictRef['mbfit'] = [20.006435534930837, 20.430453337465664, 20.79137803592443, 21.106957852340685, 21.388860253425335, 21.643240843203596, 21.875275946771513, 22.088544601198045, 22.28621053426997, 22.470159114763483, 22.642592920968404, 22.804801863419183, 22.958205505175076, 23.103367585657537, 23.241190478068184, 23.372313980135566, 23.497345309887262, 23.616965168330655, 23.73171763384528,
                            23.842071730391293, 23.948421905233754, 24.051095990149463, 24.150012532855033, 24.245404432237017, 24.337458733240545, 24.426662484018227, 24.513165417734825, 24.597371243206602, 24.679404190398955, 24.758758483003128, 24.83545473841958, 24.909948797457137, 24.982528069838313, 25.053561604365488, 25.12335255171979, 25.191183726268534, 25.257131587610104, 25.32151390441974, 25.384753330779187, -1.0]

        # check results here
        for key, vals in dictRef.items():
            assert(np.isclose(vals, selb[key].tolist()).all())

        # cleaning directory
        if os.path.isfile(simuLC):
            os.system('rm {}'.format(simuLC))

class TestSNFit(unittest.TestCase):
    def testSNFit(self):

        # get data
        # two files requested:
        # Simu_*: astropy table with a list of SN simulation parameters
        # LC_*: astropy tables with LC

        # grab LCs
        prodid = 'sncosmo_Fake_Fake_DESC_seas_-1_-2.0_0.2.hdf5'
        simu_name = 'Simu_{}'.format(prodid)
        lc_name = 'LC_{}'.format(prodid)

        # get these files from server if necessary
        getFile('unittests', simu_name)
        getFile('unittests', lc_name)
        getRefDir('SALT2_Files')

        # build configuration file
        conf = getconfig(prodid, simu_name, lc_name)

        # covmb
        covmb = MbCov('SALT2_Files', paramNames=dict(
            zip(['x0', 'x1', 'color'], ['x0', 'x1', 'c'])))

        # fit instance
        fit = Fitting(conf, covmb=covmb)

        # getting the simu file
        f = h5py.File(simu_name, 'r')
        print(f.keys())
        # reading the simu file
        for i, key in enumerate(f.keys()):
            simul = Table.read(simu_name, path=key)

        simul['z'] = np.round(simul['z'], 2)
        ll = np.round(list(np.arange(0.1, 0.9, 0.1)), 2)
        simul = simul[np.in1d(simul['z'], ll)]

        res = Table()
        for simu in simul:
            lc = None
            lc = Table.read(lc_name, path='lc_{}'.format(simu['id_hdf5']))
            resfit = fit(lc)
            res = vstack([res, resfit])

        idx = res['z'] < 0.8
        sel = res[idx]

        keychecks = ['z', 'Cov_colorcolor', 'mbfit', 'mb_recalc', 'sigma_mu']

        """
        for key in keychecks:
            print('dictRef[\'{}\']='.format(key), res[key].tolist())
        """
        dictRef = {}

        dictRef['z'] = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8]
        dictRef['Cov_colorcolor'] = [7.84149720548657e-07, 7.411367293053516e-06, 6.793756102580619e-05,
                                     0.00015816525343875836, 0.0004502363398936092, 0.001167662290074013, 0.008186441880889383, 0.01633751529759251]
        dictRef['mbfit'] = [20.005795111491146, 21.64206956980012, 22.64256913193122, 23.372165710405174,
                            23.948145769722828, 24.42636698333764, 24.835436823080904, 25.19122934728612]
        dictRef['mb_recalc'] = [20.02731344585068, 21.663516325567763, 22.66404477842765,
                                23.39363924396809, 23.969613792664823, 24.4478277942898, 24.85691319019894, 25.212677718902857]
        dictRef['sigma_mu'] = [0.0025744440149275205, 0.008148971048222674, 0.024267233165841652,
                               0.03934245428396615, 0.06908538342541443, 0.10804914411477154, 0.288358020033676, 0.4019877529973541]

        for key in keychecks:
            assert(np.isclose(dictRef[key], res[key].tolist()).all())

        # cleaning directory
        for fi in [simu_name, lc_name]:
            if os.path.isfile(fi):
                os.system('rm {}'.format(fi))
        for ddir in ['SALT2_Files']:
            if os.path.isdir(ddir):
                os.system('rm -rf {}'.format(ddir))

fit_snfit = TestSNFit
# fit_sncosmo = TestSNFitcosmo

if __name__ == "__main__":
    lsst.utils.tests.init()
    unittest.main(verbosity=5)

# unittest.main(verbosity=5)

reason: import numpy,from astropy


text:
import requests
from skimage import io
import json
import matplotlib.pyplot as plt
import numpy as np
import urllib.request

resp = requests.post("http://localhost:8001/digitreg",
                     files={"file":"hey"})

max_stars_repo_name: ashwani-rathee/Sardor

max_issues_repo_name: ashwani-rathee/Sardor

max_forks_repo_name: ashwani-rathee/Sardor


text:
import wx
import wx.lib.agw.aui as aui
import wx.py as py
import numpy as np
from glsurface import TrackingSurface
from _simxpm import run_xpm, pause_xpm

class RunApp(wx.App):
    def __init__(self):
        wx.App.__init__(self, redirect=False)

    def OnInit(self):
        frame = MainFrame()
        frame.Show(True)
        self.SetTopWindow(frame)
        self.frame = frame
        return True

class MainFrame(wx.Frame):
    ID_SHOW_CONSOLE = wx.NewId()

    def __init__(self):
        wx.Frame.__init__(self, None, -1, 'glsurface demo', size=(410, 432))
        self._mgr = aui.AuiManager()

        # tell AuiManager to manage this frame
        self._mgr.SetManagedWindow(self)
        self.panel = SurfacePanel(self)
        self._mgr.AddPane(self.panel,
                          aui.AuiPaneInfo().Caption("glsurface").CenterPane())
        ns = {}
        ns['wx'] = wx
        ns['app'] = wx.GetApp()
        ns['frame'] = self
        self.shell = py.shell.Shell(self, -1, locals=ns)
        self._mgr.AddPane(
            self.shell,
            aui.AuiPaneInfo().Caption('console').DestroyOnClose(False).Bottom(
            ).Snappable().Dockable().Layer(1).Position(1).MinimizeButton(
                True).MaximizeButton(True).CloseButton(False))
        self._mgr.ShowPane(self.shell, False)
        self._mgr.Update()

        accel = [(wx.ACCEL_CTRL, ord('T'), self.ID_SHOW_CONSOLE)]
        self.accel = wx.AcceleratorTable(accel)
        self.SetAcceleratorTable(self.accel)
        self.Bind(wx.EVT_TOOL, self.OnProcessTool)

    def OnProcessTool(self, event):
        eid = event.GetId()
        if eid == self.ID_SHOW_CONSOLE:
            self._mgr.ShowPane(self.shell, not self.shell.IsShown())
            self._mgr.Update()
        else:
            event.Skip()

class Surface(TrackingSurface):
    def __init__(self, *args, **kwargs):
        TrackingSurface.__init__(self, *args, **kwargs)

    def Initialize(self):
        super(Surface, self).Initialize()
        self.SetShowStepSurface(False)
        self.SetShowMode(mesh=True)
        self.rotate_matrix = np.array([[0.9625753, -0.21669953, 0.16275978],
                                       [0.26339024, 0.88946027, -0.3734787],
                                       [-0.06383575, 0.40237066, 0.91324866]],
                                      dtype=np.float32)

def BitmapFromXPM(xpm):
    xpm_b = [x.encode('utf-8') for x in xpm]
    if 'phoenix' in wx.version():
        return wx.Bitmap(xpm_b)
    else:
        return wx.BitmapFromXPMData(xpm_b)

class SurfacePanel(wx.Panel):
    ID_RUN = wx.NewId()
    ID_PAUSE = wx.NewId()

    def __init__(self, parent):
        wx.Panel.__init__(self, parent, -1)
        self.x = np.linspace(0, 2 * np.pi, 30).reshape((1, 30))
        z = np.cos(self.x).T * np.sin(self.x)
        sizer = wx.BoxSizer(wx.VERTICAL)

        tb = aui.AuiToolBar(self,
                            -1,
                            wx.DefaultPosition,
                            wx.DefaultSize,
                            agwStyle=aui.AUI_TB_OVERFLOW
                            | aui.AUI_TB_PLAIN_BACKGROUND)
        tb.SetToolBitmapSize(wx.Size(16, 16))
        tb.AddSimpleTool(self.ID_RUN, "Run", BitmapFromXPM(run_xpm))
        tb.AddSimpleTool(self.ID_PAUSE, "Pause", BitmapFromXPM(pause_xpm))
        tb.Realize()
        sizer.Add(tb, 0, wx.EXPAND, 0)

        self.canvas = Surface(self, {'z': z})
        sizer.Add(self.canvas, 1, wx.EXPAND | wx.ALL, 5)
        self.SetSizer(sizer)

        self.Bind(wx.EVT_UPDATE_UI, self.OnUpdateTool)
        self.Bind(wx.EVT_TOOL, self.OnProcessTool)

        self.timer = wx.Timer(self)
        self.Bind(wx.EVT_TIMER, self.OnTimer, self.timer)
        #self.timer.Start(100)
        self.is_running = False

    def OnUpdateTool(self, event):
        eid = event.GetId()
        if eid == self.ID_RUN:
            event.Enable(not self.is_running)
        elif eid == self.ID_PAUSE:
            event.Enable(self.is_running)
        else:
            event.Skip()

    def OnProcessTool(self, event):
        eid = event.GetId()
        if eid == self.ID_RUN:
            self.is_running = True
            self.timer.Start(50)
        elif eid == self.ID_PAUSE:
            self.is_running = False
            self.timer.Stop()
        else:
            event.Skip()

    def OnTimer(self, event):
        self.x += 0.01 * 2 * np.pi
        self.x %= 2 * np.pi
        z = np.cos(self.x).T * np.sin(self.x)
        self.canvas.NewFrameArrive(z, False)

def main():
    app = RunApp()
    app.MainLoop()

if __name__ == '__main__':
    main()

max_stars_repo_name: tianzhuqiao/glsurface

max_issues_repo_name: tianzhuqiao/glsurface

max_forks_repo_name: tianzhuqiao/glsurface


text:
#!/usr/bin/env python

import matplotlib.pyplot as plt
import numpy as np
import os

# enable 3d plots
from mpl_toolkits.mplot3d import Axes3D

def list_paths():
    for f in os.listdir('paths'):
        yield f.split('.')[0], np.loadtxt('paths/{}'.format(f))

def main():
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    ax.plot(*np.loadtxt('paths/circle.txt').T)
    ax.plot(*np.loadtxt('paths/triangle.txt').T)
    ax.plot(*np.loadtxt('paths/square.txt').T)
    for name, path in list_paths():
        if np.random.random() < 0.03:
            ax.plot(*path.T)
    plt.show()

if __name__ == '__main__':
    main()

max_stars_repo_name: EmbodiedCognition/tracing-experiment

max_issues_repo_name: EmbodiedCognition/tracing-experiment

max_forks_repo_name: EmbodiedCognition/tracing-experiment


text:
import numpy as np
import pandas as pd

import platform
if platform.system() == 'Windows':
    split_by = "\\"
else:
    split_by = "/"
# Note: Modern windows filesystems generally allow both / and \ in their paths and hence they can be used
# interchangeably, but for the few systems that don't; the above check is necessary.

import os
cur = os.path.abspath(__file__)
cur = cur.rsplit(split_by, 1)[0]

df = pd.read_csv(cur+split_by+"vehicle_data.csv")
car_list = list(df['VEHICLE NAME'])
range_list = list(df['DRIVING RANGE (IN KM)'])

def RemainingRange(battery_level):
    print('Please select the Model: \n')
    for car in car_list:
        print(str(car_list.index(car)+1) + "\t"+car+"\n")
    choice = int(input('Enter Here: '))
    car = car_list[choice-1]
    rang = range_list[choice-1]
    remaining_range = (rang*battery_level)/100

    return remaining_range, (car, rang)

def getRange(rem_battery):
    rem, deets = RemainingRange(rem_battery)
    return rem, deets


text:
"""
Module with the parent abstract class DataManagement. \n
Carmine Schipani, 2021
"""

from abc import ABC, abstractmethod
from OpenSeesPyAssistant.ErrorHandling import *
import numpy as np

class DataManagement(ABC):
    """
    Abstract parent class for data management.
    Using the associated MATLAB class \n
    LOAD_CLASS.m \n
    for the postprocessing in MATLAB, allowing for simpler and more reliable data management because the parameters
    from the OpenSeesPy analysis are imported automatically.
    """

    def SaveData(self, f):
        """
        Function that lists in the command window and saves in a opened file text "f" the data from the "self" class that calls it.
        Example: call this function after this line: \n
        with open(FileName, 'w') as f:

        @param f (io.TextIOWrapper): Opened file to write into

        @exception WrongDimension: The number of lists in the list self.data needs to be 2
        """
        if len(self.data[0]) != 2: raise WrongDimension()

        delimiter = "##############################" # 30 times #
        col_delimiter = "\t" # tab
        for data_line in self.data:
            f.write('\n')
            for col in data_line:
                if type(col) == np.ndarray:
                    tmp_str = np.array_str(col, max_line_width = np.inf)
                else:
                    tmp_str = str(col)
                f.write(tmp_str)
                f.write(col_delimiter)
        f.write('\n')
        f.write('NEW INFO SECTION DELIMITER \t')
        f.write(delimiter)

    @abstractmethod
    def ShowInfo(self):
        """
        Abstract method that shows the data stored in the class in the command window.
        In some cases, it's possible to plot some information (for example the curve of the material model).
        """
        pass

    @abstractmethod
    def ReInit(self):
        """
        Abstract method that computes the value of the parameters with respect of the arguments. \n
        Use after changing the value of argument inside the class (to update the values accordingly). \n
        This function can be very useful in combination with the function "deepcopy()" from the module "copy". \n
        Be careful that the parameter self.Initialized is also copied, thus it is safer to copy the class before the method that calls the actual OpenSees commands (and initialise the object).
        """
        pass

    @abstractmethod
    def UpdateStoredData(self):
        """
        Abstract method used to define and update the self.data member variable. \n
        This member variable (self.data) is a list of lists with 2 entries (info_name and info_value)
            and for each list is stored a different member variable of the class. \n
        Useful to debug the model, export data, copy object.
        """
        pass

max_stars_repo_name: DonCammne/OpenSeesPyAssistant

max_issues_repo_name: DonCammne/OpenSeesPyAssistant

max_forks_repo_name: DonCammne/OpenSeesPyAssistant


text:
import pandas as pd
from pandas.testing import assert_frame_equal
import numpy as np
from madic import interference
from scipy import signal

class TestInterference(object):

    def setup_method(self):
        # all pass_transition_ratio have been set to False for testing purposes
        self.df = pd.DataFrame({'rt_start': [3]*3,
                                'rt_end': [17]*3,
                                'times_arr': [np.arange(21)]*3,
                                'ref_tr': [0.2, 0.3, 0.4],
                                'pass_transition_ratio': [False, False, False]}
                               )

    def test_noise_should_not_be_interference(self):
        # random noise is not matrix interference
        intensities = [np.random.rand(21),
                       np.random.rand(21),
                       np.random.rand(21)]

        self.df['intensities_arr'] = intensities
        self.df['transition_ratio'] = [0.33, 0.33, 0.33]
        self.df['area'] = [3000, 3000, 3000]

        expected_frame = self.df.copy()
        expected_frame['interference'] = False

        result = interference.identify_interference(self.df)

        assert_frame_equal(result, expected_frame)

    def test_identify_interference(self):
        # 3 transitions, 1st has large interference
        intensities = [signal.gaussian(21, std=3)*1000,
                       signal.gaussian(21, std=3),
                       signal.gaussian(21, std=3)]

        self.df['intensities_arr'] = intensities
        self.df['transition_ratio'] = [0.98, 0.01, 0.01]
        self.df['area'] = [50000, 3000, 3000]

        expected_frame = self.df.copy()
        expected_frame['interference'] = [True, False, False]

        result = interference.identify_interference(self.df)

        assert_frame_equal(result, expected_frame)

max_stars_repo_name: dcroote/madic

max_issues_repo_name: dcroote/madic

max_forks_repo_name: dcroote/madic


text:
import numpy as np
import os
import sys
import csv
import json
import sys

class NumpyToJsonEncoder(json.JSONEncoder):
    def default(self, obj):
        if isinstance(obj, np.integer):
            return int(obj)
        elif isinstance(obj, np.floating):
            return float(obj)
        elif isinstance(obj, np.ndarray):
            return obj.tolist()
        else:
            return super(NumpyToJsonEncoder, self).default(obj)

def create_dir(params):

    if  params['problem_type'] == 'TP':

        if params['TP']['dimension'] == 2:
            dir_name = 'tile_planting_{}D_L_{}_p1_{}_p2_{}_p3_{}'.format( params['TP']['dimension'], params['TP']['length'], params['TP']['tile_params'][0], \
                                                                            params['TP']['tile_params'][1], params['TP']['tile_params'][2] )
        else:
            dir_name = 'tile_planting_{}D_L_{}_p2FP_{}_p4FP_{}'.format( params['TP']['dimension'], params['TP']['length'], params['TP']['tile_params'][0], \
                                                                        params['TP']['tile_params'][1] )

    elif  params['problem_type'] == 'WP':
        dir_name = 'wishart_planting_N_{}_alpha_{:.2f}'.format( params['WP']['length'], params['WP']['alpha'] )

    elif  params['problem_type'] == 'DCL':
        dir_name = 'dcl_Lx_{}_Ly_{}_alpha_{:.2f}_R_{}_lambda_{}'.format( params['DCL']['Lx'], params['DCL']['Ly'], params['DCL']['alpha'], params['DCL']['R'], \
                                                                        params['DCL']['lambda'] )

    elif params['problem_type'] == 'XORSAT':
        dir_name = 'regular_xorsat_k_{}_N_{}'.format( params['XORSAT']['k'], params['XORSAT']['N'])

    elif params['problem_type'] == 'K_LOCAL':
        dir_name = 'combined_{}_local_N_{}'.format( params['K_LOCAL']['k'], params['K_LOCAL']['total_spins'] )
    else:
        print('An unexpected error occurred.')
        sys.exit()

    if not os.path.exists(dir_name):
        os.makedirs(dir_name)

    return dir_name

def write_instance(dir_name, file_name, bonds, params):

    if params['file_format'] == 'txt':
        write_text_file(os.path.join(dir_name, file_name), bonds)
    else:
        # Write output in json format

        output_format = 'hobo' if params['convert_to_hobo'] else 'ising'

        output = { "problem": {"type": output_format, "terms": []} }

        for coupler in bonds:
            output['problem']['terms'].append( {"w":coupler[-1], "ids":coupler[:-1] } )

        json_output = json.dumps(output, cls=NumpyToJsonEncoder)

        with open(os.path.join(dir_name, file_name), "w") as the_file:
            the_file.write(json_output)

def write_gs_info(dir_name, gs_info):

    write_text_file(os.path.join(dir_name, 'gs_energies.txt'), gs_info)

def write_text_file(path_to_file, records):

    with open(path_to_file, "w") as the_file:
        csv.register_dialect("custom", delimiter="\t", skipinitialspace=True)
        writer = csv.writer(the_file, dialect="custom")

        for row in records:
            writer.writerow(['{:.7e}'.format(entry) if isinstance(entry, float) else entry for entry in row])


text:
from collections import defaultdict
import numpy as np
import random

def ZERO():
    return 0

def train_test_set(data):
    test_users = set(random.sample(list(data['user_id'].unique()), int(len(data['user_id'].unique()) * 0.2)))
    data_test = data[data['user_id'].isin(test_users)]
    data_train = data[~data['user_id'].isin(test_users)]
    return data_train, data_test

class Model:

    @staticmethod
    def from_name(model_name, *args, **kwargs):
        for model_class in Model._transitive_subclasses():
            if model_class.__name__ == model_name:
                return model_class(*args, **kwargs)
        raise Exception('There is no model with name {}.'.format(model_name))

    @staticmethod
    def _transitive_subclasses(given_class=None):
        if given_class is None:
            given_class = Model
        here = given_class.__subclasses__()
        return here + [sc for scs in [Model._transitive_subclasses(c) for c in here] for sc in scs]

class ItemAverage(Model):

    def __init__(self):
        self._nums = defaultdict(ZERO)
        self._corrects = defaultdict(ZERO)

    def train(self, trainset):
        predicted = np.zeros(len(trainset))
        for i, (p_recall, lexeme_id) in enumerate(trainset[['p_recall', 'lexeme_id']].values):
            predicted[i] = self.predict(lexeme_id)
            self._nums[lexeme_id] += 1
            self._corrects[lexeme_id] += p_recall
        return predicted

    def predict(self, lexeme_id):
        if lexeme_id not in self._nums:
            return 0.5
        return self._corrects[lexeme_id] / self._nums[lexeme_id]

class Elo(Model):

    def __init__(self, alpha=1.1, beta=0.09):
        self._skill = defaultdict(ZERO)
        self._difficulty = defaultdict(ZERO)
        self._count = defaultdict(ZERO)
        self._alpha = alpha
        self._beta = beta

    def train(self, trainset):
        predicted = np.zeros(len(trainset))
        for i, (p_recall, lexeme_id, user_id) in enumerate(trainset[['p_recall', 'lexeme_id', 'user_id']].values):
            predicted[i] = self.predict(user_id, lexeme_id, p_recall)
        return predicted

    def predict(self, user_id, lexeme_id, p_recall=None):
        prediction = _sigmoid(self._skill[user_id] - self._difficulty[lexeme_id])
        if p_recall is not None and (user_id, lexeme_id) not in self._count:
            item_alpha = self._alpha / (1 + self._beta * self._count[lexeme_id])
            user_alpha = self._alpha / (1 + self._beta * self._count[user_id])
            self._count[user_id] += 1
            self._count[lexeme_id] += 1
            self._count[user_id, lexeme_id] = 1
            self._difficulty[lexeme_id] -= item_alpha * (p_recall - prediction)
            self._skill[user_id] += user_alpha * (p_recall - prediction)
        return prediction

    def difficulty(self, lexeme_id):
        return self._difficulty[lexeme_id]

    def skill(self, user_id):
        return self._skill[user_id]

def _sigmoid(x):
    return 1.0 / (1 + np.exp(-x))

max_stars_repo_name: ihsgnef/duolingo-halflife-regression

max_issues_repo_name: ihsgnef/duolingo-halflife-regression

max_forks_repo_name: ihsgnef/duolingo-halflife-regression


text:
import pandas as pd
import numpy as np
import plotly.graph_objects as go
import plotly.express as px
import yfinance as yf
import matplotlib.pyplot as plt
from datetime import datetime
import datetime as dTime

def smaChart(df):
    tempDf = df.copy()
    tempDf['SMA5'] = tempDf.Close.rolling(5).mean()
    tempDf['SMA20'] = tempDf.Close.rolling(20).mean()
    tempDf['SMA50'] = tempDf.Close.rolling(50).mean()
    fig = go.Figure(data=[go.Ohlc(x=tempDf['Date'],
                open=tempDf['Open'],
                high=tempDf['High'],
                low=tempDf['Low'],
                close=tempDf['Close'], name = "OHLC"),
                go.Scatter(x=tempDf.Date, y=tempDf.SMA5, line=dict(color='orange', width=1), name="SMA5"),
                go.Scatter(x=tempDf.Date, y=tempDf.SMA20, line=dict(color='green', width=1), name="SMA20"),
                go.Scatter(x=tempDf.Date, y=tempDf.SMA50, line=dict(color='blue', width=1), name="SMA50")])
    newFig = emaChart(tempDf)
    return fig, newFig

def emaChart(tempDf):
    tempDf['EMA5'] = tempDf.Close.ewm(span=5, adjust=False).mean()
    tempDf['EMA20'] = tempDf.Close.ewm(span=20, adjust=False).mean()
    newFig = go.Figure(data=[go.Ohlc(x=tempDf['Date'],
                open=tempDf['Open'],
                high=tempDf['High'],
                low=tempDf['Low'],
                close=tempDf['Close'], name = "OHLC"),
                go.Scatter(x=tempDf.Date, y=tempDf.SMA5, line=dict(color='orange', width=1), name="SMA5"),
                go.Scatter(x=tempDf.Date, y=tempDf.SMA20, line=dict(color='green', width=1), name="SMA20"),
                go.Scatter(x=tempDf.Date, y=tempDf.SMA50, line=dict(color='blue', width=1), name="SMA50"),
                go.Scatter(x=tempDf.Date, y=tempDf.EMA5, line=dict(color='orange', width=1), name="EMA5"),
                go.Scatter(x=tempDf.Date, y=tempDf.EMA20, line=dict(color='blue', width=1), name="EMA50")])
    return newFig

def yesterday(currency):
    theFig, yesFig = yesterdaySma(currency)
    return theFig, yesFig

def yesterdaySma(currency):
    newDf = yf.download(str(currency),start=str(dTime.datetime.now().date() - dTime.timedelta(days = 1)),end=str(dTime.datetime.now().date()),interval="5m")
    newDf = newDf.reset_index()
    newDf['SMA5'] = newDf.Close.rolling(5).mean()
    newDf['SMA20'] = newDf.Close.rolling(20).mean()
    newDf['SMA50'] = newDf.Close.rolling(50).mean()

    theFig = go.Figure(data=[go.Ohlc(x=newDf['Datetime'],
                open=newDf['Open'],
                high=newDf['High'],
                low=newDf['Low'],
                close=newDf['Close'], name = "OHLC"),
                go.Scatter(x=newDf.Datetime, y=newDf.SMA5, line=dict(color='orange', width=1), name="SMA5"),
                go.Scatter(x=newDf.Datetime, y=newDf.SMA20, line=dict(color='green', width=1), name="SMA20"),
                go.Scatter(x=newDf.Datetime, y=newDf.SMA50, line=dict(color='blue', width=1), name="SMA50")])
    theFig.update_layout(title="SMA Chart")
    yesFig = yesterdayEma(newDf)
    return theFig, yesFig

def yesterdayEma(newDf):
    newDf['EMA5'] = newDf.Close.ewm(span=5, adjust=False).mean()
    newDf['EMA20'] = newDf.Close.ewm(span=20, adjust=False).mean()
    yesFig = go.Figure(data=[go.Ohlc(x=newDf['Datetime'],
                open=newDf['Open'],
                high=newDf['High'],
                low=newDf['Low'],
                close=newDf['Close'], name = "OHLC"),
                go.Scatter(x=newDf.Datetime, y=newDf.EMA5, line=dict(color='orange', width=1), name="EMA5"),
                go.Scatter(x=newDf.Datetime, y=newDf.EMA20, line=dict(color='blue', width=1), name="EMA50")])
    yesFig.update_layout(title="EMA Chart")
    return yesFig

def weeklyChart(dfInfo, currency):
    lstwk = datetime.today() - dTime.timedelta(days=7)
    lstwk = lstwk.strftime("%Y-%m-%d")
    weekDf = yf.download(str(currency),start=str(lstwk),interval="1h")
    weekDf = weekDf.reset_index()
    weeklyFig = px.line(weekDf, x="index", y=["Open", "High", "Adj Close"], title=f"Last Week Performance of {dfInfo['shortName']}")
    return weeklyFig

def monthlyChart(dfInfo, currency):
    lstmon = datetime.today() - dTime.timedelta(days=28)
    lstmon = lstmon.strftime("%Y-%m-%d")
    monthDf = yf.download(str(currency),start=str(lstmon))
    monthDf = monthDf.reset_index()
    monthlyFig = px.line(monthDf, x="Date", y=["Open", "High", "Adj Close"], title=f"Last Month Performance of {dfInfo['shortName']}")
    return monthlyFig

max_stars_repo_name: codervignesh/Forex-Analysis

max_issues_repo_name: codervignesh/Forex-Analysis

max_forks_repo_name: codervignesh/Forex-Analysis


text:
# Copyright 2021 Xilinx Inc.
#
# 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.

# --------------------------------------------------------
# KL loss
# Licensed under The Apache-2.0 License [see LICENSE for details]
# Written by jiashu https://github.com/jiashu-zhu/
# --------------------------------------------------------

"""
KL loss
"""

import mxnet as mx
import numpy as np

class KLLossOperator(mx.operator.CustomOp):
    def __init__(self):
        super(KLLossOperator, self).__init__()

    def forward(self, is_train, req, in_data, out_data, aux):
        bbox_diff = in_data[0].asnumpy()
        self.bbox_diff = bbox_diff
        #print(np.max(bbox_diff))
        alpha = in_data[1].asnumpy()
        self.alpha = alpha
        #print(np.max(alpha))

        loss = np.zeros((bbox_diff.shape[0],1))
        for ibatch in range(bbox_diff.shape[0]):
            _bbox_diff = bbox_diff[ibatch]
            _alpha = alpha[ibatch]
            bbox_diff_abs = np.abs(_bbox_diff)
            bbox_diff_sqr = np.power(_bbox_diff,2)
            eq10_where = np.where(bbox_diff_abs > 1)
            eq9_where = np.where(bbox_diff_abs <=1)
            eq10_label = np.zeros_like(bbox_diff_abs)
            eq10_label[eq10_where[0], eq10_where[1]] = 1
            eq9_label = np.zeros_like(bbox_diff_abs)
            eq9_label[eq9_where[0], eq9_where[1]] = 1

            eq9_result = bbox_diff_sqr * eq9_label
            eq9_result /= 2
            eq10_result = bbox_diff_abs * 0.5
            eq10_result = eq10_result * eq10_label
            result1 = eq9_result + eq10_result

            #print(_alpha)
            # idx = np.where(_alpha < -1)
            # _alpha[idx[0], idx[1]] = -1
            alpha_log = _alpha * 0.5
            alpha_neg = _alpha * -1
            alpha_exp = np.exp(alpha_neg)
            #print(alpha_exp)
            result1 = result1 * alpha_exp

            alpha_log_meanr = alpha_log - np.mean(alpha_log, axis=1,keepdims=True)
            log_loss = np.sum(alpha_log_meanr)

            result1_meanr = result1 - np.mean(result1, axis=1,keepdims=True)
            mul_loss = np.sum(result1_meanr)

            klloss = log_loss + mul_loss
            loss[ibatch] = klloss

        #print(loss)
        self.assign(out_data[0], req[0], mx.nd.array(loss))

    def backward(self, req, out_grad, in_data, out_data, in_grad, aux):
        out0 = np.zeros_like(self.bbox_diff)
        out1 = np.zeros_like(self.alpha)

        for ibatch in range(self.bbox_diff.shape[0]):
            _bbox_diff = self.bbox_diff[ibatch]

            _alpha = self.Normalize(self.alpha[ibatch])
            bbox_diff_abs = np.abs(_bbox_diff)
            bbox_diff_sqr = np.power(_bbox_diff, 2)
            eq10_where = np.where(bbox_diff_abs > 1)
            eq9_where = np.where(bbox_diff_abs <= 1)
            eq10_label = np.zeros_like(bbox_diff_abs)
            eq10_label[eq10_where[0], eq10_where[1]] = 1
            eq9_label = np.zeros_like(bbox_diff_abs)
            eq9_label[eq9_where[0], eq9_where[1]] = 1

            grad_alpha_9 = (-1 * 0.5* np.exp(_alpha * -1) * bbox_diff_sqr + 0.5) * eq9_label
            grad_alpha_10 = (-1 * np.exp(_alpha * -1)*(bbox_diff_abs - 0.5) + 0.5) * eq10_label
            grad_alpha = grad_alpha_9 + grad_alpha_10
            out1[ibatch] = grad_alpha/(grad_alpha.shape[1]*10)

            grad_bbox_9 = (-1 *  np.exp(_alpha * -1) * _bbox_diff) * eq9_label
            grad_bbox_10 = (np.exp(_alpha * -1) * bbox_diff_abs) * eq10_label
            grad_bbox = grad_bbox_9 + grad_bbox_10
            out0[ibatch] = grad_bbox/(grad_bbox.shape[1]*10)
        # print(np.max(out0))
        # print(np.max(out1))
        print(out0.shape)
        print(out1.shape)
        self.assign(in_grad[0], req[0], mx.nd.array(out0))
        self.assign(in_grad[1], req[1], mx.nd.array(out1))

    def Normalize(self,data):
        m = np.mean(data)
        mx = np.max(data)
        mn = np.min(data)
        if mx -mn == 0:
            return data - m
        else:
            return (data - m)/(mx - mn)

@mx.operator.register('KLLoss')
class KLLossProp(mx.operator.CustomOpProp):
    def __init__(self):
        super(KLLossProp, self).__init__(need_top_grad=False)

    def list_arguments(self):
        return ['bbox_diff', 'alpha']

    def list_outputs(self):
        return ['output']

    def infer_shape(self, in_shape):
        data_shape = in_shape[0]
        labels_shape = in_shape[1]
        out_shape = [in_shape[0][0],1]
        return [data_shape, labels_shape], [out_shape]

    def create_operator(self, ctx, shapes, dtypes):
        return KLLossOperator()

    def declare_backward_dependency(self, out_grad, in_data, out_data):
        return []

max_stars_repo_name: Xilinx/progressface

max_issues_repo_name: Xilinx/progressface

max_forks_repo_name: Xilinx/progressface