text:
import os
import tempfile
import numpy as np

from veros import Veros, veros_method, settings

yt_start = -39.0
yt_end = 43
yu_start = -40.0
yu_end = 42

class ACC2(Veros):
    """
    A simple global model with a Southern Ocean and Atlantic part
    """
    @veros_method
    def set_parameter(self):
        self.identifier = "acc2_restart_test"

        self.nx, self.ny, self.nz = 30, 42, 15
        self.dt_mom = 4800
        self.dt_tracer = 86400 / 2.
        self.runlen = 86400 * 365

        self.coord_degree = True
        self.enable_cyclic_x = True

        self.congr_epsilon = 1e-12
        self.congr_max_iterations = 5000

        self.enable_neutral_diffusion = True
        self.K_iso_0 = 1000.0
        self.K_iso_steep = 500.0
        self.iso_dslope = 0.005
        self.iso_slopec = 0.01
        self.enable_skew_diffusion = True

        self.enable_hor_friction = True
        self.A_h = (2 * self.degtom) ** 3 * 2e-11
        self.enable_hor_friction_cos_scaling = 1
        self.hor_friction_cosPower = 1

        self.enable_bottom_friction = True
        self.r_bot = 1e-5

        self.enable_implicit_vert_friction = True
        self.enable_tke = True
        self.c_k = 0.1
        self.c_eps = 0.7
        self.alpha_tke = 30.0
        self.mxl_min = 1e-8
        self.tke_mxl_choice = 2

        self.K_gm_0 = 1000.0

        self.enable_eke = True
        self.eke_k_max = 1e4
        self.eke_c_k = 0.4
        self.eke_c_eps = 0.5
        self.eke_cross = 2.
        self.eke_crhin = 1.0
        self.eke_lmin = 100.0
        self.enable_eke_superbee_advection = True
        self.enable_eke_isopycnal_diffusion = True

        self.enable_idemix = True
        self.enable_idemix_hor_diffusion = True
        self.enable_eke_diss_surfbot = True
        self.eke_diss_surfbot_frac = 0.2
        self.enable_idemix_superbee_advection = True

        self.eq_of_state_type = 3

    @veros_method
    def set_grid(self):
        ddz = [50., 70., 100., 140., 190., 240., 290., 340.,
               390., 440., 490., 540., 590., 640., 690.]
        self.dxt[:] = 2.0
        self.dyt[:] = 2.0
        self.x_origin = 0.0
        self.y_origin = -40.0
        self.dzt[:] = ddz[::-1]
        self.dzt[:] *= 1 / 2.5

    @veros_method
    def set_coriolis(self):
        self.coriolis_t[:, :] = 2 * self.omega * np.sin(self.yt[None, :] / 180. * self.pi)

    @veros_method
    def set_topography(self):
        (X, Y) = np.meshgrid(self.xt, self.yt)
        X = X.transpose()
        Y = Y.transpose()
        self.kbot[...] = (X > 1.0) | (Y < -20)

    @veros_method
    def set_initial_conditions(self):
        # initial conditions
        self.temp[:, :, :, 0:2] = ((1 - self.zt[None, None, :] / self.zw[0]) * 15 * self.maskT)[..., None]
        self.salt[:, :, :, 0:2] = 35.0 * self.maskT[..., None]

        # wind stress forcing
        taux = np.zeros(self.ny + 1, dtype=self.default_float_type)
        yt = self.yt[2:self.ny + 3]
        taux = (.1e-3 * np.sin(np.pi * (self.yu[2:self.ny + 3] - yu_start) / (-20.0 - yt_start))) * (yt < -20) \
            + (.1e-3 * (1 - np.cos(2 * np.pi * (self.yu[2:self.ny + 3] - 10.0) / (yu_end - 10.0)))) * (yt > 10)
        self.surface_taux[:, 2:self.ny + 3] = taux * self.maskU[:, 2:self.ny + 3, -1]

        # surface heatflux forcing
        self.t_star = 15 * np.invert((self.yt < -20) | (self.yt > 20)) \
            + 15 * (self.yt - yt_start) / (-20 - yt_start) * (self.yt < -20) \
            + 15 * (1 - (self.yt - 20) / (yt_end - 20)) * (self.yt > 20.)
        self.t_rest = self.dzt[np.newaxis, -1] / (30. * 86400.) * self.maskT[:, :, -1]

        if self.enable_tke:
            self.forc_tke_surface[2:-2, 2:-2] = np.sqrt((0.5 * (self.surface_taux[2:-2, 2:-2] + self.surface_taux[1:-3, 2:-2]))**2
                                                        + (0.5 * (self.surface_tauy[2:-2, 2:-2] + self.surface_tauy[2:-2, 1:-3]))**2)**(1.5)

        if self.enable_idemix:
            self.forc_iw_bottom[:] = 1.0e-6 * self.maskW[:, :, -1]
            self.forc_iw_surface[:] = 0.1e-6 * self.maskW[:, :, -1]

    @veros_method
    def set_forcing(self):
        self.forc_temp_surface[:] = self.t_rest * (self.t_star - self.temp[:, :, -1, self.tau])

    @veros_method
    def set_diagnostics(self):
        pass

    def after_timestep(self):
        pass

class RestartTest(object):
    timesteps = 10

    def __init__(self, backend):
        self.restart_file = tempfile.NamedTemporaryFile(suffix=".h5", delete=False).name

        self.acc_no_restart = ACC2(backend=backend)
        self.acc_no_restart.diskless_mode = True
        self.acc_restart = ACC2(backend=backend)
        self.acc_restart.restart_output_filename = self.restart_file

    def run(self):
        self.acc_no_restart.setup()
        self.acc_no_restart.runlen = self.acc_no_restart.dt_tracer * self.timesteps
        self.acc_no_restart.run()

        self.acc_restart.setup()
        self.acc_restart.runlen = (self.timesteps - 1) * self.acc_no_restart.dt_tracer
        self.acc_restart.run()

        self.acc_restart.restart_input_filename = self.restart_file
        self.acc_restart.setup()
        self.acc_restart.runlen = self.acc_no_restart.time - self.acc_restart.time
        self.acc_restart.run()

        os.remove(self.restart_file)
        return self.test_passed()

    def test_passed(self):
        passed = True
        for setting in settings.SETTINGS:
            s_1, s_2 = (getattr(obj, setting) for obj in (self.acc_no_restart, self.acc_restart))
            if s_1 != s_2:
                print(setting, s_1, s_2)
        for var in sorted(self.acc_no_restart.variables.keys()):
            if "salt" in var:
                continue
            arr_1, arr_2 = (getattr(obj, var) for obj in (self.acc_no_restart, self.acc_restart))
            try:
                arr_1 = arr_1.copy2numpy()
            except AttributeError:
                pass
            try:
                arr_2 = arr_2.copy2numpy()
            except AttributeError:
                pass
            if "psi" in var:
                arr_1 = arr_1[3:-2, 2:-2]
                arr_2 = arr_2[3:-2, 2:-2]
            np.testing.assert_allclose(*self._normalize(arr_1, arr_2), atol=1e-7)
        return passed

    def _normalize(self, *arrays):
        if any(a.size == 0 for a in arrays):
            return arrays
        norm = np.abs(arrays[0]).max()
        if norm == 0.:
            return arrays
        return (a / norm for a in arrays)

def test_restart(backend):
    RestartTest(backend=backend).run()

max_stars_repo_name: madsbk/veros

max_issues_repo_name: madsbk/veros

max_forks_repo_name: madsbk/veros


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

