text:
import tkinter
from sarpy_gui_apps.apps.wake_tool.panels.side_panel import SidePanel
from tkinter_gui_builder.panel_templates.image_canvas.image_canvas import ImageCanvas
import tkinter.colorchooser as colorchooser
from sarpy_gui_apps.supporting_classes.sarpy_canvas_image import SarpyCanvasDisplayImage
from tkinter_gui_builder.panel_templates.widget_panel.widget_panel import AbstractWidgetPanel
import sarpy.geometry.point_projection as point_projection
import sarpy.geometry.geocoords as geocoords
import numpy as np
import math

class AppVariables:
    def __init__(self):
        self.image_fname = "None"       # type: str
        self.sicd_metadata = None
        self.current_tool_selection = None      # type: str

        self.arrow_id = None             # type: int
        self.point_id = None            # type: int
        self.horizontal_line_id = None      # type: int

        self.line_color = "red"
        self.line_width = 3
        self.horizontal_line_width = 2
        self.horizontal_line_color = "green"

class WakeTool(AbstractWidgetPanel):
    side_panel = SidePanel          # type: SidePanel
    image_canvas = ImageCanvas      # type: ImageCanvas

    def __init__(self, master):
        # set the master frame
        master_frame = tkinter.Frame(master)
        AbstractWidgetPanel.__init__(self, master_frame)
        widget_list = ["image_canvas", "side_panel"]
        self.init_w_vertical_layout(widget_list)
        self.variables = AppVariables()

        self.side_panel.set_spacing_between_buttons(0)
        self.image_canvas.variables.canvas_image_object = SarpyCanvasDisplayImage()     # type: SarpyCanvasDisplayImage
        self.image_canvas.set_canvas_size(600, 400)

        # need to pack both master frame and self, since this is the main app window.
        master_frame.pack()
        self.pack()

        # set up event listeners
        self.side_panel.buttons.line_draw.on_left_mouse_click(self.callback_press_line_button)
        self.side_panel.buttons.point_draw.on_left_mouse_click(self.callback_press_point_button)
        self.side_panel.buttons.foreground_color.on_left_mouse_click(self.callback_select_color)
        self.side_panel.buttons.zoom_in.on_left_mouse_click(self.callback_set_to_zoom_in)
        self.side_panel.buttons.zoom_out.on_left_mouse_click(self.callback_set_to_zoom_out)

        self.image_canvas.variables.line_width = self.variables.line_width
        self.image_canvas.canvas.on_left_mouse_click(self.callback_handle_left_mouse_click)
        self.image_canvas.canvas.on_left_mouse_motion(self.callback_on_left_mouse_motion)

        self.side_panel.file_selector.set_fname_filters(["*.NITF", ".nitf"])
        self.side_panel.file_selector.select_file.on_left_mouse_click(self.callback_select_file)

    def callback_select_file(self, event):
        self.side_panel.file_selector.event_select_file(event)
        if self.side_panel.file_selector.fname:
            self.variables.image_fname = self.side_panel.file_selector.fname
        self.image_canvas.init_with_fname(self.variables.image_fname)

    def callback_press_line_button(self, event):
        self.side_panel.buttons.set_active_button(self.side_panel.buttons.line_draw)
        self.image_canvas.set_current_tool_to_draw_arrow_by_dragging()
        self.image_canvas.variables.current_shape_id = self.variables.arrow_id

    def callback_press_point_button(self, event):
        self.side_panel.buttons.set_active_button(self.side_panel.buttons.point_draw)
        self.image_canvas.set_current_tool_to_draw_point()
        self.image_canvas.variables.current_shape_id = self.variables.point_id

    def callback_set_to_zoom_in(self, event):
        self.side_panel.buttons.set_active_button(self.side_panel.buttons.zoom_in)
        self.image_canvas.set_current_tool_to_zoom_in()

    def callback_set_to_zoom_out(self, event):
        self.side_panel.buttons.set_active_button(self.side_panel.buttons.zoom_out)
        self.image_canvas.set_current_tool_to_zoom_out()

    def callback_select_color(self, event):
        self.side_panel.buttons.set_active_button(self.side_panel.buttons.foreground_color)
        color = colorchooser.askcolor()[1]
        self.image_canvas.change_shape_color(self.image_canvas.variables.current_shape_id, color)

    def callback_handle_left_mouse_click(self, event):
        # first do all the normal mouse click functionality of the canvas
        self.image_canvas.callback_handle_left_mouse_click(event)
        # now set the object ID's accordingly, we do this so we don't draw multiple arrows or points
        if self.image_canvas.variables.current_tool == self.image_canvas.TOOLS.DRAW_ARROW_BY_DRAGGING:
            self.variables.arrow_id = self.image_canvas.variables.current_shape_id
        if self.image_canvas.variables.current_tool == self.image_canvas.TOOLS.DRAW_POINT_BY_CLICKING:
            self.variables.point_id = self.image_canvas.variables.current_shape_id
        if self.variables.point_id is not None and self.variables.arrow_id is not None:
            self.update_distance()

    def callback_on_left_mouse_motion(self, event):
        self.image_canvas.callback_handle_left_mouse_motion(event)
        self.update_distance()

    def update_distance(self):
        if self.variables.point_id is not None and self.variables.arrow_id is not None:
            # calculate horizontal line segment
            point_x, point_y = self.image_canvas.get_shape_canvas_coords(self.variables.point_id)
            line_slope, line_intercept = self.get_line_slope_and_intercept()
            end_x = (point_y - line_intercept) / line_slope
            end_y = point_y
            horizontal_line_coords = (point_x, point_y, end_x, end_y)
            # save last object selected on canvas (either the point or the line)
            last_shape_id = self.image_canvas.variables.current_shape_id
            if self.variables.horizontal_line_id is None:
                self.variables.horizontal_line_id = \
                    self.image_canvas.create_new_line(horizontal_line_coords,
                                                      fill=self.variables.horizontal_line_color,
                                                      width=self.variables.horizontal_line_width)
            else:
                self.image_canvas.modify_existing_shape_using_canvas_coords(self.variables.horizontal_line_id, horizontal_line_coords)
            # set current object ID back to what it was after drawing the horizontal line
            self.image_canvas.variables.current_shape_id = last_shape_id
            canvas_distance = self.image_canvas.get_canvas_line_length(self.variables.horizontal_line_id)
            pixel_distance = self.image_canvas.get_image_line_length(self.variables.horizontal_line_id)
            geo_distance = self.calculate_wake_distance()
            self.side_panel.info_panel.canvas_distance_val.set_text("{:.2f}".format(canvas_distance))
            self.side_panel.info_panel.pixel_distance_val.set_text("{:.2f}".format(pixel_distance))
            self.side_panel.info_panel.geo_distance_val.set_text("{:.2f}".format(geo_distance))

    def calculate_wake_distance(self):
        horizontal_line_image_coords = self.image_canvas.canvas_shape_coords_to_image_coords(self.variables.horizontal_line_id)
        sicd_meta = self.image_canvas.variables.canvas_image_object.reader_object.sicd_meta
        points = np.asarray(np.reshape(horizontal_line_image_coords, (2, 2)))
        ecf_ground_points = point_projection.image_to_ground(points, sicd_meta)
        geo_ground_point_1 = geocoords.ecf_to_geodetic((ecf_ground_points[0, 0], ecf_ground_points[0, 1], ecf_ground_points[0, 2]))
        geo_ground_point_2 = geocoords.ecf_to_geodetic((ecf_ground_points[1, 0], ecf_ground_points[1, 1], ecf_ground_points[1, 2]))
        distance = math.sqrt( (ecf_ground_points[0, 0] - ecf_ground_points[1, 0])**2 +
                              (ecf_ground_points[0, 1] - ecf_ground_points[1, 1])**2 +
                              (ecf_ground_points[0, 2] - ecf_ground_points[1, 2])**2)
        return distance

    def get_line_slope_and_intercept(self):
        line_coords = self.image_canvas.canvas.coords(self.variables.arrow_id)

        line_x1, line_x2 = line_coords[0], line_coords[2]
        line_y1, line_y2 = line_coords[1], line_coords[3]

        line_slope = (line_y2 - line_y1)/(line_x2 - line_x1)
        line_intercept = line_y1 - line_slope * line_x1
        return line_slope, line_intercept

if __name__ == '__main__':
    root = tkinter.Tk()
    app = WakeTool(root)
    root.mainloop()

max_stars_repo_name: spowlas/sarpy

max_issues_repo_name: spowlas/sarpy

max_forks_repo_name: spowlas/sarpy


text:
# -*- coding: utf-8 -*-
"""
Created on Sun Jan  3 12:46:32 2021

@author: mclea
"""

import numpy as np

class field:
    """
    """

    def __init__(self, k):
        self.n = 2**k+2
        self.x = np.zeros((self.n, self.n))

    def set_BC(self):
        # bottom
        self.x[0, 0:(int(self.n/2))] = np.linspace(13, 5, int((self.n/2)))
        self.x[0, int(self.n/2):int(3*self.n/4)] = 5
        self.x[0, int(3*self.n/4):-1] = np.linspace(5, 13, int(self.n/4))

        # top
        self.x[-1, :] = 21

        # left
        self.x[0:int(3*self.n/8), 0] = np.linspace(13, 40, int(3*self.n/8))
        self.x[int(self.n/2):, 0] = np.linspace(40, 21, int(self.n/2))

        # right
        self.x[0:int(self.n/2), -1] = np.linspace(13, 40, int(self.n/2))
        self.x[int(5*self.n/8):-1, -1] = np.linspace(40, 21, int(3*self.n/8))

        # heaters
        self.x[int(3*self.n/8):int(self.n/2), 0:int(self.n/8+1)] = 40
        self.x[int(self.n/2):int(5*self.n/8), int(7*self.n/8):(self.n+1)] = 40

f = field(9)
f.set_BC()

ext: py

max_stars_repo_path: Fluid Dynamics/CFD/Finite Volume Method/field.py

max_stars_repo_name: mcleantom/CFDPython

max_issues_repo_name: mcleantom/CFDPython

max_forks_repo_name: mcleantom/CFDPython

reason: import numpy


text:
# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
#
#   Copyright (c) 2008 Emanuele Olivetti <emanuele@relativita.com>
#   See COPYING file distributed along with the PyMVPA package for the
#   copyright and license terms.
#
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
"""Bayesian Linear Regression (BLR)."""

__docformat__ = 'restructuredtext'

import numpy as np

from mvpa2.base.state import ConditionalAttribute
from mvpa2.clfs.base import Classifier, accepts_dataset_as_samples

if __debug__:
    from mvpa2.misc import debug

