text:
# -*- coding: utf-8 -*-
"""
Created on Tue Aug 16 13:17:59 2016

@author: as624
"""

import os
from time import strftime
from math import sqrt
import scipy.ndimage as ndi
from skimage.measure import regionprops
import numpy as np
from mmhelper.utility import logger
from mmhelper.measurements_class import BacteriaData

def find_input_filename(fname, out=None, batch=False,
                        im_area=None, debug=False):
    """
    Takes the filename input and splits it into simply the filename
    which is then used to name the csv later
    Also creates output directories.

    Parameters
    ------
    fname : string
        A string containing the input filename
    out : string, optional
        A string containing an optional output file_name (default : None)
    batch : Boolean, optional
        Whether or not the analysis is being run in batch
        mode (default : False)
    im_area : string, optional
        The area of the chip being analysed (default : None)
    debug : Boolean, optional
        If true then the filename is stamped with 'debug' instead of the time

    Returns
    ------
    new_dir : string
        The full path for the main output directory for the results
    file_name : string
        The name of the input file
    image_dir : string
        The full path to the the directory for output images

    """
    dir_name = os.path.split(os.path.realpath(fname))[0]

    timestamp = "debug" if debug else strftime("%Y%m%d_%H%M%S")

    if out is None:
        file_name = os.path.splitext(os.path.basename(fname))[0]
        if batch is True:
            new_folder = "Area" + im_area + "_results_" + timestamp
        else:
            new_folder = file_name + "_results_" + timestamp
        new_dir = os.path.join(dir_name, new_folder)
    else:
        file_name = out
        new_folder = file_name
        new_dir = os.path.join(dir_name, new_folder)

    img = "output_images"
    image_dir = os.path.join(new_dir, img)
    if not os.path.isdir(new_dir):
        os.makedirs(new_dir)
    elif not debug:
        logger.critical("Results folder already exists [{}]".format(new_dir))
    if not os.path.isdir(image_dir):
        os.makedirs(image_dir)
    elif not debug:
        logger.critical("Image folder already exists [{}]".format(image_dir))

    return new_dir, file_name, image_dir

def get_measurements(data, fluo_data, fullwellimages,
                     allwellcoords, allbacteria, bacteria_lineage):
    """Takes the original image as well as the detected outputs and
    creates an instance of the measurements class containing the
    information for each bacteria

    Parameters
    ------
    data : ndarray
        Array containing all the original data (2D if just one frame)
    fluo_data : ndarray
        Array containing all the original fluorescent data (2D if
        just one frame)
    fullwellimages : list
        A list where each entry is a labelled image showing the
        detected wells
    allwellcoords : list
        A list where each entry is a dictionary containing wells and
        their coordinates
    allbacteria : list
        A list where each entry is a dictionary containing labelled
        images showing the detected bacteria
    bacteria_lineage : dictionary
        Keys are the 'actual' bacteria label and the value
        a reference to their lineage

    Returns
    ------
    measurements : Custom class instance
        Its attribute "bacteria" is a dictionary containing information on
        each individual bacteria
    fluorescence_backgrounds : list
        A list of the background fluorescence value for each timepoint
    """
    measurements = BacteriaData()
    fluorescence_backgrounds = {}
    for tpoint, (fullwells, bacteria, coords) in enumerate(zip(
            fullwellimages, allbacteria, allwellcoords)):
        # Create bacteria labelled image
        bacteriaim = np.zeros_like(fullwells)
        for welllabel in coords:
            bacteriaim[coords[welllabel]] = bacteria[welllabel]
        for region in regionprops(bacteriaim):
            well_label = fullwells[tuple(region.coords[0])]
            measurements.add_bac_data(
                region.label, bacteria_lineage, region, tpoint,
                well_label=well_label)
            if fluo_data is None:
                continue
            timepoint_fluo = [f[tpoint] for f in fluo_data]
            fluo_values = [
                fluorescence_background(
                    fullwells,
                    bacteriaim,
                    fluo_im) for fluo_im in timepoint_fluo]
            fluorescence_backgrounds[tpoint] = [values[0]
                                                for values in fluo_values]
            measurements.measure_fluo(
                region, timepoint_fluo, fluo_values, tpoint)
    measurements.compile_results(max_tpoint=len(data))
    return measurements, fluorescence_backgrounds

def count_bacteria_in_wells(bacteria_labels):
    """
    For each well, return the number of bacteria it contains;
    NOTE: This works on lists of bacteria labels - use count_bacteria
          if working with full size labels of wells and bacteria

    Parameters:
    ------
    bacteria_labels : dictionary
        A dictionary containing labelled images showing the detected bacteria

    Returns
    ------
    numbac2 : list
        List containing the number of bacteria detected in each well
    """
    numbac2 = []
    for bac in bacteria_labels.values():
        numbac = ndi.label(bac)[1]
        numbac2.append(numbac)
    return numbac2
    # return [ndi.label(w>0)[1] for w in bacteria_labels]

def fluorescence_measurements(
    region,
    fluo_image,
    bkg,
):
    """Takes a region props area from a brightfield image and determines
    the mean fluorescence in a matching fluorescent image.
    Returns the average fluorescence, background deducted fluorescence
    (if its within 2 Standard errors of the background this is 0),
    and intergrated fluorescence

    Parameters
    ------
    region : Region props area
        A region props area with coordinates of the bacteria of interest
    fluo_image : ndarray (2D)
        The matching fluorescent image for the detected bacteria
    bkg : tuple
        tuple in the format (background fluorescence, background SEM) for the
        respective image
    """
    bkground = bkg[0]
    bg_2sem = bkg[0] + (bkg[1])
    fluo = np.mean(fluo_image[tuple(np.array(region.coords).T)])
    background_fluo = [0 if (fluo - bg_2sem) < 0 else (fluo - bkground)][0]
    integrated_fluorescence = background_fluo * region.area
    return fluo, background_fluo, integrated_fluorescence

def fluorescence_background(wells, bacteria, fluo_image):
    """
    Removes bacteria segmentation from well segmentation, then calculates
    the fluorescence_background as the mean of the remaining segmentation
    (in the fluorescece channel)
    """
    background_wells = (wells > 0) ^ (bacteria > 0)

    backg = fluo_image[background_wells]
    bkground = np.mean(backg)
    bg_stdev = np.std(backg)
    bg_sem = bg_stdev / (sqrt(len(backg)))
    return bkground, bg_sem

max_stars_repo_name: jmetz/momanalysis

max_issues_repo_name: jmetz/momanalysis

max_forks_repo_name: jmetz/momanalysis

reason: import numpy,import scipy


text:
import kivy
kivy.require('1.10.1')
import requests as req

from kivy.app import App
from kivy.uix.gridlayout import GridLayout
from kivy.uix.boxlayout import BoxLayout
from kivy.uix.popup import Popup
from kivy.uix.label import Label
from kivy.properties import StringProperty

import cv2
import numpy as np
import sys

import time

url = "https://icanhazdadjoke.com/"
choice =  True
search_q = ""
ch = ""

class YLYL(GridLayout):
	def newJoke(self):
		try:
			res = req.get(
				url,
				headers ={"Accept" : "application/json"}
			).json()
			self.display.text = str(res["joke"])
		except Exception:
			self.display.text = "Error finding the joke"

	def checkLaughter (self):
		content = ConfirmPopup(text = 'Did you laugh?')
		content.bind(on_answer = self._on_answer)
		self.popup = Popup(title = "The big question",
			content = content,
			size_hint = (None, None),
			size = (480, 480))

		self.popup.open()

	def _on_answer(self, instance, answer):
		print ("User answer: ", answer)
		if answer == "yes" :
			self.display.text = "Oh No you lose!"
		elif answer == "no" :
			self.display.text = "Congrats you win!"
		else:
			print("someerror")
		self.popup.dismiss()

class ConfirmPopup(GridLayout):
	text = StringProperty()
	def __init__(self, **kwargs):
		self.register_event_type('on_answer')
		super(ConfirmPopup,self).__init__(**kwargs)

	def on_answer(self, *args):
		pass

class CameraClick(BoxLayout):
    def capture(self):
        '''
        Function to capture the images and give them the names
        according to their captured time and date.
        '''
        camera = self.ids['camera']
        timestr = time.strftime("%Y%m%d_%H%M%S")
        camera.export_to_png(f"IMG_{timestr}.png")
        print("Captured")

class ylylApp(App):
	def build(self):
		return YLYL()

ylylApp().run()

ext: py

max_stars_repo_name: smiteshz/YLYLDadJokes

max_issues_repo_name: smiteshz/YLYLDadJokes

max_forks_repo_name: smiteshz/YLYLDadJokes

reason: import numpy


text:
from dolfin import has_petsc4py, parameters, as_backend_type
from PETScIO import arrayToVec

import petsc4py as petsc4py
import sys

petsc4py.init(sys.argv)
import petsc4py.PETSc as PETSc

import scipy.sparse as sp

def ParameterSetup(parameters, backend):
    parameters['linear_algebra_backend'] = backend

def Assemble(AA, bb=None):

    backend = AA.str(False)

    if AA.__str__().find('uBLAS') == -1:
        if AA.size(1) == AA.size(0):
            A = as_backend_type(AA).mat()

        else:
            A = as_backend_type(AA).vec()
        if bb:
            b = as_backend_type(bb).vec()
    elif AA.__str__().find('Eigen') == -1::
        if AA.size(1) == AA.size(0):
            As = AA.sparray()
            As.eliminate_zeros()
            row, col, value = As.indptr, As.indices, As.data
            A = PETSc.Mat().createAIJ(size=(AA.size(0), AA.size(1)), csr=(
                row.astype('int32'), col.astype('int32'), value))
        else:
            A = arrayToVec(AA.array())
        if bb:
            b = arrayToVec(bb.array())
    if bb:
        return A, b
    else:
        return A

def Scipy2PETSc(A):
    A = A.tocsr()
    return PETSc.Mat().createAIJ(size=A.shape, csr=(A.indptr, A.indices, A.data))

def PETSc2Scipy(A):
    row, col, value = A.getValuesCSR()
    return sp.csr_matrix((value, col, row), shape=A.size)

max_stars_repo_name: wathen/PhD

max_issues_repo_name: wathen/PhD

max_forks_repo_name: wathen/PhD

reason: import scipy


text:
# -*- coding: utf-8 -*-
# Author : Jin Kim
# e-mail : jinkim@seculayer.com
# Powered by Seculayer © 2020 Solution Development 2 Team, R&D Center.
#
import random
import time
import numpy as np
from hps.common.Common import Common
from hps.algorithms.HPOptimizationAbstract import HPOptimizationAbstract
from hps.algorithms.HPOptimizerUtil import HPOptimizerUtil