"""
from hypernets.searchers.evolution_searcher import EvolutionSearcher,Population
from hypernets.core.searcher import OptimizeDirection
from hypernets.core.ops import *
from hypernets.core.callbacks import SummaryCallback
import numpy as np

class Test_Evolution():
    def test_population(self):
        population = Population(optimize_direction=OptimizeDirection.Maximize)
        population.append('a', 0)
        population.append('b', 1)
        population.append('c', 2)
        population.append('d', 3)
        population.append('e', 4)
        population.append('f', 5)
        population.append('g', 6)
        population.append('h', 7)
        population.append('i', 8)
        population.append('i', 9)

        b1 = population.sample_best(25, np.random.RandomState(9527))
        assert b1.reward == 8

        population = Population(optimize_direction=OptimizeDirection.Minimize)
        population.append('a', 0)
        population.append('b', 1)
        population.append('c', 2)
        population.append('d', 3)
        population.append('e', 4)
        population.append('f', 5)
        population.append('g', 6)
        population.append('h', 7)
        population.append('i', 8)
        population.append('i', 9)

        b2 = population.sample_best(25, np.random.RandomState(9527))
        assert b2.reward == 0

    def test_eliminate(self):
        population = Population(optimize_direction=OptimizeDirection.Maximize)
        population.append('a', 4)
        population.append('b', 3)
        population.append('c', 2)
        population.append('d', 1)
        population.append('e', 0)
        population.append('f', 5)
        population.append('g', 6)
        population.append('h', 7)
        population.append('i', 8)
        population.append('j', 9)

        eliminates = population.eliminate(2, regularized=True)
        assert eliminates[0].space_sample == 'a' and eliminates[1].space_sample == 'b'

        eliminates = population.eliminate(2, regularized=False)
        assert eliminates[0].space_sample == 'e' and eliminates[1].space_sample == 'd'

    def test_mutate(self):
        def get_space():
            space = HyperSpace()
            with space.as_default():
                id1 = Identity(p1=Int(0, 10), p2=Choice(['a', 'b']))
                id2 = Identity(p3=Real(0., 1.), p4=Bool())(id1)
            return space

        # population = Population(optimize_direction=OptimizeDirection.Maximize)
        population = Population(optimize_direction='max')

        space1 = get_space()
        space1.random_sample()
        assert space1.all_assigned

        space2 = get_space()
        assert not space2.all_assigned
        new_space = population.mutate(space1, space2)

        pv1 = list(space1.get_assigned_param_values().values())
        pv2 = list(space2.get_assigned_param_values().values())

        assert space2.all_assigned
        assert new_space.all_assigned
        assert np.sum([v1 != v2 for v1, v2 in zip(pv1, pv2)]) == 1

    # def test_searcher_with_hp(self):
    #     def get_space():
    #         space = HyperSpace()
    #         with space.as_default():
    #             in1 = Input(shape=(10,))
    #             in2 = Input(shape=(20,))
    #             in3 = Input(shape=(1,))
    #             concat = Concatenate()([in1, in2, in3])
    #             dense1 = Dense(10, activation=Choice(['relu', 'tanh', None]), use_bias=Bool())(concat)
    #             bn1 = BatchNormalization()(dense1)
    #             dropout1 = Dropout(Choice([0.3, 0.4, 0.5]))(bn1)
    #             output = Dense(2, activation='softmax', use_bias=True)(dropout1)
    #         return space
    #
    #     rs = EvolutionSearcher(get_space, 5, 3, regularized=False, optimize_direction=OptimizeDirection.Maximize)
    #     hk = HyperKeras(rs, optimizer='adam', loss='sparse_categorical_crossentropy', metrics=['accuracy'],
    #                     callbacks=[SummaryCallback()])
    #
    #     x1 = np.random.randint(0, 10000, size=(100, 10))
    #     x2 = np.random.randint(0, 100, size=(100, 20))
    #     x3 = np.random.normal(1.0, 100.0, size=(100))
    #     y = np.random.randint(0, 2, size=(100), dtype='int')
    #     x = [x1, x2, x3]
    #
    #     hk.search(x, y, x, y, max_trails=10)
    #     assert hk.get_best_trail()
    #     best_trial = hk.get_best_trail()
    #
    #     estimator = hk.final_train(best_trial.space_sample, x, y)
    #     score = estimator.predict(x)
    #     result = estimator.evaluate(x, y)
    #     assert len(score) == 100
    #     assert result

ext: py

max_stars_repo_name: frankbaul/Hypernets

max_issues_repo_name: frankbaul/Hypernets

max_forks_repo_name: frankbaul/Hypernets

reason: import numpy


text:
#!/usr/bin/env python
# -*- encoding: utf-8 -*-
#-------------------------------------------------------------------------------
'''
@Author  :   {SEASON}
@License :   (C) Copyright 2013-2022, {OLD_IT_WANG}
@Contact :   {shiter@live.cn}
@Software:   PyCharm
@File    :   NLP_DEMO -- Word vector visualization
@Time    :   2020/9/18 0:44
@Desc    :

参考：https://zhuanlan.zhihu.com/p/194263854

tensorboard 可视化：
https://zhuanlan.zhihu.com/p/33786815
使用bert-serving生成词向量并聚类可视化：
https://zhuanlan.zhihu.com/p/110655509

word2vec词向量可视化聚类与Bert词向量可视化聚类对比：
https://zhuanlan.zhihu.com/p/114390048

spacy 插件 ：
https://spacy.io/universe/project/whatlies
'''
#-------------------------------------------------------------------------------
import gensim

from sklearn.decomposition import IncrementalPCA    # 用于最初的降维
from sklearn.manifold import TSNE                   # 用于最终的降维
import numpy as np                                  # 用于数组控制
from plotly.offline import init_notebook_mode, iplot, plot
import plotly.graph_objs as go
import matplotlib.pyplot as plt
import seaborn as sns
#%matplotlib inline
plt.rcParams['font.sans-serif']=['SimHei']  #用来正常显示中文标签
plt.rcParams['axes.unicode_minus']=False    #用来正常显示负号

def reduce_dimensions(model, plot_in_notebook=True):
    fig, ax = plt.subplots(figsize=(30, 30))
    num_dimensions = 2  # 最终的维度数（2维、3维等）
    vectors = []  # 向量空间中的位置
    labels = []  # “跟踪”词汇以便稍后再次标记我们的数据

    for word in list(model.wv.vocab)[:50]:  # 鉴于模型中的4000+词汇，进行可视化的话会耗费大量的计算时间和计算资源，所以笔者仅展示TOP 300词汇
        vectors.append(model.wv[word])
        labels.append(word)  # 将两个列表转换为numpy向量以进行降维

    vectors = np.asarray(vectors)
    labels = np.asarray(labels)  # 使用t-SNE进行降维
    vectors = np.asarray(vectors)
    # logging.info - 这个网站可出售。 - 最佳的Logging 来源和相关信息。('即将开始t-SNE降维处理。这可能会花费些时间...')
    tsne = TSNE(n_components=num_dimensions, n_iter=250, metric='euclidean', random_state=2019)
    vectors = tsne.fit_transform(vectors)
    x_vals = [v[0] for v in vectors]
    y_vals = [v[1] for v in vectors]  # 创建一个 trace
    trace = go.Scatter(x=x_vals, y=y_vals, mode='text', text=labels)
    data = [trace]
  # logging.info - 这个网站可出售。 - 最佳的Logging 来源和相关信息。('词嵌入可视化图谱绘制已经完成！')
    if plot_in_notebook:
      init_notebook_mode(connected=True)
      iplot(data, filename='word-embedding-plot')
    else:
      plot(data, filename='word-embedding-plot.html')

def main():
    model = gensim.models.Word2Vec.load('../word2vector/gensim_temp.txt')
    print(len(model.wv.vocab))

    #reduce_dimensions(model, plot_in_notebook=False)

if __name__=="__main__":
    main()


text:
"""Problem Set 4: Motion Detection"""

import numpy as np
import cv2

# Utility function
def normalize_and_scale(image_in, scale_range=(0, 255)):
    """Normalizes and scales an image to a given range [0, 255].

    Utility function. There is no need to modify it.

    Args:
        image_in (numpy.array): input image.
        scale_range (tuple): range values (min, max). Default set to
                             [0, 255].

    Returns:
        numpy.array: output image.
    """
    image_out = np.zeros(image_in.shape)
    cv2.normalize(image_in, image_out, alpha=scale_range[0],
                  beta=scale_range[1], norm_type=cv2.NORM_MINMAX)

    return image_out

# Assignment code
def gradient_x(image):
    """Computes image gradient in X direction.

    Use cv2.Sobel to help you with this function. Additionally you
    should set cv2.Sobel's 'scale' parameter to one eighth and ksize
    to 3.

    Args:
        image (numpy.array): grayscale floating-point image with
                             values in [0.0, 1.0].

    Returns:
        numpy.array: image gradient in the X direction. Output
                     from cv2.Sobel.
    """
    image_cpy = image.copy()
    # sobel scale = 1/8, ksize = 3 default border type
    sobel_x = cv2.Sobel(image_cpy, cv2.CV_64F, 1, 0, scale=0.125, ksize=3)
    return sobel_x

def gradient_y(image):
    """Computes image gradient in Y direction.

    Use cv2.Sobel to help you with this function. Additionally you
    should set cv2.Sobel's 'scale' parameter to one eighth and ksize
    to 3.

    Args:
        image (numpy.array): grayscale floating-point image with
                             values in [0.0, 1.0].

    Returns:
        numpy.array: image gradient in the Y direction.
                     Output from cv2.Sobel.
    """
    image_cpy = image.copy()
    # sobel scale = 1/8, ksize = 3 default border type
    sobel_y = cv2.Sobel(image_cpy, cv2.CV_64F, 0, 1, scale=0.125, ksize=3)
    return sobel_y

def optic_flow_lk(img_a, img_b, k_size, k_type, sigma=1):
    """Computes optic flow using the Lucas-Kanade method.

    For efficiency, you should apply a convolution-based method.

    Note: Implement this method using the instructions in the lectures
    and the documentation.

    You are not allowed to use any OpenCV functions that are related
    to Optic Flow.

    Args:
        img_a (numpy.array): grayscale floating-point image with
                             values in [0.0, 1.0].
        img_b (numpy.array): grayscale floating-point image with
                             values in [0.0, 1.0].
        k_size (int): size of averaging kernel to use for weighted
                      averages. Here we assume the kernel window is a
                      square so you will use the same value for both
                      width and height.
        k_type (str): type of kernel to use for weighted averaging,
                      'uniform' or 'gaussian'. By uniform we mean a
                      kernel with the only ones divided by k_size**2.
                      To implement a Gaussian kernel use
                      cv2.getGaussianKernel. The autograder will use
                      'uniform'.
        sigma (float): sigma value if gaussian is chosen. Default
                       value set to 1 because the autograder does not
                       use this parameter.

    Returns:
        tuple: 2-element tuple containing:
            U (numpy.array): raw displacement (in pixels) along
                             X-axis, same size as the input images,
                             floating-point type.
            V (numpy.array): raw displacement (in pixels) along
                             Y-axis, same size and type as U.
    """

    # fix to make sure the kernel size is an odd number.
    if k_size % 2 == 0 and k_size > 2:
        k_size = k_size - 1

    # copy images for use.
    image_a = img_a.copy()
    image_b = img_b.copy()

    k = np.ones((k_size, k_size)) / (k_size ** 2)
    # base gradients.
    gx = gradient_x(image_a)
    gy = gradient_y(image_a)
    gt = image_b - image_a

    # calculate the weighted gradients using blurs to prevent having to loop. Makes this much easier to manage.
    if k_type == 'uniform':
        gx_x = cv2.filter2D(gx * gx, ddepth=-1, kernel=k)
        gy_y = cv2.filter2D(gy * gy, ddepth=-1, kernel=k)
        gx_y = cv2.filter2D(gx * gy, ddepth=-1, kernel=k)
        gy_x = cv2.filter2D(gy * gx, ddepth=-1, kernel=k)
        gx_t = cv2.filter2D(gx * gt, ddepth=-1, kernel=k)
        gy_t = cv2.filter2D(gy * gt, ddepth=-1, kernel=k)
    else:
        blur_kernel = (k_size, k_size)
        gx_x = cv2.GaussianBlur(gx * gx, blur_kernel, sigma)
        gy_y = cv2.GaussianBlur(gy * gy, blur_kernel, sigma)
        gx_y = cv2.GaussianBlur(gx * gy, blur_kernel, sigma)
        gy_x = cv2.GaussianBlur(gy * gx, blur_kernel, sigma)
        gx_t = cv2.GaussianBlur(gx * gt, blur_kernel, sigma)
        gy_t = cv2.GaussianBlur(gy * gt, blur_kernel, sigma)

    # create a new matrix from the gradients in the x, the xy, the yx, and the yy directions. To allow the use
    # of matrix multiplier provided by numpy, we flatten, then transpose to the correct shape of a 2x2 matrix.
    # np.matmul states that if either argument is N-D, N > 2, it is treated as a stack of matrices. We want a stack
    # of 2x2 matrices. residing in the last two indexes and broadcast accordingly.

    image_a_transform = np.array([
        gx_x.flatten(), gx_y.flatten(),
        gy_x.flatten(), gy_y.flatten()
    ]).T.reshape(-1, 2, 2)  # transpose and take from last a 2 x 2 matrix in the reshape.

    # 2 x 1 matrix for matmul
    time_transform = np.array([
        -gx_t.flatten(),
        -gy_t.flatten()
    ]).T.reshape(-1, 2, 1)

    try:
        image_a_transform_inv = np.linalg.inv(image_a_transform)
    except np.linalg.linalg.LinAlgError:
        # if this is a singular matrix (happens when the shift is the same as the original)
        return np.zeros(image_a.shape), np.zeros(image_a.shape)

    transform = np.matmul(image_a_transform_inv, time_transform)
    U = transform[:, 0, 0].reshape(image_a.shape)  # get hte U matrix
    V = transform[:, 1, 0].reshape(image_a.shape)  # get the V
    return U, V

def reduce_image(image):
    """Reduces an image to half its shape.

    The autograder will pass images with even width and height. It is
    up to you to determine values with odd dimensions. For example the
    output image can be the result of rounding up the division by 2:
    (13, 19) -> (7, 10)

    For simplicity and efficiency, implement a convolution-based
    method using the 5-tap separable filter.

    Follow the process shown in the lecture 6B-L3. Also refer to:
    -  Burt, P. J., and Adelson, E. H. (1983). The Laplacian Pyramid
       as a Compact Image Code
    You can find the link in the problem set instructions.

    Args:
        image (numpy.array): grayscale floating-point image, values in
                             [0.0, 1.0].

    Returns:
        numpy.array: output image with half the shape, same type as the
                     input image.
    """

    image = image.copy()
    # only select odd columns out.
    image = image.copy()
    # 1/16 4/16  6/16 4/12 1/16
    # k = np.array([1, 4, 6, 4, 1]) / 16.0
    k = np.array([np.array([1, 4, 6, 4, 1]) / 16.0])
    r_k = np.dot(k.T, k)
    filter_img = cv2.filter2D(image, -1, r_k)
    return filter_img[::2, ::2]

def gaussian_pyramid(image, levels):
    """Creates a Gaussian pyramid of a given image.

    This method uses reduce_image() at each level. Each image is
    stored in a list of length equal the number of levels.

    The first element in the list ([0]) should contain the input
    image. All other levels contain a reduced version of the previous
    level.

    All images in the pyramid should floating-point with values in

    Args:
        image (numpy.array): grayscale floating-point image, values
                             in [0.0, 1.0].
        levels (int): number of levels in the resulting pyramid.

    Returns:
        list: Gaussian pyramid, list of numpy.arrays.
    """
    image = image.copy()
    images = [image]   # the list for the pyramid. has base image first.
    for i in range(levels-1):
        image = reduce_image(image)  # reduce it then append. n times.
        images.append(image)
    return images

def create_combined_img(img_list):
    """Stacks images from the input pyramid list side-by-side.

    Ordering should be large to small from left to right.

    See the problem set instructions for a reference on how the output
    should look like.

    Make sure you call normalize_and_scale() for each image in the
    pyramid when populating img_out.

    Args:
        img_list (list): list with pyramid images.

    Returns:
        numpy.array: output image with the pyramid images stacked
                     from left to right.
    """

    # height is constant. image will never be larger than the first.
    h = img_list[0].shape[0]
    # width will be sum of all images size in the x direction
    w = sum([i.shape[1] for i in img_list])

    output = np.zeros((h, w))  # empty image.
    curr_x = 0
    for image in img_list:
        ih, iw = image.shape  # use this to determine where to place.
        output[: ih, curr_x: curr_x + iw] = normalize_and_scale(image)
        curr_x += iw

    return output

def expand_image(image):
    """Expands an image doubling its width and height.

    For simplicity and efficiency, implement a convolution-based
    method using the 5-tap separable filter.

    Follow the process shown in the lecture 6B-L3. Also refer to:
    -  Burt, P. J., and Adelson, E. H. (1983). The Laplacian Pyramid
       as a Compact Image Code

    You can find the link in the problem set instructions.

    Args:
        image (numpy.array): grayscale floating-point image, values
                             in [0.0, 1.0].

    Returns:
        numpy.array: same type as 'image' with the doubled height and
                     width.
    """
    image = image.copy()
    dr = 2 * image.shape[0]  # double rows.
    dc = 2 * image.shape[1]  # double columns

    output = np.zeros((dr, dc))
    # # filling alternate rows
    output[::2, ::2] = image
    # 2/16, 8/16, 12/16, 8/16, 2/16
    k = np.array([np.array([2, 8, 12, 8, 2]) / 16.0])
    r_k = np.dot(k.T, k)
    output = cv2.filter2D(output, -1, r_k)
    return output

def laplacian_pyramid(g_pyr):
    """Creates a Laplacian pyramid from a given Gaussian pyramid.

    This method uses expand_image() at each level.

    Args:
        g_pyr (list): Gaussian pyramid, returned by gaussian_pyramid().

    Returns:
        list: Laplacian pyramid, with l_pyr[-1] = g_pyr[-1].
    """
    images = []
    for idx in range(len(g_pyr)):
        if idx == len(g_pyr) - 1:
            image = g_pyr[idx]
        else:
            h = g_pyr[idx].shape[0]
            w = g_pyr[idx].shape[1]
            # expand the image that is found to the next images size.
            image = g_pyr[idx] - expand_image(g_pyr[idx + 1])[:h, :w]
        images.append(image)

    return images

def warp(image, U, V, interpolation, border_mode):
    """Warps image using X and Y displacements (U and V).

    This function uses cv2.remap. The autograder will use cubic
    interpolation and the BORDER_REFLECT101 border mode. You may
    change this to work with the problem set images.

    See the cv2.remap documentation to read more about border and
    interpolation methods.

    Args:
        image (numpy.array): grayscale floating-point image, values
                             in [0.0, 1.0].
        U (numpy.array): displacement (in pixels) along X-axis.
        V (numpy.array): displacement (in pixels) along Y-axis.
        interpolation (Inter): interpolation method used in cv2.remap.
        border_mode (BorderType): pixel extrapolation method used in
                                  cv2.remap.

    Returns:
        numpy.array: warped image, such that
                     warped[y, x] = image[y + V[y, x], x + U[y, x]]
    """

    image = image.copy()

    # displacements in x-axis and y-axis as U and V for image as float32
    U = U.astype(np.float32)
    V = V.astype(np.float32)

    new_x, new_y = np.meshgrid(range(image.shape[1]), range(image.shape[0]))
    new_x = new_x.astype(np.float32) + U
    new_y = new_y.astype(np.float32) + V

    output = cv2.remap(src=image, map1=new_x, map2=new_y, interpolation=interpolation, borderMode=border_mode)
    return output

def hierarchical_lk(img_a, img_b, levels, k_size, k_type, sigma, interpolation,
                    border_mode):
    """Computes the optic flow using Hierarchical Lucas-Kanade.

    This method should use reduce_image(), expand_image(), warp(),
    and optic_flow_lk().

    Args:
        img_a (numpy.array): grayscale floating-point image, values in
                             [0.0, 1.0].
        img_b (numpy.array): grayscale floating-point image, values in
                             [0.0, 1.0].
        levels (int): Number of levels.
        k_size (int): parameter to be passed to optic_flow_lk.
        k_type (str): parameter to be passed to optic_flow_lk.
        sigma (float): parameter to be passed to optic_flow_lk.
        interpolation (Inter): parameter to be passed to warp.
        border_mode (BorderType): parameter to be passed to warp.

    Returns:
        tuple: 2-element tuple containing:
            U (numpy.array): raw displacement (in pixels) along X-axis,
                             same size as the input images,
                             floating-point type.
            V (numpy.array): raw displacement (in pixels) along Y-axis,
                             same size and type as U.
    """

    image_a = img_a.copy()
    image_b = img_b.copy()

    gaussian_pyramid_a = gaussian_pyramid(image_a, levels)
    gaussian_pyramid_b = gaussian_pyramid(image_b, levels)

    # init U and V
    U = np.zeros(gaussian_pyramid_a[-1].shape)
    V = np.zeros(gaussian_pyramid_a[-1].shape)

    for aa, bb in reversed(zip(gaussian_pyramid_a, gaussian_pyramid_b)):
        h, w = aa.shape
        # expand the flow field and double it to get the the next image.
        U = (expand_image(U) * 2)[:h, :w]
        V = (expand_image(V) * 2)[:h, :w]
        # warp it
        c = warp(bb, U, V, interpolation=interpolation, border_mode=border_mode)

        # perform Lucas Kanade
        nx, ny = optic_flow_lk(aa, c, k_size, k_type, sigma)

        # add to the orignal flow.
        U = U + nx
        V = V + ny

    return U, V

max_stars_repo_name: jperuggia/ComputerVision

max_issues_repo_name: jperuggia/ComputerVision

max_forks_repo_name: jperuggia/ComputerVision


text:
# --------------------------------------------------------
# FCN
# Copyright (c) 2016 RSE at UW
# Licensed under The MIT License [see LICENSE for details]
# Written by Yu Xiang
# --------------------------------------------------------

"""Test a FCN on an imdb (image database)."""

from fcn.config import cfg, get_output_dir
import argparse
from utils.timer import Timer
from utils.blob import im_list_to_blob, pad_im, unpad_im
from utils.voxelizer import Voxelizer, set_axes_equal
from utils.se3 import *
import numpy as np
import cv2
import cPickle
import os
import math
import tensorflow as tf
import scipy.io
import time
from normals import gpu_normals
# from pose_estimation import ransac
#from kinect_fusion import kfusion
import cv2
import scipy.ndimage

import triplet_flow_loss.run_slow_flow_calculator_process
# import pyximport; pyximport.install()
import utils.plotting_tools

#output/lov_synthetic_features_drill_crackerbox/batch_size_4_loss_L2_optimizer_ADAM_network_net_labeled_concat_features_2017-08-06_2017-08-06_19-35-28_909447/vgg16_flow_lov_features_iter_9000.ckpt
#output/lov_synthetic_features/batch_size_1_loss_L2_optimizer_ADAM_network_net_labeled_concat_features_2017-08-06_2017-08-06_19:35:48_853842/vgg16_flow_lov_features_iter_69000.ckpt

def get_network_output(sess,net, blob):
    left_blob = blob['left_image']
    depth_blob = blob['depth']
    left_labels_blob = blob['left_labels']

    network_inputs = {net.data_left: left_blob, net.keep_prob: 1.0}

    network_outputs = [net.get_output('features_1x_l'), net.get_output("features_2x_l"),
                       net.get_output("features_4x_l"), net.get_output("features_8x_l")]
    results = siphon_outputs_single_frame(sess, net, network_inputs, network_outputs)

    features_l = np.concatenate([
        scipy.ndimage.zoom(np.squeeze(results[7]), (8, 8, 1), order=1),
        scipy.ndimage.zoom(np.squeeze(results[6]), (4, 4, 1), order=1),
        scipy.ndimage.zoom(np.squeeze(results[5]), (2, 2, 1), order=1),
        np.squeeze(results[0])], axis=2)

    return features_l

def create_voxel(sess, net, imdb, weights_filename, image_set):

    if weights_filename is not None:
        output_dir = get_output_dir(imdb, weights_filename)
        if not os.path.exists(output_dir):
            os.makedirs(output_dir)

    # get a set of all the images that will go into the voxel
    images = {}
    for pair in imdb.roidb:
        if image_set in pair['video_id']:
            if pair['image'] not in images:
                images[pair['image']] = pair

    voxelizer = Voxelizer(1, 1)

    for image in images.values():
        minibatch = get_batch(image, voxelizer)
        features = get_network_output(sess, net, minibatch)
        pass

def display_img(img, title, right=False, ax=None):
    global iiiiii, x_plots, y_plots, axes_left_list, axes_right_list, fig
    if ax is None:
        ax = fig.add_subplot(y_plots, x_plots, iiiiii)
        iiiiii += 1
        if right:
            axes_right_list.append(ax)
        else:
            axes_left_list.append(ax)
    ax.imshow(img)
    ax.set_title(title)
    return ax

def fix_rgb_image(image_in):
    image = image_in.copy() + cfg.PIXEL_MEANS
    image = image[:, :, (2, 1, 0)]
    image = image.astype(np.uint8)
    return image

def siphon_outputs_single_frame(sess, net, data_feed_dict, outputs):
    # compute image blob

    sess.run(net.enqueue_op, feed_dict=data_feed_dict)
    output = sess.run(outputs)

    # assert sess.run(net.queue_size_op) == queue_start_size, "data queue size changed"
    return output

def get_batch(roidb, voxelizer):
    """Given a roidb, construct a minibatch sampled from it."""

    # Get the input image blob, formatted for tensorflow
    image_left_blob, left_label_blob, depth_blob, rt_matrix, intrinsic_matrix = _get_image_blob(roidb)

    blobs = {'left_image': image_left_blob,
             'left_labels': left_label_blob,
             'depth': depth_blob,
             'rt_matrix': rt_matrix,
             'intrinsic_matrix': intrinsic_matrix}

    return blobs

def _get_image_blob(roidb):
    """Builds an input blob from the images in the roidb at the specified
    scales.
    """
    num_images = len(roidb)
    processed_left = []
    processed_depth = []
    processed_left_labels = []

    im_left_raw = cv2.imread(roidb['image'], cv2.IMREAD_UNCHANGED)
    im_left = pad_im(im_left_raw[:, :, :3], 16).astype(np.float32)
    im_left -= cfg.PIXEL_MEANS
    im_left_processed = cv2.resize(im_left, None, None, fx=1, fy=1,
                                   interpolation=cv2.INTER_LINEAR)
    processed_left.append(im_left_processed)

    # meta data
    # meta_data1 = scipy.io.loadmat(roidb[i*2]['meta_data'])
    K1 = [[1.06677800e+03, 0.00000000e+00, 3.12986900e+02],
           [0.00000000e+00, 1.06748700e+03, 2.41310900e+02],
           [0.00000000e+00, 0.00000000e+00, 1.00000000e+00]]
    K1 = np.matrix(K1)
    rt_mat = np.genfromtxt(roidb['pose'], skip_header=1)

    # depth
    depth_raw = pad_im(cv2.imread(roidb['depth'], cv2.IMREAD_UNCHANGED), 16).astype(np.float32)
    processed_depth.append(depth_raw / 10000)

    if cfg.TRAIN.USE_MASKS:
        # read left label image
        alpha_channel_l = im_left_raw[:, :, 3] / 255.
        im_cls = (np.dstack([np.zeros(im_left_raw.shape[:2]), alpha_channel_l]))
        processed_left_labels.append(im_cls)
    else:
        processed_left_labels.append(np.zeros([im_left_processed.shape[0], im_left_processed.shape[1],2], dtype=np.int32) + [0, 1])

    # Create a blob to hold the input images
    left = im_list_to_blob(processed_left, 3)
    depth = im_list_to_blob(processed_depth, 1)
    if cfg.TRAIN.USE_MASKS:
        left_labels = im_list_to_blob(processed_left_labels, processed_left_labels[0].shape[2])
    else:
        left_labels = im_list_to_blob(processed_left_labels, 2, mess_up_single_channel=True)

    return left, left_labels, depth, K1, rt_mat

# RT is a 3x4 matrix
def se3_inverse(RT):
    R = RT[0:3, 0:3]
    T = RT[0:3, 3].reshape((3, 1))
    RT_new = np.zeros((3, 4), dtype=np.float32)
    RT_new[0:3, 0:3] = R.transpose()
    RT_new[0:3, 3] = -1 * np.dot(R.transpose(), T).reshape((3))
    return RT_new

def se3_mul(RT1, RT2):
    R1 = RT1[0:3, 0:3]
    T1 = RT1[0:3, 3].reshape((3, 1))

    R2 = RT2[0:3, 0:3]
    T2 = RT2[0:3, 3].reshape((3, 1))

    RT_new = np.zeros((3, 4), dtype=np.float32)
    RT_new[0:3, 0:3] = np.dot(R1, R2)
    T_new = np.dot(R1, T2) + T1
    RT_new[0:3, 3] = T_new.reshape((3))
    return RT_new

reason: import numpy,import scipy


text:
import os
import sys
import zipfile
import requests
import shutil
from shapely.geometry.polygon import Polygon
from shapely.ops import transform
import pyproj
from osgeo import gdal
from pandas import DataFrame
import pandas as pd
import numpy as np
from osgeo import gdal
import pickle
from tqdm import tqdm
import click
from threading import Thread
from multiprocessing.dummy import Pool as ThreadPool
import json
import subprocess

# pixel area only depends on latitude (not longitude)
# we re-project WGS84 to cylindrical equal area
def pixel_area(pix_deg):
    project = lambda x, y: pyproj.transform(pyproj.Proj(init='epsg:4326'), pyproj.Proj(proj='cea'), x, y)
    offset = pix_deg / 2
    lts = np.arange(90-offset, -90, -pix_deg)
    area = np.empty_like(lts)
    lon = 0
    for y, lat in enumerate(lts):
        pixel1 = Polygon([(lon - offset, lat + offset), (lon + offset, lat + offset), (lon + offset, lat - offset), (lon - offset, lat - offset)])
        pixel2 = transform(project, pixel1)
        area[y] = pixel2.area
    return area

def get_flow_dir(row):
    if not os.path.exists(f'tiles/dir/3s/{row.tile}'):
        tqdm.write(f'Downloading {row.tile}...')
        r = requests.get(row.url + row.tile)
        with open(f'tiles/dir/3s/{row.tile}', 'wb') as f:
            f.write(r.content)
    try:
        with zipfile.ZipFile(f'tiles/dir/3s/{row.tile}', 'r') as z:
            z.extractall(path = 'tmp/')
        flow_dir = gdal.Open(f'tmp/{row.tile[:-9]}/{row.tile[:-9]}/w001001.adf')
        geo = flow_dir.GetGeoTransform()
        ySize, xSize = flow_dir.RasterYSize, flow_dir.RasterXSize
        flow_dir = flow_dir.ReadAsArray()
        shutil.rmtree(f'tmp/{row.tile[:-9]}')
        # data is padded into a 6000x6000 array (some tiles may be smaller):
        array_5x5 = np.ones((6000, 6000), dtype='uint8') * 255
        y0 = int(round((geo[3] - row.lat) / geo[5]))
        y1 = 6000 - int(round(((row.lat - 5) - (geo[3] + geo[5] * ySize)) / geo[5]))
        x0 = int(round((geo[0] - row.lon) / geo[1]))
        x1 = 6000 - int(round(((row.lon + 5) - (geo[0] + geo[1] * xSize)) / geo[1]))
        array_5x5[y0:y1, x0:x1] = flow_dir
    except:
        tqdm.write('Not a ZIP file!')
        array_5x5 = np.zeros((6000, 6000), dtype='uint8')
    return array_5x5

def pass1(cpu, drop_pixel, pix_area, df):
    for row in df.iterrows():
        process_tile(cpu, drop_pixel, pix_area, row[1], df, True)

def pass2(parallel, drop_pixel, pix_area, df):
    i = 0
    acc_dict = {}
    udlr_dict = {}
    while not np.all(df.done):
        row = df[df.done==False].iloc[0]
        process_tile(0, drop_pixel, pix_area, row, df, False, acc_dict, udlr_dict)
        i += 1
        if (i % 100) == 0:
            print('Compressing tiles...')
            compress_tiles(parallel)
    for k, v in acc_dict.items():
        np.savez_compressed(f'tiles/acc/3s/{k}', a=v)
    compress_tiles(parallel)
    # create TIFF files
    print('Creating TIFF files')
    for row in tqdm(df.iterrows()):
        name = row[1].tile.replace('_dir_grid.zip', '_acc')
        a = np.load(f'tiles/acc/3s/{name}.npz')['a']
        driver = gdal.GetDriverByName('GTiff')
        ds = driver.Create(f'tiles/acc/3s/{name}.tif', a.shape[1], a.shape[0], 1, gdal.GDT_Float64, ['COMPRESS=LZW'])
        ds.SetGeoTransform((row[1].lon, 5/6000, 0.0, row[1].lat, 0.0, -5/6000))
        ds.SetProjection('GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],AUTHORITY[\"EPSG\",\"4326\"]]')
        band = ds.GetRasterBand(1)
        band.WriteArray(a)
        band.SetNoDataValue(0.)
        ds = None
        os.remove(f'tiles/acc/3s/{name}.npz')
    print('Creating VRT (Virtual Dataset)')
    subprocess.check_call('gdalbuildvrt acc.vrt tiles/acc/3s/*.tif', shell=True)

def compress_tiles(parallel):
    paths = []
    for path in ['tmp/udlr', 'tiles/acc/3s']:
        for fname in os.listdir(path):
            if fname.endswith('.npy'):
                paths.append(f'{path}/{fname}')
    pool = ThreadPool(parallel)
    pool.map(compress_tile, paths)

def compress_tile(path):
    a = np.load(path)
    np.savez_compressed(path[:-4], a=a)
    os.remove(path)

def process_tile(cpu, drop_pixel, pix_area, row, df, first_pass, acc_dict={}, udlr_dict={}):
    lat, lon = row['lat'], row['lon']
    flow_dir = get_flow_dir(row)
    name = row['tile'][:-len('_dir_grid.zip')]
    if first_pass:
        udlr_in = np.zeros((4, 6000), dtype='float64')
    else:
        df.loc[(df.lat==lat) & (df.lon==lon), 'done'] = True
        if f'udlr_{lat}_{lon}' in udlr_dict:
            udlr_in = udlr_dict[f'udlr_{lat}_{lon}']
        elif os.path.exists(f'tmp/udlr/udlr_{lat}_{lon}.npz'):
            udlr_in = np.load(f'tmp/udlr/udlr_{lat}_{lon}.npz')['a']
        elif os.path.exists(f'tmp/udlr/udlr_{lat}_{lon}.npy'):
            udlr_in = np.load(f'tmp/udlr/udlr_{lat}_{lon}.npy')
        else:
            print(f'Nothing to do for {name}')
            df.to_pickle('tmp/df.pkl')
            return
    if (not first_pass) and (f'{name}_acc' in acc_dict):
        flow_acc = acc_dict[f'{name}_acc']
    elif os.path.exists(f'tiles/acc/3s/{name}_acc.npz'):
        flow_acc = np.load(f'tiles/acc/3s/{name}_acc.npz')['a']
    elif os.path.exists(f'tiles/acc/3s/{name}_acc.npy'):
        flow_acc = np.load(f'tiles/acc/3s/{name}_acc.npy')
    else:
        flow_acc = np.zeros((6000, 6000), dtype='float64')
    udlr_out = np.zeros((4, 6000+2), dtype='float64')
    do_inside = first_pass
    tqdm.write(f'Processing {name} (inside: {do_inside})')
    row_nb = 60
    for row_i in tqdm(range(0, 6000, row_nb)):
        pix_area2 = pix_area[6000*((90-lat)//5):]
        drop_pixel(flow_dir, flow_acc, udlr_in, udlr_out, pix_area2, do_inside, row_i, row_nb)
    if first_pass:
        np.savez_compressed(f'tiles/acc/3s/{name}_acc', a=flow_acc)
    else:
        if (f'{name}_acc' in acc_dict) or (len(acc_dict) < 10):
            acc_dict[f'{name}_acc'] = flow_acc
        else:
            first = list(acc_dict.keys())[0]
            print(f'Saving {first}')
            np.save(f'tiles/acc/3s/{first}', acc_dict[first])
            del acc_dict[first]
            acc_dict[f'{name}_acc'] = flow_acc
        if f'udlr_{lat}_{lon}' in udlr_dict:
            del udlr_dict[f'udlr_{lat}_{lon}']
        if os.path.exists(f'tiles/acc/3s/{name}_acc.npz'):
            os.remove(f'tiles/acc/3s/{name}_acc.npz')
        if os.path.exists(f'tmp/udlr/udlr_{lat}_{lon}.npz'):
            os.remove(f'tmp/udlr/udlr_{lat}_{lon}.npz')
        if os.path.exists(f'tmp/udlr/udlr_{lat}_{lon}.npy'):
            os.remove(f'tmp/udlr/udlr_{lat}_{lon}.npy')
    var = [[5, 0, 1, 0, (0, 0), (1, -1), 5, -5], [-5, 0, 0, 1, (0, -1), (1, 0), 5, 5], [0, -5, 3, 2, (1, 0), (0, -1), -5, -5], [0, 5, 2, 3, (1, -1), (0, 0), -5, 5]]
    for i in range(4):
        # do the sides
        if np.any(udlr_out[i][1:-1]):
            lat2 = lat + var[i][0]
            lon2 = lon + var[i][1]
            if first_pass:
                udlr_name = f'tmp/udlr{cpu}/udlr_{lat2}_{lon2}'
            else:
                df.loc[(df.lat==lat2) & (df.lon==lon2), 'done'] = False
                udlr_name = f'tmp/udlr/udlr_{lat2}_{lon2}'
            if (not first_pass) and (os.path.basename(udlr_name) in udlr_dict):
                udlr = udlr_dict[os.path.basename(udlr_name)]
            elif os.path.exists(f'{udlr_name}.npz'):
                udlr = np.load(f'{udlr_name}.npz')['a']
            elif os.path.exists(f'{udlr_name}.npy'):
                udlr = np.load(f'{udlr_name}.npy')
            else:
                udlr = np.zeros((4, 6000), dtype='float64')
            udlr[var[i][2]] += udlr_out[var[i][3]][1:-1]
            if first_pass:
                np.savez_compressed(udlr_name, a=udlr)
            else:
                if (os.path.basename(udlr_name) in udlr_dict) or (len(udlr_dict) < 10):
                    udlr_dict[os.path.basename(udlr_name)] = udlr
                else:
                    first = list(udlr_dict.keys())[0]
                    print(f'Saving {first}')
                    np.save(f'tmp/udlr/{first}', udlr_dict[first])
                    del udlr_dict[first]
                    udlr_dict[os.path.basename(udlr_name)] = udlr
                if os.path.exists(f'{udlr_name}.npz'):
                    os.remove(f'{udlr_name}.npz')
        # do the corners
        if udlr_out[var[i][4]] != 0:
            lat2 = lat + var[i][6]
            lon2 = lon + var[i][7]
            if first_pass:
                udlr_name = f'tmp/udlr{cpu}/udlr_{lat2}_{lon2}'
            else:
                df.loc[(df.lat==lat2) & (df.lon==lon2), 'done'] = False
                udlr_name = f'tmp/udlr/udlr_{lat2}_{lon2}'
            if (not first_pass) and (os.path.basename(udlr_name) in udlr_dict):
                udlr = udlr_dict[os.path.basename(udlr_name)]
            elif os.path.exists(f'{udlr_name}.npz'):
                udlr = np.load(f'{udlr_name}.npz')['a']
            elif os.path.exists(f'{udlr_name}.npy'):
                udlr = np.load(f'{udlr_name}.npy')
            else:
                udlr = np.zeros((4, 6000), dtype='float64')
            udlr[var[i][5]] += udlr_out[var[i][4]]
            if first_pass:
                np.savez_compressed(udlr_name, a=udlr)
            else:
                if (os.path.basename(udlr_name) in udlr_dict) or (len(udlr_dict) < 10):
                    udlr_dict[os.path.basename(udlr_name)] = udlr
                else:
                    first = list(udlr_dict.keys())[0]
                    print(f'Saving {first}')
                    np.save(f'tmp/udlr/{first}', udlr_dict[first])
                    del udlr_dict[first]
                    udlr_dict[os.path.basename(udlr_name)] = udlr
                if os.path.exists(f'{udlr_name}.npz'):
                    os.remove(f'{udlr_name}.npz')
    if first_pass:
        df.to_pickle(f'tmp/df{cpu}.pkl')
    else:
        df.to_pickle('tmp/df.pkl')

@click.command()
@click.option('-n', '--numba', is_flag=True, help='Use Numba as the computing backend (otherwise, use Cython).')
@click.option('-p1', '--parallel1', default=1, help='Number of CPU cores to use for first pass.')
@click.option('-p2', '--parallel2', default=1, help='Number of CPU cores to use for second pass.')
@click.option('-r', '--reset', is_flag=1, help="Start the processing from scratch (don't download tiles if already downloaded).")
def acc_flow(numba, parallel1, parallel2, reset):
    if reset:
        shutil.rmtree('tmp', ignore_errors=True)
        shutil.rmtree('tiles/acc', ignore_errors=True)
    if numba:
        sys.path.append('numba')
        print('Using Numba')
    else:
        sys.path.append('cython')
        print('Using Cython')
    from drop_pixel import drop_pixel
    os.makedirs('tiles/dir/3s', exist_ok=True)
    os.makedirs('tiles/acc/3s', exist_ok=True)
    for cpu in range(parallel1):
        os.makedirs(f'tmp/udlr{cpu}', exist_ok=True)
    if os.path.exists('tmp/pix_area.npy'):
        pix_area = np.load('tmp/pix_area.npy')
    else:
        pix_area = pixel_area(5 / 6000)
        pix_area = pix_area / np.max(pix_area)
        np.save('tmp/pix_area', pix_area)

    if os.path.exists('tmp/df.pkl'):
        df = pd.read_pickle(f'tmp/df.pkl')
    else:
        # first pass
        print('1st pass')
        df = []
        df_ok = True
        for cpu in range(parallel1):
            if os.path.exists(f'tmp/df{cpu}.pkl'):
                df.append(pd.read_pickle(f'tmp/df{cpu}.pkl'))
            else:
                df_ok = False
        if not df_ok:
            with open('tiles.json') as f:
                dire = json.load(f)

            urls, tiles, lats, lons = [], [], [], []
            for continent in dire:
                for tile in dire[continent][1]:
                    lat = int(tile[1:3])
                    if tile[0] == 's':
                        lat = -lat
                    lon = int(tile[4:7])
                    if tile[3] == 'w':
                        lon = -lon
                    if tile not in tiles:
                        lats.append(lat + 5) # upper left
                        lons.append(lon)
                        tiles.append(tile)
                        urls.append(dire[continent][0])
            df = DataFrame({'lat': lats, 'lon': lons, 'tile': tiles, 'url': urls}).sort_values(by=['lat', 'lon'], ascending = [0, 1])  # top-down, left-right
            df.index = range(len(df))
            df_keep = df
            df = []
            i0 = 0
            size = len(df_keep) // parallel1
            i1 = size
            for cpu in range(parallel1):
                if cpu == parallel1 - 1:
                    i1 = None
                df.append(df_keep.iloc[i0:i1].copy(deep=True))
                df[cpu].to_pickle(f'tmp/df{cpu}.pkl')
                if cpu != parallel1 - 1:
                    i0 += size
                    i1 += size
        threads = []
        for cpu in range(parallel1):
            t = Thread(target=pass1, args=(cpu, drop_pixel, pix_area, df[cpu]))
            threads.append(t)
            t.start()
        for t in threads:
            t.join()
        df = pd.concat(df)
        df['done'] = False
        df.to_pickle('tmp/df.pkl')
        if parallel1 == 1:
            shutil.copytree('tmp/udlr0', 'tmp/udlr')
        else:
            os.makedirs(f'tmp/udlr', exist_ok=True)
            for cpu in range(parallel1):
                for fname in os.listdir(f'tmp/udlr{cpu}'):
                    if os.path.exists(f'tmp/udlr/{fname}'):
                        udlr = np.load(f'tmp/udlr{cpu}/{fname}')['a']
                        udlr += np.load(f'tmp/udlr/{fname}')['a']
                        np.savez_compressed(f'tmp/udlr/{fname[:-4]}', a=udlr)
                    else:
                        shutil.copyfile(f'tmp/udlr{cpu}/{fname}', f'tmp/udlr/{fname}')

    # second pass
    print('2nd pass')
    pass2(parallel2, drop_pixel, pix_area, df)

if __name__ == '__main__':
    acc_flow()


text:
#!/usr/bin/env python
# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""
Tests for pyvo.nameresolver.sesame requiring the network
"""
from __future__ import print_function, division