class BLR(Classifier):
    """Bayesian Linear Regression (BLR).

    """

    predicted_variances = ConditionalAttribute(enabled=False,
        doc="Variance per each predicted value")

    log_marginal_likelihood = ConditionalAttribute(enabled=False,
        doc="Log Marginal Likelihood")

    __tags__ = [ 'blr', 'regression', 'linear' ]

    def __init__(self, sigma_p = None, sigma_noise=1.0, **kwargs):
        """Initialize a BLR regression analysis.

        Parameters
        ----------
        sigma_noise : float
          the standard deviation of the gaussian noise.
          (Defaults to 0.1)

        """
        # init base class first
        Classifier.__init__(self, **kwargs)

        # pylint happiness
        self.w = None

        # It does not make sense to calculate a confusion matrix for a
        # BLR:
        self.ca.enable('training_stats', False)

        # set the prior on w: N(0,sigma_p) , specifying the covariance
        # sigma_p on w:
        self.sigma_p = sigma_p

        # set noise level:
        self.sigma_noise = sigma_noise

        self.ca.predicted_variances = None
        self.ca.log_marginal_likelihood = None
        # Yarik: what was those about??? just for future in
        #        compute_log_marginal_likelihood ?
        # self.targets = None
        pass

    def __repr__(self):
        """String summary of the object
        """
        return """BLR(w=%s, sigma_p=%s, sigma_noise=%f, enable_ca=%s)""" % \
               (self.w, self.sigma_p, self.sigma_noise, str(self.ca.enabled))

    def compute_log_marginal_likelihood(self):
        """
        Compute log marginal likelihood using self.train_fv and self.targets.
        """
        # log_marginal_likelihood = None
        # return log_marginal_likelihood
        raise NotImplementedError

    def _train(self, data):
        """Train regression using `data` (`Dataset`).
        """
        # BLR relies on numerical labels
        train_labels = self._attrmap.to_numeric(data.sa[self.get_space()].value)
        # provide a basic (i.e. identity matrix) and correct prior
        # sigma_p, if not provided before or not compliant to 'data':
        if self.sigma_p == None: # case: not provided
            self.sigma_p = np.eye(data.samples.shape[1]+1)
        elif self.sigma_p.shape[1] != (data.samples.shape[1]+1): # case: wrong dimensions
            self.sigma_p = np.eye(data.samples.shape[1]+1)
        else:
            # ...then everything is OK :)
            pass

        # add one fake column of '1.0' to model the intercept:
        self.samples_train = np.hstack([data.samples,np.ones((data.samples.shape[0],1))])
        if type(self.sigma_p)!=type(self.samples_train): # if sigma_p is a number...
            self.sigma_p = np.eye(self.samples_train.shape[1])*self.sigma_p # convert in matrix
            pass

        self.A_inv = np.linalg.inv(1.0/(self.sigma_noise**2) *
                                  np.dot(self.samples_train.T,
                                        self.samples_train) +
                                  np.linalg.inv(self.sigma_p))
        self.w = 1.0/(self.sigma_noise**2) * np.dot(self.A_inv,
                                                   np.dot(self.samples_train.T,
                                                         train_labels))
        pass

    @accepts_dataset_as_samples
    def _predict(self, data):
        """
        Predict the output for the provided data.
        """

        data = np.hstack([data,np.ones((data.shape[0],1),dtype=data.dtype)])
        predictions = np.dot(data,self.w)

        if self.ca.is_enabled('predicted_variances'):
            # do computation only if conditional attribute was enabled
            self.ca.predicted_variances = np.dot(data, np.dot(self.A_inv, data.T)).diagonal()[:,np.newaxis]
        self.ca.estimates = predictions
        return predictions

    def set_hyperparameters(self,*args):
        """
        Set hyperparameters' values.

        Note that this is a list so the order of the values is
        important.
        """
        args=args[0]
        self.sigma_noise = args[0]
        if len(args)>1:
            self.sigma_p = np.array(args[1:]) # XXX check if this is ok
            pass
        return

    pass

if __name__ == "__main__":
    import pylab
    pylab.close("all")
    pylab.ion()

    from mvpa2.misc.data_generators import linear_awgn

    train_size = 10
    test_size = 100
    F = 1 # dimensions of the dataset

    # np.random.seed(1)

    slope = np.random.rand(F)
    intercept = np.random.rand(1)
    print "True slope:",slope
    print "True intercept:",intercept

    dataset_train = linear_awgn(train_size, intercept=intercept, slope=slope)

    dataset_test = linear_awgn(test_size, intercept=intercept, slope=slope, flat=True)

    regression = True
    logml = False

    b = BLR(sigma_p=np.eye(F+1), sigma_noise=0.1)
    b.ca.enable("predicted_variances")
    b.train(dataset_train)
    predictions = b.predict(dataset_test.samples)
    print "Predicted slope and intercept:",b.w

    if F==1:
        pylab.plot(dataset_train.samples,
                   dataset_train.sa[b.get_space()].value,
                   "ro", label="train")

        pylab.plot(dataset_test.samples, predictions, "b-", label="prediction")
        pylab.plot(dataset_test.samples,
                   predictions+np.sqrt(b.ca.predicted_variances),
                   "b--", label="pred(+/-)std")
        pylab.plot(dataset_test.samples,
                   predictions-np.sqrt(b.ca.predicted_variances),
                   "b--", label=None)
        pylab.legend()
        pylab.xlabel("samples")
        pylab.ylabel("labels")
        pylab.title("Bayesian Linear Regression on dataset 'linear_AWGN'")
        pass

max_issues_repo_name: mih/PyMVPA

max_forks_repo_name: mih/PyMVPA


text:
##################################################################################
#
#   pysvo.extcurves: rudimentary extinction curve module
#
#   Main Functions:
#
#       - Cardelli1989:    Cardelli et al. (1989) extinction curve
#       - Schlafly2016:    Schlafly et al. (2016) extinction curve
#
#   For more options, use the extinction module available on pip.
#
###################################################################################
import numpy as np
from pysvo import fitting_tools

"""
Defining the transformation constants between a set of filters
relative to to the V-band extinction.
"""

def Cardelli1989( wavelength, RV=3.1 ):
    """
    Get A(lambda) / AV using recipe of Cardelli et al (1989).

    Input:
        wavelength - in Angstroms
    Optional:
        RV         - the preferred Rv value
    Output:
        Extinction relative to V-band: A(lambda) / AV

    History:
        2016       - Written             - F. Anders (AIP)
        2019-06-10 - Renamed for starsed - F. Anders (ICCUB)
    """
    x = 10000. / wavelength     # x is in (mum)^-1
    # Different coefficients for UV, optical, and infrared waveleghts:
    # opt: 0.25 - 0.9 micron, IR: 0.9 - 3.3 micron
    a = np.piecewise( x, [x<1.1, x>=1.1],
                      [lambda x: 0.574 * x**1.61,
                       lambda x: 1 + 0.17699*(x-1.82)  - 0.50447*(x-1.82)**2
                                 - 0.02427*(x-1.82)**3 + 0.72085*(x-1.82)**4
                                 + 0.01979*(x-1.82)**5 - 0.77530*(x-1.82)**6
                                 + 0.32999*(x-1.82)**7  ] )
    b = np.piecewise( x, [x<1.1, x>=1.1],
                      [lambda x: -0.527 * x**1.61,
                       lambda x:   1.41338*(x-1.82)    + 2.28305*(x-1.82)**2
                                 + 1.07233*(x-1.82)**3 - 5.38434*(x-1.82)**4
                                 - 0.62251*(x-1.82)**5 + 5.30260*(x-1.82)**6
                                 - 2.09002*(x-1.82)**7  ] )
    return a + b / RV

def Schlafly2016( wavelength, RV=3.32 ):
    """
    Based on E. Schlafly's code: http://e.schlaf.ly/apored/extcurve.html

    Returns the extinction curve, A(lambda)/A(5420 A), according to
    Schlafly+2016. RV is transformed to the parameter "x," which controls the overall shape of
    the extinction curve in an R(V)-like way. The extinction curve is based on broad band photometry between the PS1 g
    band and the WISE W2 band, which have effective wavelengths between 5000
    and 45000 A, and is blindly extrapolated outside that
    range.  The gray component of the extinction curve is fixed by enforcing
    A(H)/A(K) = 1.55 (Indebetouw+2005).

    Args:
        wavelength - in Angstrom
    Optional:
        RV         - default is 3.32 (Schlafly+2016, stdev 0.18)

    Returns:
        A(lambda)/A(5420 A)

    History:
        2017       - adapted for starhorse - F. Anders (AIP)
        2019-06-10 - ported to starsed     - F. Anders (ICCUB)
    """

    # Transform RV' to Schlafly's x parameter
    x = (RV - 3.3) / 9.1
    #   ra: extinction vector at anchor wavelengths
    #   dra: derivative of extinction vector at anchor wavelengths
    #   lam: anchor wavelengths (angstroms)
    ra = np.array([ 0.65373283,  0.39063843,  0.20197893,  0.07871701, -0.00476316,
                   -0.14213929, -0.23660605, -0.28522577, -0.321301  , -0.33503192])
    dra = np.array([-0.54278669,  0.03404903,  0.36841725,  0.42265873,  0.38247769,
                     0.14148814, -0.04020524, -0.13457319, -0.26883343, -0.36269229])
    lam = np.array([  5032.36441067,   6280.53335141,   7571.85928312,   8690.89321059,
                          9635.52560909,  12377.04268274,  16381.78146718,  21510.20523237,
                         32949.54009328,  44809.4919175 ])
    # Indebetouw (2005) A_H / A_K
    rhk = 1.55

    anchors = ra + x*dra
    # fix gray component so that A(H)/A(K) = 1.55
    anchors += (-anchors[6] + rhk*anchors[7])/(1 - rhk)
    cs0 = fitting_tools.CubicSpline(lam, anchors, yp='3d=0')
    # normalize at 5420 angstroms
    ec  = fitting_tools.CubicSpline(lam, anchors/cs0(5420.), yp='3d=0')
    return ec( wavelength )

max_stars_repo_name: fjaellet/pysvo

max_issues_repo_name: fjaellet/pysvo

max_forks_repo_name: fjaellet/pysvo


text:
# -*- coding: utf-8 -*-
"""
Created on Thu Apr  2 13:19:58 2015
@author: ajaver
"""
import os

import cv2
import tables
import numpy as np

from tierpsy.analysis.compress.BackgroundSubtractor import BackgroundSubtractorVideo
from tierpsy.analysis.compress.extractMetaData import store_meta_data, read_and_save_timestamp
from tierpsy.analysis.compress.selectVideoReader import selectVideoReader
from tierpsy.helper.params import compress_defaults, set_unit_conversions
from tierpsy.helper.misc import TimeCounter, print_flush, TABLE_FILTERS
from tierpsy.analysis.split_fov.helper import parse_camera_serial
from tierpsy.analysis.split_fov.FOVMultiWellsSplitter import FOVMultiWellsSplitter