class GridSearch(HPOptimizationAbstract):

    def __init__(self, **kwargs):
        # inheritance init
        super(GridSearch, self).__init__(**kwargs)
        self._check_hpo_params()
        self.DUP_CHECK = False
        self._categorical_param_dict = dict()

        self.categorical_boundary()
        self.LOGGER = Common.LOGGER.getLogger()

    ##### implement methods
    def _check_hpo_params(self):
        self._n_pop = self._n_params
        self._n_steps = self._n_steps

    ## grid search main function

    def _generate(self, param_list, score_list):

        best_param_list = self._population(param_list)
        result_param_list = HPOptimizerUtil.remove_duplicated_params(self.unique_param_dict, best_param_list)

        num_result_params = len(result_param_list)

        ## leak
        if  num_result_params < self._n_pop:
            result_param_list += self._generate_categorical_param_dict_list(self._n_pop - num_result_params)
        ## over
        elif num_result_params > self._n_pop :
            random.shuffle(result_param_list)
            result_param_list = result_param_list[:self._n_pop]
        return result_param_list

    #############################################################################
    ### grid search private functions
    def _population(self, param_list):
        if len(param_list) == 0:
            return self._generate_categorical_param_dict_list(self._n_pop)
        else :
            return param_list

    def _generate_categorical_param_dict_list(self, n_pop):
        temp_list = list()

        for i in range(n_pop):
            temp_dict = self._generate_catecorical_param_dict()
            temp_list.append(temp_dict)
        return temp_list

    def categorical_boundary(self):
        for k in self._pbounds.keys():
            if k == "optimizer_fn" or k == "act_fn":
                self._categorical_param_dict[k] = self._pbounds[k]
            elif k=="hidden_units":
                hidden_list=list()
                for i in range(np.random.randint(1, self._pbounds[k][0])):
                    hidden_list.append(int(np.random.randint(1, self._pbounds[k][1])))
                self._categorical_param_dict[k] = hidden_list
            else:
                temp_list = list()
                for i in range (5):
                    temp_list.append(np.random.uniform(self._pbounds[k][0], self._pbounds[k][1]))
                self._categorical_param_dict[k] = temp_list

    def _generate_catecorical_param_dict(self):
        param_dict = dict()
        for k in self._categorical_param_dict.keys():
            param_dict[k] = np.random.choice(self._categorical_param_dict[k])

        tmp_hash = HPOptimizerUtil.param_dict_to_hash(param_dict)
        if not HPOptimizerUtil.is_duplicate(self.unique_param_dict, tmp_hash):
            self.unique_param_dict[tmp_hash] = {"param_dict": param_dict}
            return param_dict
        else:
            return self._generate_catecorical_param_dict()

max_stars_repo_name: GyunHyukLee/GoldMine

max_issues_repo_name: GyunHyukLee/GoldMine

max_forks_repo_name: GyunHyukLee/GoldMine


text:
import torch
import torch.nn as nn
import torch.nn.functional as F
import pickle
import numpy as np

from mmdet.core import multiclass_nms
from mmcv.runner import force_fp32
from .convfc_bbox_head import Shared2FCBBoxHead
from mmdet.models.builder import HEADS, build_loss