import os, sys, shutil, re, imp
import unittest, pdb
from urllib2 import URLError, HTTPError

from .. import sesame
import xml.etree.ElementTree as ET
from astropy.tests.helper import pytest, remote_data

mirrors = ["cfa", "cds"]
default_mirror = "cds"

if default_mirror in mirrors:
    mirrors.remove(default_mirror)
alt_mirror = mirrors[0]
sesame.set_default_endpoint(default_mirror)

do_remote = False

@remote_data
class SesameQueryTest(unittest.TestCase):

    def testExecuteStream(self):
        q = sesame.SesameQuery()
        q.names = ["ngc4258", "goob's star", "m51"]
        strm = q.execute_stream()
        doc = ET.parse(strm)
        strm.close()
        self.assertEquals("Sesame", doc.getroot().tag)

        strm = q.execute_stream(astext=True)
        doc = ET.parse(strm)
        strm.close()
        self.assertEquals("Sesame", doc.getroot().tag)

        strm = q.execute_stream(format="pc")
        lines = []
        for line in strm:
            lines.append(line)
        strm.close()
        self.assertTrue(lines[0].startswith(b'#'),
                        b"unexpected first line: " + lines[0])
        self.assertTrue(lines[5].startswith(b'%'),
                        b"unexpected mid-line: " + lines[5])

        q.dbs = "GB"
        self.assertRaises(sesame.DALQueryError, q.execute_stream)
        strm = q.execute_stream(lax=True)
        strm.close()

        del q.dbs
        self.assertEquals("", q.dbs)
        q.opts = "bo"
        self.assertRaises(sesame.DALQueryError, q.execute_stream)
        strm = q.execute_stream(lax=True)
        strm.close()

    def testExecute(self):
        q = sesame.SesameQuery()
        q.names = ["ngc4258"]

        r = q.execute()
        self.assertTrue(isinstance(r, list))
        self.assertEquals(1, len(r))

        q.names = ["ngc4258", "goob's star", "m51"]
        r = q.execute()
        self.assertTrue(isinstance(r, list))
        self.assertEquals(3, len(r))

        target = r[0]
        self.assertEquals(1, len(target.responses))
        self.assertTrue(target.resolved)
        odata = target.responses[0]
        self.assertEquals("M 106", odata.oname)
        odata = target.according_to("sim")
        self.assertTrue(odata is not None)
        self.assertEquals("M 106", odata.oname)

        target = r[1]
        self.assertEquals(0, len(target.responses))
        self.assertFalse(target.resolved)
        self.assertTrue(target.name.startswith("goob's"))
        odata = target.according_to("sim")
        self.assertTrue(odata is None)

        q.useDatabases("Simbad", "Vizier", "all")
        r = q.execute()
        self.assertTrue(isinstance(r, list))
        self.assertEquals(3, len(r))

        target = r[0]
        self.assertEquals(2, len(target.responses))
        self.assertTrue(target.resolved)
        odata = target.according_to("sim")
        self.assertTrue(odata is not None)
        self.assertEquals("M 106", odata.oname)
        odata = target.according_to("viz")
        self.assertTrue(odata is not None)
        self.assertEquals("{NGC} 4258", odata.oname)