def getROIMask(
        image,
        min_area,
        max_area,
        thresh_block_size,
        thresh_C,
        dilation_size,
        keep_border_data,
        is_light_background,
        wells_mask=None):
    '''
    Calculate a binary mask to mark areas where it is possible to find worms.
    Objects with less than min_area or more than max_area pixels are rejected.
        > min_area -- minimum blob area to be considered in the mask
        > max_area -- max blob area to be considered in the mask
        > thresh_C -- threshold used by openCV adaptiveThreshold
        > thresh_block_size -- block size used by openCV adaptiveThreshold
        > dilation_size -- size of the structure element to dilate the mask
        > keep_border_data -- (bool) if false it will reject any blob that touches the image border
        > is_light_background -- (bool) true if bright field, false if fluorescence
        > wells_mask -- (bool 2D) mask that covers (with False) the edges of wells in a MW plate
    '''
    # Objects that touch the limit of the image are removed. I use -2 because
    # openCV findCountours remove the border pixels
    IM_LIMX = image.shape[0] - 2
    IM_LIMY = image.shape[1] - 2

    #this value must be at least 3 in order to work with the blocks
    thresh_block_size = max(3, thresh_block_size)
    if thresh_block_size % 2 == 0:
        thresh_block_size += 1  # this value must be odd

    #let's add a median filter, this will smooth the image, and eliminate small variations in intensity
    # now done with opencv instead of scipy
    image = cv2.medianBlur(image, 5)

    # adaptative threshold is the best way to find possible worms. The
    # parameters are set manually, they seem to work fine if there is no
    # condensation in the sample
    if not is_light_background:  # invert the threshold (change thresh_C->-thresh_C and cv2.THRESH_BINARY_INV->cv2.THRESH_BINARY) if we are dealing with a fluorescence image
        mask = cv2.adaptiveThreshold(
            image,
            255,
            cv2.ADAPTIVE_THRESH_MEAN_C,
            cv2.THRESH_BINARY,
            thresh_block_size,
            -thresh_C)
    else:
        mask = cv2.adaptiveThreshold(
            image,
            255,
            cv2.ADAPTIVE_THRESH_MEAN_C,
            cv2.THRESH_BINARY_INV,
            thresh_block_size,
            thresh_C)

    # find the contour of the connected objects (much faster than labeled
    # images)

    contours, hierarchy = cv2.findContours(
        mask.copy(), cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)[-2:]

    # find good contours: between max_area and min_area, and do not touch the
    # image border
    goodIndex = []
    for ii, contour in enumerate(contours):
        if not keep_border_data:
            if wells_mask is None:
                # eliminate blobs that touch a border
                # TODO: double check this next line. I suspect contour is in
                # x,y and not row columns
                keep = not np.any(contour == 1) and \
                    not np.any(contour[:, :, 0] ==  IM_LIMY)\
                    and not np.any(contour[:, :, 1] == IM_LIMX)
            else:
                # keep if no pixel of contour is in the 0 part of the mask
                keep = not np.any(wells_mask[contour[:, :, 1],
                                             contour[:, :, 0]] == 0)
        else:
            keep = True

        if keep:
            area = cv2.contourArea(contour)
            if (area >= min_area) and (area <= max_area):
                goodIndex.append(ii)

    # typically there are more bad contours therefore it is cheaper to draw
    # only the valid contours
    mask = np.zeros(image.shape, dtype=image.dtype)
    for ii in goodIndex:
        cv2.drawContours(mask, contours, ii, 1, cv2.FILLED)

    # drawContours left an extra line if the blob touches the border. It is
    # necessary to remove it
    mask[0, :] = 0
    mask[:, 0] = 0
    mask[-1, :] = 0
    mask[:, -1] = 0

    # dilate the elements to increase the ROI, in case we are missing
    # something important
    struct_element = cv2.getStructuringElement(
        cv2.MORPH_ELLIPSE, (dilation_size, dilation_size))
    mask = cv2.dilate(mask, struct_element, iterations=3)

    return mask

def normalizeImage(img):
    # normalise image intensities if the data type is other
    # than uint8
    image = image.astype(np.double)

    imax = img.max()
    imin = img.min()
    factor = 255/(imax-imin)

    imgN = ne.evaluate('(img-imin)*factor')
    imgN = imgN.astype(np.uint8)

    return imgN, (imin, imax)

def reduceBuffer(Ibuff, is_light_background):
    if is_light_background:
        return np.min(Ibuff, axis=0)
    else:
        return np.max(Ibuff, axis=0)

def createImgGroup(fid, name, tot_frames, im_height, im_width, is_expandable=True):
    parentnode, _, name = name.rpartition('/')
    parentnode += '/'

    if is_expandable:
        img_dataset = fid.create_earray(
                        parentnode,
                        name,
                        atom=tables.UInt8Atom(),
                        shape =(0,
                             im_height,
                             im_width),
                        chunkshape=(1,
                             im_height,
                             im_width),
                        expectedrows=tot_frames,
                        filters=TABLE_FILTERS
                        )
    else:
        img_dataset = fid.create_carray(
                        parentnode,
                        name,
                        atom=tables.UInt8Atom(),
                        shape =(tot_frames,
                             im_height,
                             im_width),
                        filters=TABLE_FILTERS
                        )

    img_dataset._v_attrs["CLASS"] = np.string_("IMAGE")
    img_dataset._v_attrs["IMAGE_SUBCLASS"] = np.string_("IMAGE_GRAYSCALE")
    img_dataset._v_attrs["IMAGE_WHITE_IS_ZERO"] = np.array(0, dtype="uint8")
    img_dataset._v_attrs["DISPLAY_ORIGIN"] = np.string_("UL")  # not rotated
    img_dataset._v_attrs["IMAGE_VERSION"] = np.string_("1.2")

    return img_dataset