@HEADS.register_module()
class GSBBoxHeadWith0(Shared2FCBBoxHead):

    def __init__(self,
                 fc_out_channels=1024,
                 gs_config=None,
                 *args,
                 **kwargs):
        super(GSBBoxHeadWith0, self).__init__(
                fc_out_channels=fc_out_channels,
                *args, **kwargs)
        self.fc_cls = nn.Linear(self.cls_last_dim,
                                self.num_classes + gs_config.num_bins)

        # self.loss_bg = build_loss(gs_config.loss_bg)

        self.loss_bins = []
        for i in range(gs_config.num_bins):
            self.loss_bins.append(build_loss(gs_config.loss_bin))

        # self.label2binlabel = torch.load(gs_config.label2binlabel).cuda()
        # self.pred_slice = torch.load(gs_config.pred_slice).cuda()
        with open(gs_config.label2binlabel, 'rb') as f:
            self.label2binlabel = torch.from_numpy(
                pickle.load(f)
            ).cuda()
        with open(gs_config.pred_slice, 'rb') as f:
            self.pred_slice = torch.from_numpy(
                pickle.load(f)
            ).cuda()

        # TODO: update this ugly implementation. Save fg_split to a list and
        #  load groups by gs_config.num_bins
        with open(gs_config.fg_split, 'rb') as fin:
            fg_split = pickle.load(fin)

        self.fg_splits = []
        # self.fg_splits.append(torch.from_numpy(fg_split['(0, 10)']).cuda())
        # self.fg_splits.append(torch.from_numpy(fg_split['[10, 100)']).cuda())
        # self.fg_splits.append(torch.from_numpy(fg_split['[100, 1000)']).cuda())
        # self.fg_splits.append(torch.from_numpy(fg_split['[1000, ~)']).cuda())
        self.fg_splits.append(torch.from_numpy(fg_split['(0, 1e2)']).cuda())
        self.fg_splits.append(torch.from_numpy(fg_split['[1e2, 1e3)']).cuda())
        self.fg_splits.append(torch.from_numpy(fg_split['[1e3, 1e4)']).cuda())
        self.fg_splits.append(torch.from_numpy(fg_split['[1e4, ~)']).cuda())

        self.others_sample_ratio = gs_config.others_sample_ratio

    def _sample_others(self, label):

        # only works for non bg-fg bins

        fg = torch.where(label > 0, torch.ones_like(label),
                         torch.zeros_like(label))
        fg_idx = fg.nonzero(as_tuple=True)[0]
        fg_num = fg_idx.shape[0]
        if fg_num == 0:
            return torch.zeros_like(label)

        bg = 1 - fg
        bg_idx = bg.nonzero(as_tuple=True)[0]
        bg_num = bg_idx.shape[0]

        bg_sample_num = int(fg_num * self.others_sample_ratio)

        if bg_sample_num >= bg_num:
            weight = torch.ones_like(label)
        else:
            sample_idx = np.random.choice(bg_idx.cpu().numpy(),
                                          (bg_sample_num, ), replace=False)
            sample_idx = torch.from_numpy(sample_idx).cuda()
            fg[sample_idx] = 1
            weight = fg

        return weight

    def _remap_labels(self, labels_):
        # BAGS: mmdet now use self.num_clases represent bg
        labels = labels_.clone()
        bg = (labels == self.num_classes)
        labels += 1
        labels[bg] = 0

        num_bins = self.label2binlabel.shape[0]
        new_labels = []
        new_weights = []
        new_avg = []
        for i in range(num_bins):
            mapping = self.label2binlabel[i]
            new_bin_label = mapping[labels]

            if i < 1:
                weight = torch.ones_like(new_bin_label)
                # weight = torch.zeros_like(new_bin_label)
            else:
                weight = self._sample_others(new_bin_label)
            new_labels.append(new_bin_label)
            new_weights.append(weight)

            avg_factor = max(torch.sum(weight).float().item(), 1.)
            new_avg.append(avg_factor)

        return new_labels, new_weights, new_avg

    def _slice_preds(self, cls_score):

        new_preds = []

        num_bins = self.pred_slice.shape[0]
        for i in range(num_bins):
            start = self.pred_slice[i, 0]
            length = self.pred_slice[i, 1]
            sliced_pred = cls_score.narrow(1, start, length)
            new_preds.append(sliced_pred)

        return new_preds

    @force_fp32(apply_to=('cls_score', 'bbox_pred'))
    def loss(self,
             cls_score,
             bbox_pred,
             rois,
             labels,
             label_weights,
             bbox_targets,
             bbox_weights,
             reduction_override=None):
        losses = dict()

        if cls_score is not None:
            # Original label_weights is 1 for each roi.
            new_labels, new_weights, new_avgfactors = self._remap_labels(labels)
            new_preds = self._slice_preds(cls_score)

            num_bins = len(new_labels)
            for i in range(num_bins):
                losses['loss_cls_bin{}'.format(i)] = self.loss_bins[i](
                    new_preds[i],
                    new_labels[i],
                    new_weights[i],
                    avg_factor=new_avgfactors[i],
                    reduction_override=reduction_override
                )

        if bbox_pred is not None:
            bg_class_ind = self.num_classes
            # 0~self.num_classes-1 are FG, self.num_classes is BG
            pos_inds = (labels >= 0) & (labels < bg_class_ind)
            # do not perform bounding box regression for BG anymore.
            if pos_inds.any():
                if self.reg_decoded_bbox:
                    # When the regression loss (e.g. `IouLoss`,
                    # `GIouLoss`, `DIouLoss`) is applied directly on
                    # the decoded bounding boxes, it decodes the
                    # already encoded coordinates to absolute format.
                    bbox_pred = self.bbox_coder.decode(rois[:, 1:], bbox_pred)
                if self.reg_class_agnostic:
                    pos_bbox_pred = bbox_pred.view(
                        bbox_pred.size(0), 4)[pos_inds.type(torch.bool)]
                else:
                    pos_bbox_pred = bbox_pred.view(
                        bbox_pred.size(0), -1,
                        4)[pos_inds.type(torch.bool),
                           labels[pos_inds.type(torch.bool)]]
                losses['loss_bbox'] = self.loss_bbox(
                    pos_bbox_pred,
                    bbox_targets[pos_inds.type(torch.bool)],
                    bbox_weights[pos_inds.type(torch.bool)],
                    avg_factor=bbox_targets.size(0),
                    reduction_override=reduction_override)
            else:
                losses['loss_bbox'] = bbox_pred[pos_inds].sum()
        return losses

    @force_fp32(apply_to=('cls_score'))
    def _merge_score(self, cls_score):
        '''
        Do softmax in each bin. Decay the score of normal classes
        with the score of fg.
        From v1.
        '''
        batch_inference = False
        if cls_score.dim() == 3:
            cls_score.squeeze_(0)
            batch_inference = True

        num_proposals = cls_score.shape[0]

        new_preds = self._slice_preds(cls_score)    # split into each bins
        new_scores = [F.softmax(pred, dim=1) for pred in new_preds]

        bg_score = new_scores[0]
        fg_score = new_scores[1:]

        fg_merge = torch.zeros((num_proposals, self.num_classes+1)).cuda()
        merge = torch.zeros((num_proposals, self.num_classes+1)).cuda()

        # import pdb
        # pdb.set_trace()
        for i, split in enumerate(self.fg_splits):
            fg_merge[:, split] = fg_score[i][:, 1:]

        weight = bg_score.narrow(1, 1, 1)

        # Whether we should add this? Test
        fg_merge = weight * fg_merge

        # merge[:, 0] = bg_score[:, 0]
        # merge[:, 1:] = fg_merge[:, 1:]

        # BAGS: mmdet has moved bg to the last dim
        merge[:, :self.num_classes] = fg_merge[:, 1:]
        merge[:, self.num_classes] = bg_score[:, 0]

        # fg_idx = (bg_score[:, 1] > 0.5).nonzero(as_tuple=True)[0]
        # erge[fg_idx] = fg_merge[fg_idx]
        if batch_inference:
            return merge.unsqueeze(0)
        else:
            return merge

    @force_fp32(apply_to=('cls_score', 'bbox_pred'))
    # def get_det_bboxes(self,
    def get_bboxes(self,
                    rois,
                    cls_score,
                    bbox_pred,
                    img_shape,
                    scale_factor,
                    rescale=False,
                    cfg=None):
        # BAGS: F.softmax -> self._merge_score
        if isinstance(cls_score, list):
            cls_score = sum(cls_score) / float(len(cls_score))

        # if self.custom_cls_channels:
        #     scores = self.loss_cls.get_activation(cls_score)
        # else:
        #     scores = F.softmax(
        #         cls_score, dim=-1) if cls_score is not None else None
        assert not self.custom_cls_channels
        scores = self._merge_score(cls_score)

        batch_mode = True
        if rois.ndim == 2:
            # e.g. AugTest, Cascade R-CNN, HTC, SCNet...
            batch_mode = False

            # add batch dimension
            if scores is not None:
                scores = scores.unsqueeze(0)
            if bbox_pred is not None:
                bbox_pred = bbox_pred.unsqueeze(0)
            rois = rois.unsqueeze(0)

        if bbox_pred is not None:
            bboxes = self.bbox_coder.decode(
                rois[..., 1:], bbox_pred, max_shape=img_shape)
        else:
            bboxes = rois[..., 1:].clone()
            if img_shape is not None:
                max_shape = bboxes.new_tensor(img_shape)[..., :2]
                min_xy = bboxes.new_tensor(0)
                max_xy = torch.cat(
                    [max_shape] * 2, dim=-1).flip(-1).unsqueeze(-2)
                bboxes = torch.where(bboxes < min_xy, min_xy, bboxes)
                bboxes = torch.where(bboxes > max_xy, max_xy, bboxes)

        if rescale and bboxes.size(-2) > 0:
            if not isinstance(scale_factor, tuple):
                scale_factor = tuple([scale_factor])
            # B, 1, bboxes.size(-1)
            scale_factor = bboxes.new_tensor(scale_factor).unsqueeze(1).repeat(
                1, 1,
                bboxes.size(-1) // 4)
            bboxes /= scale_factor

        # Replace multiclass_nms with ONNX::NonMaxSuppression in deployment
        if torch.onnx.is_in_onnx_export():
            from mmdet.core.export import add_dummy_nms_for_onnx
            batch_size = scores.shape[0]
            # ignore background class
            scores = scores[..., :self.num_classes]
            labels = torch.arange(
                self.num_classes, dtype=torch.long).to(scores.device)
            labels = labels.view(1, 1, -1).expand_as(scores)
            labels = labels.reshape(batch_size, -1)
            scores = scores.reshape(batch_size, -1)
            bboxes = bboxes.reshape(batch_size, -1, 4)

            max_size = torch.max(img_shape)
            # Offset bboxes of each class so that bboxes of different labels
            #  do not overlap.
            offsets = (labels * max_size + 1).unsqueeze(2)
            bboxes_for_nms = bboxes + offsets
            max_output_boxes_per_class = cfg.nms.get(
                'max_output_boxes_per_class', cfg.max_per_img)
            iou_threshold = cfg.nms.get('iou_threshold', 0.5)
            score_threshold = cfg.score_thr
            nms_pre = cfg.get('deploy_nms_pre', -1)
            batch_dets, labels = add_dummy_nms_for_onnx(
                bboxes_for_nms,
                scores.unsqueeze(2),
                max_output_boxes_per_class,
                iou_threshold,
                score_threshold,
                pre_top_k=nms_pre,
                after_top_k=cfg.max_per_img,
                labels=labels)
            # Offset the bboxes back after dummy nms.
            offsets = (labels * max_size + 1).unsqueeze(2)
            # Indexing + inplace operation fails with dynamic shape in ONNX
            # original style: batch_dets[..., :4] -= offsets
            bboxes, scores = batch_dets[..., 0:4], batch_dets[..., 4:5]
            bboxes -= offsets
            batch_dets = torch.cat([bboxes, scores], dim=2)
            return batch_dets, labels
        det_bboxes = []
        det_labels = []
        for (bbox, score) in zip(bboxes, scores):
            if cfg is not None:
                det_bbox, det_label = multiclass_nms(bbox, score,
                                                     cfg.score_thr, cfg.nms,
                                                     cfg.max_per_img)
            else:
                det_bbox, det_label = bbox, score
            det_bboxes.append(det_bbox)
            det_labels.append(det_label)

        if not batch_mode:
            det_bboxes = det_bboxes[0]
            det_labels = det_labels[0]
        return det_bboxes, det_labels

    # BAGS: not used codes
    # def _remap_labels1(self, labels):

    #     num_bins = self.label2binlabel.shape[0]
    #     new_labels = []
    #     new_weights = []
    #     new_avg = []
    #     for i in range(num_bins):
    #         mapping = self.label2binlabel[i]
    #         new_bin_label = mapping[labels]

    #         weight = torch.ones_like(new_bin_label)

    #         new_labels.append(new_bin_label)
    #         new_weights.append(weight)

    #         avg_factor = max(torch.sum(weight).float().item(), 1.)
    #         new_avg.append(avg_factor)

    #     return new_labels, new_weights, new_avg

    # @force_fp32(apply_to=('cls_score'))
    # def _merge_score1(self, cls_score):
    #     '''
    #     Do softmax in each bin. Merge the scores directly.
    #     '''
    #     num_proposals = cls_score.shape[0]

    #     new_preds = self._slice_preds(cls_score)
    #     new_scores = [F.softmax(pred, dim=1) for pred in new_preds]

    #     bg_score = new_scores[0]
    #     fg_score = new_scores[1:]

    #     fg_merge = torch.zeros((num_proposals, 1231)).cuda()
    #     merge = torch.zeros((num_proposals, 1231)).cuda()

    #     for i, split in enumerate(self.fg_splits):
    #         fg_merge[:, split] = fg_score[i]

    #     merge[:, 0] = bg_score[:, 0]
    #     fg_idx = (bg_score[:,1] > 0.5).nonzero(as_tuple=True)[0]
    #     merge[fg_idx] = fg_merge[fg_idx]

    #     return merge

    # @force_fp32(apply_to=('cls_score'))
    # def _merge_score2(self, cls_score):
    #     '''
    #     Do softmax in each bin. Softmax again after merge.
    #     '''
    #     num_proposals = cls_score.shape[0]

    #     new_preds = self._slice_preds(cls_score)
    #     new_scores = [F.softmax(pred, dim=1) for pred in new_preds]

    #     bg_score = new_scores[0]
    #     fg_score = new_scores[1:]

    #     fg_merge = torch.zeros((num_proposals, 1231)).cuda()
    #     merge = torch.zeros((num_proposals, 1231)).cuda()

    #     for i, split in enumerate(self.fg_splits):
    #         fg_merge[:, split] = fg_score[i]

    #     merge[:, 0] = bg_score[:, 0]
    #     fg_idx = (bg_score[:,1] > 0.5).nonzero(as_tuple=True)[0]
    #     merge[fg_idx] = fg_merge[fg_idx]
    #     merge = F.softmax(merge)

    #     return merge

    # @force_fp32(apply_to=('cls_score'))
    # def _merge_score4(self, cls_score):
    #     '''
    #     Do softmax in each bin.
    #     Do softmax on merged fg classes.
    #     Decay the score of normal classes with the score of fg.
    #     From v2 and v3
    #     '''
    #     num_proposals = cls_score.shape[0]

    #     new_preds = self._slice_preds(cls_score)
    #     new_scores = [F.softmax(pred, dim=1) for pred in new_preds]

    #     bg_score = new_scores[0]
    #     fg_score = new_scores[1:]

    #     fg_merge = torch.zeros((num_proposals, 1231)).cuda()
    #     merge = torch.zeros((num_proposals, 1231)).cuda()

    #     for i, split in enumerate(self.fg_splits):
    #         fg_merge[:, split] = fg_score[i]

    #     fg_merge = F.softmax(fg_merge, dim=1)
    #     weight = bg_score.narrow(1, 1, 1)
    #     fg_merge = weight * fg_merge

    #     merge[:, 0] = bg_score[:, 0]
    #     merge[:, 1:] = fg_merge[:, 1:]
    #     # fg_idx = (bg_score[:, 1] > 0.5).nonzero(as_tuple=True)[0]
    #     # erge[fg_idx] = fg_merge[fg_idx]

    #     return merge

    # @force_fp32(apply_to=('cls_score'))
    # def _merge_score5(self, cls_score):
    #     '''
    #     Do softmax in each bin.
    #     Pick the bin with the max score for each box.
    #     '''
    #     num_proposals = cls_score.shape[0]

    #     new_preds = self._slice_preds(cls_score)
    #     new_scores = [F.softmax(pred, dim=1) for pred in new_preds]

    #     bg_score = new_scores[0]
    #     fg_score = new_scores[1:]
    #     max_scores = [s.max(dim=1, keepdim=True)[0] for s in fg_score]
    #     max_scores = torch.cat(max_scores, 1)
    #     max_idx = max_scores.argmax(dim=1)

    #     fg_merge = torch.zeros((num_proposals, 1231)).cuda()
    #     merge = torch.zeros((num_proposals, 1231)).cuda()

    #     for i, split in enumerate(self.fg_splits):
    #         tmp_merge = torch.zeros((num_proposals, 1231)).cuda()
    #         tmp_merge[:, split] = fg_score[i]
    #         roi_idx = torch.where(max_idx == i,
    #                               torch.ones_like(max_idx),
    #                               torch.zeros_like(max_idx)).nonzero(
    #             as_tuple=True)[0]
    #         fg_merge[roi_idx] = tmp_merge[roi_idx]

    #     merge[:, 0] = bg_score[:, 0]
    #     fg_idx = (bg_score[:, 1] > 0.5).nonzero(as_tuple=True)[0]
    #     merge[fg_idx] = fg_merge[fg_idx]

    #     return merge


text:
# -*- coding: utf-8 -*-
"""
% File:             med_to_tec.py
% File type:        Function
% Created on:       Jan 5, 2015
% Created by:       Kyra Swanson
% Last revised on:  Jun 5, 2018
% Last revised by:  Kyra Swanson

Converts MedPC file into list of [time,event] codes

Based on med_to_tec.m
Created by Marcelo S. Caetano

Requires Time-Event pairs to be saved as:
TIME.EVENT
For example, an event 1 at exactly 536 seconds from start will be recorded as
268000.001
Currently, the system clocks at 2ms, so time codes are divided by 500
to convert them to seconds.
Events are converted to integers.
The final product for this event is
[536,1]

"""
import numpy as np
import re

#convert time to seconds (BTIME at 2ms, 500/sec)
timeslice = 500
# timeslice = 200 (old files timeslice at 5ms = 200)

#convert a .med file into a readable format
def cleanLines(item):
    start = item.index(":") + 1
    item = item[start:]
    item = item.replace("\n","")
    item = item.split()
    return item

def breakLine(lines):
    Time_Event = []
    for item in lines:
        for el in item:
            temp = el.split(".")
            if temp[0] != "0":
                temp[0] = float(temp[0])/timeslice
                temp[1] = float(temp[1])
                Time_Event.append(temp)
    return np.asarray(Time_Event)

def convert(fileName,variable = "A"):
    fid = open(fileName) #open file
    lines = fid.readlines() #store lines in list
    fid.close()

    flag = variable + ":\n" #find variable flag (eg A:)
    start_index = next(i for i,x in enumerate(lines) if x == flag) + 1
    lines = lines[start_index:]
    stop_index = 0

    while not (re.search('[a-zA-Z]', lines[stop_index])): #stops at the next variable or comment
        lines[stop_index] = cleanLines(lines[stop_index]) #clean lines
        stop_index +=  1

    lines = lines[:stop_index]
    TEC = breakLine(lines) #break lines, store time,event in list Time_Event
    return TEC

#****************************************************************************
#   PULLING OUT DATA

#returns times for a given code
def findTimes(EV,code):
    return [val[0] for (i,val) in enumerate(EV) if val[1] == code]

#Returns the value for a given variable
def findVariable(fileName,variable):
    fid = open(fileName)
    lines = fid.readlines()
    fid.close()

    for i,j in enumerate(lines):
        if j[0] == variable:
            answer = lines[i]

    answer = answer[:-5]
    start = re.search("\d", answer).start()
    return int(answer[start:])

max_stars_repo_name: LaubachLab/MedParse

max_issues_repo_name: LaubachLab/MedParse

max_forks_repo_name: LaubachLab/MedParse


text:
'''
Segment
=======

This file will segment all the recordings into files containing only one
specific audio sample recording (sweep, noise, speech, etc).

Author: 2018 (c) Robin Scheibler
License: MIT License
'''

from __future__ import division, print_function

import argparse, os, re, sys
import numpy
from scipy.io import wavfile
from samplerate import resample

basedir = os.path.split(os.path.abspath(__file__))[0]
sys.path.append(basedir)
import ffmpeg_audio

# need that when no display is available
if 'DISPLAY' not in os.environ:
    import matplotlib as mpl
    mpl.use('Agg')
    need_agg = True
else:
    need_agg = False

from matplotlib import pyplot

# import reference file globally
path = os.path.split(os.path.dirname(os.path.realpath(__file__)))[0]
all_samples_file = path + '/all_samples.wav'
fs, all_samples = wavfile.read(all_samples_file)

# these are the lengths in samples of all the files concatenated during playback
# This list is printed by the stitch script when 'all_samples.wav' is created.
lengths= [48000, 144000, 1500, 144000, 1500, 144000, 1500, 151451,
        1500, 113357, 1500, 106598, 1500, 152987, 1500, 190463, 48000]

# these are the names of the individual segments in the order they were played
labels = ['silence', 'sweep_lin', 'sweep_exp', 'noise', 'fq_sample4',
        'fq_sample3', 'fq_sample2', 'fq_sample1', 'fq_sample0']

# filename pattern
pattern = re.compile('pyramic_spkr[012]_all_samples_\d+.[wav|tta]')

def compute_cost(audio, noise_mean, noise_est, noise_thresh):
    import numpy
    return numpy.mean(numpy.abs(audio - noise_mean[None,:]), axis=1) - noise_thresh * noise_est

def find_segments(audio, lengths, noise_mean, noise_est, noise_thresh):
    import numpy

    # end recursion
    if len(lengths) < 2:
        return []

    # find beginning of segment
    cost = compute_cost(audio[:lengths[0]+lengths[1],:], noise_mean, noise_est, noise_thresh)
    start_offset = int(numpy.min( numpy.where(cost > 0)[0] ))

    # find end of segment
    L = start_offset + lengths[1]
    cost = compute_cost(audio[L:L+lengths[2]*5,:], noise_mean, noise_est, noise_thresh)
    end_offset = 0
    q = cost > 0
    chunk = 50
    jump = 25
    while numpy.sum(q[end_offset:end_offset+chunk]) > 0:
        end_offset += jump

    end_offset += L

    # recursively search segments
    beyond = [x + end_offset for x in find_segments(audio[end_offset:,:], lengths[2:], noise_mean, noise_est, noise_thresh)]

    return [start_offset, end_offset] + beyond

def open_segment(filename, noise_thresh=3, off_lo=50, off_hi=50, plot=False):
    global lengths, labels

    import re, os
    import numpy
    from scipy.io import wavfile

    from matplotlib import pyplot

    # the file to segment
    rate, audio = wavfile.read(filename)

    # find offset here
    noise_mean = numpy.mean(audio[:500, :], axis=0)
    noise_est = numpy.mean(numpy.std(audio[:500,:] - noise_mean[None,:], axis=0))

    # find the boundary of first sweep, with length[1]
    boundaries = [0] + find_segments(audio, lengths, noise_mean, noise_est, noise_thresh)

    # now extract
    signals = [audio[:boundaries[1]-off_lo,:]]

    # list of views
    for i in range(1, len(boundaries)-1, 2):
        b_lo, b_hi = boundaries[i], boundaries[i+1]
        signals.append(audio[b_lo-off_lo:b_hi+off_hi])

    # make a dictionary
    d = dict(zip(labels, signals))
    if plot:
        for label, signal in d.items():
            pyplot.figure()
            pyplot.plot(signal[:,0])
            pyplot.title(label)
        pyplot.show()

    # add the rate
    d['rate'] = rate

    return d

def open_segment_rigid(filename, noise_thresh=3, off_lo=150, off_hi=150, plot=False):
    global lengths, labels, basedir

    import re, os, sys
    import numpy
    from scipy.io import wavfile

    from matplotlib import pyplot

    sys.path.append(basedir)
    import ffmpeg_audio

    # the file to segment
    if os.path.splitext(filename)[1] == '.wav':
        rate, audio = wavfile.read(filename)
    else:
        rate, audio = ffmpeg_audio.read(filename)

    # find offset here
    noise_mean = numpy.mean(audio[:500, :], axis=0)
    noise_est = numpy.mean(numpy.std(audio[:500,:] - noise_mean[None,:], axis=0))

    # recursive code (method 2)
    boundaries = [0] + find_segments(audio, lengths[:3], noise_mean, noise_est, noise_thresh)
    boundaries[2] = boundaries[1] + lengths[1]

    for i in range(2, len(lengths)-1):
        boundaries.append(boundaries[-1] + lengths[i])

    # now extract
    signals = [audio[:boundaries[1]-off_lo,:]]

    # list of views
    for i in range(1, len(boundaries)-1, 2):
        b_lo, b_hi = boundaries[i], boundaries[i+1]
        signals.append(audio[b_lo-off_lo:b_hi+off_hi])

    # make a dictionary
    d = dict(zip(labels, signals))
    if plot:
        for label, signal in d.items():
            pyplot.figure()
            pyplot.plot(signal[:,0])
            pyplot.title(label)
        pyplot.show()

    # add the rate
    d['rate'] = rate

    return d

def save_samples(filename):
    global pattern, output_dir, qc_images

    import re, os
    import numpy
    from scipy.io import wavfile
    from samplerate import resample

    from matplotlib import pyplot

    _, fname = os.path.split(filename)

    if not pattern.match(fname):
        return None

    spkr, angle = [int(i) for i in re.findall(r'\d+', fname)]
    signals = open_segment_rigid(filename, noise_thresh=3, off_lo=150, off_hi=150, plot=False)

    out_name = '{}_spkr{}_angle{}.{}'

    for name, signal in signals.items():

        if name == 'rate':
            continue

        if 'fq_sample' in name:
            # we resample the speech signals to 16kHz because this was the
            # orignal rate of the samples played
            # this saves some storage and load speed is faster
            signal = resample(signal, 16000. / signals['rate'], 'sinc_best')
            rate = 16000
        else:
            rate = signals['rate']

        # save in folder with name of the sample
        folder = os.path.join(output_dir, name)
        if not os.path.exists(folder):
            os.mkdir(folder)

        # format filename with angle and speaker location
        filename = out_name.format(name, spkr, (360 - angle) % 360, 'wav')

        signal_float = signal.astype(numpy.int16)
        wavfile.write(os.path.join(folder, filename), rate, signal_float)

        if qc_images:
            # save a spectrogram for later inspection

            folder = os.path.join(output_dir, 'qc_images', name)
            if not os.path.exists(folder):
                os.mkdir(folder)

            filename = out_name.format(name, spkr, (360 - angle) % 360, 'png')

            pyplot.specgram(signal[:,0].astype(numpy.float32) / (2**15+1), Fs=rate, NFFT=1024, noverlap=512)
            pyplot.savefig(os.path.join(folder, filename))
            pyplot.clf()

if __name__ == "__main__":

    parser = argparse.ArgumentParser(description='Segments files')
    parser.add_argument('path', type=str,
                        help='path to file to segment, or directory containing all files when using the option --all')
    parser.add_argument('-o', '--output', type=str, default='segmented',
                        help='output folder for the segemented recordings')
    parser.add_argument('-a', '--all', action='store_true',
                        help='segment all files in a directory')
    parser.add_argument('-q', '--qc', action='store_true',
                        help='produce spectrograms of all samples for quality control')

    args = parser.parse_args()

    # create necessary directories if not existing
    if not os.path.exists(args.output):
        os.mkdir(args.output)
    if args.qc and not os.path.exists(args.output + '/qc_images'):
        os.mkdir(args.output + '/qc_images')

    if args.all:

        if not os.path.isdir(args.path):
            raise ValueError('When segmenting multiple files, the argument should be a folder name')
            sys.exit()

        print('All files in {} will be segmented'.format(args.path))

        import ipyparallel as ip

        # Start the parallel processing
        c = ip.Client()
        NC = len(c.ids)
        print(NC,'workers on the job')

        # need that because no display is available
        if need_agg:
            c[:].map_sync(mpl.use, ['Agg'] * len(c))

        # push necessary global variables
        c[:].push(
                dict(
                    pattern=pattern,
                    output_dir=args.output,
                    qc_images=args.qc,
                    lengths=lengths,
                    labels=labels,
                    open_segment=open_segment,
                    open_segment_rigid=open_segment_rigid,
                    find_segments=find_segments,
                    compute_cost=compute_cost,
                    basedir=basedir,
                    ),
                block=True
                )
        print('pushing global variables to engines.')

        # all the files to process
        all_filenames = [fname for fname in os.listdir(args.path) if fname[0] != '.']

        # delegate to workers
        print('Dispatching {} files to engines.'.format(len(all_filenames)))
        try:
            c[:].map_sync(save_samples, [args.path + '/' + fname for fname in all_filenames])
        except:
            c.abort(block=True)

    else:
        open_segment_rigid(args.path, plot=True)
        output_dir = 'segmented'
        qc_images = True
        save_samples(args.path)

max_stars_repo_name: fakufaku/pyramic-dataset

max_issues_repo_name: fakufaku/pyramic-dataset

max_forks_repo_name: fakufaku/pyramic-dataset

reason: import numpy,from scipy


text:
#! /usr/bin/env python

import numpy as np
import os
import sys

names = []
weights = []

threshold = 0.01
best_mix = sys.stdin.readline()
ss = best_mix.split()
for i in range(len(ss)//4):
  name = ss[i*4+1]
  weight = float(ss[i*4+3])
  if weight >= threshold:
    names.append(os.path.basename(name).partition(".")[0])
    weights.append(weight)

weights = np.array(weights)
weights = weights / sum(weights)

for i in range(len(names)):
  print("%s 1 %f" % (names[i], weights[i]))

max_stars_repo_name: alumae/kaldi-estonian-recipe

max_issues_repo_name: alumae/kaldi-estonian-recipe

max_forks_repo_name: alumae/kaldi-estonian-recipe


text:
"""Tests for PhonopyAtoms."""
from io import StringIO

import numpy as np
import yaml

from phonopy.structure.atoms import PhonopyAtoms

symbols_SiO2 = ["Si"] * 2 + ["O"] * 4
symbols_AcO2 = ["Ac"] * 2 + ["O"] * 4
_lattice = [[4.65, 0, 0], [0, 4.75, 0], [0, 0, 3.25]]
_points = [
    [0.0, 0.0, 0.0],
    [0.5, 0.5, 0.5],
    [0.3, 0.3, 0.0],
    [0.7, 0.7, 0.0],
    [0.2, 0.8, 0.5],
    [0.8, 0.2, 0.5],
]
cell_SiO2 = PhonopyAtoms(cell=_lattice, scaled_positions=_points, symbols=symbols_SiO2)
cell_AcO2 = PhonopyAtoms(cell=_lattice, scaled_positions=_points, symbols=symbols_AcO2)

def test_SiO2():
    """Test of attributes by SiO2."""
    _test_cell(cell_SiO2, _lattice, _points, symbols_SiO2)

def test_AcO2():
    """Test of attributes by AcO2."""
    _test_cell(cell_AcO2, _lattice, _points, symbols_AcO2)

def _test_cell(cell, lattice, points, symbols):
    np.testing.assert_allclose(cell.cell, lattice, atol=1e-8)
    for s1, s2 in zip(cell.symbols, symbols):
        assert s1 == s2
    diff = cell.scaled_positions - points
    diff -= np.rint(diff)
    dist = np.linalg.norm(np.dot(diff, cell.cell), axis=1)
    np.testing.assert_allclose(dist, np.zeros(len(dist)), atol=1e-8)

def test_phonopy_atoms_SiO2():
    """Test of PhonopyAtoms __str__ by SiO2."""
    _test_phonopy_atoms(cell_SiO2)

def test_phonopy_atoms_AcO2():
    """Test of PhonopyAtoms __str__ by AcO2."""
    _test_phonopy_atoms(cell_AcO2)

def _test_phonopy_atoms(cell):
    with StringIO(str(PhonopyAtoms(atoms=cell))) as f:
        data = yaml.safe_load(f)
        np.testing.assert_allclose(cell.cell, data["lattice"], atol=1e-8)
        positions = []
        for atom, symbol in zip(data["points"], cell.symbols):
            positions.append(atom["coordinates"])
            assert atom["symbol"] == symbol
        diff = cell.scaled_positions - positions
        diff -= np.rint(diff)
        dist = np.linalg.norm(np.dot(diff, cell.cell), axis=1)
        np.testing.assert_allclose(dist, np.zeros(len(dist)), atol=1e-8)

max_stars_repo_name: phonopy/phonopy

max_issues_repo_name: phonopy/phonopy

max_forks_repo_name: phonopy/phonopy


text:
import numpy    as np
import warnings
from   numba    import njit

class SSECalculator():
    ###########################################################################
    #SSE Calculation
    ###########################################################################
    def reconstruct(self, free_params_array):
        sampled_index = self.sampled_index
        nominal_dct   = self.nominal
        p             = {}

        for scenario in self.init:
            nominal     = self.nominal[scenario]
            recon_array = self._reconstruct(nominal_dct[scenario], sampled_index, free_params_array)
            p[scenario] = recon_array

        return p

    @staticmethod
    @njit
    def _reconstruct(nominal, sampled_index, free_params_array):
        params                = nominal.copy()
        params[sampled_index] = free_params_array

        return params

    @staticmethod
    @njit
    def get_error(ym, yd, sd, idx):
        ym_ = ym if idx is None else ym[idx]
        return np.sum(np.abs(ym_-yd)**2)/sd

    @staticmethod
    def sort_params(params, states):
        try:
            return params.loc[states.index]
        except:
            raise ValueError('Parameters are missing one or more indices.')

    ###########################################################################
    #Instantiators
    ###########################################################################
    def __init__(self, model, dataset):
        free_params                   = model.optim_args['free_params']
        tspan, t_data, y_data, s_data = split_dataset(model, dataset)
        init                          = model._states

        nominal_vals  = model._params
        param_names   = model.get_param_names()
        sampled_index = [i for i, p in enumerate(param_names) if p in free_params]

        #Check
        param_check = set(param_names)
        if len(param_check.intersection(free_params)) != len(free_params):
            raise ValueError('Unexpected parameter(s) in free params.')

        #Assign
        self.tspan         = tspan
        self.t_data        = t_data
        self.y_data        = y_data
        self.s_data        = s_data
        self.init          = init
        self.model         = model
        self.nominal       = nominal_vals
        self.sampled_index = np.array(sampled_index)

    ###########################################################################
    #Integration
    ###########################################################################
    def __call__(self, free_params_array):
        SSE     = 0
        t_data  = self.t_data
        y_data  = self.y_data
        s_data  = self.s_data
        model   = self.model
        p       = self.reconstruct(free_params_array)

        for ir in model(p=p, tspan=self.tspan, overlap=False, include_events=False):
            scenario = ir.scenario

            for var, yd in t_data.get(scenario, {}).items():
                ym  = ir[var]
                idx = t_data[scenario][var]
                yd  = y_data[scenario][var]
                sd  = s_data[var]

                if type(ym) == dict:
                    ym = ym['y']

                SSE += self.get_error(ym, yd, sd, idx) if hasattr(ym, '__iter__') else self.get_error(ym, yd, sd, None)

        return SSE

    ###########################################################################
    #Miscellaneous
    ###########################################################################
    def __repr__(self):
        return f'{type(self).__name__}(model: {self.model.model_key})'

    def __str__(self):
        return self.__repr__()

###############################################################################
#Preprocessing
###############################################################################
def check_scenarios_match(model, dataset):

    d_scenarios = set([k[1] for k in dataset.keys()])
    m_scenarios = set(model.states.index)
    diff1       = d_scenarios.difference(m_scenarios)
    diff2       = m_scenarios.difference(d_scenarios)

    if diff1:
        msg = f'There are one or scenarios in dataset that are not in model.states for model {model.model_key}'
        warnings.warn(msg)

    if diff2:
        msg = f'There are one or scenarios in dataset that are missing from model.states for model {model.model_key}'
        warnings.warn(msg)

def split_dataset(model, dataset):
    mod_vars  = [*model.get_state_names(), *model.vrbs] if model.vrbs else model.get_state_names()
    model_exvs = model.exvs if model.exvs else []
    y_data     = {}
    t_data     = {}
    s_data     = {}
    y_set      = set()
    x_set      = set()
    t_set      = set()
    s_set      = set()
    t_len      = {}
    tpoints    = {}
    max_vals   = {}

    for key, value in dataset.items():
        if key[0] == 'Data':
            _, scenario, var = key

            #Update y_data
            y_data.setdefault(scenario, {})[var] = np.array(value)

            #Check if variable is a state or exv
            if var in mod_vars :
                y_set.add((scenario, var))
            elif var not in model_exvs:
                raise ValueError(f'exp_data contains an exv {var} which is not in the model.')

            #Track variable to determine if sd is provided
            x_set.add(var)

            #Track max value for sd calculations
            max_vals.setdefault(var, [])
            max_vals[var] += [*value]

        elif key[0] == 'Time':
            _, scenario, var = key

            #Create a list to store time points
            tpoints.setdefault(scenario, [])

            #Update t_data
            #Check length of tpoints[scenario]
            #Determine the start and stop indices
            t_data.setdefault(scenario, {})
            t_data[scenario][var] = len(tpoints[scenario]), len(tpoints[scenario]) + len(value)

            #Track to see if the corresponding y_data is availabe
            t_set.add((scenario, var))

            #Update tpoints AFTER update t_data
            tpoints[scenario] += [*value]
            t_len[(scenario, var)] = len(value)

        elif key[0] == 'Std':
            _, scenario, var = key

            #Update s_data
            s_data[var] = max(value, s_data.get(value, 0))

            #Track variables for which the user has provided sd
            s_set.add(var)

    #Check that all states have an associated time
    #Does not apply to exvs
    if len(y_set.intersection(t_set)) != len(y_set):
        raise ValueError()

    #Check that t and y arrays have the same length
    mismatched = [k for k, v in t_len.items() if len(y_data[k[0]][k[1]]) != v]
    if mismatched:
        raise ValueError('Mismatched data lengths.')

    #Make sd data if absent
    for state in x_set.difference(s_set):
        s_data[state] = np.percentile(max_vals[state], 75)/20

    #Make tspan for SSE calculation
    #Get the indices of tspan for each state/scenario
    tspan   = {}
    for scenario, tps in tpoints.items():
        tspan_, indices = np.unique(tps, return_inverse=True)
        if not len(tspan_):
            raise ValueError('tspan cannot have only 1 time point.')

        if tspan_[0] <= 0:
            delta           = 0
            tspan[scenario] = tspan_
        else:
            delta           = 1
            tspan[scenario] = np.array([0, *tspan_])

        for var in t_data[scenario]:
            start, stop = t_data[scenario][var]
            indices_    = indices[start:stop] + delta

            #Reassign with the indices
            t_data[scenario][var] = indices_

    #Check that scenarios match
    check_scenarios_match(model, dataset)

    return tspan, t_data, y_data, s_data

max_stars_repo_name: russellnjk/dunlin

max_issues_repo_name: russellnjk/dunlin

max_forks_repo_name: russellnjk/dunlin


text:
import itertools
import numpy as np

class Graph(object):

  def __init__(self, name, sub_graphs):
    self.name = name
    self.sub_graphs = list(sub_graphs)

  def vertex_count(self):
    max_vertexes = map(lambda sub_graph: max(
        sub_graph.vertexes), self.sub_graphs)
    return max(max_vertexes)+1

  def neighbours(self, vertex):
    return [neighbour
            for sub_graph in self.sub_graphs
            for neighbour in sub_graph.neighbours(vertex)]

  def adjacency_matrix(self):
    N = self.vertex_count()
    adjacencies = list(map(lambda sg: sg.adjacency_matrix(N),
                           self.sub_graphs))

    return np.sum(adjacencies, 0)

  def coin_faces(self):
    N = self.vertex_count()
    r = self.max_degree()
    coin_faces = [np.zeros([N, N]) for _ in range(r)]

    for i in range(N):
      neighbours = self.neighbours(i)
      for (index, j) in enumerate(neighbours):
        coin_faces[index][j, i] += 1

    return coin_faces

  def max_degree(self):
    return int(self.adjacency_matrix().sum(axis=0).max())

max_stars_repo_name: nemkin/quantum

max_issues_repo_name: nemkin/quantum

max_forks_repo_name: nemkin/quantum


text:
import numpy as np
import tensorflow as tf
from keras import keras_parameterized, testing_utils
from keras.mixed_precision import policy as mixed_precision
from tensorflow.python.training.tracking import util as trackable_util
from tensorflow.python.util import object_identity
from ..model import build_tri_trans_net, TriTransNet
from ....testing_utils import layer_multi_io_test

@keras_parameterized.run_all_keras_modes
class TestTriTransNet(keras_parameterized.TestCase):
    def setUp(self):
        super(TestTriTransNet, self).setUp()
        self.default_policy = mixed_precision.global_policy()

    def tearDown(self):
        super(TestTriTransNet, self).tearDown()
        mixed_precision.set_policy(self.default_policy)

    def test_layer(self):
        layer_multi_io_test(
            TriTransNet,
            kwargs={},
            input_shapes=[(2, 256, 256, 3), (2, 256, 256, 1)],
            input_dtypes=['uint8', 'uint16'],
            expected_output_shapes=[(None, 256, 256, 1)] * 4,
            expected_output_dtypes=['float32'] * 4
        )

        # TODO: Too slow
        # mixed_precision.set_policy('mixed_float16')
        # layer_multi_io_test(
        #     TriTransNet,
        #     kwargs={},
        #     input_shapes=[(2, 128, 128, 3), (2, 128, 128, 1)],
        #     input_dtypes=['uint8', 'uint16'],
        #     expected_output_shapes=[(None, 128, 128, 1)] * 4,
        #     expected_output_dtypes=['float32'] * 4
        # )

    def test_model(self):
        model = build_tri_trans_net(image_size=128)
        model.compile(
            optimizer='sgd', loss='binary_crossentropy',
            run_eagerly=testing_utils.should_run_eagerly())
        model.fit(
            [np.random.random((2, 128, 128, 3)).astype(np.uint8), np.random.random((2, 128, 128, 1)).astype(np.uint16)],
            np.random.randint(0, 1, (2, 128, 128)),
            epochs=1, batch_size=1)

        # test config
        model.get_config()

        # check whether the model variables are present
        # in the trackable list of objects
        checkpointed_objects = object_identity.ObjectIdentitySet(trackable_util.list_objects(model))
        for v in model.variables:
            self.assertIn(v, checkpointed_objects)

if __name__ == '__main__':
    tf.test.main()

max_stars_repo_name: shkarupa-alex/segme

max_issues_repo_name: shkarupa-alex/segme

max_forks_repo_name: shkarupa-alex/segme


text:
import numpy as np
import tensorflow as tf
from .tensorflow_interp import interpolate3d
from .utils import unit_interval

def construct_beam_grid(uvb, Nza, Naz, freq=None, axis=0, feed=0, spw=0,
                        dtype=tf.float32, **uvbeam_interp_opts):
    """
    Construct a 3D grid of beam values from a UVBeam object, with axes zenith
    angle, azimuth, and frequency. The resulting grid of values can then be
    used by a 3D interpolator.

    The output is a pair of TensorFlow tensors for the real and imaginary parts
    respectively.

    Parameters
    ----------
    uvb : UVBeam object
        A UVBeam object that has been correctly initialised so that the
        `interp()` method can be used.

    Nza, Naz : int
        Number of grid points in the zenith angle and azimuth directions. Note
        that the range of `za` is [0, pi/2] and the range of `az` is [0, 2*pi].

    freq : array_like
        1D arrays of grid points in frequency (in Hz). The beam function is
        evaluated on a regular 3D grid constructed from these 1D arrays.

        If `freq` is `None`, it will be taken from the UVBeam object,
        `uvb.freq_array`.

    axis, feed : int, optional
        Indices of the direction vector (axis) on the sky, and the feed of the
        receiver. How these are defined depends on the input UVBeam object, but
        the axes will normally be 0 or 1 (theta and phi unit vectors), and the
        feeds will normally be 0 or 1 (n or e feeds).
        Defaults: 0.

    spw : int, optional
        Index of the spectral window to use in the UVBeam object. Default: 0.

    **uvbeam_interp_opts : kwargs, optional
        Options to pass to the UVBeam `uvb.interp` method.

    Returns
    -------
    grid_re, grid_im : Tensor, complex
        3D arrays of beam values, with shape (za.size, az.size, uvb.freq_array.size).
        Returned as real and imaginary parts.
    """
    # Define coordinate grids
    az = np.linspace(0., 2.*np.pi, Naz)
    za = np.linspace(0., np.pi/2., Nza)

    # Interpolate beam on 3D grid
    # Returns array of shape: Naxes_vec, Nspws, Nfeeds, Nfreqs, Naz * Nza
    az2d, za2d = np.meshgrid(az, za)
    beam, basis_vec = uvb.interp(az2d.flatten(), za2d.flatten(),
                                 freq_array=np.atleast_1d(freq),
                                 **uvbeam_interp_opts)

    # Reshape grid
    Nfreq = beam.shape[3]
    beam = beam[axis,spw,feed,:,:].reshape(Nfreq, za.size, az.size)
    beam = np.swapaxes(beam, 0, -1) # az, za, freq

    # Output tensors
    grid_re = tf.convert_to_tensor(beam[np.newaxis,:,:,:,np.newaxis].real,
                                   dtype=dtype)
    grid_im = tf.convert_to_tensor(beam[np.newaxis,:,:,:,np.newaxis].imag,
                                   dtype=dtype)
    return grid_re, grid_im,

@tf.function
def coords_for_interp(za, az, freqs, freq_range, grid_shape, dtype=tf.float64):
    """
    Construct a Tensor with coordinates within the unit cube for each zenith
    angle and azimuth, at each frequency. These are the coordinates that the
    beam interpolation will be evaluated at.

    Parameters
    ----------
    za, az : array_like or Tensor
        1D numpy arrays or Tensors containing zenith angle and azimuth
        coordinates of point sources, in radians.

    freqs : array_like or Tensor
        1D numpy array of frequency values (in Hz) that the beam interpolation
        will be evaluated at.

    freq_range : tuple or list
        Tuple or list of len(2), with the maximum and minimum of the frequency
        range that will be covered by the interpolation grid.

    grid_shape : tuple
        Shape of interpolation data grid, (Ngrid_za, Ngrid_az, Ngrid_freq).

    dtype : tf.dtype
        Data type to use for coordinate Tensors. Default: tf.float64.

    Returns
    -------
    sample_pts : Tensor
        Tensor with the correct shape for evaluating a trilinear spline with
        `interpolate3d()`.
    """
    # Unit interval datatype (can be reduced precision)
    unit_type = tf.float64

    # Interpolation grid size
    Ngrid_za, Ngrid_az, Ngrid_freq = grid_shape

    # Convert frequency array to right shape/type
    freqs = tf.squeeze(tf.convert_to_tensor(freqs, dtype=dtype))
    Nfreqs = freqs.shape[0]
    Nptsrc = tf.convert_to_tensor([za.shape[0],], tf.int32) # int32/64 for tf.tile

    # Expand arrays to get full set of sample points
    az_rpt = tf.repeat(tf.convert_to_tensor(az, dtype=dtype), Nfreqs)
    za_rpt = tf.repeat(tf.convert_to_tensor(za, dtype=dtype), Nfreqs)
    freq_rpt = tf.tile(freqs, Nptsrc)

    # Convert arrays to unit interval and cast to different dtype
    az_unit = tf.cast(unit_interval(az_rpt, 0., 2.*np.pi, scale_factor=Ngrid_az),
                      dtype=unit_type)
    za_unit = tf.cast(unit_interval(za_rpt, 0., np.pi/2., scale_factor=Ngrid_za),
                      dtype=unit_type)
    freq_unit = tf.cast(unit_interval(freq_rpt, freq_range[0], freq_range[1],
                                      Ngrid_freq),
                        dtype=unit_type)

    # Combine coords into a structure of the expected shape
    sample_pts = tf.expand_dims(tf.transpose(
                                    tf.stack([az_unit, za_unit, freq_unit]) ),
                                axis=0)
    return sample_pts

@tf.function
def interpolate_beam(grid_re, grid_im, za, az, freqs, freq_range, dtype=tf.float64):
    """
    Interpolate beam on a 3D grid in zenith angle, azimuth, and frequency,
    using TensorFlow routines.

    Parameters
    ----------
    grid_re, grid_im : array_like, complex
        3D array of beam values, with shape (za.size, az.size, uvb.freq_array.size).
        Assumed to be a complex numpy ndarray.

    za, az : array_like or Tensor
        1D numpy arrays or Tensors containing zenith angle and azimuth
        coordinates of point sources, in radians.

    freqs : array_like or Tensor
        1D numpy array of frequency values (in Hz) that the beam interpolation
        will be evaluated at.

    freq_range : tuple or list
        Tuple or list of len(2), with the maximum and minimum of the frequency
        range that will be covered by the interpolation grid.

    dtype : tf.dtype
        Data type to use for coordinate Tensors. Default: tf.float64.

    Returns
    -------
    beam_vals : Tensor, complex
        Complex-valued tensor with beam value at the position of each source,
        for each requested frequency.
    """
    # Put sample points in appropriately-scaled unit cube
    grid_shape = (grid_re.shape[2], grid_re.shape[1], grid_re.shape[3]) # za, az, freq
    sample_pts = coords_for_interp(za, az, freqs, freq_range=freq_range,
                                   grid_shape=grid_shape, dtype=dtype)

    # Interpolate on grid
    beam_re = interpolate3d(grid_re, sample_pts)
    beam_im = interpolate3d(grid_im, sample_pts)
    beam = tf.complex(beam_re, beam_im)

    # Reshape and return
    # FIXME: Need to check that this is being shaped correctly
    return tf.reshape(beam, (za.shape[0], -1))

max_stars_repo_name: philbull/TensorVis

max_issues_repo_name: philbull/TensorVis

max_forks_repo_name: philbull/TensorVis


text:
import numpy as np
import scipy.stats as spstat

"""Metropolis Hastings with "known" covariance"""

def sampler(logpostfunc, options={}):
    '''
    Parameters
    ----------
    logpostfunc : function
        a function returns the log of the posterior for a given theta.
    options : dict, optional
        a dictionary providing the sampler options. The default is {}.
        The possible parameters for the sampler:
        - numsamp (int) : number of samples to draw
        - stepType : either 'uniform' or 'normal'
        - stepParam (float) : scaling parameter
        - theta0 : initial theta value
        - covMat : covariance matrix for selecting candidates #EOW ADD
    Raises
    ------
    ValueError
        If a stepParam or theta0 is not provided.
    Returns
    -------
    sampler_info : dict
        a dictionary contains the output of the sampler.
    '''

    # Initialize
    if 'numsamp' in options.keys():
        n = options['numsamp']
    else:
        n = 2000

    if 'stepType' in options.keys():
        stepType = options['stepType']
    else:
        # default is normal
        stepType = 'normal'

    # scaling parameter
    if 'stepParam' in options.keys():
        stepParam = options['stepParam']
    else:
        raise ValueError('Unknown stepParam')

    # intial theta to start the chain
    if 'theta0' in options.keys():
        thetastart = options['theta0'][0]
    else:
        raise ValueError('Unknown theta0')

    # EOW ADD - number of burn-in iterations
    if 'adaptNum' in options.keys():
        adaptNum = options['adaptNum']
    else:
        adaptNum = 1000
        raise ValueError('Unknown adaptNum')

    p = thetastart.shape[0]
    lposterior = np.zeros(adaptNum + n)
    theta = np.zeros((adaptNum + n, p))
    lposterior[0] = logpostfunc(thetastart)
    theta[0, :] = thetastart
    n_acc = 0
    MVN = spstat.multivariate_normal(mean=None, cov=options['covMat'])

    for i in range(1, adaptNum + n):
        # Candidate theta
        if stepType == 'normal':
            theta_cand = theta[i-1, :] + stepParam[0] * MVN.rvs()
            #theta_cand = [theta[i-1, :][k] + stepParam[k] *
                          #spstat.multivariate_normal(mean=None, cov=options['covMat'],
                                                     #size=1) for k in range(p)]
                          #spstat.norm.rvs(0, 1, size=1) for k in range(p)]
        elif stepType == 'uniform':
            theta_cand = [theta[i-1, :][k] + stepParam[k] *
                          spstat.uniform.rvs(-0.5, 0.5, size=1) for k in range(p)]

        theta_cand = np.reshape(np.array(theta_cand), (1, p))

        # Compute loglikelihood
        logpost = logpostfunc(theta_cand)

        if np.isfinite(logpost):
            p_accept = min(1, np.exp(logpost - lposterior[i-1]))
            accept = np.random.uniform() < p_accept
        else:
            accept = False

        # Accept candidate?
        if accept:
            # Update position
            theta[i, :] = theta_cand
            lposterior[i] = logpost
            if i >= adaptNum:
                n_acc += 1
        else:
            theta[i, :] = theta[i-1, :]
            lposterior[i] = lposterior[i-1]

    theta = theta[(adaptNum):(adaptNum + n), :]
    sampler_info = {'theta': theta, 'acc_rate': n_acc/n}
    return sampler_info


text:
import numpy as np
import isingify
import pandas as pd
import minorminer as mm
import dimod
from dwave.system.composites import FixedEmbeddingComposite
from dwave.system.samplers import DWaveSampler
import os

class Ising_SVP_Experiment():

    def __init__(self, dimension, k, times, n_reads, qudit_mapping):
        self.dim = dimension
        self.k = k
        self.times = times
        self.n_reads = n_reads
        self.qudit_mapping = qudit_mapping
        if self.qudit_mapping == 'bin':
            self.qubits_per_qudit = np.ceil(np.log2(2 * self.k))
        else:  # 'ham'
            self.qubits_per_qudit = 2 * self.k
        self.nqubits = self.dim * self.qubits_per_qudit
        self.t = len(self.times.values())
        """
        self.scaling_coefficient should be chosen so that big enough that you achieve almost no chain breaks,
        but at the same time small enough that you are still expanding the eigenspectrum versus the binary mapping
        # i.e. it should be less than log2(2k) times the scaling coefficient used for log, which here is
        """
        self.scaling_coefficient = 8

    def bitstr_to_coeff_vector(self, bitstr):
        """
        Parameters
        ----------
        bitstr - a list of int Pauli-Z eigenstates (+/- 1) which can be interpretted as a bitstring through
         application of a simple operator 0.5(1-Z)

        Returns
        -------
        str - comma-delimited coefficient vector as a string
        """

        if self.qudit_mapping == "bin":
            ind = 0
            vect = []
            for i in range(self.dim):
                num = -2 ** (self.qubits_per_qudit - 1)
                for j in range(int(self.qubits_per_qudit)):
                    num += bitstr[ind] * (2 ** (j))
                    ind += 1
                vect.append(num)
            return (','.join(list(map(str, [int(y) for y in vect]))))

        else:  # 'ham'
            ind = 0
            vect = []
            for i in range(self.dim):
                num = 0
                for j in range(int(self.qubits_per_qudit)):
                    num += bitstr[ind] / 2
                    ind += 1
                vect.append(num)
            return (','.join(list(map(str, [int(y) for y in vect]))))

    def preprocess(self, lattice):
        """
        Used to map the input lattice to properties used by the D-Wave API, i.e. the interaction strengths
        Jmat, hvec, ic are explained in detail in isingify.py

        Parameters
        ----------
        lattice - numpy 2D array of type int representing a row-basis to input to the algorithm

        Returns
        -------
        """

        self.Gram = lattice@lattice.T

        if self.qudit_mapping == "bin":
            self.Jmat, self.hvec, self.ic = isingify.svp_isingcoeffs_bin(self.Gram, self.k)
        else:
            self.Jmat, self.hvec, self.ic = isingify.svp_isingcoeffs_ham(self.Gram, self.k)

        self.scale = np.max((np.abs(self.Jmat).max(), np.abs(self.hvec / 2).max())) * self.scaling_coefficient
        self.l = {i: self.hvec[i] / self.scale for i in range(self.hvec.shape[0])}
        self.q = {(i, j): self.Jmat[i, j] / self.scale for i in range(self.Jmat.shape[0]) for j in
                  range(self.Jmat.shape[1])
                  if
                  not i == j}

    def embed_graph(self):
        """
        Separate function as graph embedding can be reused so long as self.dim, self.k remain constant
        Embed graph maps the fully connected Ising model we have defined to the D-Wave Chimera topology
        """
        source_edgelist = [(i, j) for i in range(int(self.nqubits)) for j in range(i, int(self.nqubits))]
        dws = DWaveSampler(qpu=True)  # instantiate DWaveSampler
        dws_structure = dimod.child_structure_dfs(dws)
        target_edgelist = dws_structure.edgelist
        self.embedding = mm.find_embedding(source_edgelist, target_edgelist)
        physical_qubits = []
        for qubit_list in self.embedding.values():  # this loop just puts all the physical qubit indices into one list
            physical_qubits += qubit_list
        print("Number of physical qubits: ", len(physical_qubits))

    def execute_dwave(self):
        """
        Executes request to D-Wave API using predefined Ising_SVP_Experiment attributes (from preprocess())
        Saves results to response_dict attribute
        """
        self.response_dict = {}
        for i in range(self.t):
            response = FixedEmbeddingComposite(DWaveSampler(qpu=True),
                                               embedding=self.embedding).sample_ising(self.l,
                                                                                self.q,
                                                                                num_reads=self.n_reads,
                                                                                annealing_time=
                                                                                self.times[i],
                                                                                auto_scale=False,
                                                                                chain_strength=2.0)
            self.response_dict[self.times[i]] = response

    def postprocess(self, lattice_index):
        """
        Parameters
        ----------
        lattice_index - int index of the lattice experiment

        Returns
        -------
        N/A - processes D-Wave results and saves to CSV in file 'data'
        """
        for i in range(self.t):
            raw = list(self.response_dict[self.times[i]].data())
            if self.qudit_mapping == 'bin':
                bitlists_int = [[int(y) for y in [(1 - x) / 2 for x in list(samp[0].values())]] for samp in raw]
                bitlists_str = [''.join(list(map(str, [int(y) for y in [(1 - x) / 2 for x in list(samp[0].values())]])))
                                for
                                samp in raw]
            else:  # 'ham'
                bitlists_int = [[int(y) for y in list(samp[0].values())] for samp in raw]
                bitlists_str = [''.join(list(map(str, [int(y) for y in list(samp[0].values())])))
                                for samp in raw]
            vects = [self.bitstr_to_coeff_vector(bitlist) for bitlist in bitlists_int]
            eigenenergies = [samp[1] for samp in raw]
            occurrences = [samp[2] for samp in raw]
            chain_break = [samp[3] for samp in raw]
            df = pd.DataFrame(data=list(zip(eigenenergies, occurrences, bitlists_str, vects, chain_break)),
                              columns=['eigenenergy', 'occurrences', 'bitstring', 'coefficient_vectors', 'chain_break'])
            df['vector_lengths_sq'] = df['eigenenergy'].apply(lambda x: int(round((x * self.scale) + self.ic)))
            path = os.path.join('data',
                                '{}D_sample{}_T{}_k{}.csv'.format(self.dim, lattice_index, self.times[i], self.k))
            df.to_csv(path)

    def run(self):
        """ Preprocess, Execute, Postprocess D-Wave experiment """

        # use same embedding for all lattices of same dimension
        self.embed_graph()

        # run experiment and save results for each lattice, for a range of times
        # each lattice_index corresponds to a different random lattice to experiment on
        for lattice_index in range(2):
            lattice = np.random.randint(low=0, high=5, size=(self.dim, self.dim), dtype=int)
            self.preprocess(lattice_index, lattice)
            self.execute_dwave()
            self.postprocess(lattice_index, lattice)

def main():

    """
    Set experiment parameters

    dimension - int: lattice dimension
    k - int:  qubits per qudit
    reps - int: number of D-Wave samples collected
    qudit_mapping - str: binary or Hamming encoded qudits
    times - dict of key/value type int: keys are sweeptime indices, values are sweeptimes in microseconds
    """

    dimension = 3
    k = 4
    reps = 90
    qudit_mapping = "ham"
    times = {0: 1}

    """ Initialise Ising_SVP_Experiment object with parameters set above """
    experiment = Ising_SVP_Experiment(dimension, k, times, reps, qudit_mapping)
    experiment.run()

if __name__ == "__main__":
    main()

max_stars_repo_name: adambrzosko/D-Wave-Ising-SVP

max_issues_repo_name: adambrzosko/D-Wave-Ising-SVP

max_forks_repo_name: adambrzosko/D-Wave-Ising-SVP


text:
"""
Comparison of a UW model with a GC analysis
$Header: /nfs/slac/g/glast/ground/cvs/pointlike/python/uw/like2/analyze/gc_comparison.py,v 1.3 2017/11/17 22:43:17 burnett Exp $

"""

import os
from astropy.io import fits as pyfits
import numpy as np
import pylab as plt
import pandas as pd
from astropy.io import fits

from skymaps import SkyDir
# FIX LATER from matplotlib_venn import venn2
from . import sourceinfo
from . analysis_base import FloatFormat, html_table

_filenames=dict(
    GC300='GCPass8_listsources_cleaned07deg_TS10_tsmap_300MeVref_fit3TS10.fits',
    GC500='GCPass8_listsources_all_TS10_tsmap_500MeV_cleaned02deg_fit3TS10_ref.fits',
)

class GCcomparison(sourceinfo.SourceInfo):
    """Comparison with a GC analysis

    %(setup_info)s
    """

    def setup(self, gc_root = '/nfs/farm/g/glast/u/burnett/analysis/GC_check/',
        cols = 'ra dec glat glon ts pos_sigma Spectral_Index'.split(),
        gcname = 'GC300', **kw):

        super(GCcomparison, self).setup(**kw)
        self.plotfolder='comparison_%s' % gcname
        filename = _filenames[gcname]

        # load the GC analysis
        assert os.path.exists(gc_root + filename)
        hdulist =fits.open(gc_root + filename)
        data = hdulist[1].data
        df = pd.DataFrame([data.field(x) for x in cols], index=cols).T
        df.index=data.field('Source_name');
        df.index.name='name'
        df['sd'] = map(SkyDir, df.ra,df.dec)
        print 'Created a DataFrame from {} with {} rows'.format(filename, len(df))
        self.dfgc = df
        self.gcname=gcname
        self.setup_info='Comparing sources in file {} with {}'\
            .format(filename, self.skymodel)
        self.select_uw()

    def gc_plots(self):
        """GC plots

        Using all entries
        """
        df = self.dfgc
        fig,axx =plt.subplots(2,2, figsize=(12,12))
        ax = axx[0,0]
        ax.hist(df.ts, np.logspace(1,5,41), histtype='step');
        plt.setp(ax, xscale='log', xlabel='ts');
        ax.grid(alpha=0.5)
        ax=axx[0,1]
        df.glon[df.glon>180]-=360
        ax.plot(df.glon, df.glat, '.')
        plt.setp(ax, xlabel='glon', ylabel='glat', xlim=(25,-25), ylim=(-25,25))
        ax.grid(alpha=0.5)
        ax=axx[1,0]
        ax.loglog(df.ts, df.pos_sigma, '.')
        plt.setp(ax, xlabel='ts', ylabel='pos_sigma');
        ax.grid(alpha=0.5)
        ax=axx[1,1]
        ax.hist(df.Spectral_Index, np.linspace(1,4,31), histtype='step');
        plt.setp(ax, xlabel='Spectral Index')
        ax.grid(alpha=0.5)
        fig.suptitle('{} plots'.format(self.gcname), fontsize=16)
        return fig

    def select_uw(self, angle_cut=21):
        dfuwall = self.df
        glon=dfuwall.glon
        glon[glon>180]-=360
        glat=dfuwall.glat
        dfuw = dfuwall[(abs(glon)<angle_cut) & (abs(glat)<angle_cut)]
        dfuw['sd'] = map(SkyDir, dfuw.ra, dfuw.dec)
        dfuw['pos_sigma'] = (dfuw.a * dfuw.b) ** 0.5
        text = 'Selected {} uw sources within box of {} deg about GC'.format(len(dfuw),angle_cut)
        print text; self.setup_info +='<br>'+text
        self.dfuw = dfuw

    def correlate(self, angle_cut=0.1):
        """Correlate the two lists
        """
        def differences(a,b):
            matrix = np.array([[x.difference(y) for y in b] for x in a], np.float32)
            return matrix
        def closest(t):
            n = t.argmin()
            return (n, (np.degrees(t[n])*3600).round())

        diff_array =differences(self.dfgc.sd, self.dfuw.sd)
        cl_gc = np.array([closest(r) for r in diff_array[:,]], int)
        cl_uw = np.array([closest(c) for c in diff_array[:,].T], int)

        # add correlation info to the GC dataframe
        close_cut = 3600*angle_cut
        df = self.dfgc
        dfuw = self.dfuw
        df['uw_id'] = cl_gc[:,0]
        df['dist'] = cl_gc[:,1]
        df['uw_ts'] = [dfuw.ix[i].ts for i in cl_gc[:,0]]
        df['uw_pos_sigma'] = [dfuw.ix[i].pos_sigma for i in cl_gc[:,0]]
        df['uw_pindex'] = [dfuw.ix[i].pindex for i in cl_gc[:,0]]
        df['uw_roi'] = [int((dfuw.ix[i].roiname)[-4:]) for i in cl_gc[:,0]]
        df['uwok'] = df.dist<close_cut
        df['nouw'] = np.logical_not(df.uwok)
        print 'GC sources associated with uw: {}/{}'.format(sum(df.uwok), len(df))

        #add correlation info to the UW dataframe
        dfuw['gc_id'] = np.array(cl_uw[:,0])
        dfuw['dist'] = np.array(cl_uw[:,1])

        dfuw['gc_ts'] = [df.ix[i].ts for i in cl_uw[:,0]]
        dfuw['gc_pos_sigma'] = [df.ix[i].pos_sigma for i in cl_uw[:,0]]
        dfuw['gc_pindex'] = [df.ix[i].Spectral_Index for i in cl_uw[:,0]]
        dfuw['gcok'] = dfuw.dist<close_cut
        dfuw['nogc'] = np.logical_not(dfuw.gcok)
        print 'UW sources associated with GC: {}/{}'.format(sum(dfuw.gcok), len(dfuw))

        fig, axx =plt.subplots(1,2, figsize=(12,6))
        ax=axx[0]
        hist_kw=dict(bins=np.linspace(0,0.5, 26), histtype='step', log=True)
        ax.hist(cl_gc[:,1].clip(0,1800)/3600.,
            label='{} [{}]'.format(self.gcname,len(cl_gc)), **hist_kw);
        ax.hist(cl_uw[:,1].clip(0,1800)/3600.,
            label='{} [{}]'.format(self.skymodel,len(cl_uw)), color='orange',**hist_kw);
        ax.axvline(angle_cut, color='red', ls='--', label='angle cut')
        plt.setp(ax, ylim=(0.8,None), xlabel='distance (deg)', title='minimum distance')
        ax.legend(loc='upper left');
        ax.grid(alpha=0.5)

        # make a Venn diagram
        common = sum(cl_uw[:,1]<close_cut)
        venn2(ax=axx[1], subsets=(len(self.dfuw)-common, len(self.dfgc)-common,common)
        ,set_labels=(self.skymodel,self.gcname))
        return fig

    def common_plots(self):
        """Common plots

        """
        df = self.dfgc
        fig, axx = plt.subplots(2,2, figsize=(12,12))
        ax = axx[0,0]
        dfok = df[df.uwok]
        ax.loglog(dfok.ts,dfok.uw_ts, '+b');
        ax.plot([10,1e6], [10,1e6], '--g');
        plt.setp(ax, xlabel='GC ts', ylabel='UW ts')
        ax.grid(alpha=0.5)
        ax = axx[0,1]
        lim = (1e-3,0.1)
        ax.loglog(dfok.pos_sigma.clip(*lim), dfok.uw_pos_sigma.clip(*lim), '+b');
        ax.plot(lim, lim, '--g')
        plt.setp(ax, xlim=lim, ylim=lim, xlabel='GC pos_sigma', ylabel='UW pos_sigma');
        ax.grid(alpha=0.5);
        ax = axx[1,0]
        lim=(1.5,3.5)
        ax.plot(dfok.Spectral_Index.clip(*lim), dfok.uw_pindex.clip(*lim), '+b')

        plt.setp(ax, xlim=lim, ylim=lim, xlabel='GC Spectral Index', ylabel='UW pindex')
        ax.plot(lim, lim, '--g')
        ax.grid(alpha=0.5)
        axx[1,1].set_visible(False)
        fig.suptitle('Comparison of values for common sources', fontsize=14)
        return fig

    def nouwcounterpart(self):
        """Plots for GC sources with no UW counterpart
        """
        df=self.dfgc
        fig,axx =plt.subplots(2,2, figsize=(12,12))
        ax = axx[0,0]
        ax.hist(df.ts, np.logspace(1,5,21), histtype='step', label='all');
        ax.hist(df[df.nouw].ts, np.logspace(1,5,21), histtype='step', color='red', label='no UW');
        plt.setp(ax, xscale='log', xlabel='ts');
        ax.legend()
        ax.grid(alpha=0.5)
        ax=axx[0,1]
        df.glon[df.glon>180]-=360
        ax.plot(df.glon, df.glat, '.')
        ax.plot(df[df.nouw].glon, df[df.nouw].glat, '.r')
        plt.setp(ax, xlabel='glon', ylabel='glat', xlim=(25,-25), ylim=(-25,25))
        ax.grid(alpha=0.5)
        ax=axx[1,0]
        ax.plot(df.glon, df.glat, '.')
        ax.plot(df[df.nouw].glon, df[df.nouw].glat, '.r')
        plt.setp(ax, xlabel='glon', ylabel='glat', xlim=(25,-25), ylim=(-5,5))
        ax.grid(alpha=0.5)
        ax=axx[1,1]
        ax.hist(df.Spectral_Index, np.linspace(1,4,31), histtype='step');
        ax.hist(df[df.nouw].Spectral_Index, np.linspace(1,4,31), histtype='step', color='red');
        ax.grid(alpha=0.5)
        plt.setp(ax, xlabel='Specral Index')
        fig.suptitle('{} distributions with no {} counterpart'.format(self.gcname, self.skymodel), fontsize=16);
        return fig

    def nogccounterpart(self):
        """Plots for UW sources with no GC counterpart
        """
        dfuw=self.dfuw
        fig,axx =plt.subplots(2,2, figsize=(12,12))
        ax = axx[0,0]
        ax.hist(dfuw.ts, np.logspace(1,5,41), histtype='step');
        ax.hist(dfuw[dfuw.nogc].ts, np.logspace(1,5,41), histtype='step', color='red', label='no GC');
        plt.setp(ax, xscale='log', xlabel='ts');
        ax.legend()
        ax.grid(alpha=0.5)
        ax=axx[0,1]
        dfuw.glon[dfuw.glon>180]-=360
        ax.plot(dfuw.glon, dfuw.glat, '.')
        ax.plot(dfuw[dfuw.nogc].glon, dfuw[dfuw.nogc].glat, '.r')
        plt.setp(ax, xlabel='glon', ylabel='glat', xlim=(25,-25), ylim=(-25,25))
        ax.grid(alpha=0.5)
        ax=axx[1,0]
        ax.plot(dfuw.glon, dfuw.glat, '.')
        ax.plot(dfuw[dfuw.nogc].glon, dfuw[dfuw.nogc].glat, '.r')
        plt.setp(ax, xlabel='glon', ylabel='glat', xlim=(25,-25), ylim=(-5,5))
        ax.grid(alpha=0.5)
        ax=axx[1,1]
        ax.hist(dfuw.pindex, np.linspace(1,4,31), histtype='step');
        ax.hist(dfuw[dfuw.nogc].pindex, np.linspace(1,4,31), histtype='step', color='red');
        ax.grid(alpha=0.5)
        plt.setp(ax, xlabel='Spectral Index')
        fig.suptitle('{} sources with no {} counterpart'.format(self.skymodel, self.gcname),
            fontsize=16);
        return fig

    def gc_subset_table(self, tsmin=25, bmin=0,
        cols='ra dec glon glat ts pos_sigma Spectral_Index dist uw_roi'.split()
        ):
        """Table of selected strong GC sources not in the UW list

        %(selected_sources)s
        """
        df = self.dfgc
        selection = df.nouw & (df.ts>tsmin) & (np.abs(df.glat)>bmin)
        dfsel = df[selection][cols]
        # replace all prefixes
        split_names = [n.split()[1] for n in dfsel.index]
        dfsel.index = [self.gcname+' '+sn for sn in split_names]
        dfsel.index.name='name'
        self.selected_sources=\
            html_table(
                dfsel,
                 href=False,
                 heading = '<p>%d GC sources with TS>%.0f, |b|>%.1f' % (sum(selection),tsmin, bmin),
                    name=self.plotfolder+'/selected', maxlines=20,
                    float_format=(FloatFormat(3))
            )
        outfile = self.plotfolder+'/selected.csv'
        dfsel.to_csv(outfile)
        print 'Wrote {} seeds to file {}'.format(sum(selection), outfile)

    def all_plots(self):
        """
        %(setup_info)s
        """
        self.runfigures([
                self.gc_plots, self.correlate, self.common_plots, self.nouwcounterpart, self.nogccounterpart,self.gc_subset_table,
        ])

max_stars_repo_name: tburnett/pointlike

max_issues_repo_name: tburnett/pointlike

max_forks_repo_name: tburnett/pointlike

reason: import numpy,from astropy


text:
from UpDownMethods import CORRECT, INCORRECT
import UpDownMethods as ud
import numpy as np
import matplotlib.pyplot as plt
import unittest

#
# Simulation parameters

responses = [CORRECT, CORRECT, INCORRECT, INCORRECT, INCORRECT, INCORRECT,
             CORRECT, CORRECT, INCORRECT, INCORRECT, INCORRECT, CORRECT,
             CORRECT, CORRECT, CORRECT, CORRECT, INCORRECT, INCORRECT,
             INCORRECT, CORRECT, CORRECT, INCORRECT, INCORRECT, INCORRECT]

initalValue = 0.0
stepSize = 1.0
down = 1
up = 1

#
# Test code

class TestSequenceFunctions(unittest.TestCase):

    def setUp(self):
        self.results = ud.initiate_procedure()
        nextValue, self.results = ud.append_result(self.results, responses[0],
                                                   down, up, stepSize,
                                                   initalValue)
        for resp in responses[1:]:
            nextValue, self.results = ud.append_result(self.results, resp,
                                                       down, up, stepSize,
                                                       nextValue)
            ud.midpoints(self.results)
            ud.reversals(self.results)
            ud.estimate_reversals(self.results)
            ud.runs(self.results)

    def test_initiateResults(self):
        self.results = ud.initiate_procedure()
        self.assertIs(len(self.results), 0)

    def test_calculateMidpoints(self):
        mids = ud.midpoints(self.results)
        mids = mids["Midpoint"]
        mids = mids[[1, 3, 5]].values   # Runs 2, 4, 6
        self.assertIsNone(np.testing.assert_array_equal(mids, [0., 1.5, -0.5]))

    def test_plotResults(self):
        ud.plot_results(self.results, midpoints=True)
        plt.savefig('doc/images/Levitt-Fig4.png', bbox_inches='tight')

    def test_reversals(self):
        revs = ud.reversals(self.results)
        self.assertIsNone(np.testing.assert_array_equal(
           revs["Reversal"].values, [1, 2, 3, 4, 5, 6, 7]))
        self.assertIsNone(np.testing.assert_array_equal(
           revs["Trial"].values, [3, 7, 9, 12, 17, 20, 22]))
        self.assertIsNone(np.testing.assert_array_equal(
           revs["Value"].values, [-2, 2, 0, 3, -2, 1, -1]))

    def test_estimate_reversals(self):
        est = ud.estimate_reversals(self.results)
        self.assertEqual(est, 0)
        est = ud.estimate_reversals(self.results, num=4)
        self.assertEqual(est, 0.25)

    def test_runs(self):
        runs = ud.runs(self.results)
        self.assertIsNone(np.testing.assert_array_equal(
           runs["Start"].values, [1, 3, 7, 9, 12, 17, 20, 22]))
        self.assertIsNone(np.testing.assert_array_equal(
           runs["Finish"].values, [3, 7, 9, 12, 17, 20, 22, 24]))

if __name__ == '__main__':
    suite = unittest.TestLoader().loadTestsFromTestCase(TestSequenceFunctions)
    unittest.TextTestRunner(verbosity=2).run(suite)

max_stars_repo_name: codles/UpDownMethods

max_issues_repo_name: codles/UpDownMethods

max_forks_repo_name: codles/UpDownMethods