@remote_data
class ResolveTest(unittest.TestCase):

    def testDefault(self):
        odata = sesame.resolve("NGC4258")
        self.assertTrue(odata is not None)
        self.assertFalse(isinstance(odata, list))
        self.assertTrue(odata.resolver_name.startswith("S=Simbad"))
        self.assertEquals("M 106", odata.oname)
        self.assertTrue("otype" in odata.keys())

        odata = sesame.resolve(["NGC4258", "M51"])
        self.assertTrue(odata is not None)
        self.assertTrue(isinstance(odata, list))
        self.assertTrue(odata[1].resolver_name.startswith("S=Simbad"))
        self.assertEquals("M 106", odata[0].oname)
        self.assertEquals("M  51", odata[1].oname)
        self.assertTrue("otype" in odata[0].keys())
        self.assertTrue("otype" in odata[1].keys())

    def testDb(self):
        odata = sesame.resolve("NGC4258", "NED")
        self.assertTrue(odata is not None)
        self.assertFalse(isinstance(odata, list))
        self.assertTrue(odata.resolver_name.startswith("N=NED"))
        self.assertEquals("MESSIER 106", odata.oname)
        self.assertTrue("MType" in odata.keys())

        odata = sesame.resolve(["NGC4258", "M51"], "NED")
        self.assertTrue(odata is not None)
        self.assertTrue(isinstance(odata, list))
        self.assertTrue(odata[1].resolver_name.startswith("N=NED"))
        self.assertEquals("MESSIER 106", odata[0].oname)
        self.assertEquals("MESSIER 051", odata[1].oname)
        self.assertTrue("nrefs" in odata[0].keys())
        self.assertTrue("nrefs" in odata[1].keys())

        odata = sesame.resolve("NGC4258", "Vizier")
        self.assertTrue(odata is not None)
        self.assertFalse(isinstance(odata, list))
        self.assertTrue(odata.resolver_name.startswith("V=VizieR"))
        self.assertEquals("{NGC} 4258", odata.oname)
        self.assertTrue("jpos" in odata.keys())

    def testInclude(self):
        odata = sesame.resolve("NGC4258", include=" aliases fluxes")
        self.assertTrue(odata is not None)
        self.assertFalse(isinstance(odata, list))
        self.assertTrue("alias" in odata.keys())
        self.assertTrue(len(odata.aliases) > 0)

        odata = sesame.resolve("NGC4258", include=["fluxes", "aliases"])
        self.assertTrue(odata is not None)
        self.assertFalse(isinstance(odata, list))
        self.assertTrue("alias" in odata.keys())
        self.assertTrue(len(odata.aliases) > 0)

        odata = sesame.resolve("NGC4258", include="aliases ")
        self.assertTrue(odata is not None)
        self.assertFalse(isinstance(odata, list))
        self.assertTrue("alias" in odata.keys())
        self.assertTrue(len(odata.aliases) > 0)

    def testMirror(self):
        odata = sesame.resolve("NGC4258", mirror="cds")
        self.assertTrue(odata is not None)
        self.assertTrue(odata.resolver_name.startswith("S=Simbad"))
        self.assertEquals("M 106", odata.oname)
        self.assertTrue("otype" in odata.keys())

        odata = sesame.resolve("NGC4258", mirror="cfa")
        self.assertTrue(odata is not None)
        self.assertTrue(odata.resolver_name.startswith("S=Simbad"))
        self.assertEquals("M 106", odata.oname)
        self.assertTrue("otype" in odata.keys())

        self.assertRaises(LookupError, sesame.resolve, "NGC4258", mirror="ncsa")