def initMasksGroups(fid, expected_frames, im_height, im_width,
    attr_params, save_full_interval, is_expandable=True):

    # open node to store the compressed (masked) data
    mask_dataset = createImgGroup(fid, "/mask", expected_frames, im_height, im_width, is_expandable)

    tot_save_full = (expected_frames // save_full_interval) + 1
    full_dataset = createImgGroup(fid, "/full_data", tot_save_full, im_height, im_width, is_expandable)
    full_dataset._v_attrs['save_interval'] = save_full_interval

    assert all(x in ['expected_fps', 'is_light_background', 'microns_per_pixel'] for x in attr_params)
    set_unit_conversions(mask_dataset, **attr_params)
    set_unit_conversions(full_dataset, **attr_params)

    if is_expandable:
        mean_intensity = fid.create_earray('/',
                                        'mean_intensity',
                                        atom=tables.Float32Atom(),
                                        shape=(0,),
                                        expectedrows=expected_frames,
                                        filters=TABLE_FILTERS)
    else:
        mean_intensity = fid.create_carray('/',
                                        'mean_intensity',
                                        atom=tables.Float32Atom(),
                                        shape=(expected_frames,),
                                        filters=TABLE_FILTERS)

    return mask_dataset, full_dataset, mean_intensity

def compressVideo(video_file, masked_image_file, mask_param,  expected_fps=25,
                  microns_per_pixel=None, bgnd_param ={}, buffer_size=-1,
                  save_full_interval=-1, max_frame=1e32, is_extract_timestamp=False,
                  fovsplitter_param={}):
    '''
    Compresses video by selecting pixels that are likely to have worms on it and making the rest of
    the image zero. By creating a large amount of redundant data, any lossless compression
    algorithm will dramatically increase its efficiency. The masked images are saved as hdf5 with gzip compression.
    The mask is calculated over a minimum projection of an image stack. This projection preserves darker regions
    (or brighter regions, in the case of fluorescent labelling)
    where the worm has more probability to be located. Additionally it has the advantage of reducing
    the processing load by only requiring to calculate the mask once per image stack.
     video_file --  original video file
     masked_image_file --
     buffer_size -- size of the image stack used to calculate the minimal projection and the mask
     save_full_interval -- have often a full image is saved
     max_frame -- last frame saved (default a very large number, so it goes until the end of the video)
     mask_param -- parameters used to calculate the mask
    '''

    #get the default values if there is any bad parameter
    output = compress_defaults(masked_image_file,
                                expected_fps,
                                buffer_size = buffer_size,
                                save_full_interval = save_full_interval)

    buffer_size = output['buffer_size']
    save_full_interval = output['save_full_interval']

    if len(bgnd_param) > 0:
        is_bgnd_subtraction = True
        assert bgnd_param['buff_size']>0 and bgnd_param['frame_gap']>0
    else:
        is_bgnd_subtraction = False

    if len(fovsplitter_param) > 0:
        is_fov_tosplit = True
        assert all(key in fovsplitter_param for key in ['total_n_wells', 'whichsideup', 'well_shape'])
        assert fovsplitter_param['total_n_wells']>0
    else:
        is_fov_tosplit = False

    # processes identifier.
    base_name = masked_image_file.rpartition('.')[0].rpartition(os.sep)[-1]

    # select the video reader class according to the file type.
    vid = selectVideoReader(video_file)

    # delete any previous  if it existed
    with tables.File(masked_image_file, "w") as mask_fid:
        pass

    #Extract metadata
    if is_extract_timestamp:
        # extract and store video metadata using ffprobe
        #NOTE: i cannot calculate /timestamp until i am sure of the total number of frames
        print_flush(base_name + ' Extracting video metadata...')
        expected_frames = store_meta_data(video_file, masked_image_file)

    else:
        expected_frames = 1

    # Initialize background subtraction if required

    if is_bgnd_subtraction:
        print_flush(base_name + ' Initializing background subtraction.')
        bgnd_subtractor = BackgroundSubtractorVideo(video_file, **bgnd_param)

    # intialize some variables
    max_intensity, min_intensity = np.nan, np.nan
    frame_number = 0
    full_frame_number = 0
    image_prev = np.zeros([])

    # Initialise FOV splitting if needed
    if is_bgnd_subtraction:
        img_fov = bgnd_subtractor.bgnd.astype(np.uint8)
    else:
        ret, img_fov = vid.read()
        # close and reopen the video, to restart from the beginning
        vid.release()
        vid = selectVideoReader(video_file)

    if is_fov_tosplit:
        # TODO: change class creator so it only needs the video name? by using
        # Tierpsy's functions such as selectVideoReader it can then read the first image by itself

        camera_serial = parse_camera_serial(masked_image_file)

        fovsplitter = FOVMultiWellsSplitter(img_fov,
                                            camera_serial=camera_serial,
                                            px2um=microns_per_pixel,
                                            **fovsplitter_param)
        wells_mask = fovsplitter.wells_mask
    else:
        wells_mask = None

    # initialize timers
    print_flush(base_name + ' Starting video compression.')

    if expected_frames == 1:
        progressTime = TimeCounter('Compressing video.')
    else:
        #if we know the number of frames display it in the progress
        progressTime = TimeCounter('Compressing video.', expected_frames)

    with tables.File(masked_image_file, "r+") as mask_fid:

        #initialize masks groups
        attr_params = dict(
            expected_fps = expected_fps,
            microns_per_pixel = microns_per_pixel,
            is_light_background = int(mask_param['is_light_background'])
            )
        mask_dataset, full_dataset, mean_intensity = initMasksGroups(mask_fid,
            expected_frames, vid.height, vid.width,
            attr_params, save_full_interval)

        if is_bgnd_subtraction:
            bg_dataset = createImgGroup(mask_fid, "/bgnd", 1, vid.height, vid.width, is_expandable=False)
            # because we only save the one background:
            bg_dataset._v_attrs['save_interval'] = len(vid)
            # except that if reading with ffmpeg, this could be not accurate.
            # call it again after reading the whole video!
            bg_dataset[0,:,:] = img_fov

        if vid.dtype != np.uint8:
            # this will worm as flags to be sure that the normalization took place.
            normalization_range = mask_fid.create_earray('/',
                                        'normalization_range',
                                        atom=tables.Float32Atom(),
                                        shape=(0, 2),
                                        expectedrows=expected_frames,
                                        filters=TABLE_FILTERS
                                        )

        while frame_number < max_frame:

            ret, image = vid.read()
            if ret != 0:
                # increase frame number
                frame_number += 1

                # opencv can give an artificial rgb image. Let's get it back to
                # gray scale.
                if image.ndim == 3:
                    image = cv2.cvtColor(image, cv2.COLOR_RGB2GRAY)

                if image.dtype != np.uint8:
                    # normalise image intensities if the data type is other
                    # than uint8
                    image, img_norm_range = normalizeImage(image)
                    normalization_range.append(img_norm_range)

                #limit the image range to 1 to 255, 0 is a reserved value for the background
                assert image.dtype == np.uint8
                image = np.clip(image, 1,255)

                # Add a full frame every save_full_interval
                if frame_number % save_full_interval == 1:
                    full_dataset.append(image[np.newaxis, :, :])
                    full_frame_number += 1

                # buffer index
                ind_buff = (frame_number - 1) % buffer_size

                # initialize the buffer when the index correspond to 0
                if ind_buff == 0:
                    Ibuff = np.zeros(
                        (buffer_size, vid.height, vid.width), dtype=np.uint8)

                # add image to the buffer
                Ibuff[ind_buff, :, :] = image.copy()
                mean_int = np.mean(image)
                assert mean_int >= 0
                mean_intensity.append(np.array([mean_int]))

            else:
                # sometimes the last image is all zeros, control for this case
                if np.all(Ibuff[ind_buff] == 0):
                    frame_number -= 1
                    ind_buff -= 1

                # close the buffer
                Ibuff = Ibuff[:ind_buff + 1]

            # mask buffer and save data into the hdf5 file
            if (ind_buff == buffer_size - 1 or ret == 0) and Ibuff.size > 0:
                if is_bgnd_subtraction:
                    Ibuff_b  = bgnd_subtractor.apply(Ibuff, frame_number)
                else:
                    Ibuff_b = Ibuff

                #calculate the max/min in the of the buffer
                img_reduce = reduceBuffer(Ibuff_b, mask_param['is_light_background'])

                mask = getROIMask(img_reduce, wells_mask=wells_mask, **mask_param)

                Ibuff *= mask

                # now apply the well_mask if is MWP
                if is_fov_tosplit:
                    fovsplitter.apply_wells_mask(Ibuff) # Ibuff will be modified after this

                # add buffer to the hdf5 file
                frame_first_buff = frame_number - Ibuff.shape[0]
                mask_dataset.append(Ibuff)

            if frame_number % 500 == 0:
                # calculate the progress and put it in a string
                progress_str = progressTime.get_str(frame_number)
                print_flush(base_name + ' ' + progress_str)

            # finish process
            if ret == 0:
                break

        # now that the whole video is read, we definitely have a better estimate
        # for its number of frames. so set the save_interval again
        if is_bgnd_subtraction:
            # bg_dataset._v_attrs['save_interval'] = len(vid)
            # so that didn't work. Either I have an off by one,
            # or if the video is corrupted it's just safer to do:
            bg_dataset._v_attrs['save_interval'] = mask_dataset.shape[0]

        # close the video
        vid.release()

    # save fovsplitting data
    if is_fov_tosplit:
        fovsplitter.write_fov_wells_to_file(masked_image_file)
        if fovsplitter.is_dubious:
            print(f'Check {masked_image_file} for plate alignment')

    read_and_save_timestamp(masked_image_file)
    print_flush(base_name + ' Compressed video done.')


text:
import numpy as np
from tectosaur.ops.dense_op import DenseOp
from tectosaur.ops.composite_op import CompositeOp
from tectosaur.ops.sum_op import SumOp

def test_composite_op():
    n = 20
    split_n = 10
    mat = np.random.rand(n, n)
    v = np.random.rand(n)
    comp_op = CompositeOp(
        (DenseOp(mat[:split_n,:split_n]), 0, 0),
        (DenseOp(mat[split_n:,:split_n]), split_n, 0),
        (DenseOp(mat[:split_n,split_n:]), 0, split_n),
        (DenseOp(mat[split_n:,split_n:]), split_n, split_n)
    )

    correct = mat.dot(v)
    check = comp_op.dot(v)
    np.testing.assert_almost_equal(correct, check)

def test_sum_op():
    mat = np.random.rand(2,2)
    v = np.random.rand(2)
    sum_op = SumOp([DenseOp(mat), DenseOp(mat)])
    np.testing.assert_almost_equal(sum_op.dot(v), 2 * mat.dot(v))

max_stars_repo_name: jlmaurer/tectosaur

max_issues_repo_name: jlmaurer/tectosaur

max_forks_repo_name: jlmaurer/tectosaur


text:
# import itertools
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from postladim import ParticleFile
from gridforce_discrete import Grid

# Input file
particle_file = "obstacle.nc"

# Get the grid information
grid = Grid(None)

# Open the particle_file
pf = ParticleFile(particle_file)
num_times = pf.num_times

# Set up the plot area
fig = plt.figure(figsize=(12, 8))
ax = plt.axes(xlim=(0, grid.imax), ylim=(0, grid.jmax), aspect="equal")

# Landmask
Xb = np.arange(-0.5, grid.imax)
Yb = np.arange(-0.5, grid.jmax)
constmap = plt.matplotlib.colors.ListedColormap([0.2, 0.6, 0.4])
M = np.ma.masked_where(grid.M > 0, grid.M)
plt.pcolormesh(Xb, Yb, M, cmap=constmap)
# Draw the cirular boundary
T = np.linspace(0, np.pi)
plt.plot(grid.X0 + grid.R * np.cos(T), grid.R * np.sin(T), color="black")

# Plot initial particle distribution
X, Y = pf.position(0)
particle_dist, = ax.plot(X, Y, ".", color="red", markeredgewidth=0.5, lw=0.5)
time0 = pf.time(0)  # Save start-time
timestr = "00:00"
timestamp = ax.text(0.02, 0.93, timestr, fontsize=16, transform=ax.transAxes)

# Update function
def animate(t):
    X, Y = pf.position(t)
    particle_dist.set_data(X, Y)
    # Time since start in minutes
    # dtime = int((pf.time(t) - time0).total_seconds() / 60)
    dtime = (pf.time(t) - time0) / 60
    dtimestr = str(dtime)
    timestamp.set_text(dtimestr)
    return particle_dist, timestamp

# Do the animation
anim = FuncAnimation(
    fig,
    animate,
    frames=num_times,
    interval=20,
    repeat=True,
    repeat_delay=500,
    blit=True,
)

plt.show()

max_stars_repo_name: pnsaevik/ladim

max_issues_repo_name: pnsaevik/ladim

max_forks_repo_name: pnsaevik/ladim


text:
#Copyright (c) 2018-2020 Gavin W. Wilson

#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 .snvs import calc
from collections import Counter
from .plt import format_subplots, rand_jitter
from scipy.sparse import csr_matrix
import pylab as plt
import pandas as pd
from matplotlib.lines import Line2D
import numpy as NP
from sklearn.metrics import roc_curve, auc
import matplotlib
from matplotlib.colors import ListedColormap, Normalize
from matplotlib.collections import PatchCollection
from scipy.stats import binom
from cyvcf2 import VCF
import operator
from mpl_toolkits.axes_grid1 import make_axes_locatable
import scipy.cluster.hierarchy as sch
import scipy.spatial.distance as dist
from mpl_toolkits.axes_grid1.axes_divider import make_axes_locatable

def iter_snvs(ref, alt):
    for i, (s, e) in enumerate(zip(ref.indptr, ref.indptr[1:])):
        yield i, ref.indices[s:e], s, e, ref.data[s:e], alt.data[s:e]

class FindSNVs(object):
    def __init__(self, ref, alt, snvs, barcodes, labels, close_fig=True):
        self._ref = ref
        self._alt = alt
        self._snvs = snvs
        self._cells = ref.shape[1]
        self._close_fig = close_fig
        self._barcodes = barcodes
        self._labels = labels

    def score_snvs(self, snv_min_cells_per_genotype=2, snv_total_cells=5, snv_AF_diff = 0.3, top_snvs=5):
        self._params = {
            'snv_min_cells_per_genotype':snv_min_cells_per_genotype,
            'snv_total_cells':snv_total_cells,
            'snv_AF_diff':snv_AF_diff,
            #'eps':eps,
            'top_snvs':top_snvs,
            #'snv_score_diff':score_diff,
            #'snv_score_perc':score_perc
        }
        labels = self._labels
        ogts = calc.calculate_genotypes(self._ref, self._alt, labels.astype('int32'))
        x, y = ogts['cells_0'], ogts['cells_1']
        p = self._params
        idx = NP.where((x >= p['snv_min_cells_per_genotype']) & (y >= p['snv_min_cells_per_genotype']) & ((x + y) >= p['snv_total_cells']))[0]
        keep = idx[(NP.abs(ogts['mean_AF_0'][idx] - ogts['mean_AF_1'][idx]) >= p['snv_AF_diff'])]
        snv_scores, predictive, snv_index = self._scores(ogts, self._ref, self._alt, keep, labels=labels)
        self._keep = keep
        self._snv_scores = snv_scores
        self._predictive = predictive
        self._snv_index = snv_index

    def find_predictive(self, N = 5, max_snvs = 50, use_genotypes=False):
        totals = self._snv_scores['total_correct'].values
        gts = None
        if use_genotypes:
            gts = self._snv_scores['exome_gt1'].values
        ai = NP.argsort(totals)
        C = self._predictive.shape[1]
        cell_counts = NP.zeros(C, dtype='int')
        labs = self._labels
        glabs = NP.where(labs < 2)[0]
        pp = self._predictive
        snvs = []
        for i in ai[::-1]:
            if len(snvs) > max_snvs:
                break
            if use_genotypes and gts[i] == '':
                continue
            s, e = pp.indptr[i], pp.indptr[i + 1]
            ci = pp.indices[s:e]
            lidx = ci[labs[ci] < 2]
            updated = (cell_counts[lidx] < N).sum()
            if updated > 0:
                cell_counts[ci] += 1
                finished = (cell_counts[glabs] >= N).sum()
                snvs.append(i)
                if finished >= len(glabs):
                    break

        self._predictive_snvs = self._snv_scores.iloc[snvs]
        self._predictive_counts = cell_counts
        print(f'Needed {len(snvs)} SNVs to cover each singlet {N} times')

    def dotplot(self, label_map, out = None, close=False, use_genotypes=False):
        snvs = self._predictive_snvs
        idx = snvs['snv_idx'].values
        ref, alt = self._ref[idx], self._alt[idx]

        snames = NP.array([f'{r}: {p + 1}' for r, p in zip(snvs.chrom, snvs.pos)])

        good_index = [i for i, (b, l) in enumerate(zip(self._barcodes, self._labels)) if b in label_map and l < 3]
        print(f'Using {len(good_index)} out of {len(self._barcodes)} barcodes')
        labels_raw = [label_map[b] for b in self._barcodes[good_index]]
        label_names = NP.array(list(set(labels_raw)) + ['Doublets'])
        label_imap = {l:i for i,l in enumerate(label_names)}
        labels = NP.fromiter((label_imap[l] for l in labels_raw), dtype='int')
        ref, alt = ref[:,good_index], alt[:,good_index]

        NL = len(label_imap)
        gts = self._labels[good_index]

        labels[gts > 1] = NL - 1
        groups = sorted(Counter(zip(gts, labels)).items(), key=lambda k: (k[0][0], -k[1]), reverse=True)
        print(groups)
        gcounts = NP.fromiter((g[1] for g in groups), dtype='int')
        groups = NP.array([g[0] for g in groups])

        gmap = {tuple(d):i for i, d in enumerate(groups)}

        cafs = NP.zeros((len(groups), len(idx)), dtype='double')
        ctots = NP.zeros((len(groups), len(idx)), dtype='int')

        cell_rows = NP.fromiter((gmap[(g, l)] for g, l in zip(gts, labels)), dtype='int')

        for j, barcodes, start, end, r, a in iter_snvs(ref, alt):
            #print(cell_rows[barcodes])
            NP.add.at(ctots[:,j], cell_rows[barcodes], 1)
            AF = a / (r + a)
            NP.add.at(cafs[:,j], cell_rows[barcodes], AF)

        cperc = ctots / gcounts[:,NP.newaxis]
        R = (cperc * 0.90)  / 2
        A = NP.zeros(R.shape, dtype='double')
        A[ctots > 0] = cafs[ctots > 0] / ctots[ctots > 0]

        d1 = dist.pdist(A.T, metric='euclidean')
        rlink = sch.linkage(d1, method='ward') ### gene-clustering metric - 'average', 'single', 'centroid', 'complete'
        cidx = sch.leaves_list(rlink)

        A = A[:,cidx]
        R = R[:,cidx]

        fig, ax = format_subplots(1, 1, figsize=(8, 5), dpi=1000)
        ax.set_aspect('equal')

        # SNV's along the X axis
        x, y = NP.meshgrid(NP.arange(R.shape[1]), NP.arange(R.shape[0]))
        circles = NP.array([plt.Circle((j,i), radius=r) for r, j, i in zip(R.flat, x.flat, y.flat)])
        norm = matplotlib.colors.Normalize(vmin=0.,vmax=1.)
        col = PatchCollection(circles, array=A.flatten(), cmap=plt.cm.Spectral_r, norm=norm, linewidth=0.5, edgecolor='k')
        ax.add_collection(col)
        ax.set_ylim(-0.5, R.shape[0] - 0.5)
        ax.set_xlim(-0.5, R.shape[1] - 0.5)
        ax.set_yticks(NP.arange(R.shape[0]))
        ax.set_xticks(NP.arange(R.shape[1]))
        ax.set_xticklabels(snames[cidx], fontsize=10, rotation=90)
        labs = [f'{l} N = {c}' for l, c in zip(label_names[groups[:,1]], gcounts)]
        ax.set_yticklabels(labs)
        ax.axhline(NP.where(groups[:,0] == 1)[0][0] - 0.5, color='k')
        ax.axhline(NP.where(groups[:,0] == 0)[0][0] - 0.5, color='k')
        ax.spines['right'].set_visible(True)
        ax.spines['top'].set_visible(True)

        #for g, c in zip(groups, gcounts):
        #    print(g[0], label_names[g[1]], c)

        g1 = NP.where(groups[:,0] == 0)[0]
        g2 = NP.where(groups[:,0] == 1)[0]
        print(g1, g2)
        trans = ax.get_yaxis_transform()
        ax.text(1.01, NP.median(g1), 'GT1', va='center', transform=trans, rotation=270, fontsize=15)
        ax.text(1.01, NP.median(g2), 'GT2', va='center', transform=trans, rotation=270, fontsize=15)

        pos = ax.get_position()

        W = pos.xmax - pos.xmin
        H = pos.ymax - pos.ymin
        WR = W / A.shape[1]
        HR = H / A.shape[0]
        cb_lgd = fig.add_axes([pos.xmax * 1.15, pos.ymax * 0.6, W / A.shape[1], H / A.shape[0] * 4 ])
        lcircles = [plt.Circle((0, i), radius=((r * 0.95) / 2), fill='k') for i, r in enumerate([0.25, 0.50, 0.75, 1])]
        lcol = PatchCollection(lcircles, array=NP.array([1] * len(lcircles)), cmap="binary", norm=norm, linewidth=0.5)
        cb_lgd.add_collection(lcol)
        cb_lgd.set_ylim(-0.5, 3.5)
        cb_lgd.set_xlim(-0.5, 0.5)
        cb_lgd.set_aspect('equal')
        cb_lgd.set_xticks([])
        cb_lgd.set_yticks([0, 1, 2, 3])
        cb_lgd.set_yticklabels([25, 50, 75, 100])
        cb_lgd.set_ylabel('Cells Expressing SNV (%)')

        pos = cb_lgd.get_position()
        cb_ax = fig.add_axes([pos.xmax * 1.1, pos.ymin, 0.05, 0.25 ])
        fig.colorbar(col, cax=cb_ax)
        cb_ax.yaxis.tick_left()
        cb_ax.set_ylabel('Mean Allele Fraction')
        cb_ax.yaxis.set_label_position("left")

        if use_genotypes:
            axd = make_axes_locatable(ax)
            gt_ax = axd.append_axes("top", size="20%", pad="1%")
            #gt_ax = fig.add_axes([pos.xmin, pos.ymax, W, 0.25 ])

            gt_labs = NP.zeros((2, len(snvs)), dtype='double')
            gt_labs[0] = [float(s) for s in snvs['exome_gt1'].values]
            gt_labs[1] = [float(s) for s in snvs['exome_gt2'].values]

            gt_ax.imshow(gt_labs[:,cidx], vmin=0, vmax=1, cmap=plt.cm.Spectral_r, aspect='auto')
            gt_ax.set_xticks([])
            gt_ax.set_yticks([0, 1])
            gt_ax.set_yticklabels(['GT1 Exome GT', 'GT2 Exome GT'], fontsize=8)

        if out is not None:
            fig.savefig(out, bbox_inches='tight')

        if close:
            plt.close(fig)

    def add_genotypes(self, vcf, gt1, gt2):
        fo = VCF(vcf, samples=[gt1, gt2])
        s1 = fo.samples.index(gt1)
        s2 = fo.samples.index(gt2)
        gmap = {0:0.0, 1:0.5, 3:1.0, 2:None}
        exome = {}
        for r in fo:
            if r.FILTER or r.is_indel or len(r.ALT) > 1:
                continue
            pos = r.POS - 1
            ref = r.REF
            alt = r.ALT[0]
            chrom = r.CHROM
            g1 = r.gt_types[s1]
            g2 = r.gt_types[s2]
            if g1 == g2 or (g1 == 2 or g2 == 2):
                continue
            exome[chrom, pos] = (ref, alt, gmap[g1], gmap[g2])
        dgt1 = []
        dgt2 = []
        for i, r in self._snv_scores.iterrows():
            t = (r.chrom, r.pos)
            d = exome.get(t)
            if d is not None and d[0] == r['ref'] and d[1] == r['alt']:
                dgt1.append(d[2])
                dgt2.append(d[3])
            else:
                dgt1.append('')
                dgt2.append('')
        self._snv_scores['exome_gt1'] = dgt1
        self._snv_scores['exome_gt2'] = dgt2

    def _scores(self, gts, ref, alt, keep, labels, NL = 2):
        p = self._params
        #rr = ref[keep].tocsc()
        #aa = alt[keep].tocsc()
        rr = ref[keep]
        aa = alt[keep]
        snvs = self._snvs.iloc[keep]
        N = p['top_snvs']
        means = [gts[f'mean_AF_{l}'][keep] for l in range(NL)]
        #eps = p['eps']
        #eps1 = NP.log(eps)
        #eps2 = NP.log(1 - eps)
        #smin = p['snv_score_diff']
        #pmin = p['snv_score_perc']
        #predictive = NP.full(len(rr.data), -1, dtype='int8')
        rows = []

        data = []
        indptr = [0]
        indices = []
        snv_index = []

        for i, barcodes, start, end, r, a in iter_snvs(rr, aa):
            if NP.isnan(means[0][i]) or NP.isnan(means[1][i]):
                continue
            snv = snvs.iloc[i]
            labs = labels[barcodes]
            lidx = NP.where(labs < 2)[0]
            ll = labels[barcodes[lidx]]

            AF = a[lidx] / (a + r)[lidx]
            m1 = means[0][i]
            m2 = means[1][i]

            plab = 1 if (m1 < m2) else 0
            fpr, tpr, thresholds = roc_curve(ll, AF, pos_label=plab) #, drop_intermediate=False)
            #a = NP.abs(auc(fpr, tpr) - 0.5) * 2
            a = auc(fpr, tpr)
            thresh = thresholds[NP.argmax(tpr - fpr)]

            g1c = (ll == plab) & (AF >= thresh)
            gt1c = g1c.sum()
            g2c = (ll != plab) & (AF < thresh)
            gt2c = g2c.sum()
            #predictive[lidx + start] = (g2c | g1c)

            kidx = NP.where(g2c | g1c)[0]
            indices.extend(barcodes[lidx[kidx]])
            #indices.extend(barcodes[lidx])
            #data.extend((g2c | g1c))
            data.extend(NP.ones(len(kidx), dtype='bool'))
            indptr.append(indptr[-1] + len(kidx))
            #indptr.append(indptr[-1] + len(lidx))
            snv_index.append(i)

            if plab == 1:
                gt1c, gt2c = gt2c, gt1c

            rows.append((keep[i], snvs.index[i][0], snvs.index[i][1], snv['ref'], snv['alt'], (ll == 0).sum(), (ll == 1).sum(), m1, m2, a, thresh, gt1c, gt2c, gt1c + gt2c))

        df = pd.DataFrame(rows, columns=('snv_idx', 'chrom', 'pos', 'ref', 'alt', 'gt1', 'gt2', 'mgt1', 'mgt2', 'auc', 'threshold', 'gt1_correct', 'gt2_correct', 'total_correct'))

        predictive = csr_matrix((data, indices, indptr), shape=(len(snv_index), ref.shape[1]))
        print(predictive.shape, len(predictive.data), len(predictive.indptr), len(snv_index))
        return df, predictive, NP.array(snv_index)

max_stars_repo_name: GWW/sctx

max_issues_repo_name: GWW/sctx

max_forks_repo_name: GWW/sctx

reason: import numpy,import scipy,from scipy


text:
"""
ace.solvers
~~~~~~~~~~~

Solve a simple linear equation Ax = b in a particular way
"""
from numpy import append, array, diagonal, tril, triu
from numpy.linalg import inv, solve
from scipy.linalg import lu

def lu_decomposition(a, b):
    """
    Solve a linear equation by LU-decomposition

    Comes from LU decomposition of a matrix A s.t. A = LU

    Then

        LUx = b => Ux = y => Ly = b
    """
    _, l, u = lu(a)
    y = solve(l, b)
    return solve(u, y)

def gauss_seidel(a, b, iterations, x=None):
    """
    Solve a linear equation by the gauss seidel iteration outlined in the book

    Follows the eq:

        x = inv(L)*(b - U*x)
    """
    l = tril(a)
    upper_plus_diagonal = triu(a)
    u = upper_plus_diagonal - _diagonal_matrix(a)
    x = _check_for_initial_guess(a, x)
    for _ in xrange(iterations):
        x = inv(l).dot(b - u.dot(x))
    return x

def gauss_jacobi(a, b, iterations, x=None):
    """
    Solve a linear equation by the gauss jacobi iteration outlined in the book.

    Follows the eq:

        x = inv(D)(b - Rx)

    Where D is the diagonal matrix of A and R is the remainder s.t D + R = A
    """
    d = _diagonal_matrix(a)
    # Calculate the remainder matrix
    r = a - d
    x = _check_for_initial_guess(a, x)
    for _ in xrange(iterations):
        x = inv(d).dot(b - r.dot(x))
    return x

def _diagonal_matrix(a):
    """
    Given a square, 2D matrix a, create a diagonal matrix from it
    """
    diag = diagonal(a)
    # Create first row to initialize the correct shape
    first_row = [diag[0]] + [0 for _ in range(len(diag) - 1)]
    diag_matrix = array([first_row])
    # Construct the remaining rows in the diagonal matrix
    for index in range(1, len(diag)):
        row = [0 if index != i else diag[i] for i in range(len(diag))]
        diag_matrix = append(diag_matrix, [row], axis=0)
    return diag_matrix

def _check_for_initial_guess(a, x):
    """
    If we have not provided an initial array for x make a new one
    """
    if not x:
        x = array([1 for _ in range(a.shape[1])])
    return x

max_stars_repo_name: hahnicity/ace

max_issues_repo_name: lnsongxf/ace

max_forks_repo_name: lnsongxf/ace

reason: from numpy,from scipy


text:
import tensorflow as tf
import numpy as np
from models import utils
import json

class CDAE():
    def __init__(self, input_class, sparse=True):
        # Input Data Properties
        print 'Importing Data ...'
        self.input_data = input_class(sparse=sparse)
        print 'Done.'
        self.sparse = sparse
        self.n_items = self.input_data.get_n_items()
        self.n_users = self.input_data.get_n_users()

    def _set_up_training_ops(self, **kwargs):
        '''
        '''
        print 'Set up Training Operations...'
        # Init parameters
        dropout_prob = kwargs.get('dropout_prob') if kwargs.get('dropout_prob') else 0.2
        transfer_fn = kwargs.get('transfer_fn') if kwargs.get('transfer_fn') else tf.sigmoid
        output_fn = kwargs.get('output_fn') if kwargs.get('output_fn') else tf.identity
        n_hidden = kwargs.get('n_hidden') if kwargs.get('n_hidden') else 50

        # AdaGrad Parameters
        ada_beta = kwargs.get('ada_beta') if kwargs.get('ada_beta') else 1.0
        ada_learning_rate = kwargs.get('ada_learning_rate') if kwargs.get('ada_learning_rate') else 0.01
        #optimizer = tf.train.AdagradOptimizer(ada_learning_rate, ada_beta)
        optimizer = tf.train.GradientDescentOptimizer(ada_learning_rate)

        # Regularization Parameter
        reg_lambda = kwargs.get('reg_lambda') if kwargs.get('reg_lambda') else 1.0

        # Model Ops
        self.y = tf.placeholder(tf.float32, [None, self.n_items])
        self.z = lambda u: transfer_fn(tf.add_n( [tf.matmul(self.y, self.weight['W']),
                                                  tf.reshape(self.weight['V'][u], [1, n_hidden]),
                                                  tf.reshape(self.weight['b'], [1, n_hidden])]
                                                ))
        self.y_hat = lambda u: output_fn(tf.add(tf.matmul(self.z(u), tf.transpose(self.weight['W_prime'])),
                                                tf.reshape(self.weight['b_prime'], [1, self.n_items]))
                                        )
        # Loss Function
        # TODO - Figure out how to include logistic loss efficiently
        self.loss = lambda u: tf.nn.l2_loss(tf.sub(self.y, self.y_hat(u)))

        # Gradient Ops
        #   - To my knowledge this when cost_gradient_apply gets called per user per item
        #     we are simulating online training. However, if we dont call apply_gradients
        #     we may be able to accumalate the gradients and apply them on a per user basis
        #     to simulate something like batch training.
        #   - TODO compare results of both online and batch methods

        self.dLoss_d = lambda u, var: optimizer.compute_gradients(self.loss(u), var_list=[var])[0][0]
        self.dCost_d = lambda u, var: tf.mul(1.0 / self.n_users, self.dLoss_d(u, var)) - tf.mul(reg_lambda, var)
        self.cost_gradient_apply = lambda u, var: optimizer.apply_gradients([[self.dCost_d(u, var), var]])
        self.init_op = tf.initialize_all_variables()
        print 'Done.'

    def _create_augmented_O_set(self, ith_data, neg_sample_ratio = 5):
        '''
        '''
        pos_examples_indices = np.where(ith_data > 0)[0]
        neg_examples_indices = np.where(ith_data == 0)[0]
        aug_set = {index for index in pos_examples_indices}
        if len(neg_examples_indices) <= neg_sample_ratio * len(pos_examples_indices):
            for index in neg_examples_indices:
                aug_set.add(index)
        else:
            for index in np.random.permutation(neg_examples_indices)[0:neg_sample_ratio * len(pos_examples_indices)]:
                aug_set.add(index)
        return aug_set

    def train(self, **kwargs):
        '''
            Method
        '''
        MAX_ITERATION = kwargs.get('max_iteration') if kwargs.get('max_iteration') else 1000
        dropout_prob = kwargs.get('dropout_prob') if kwargs.get('dropout_prob') else 0.2
        n_hidden = kwargs.get('n_hidden') if kwargs.get('n_hidden') else 50
        dropout_prob = tf.constant(dropout_prob, dtype=tf.float32)
        iteration = 0

        self.weight = self._init_variables(n_hidden)
        self._set_up_training_ops(**kwargs)

        old_avg_cost = None
        model_improved = True
        input_data = self.input_data.get_data()
        sess = tf.Session()
        print 'Initialing all Variables...'
        sess.run(self.init_op)
        print 'Done.'
        print 'Start Training... '
        while iteration < MAX_ITERATION or model_improved:
            old_avg_cost = avg_cost if iteration >= 1 else None
            sum_cost = 0
            for u in xrange(self.n_users):
                val_data = input_data[u].get_full_data()
                train_data = input_data[u].get_train_data()

                if self.sparse:
                    try:
                        val_data = sess.run(tf.sparse_tensor_to_dense(val_data))
                        train_data = sess.run(tf.sparse_tensor_to_dense(train_data))
                    except:
                        print 'user: {}'.format(u)
                y_tilde = sess.run(tf.nn.dropout(tf.constant(train_data, dtype=tf.float32, shape=[1, self.n_items]), dropout_prob))
                augmented_O_set = self._create_augmented_O_set(np.array(val_data))

                for item in augmented_O_set:
                    sess.run(self.cost_gradient_apply(u, self.weight['W_prime'][item]), feed_dict={self.y: y_tilde})
                    sess.run(self.cost_gradient_apply(u, self.weight['b_prime'][item]), feed_dict={self.y: y_tilde})
                for item in np.where(x_tilde == 1):
                    sess.run(self.cost_gradient_apply(u, self.weight['W'][item]), feed_dict={self.y: y_tilde})

                sess.run(self.cost_gradient_apply(u, self.weight['V'][u]), feed_dict={self.y: y_tilde})
                sess.run(self.cost_gradient_apply(u, self.weight['b']), feed_dict={self.y: y_tilde})
                sum_cost += sess.run(self.loss(u), feed_dict={self.y: val_data})

            avg_cost = sum_cost / self.n_users
            model_improved = True if old_avg_cost == None or old_avg_cost > avg_cost else False
            print 'Iteration: {}                    Average Cost per User: {}'.format(i, cost)
            iteration += 1
        sess.close()

    def recommend(self, u, Y, N=5, test_set=False):
        '''
           Method recommend outputs the top-N item recommendations for
           user u

           Inputs:
               u (Int): id for user to be recommended
               Y (AbstractInputSource): user data class that encapsulates
                                        training, test, and full data
               N (Int): number of items to recommend user
               test_set (Boolean): necessary flag for checking the accuracy
                                   of the training set to the test set
           Output:
               List: Top N list of items recommend to user u
        '''

        if test_set:
            if not isinstance(Y, AbstractInputClass):
                raise ValueError
            recommendation_candidates = Y.get_neg_indices()
            if self.sparse:
                train_data = self.sess.run(tf.sparse_tensor_to_dense(Y.get_train_data()))
            else:
                train_data = Y.get_train_data()

            feed_dict = {self.y: train_data}
        else:
            recommendation_candidates = np.where(Y == 0)
            feed_dict = {self.y: Y}
        y_hat = self.sess.run(self.y_hat(u), feed_dict=feed_dict)
        top_recommended_items = sorted(recommendation_candidates, key=lambda i: y_hat[i], reverse=True)
        return top_items_sorted[:N] if len(top_items_sorted) > N else top_items_sorted

    def _save_weights(self, with_error=False):
        output_dict = {}
        if with_error:
            # Mark the file as distinct
            pass
        for weight in self.weight.items():
            weight = self.sess.run(weight)

    def _load_weights(self):
        pass
    def _count_num_accurate_classification(self, test_set_items, recommendations):
        '''
            | Intersection(C-adopted, C-NRecommended) |
            Inputs:
                test_set_items (List): Actual items that user liked
                recommendation (List): Items the system recommended to the user
            Output:
                Int: Length of intersection between system recommended items and
                     items the user actually liked
        '''
        num_accurate_predictions = 0
        for item in recommendations:
            if item in test_set_items:
                num_accurate_predictions += 1
        return num_accurate_predictions

    def precision(self, actual, predicted, N=5):
        '''
            Calculates precision defined by:
                | Intersection(C-adopted, C-NRecommended) | / N

            Inputs:
                actual (List): Actual items that were liked
                predicted (List): Prediction of actual items
            Output:
                Float: Number b/w 0 and 1 corresponding to the precision
                       of the recommendation list
        '''

        if len(predicted) > N:
            predicted = predict[:N]

        num_accurate_predictions = self._count_num_accurate_classification(actual, predicted)
        return float(num_accurate_predictions) / N

    def recall(self, actual, predicted, N=5):
        '''
            Calculates recall defined by:
                | Intersection(C-adopted, C-NRecommended) | / | C-adopted |

            Inputs:
                actual (List): Actual items that were liked
                predicted (List): Prediction of actual items
                N (Int): number of items to be recommended
            Output:
                Float: Number b/w 0 to 1 corresponding to the recall
                       of the recommendation list
        '''

        num_accurate_predictions = self._count_num_accurate_classification(test_set_positive_items, recommendations)
        return float(num_accurate_predictions) / len(test_set_positive_items)

    def average_precision(self, u, Y, N):
        '''
            Average Precision at N is defined as:
                sum(Precision@k * Rel(K), k, 1, N) / min(N, |C-adopted|)
            Inputs:
                u (Int): Represents user u
                Y (AbstractInputSource): user data class that encapsulates
                                         training, test, and full data
                N (Int): number of items to be recommended
            Output:
                Float; Number b/w 0 to 1 corresponding to the Average precision
                       of the recommendation list

        '''
        min_u_adopted = min(N, len(test_set_positive_items))
        actual = Y.get_test_data()
        predicted = self.recommend(u, Y, N=N, test_set=True)
        return (sum([self.precision(actual, predicted, N=k) if predicted[k] in actual else 0.0 for k in xrange(1, N+1)]) /
                float(min_N_adopted))

    def mean_average_precision(self, N):
        '''
            Average of AP for all users

            Inputs:
                N (Int): Number of items to be recommended
            Output:
                Float: Number b/w 0 to 1 corresponding to the Mean Average Precision@N
        '''
        return (sum([self.average_precision(u, user, N) for u, user in enumerate(self.input_data.get_data())]) /
               float(len(self.users)))

    def _init_variables(self, n_hidden=50):
        '''

            Initialize Tensorflow Variables
            Inputs:

            Outputs:

        '''
        # TODO-Look into slicing Big Variables and see if they are updatable.. Check performance
        print 'Creating Tensorflow Variables...'
        W = utils.xavier_init(self.n_items, n_hidden)
        W_prime = utils.xavier_init(self.n_items, n_hidden)
        all_weights = dict()
        all_weights['W'] = [tf.Variable(W[i]) for i in xrange(self.n_items)]
        all_weights['W_prime'] = [tf.Variable(W_prime[i]) for i in xrange(self.n_items)]
        all_weights['b'] = tf.Variable(tf.zeros([1, n_hidden], dtype=tf.float32))
        all_weights['b_prime'] = [tf.Variable(tf.constant(0, dtype=tf.float32), dtype=tf.float32) for _ in xrange(self.n_items)]
        all_weights['V'] = [tf.Variable(tf.zeros([n_hidden])) for _ in xrange(self.n_users)]
        print 'Done.'
        return all_weights


text:
# NVIDIA

import numpy as np
import pickle
import unittest

from parameterized import parameterized, parameterized_class

import test_data_path
import unit_test_utils
import global_vars as gv

from modeling import BertConfig, BertLayer

import torch
from torch import nn

class TestEncoders(unittest.TestCase):
  @classmethod
  def setUpClass(self):
    self.device = 'cuda:0'
    self.data_path = test_data_path.get_path() + '/gold/batch_1/'
    self.config_path = '/workspace/phase1/bert_config.json'
    self.bert_config = BertConfig(vocab_size_or_config_json_file=self.config_path)

    # Additional config options
    self.bert_config.fused_mha = False
    self.bert_config.fused_gelu_bias = True
    self.bert_config.dense_seq_output = True
    self.bert_config.unpad = False
    self.bert_config.pad = False
    self.bert_config.fuse_qkv = True
    self.bert_config.fuse_scale = True
    self.bert_config.fuse_mask = True
    self.bert_config.fuse_dropout = False
    self.bert_config.apex_softmax = True
    self.bert_config.enable_stream = False
    self.bert_config.fused_dropout_add = False

    self.n_encoder_layers = 24
    self.encoder_layers = [BertLayer(self.bert_config) for x in range(self.n_encoder_layers)]

    # Disable things like dropout
    for layer in self.encoder_layers:
      layer.eval()

    self.tf_input_mask_path = self.data_path + 'input_tensors/input_mask.np'
    self.tf_input_mask = torch.from_numpy(np.load(self.tf_input_mask_path)).to(self.device)

    tf_tensors_path = self.data_path + 'intermediate_tensors/all_intermediate_tensors.pkl'
    if not gv.tf_tensors:
      with open(tf_tensors_path, mode='rb') as f:
        gv.tf_tensors = pickle.load(f)

    # Tensor keys needed for retrieval
    self.encoder_key_template = {
      'bert/encoder/layer_%/attention/self/from_tensor' : None,
      'bert/encoder/layer_%/attention/self/query_output' : None,
      'bert/encoder/layer_%/attention/self/key_output' : None,
      'bert/encoder/layer_%/attention/self/value_output' : None,
      'bert/encoder/layer_%/attention/self/attention_score_output' : None,
      'bert/encoder/layer_%/attention/self/attention_score_scaled_output' : None,
      'bert/encoder/layer_%/attention/self/attention_score_additive_output' : None,
      'bert/encoder/layer_%/attention/self/attention_probs_output' : None,
      'bert/encoder/layer_%/attention/self/context_layer' : None
    }

    # Check to see if this is the first setup to load the TF weights, otherwise just use the already loaded data
    if not gv.tf_weights:
      with open(self.data_path + 'tf_checkpoint.pkl', mode='rb') as f:
        gv.tf_weights = pickle.load(f)

    self.tf_encoder_tensors = []

    for layer_idx in range(self.n_encoder_layers):
      self.tf_encoder_tensors.append({key.replace('%', str(layer_idx)):torch.from_numpy(gv.tf_tensors[key.replace('%', str(layer_idx))]).to(self.device) for (key,value) in self.encoder_key_template.items()})

    # Weight keys needed for retrieval
    for layer_idx, encoder_layer in enumerate(self.encoder_layers):
      state_dict = encoder_layer.state_dict()
      pyt_tf_map = unit_test_utils.pyt_tf_mapping(state_dict, 'bert/encoder/layer_' + str(layer_idx) + '/')
      for key in pyt_tf_map:
        mapped_key = pyt_tf_map[key]
        if 'dense' in mapped_key or 'query' in mapped_key or 'key' in mapped_key or 'value' in mapped_key:
          state_dict[key] = torch.from_numpy(np.transpose(gv.tf_weights[mapped_key]))
        else:
          state_dict[key] = torch.from_numpy(gv.tf_weights[mapped_key])
      self.encoder_layers[layer_idx].load_state_dict(state_dict)

      self.encoder_layers[layer_idx] = encoder_layer.to(self.device)

    self.tolerance = 1e-5

  def test_encoders_initialized_properly(self):
    test_result = True

    for layer_idx, encoder_layer in enumerate(self.encoder_layers):
      state_dict = encoder_layer.state_dict()
      pyt_tf_map = unit_test_utils.pyt_tf_mapping(state_dict, 'bert/encoder/layer_' + str(layer_idx) + '/')
      for key in pyt_tf_map:
        mapped_key = pyt_tf_map[key]
        if 'dense' in mapped_key or 'query' in mapped_key or 'key' in mapped_key or 'value' in mapped_key:
          new_result = torch.equal(state_dict[key], torch.from_numpy(np.transpose(gv.tf_weights[mapped_key])).to(self.device))
          test_result = test_result and new_result
          if not new_result:
            print("Failure:", mapped_key)
        else:
          new_result = torch.equal(state_dict[key], torch.from_numpy(gv.tf_weights[mapped_key]).to(self.device))
          test_result = test_result and new_result
          if not new_result:
            print("Failure:", mapped_key)

    self.assertTrue(test_result)
  @parameterized.expand([(x,) for x in range(24)])
  def test_encoder_attention_query(self, layer_idx=0):
    pyt_result = self.encoder_layers[layer_idx].attention.self.query(self.tf_encoder_tensors[layer_idx]['bert/encoder/layer_' + str(layer_idx) + '/attention/self/from_tensor'])
    pyt_result = pyt_result.detach().cpu().numpy()

    counts, bins_absolute, bins_relative, max_idx = unit_test_utils.max_abs_diff_binning(pyt_result, self.tf_encoder_tensors[layer_idx]['bert/encoder/layer_' + str(layer_idx) + '/attention/self/query_output'].view(1,512,1024).detach().cpu().numpy())

    self.assertTrue(bins_absolute[-1] < self.tolerance)

  @parameterized.expand([(x,) for x in range(24)])
  def test_encoder_attention_key(self, layer_idx=0):
    pyt_result = self.encoder_layers[layer_idx].attention.self.key(self.tf_encoder_tensors[layer_idx]['bert/encoder/layer_' + str(layer_idx) + '/attention/self/from_tensor'])
    pyt_result = pyt_result.detach().cpu().numpy()

    counts, bins_absolute, bins_relative, max_idx = unit_test_utils.max_abs_diff_binning(pyt_result, self.tf_encoder_tensors[layer_idx]['bert/encoder/layer_' + str(layer_idx) + '/attention/self/key_output'].view(1,512,1024).detach().cpu().numpy())

    self.assertTrue(bins_absolute[-1] < self.tolerance)

  @parameterized.expand([(x,) for x in range(24)])
  def test_encoder_attention_value(self, layer_idx=0):
    pyt_result = self.encoder_layers[layer_idx].attention.self.value(self.tf_encoder_tensors[layer_idx]['bert/encoder/layer_' + str(layer_idx) + '/attention/self/from_tensor'])
    pyt_result = pyt_result.detach().cpu().numpy()

    counts, bins_absolute, bins_relative, max_idx = unit_test_utils.max_abs_diff_binning(pyt_result, self.tf_encoder_tensors[layer_idx]['bert/encoder/layer_' + str(layer_idx) + '/attention/self/value_output'].view(1,512,1024).detach().cpu().numpy())

    self.assertTrue(bins_absolute[-1] < self.tolerance)

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


text:
# Copyright 2019 The dm_control Authors.
#
# 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.
# ============================================================================

"""Tests for locomotion.arenas.floors."""

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

from absl.testing import absltest
from dm_control import mjcf
from dm_control.locomotion.arenas import floors
import numpy as np

class FloorsTest(absltest.TestCase):

  def test_can_compile_mjcf(self):
    arena = floors.Floor()
    mjcf.Physics.from_mjcf_model(arena.mjcf_model)

  def test_size(self):
    floor_size = (12.9, 27.1)
    arena = floors.Floor(size=floor_size)
    self.assertEqual(tuple(arena.ground_geoms[0].size[:2]), floor_size)

  def test_top_camera(self):
    floor_width, floor_height = 12.9, 27.1
    arena = floors.Floor(size=[floor_width, floor_height])

    self.assertGreater(arena._top_camera_y_padding_factor, 1)
    np.testing.assert_array_equal(arena._top_camera.quat, (1, 0, 0, 0))

    expected_camera_y = floor_height * arena._top_camera_y_padding_factor
    np.testing.assert_allclose(
        np.tan(np.deg2rad(arena._top_camera.fovy / 2)),
        expected_camera_y / arena._top_camera.pos[2])

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


text:
# encoding: utf-8
# ******************************************************
# requirement:python3
# Author: chyb
# Last modified: 20181023 14:00
# Email:chyb3.14@gmail.com
# Filename:train.py
# Description:
# ******************************************************
import numpy as np
import math
from input_data import input_data
import matplotlib.pyplot as plt

import time
from tqdm import trange
from random import random, randint

class BPNetwork(object):
    def __init__(self,feature_length,input_numbers,hiden_numbers,output_kinds,base_lr1,base_lr2,base_lr3):
        self.feature_length=feature_length
        self.input_numbers=input_numbers
        self.hiden_numbers=hiden_numbers
        self.output_kinds=output_kinds
        self.hiden_lr=base_lr1
        self.out_lr=base_lr2
        self.input_lr=base_lr3
        """
        the initialization for weights is important
        """
        self.w1=0.5*(np.random.random((self.feature_length, self.input_numbers))-0.5)
        self.w2=0.5*(np.random.random((self.input_numbers, self.hiden_numbers))-0.5)
        self.w3=0.5*(np.random.random((self.hiden_numbers, self.output_kinds))-0.5)

        self.input_offset=np.zeros(self.input_numbers)
        self.hiden_offset=np.zeros(self.hiden_numbers)
        self.output_offset=np.zeros(self.output_kinds)
    def sigmoid(self,x):
        output=[]
        for i in x:
            output.append(1/(1+math.exp(-i)))
        output=np.array(output)
        return output
    def forward(self,input_feature):
        self.input_feature=input_feature
        self.input_val=np.dot(self.input_feature,self.w1)+self.input_offset
        self.input_out=self.sigmoid(self.input_val)
        self.hiden_val=np.dot(self.input_out,self.w2)+self.hiden_offset
        self.hiden_out=self.sigmoid(self.hiden_val)
        self.out_val=np.dot(self.hiden_val,self.w3)+self.output_offset
        self.out_out=self.sigmoid(self.out_val)

    def backward(self,label):
        self.erro=label-self.out_out
        delta_out=self.erro*self.out_out*(1-self.out_out)
        delta_hiden=self.hiden_out*(1-self.hiden_out)*np.dot(self.w3,delta_out)
        delta_input=self.input_out*(1-self.input_out)*np.dot(self.w2,delta_hiden)
        for i in range(0,self.output_kinds):
            self.w3[:,i]+=self.hiden_lr*delta_out[i]*self.hiden_out
        for i in range(0,self.hiden_numbers):
            self.w2[:,i]+=self.out_lr*delta_hiden[i]*self.input_out
        for i in range(0,self.input_numbers):
            self.w1[:,i]+=self.input_lr*delta_input[i]*self.input_feature
        self.output_offset+=self.out_lr*delta_out
        self.hiden_offset+=self.hiden_lr*delta_hiden
        self.input_offset+=self.input_lr*delta_input
    def reduce_lr(self,lr1,lr2,lr3):
        """
        reduce learning rate during training
        """
        self.hiden_lr=lr1
        self.out_lr=lr2
        self.input_lr=lr3
if __name__=="__main__":
    print ("*******Data loading*******")
    sample, label, test_s, test_l = input_data()
    sample = np.array(sample,dtype='float')
    sample=(sample)/256.0
    test_s = np.array(test_s,dtype='float')
    test_s =(test_s)/256.0

    samp_num = len(sample)
    inp_num = len(sample[0])
    out_num = 10
    hid_num = 15
    BP=BPNetwork(inp_num,30,hid_num,out_num,0.1,0.1,0.1)
    X=[]
    Loss=[]
    Loss_temp=[]
    loss=0
    logtext = open("./saved/train_test_log.txt", "w")
    print("*******Model trainng*******")
    logtext.write("*******Model trainng*******" + "\n")
    for step in range(0,1):
        if step==1:
            BP.reduce_lr(0.01,0.01,0.01)
        elif step==2:
            BP.reduce_lr(0.001,0.001,0.001)
        with trange(samp_num) as t:
            # for i in range(0,samp_num):
            for i in t:
                # 设置进度条左边显示的信息
                t.set_description("step %i" % i)

                train_label = np.zeros(out_num)
                train_label[label[i]] = 1

                BP.forward(sample[i])
                BP.backward(np.array(train_label))
                error=BP.erro
                for j in range(0,len(error)):
                        loss=loss+abs(error[j])
                Loss_temp.append(loss)
                if i%100==99 and step==0:
                    Loss.append(sum(Loss_temp)/len(Loss_temp))
                    Loss_temp=[]
                    X.append(i+step*60000)
                    # 设置进度条右边显示的信息
                    t.set_postfix(loss=loss, step=i)
                    logtext.write("step:" + str(i) + " " + "loss:" + str(loss) + "\n")
                    loss=0

    print ("*******Model Testing*******")
    logtext.write("*******Model Testing*******" + "\n")
    correct=0
    for i in range(0,len(test_s)):
        BP.forward(test_s[i])
        result=BP.out_out
        if np.argmax(result) == test_l[i]:
            correct+=1
    print ("Predicted Correct:", correct)
    print ("Predicted Correct Rate:", float(correct)/float(len(test_s)))
    logtext.write("Predicted Correct:" + str(correct) + "\n"+ "Predicted Correct Rate:"+str(float(correct)/float(len(test_s))))
    logtext.close()
    plt.figure()
    plt.plot(X,Loss)
    plt.savefig("./saved/train_loss.jpg")
    plt.show()


text:
import numpy as np
from layer import Layer

class DenseLayer(Layer):
    """ Add dense layer.

    Parameters
    ----------
    nodes_prev : int
        number of hidden units in previous layer
    nodes_curr : int
        number of hidden units in current layer
    init : obj
        how to initialize weights. Methods are found in initialize.py
    activation : obj
        activation function. Functions are found in activation.py
    optimizer : obj
        optimizer function. Methods are found in optimizer.py
    bias : bool
        bias on previous layer on (True) / off (False)
    """
    def __init__(self, nodes_prev, nodes_curr, init=None, activation=None, optimizer=None, bias=None):
        '''
        if init is None:
            super(Layer, self).__init__()
        else:
            super(Layer, self).__init__(init)
        if activation is None:
            super(Layer, self).__init__()
        else:
            super(Layer, self).__init__(activation)
        if optimizer is None:
            super(Layer, self).__init__()
        else:
            super(Layer, self).__init__(optimizer)
        if bias is None:
            super(Layer, self).__init__()
        else:
            super(Layer, self).__init__(bias)
        '''
        #super().__init__()
        self.bias = bias
        if self.bias:
            nodes_prev += 1   # Adding bias node
        self.weight = init(size=(nodes_prev, nodes_curr))
        self.activation = activation
        self.optimizer = optimizer

    def forward(self, input_layer):
        """Forward propagation. Multiply input_layer with weight matrix.

        Parameters
        ----------
        input_layer : ndarray
            Output from previous layer.
        """
        return input_layer.dot(self.weight)

    def __call__(self, input_layer):
        """Activation of output from forward propagation.

        Parameters
        ----------
        input_layer : ndarray
            Output from previous layer.
        """
        if self.bias:
            input_layer = np.c_[input_layer, np.ones((len(input_layer),1))]
        self.input_layer = input_layer
        z = self.forward(input_layer)
        self.output_layer = self.activation(z)
        return self.output_layer

    def backward(self, dcost):
        """ Backward propagation.

        Parameters
        ----------
        dcost : ndarray
            derivative of cost function with respect to the activation array
        """
        df = self.activation.derivate()
        self.delta = dcost * df
        dcost_new = np.einsum('ik,jk->ij',self.delta,self.weight)
        if self.bias:
            return dcost_new[:,:-1]
        else:
            return dcost_new

    def update_weights(self, step):
        gradient = self.input_layer.T.dot(self.delta)
        self.weight -= self.optimizer(step+1, gradient)
        return self.weight

max_stars_repo_name: evenmn/tensornet

max_issues_repo_name: evenmn/tensornet

max_forks_repo_name: evenmn/tensornet


text:
#! /usr/bin/env python
import sys
import argparse
import screed
import math
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import os.path
from scipy.stats import t

def getstats(infile1, infile2):

    sample1 = []
    sample2 = []
    count = 0
    for line in open(infile1):
        if count == 0:
            count+=1
            continue
        dr1, _, _, _, SE1, _, _, _, n1 = line.split('\t')
    count = 0
    for line in open(infile2):
        if count == 0:
             count+=1
             continue
        dr2, _, _, _, SE2, _, _, _, n2 = line.split('\t')
    df = ( float(n1) + float(n2) )-4
    tscore = (float(dr1) - float(dr2)) / np.sqrt(np.square(float(SE1)) +np.square(float(SE2)) )
    pvalue =  2 *  (t.sf(abs(tscore),df= df))
    pvalue = np.round(pvalue, decimals=4)
    print(infile1,'\t', infile2)
    print("tscore",'\t', "pvalue")
    print(tscore,'\t', pvalue)

def main():
    parser = argparse.ArgumentParser()
    parser.add_argument('infile1', default=False)
    parser.add_argument('infile2', default=False)
    args = parser.parse_args()
    getstats(args.infile1, args.infile2)

if __name__ == '__main__':
    main()

max_stars_repo_name: SherineAwad/ribofilio

max_issues_repo_name: SherineAwad/ribofilio

max_forks_repo_name: SherineAwad/ribofilio

reason: import numpy,from scipy