@remote_data
class Object2posTest(unittest.TestCase):

    def testDefault(self):
        pos = sesame.object2pos("NGC4258")
        self.assertTrue(pos is not None)
        self.assertFalse(isinstance(pos, list))
        self.assertTrue(isinstance(pos, tuple))
        self.assertEquals(2, len(pos))
        self.assertTrue(all(map(lambda p: isinstance(p,float), pos)))
        self.assertAlmostEquals(184.74008333, pos[0])
        self.assertAlmostEquals(47.30371944, pos[1])

        pos = sesame.object2pos(["NGC4258", "M51"])
        self.assertTrue(pos is not None)
        self.assertTrue(isinstance(pos, list))
        self.assertTrue(isinstance(pos[0], tuple))
        self.assertTrue(isinstance(pos[1], tuple))
        self.assertEquals(2, len(pos[1]))
        self.assertAlmostEquals(184.74008333, pos[0][0])
        self.assertAlmostEquals(47.30371944, pos[0][1])

        pos = sesame.object2pos(["NGC4258"])
        self.assertTrue(pos is not None)
        self.assertTrue(isinstance(pos, list))

    def testDb(self):
        pos = sesame.object2pos("NGC4258", "V")
        self.assertTrue(pos is not None)
        self.assertAlmostEquals(184.73, pos[0], places=2)
        self.assertAlmostEquals(47.31, pos[1], places=2)

    def testMirror(self):
        pos = sesame.object2pos("NGC4258", mirror="cds")
        self.assertTrue(pos is not None)
        self.assertFalse(isinstance(pos, list))
        self.assertTrue(isinstance(pos, tuple))
        self.assertEquals(2, len(pos))
        self.assertTrue(all(map(lambda p: isinstance(p,float), pos)))
        self.assertAlmostEquals(184.74008333, pos[0])
        self.assertAlmostEquals(47.30371944, pos[1])

        pos = sesame.object2pos("NGC4258", mirror="cfa")
        self.assertTrue(pos is not None)
        self.assertFalse(isinstance(pos, list))
        self.assertTrue(isinstance(pos, tuple))
        self.assertEquals(2, len(pos))
        self.assertTrue(all(map(lambda p: isinstance(p,float), pos)))
        self.assertAlmostEquals(184.74008333, pos[0])
        self.assertAlmostEquals(47.30371944, pos[1])

        self.assertRaises(LookupError, sesame.object2pos, "NGC4258",
                          mirror="ncsa")

@remote_data
class Object2sexaposTest(unittest.TestCase):

    def testDefault(self):
        pos = sesame.object2sexapos("NGC4258")
        self.assertTrue(pos is not None)
        self.assertFalse(isinstance(pos, list))
        self.assertTrue(isinstance(pos, str))
        self.assertEquals("12:18:57.61 +47:18:13.3", pos)

        pos = sesame.object2sexapos(["NGC4258", "M51"])
        self.assertTrue(pos is not None)
        self.assertTrue(isinstance(pos, list))
        self.assertEquals(2, len(pos))
        self.assertTrue(isinstance(pos[0], str))
        self.assertTrue(isinstance(pos[1], str))
        self.assertEquals("12:18:57.61 +47:18:13.3", pos[0])

        pos = sesame.object2sexapos(["NGC4258"])
        self.assertTrue(pos is not None)
        self.assertTrue(isinstance(pos, list))
        self.assertEquals(1, len(pos))
        self.assertTrue(isinstance(pos[0], str))
        self.assertEquals("12:18:57.61 +47:18:13.3", pos[0])

    def testDb(self):
        pos = sesame.object2sexapos("NGC4258", "V")
        self.assertTrue(pos is not None)
        self.assertEquals("12:18.9     +47:19", pos)

__all__ = "SesameQueryTest".split()
def suite():
    tests = []
    for t in __all__:
        tests.append(unittest.makeSuite(globals()[t]))
    return unittest.TestSuite(tests)

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

max_stars_repo_name: bsipocz/pyvo

max_issues_repo_name: bsipocz/pyvo

max_forks_repo_name: bsipocz/pyvo

reason: from astropy


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
#    Project: Fable Input Output
#             https://github.com/silx-kit/fabio
#
#    Copyright (C) European Synchrotron Radiation Facility, Grenoble, France
#
#    Principal author:       Jérôme Kieffer (Jerome.Kieffer@ESRF.eu)
#
#  Permission is hereby granted, free of charge, to any person
#  obtaining a copy of this software and associated documentation files
#  (the "Software"), to deal in the Software without restriction,
#  including without limitation the rights to use, copy, modify, merge,
#  publish, distribute, sublicense, and/or sell copies of the Software,
#  and to permit persons to whom the Software is furnished to do so,
#  subject to the following conditions:
#
#  The above copyright notice and this permission notice shall be
#  included in all copies or substantial portions of the Software.
#
#  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
#  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
#  OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
#  NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
#  HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
#  WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
#  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
#  OTHER DEALINGS IN THE SOFTWARE.

from __future__ import with_statement, print_function, division

__authors__ = ["Jérôme Kieffer"]
__contact__ = "Jerome.Kieffer@esrf.fr"
__license__ = "MIT"
__copyright__ = "2011-2016 ESRF"
__date__ = "29/10/2018"

import unittest
import numpy
import logging

logger = logging.getLogger(__name__)

from fabio import compression

class TestByteOffset(unittest.TestCase):
    """
    test the byte offset compression and decompression
    """
    def setUp(self):
        self.ds = numpy.array([0, 1, 2, 127, 0, 1, 2, 128, 0, 1, 2, 32767, 0, 1, 2, 32768, 0, 1, 2, 2147483647, 0, 1, 2, 2147483648, 0, 1, 2, 128, 129, 130, 32767, 32768, 128, 129, 130, 32768, 2147483647, 2147483648])
        self.ref = b'\x00\x01\x01}\x81\x01\x01~\x80\x80\xff\x01\x01\x80\xfd\x7f\x80\x01\x80\x01\x01\x80\xfe\x7f\x80\x00\x80\x00\x80\xff\xff\x01\x01\x80\x00\x80\xfd\xff\xff\x7f\x80\x00\x80\x01\x00\x00\x80\x01\x01\x80\x00\x80\xfe\xff\xff\x7f\x80\x00\x80\x00\x00\x00\x80\x00\x00\x00\x80\xff\xff\xff\xff\x01\x01~\x01\x01\x80}\x7f\x01\x80\x80\x80\x01\x01\x80~\x7f\x80\x00\x80\xff\x7f\xff\x7f\x01'

    def tearDown(self):
        unittest.TestCase.tearDown(self)
        self.ds = self.ref = None

    def testComp(self):
        """
        """
        # first with numpy
        ds = numpy.array([0, 128])
        ref = b"\x00\x80\x80\00"
        self.assertEqual(ref, compression.compByteOffset_numpy(ds), "test +128")
        ds = numpy.array([0, -128])
        ref = b'\x00\x80\x80\xff'
        self.assertEqual(ref, compression.compByteOffset_numpy(ds), "test -128")
        ds = numpy.array([10, -128])
        ref = b'\n\x80v\xff'
        self.assertEqual(ref, compression.compByteOffset_numpy(ds), "test +10 -128")
        self.assertEqual(self.ref, compression.compByteOffset_numpy(self.ds), "test larger")

        # Then with cython 32 bits
        ds = numpy.array([0, 128], dtype="int32")
        ref = b"\x00\x80\x80\00"
        self.assertEqual(ref, compression.compByteOffset_cython(ds), "test +128")
        ds = numpy.array([0, -128], dtype="int32")
        ref = b'\x00\x80\x80\xff'
        self.assertEqual(ref, compression.compByteOffset_cython(ds), "test -128")
        ds = numpy.array([10, -128], dtype="int32")
        ref = b'\n\x80v\xff'
        self.assertEqual(ref, compression.compByteOffset_cython(ds), "test +10 -128")
        self.assertEqual(self.ref, compression.compByteOffset_cython(self.ds), "test larger")

        # Then with cython 64bits
        ds = numpy.array([0, 128], dtype="int64")
        ref = b"\x00\x80\x80\00"
        self.assertEqual(ref, compression.compByteOffset_cython(ds), "test +128")
        ds = numpy.array([0, -128], dtype="int64")
        ref = b'\x00\x80\x80\xff'
        self.assertEqual(ref, compression.compByteOffset_cython(ds), "test -128")
        ds = numpy.array([10, -128], dtype="int64")
        ref = b'\n\x80v\xff'
        self.assertEqual(ref, compression.compByteOffset_cython(ds), "test +10 -128")
        self.assertEqual(self.ref, compression.compByteOffset_cython(self.ds), "test larger")

    def testSC(self):
        """test that datasets are unchanged after various compression/decompressions"""

        obt_np = compression.decByteOffset_numpy(compression.compByteOffset_numpy(self.ds))
        self.assertEqual(abs(self.ds - obt_np).max(), 0.0, "numpy-numpy algo")
        obt_cy = compression.decByteOffset_cython(compression.compByteOffset_numpy(self.ds))
        self.assertEqual(abs(self.ds - obt_cy).max(), 0.0, "cython-numpy algo")
        obt_cy2 = compression.decByteOffset_cython(compression.compByteOffset_numpy(self.ds), self.ds.size)
        self.assertEqual(abs(self.ds - obt_cy2).max(), 0.0, "cython2-numpy algo_orig")

        obt_np = compression.decByteOffset_numpy(compression.compByteOffset_cython(self.ds))
        self.assertEqual(abs(self.ds - obt_np).max(), 0.0, "numpy-numpy algo")
        obt_cy = compression.decByteOffset_cython(compression.compByteOffset_cython(self.ds))
        self.assertEqual(abs(self.ds - obt_cy).max(), 0.0, "cython-numpy algo")
        obt_cy2 = compression.decByteOffset_cython(compression.compByteOffset_cython(self.ds), self.ds.size)
        self.assertEqual(abs(self.ds - obt_cy2).max(), 0.0, "cython2-numpy algo_orig")

def suite():
    loadTests = unittest.defaultTestLoader.loadTestsFromTestCase
    testsuite = unittest.TestSuite()
    testsuite.addTest(loadTests(TestByteOffset))
    return testsuite

if __name__ == '__main__':
    runner = unittest.TextTestRunner()
    runner.run(suite())

max_stars_repo_name: picca/fabio

max_issues_repo_name: picca/fabio

max_forks_repo_name: boesecke/fabio


text:
"""Supporting functions for the 'fix' command."""
import logging

import numpy as np
import pandas as pd

from . import descriptives, params, smoothing

def do_fix(target_raw, antitarget_raw, reference,
           do_gc=True, do_edge=True, do_rmask=True, do_cluster=False):
    """Combine target and antitarget coverages and correct for biases."""
    # Load, recenter and GC-correct target & antitarget probes separately
    logging.info("Processing target: %s", target_raw.sample_id)
    cnarr, ref_matched = load_adjust_coverages(target_raw, reference,
                                               True, do_gc, do_edge, False)
    logging.info("Processing antitarget: %s", antitarget_raw.sample_id)
    anti_cnarr, ref_anti = load_adjust_coverages(antitarget_raw, reference,
                                                 False, do_gc, False, do_rmask)
    if len(anti_cnarr):
        # Combine target and antitarget bins
        cnarr.add(anti_cnarr)
        ref_matched.add(ref_anti)

    # Find reference clusters, if requested
    log2_key = 'log2'
    spread_key = 'spread'
    if do_cluster:
        ref_log2_cols = [col for col in ref_matched.data.columns
                         if col == "log2"
                         or col.startswith("log2")]
        if len(ref_log2_cols) == 1:
            logging.info("Reference does not contain any sub-clusters; "
                         "using %s", log2_key)
        else:
            # Get correlations between test sample and each reference cluster
            corr_coefs = np.array([cnarr.log2.corr(ref_matched[ref_col])
                                   for ref_col in ref_log2_cols])
            ordered = [(k, r) for r, k in sorted(zip(corr_coefs, ref_log2_cols),
                                                 reverse=True)]
            logging.info("Correlations with each cluster:\n\t%s",
                         "\n\t".join(["{}\t: {}".format(k, r)
                                      for k, r in ordered]))
            log2_key = ordered[0][0]
            if log2_key.startswith('log2_'):
                suffix = log2_key.split('_', 1)[1]
                spread_key = 'spread_' + suffix
            logging.info(" -> Choosing columns %r and %r", log2_key, spread_key)

    # Normalize coverages according to the reference
    # (Subtract the reference log2 copy number to get the log2 ratio)
    cnarr.data['log2'] -= ref_matched[log2_key]
    cnarr = apply_weights(cnarr, ref_matched, log2_key, spread_key)
    cnarr.center_all(skip_low=True)
    return cnarr

def load_adjust_coverages(cnarr, ref_cnarr, skip_low,
                          fix_gc, fix_edge, fix_rmask):
    """Load and filter probe coverages; correct using reference and GC."""
    if 'gc' in cnarr:
        # Don't choke on Picard-derived files that have the GC column
        cnarr = cnarr.keep_columns(cnarr._required_columns + ('depth',))

    # No corrections needed if there are no data rows (e.g. no antitargets)
    if not len(cnarr):
        return cnarr, ref_cnarr[:0]

    ref_matched = match_ref_to_sample(ref_cnarr, cnarr)

    # Drop bins that had poor coverage in the pooled reference
    ok_cvg_indices = ~mask_bad_bins(ref_matched)
    logging.info("Keeping %d of %d bins", sum(ok_cvg_indices), len(ref_matched))
    cnarr = cnarr[ok_cvg_indices]
    ref_matched = ref_matched[ok_cvg_indices]

    # Apply corrections for known systematic biases in coverage
    cnarr.center_all(skip_low=skip_low)
    # Skip bias corrections if most bins have no coverage (e.g. user error)
    if (cnarr['log2'] > params.NULL_LOG2_COVERAGE - params.MIN_REF_COVERAGE
        ).sum() <= len(cnarr) // 2:
        logging.warning("WARNING: most bins have no or very low coverage; "
                        "check that the right BED file was used")
    else:
        cnarr_index_reset = False
        if fix_gc:
            if 'gc' in ref_matched:
                logging.info("Correcting for GC bias...")
                cnarr = center_by_window(cnarr, .1, ref_matched['gc'])
                cnarr_index_reset = True
            else:
                logging.warning("WARNING: Skipping correction for GC bias")
        if fix_edge:
            logging.info("Correcting for density bias...")
            edge_bias = get_edge_bias(cnarr, params.INSERT_SIZE)
            cnarr = center_by_window(cnarr, .1, edge_bias)
            cnarr_index_reset = True
        if fix_rmask:
            if 'rmask' in ref_matched:
                logging.info("Correcting for RepeatMasker bias...")
                cnarr = center_by_window(cnarr, .1, ref_matched['rmask'])
                cnarr_index_reset = True
            else:
                logging.warning("WARNING: Skipping correction for "
                                "RepeatMasker bias")
        if cnarr_index_reset:
            ref_matched.data.reset_index(drop=True, inplace=True)
    return cnarr, ref_matched

def mask_bad_bins(cnarr):
    """Flag the bins with excessively low or inconsistent coverage.

    Returns
    -------
    np.array
        A boolean array where True indicates bins that failed the checks.
    """
    mask = ((cnarr['log2'] < params.MIN_REF_COVERAGE) |
            (cnarr['log2'] > -params.MIN_REF_COVERAGE) |
            (cnarr['spread'] > params.MAX_REF_SPREAD))
    if 'depth' in cnarr:
        mask |= cnarr['depth'] == 0
    if 'gc' in cnarr:
        mask |= (cnarr['gc'] > .7) | (cnarr['gc'] < .3)
    return mask

def match_ref_to_sample(ref_cnarr, samp_cnarr):
    """Filter the reference bins to match the sample (target or antitarget)."""
    # Assign each bin a unique string ID based on genomic coordinates
    samp_labeled = samp_cnarr.data.set_index(pd.Index(samp_cnarr.coords()))
    ref_labeled = ref_cnarr.data.set_index(pd.Index(ref_cnarr.coords()))
    for dset, name in ((samp_labeled, "sample"),
                       (ref_labeled, "reference")):
        dupes = dset.index.duplicated()
        if dupes.any():
            raise ValueError(("Duplicated genomic coordinates in {} set. Total duplicated regions: {}, starting with:\n"
                              "{}.").format(name, len(dset.index[dupes]), "\n".join(map(str, dset.index[dupes][:10]))))
    # Take the reference bins with IDs identical to those in the sample
    ref_matched = ref_labeled.reindex(index=samp_labeled.index)
    # Check for signs that the wrong reference was used
    num_missing = pd.isnull(ref_matched.start).sum()
    if num_missing > 0:
        raise ValueError("Reference is missing %d bins found in %s"
                         % (num_missing, samp_cnarr.sample_id))
    x = ref_cnarr.as_dataframe(ref_matched.reset_index(drop=True)
                               .set_index(samp_cnarr.data.index))
    return x

def center_by_window(cnarr, fraction, sort_key):
    """Smooth out biases according to the trait specified by sort_key.

    E.g. correct GC-biased bins by windowed averaging across similar-GC
    bins; or for similar interval sizes.
    """
    # Separate neighboring bins that could have the same key
    # (to avoid re-centering actual CNV regions -- only want an independently
    # sampled subset of presumably overall-CN-neutral bins)
    df = cnarr.data.reset_index(drop=True)
    np.random.seed(0xA5EED)
    shuffle_order = np.random.permutation(df.index)
    #df = df.reindex(shuffle_order)
    df = df.iloc[shuffle_order]
    # Apply the same shuffling to the key array as to the target probe set
    if isinstance(sort_key, pd.Series):
        # XXX why
        sort_key = sort_key.values
    sort_key = sort_key[shuffle_order]
    # Sort the data according to the specified parameter
    order = np.argsort(sort_key, kind='mergesort')
    df = df.iloc[order]
    biases = smoothing.rolling_median(df['log2'], fraction)
    # biases = smoothing.savgol(df['log2'], fraction)
    df['log2'] -= biases
    fixarr = cnarr.as_dataframe(df)
    fixarr.sort()
    return fixarr

def get_edge_bias(cnarr, margin):
    """Quantify the "edge effect" of the target tile and its neighbors.

    The result is proportional to the change in the target's coverage due to
    these edge effects, i.e. the expected loss of coverage near the target
    edges and, if there are close neighboring tiles, gain of coverage due
    to "spill over" reads from the neighbor tiles.

    (This is not the actual change in coverage. This is just a tribute.)
    """
    output_by_chrom = []
    for _chrom, subarr in cnarr.by_chromosome():
        tile_starts = subarr['start'].values
        tile_ends = subarr['end'].values
        tgt_sizes = tile_ends - tile_starts
        # Calculate coverage loss at (both edges of) each tile
        losses = edge_losses(tgt_sizes, margin)
        # Find tiled intervals within a margin (+/- bp) of the given probe
        # (excluding the probe itself), then calculate the relative coverage
        # "gain" due to the neighbors, if any
        gap_sizes = tile_starts[1:] - tile_ends[:-1]
        ok_gaps_mask = (gap_sizes < margin)
        ok_gaps = gap_sizes[ok_gaps_mask]
        left_gains = edge_gains(tgt_sizes[1:][ok_gaps_mask], ok_gaps, margin)
        right_gains = edge_gains(tgt_sizes[:-1][ok_gaps_mask], ok_gaps, margin)
        gains = np.zeros(len(subarr))
        gains[np.concatenate([[False], ok_gaps_mask])] += left_gains
        gains[np.concatenate([ok_gaps_mask, [False]])] += right_gains
        output_by_chrom.append(gains - losses)
    return pd.Series(np.concatenate(output_by_chrom), index=cnarr.data.index)

def edge_losses(target_sizes, insert_size):
    """Calculate coverage losses at the edges of baited regions.

    Letting i = insert size and t = target size, the proportional loss of
    coverage near the two edges of the baited region (combined) is:

    .. math :: i/2t

    If the "shoulders" extend outside the bait $(t < i), reduce by:

    .. math :: (i-t)^2 / 4it

    on each side, or (i-t)^2 / 2it total.
    """
    losses = insert_size / (2 * target_sizes)
    # Drop the shoulder part that would extend past the bait
    small_mask = (target_sizes < insert_size)
    t_small = target_sizes[small_mask]
    losses[small_mask] -= ((insert_size - t_small)**2
                           / (2 * insert_size * t_small))
    return losses

def edge_gains(target_sizes, gap_sizes, insert_size):
    """Calculate coverage gain from neighboring baits' flanking reads.

    Letting i = insert size, t = target size, g = gap to neighboring bait,
    the gain of coverage due to a nearby bait, if g < i, is::

    .. math :: (i-g)^2 / 4it

    If the neighbor flank extends beyond the target (t+g < i), reduce by::

    .. math :: (i-t-g)^2 / 4it

    If a neighbor overlaps the target, treat it as adjacent (gap size 0).
    """
    if not (gap_sizes <= insert_size).all():
        raise ValueError("Gaps greater than insert size:\n" +
                         gap_sizes[gap_sizes > insert_size].head())
    gap_sizes = np.maximum(0, gap_sizes)
    gains = ((insert_size - gap_sizes)**2
             / (4 * insert_size * target_sizes))
    # Drop the flank part that extends past this baited region
    past_other_side_mask = (target_sizes + gap_sizes < insert_size)
    g_past = gap_sizes[past_other_side_mask]
    t_past = target_sizes[past_other_side_mask]
    gains[past_other_side_mask] -= ((insert_size - t_past - g_past)**2
                                    / (4 * insert_size * t_past))
    return gains

def apply_weights(cnarr, ref_matched, log2_key, spread_key, epsilon=1e-4):
    """Calculate weights for each bin.

    Bin weight is an estimate of (1 - variance) and within the range
    ``(0, 1]``.

    Weights are derived from:

    - Each bin's size
    - Sample's genome-wide average (on/off-target) coverage depth
    - Sample's genome-wide observed (on/off-target) bin variances

    And with a pooled reference:

    - Each bin's coverage depth in the reference
    - The "spread" column of the reference (approx. stdev)

    These estimates of variance assume the number of aligned reads per bin
    follows a Poisson distribution, approximately log-normal.

    Parameters
    ----------
    cnarr : CopyNumArray
        Sample bins.
    ref_match : CopyNumArray
        Reference bins.
    log2_key : string
        The 'log2' column name in the reference to use. A clustered reference
        may have a suffix indicating the cluster, e.g. "log2_1".
    spread_key : string
        The 'spread' or 'spread_<cluster_id>' column name to use.
    epsilon : float
        Minimum value for bin weights, to avoid 0-weight bins causing errors
        later during segmentation. (CBS doesn't allow 0-weight bins.)

    Returns: The input `cnarr` with a `weight` column added.
    """
    # Weight by sample-level features -- works for flat reference, too
    logging.debug("Weighting bins by size and overall variance in sample")
    simple_wt = np.zeros(len(cnarr))
    # Calculate separately for on-, off-target bins
    is_anti = cnarr['gene'].isin(params.ANTITARGET_ALIASES)
    tgt_cna = cnarr[~is_anti]
    tgt_var = descriptives.biweight_midvariance(tgt_cna
                                                .drop_low_coverage()
                                                .residuals()) ** 2
    bin_sz = np.sqrt(tgt_cna['end'] - tgt_cna['start'])
    tgt_simple_wts = 1 - tgt_var / (bin_sz / bin_sz.mean())
    simple_wt[~is_anti] = tgt_simple_wts

    if is_anti.any():
        # Check for a common user error
        anti_cna = cnarr[is_anti]
        anti_ok = anti_cna.drop_low_coverage()
        frac_anti_low = 1 - (len(anti_ok) / len(anti_cna))
        if frac_anti_low > .5:
            # Off-target bins are mostly garbage -- skip reweighting
            logging.warning("WARNING: Most antitarget bins ({:.2f}%, {:d}/{:d})"
                            " have low or no coverage; is this amplicon/WGS?"
                            .format(100 * frac_anti_low,
                                    len(anti_cna) - len(anti_ok),
                                    len(anti_cna)))

        anti_var = descriptives.biweight_midvariance(anti_ok.residuals()) ** 2
        anti_bin_sz = np.sqrt(anti_cna['end'] - anti_cna['start'])
        anti_simple_wts = 1 - anti_var / (anti_bin_sz / anti_bin_sz.mean())
        simple_wt[is_anti] = anti_simple_wts

        # Report any difference in bin set variability
        var_ratio = max(tgt_var, .01) / max(anti_var, .01)
        if var_ratio > 1:
            logging.info("Targets are %.2f x more variable than antitargets",
                         var_ratio)
        else:
            logging.info("Antitargets are %.2f x more variable than targets",
                         1. / var_ratio)

    if ((ref_matched[spread_key] > epsilon).any() and
        (np.abs(np.mod(ref_matched[log2_key], 1)) > epsilon).any()):
        # Pooled/paired reference only
        logging.debug("Weighting bins by coverage spread in reference")
        # NB: spread ~= SD, so variance ~= spread^2
        fancy_wt = 1.0 - ref_matched[spread_key] ** 2
        # Average w/ simple weights, giving this more emphasis
        x = .9
        weights = (x * fancy_wt
                   + (1 - x) * simple_wt)
    else:
        # Flat reference, only 1 weight estimate
        weights = simple_wt

    return cnarr.add_columns(weight=weights.clip(epsilon, 1.0))


text:
'''
Python wrapper for the library correlation.so
'''

import ctypes as ct
import numpy as np

# C library with the bootstrap method
libc = ct.CDLL('../C/correlation.so')

# types of the arguments and return value for the function to calculate the
# correlation coefficient
libc.corr.argtypes = [np.ctypeslib.ndpointer(dtype='double'), np.ctypeslib.ndpointer(dtype='double'), ct.c_int]
libc.corr.restype = ct.c_double

# types of the arguments and return value for the function to calculate the
# standard deviation of the correlation coefficient
libc.std_corr.argtypes = [np.ctypeslib.ndpointer(dtype='double'), np.ctypeslib.ndpointer(dtype='double'), ct.c_int, ct.c_int]
libc.std_corr.restype = ct.c_double

# Python wrapper: correlation coefficient
def corr_coeff(x,y):
    '''
    computes the correlation coefficient between x and y

    return type: float

    x: numpy array of float64
    y: numpy array of float64
    '''
    length = min([len(x), len(y)])
    return libc.corr(x, y, length)

# Python wrapper: standard deviation of the correlation coefficient
def std_corr_coeff(x,y):
    '''
    computes the standard deviation for the correlation coefficient between
    x and y

    return type: float

    x: numpy array of float64
    y: numpy array of float64
    '''
    length = min([len(x), len(y)])
    return libc.std_corr(x, y, length, length)

def corr_and_error(dframe, column1, column2):
    '''
    computes the correlation coefficient between the values in columns column1
    and column2 of the dataframe dframe and its standard deviation
    The two columns are assumed to contain numeric values only.

    dframe: pandas dataframe
    column1: string
    column2: string
    '''
	# remove nan values and ensure the type is float64
    dframe = dframe.loc[-dframe[column1].isnull() & -dframe[column2].isnull(), [column1,column2]].astype('float64')
    x = np.array(dframe[column1], dtype='float64')
    y = np.array(dframe[column2], dtype='float64')
    return corr_coeff(x,y), std_corr_coeff(x,y)

# function to print the result in a self-explanatory way
def print_corr_and_error(dframe, column1, column2):
    phrase = 'The correlation coefficient between columns {} and {} is {:.2g} ± {:.2g}.'
    corr, std_corr = corr_and_error(dframe, column1, column2)
    print(phrase.format(column1, column2, corr, std_corr))


text:
# coding=utf-8
# Copyright 2021 Google LLC.
#
# 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.

# Lint as: python3
"""Tests for the Sinkhorn divergence."""

from absl.testing import absltest
from absl.testing import parameterized
import jax
import jax.numpy as jnp
import jax.test_util
from ott.geometry import geometry
from ott.geometry import pointcloud
from ott.tools import sinkhorn_divergence

class SinkhornDivergenceTest(jax.test_util.JaxTestCase):

  def setUp(self):
    super().setUp()
    self.rng = jax.random.PRNGKey(0)
    self._dim = 4
    self._num_points = 30, 37
    self.rng, *rngs = jax.random.split(self.rng, 3)
    a = jax.random.uniform(rngs[0], (self._num_points[0],))
    b = jax.random.uniform(rngs[1], (self._num_points[1],))
    self._a = a / jnp.sum(a)
    self._b = b / jnp.sum(b)

  def test_euclidean_point_cloud(self):
    rngs = jax.random.split(self.rng, 2)
    x = jax.random.uniform(rngs[0], (self._num_points[0], self._dim))
    y = jax.random.uniform(rngs[1], (self._num_points[1], self._dim))
    geometry_xx = pointcloud.PointCloud(x, x, epsilon=0.01)
    geometry_xy = pointcloud.PointCloud(x, y, epsilon=0.01)
    geometry_yy = pointcloud.PointCloud(y, y, epsilon=0.01)
    div = sinkhorn_divergence._sinkhorn_divergence(
        geometry_xy,
        geometry_xx,
        geometry_yy,
        self._a,
        self._b,
        threshold=1e-2)
    self.assertGreater(div.divergence, 0.0)
    self.assertLen(div.potentials, 3)

    # Test symmetric setting,
    # test that symmetric evaluation converges earlier/better.
    div = sinkhorn_divergence.sinkhorn_divergence(
        pointcloud.PointCloud, x, x, epsilon=1e-1,
        sinkhorn_kwargs={'inner_iterations': 1})
    self.assertAllClose(div.divergence, 0.0, rtol=1e-5, atol=1e-5)
    iters_xx = jnp.sum(div.errors[0] > 0)
    iters_xx_sym = jnp.sum(div.errors[1] > 0)
    self.assertGreater(iters_xx, iters_xx_sym)

  def test_euclidean_autoepsilon(self):
    rngs = jax.random.split(self.rng, 2)
    cloud_a = jax.random.uniform(rngs[0], (self._num_points[0], self._dim))
    cloud_b = jax.random.uniform(rngs[1], (self._num_points[1], self._dim))
    div = sinkhorn_divergence.sinkhorn_divergence(
        pointcloud.PointCloud,
        cloud_a, cloud_b,
        a=self._a, b=self._b,
        sinkhorn_kwargs=dict(threshold=1e-2))
    self.assertGreater(div.divergence, 0.0)
    self.assertLen(div.potentials, 3)
    self.assertLen(div.geoms, 3)
    self.assertAllClose(div.geoms[0].epsilon, div.geoms[1].epsilon)

  def test_euclidean_autoepsilon_not_share_epsilon(self):
    rngs = jax.random.split(self.rng, 2)
    cloud_a = jax.random.uniform(rngs[0], (self._num_points[0], self._dim))
    cloud_b = jax.random.uniform(rngs[1], (self._num_points[1], self._dim))
    div = sinkhorn_divergence.sinkhorn_divergence(
        pointcloud.PointCloud,
        cloud_a, cloud_b,
        a=self._a, b=self._b,
        sinkhorn_kwargs=dict(threshold=1e-2), share_epsilon=False)
    self.assertGreater(jnp.abs(div.geoms[0].epsilon - div.geoms[1].epsilon), 0)

  def test_euclidean_point_cloud_wrapper(self):
    rngs = jax.random.split(self.rng, 2)
    cloud_a = jax.random.uniform(rngs[0], (self._num_points[0], self._dim))
    cloud_b = jax.random.uniform(rngs[1], (self._num_points[1], self._dim))
    div = sinkhorn_divergence.sinkhorn_divergence(
        pointcloud.PointCloud,
        cloud_a, cloud_b, epsilon=0.1,
        a=self._a, b=self._b,
        sinkhorn_kwargs=dict(threshold=1e-2))
    self.assertGreater(div.divergence, 0.0)
    self.assertLen(div.potentials, 3)
    self.assertLen(div.geoms, 3)

  def test_euclidean_point_cloud_wrapper_no_weights(self):
    rngs = jax.random.split(self.rng, 2)
    cloud_a = jax.random.uniform(rngs[0], (self._num_points[0], self._dim))
    cloud_b = jax.random.uniform(rngs[1], (self._num_points[1], self._dim))
    div = sinkhorn_divergence.sinkhorn_divergence(
        pointcloud.PointCloud,
        cloud_a, cloud_b, epsilon=0.1,
        sinkhorn_kwargs=dict(threshold=1e-2))
    self.assertGreater(div.divergence, 0.0)
    self.assertLen(div.potentials, 3)
    self.assertLen(div.geoms, 3)

  def test_euclidean_point_cloud_unbalanced_wrapper(self):
    rngs = jax.random.split(self.rng, 2)
    cloud_a = jax.random.uniform(rngs[0], (self._num_points[0], self._dim))
    cloud_b = jax.random.uniform(rngs[1], (self._num_points[1], self._dim))
    div = sinkhorn_divergence.sinkhorn_divergence(
        pointcloud.PointCloud,
        cloud_a, cloud_b, epsilon=0.1,
        a=self._a +.001, b=self._b +.002,
        sinkhorn_kwargs=dict(threshold=1e-2, tau_a=0.8, tau_b=0.9))
    self.assertGreater(div.divergence, 0.0)
    self.assertLen(div.potentials, 3)
    self.assertLen(div.geoms, 3)

  def test_generic_point_cloud_wrapper(self):
    rngs = jax.random.split(self.rng, 2)
    x = jax.random.uniform(rngs[0], (self._num_points[0], self._dim))
    y = jax.random.uniform(rngs[1], (self._num_points[1], self._dim))

    # Tests with 3 cost matrices passed as args
    cxy = jnp.sum(jnp.abs(x[:, jnp.newaxis] - y[jnp.newaxis, :])**2, axis=2)
    cxx = jnp.sum(jnp.abs(x[:, jnp.newaxis] - x[jnp.newaxis, :])**2, axis=2)
    cyy = jnp.sum(jnp.abs(y[:, jnp.newaxis] - y[jnp.newaxis, :])**2, axis=2)
    div = sinkhorn_divergence.sinkhorn_divergence(
        geometry.Geometry,
        cxy, cxx, cyy, epsilon=0.1,
        a=self._a, b=self._b,
        sinkhorn_kwargs=dict(threshold=1e-2))
    self.assertIsNotNone(div.divergence)
    self.assertLen(div.potentials, 3)
    self.assertLen(div.geoms, 3)

    # Tests with 2 cost matrices passed as args
    div = sinkhorn_divergence.sinkhorn_divergence(
        geometry.Geometry,
        cxy, cxx, epsilon=0.1,
        a=self._a, b=self._b,
        sinkhorn_kwargs=dict(threshold=1e-2))
    self.assertIsNotNone(div.divergence)
    self.assertLen(div.potentials, 3)
    self.assertLen(div.geoms, 3)

    # Tests with 3 cost matrices passed as kwargs
    div = sinkhorn_divergence.sinkhorn_divergence(
        geometry.Geometry,
        cost_matrix=(cxy, cxx, cyy), epsilon=0.1,
        a=self._a, b=self._b,
        sinkhorn_kwargs=dict(threshold=1e-2))
    self.assertIsNotNone(div.divergence)
    self.assertLen(div.potentials, 3)
    self.assertLen(div.geoms, 3)

  @parameterized.parameters([True, False])
  def test_segment_sinkhorn_result(self, shuffle):
    # Test that segmented sinkhorn gives the same results:
    rngs = jax.random.split(self.rng, 4)
    x = jax.random.uniform(rngs[0], (self._num_points[0], self._dim))
    y = jax.random.uniform(rngs[1], (self._num_points[1], self._dim))
    geom_kwargs = dict(epsilon=0.01)
    sinkhorn_kwargs = dict(threshold=1e-2)
    true_divergence = sinkhorn_divergence.sinkhorn_divergence(
        pointcloud.PointCloud,
        x,
        y,
        a=self._a,
        b=self._b,
        sinkhorn_kwargs=sinkhorn_kwargs,
        **geom_kwargs).divergence

    if shuffle:
      # Now, shuffle the order of both arrays, but
      # still maintain the segment assignments:
      idx_x = jax.random.shuffle(rngs[2], jnp.arange(x.shape[0] * 2))
      idx_y = jax.random.shuffle(rngs[3], jnp.arange(y.shape[0] * 2))
    else:
      idx_x = jnp.arange(x.shape[0] * 2)
      idx_y = jnp.arange(y.shape[0] * 2)

    # Duplicate arrays:
    x_copied = jnp.concatenate((x, x))[idx_x]
    a_copied = jnp.concatenate((self._a, self._a))[idx_x]
    segment_ids_x = jnp.arange(2).repeat(x.shape[0])[idx_x]

    y_copied = jnp.concatenate((y, y))[idx_y]
    b_copied = jnp.concatenate((self._b, self._b))[idx_y]
    segment_ids_y = jnp.arange(2).repeat(y.shape[0])[idx_y]

    segmented_divergences = sinkhorn_divergence.segment_sinkhorn_divergence(
        x_copied,
        y_copied,
        segment_ids_x=segment_ids_x,
        segment_ids_y=segment_ids_y,
        indices_are_sorted=False,
        weights_x=a_copied,
        weights_y=b_copied,
        sinkhorn_kwargs=sinkhorn_kwargs,
        **geom_kwargs)

    self.assertArraysAllClose(
        true_divergence.repeat(2), segmented_divergences)

  def test_segment_sinkhorn_different_segment_sizes(self):

    # Test other array sizes
    x1 = jnp.arange(10)[:, None].repeat(2, axis=1)
    y1 = jnp.arange(11)[:, None].repeat(2, axis=1) + 0.1

    # Should have larger divergence since further apart:
    x2 = jnp.arange(12)[:, None].repeat(2, axis=1)
    y2 = 2 * jnp.arange(13)[:, None].repeat(2, axis=1) + 0.1

    segmented_divergences = sinkhorn_divergence.segment_sinkhorn_divergence(
        jnp.concatenate((x1, x2)),
        jnp.concatenate((y1, y2)),
        num_per_segment_x=jnp.array([10, 12]),
        num_per_segment_y=jnp.array([11, 13]),
        epsilon=0.01)

    self.assertEqual(segmented_divergences.shape[0], 2)
    self.assertGreater(segmented_divergences[1], segmented_divergences[0])

    true_divergences = jnp.array([
        sinkhorn_divergence.sinkhorn_divergence(
            pointcloud.PointCloud, x, y, epsilon=0.01).divergence
        for x, y in zip((x1, x2), (y1, y2))
    ])
    self.assertArraysAllClose(segmented_divergences, true_divergences)

  @parameterized.parameters(
      [dict(anderson_acceleration=3), 1e-2],
      [dict(anderson_acceleration=6), None],
      [dict(chg_momentum_from=20), 1e-3],
      [dict(chg_momentum_from=30), None],
      [dict(momentum=1.05), 1e-3],
      [dict(momentum=1.01), None])
  def test_euclidean_momentum_params(self, sinkhorn_kwargs, epsilon):
    # check if sinkhorn divergence sinkhorn_kwargs parameters used for
    # momentum/Anderson are properly overriden for the symmetric (x,x) and
    # (y,y) parts.
    rngs = jax.random.split(self.rng, 2)
    threshold = 3.2e-3
    cloud_a = jax.random.uniform(rngs[0], (self._num_points[0], self._dim))
    cloud_b = jax.random.uniform(rngs[1], (self._num_points[1], self._dim))
    div = sinkhorn_divergence.sinkhorn_divergence(
        pointcloud.PointCloud,
        cloud_a,
        cloud_b,
        epsilon=epsilon,
        a=self._a,
        b=self._b,
        sinkhorn_kwargs=sinkhorn_kwargs.update({'threshold': threshold}))
    self.assertGreater(threshold, div.errors[0][-1])
    self.assertGreater(threshold, div.errors[1][-1])
    self.assertGreater(threshold, div.errors[2][-1])
    self.assertGreater(div.divergence, 0.0)

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

max_stars_repo_name: MUCDK/ott

max_issues_repo_name: RamonYeung/ott

max_forks_repo_name: RamonYeung/ott

reason: import jax


text:
# Simple Linear Regression

# Importing the libraries
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

# Importing the dataset
dataset = pd.read_csv('Salary_Data.csv')
X = dataset.iloc[:, :-1].values
Y = dataset.iloc[:, -1].values

# Splitting the dataset into the training set and test set
from sklearn.cross_validation import train_test_split
X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size= 1.0/3, random_state= 0)

"""# Feature scaling
from sklearn.preprocessing import StandardScaler
sc_X = StandardScaler()
X_train = sc_X.fit_transform(X_train)
X_test = sc_X.transform(X_test)"""

# Fitting simple linear regression to the training set
from sklearn.linear_model import LinearRegression
regressor = LinearRegression()
regressor.fit(X_train,Y_train)

# Predicting the test set results
Y_pred = regressor.predict(X_test)

# Visualising the training set results
plt.scatter(X_train,Y_train, color = 'red')
plt.plot(X_train, regressor.predict(X_train), color = 'blue')
plt.title('Salary vs Experience (Training set)')
plt.xlabel('Years of experience')
plt.ylabel('Salary ($)')

# Visualising the test set results
plt.figure()
plt.scatter(X_test,Y_test, color = 'red')
plt.plot(X_train, regressor.predict(X_train), color = 'blue')
plt.title('Salary vs Experience (Test set)')
plt.xlabel('Years of experience')
plt.ylabel('Salary ($)')
plt.show()

max_stars_repo_name: BigMani/Machine-learning-projects

max_issues_repo_name: BigMani/Machine-learning-projects

max_forks_repo_name: BigMani/Machine-learning-projects


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
############
# Standard #
############
import time
from collections import OrderedDict
###############
# Third Party #
###############
import numpy as np

##########
# Module #
##########
from pswalker.examples import (patch_pims, _calc_cent_x, _m1_calc_cent_x,
                               _m1_m2_calc_cent_x)
from pswalker.sim import (source, mirror, pim)

def test_patch_pims_with_no_bounce_func(simple_two_bounce_system):
    s, m1, m2 = simple_two_bounce_system
    pim_1 = [pim.PIM('test_pim', name='pim_1', z=3)]
    pim_2 = [pim.PIM('test_pim', name='pim_2', z=5)]
    pim_3 = [pim.PIM('test_pim', name='pim_3', z=7)]

    pim_1 = patch_pims(pim_1)
    assert pim_1.detector.stats2.centroid.x.get() == _calc_cent_x(s, pim_1)

    pim_2 = patch_pims(pim_2, m1)
    assert pim_2.detector.stats2.centroid.x.get() == _calc_cent_x(s, pim_2)

    pim_3 = patch_pims(pim_3, [m1,m2])
    assert pim_3.detector.stats2.centroid.x.get() == _calc_cent_x(s, pim_3)

def test_patch_pims_with_one_bounce_func(simple_two_bounce_system):
    s, m1, m2 = simple_two_bounce_system
    pim_1 = [pim.PIM('test_pim', name='pim_1', z=13)]
    pim_2 = [pim.PIM('test_pim', name='pim_1', z=15)]

    pim_1 = patch_pims(pim_1, [m1,m2])
    assert pim_1.detector.stats2.centroid.x.get() == _m1_calc_cent_x(s, m1,
                                                                     pim_1)

    pim_2 = patch_pims(pim_2, m1)
    assert pim_2.detector.stats2.centroid.x.get() == _m1_calc_cent_x(s, m1,
                                                                     pim_2)

def test_patch_pims_with_two_bounce_func(simple_two_bounce_system):
    s, m1, m2 = simple_two_bounce_system
    pim_1 = [pim.PIM('test_pim', name='pim_1', z=25)]
    pim_1 = patch_pims(pim_1, [m1,m2])
    assert pim_1.detector.stats2.centroid.x.get() == _m1_m2_calc_cent_x(
        s, m1, m2, pim_1)

def test_pim_set_and_read(one_bounce_system):
    s, mot, pim = one_bounce_system
    set_x = 100
    set_z = 10
    pim.sim_x.put(set_x)
    pim.sim_z.put(set_z)
    pim = patch_pims(pim, [mot])
    assert pim.sim_x.get() == set_x
    assert pim.sim_z.get() == set_z
    assert pim.detector.stats2.centroid.x.get() == _m1_calc_cent_x(s, mot, pim)

max_stars_repo_name: ZLLentz/pswalker

max_issues_repo_name: ZLLentz/pswalker

max_forks_repo_name: ZLLentz/pswalker