text:
from abc import ABCMeta, abstractmethod

import numpy as np
import tensorflow as tf
from skimage.draw import circle
import Constants as Constants
import datasets.Util.Util as Util
from Log import log
from datasets.Augmentors import parse_augmentors
from datasets.Util.Batch import create_batch_dict
from datasets.Util.Normalization import unnormalize
from datasets.Util.Reader import read_images_from_disk, load_label_default, load_img_default
from datasets.Util.Resize import parse_resize_mode, ResizeMode

class Dataset(object, metaclass=ABCMeta):
    def __init__(self, subset):
        self.summaries = []
        self.subset = subset

    def _get_resize_params(self, subset, input_size_default, resize_mode_default=ResizeMode.Unchanged):
        if subset == "valid":
            subset = "val"
        resize_mode_string = self.config.str("resize_mode_" + subset, "")
        # use same as for train if nothing specified
        if resize_mode_string == "":
            resize_mode_string = self.config.str("resize_mode_train", "")
        if resize_mode_string == "":
            resize_mode = resize_mode_default
        else:
            resize_mode = parse_resize_mode(resize_mode_string)
        input_size = self.config.int_list("input_size_" + subset, [])
        # use same as for train if nothing specified
        if len(input_size) == 0:
            if resize_mode == ResizeMode.Unchanged:
                input_size = [None, None]
            else:
                input_size = self.config.int_list("input_size_train", [])
        if len(input_size) == 0:
            input_size = input_size_default
        if resize_mode == ResizeMode.RandomResize:
            input_size = [None, None]
        return resize_mode, input_size

    def _parse_augmentors_and_shuffle(self):
        if self.subset == "train":
            shuffle = True
            augmentor_strs = self.config.unicode_list("augmentors_train", [])
            augmentors = parse_augmentors(augmentor_strs, self.void_label())
            if len(augmentors) == 0:
                print("warning, no data augmentors used on train", file=log.v1)
        else:
            shuffle = False
            augmentor_strs = self.config.unicode_list("augmentors_val", [])
            augmentors = parse_augmentors(augmentor_strs, self.void_label())
        return augmentors, shuffle

    @abstractmethod
    def num_classes(self):
        pass

    @abstractmethod
    def num_examples_per_epoch(self):
        pass

    # should contain inputs, labels, tags, maybe raw_labels, index_img
    @abstractmethod
    def create_input_tensors_dict(self, batch_size):
        pass

    @abstractmethod
    def void_label(self):
        pass

class ImageDataset(Dataset):
    def __init__(self, dataset_name, default_path, num_classes, config, subset, coord, image_size, void_label=255,
                 fraction=1.0, label_postproc_fn=lambda x: x, label_load_fn=load_label_default,
                 img_load_fn=load_img_default, ignore_classes=[]):
        super(ImageDataset, self).__init__(subset)
        self._num_classes = num_classes
        self._void_label = void_label
        assert subset in ("train", "valid"), subset
        if subset == "train":
            self.epoch_length = config.int("epoch_length", - 1)
        else:
            self.epoch_length = -1
        self.config = config
        self.data_dir = config.str(dataset_name + "_data_dir", default_path)
        self.coord = coord
        self.image_size = image_size
        self.inputfile_lists = None
        self.fraction = fraction
        self.label_postproc_fn = label_postproc_fn
        self.label_load_fn = label_load_fn
        self.img_load_fn = img_load_fn
        self.use_summaries = self.config.bool("use_summaries", False)
        self.ignore_classes = ignore_classes

    def _load_inputfile_lists(self):
        if self.inputfile_lists is not None:
            return
        self.inputfile_lists = self.read_inputfile_lists()
        # make sure all lists have the same length
        assert all([len(l) == len(self.inputfile_lists[0])
                    for l in self.inputfile_lists])
        if self.fraction < 1.0:
            n = int(self.fraction * len(self.inputfile_lists[0]))
            self.inputfile_lists = [l[:n] for l in self.inputfile_lists]

    def create_input_tensors_dict(self, batch_size):
        self._load_inputfile_lists()
        resize_mode, input_size = self._get_resize_params(
            self.subset, self.image_size, ResizeMode.Unchanged)
        augmentors, shuffle = self._parse_augmentors_and_shuffle()

        #inputfile_tensors = [tf.convert_to_tensor(l, dtype=tf.string) for l in self.inputfile_lists]
        inputfile_tensors = [tf.convert_to_tensor(
            l) for l in self.inputfile_lists]
        queue = tf.train.slice_input_producer(
            inputfile_tensors, shuffle=shuffle)

        tensors_dict, summaries = self._read_inputfiles(
            queue, resize_mode, input_size, augmentors)
        tensors_dict, summ = create_batch_dict(batch_size, tensors_dict)
        if summ is not None:
            summaries.append(summ)
        if Constants.BBOXES in tensors_dict and Constants.IDS in tensors_dict and Constants.CLASSES in tensors_dict:
            assert "labels" not in tensors_dict
            tensors_dict["labels"] = (tensors_dict[Constants.BBOXES], tensors_dict[Constants.IDS],
                                      tensors_dict[Constants.CLASSES])
            if Constants.IGNORE_REGIONS in tensors_dict:
                tensors_dict["labels"] += (
                    tensors_dict[Constants.IGNORE_REGIONS],)
                if Constants.SCENE_INFOS in tensors_dict:
                    tensors_dict["labels"] += (
                        tensors_dict[Constants.SCENE_INFOS],)

        if self.use_summaries:
            summaries += self._create_summaries(tensors_dict)

        self.summaries += summaries
        return tensors_dict

    def _create_summaries(self, tensors_dict):
        inputs = tensors_dict["inputs"]
        input_imgs = unnormalize(inputs[:, :, :, :3])
        # count is incremented after each summary creation. This helps us to keep track of the number of channels
        start = 3
        summaries = []
        # Old label would be present either as a single channel, or along with 2 other distance transform channels.
        if inputs.get_shape()[-1] in [4, 6, 8]:
            if Constants.OLD_LABEL_AS_DT in tensors_dict:
                [input_imgs] = tf.py_func(self.add_clicks, [tf.concat([input_imgs, inputs[:, :, :, start:start + 1]],
                                                                      axis=3), 'g'], [tf.float32])
            else:
                [old_label, input_imgs] = tf.py_func(self.get_masked_image, [input_imgs, inputs[:, :, :, start:start + 1]],
                                                     [tf.float32, tf.float32])

            summ = tf.summary.image(
                "old_labels", inputs[:, :, :, start:start + 1])
            summaries.append(summ)
            start += 1

        # Add clicks to the image so that they can be viewed in tensorboard
        if inputs.get_shape()[-1] > 6:
            [input_imgs] = tf.py_func(self.add_clicks, [tf.concat([input_imgs, inputs[:, :, :, start:start + 1]],
                                                                  axis=3), 'r'], [tf.float32])
            [input_imgs] = tf.py_func(self.add_clicks, [tf.concat([input_imgs, inputs[:, :, :, start + 1:start + 2]],
                                                                  axis=3), 'r'], [tf.float32])
            [input_imgs] = tf.py_func(self.add_clicks, [tf.concat([input_imgs, inputs[:, :, :, start + 2:start + 3]],
                                                                  axis=3), 'g'], [tf.float32])
            [input_imgs] = tf.py_func(self.add_clicks, [tf.concat([input_imgs, inputs[:, :, :, start + 3:start + 4]],
                                                                  axis=3), 'g'], [tf.float32])
        elif inputs.get_shape()[-1] > 4:
            [input_imgs] = tf.py_func(self.add_clicks, [tf.concat([input_imgs, inputs[:, :, :, start:start + 1]],
                                                                  axis=3), 'r'], [tf.float32])
            [input_imgs] = tf.py_func(self.add_clicks, [tf.concat([input_imgs, inputs[:, :, :, start + 1:start + 2]],
                                                                  axis=3), 'g'], [tf.float32])

        # bounding boxes
        if Constants.BBOXES in tensors_dict:
            bboxes = tensors_dict[Constants.BBOXES]
            # permute y1, y2, x1, x2 -> y1, x1, y2, x1
            bboxes = tf.stack([bboxes[..., 0], bboxes[..., 2],
                               bboxes[..., 1], bboxes[..., 3]], axis=-1)
            # normalize bboxes to [0..1]
            height = tf.shape(input_imgs)[1]
            width = tf.shape(input_imgs)[2]
            bboxes = tf.cast(bboxes, tf.float32) / tf.cast(
                tf.stack([height, width, height, width], axis=0), tf.float32)
            imgs_with_bboxes = tf.image.draw_bounding_boxes(input_imgs, bboxes)
            summ = tf.summary.image("inputs", imgs_with_bboxes)
            summaries.append(summ)
        else:
            summ = tf.summary.image("inputs", input_imgs)
            summaries.append(summ)

        if not isinstance(tensors_dict["labels"], tuple):
            # will only work well for binary segmentation
            summ = tf.summary.image("labels", tensors_dict["labels"] * 255)
            summaries.append(summ)

        # Append the distance transforms, if they are available.
        if inputs.get_shape()[-1] > 6:
            # Get negative distance transform from the extra input channels.
            summ = tf.summary.image(
                Constants.DT_NEG, inputs[:, :, :, start:start + 1])
            summaries.append(summ)
            summ = tf.summary.image(
                Constants.DT_NEG, inputs[:, :, :, start+1:start + 2])
            summaries.append(summ)
            summ = tf.summary.image(
                Constants.DT_POS, inputs[:, :, :, start + 2:start + 3])
            summaries.append(summ)
            summ = tf.summary.image(
                Constants.DT_POS, inputs[:, :, :, start + 3:start + 4])
            summaries.append(summ)
        elif inputs.get_shape()[-1] > 4:
            # Get negative distance transform from the extra input channels.
            summ = tf.summary.image(
                Constants.DT_NEG, inputs[:, :, :, start:start + 1])
            summaries.append(summ)
            summ = tf.summary.image(
                Constants.DT_POS, inputs[:, :, :, start + 1:start + 2])
            summaries.append(summ)

        return summaries

    # default implementation, should in many cases be overwritten
    def _read_inputfiles(self, queue, resize_mode, input_size, augmentors):
        tensors, summaries = read_images_from_disk(queue, input_size, resize_mode, label_postproc_fn=self.label_postproc_fn,
                                                   augmentors=augmentors, label_load_fn=self.label_load_fn,
                                                   img_load_fn=self.img_load_fn)
        return tensors, summaries

    def add_clicks(self, inputs, c):
        out_imgs = None

        for input in inputs:
            img = input[:, :, :3]
            # Radius of the point to be diplayed
            r = 3
            pts = np.where(input[:, :, 3] == 0)
            pts_zipped = list(zip(pts[0], pts[1]))
            if len(pts[0]) > 0:
                for pt in pts_zipped:
                    if r < pt[0] < img.shape[0] - r and r < pt[1] < img.shape[1] - r:
                        rr, cc = circle(pt[0], pt[1], 5, img.shape)
                        img[rr, cc, :] = [np.max(img), np.min(img), np.min(img)] if c == 'r' \
                            else [np.min(img), np.min(img), np.max(img)]

            img = img[np.newaxis, :, :, :]
            if out_imgs is None:
                out_imgs = img
            else:
                out_imgs = np.concatenate((out_imgs, img), axis=0)

        return out_imgs.astype(np.float32)

    def get_masked_image(self, inputs, old_labels):
        out_imgs = None
        out_labels = None
        for input_img, old_label_unmodified in zip(inputs, old_labels):
            if len(np.unique(old_label_unmodified)) > 2:
                old_label = np.where(old_label_unmodified*255 <= 128, 1, 0)
            else:
                old_label = old_label_unmodified

            out_img = Util.get_masked_image(input_img, old_label[:, :, 0])
            out_img = out_img[np.newaxis, :, :, :]
            out_label = old_label[np.newaxis, :, :, :]
            if out_imgs is None:
                out_imgs = out_img
            else:
                out_imgs = np.concatenate((out_imgs, out_img), axis=0)

            if out_labels is None:
                out_labels = out_label
            else:
                out_labels = np.concatenate((out_labels, out_label), axis=0)

        return out_labels.astype(np.float32), out_imgs.astype(np.float32)

    @abstractmethod
    def read_inputfile_lists(self):
        pass

    def num_examples_per_epoch(self):
        if self.epoch_length != -1:
            return self.epoch_length
        else:
            self._load_inputfile_lists()
            return len(self.inputfile_lists[0])

    def num_classes(self):
        return self._num_classes

    def void_label(self):
        return self._void_label

    def ignore_classes(self):
        return self.ignore_classes

max_stars_repo_name: GeniusGaryant/PReMVOS

max_issues_repo_name: GeniusGaryant/PReMVOS

max_forks_repo_name: GeniusGaryant/PReMVOS


text:
#!/usr/bin/env python

'''Correlate measurement time series with a population of objects to find the best match.

Currently only works for Mono-static measurements.

# TODO: Assume a uniform prior distribution over population index, posterior distribution is the probability of what object generated the data. Probability comes from measurement covariance.
'''

import sys
import os
import time
import glob
from tqdm import tqdm

try:
    from mpi4py import MPI
    comm = MPI.COMM_WORLD
except ImportError:
    comm = None

import numpy as np
import scipy
import h5py
import scipy.constants as constants

def residual_distribution_metric(t, r, v, r_ref, v_ref):
    '''Using the simulated and the measured ranges and rage-rates calculate a de-correlation metric.

    :param numpy.ndarray t: Times in seconds corresponding to measurement and object simulated data.
    :param numpy.ndarray r: Measured ranges in meters
    :param numpy.ndarray v: Measured rage-rates in meters per second
    :param numpy.ndarray r_ref: Object simulated ranges in meters
    :param numpy.ndarray v_ref: Object simulated rage-rates in meters per second
    :return: Metric value, smaller values indicate better match.
    :rtype: float
    '''
    residual_r_mu = np.mean(r_ref - r)
    residual_v_mu = np.mean(v_ref - v)

    metric = np.abs(residual_r_mu) + np.abs(residual_v_mu)

    return metric

def generate_measurements(state_ecef, rx_ecef, tx_ecef):

    r_tx = tx_ecef[:,None] - state_ecef[:3,:]
    r_rx = rx_ecef[:,None] - state_ecef[:3,:]

    r_tx_n = np.linalg.norm(r_tx, axis=0)
    r_rx_n = np.linalg.norm(r_rx, axis=0)

    r_sim = r_tx_n + r_rx_n

    v_tx = -np.sum(r_tx*state_ecef[3:,:], axis=0)/r_tx_n
    v_rx = -np.sum(r_rx*state_ecef[3:,:], axis=0)/r_rx_n

    v_sim = v_tx + v_rx

    return r_sim, v_sim

def correlate(
        measurements,
        population,
        metric=residual_distribution_metric,
        metric_reduce=lambda x,y: x+y,
        forward_model=generate_measurements,
        variables=['r','v'],
        n_closest=1,
        profiler=None,
        logger=None,
        MPI=False,
    ):
    '''Given a mono-static measurement of ranges and rage-rates, a radar model and a population: correlate measurements with population.

    :param list data: List of dictionaries that contains measurement data. Contents are described below.
    :param station station: Model of receiver station that performed the measurement.
    :param Population population: Population to correlate against.
    :param function metric: Metric used to correlate measurement and simulated object measurement.
    :param function forward_model: A pointer to a function that takes in the ecef-state, the rx and tx station ecefs and calculates the observed variables return as a tuple.
    :param list variables: The data variables recorded by the system.
    :param Profiler profiler:
    :param logging.Logger logger:
    :param bool MPI: If True use internal parallelization with MPI to calculate correlation. Turn to False to externally parallelize with MPI.

    **Measurement data:**

      Each entry in the input :code:`measurements` list must be a dictionary that contains the following fields:
        * 't': [numpy.ndarray] Times relative epoch in seconds
        * 'epoch': [astropy.Time] epoch for measurements
        * 'tx': [sorts.TX] Pointer to the TX station
        * 'rx': [sorts.RX] Pointer to the RX station

      Then it will contain a entry for each name in the `variables` list. By default this is
        * 'r': [numpy.ndarray] Two-way ranges in meters
        * 'v': [numpy.ndarray] Two-way range-rates in meters per second

    '''

    correlation_data = [None]*len(population)
    match_pop = np.empty((len(population),), dtype=np.float64)

    if MPI and comm is not None:
        step = comm.size
        next_check = comm.rank
    else:
        step = 1
        next_check = 0

    pbars = []
    for pbar_id in range(step):
        pbars.append(tqdm(range(len(population)//step), ncols=100))
    pbar = pbars[next_check]

    for ind, obj in enumerate(population):
        if ind != next_check:
            logger.debug('skipping {}/{}'.format(ind+1, len(population)))
            continue

        pbar.set_description('Correlating object {} '.format(ind))
        pbar.update(1)

        if logger is not None:
            logger.debug('correlating {}/{}'.format(ind+1, len(population)))
            logger.debug(obj)

        correlation_data[ind] = []

        for di, data in enumerate(measurements):
            t = data['t']
            t_sort = t.argsort()
            t = t[t_sort]
            var = tuple(data[x][t_sort] for x in variables)
            tx = data['tx']
            rx = data['rx']
            epoch = data['epoch']

            tx_ecef = tx.ecef.copy()
            tx_ecef_norm = tx_ecef/np.linalg.norm(tx_ecef)
            rx_ecef = rx.ecef.copy()
            rx_ecef_norm = rx_ecef/np.linalg.norm(rx_ecef)

            t_prop = (epoch - obj.epoch).sec + t

            states = obj.get_state(t_prop)

            ref_var = forward_model(states, rx_ecef, tx_ecef)

            match = metric(t, *(var + ref_var))

            if di == 0:
                match_pop[ind] = match
            else:
                match_pop[ind] = metric_reduce(match_pop[ind], match)

            cdat = {f'{x}_ref': ref_var[i].copy()  for i,x in enumerate(variables)}
            cdat['match'] = match

            correlation_data[ind].append(cdat)

        next_check += step

    if MPI and comm is not None and step > 1:

        if comm.rank == 0:
            if logger is not None:
                logger.debug('Receiving all results <barrier>')

            for T in range(1,comm.size):
                for ID in range(T,len(population),comm.size):
                    match_pop[ID] = comm.recv(source=T, tag=ID)
                    if logger is not None:
                        logger.debug('received packet {} from PID{}'.format(ID, T))
        else:
            if logger is not None:
                logger.debug('Distributing all correlation results to process 0 <barrier>')

            for ID in range(comm.rank,len(population),comm.size):
                comm.send(match_pop[ID], dest=0, tag=ID)

        if logger is not None:
            logger.debug('---> Distributing done </barrier>')

    if MPI and comm is not None and step > 1:
        if comm.rank == 0:
            best_matches = np.argsort(match_pop)
        else:
            best_matches = None
        best_matches = comm.bcast(best_matches, root=0)
    else:
        best_matches = np.argsort(match_pop)

    if len(best_matches) > n_closest:
        best_matches = best_matches[:n_closest]

    for pbar in pbars:
        pbar.close()

    best_cdata = [None]*len(best_matches)

    for cind in range(len(best_cdata)):
        if MPI and comm is not None and step > 1:
            source_rank = best_matches[cind] % comm.size
            best_cdata[cind] = comm.bcast(correlation_data[best_matches[cind]], root=source_rank)
        else:
            best_cdata[cind] = correlation_data[best_matches[cind]]

    return best_matches, match_pop[best_matches], best_cdata

ext: py

reason: import numpy,import scipy


text:
from __future__ import absolute_import, division, print_function

import os

import matplotlib.pyplot as plt
from matplotlib.patches import Patch
import networkx as nx

from vivarium.core.experiment import Compartment
from vivarium.core.composition import (
    COMPARTMENT_OUT_DIR,
    simulate_compartment_in_experiment,
    plot_simulation_output,
)
from vivarium.library.make_network import save_network
from vivarium.library.units import units
# processes
from vivarium.data.amino_acids import amino_acids
from vivarium.processes.tree_mass import TreeMass
from vivarium.processes.transcription import Transcription, UNBOUND_RNAP_KEY
from vivarium.processes.translation import Translation, UNBOUND_RIBOSOME_KEY
from vivarium.processes.degradation import RnaDegradation
from vivarium.processes.complexation import Complexation
from vivarium.processes.division_volume import DivisionVolume
from vivarium.data.amino_acids import amino_acids
from vivarium.data.nucleotides import nucleotides
from vivarium.states.chromosome import toy_chromosome_config

NAME = 'gene_expression'

class GeneExpression(Compartment):

    defaults = {
        'global_path': ('global',),
        'initial_mass': 1339.0 * units.fg}

    def __init__(self, config):
        self.config = config
        self.global_path = config.get('global_path', self.defaults['global_path'])
        self.initial_mass = config.get('initial_mass', self.defaults['initial_mass'])

    def generate_processes(self, config):
        transcription = Transcription(config.get('transcription', self.config.get('transcription')))
        translation = Translation(config.get('translation', self.config.get('translation')))
        degradation = RnaDegradation(config.get('degradation', self.config.get('degradation')))
        complexation = Complexation(config.get('complexation', self.config.get('complexation')))
        mass_deriver = TreeMass(config.get('mass_deriver', {
            'initial_mass': config.get('initial_mass', self.initial_mass)}))
        division = DivisionVolume(config)

        return {
            'mass_deriver': mass_deriver,
            'transcription': transcription,
            'translation': translation,
            'degradation': degradation,
            'complexation': complexation,
            'division': division}

    def generate_topology(self, config):
        global_path = config.get('global_path', self.global_path)

        return {
            'mass_deriver': {
                'global': global_path},

            'transcription': {
                'chromosome': ('chromosome',),
                'molecules': ('molecules',),
                'proteins': ('proteins',),
                'transcripts': ('transcripts',),
                'factors': ('concentrations',),
                'global': global_path},

            'translation': {
                'ribosomes': ('ribosomes',),
                'molecules': ('molecules',),
                'transcripts': ('transcripts',),
                'proteins': ('proteins',),
                'concentrations': ('concentrations',),
                'global': global_path},

            'degradation': {
                'transcripts': ('transcripts',),
                'proteins': ('proteins',),
                'molecules': ('molecules',),
                'global': global_path},

            'complexation': {
                'monomers': ('proteins',),
                'complexes': ('proteins',),
                'global': global_path},

            'division': {
                'global': global_path}}

# analysis
def gene_network_plot(data, out_dir, filename='gene_network'):
    '''
    Make a gene network plot from configuration data

    - data (dict):
        {
        'operons': operons, the "genes" in a chromosome config with {operon: [genes list]}
        'templates': promoters, the "promoters" in a chromosome config with {promoter: {sites: [], thresholds: []}}
        'complexes': complexes, stoichiometry from a complexation process.
        }

    '''
    node_size = 400
    node_distance = 30
    edge_weight = 1
    iterations = 1000

    operon_suffix = '_o'
    tf_suffix = '_tf'
    promoter_suffix = '_p'
    gene_suffix = '_g'
    complex_suffix = '_cxn'

    # plotting parameters
    color_legend = {
        'operon': [x/255 for x in [199,164,53]],
        'gene': [x / 255 for x in [181,99,206]],
        'promoter': [x / 255 for x in [110,196,86]],
        'transcription_factor': [x / 255 for x in [222,85,80]],
        'complex': [x / 255 for x in [153, 204, 255]],
    }

    # get data
    operons = data.get('operons', {})
    templates = data.get('templates', {})
    complexes = data.get('complexes', {})

    # make graph from templates
    G = nx.Graph()

    # initialize node lists for graphing
    operon_nodes = set()
    gene_nodes = set()
    promoter_nodes = set()
    tf_nodes = set()
    complex_nodes = set()

    edges = []

    # add operon --> gene connections
    for operon, genes in operons.items():
        operon_name = operon + operon_suffix
        gene_names = [gene + gene_suffix for gene in genes]

        operon_nodes.add(operon_name)
        gene_nodes.update(gene_names)

        G.add_node(operon_name)
        G.add_nodes_from(gene_names)

        # set node attributes
        operon_attrs = {operon_name: {'type': 'operon'}}
        gene_attrs = {gene: {'type': 'gene'} for gene in gene_names}
        nx.set_node_attributes(G, operon_attrs)
        nx.set_node_attributes(G, gene_attrs)

        # add operon --> gene edge
        for gene_name in gene_names:
            edge = (operon_name, gene_name)
            edges.append(edge)
            G.add_edge(operon_name, gene_name)

    # add transcription factor --> promoter --> operon connections
    for promoter, specs in templates.items():
        promoter_name = promoter + promoter_suffix

        promoter_nodes.add(promoter_name)
        G.add_node(promoter_name)

        # set node attributes
        promoter_attrs = {promoter_name: {'type': 'promoter'}}
        nx.set_node_attributes(G, promoter_attrs)

        # get sites and terminators
        promoter_sites = specs['sites']
        terminators = specs['terminators']

        # add transcription factors
        for site in promoter_sites:
            thresholds = site['thresholds']
            for threshold in thresholds:
                tf_name = threshold + tf_suffix

                tf_nodes.add(tf_name)
                G.add_node(tf_name)

                # set node attributes
                tf_attrs = {tf_name: {'type': 'transcription_factor'}}
                nx.set_node_attributes(G, tf_attrs)

                # connect transcription_factor --> promoter
                edge = (tf_name, promoter_name)
                edges.append(edge)
                G.add_edge(tf_name, promoter_name)

        # add gene products
        for site in terminators:
            products = site['products']
            for product in products:
                operon_name = product + operon_suffix

                operon_nodes.add(operon_name)
                G.add_node(operon_name)

                # set node attributes
                operon_attrs = {operon_name: {'type': 'operon'}}
                nx.set_node_attributes(G, operon_attrs)

                # connect promoter --> operon
                edge = (promoter_name, operon_name)
                edges.append(edge)
                G.add_edge(promoter_name, operon_name)

    ## add gene product --> complex
    # first get the sets of complexes and reactants.
    complex_set = set()
    monomer_set = set()
    for complex, stoichiometry in complexes.items():
        complex_list = [mol_id for mol_id, coeff in stoichiometry.items() if coeff>0]
        reactant_list = [mol_id for mol_id, coeff in stoichiometry.items() if coeff<0]

        complex_set.update(complex_list)
        monomer_set.update(reactant_list)

    # complexes are removed from reactants
    for complex in complex_set:
        if complex in monomer_set:
            monomer_set.remove(complex)

    # add nodes and edges
    for complex, stoichiometry in complexes.items():
        complexes = [mol_id for mol_id, coeff in stoichiometry.items() if coeff>0]
        reactants = {mol_id: coeff for mol_id, coeff in stoichiometry.items() if coeff<0}

        assert len(complexes) == 1, 'too many complexes'
        complex = complexes[0]
        complex_name = complex + complex_suffix

        complex_nodes.add(complex_name)
        G.add_node(complex_name)

        # set node attributes
        complex_attrs = {complex_name: {'type': 'complex'}}
        nx.set_node_attributes(G, complex_attrs)

        # make edge
        for reactant, coeff in reactants.items():

            # is this reactant a monomer or a complex?
            if reactant in monomer_set:
                # TODO -- check that it is actually included in the genes?
                reactant_name = reactant + gene_suffix
            elif reactant in complex_set:
                reactant_name = reactant + complex_suffix

            # connect reactant --> complex
            edge = (reactant_name, complex_name)
            edges.append(edge)
            G.add_edge(reactant_name, complex_name)

    # add edges between proteins/complexes and transcription factors that share the same name
    tf_names = [node.replace(tf_suffix, '') for node in tf_nodes]
    for node in gene_nodes:
        gene_name = node.replace(gene_suffix, '')
        if gene_name in tf_names:
            tf_node = gene_name + tf_suffix
            edge = (node, tf_node)
            edges.append(edge)
            G.add_edge(node, tf_node)

    for node in complex_nodes:
        complex_name = node.replace(complex_suffix, '')
        if complex_name in tf_names:
            tf_node = complex_name + tf_suffix
            edge = (node, tf_node)
            edges.append(edge)
            G.add_edge(node, tf_node)

    # separate out the disconnected graphs
    subgraphs = sorted(nx.connected_components(G), key = len, reverse=True)

    # make node labels by removing suffix
    node_labels = {}
    o_labels = {node: node.replace(operon_suffix,'') for node in operon_nodes}
    g_labels = {node: node.replace(gene_suffix, '') for node in gene_nodes}
    p_labels = {node: node.replace(promoter_suffix, '') for node in promoter_nodes}
    tf_labels = {node: node.replace(tf_suffix, '') for node in tf_nodes}
    cxn_labels = {node: node.replace(complex_suffix, '') for node in complex_nodes}
    node_labels.update(o_labels)
    node_labels.update(g_labels)
    node_labels.update(p_labels)
    node_labels.update(tf_labels)
    node_labels.update(cxn_labels)

    # save network for use in gephi
    gephi_nodes = {}
    gephi_edges = [[node1, node2, edge_weight] for (node1, node2) in edges]
    for node_id in operon_nodes:
        gephi_nodes[node_id] = {
            'label': o_labels[node_id],
            'type': 'operon',
            'size': 1}
    for node_id in gene_nodes:
        gephi_nodes[node_id] = {
            'label': g_labels[node_id],
            'type': 'gene',
            'size': 1}
    for node_id in promoter_nodes:
        gephi_nodes[node_id] = {
            'label': p_labels[node_id],
            'type': 'promoter',
            'size': 1}
    for node_id in tf_nodes:
        gephi_nodes[node_id] = {
            'label': tf_labels[node_id],
            'type': 'transcription factor',
            'size': 1}
    for node_id in complex_nodes:
        gephi_nodes[node_id] = {
            'label': cxn_labels[node_id],
            'type': 'complex',
            'size': 1}
    save_network(gephi_nodes, gephi_edges, out_dir)

    # plot
    n_rows = len(list(subgraphs))
    n_cols = 1
    fig = plt.figure(3, figsize=(6*n_cols, 6*n_rows))
    grid = plt.GridSpec(n_rows, n_cols, wspace=0.2, hspace=0.2)

    for idx, subgraph_nodes in enumerate(subgraphs):
        ax = fig.add_subplot(grid[idx, 0])

        subgraph = G.subgraph(list(subgraph_nodes))
        subgraph_node_labels = {node: node_labels[node] for node in subgraph_nodes}

        # get positions
        dist = node_distance / (len(subgraph)**0.5)
        # pos = nx.spring_layout(subgraph, k=dist, iterations=iterations)
        pos = nx.spring_layout(subgraph,
                               k=dist,
                               scale=20,
                               iterations=iterations)

        color_map = []
        for node in subgraph:
            type = subgraph.nodes[node]['type']
            node_color = color_legend[type]
            color_map.append(node_color)

        nx.draw(subgraph, pos, node_size=node_size, node_color=color_map)

        # edges
        # colors = list(range(1,len(edges)+1))
        nx.draw_networkx_edges(subgraph, pos,
                               # edge_color=colors,
                               width=1.5)

        nx.draw_networkx_labels(subgraph, pos,
                                labels=subgraph_node_labels,
                                font_size=8)

        if idx == 0:
            # make legend
            legend_elements = [Patch(facecolor=color, label=name) for name, color in color_legend.items()]
            plt.legend(handles=legend_elements, bbox_to_anchor=(1.5, 1.0))

    # save figure
    fig_path = os.path.join(out_dir, filename)
    plt.figure(3, figsize=(12, 12))
    plt.axis('off')
    plt.savefig(fig_path, bbox_inches='tight')
    plt.close()

def plot_gene_expression_output(timeseries, config, out_dir='out'):

    name = config.get('name', 'gene_expression')
    ports = config.get('ports', {})
    molecules = timeseries[ports['molecules']]
    transcripts = timeseries[ports['transcripts']]
    proteins = timeseries[ports['proteins']]
    time = timeseries['time']

    # make figure and plot
    n_cols = 1
    n_rows = 5
    plt.figure(figsize=(n_cols * 6, n_rows * 1.5))

    # define subplots
    ax1 = plt.subplot(n_rows, n_cols, 1)
    ax2 = plt.subplot(n_rows, n_cols, 2)
    ax3 = plt.subplot(n_rows, n_cols, 3)
    ax4 = plt.subplot(n_rows, n_cols, 4)
    ax5 = plt.subplot(n_rows, n_cols, 5)

    polymerase_ids = [
        UNBOUND_RNAP_KEY,
        UNBOUND_RIBOSOME_KEY]
    amino_acid_ids = list(amino_acids.values())
    nucleotide_ids = list(nucleotides.values())

    # plot polymerases
    for poly_id in polymerase_ids:
        ax1.plot(time, proteins[poly_id], label=poly_id)
    ax1.legend(loc='center left', bbox_to_anchor=(1.0, 0.5))
    ax1.title.set_text('polymerases')

    # plot nucleotides
    for nuc_id in nucleotide_ids:
        ax2.plot(time, molecules[nuc_id], label=nuc_id)
    ax2.legend(loc='center left', bbox_to_anchor=(1.0, 0.5))
    ax2.title.set_text('nucleotides')

    # plot molecules
    for aa_id in amino_acid_ids:
        ax3.plot(time, molecules[aa_id], label=aa_id)
    ax3.legend(loc='center left', bbox_to_anchor=(2.0, 0.5))
    ax3.title.set_text('amino acids')

    # plot transcripts
    for transcript_id, series in transcripts.items():
        ax4.plot(time, series, label=transcript_id)
    ax4.legend(loc='center left', bbox_to_anchor=(1.0, 0.5))
    ax4.title.set_text('transcripts')

    # plot proteins
    for protein_id in sorted(proteins.keys()):
        if protein_id != UNBOUND_RIBOSOME_KEY:
            ax5.plot(time, proteins[protein_id], label=protein_id)
    ax5.legend(loc='center left', bbox_to_anchor=(1.5, 0.5))
    ax5.title.set_text('proteins')

    # adjust axes
    for axis in [ax1, ax2, ax3, ax4, ax5]:
        axis.spines['right'].set_visible(False)
        axis.spines['top'].set_visible(False)

    ax1.set_xticklabels([])
    ax2.set_xticklabels([])
    ax3.set_xticklabels([])
    ax4.set_xticklabels([])
    ax5.set_xlabel('time (s)', fontsize=12)

    # save figure
    fig_path = os.path.join(out_dir, name)
    plt.subplots_adjust(wspace=0.3, hspace=0.5)
    plt.savefig(fig_path, bbox_inches='tight')

# test
def run_gene_expression(total_time=10, out_dir='out'):
    timeseries = test_gene_expression(total_time)
    plot_settings = {
        'name': 'gene_expression',
        'ports': {
            'transcripts': 'transcripts',
            'molecules': 'molecules',
            'proteins': 'proteins'}}
    plot_gene_expression_output(timeseries, plot_settings, out_dir)

    sim_plot_settings = {'max_rows': 25}
    plot_simulation_output(timeseries, sim_plot_settings, out_dir)

def test_gene_expression(total_time=10):
    # load the compartment
    compartment_config = {
        'external_path': ('external',),
        'exchange_path': ('exchange',),
        'global_path': ('global',),
        'agents_path': ('..', '..', 'cells',),
        'transcription': {
            'sequence': toy_chromosome_config['sequence'],
            'templates': toy_chromosome_config['promoters'],
            'genes': toy_chromosome_config['genes'],
            'promoter_affinities': toy_chromosome_config['promoter_affinities'],
            'transcription_factors': ['tfA', 'tfB'],
            'elongation_rate': 10.0}}
    compartment = GeneExpression(compartment_config)

    molecules = {
        nt: 1000
        for nt in nucleotides.values()}
    molecules.update({
        aa: 1000
        for aa in amino_acids.values()})

    proteins = {
        polymerase: 100
        for polymerase in [
                UNBOUND_RNAP_KEY,
                UNBOUND_RIBOSOME_KEY]}

    proteins.update({
        factor: 1
        for factor in [
                'tfA',
                'tfB']})

    # simulate
    settings = {
        'timestep': 1,
        'total_time': total_time,
        'initial_state': {
            'proteins': proteins,
            'molecules': molecules}}
    return simulate_compartment_in_experiment(compartment, settings)

if __name__ == '__main__':
    out_dir = os.path.join(COMPARTMENT_OUT_DIR, NAME)
    if not os.path.exists(out_dir):
        os.makedirs(out_dir)

    run_gene_expression(600, out_dir)

reason: import networkx


text:
#
#       Yaroslav Bershatsky 2017
#
#       dx/dt = t
#
#

import numpy as np
import matplotlib.pyplot as plt

def main():

    t_n = np.arange(0, 20, h)
    n = int(len(t_n))

    x = t_n * t_n / 2

    #print('x\n', x)

    t = t_n

    a_implicit = np.diag(3 * np.ones(n - 1), k = -1) + np.diag(-4 * np.ones(n), k = 0) + np.diag(np.ones(n - 1), k = 1) # implicit solution k = 2
    #a_implicit = np.diag(25/12 * np.ones(n - 1), k = -1) + np.diag(-4 * np.ones(n), k = 0) + np.diag(3 * np.ones(n - 1), k = 1) + np.diag(-4 / 3 * np.ones(n - 2), k = 2) + np.diag(1 / 4 * np.ones(n - 3), k = 3) # implicit solution k = 4

    a_explicit = np.diag(np.ones(n - 1), k = -1) + np.diag( -1 * np.ones(n), k = 0)  # explicit solution k = 1
    #a_explicit = np.diag(np.ones(n - 1), k = -1) + np.diag( -1 * np.ones(n - 1), k = 1)  # explicit solution k = 2
    #a_explicit = np.diag(1 / 3 * np.ones(n - 1), k = -1) + np.diag( 1/2 * np.ones(n), k = 0)+ np.diag(-1 * np.ones(n - 1), k = 1) + np.diag(1/6 * np.ones(n - 2), k = 2)  # explicit solution k = 3

    #print('a_implicit\n', a_implicit)
    #print('a_explicit\n', a_explicit)

    a_implicit_inv = np.linalg.inv(a_implicit)
    a_explicit_inv = np.linalg.inv(a_explicit)

    #print('a_inv\n', a_inv)
    #print('a_inv * a\n', np.dot(a_inv,a))

    x_implicit_n = -2 * h * np.dot(a_implicit_inv, t_n)
    x_explicit_n = -1 * h * np.dot(a_explicit_inv, t)

    #print(t,'\n', b,'\n', t + b)
    #print('x_n\n', x_n)

    fig = plt.figure(facecolor = 'white')
    plt.plot(x_implicit_n, t_n, linewidth = 1, label = 'Implicit solution')
    plt.plot(x_explicit_n, t_n, linewidth = 1, label = 'Explicit solution')

    plt.legend(fontsize = 12)
    plt.grid(True)
    plt.show()
if __name__ == '__main__':
    main()

max_stars_repo_name: BershackyJaroslav/brusselator

max_issues_repo_name: BershackyJaroslav/brusselator

max_forks_repo_name: BershackyJaroslav/brusselator

reason: import numpy


text:
import numpy as _np
import matplotlib.pyplot as _plt
from scipy.signal import decimate as _decimate
from johnspythonlibrary2.Plot import finalizeSubplot as _finalizeSubplot
import xarray as _xr
from deprecated import deprecated as _deprecated

# %% windowing

def hann_window_1D(array, hann_width, plot=False):
	"""

	Example
	-------
	Example 1 ::

		array = _np.random.rand(1000)
		hann_width = 500
		array_hann = hann_window_1D(array, hann_width, plot=True)

	Example 2 ::

		# xarray

		array = _xr.DataArray(_np.random.rand(1000), dims='t', coords=[np.arange(1000)*0.1])
		hann_width = 500
		array_hann = hann_window_1D(array, hann_width, plot=True)

	"""
	if type(array) == _xr.core.dataarray.DataArray:
		xr_array = array.copy()
		array = array.data
		data_xr = True
	else:
		data_xr = False

	from scipy.signal import hann
	if hann_width == array.shape[0]:
		y_hann = hann(hann_width)
	elif hann_width < array.shape[0]:
		y_hann = zero_pad(hann(hann_width), array.shape[0])
	else:
		raise Exception('Array should be longer than hann_width')

	y_out = y_hann * array

	if plot is True:
		fig, ax = _plt.subplots()
		ax.plot(array, label='Original data')
		ax.plot(y_hann, label='Hann window', linestyle='--', linewidth=2)
		ax.plot(y_out, label='Hann applied to data')
		ax.legend()

	if data_xr is True:
		xr_array.data = y_out
		y_out = xr_array

	return y_out

def hann_window_1D_to_2D_data(array, hann_width, axis=0, plot=False):
	"""

	Example
	-------
	Example 1 ::

		array = _np.random.rand(1000, 1200)
		hann_width = 500
		axis = 0
		array_hann = hann_window_1D_to_2D_data(array, hann_width, axis=axis, plot=True)

	Example 2 ::

		array = _np.random.rand(1000, 1200)
		hann_width = 500
		axis = 1
		array_hann = hann_window_1D_to_2D_data(array, hann_width, axis=axis, plot=True)

	Example 3 ::

		# xr example

		array = _xr.DataArray(_np.random.rand(1000, 1200), dims=['x', 'y'], coords=[_np.arange(1000), _np.arange(1200)])
		hann_width = 500
		axis = 1
		array_hann = hann_window_1D_to_2D_data(array, hann_width, axis=axis, plot=True)

	"""
	if type(array) == _xr.core.dataarray.DataArray:
		xr_array = array.copy()
		array = array.data
		data_xr = True
	else:
		data_xr = False

	from scipy.signal import hann
	if hann_width == array.shape[axis]:
		y_hann = hann(hann_width)
	elif hann_width < array.shape[axis]:
		y_hann = zero_pad(hann(hann_width), array.shape[axis])
	else:
		raise Exception('Array should be longer than hann_width')

	if axis == 0:
		y_hann = _np.repeat(y_hann[..., _np.newaxis], array.shape[1], axis=1)
	elif axis == 1:
		y_hann = _np.repeat(y_hann[_np.newaxis, ...], array.shape[0], axis=0)

	y_out = y_hann * array

	if plot is True:
		fig, ax = _plt.subplots(3, sharex=True)
		ax[0].contourf(array, label='Original data')
		ax[1].contourf(y_hann, label='Hann window', linestyle='--', linewidth=2)
		ax[2].contourf(y_out, label='Hann applied to data')

	if data_xr is True:
		xr_array.data = y_out
		y_out = xr_array

	return y_out

def zero_pad(array: _np.ndarray, target_length: int, axis: int = 0):
	"""
	Pads an array along a specified axis with zeros

	Example
	-------
	Example 1 ::

		# a 2D case where np.remainder(pad_size, 2) == 0

		array = _np.arange(5*7).reshape(5,7)
		axis = 1
		target_length = 11
		b = zero_pad(array, target_length=target_length, axis=axis)

	Example 2 ::

		# a 2D case where np.remainder(pad_size, 2) != 0

		array = _np.arange(5*7).reshape(5,7)
		axis = 1
		target_length = 12
		b = zero_pad(array, target_length=target_length, axis=axis)

	Example 3 ::

		# check the 1D case

		array = _np.arange(5) + 1
		axis = 0
		target_length = 10
		b = zero_pad(array, target_length=target_length, axis=axis)

	References
	----------
	 *  https://stackoverflow.com/questions/19349410/how-to-pad-with-zeros-a-tensor-along-some-axis-python
	"""

	pad_size = (target_length - array.shape[axis]) / 2
	if pad_size <= 0:
		raise Exception('Target length must be larger than the array')

	if _np.remainder(pad_size, 1) == 0:
		pad_sizes = (int(pad_size), int(pad_size))
	else:
		pad_sizes = (int(pad_size), int(pad_size+1))

	npad = [(0, 0)] * array.ndim
	npad[axis] = pad_sizes

	return _np.pad(array, pad_width=npad, mode='constant', constant_values=0)

# %% downsampling

@_deprecated(reason='Replaced by downsample().')
def downsample_and_antialiasing(da, downsample_factor=10, plot=False):

	from scipy.signal import decimate #, resample

	t_orig = da.coords[da.dims[0]].data
	t_new = t_orig[::downsample_factor]
	da_new = _xr.DataArray(decimate(da.data, q=downsample_factor, ftype='fir'), dims=da.dims, coords=[t_new])

	if plot is True:
		fig, ax = _plt.subplots()
		da.plot(ax=ax, label='orig')
		da_new.plot(ax=ax, label='downsampled')
		ax.legend()

	return da_new

def downsample(da, downsample_factor=10, antialiasing=True, plot=False):
	"""

	Parameters
	----------
	da : xr.DataArray
		1D dataarray to be downsampled
	downsample_factor : int
		The factor to downsampled by.
	antialiasing : bool
		Downsampling can be done with and without an antialiasing filter. The default is True (i.e. with antialiasing)
	plot : bool
		Optional plot of the results

	Returns
	-------
	da_new : xr.DataArray
		The downsampled dataarray

	Example
	-------
	Example 1 ::

		dt = 1e-6
		t = _np.arange(0, 1e4) * dt
		da = _xr.DataArray(_np.sin(2 * _np.pi * 1e3 * t) + _np.random.rand(len(t)))
		da_downsampled = downsample(da, downsample_factor=20, plot=True)
		da_downsampled = downsample(da, downsample_factor=20, plot=True, antialiasing=False)

	"""
	# downsample time
	t_orig = da.coords[da.dims[0]].data
	t_new = t_orig[::downsample_factor]

	# downsample data (with or without antialiasing)
	if antialiasing is True:
		da_new = _xr.DataArray(_decimate(da.data, q=downsample_factor, ftype='fir'), dims=da.dims, coords=[t_new])
	else:
		da_new = _xr.DataArray(da.data[::downsample_factor], dims=da.dims, coords=[t_new])

	# optional plot of results
	if plot is True:
		fig, ax = _plt.subplots()
		da.plot(ax=ax, label='orig')
		da_new.plot(ax=ax, label='downsampled')
		ax.legend()

	return da_new

# %% low and highpass filtering

def filtfilt(	da,
				cornerFreq,
				filterType='low',
				filterOrder=1,
				plot=False,
				axis='t',
				time=None):
	"""
	Forward-backwards filter using a butterworth filter

	Parameters
	----------
	da : xarray.core.dataarray.DataArray
		Signal.  Index is time with units in seconds.
	cornerFreq : float or numpy.array of floats
		Filter's corner frequency
		If bandpass or bandstop, cornerFreq is a numpy.array with two elements
	filterType : str
		* 'low' - low-pass filter
		* 'high' - high-pass filter
		* 'bandpass'
		* 'bandstop'
	filterOrder : int
		Order of the butterworth filter.  8 is default.
	plot : bool
		Optional plot
	time : numpy array, optional
		if da is a numpy array, then time needs to be provided

	Returns
	-------
	dfOut : xarray.core.dataarray.DataArray
		Filtered signal.  Index is time with units in seconds.

	Notes
	-----
	* Be careful with high filterOrders.  I've seen them produce nonsense results.  I recommend starting low and then turning up the order

	Examples
	--------
	Example 1::

		import numpy as np
		import xarray as xr
		t = np.linspace(0, 1.0, 2001)
		dt=t[1]-t[0]
		fs=1.0/dt
		xlow = np.sin(2 * np.pi * 5 * t)
		xmid = np.sin(2 * np.pi * 50 * t)
		xhigh = np.sin(2 * np.pi * 500 * t)
		x = xlow + xhigh + xmid
		da=xr.DataArray(x,dims=['t'],
				  coords={'t':t})
		dfOut=filtfilt( 	da,
							cornerFreq=30.0,
							filterType='low',
							filterOrder=8,
							plot=True)
		dfOut=filtfilt( 	da,
							cornerFreq=200.0,
							filterType='high',
							filterOrder=8,
							plot=True)
		dfOut=filtfilt( 	da,
							cornerFreq=np.array([25,150]),
							filterType='bandpass',
							filterOrder=4,
							plot=True)
		dfOut=filtfilt( 	x,
							cornerFreq=np.array([25,150]),
							filterType='bandpass',
							filterOrder=4,
							plot=True,
							time=t)

 	Example 2::

		t=_np.arange(0,100e-3-4e-6,2e-6)
		fStart=2e2
		fStop=2.0e3
		y1=_chirp(t,[10e-3,90e-3],[fStart,fStop])#[1e-3,19.46e-3]
		da=_xr.DataArray(y1,
 					  dims=['t'],
					   coords={'t':t})

		filterOrder=2
		cornerFreq=1000.0
		daHP=filtfilt(da,cornerFreq=cornerFreq,filterType='high',plot=False,filterOrder=filterOrder)
		daLP=filtfilt(da,cornerFreq=cornerFreq,filterType='low',plot=False,filterOrder=filterOrder)

		fig,ax=_plt.subplots()
		da.plot(ax=ax,label='Original')
		daHP.plot(ax=ax,label='Highpass')
		daLP.plot(ax=ax,label='Lowpass')
		_finalizeSubplot(ax,
					   xlabel='Time',
					   ylabel='Signal amplitude',)

 	Example 3::

		import numpy as _np
		import xarray as _xr

		t=_np.arange(0,100e-3-4e-6,2e-6)
		y=jpl2.Process.SigGen.gaussianNoise(t.shape,plot=False)
		da=_xr.DataArray(y,
 					  dims=['t'],
					   coords={'t':t})

		daOut=filtfilt( 	da,
										cornerFreq=np.array([1e5,1.5e5]),
										filterType='bandpass',
										filterOrder=4,
										plot=False)
		jpl2.Process.Spectral.fft(daOut,plot=True,
									trimNegFreqs=True)

	References
	----------
	* https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.signal.filtfilt.html

	"""
	if type(da) not in [_xr.core.dataarray.DataArray]:
		if type(da) in [_np.ndarray]:
			da=_xr.DataArray(	da,
								dims=['t'],
								coords={'t':time})
		else:
			raise Exception('Invalid input type')
	if 'time' in da.dims:
		da = da.rename({'time': 't'})
	if 't' not in da.dims:
		raise Exception('time not present')

	from scipy.signal import butter, filtfilt

	# construct butterworth filter
	samplingFreq = float(1.0 / (da.t[1] - da.t[0]))
	Wn = _np.array(cornerFreq).astype(float) / samplingFreq * 2  # I don't know why this factor of 2 needs to be here
	b, a = butter(	filterOrder,
					Wn=Wn,
					btype=filterType,
					analog=False,
					)

	# perform forwards-backwards filter
	daOut = _xr.DataArray(	filtfilt(	b, a, da.values.reshape(-1),
										),
							   dims='t',
							   coords={'t': da.t})

	if plot is True:
		fig, ax = _plt.subplots()
		da.plot(ax=ax, label='original')
		daOut.plot(ax=ax, label='filtered')
		_finalizeSubplot(	ax,
							xlabel='Time',
							ylabel='Signal amplitude')

	return daOut

def boxcar_convolution_filter(	da,
								width_in_time,
								filter_type='low',
								mode='reflect',
								plot=False):
	from scipy.ndimage import uniform_filter1d

	dt = _np.mean(da.t[1:].data - da.t[:-1].data)
	width = _np.round(width_in_time / dt).astype(int)
	da_lp = _xr.DataArray(uniform_filter1d(da.data, size=width, mode=mode), dims=da.dims, coords=da.coords)
	da_hp = da - da_lp.data

	if plot is True:
		fig, ax = _plt.subplots()
		da.plot(ax=ax, label='Original')
		da_lp.plot(ax=ax, label='Lowpass')
		da_hp.plot(ax=ax, label='Highpass')
		ax.legend()

	if filter_type == 'high':
		return da_hp
	else:
		return da_lp

def gaussianFilter(	da,
				    timeFWHM,
					filterType='high',
					plot=False):
	"""
	Low and highpass filters using scipy's gaussian convolution filter

	Parameters
	----------
	da  :  xarray.core.dataarray.DataArray
		unfiltered (raw) signal
	timeFWHM : float
		full width at half maximum of the gaussian with units in time.  this
		effectively sets the corner frequency of the filter
	filterType : str
		'high' - high-pass filter
		'low' - low-pass filter
	plot : bool
		plots the results
	plotGaussian : bool
		plots the gaussian distribution used for the filter

	Returns
	-------
	xarray.core.dataarray.DataArray
		filtered signal

	Examples
	--------
	Example 1::

		import xarray as xr
		t=_np.arange(0,100e-3-4e-6,2e-6)
		fStart=2e2
		fStop=2.0e3
		y1=_chirp(t,[10e-3,90e-3],[fStart,fStop])#[1e-3,19.46e-3]
		da=xr.DataArray(y1,
					  dims=['t'],
					  coords={'t':t})

		fwhm=0.5e-3
		daHP=gaussian(da,fwhm,'high',plot=False)
		daLP=gaussian(da,fwhm,'low',plot=False)

		fig,ax=_plt.subplots()
		ax.plot(da,label='Original')
		ax.plot(daHP,label='Highpass')
		ax.plot(daLP,label='Lowpass')
		_finalizeSubplot(ax,
					   xlabel='Time',
					   ylabel='Signal amplitude',)

	References
	----------
	* https://en.wikipedia.org/wiki/Full_width_at_half_maximum
	* https://docs.scipy.org/doc/scipy-0.19.0/reference/generated/scipy.signal.gaussian.html
	* https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.gaussian_filter1d.html

	"""

	# check input
	if type(da) not in [_xr.core.dataarray.DataArray]:
		raise Exception('Invalid input type')
	if 'time' in da.dims:
		da = da.rename({'time':'t'})
	if 't' not in da.dims:
		raise Exception('time not present')

	# import
	from scipy.ndimage import gaussian_filter1d

	# convert FWHM to standard deviation
	dt = float(da.t[1] - da.t[0])
	def fwhmToGaussFilterStd(fwhm, dt=dt):
		std = 1.0 / _np.sqrt(8 * _np.log(2)) * fwhm / dt
		return std
	std = fwhmToGaussFilterStd(timeFWHM, dt)

	# perform gaussian filter
	daFiltered=_xr.DataArray(	gaussian_filter1d(	da,
													std,
													axis=0,
													mode='nearest'),
								dims=['t'],
								coords={'t': da.t})

	# optional plot of results
	if plot is True:
		fig, ax = _plt.subplots()
		da.plot(ax=ax,label='Raw')
		daFiltered.plot(ax=ax, label='Low-pass')
		(da - daFiltered).plot(ax=ax, label='High-pass')
		_finalizeSubplot(	ax,
							xlabel='Time',
						    ylabel='Signal amplitude',)

	if filterType == 'low':
		return daFiltered
	elif filterType == 'high':
		return da-daFiltered
	else:
		raise Exception('Bad filter type')

def butterworth(	da,
					cornerFreq,
					filterType='low',
					filterOrder=1,
					plot=False,
					time=None):
	"""
	Apply a butterworth filter to an input signal
	Introduces a frequency-dependent phase-shift on the output.
	If no phase-shift is desired, use filtfilt() instead.

	Parameters
	----------
	da : xarray.core.dataarray.DataArray
		Signal.  Index is time with units in seconds.
	cornerFreq : float or numpy.array of floats
		Filter's corner frequency
		If bandpass or bandstop, cornerFreq is a numpy.array with two elements
	filterType : str
		* 'low' - low-pass filter
		* 'high' - high-pass filter
		* 'bandpass'
		* 'bandstop'
	filterOrder : int
		Order of the butterworth filter.  8 is default.
	plot : bool
		Optional plot
	time : numpy array, optional
		if da is a numpy array, then time needs to be provided

	Returns
	-------
	daOut : xarray.core.dataarray.DataArray
		Filtered signal.  Index is time with units in seconds.

	Examples
	--------
 	Example 2::

		t=_np.arange(0,100e-3-4e-6,2e-6)
		fStart=2e2
		fStop=2.0e3
		from johnspythonlibrary2.Process.SigGen import chirp as _chirp
		y1=_chirp(t,[10e-3,90e-3],[fStart,fStop])#[1e-3,19.46e-3]
		da=_xr.DataArray(y1,
 					  dims=['t'],
					   coords={'t':t})

		filterOrder=2
		cornerFreq=600.0
		daHP=butterworth(da,cornerFreq=cornerFreq,filterType='high',plot=False,filterOrder=filterOrder)
		daLP=butterworth(da,cornerFreq=cornerFreq,filterType='low',plot=False,filterOrder=filterOrder)

		fig,ax=_plt.subplots()
		da.plot(ax=ax,label='Original')
		daHP.plot(ax=ax,label='Highpass')
		daLP.plot(ax=ax,label='Lowpass')
		_finalizeSubplot(ax,
					   xlabel='Time',
					   ylabel='Signal amplitude',)

 	Example 3::

		import numpy as _np
		import xarray as _xr

		t=_np.arange(0,100e-3-4e-6,2e-6)
		from johnspythonlibrary2.Process.SigGen import gaussianNoise
		from johnspythonlibrary2.Process.Spectral import fft
		y=gaussianNoise(t.shape,plot=False)
		da=_xr.DataArray( 	y,
						    dims=['t'],
							coords={'t':t})

		daOut=butterworth( 	da,
							cornerFreq=_np.array([1e5,1.5e5]),
							filterType='bandpass',
							filterOrder=4,
							plot=False)
		fft(	daOut,
				plot=True,
				trimNegFreqs=True)

	"""

	# check input signal
	if type(da) not in [_xr.core.dataarray.DataArray]:
		if type(da) in [_np.ndarray]:
 			da=_xr.DataArray(	da,
								dims=['t'],
								coords={'t': time})
		else:
 			raise Exception('Invalid input type')
	if 'time' in da.dims:
		da = da.rename({'time': 't'})
	if 't' not in da.dims:
		raise Exception('time not present')

	from scipy.signal import butter, lfilter

 	# construct butterworth filter
	samplingFreq = float(1.0 / (da.t[1] - da.t[0]))
	Wn = _np.array(cornerFreq).astype(float) / samplingFreq * 2  # I don't know why this factor of 2 needs to be here
	b, a = butter(	filterOrder,
 					Wn=Wn,
 					btype=filterType,
 					analog=False,
 					)

	# perform forwards-backwards filter
	daOut = _xr.DataArray(	lfilter(b, a, da.values.reshape(-1)),
							dims='t',
							coords={'t': da.t})

	if plot is True:
		fig ,ax = _plt.subplots()
		da.plot(ax=ax, label='original')
		daOut.plot(ax=ax, label='filtered')
		_finalizeSubplot(	ax,
							xlabel='Time',
							ylabel='Signal amplitude',)

	return daOut

reason: import numpy,from scipy


text:
from numpy import all, asarray, array, where, exp
from pandas import DataFrame
from skimage.filters import gaussian_filter
from skimage.feature import peak_local_max
from scipy.optimize import curve_fit
from scipy.stats import multivariate_normal
import matplotlib.pyplot as plt

def compute(im, options):
    beads, maxima, centers, smoothed = getCenters(im, options)
    return [getPSF(x, options) for x in beads], beads, maxima, centers, smoothed

def inside(shape, center, window):
    """
    Returns boolean if a center and its window is fully contained
    within the shape of the image on all three axes
    """
    return all([(center[i]-window[i] >= 0) & (center[i]+window[i] <= shape[i]) for i in range(0,3)])

def volume(im, center, window):
    if inside(im.shape, center, window):
        volume = im[(center[0]-window[0]):(center[0]+window[0]), (center[1]-window[1]):(center[1]+window[1]), (center[2]-window[2]):(center[2]+window[2])]
        volume = volume.astype('float64')
        baseline = volume[[0,-1],[0,-1],[0,-1]].mean()
        volume = volume - baseline
        volume = volume/volume.max()
        return volume

def findBeads(im, window, thresh):
    smoothed = gaussian_filter(im, 1, output=None, mode='nearest', cval=0, multichannel=None)
    centers = peak_local_max(smoothed, min_distance=3, threshold_rel=thresh, exclude_border=True)
    return centers, smoothed.max(axis=0)

def keepBeads(im, window, centers, options):
    centersM = asarray([[x[0]/options['pxPerUmAx'], x[1]/options['pxPerUmLat'], x[2]/options['pxPerUmLat']] for x in centers])
    centerDists = [nearest(x,centersM) for x in centersM]
    keep = where([x>3 for x in centerDists])
    centers = centers[keep[0],:]
    keep = where([inside(im.shape, x, window) for x in centers])
    return centers[keep[0],:]

def getCenters(im, options):
    window = [options['windowUm'][0]*options['pxPerUmAx'], options['windowUm'][1]*options['pxPerUmLat'], options['windowUm'][2]*options['pxPerUmLat']]
    window = [round(x) for x in window]
    centers, smoothed = findBeads(im, window, options['thresh'])
    centers = keepBeads(im, window, centers, options)
    beads = [volume(im, x, window) for x in centers]
    maxima = [im[x[0], x[1], x[2]] for x in centers]
    return beads, maxima, centers, smoothed

def getPSF(bead, options):
    latProfile, axProfile = getSlices(bead)
    latFit = fit(latProfile,options['pxPerUmLat'])
    axFit = fit(axProfile,options['pxPerUmAx'])
    data = DataFrame([latFit[3], axFit[3]],index = ['FWHMlat', 'FWHMax']).T
    return data, latFit, axFit

def getSlices(average):
    latProfile = (average.mean(axis=0).mean(axis=1) + average.mean(axis=0).mean(axis=1))/2
    axProfile = (average.mean(axis=1).mean(axis=1) + average.mean(axis=2).mean(axis=1))/2
    return latProfile, axProfile

def fit(yRaw,scale):
    y = yRaw - (yRaw[0]+yRaw[-1])/2
    x = (array(range(y.shape[0])) - y.shape[0]/2)
    x = (array(range(y.shape[0])) - y.shape[0]/2)
    popt, pcov = curve_fit(gauss, x, y, p0 = [1, 0, 1, 0])
    FWHM = 2.355*popt[2]/scale
    yFit = gauss(x, *popt)
    return x, y, yFit, FWHM

def plotPSF(x,y,yFit,FWHM,scale,Max):
    plt.plot(x.astype(float)/scale,yFit/yFit.max(), lw=2);
    plt.plot(x.astype(float)/scale,y/yFit.max(),'ok');
    plt.xlim([-x.shape[0]/2/scale, x.shape[0]/2/scale])
    plt.ylim([0, 1.1])
    plt.xlabel('Distance (um)')
    plt.ylabel('Norm. intensity')
    plt.annotate('FWHM %.2f um' % FWHM,xy=(x.shape[0]/4/scale, .6), size=14)
    plt.annotate('Brightness %.2f' % Max,xy=(x.shape[0]/4/scale, .5), size=14)

def plotAvg(i):
    plt.figure(figsize=(5,5));
    plt.imshow(average[i], vmin=0, vmax=.9);
    if i==average.shape[0]/2:
        plt.plot(average.shape[1]/2, average.shape[2]/2, 'r.', ms=10);
    plt.xlim([0, average.shape[1]])
    plt.ylim([average.shape[2], 0])
    plt.axis('off');

def plotAvg(i):
    plt.figure(figsize=(5,5));
    plt.imshow(average[i], vmin=0, vmax=.9);
    if i==average.shape[0]/2:
        plt.plot(average.shape[1]/2, average.shape[2]/2, 'r.', ms=10);
    plt.xlim([0, average.shape[1]])
    plt.ylim([average.shape[2], 0])
    plt.axis('off');

def dist(x,y):
    return ((x - y)**2)[1:].sum()**(.5)

def nearest(x,centers):
    z = [dist(x,y) for y in centers if not (x == y).all()]
    return abs(array(z)).min(axis=0)

def gauss(x, a, mu, sigma, b):
    return a*exp(-(x-mu)**2/(2*sigma**2))+b

max_stars_repo_name: sofroniewn/psf

max_issues_repo_name: sofroniewn/psf

max_forks_repo_name: sofroniewn/psf

reason: from numpy,from scipy


text:
# ==============================================================================
"""Additional Numpy methods. Big mess of many things!
"""
import numpy as np

# =========================================================================== #
# Numpy implementations of SSD boxes functions.
# =========================================================================== #
def ssd_bboxes_decode(feat_localizations,
                      anchor_bboxes,
                      prior_scaling=[0.1, 0.1, 0.2, 0.2]):
    """Compute the relative bounding boxes from the layer features and
    reference anchor bounding boxes.

    Return:
      numpy array Nx4: ymin, xmin, ymax, xmax
    """
    # Reshape for easier broadcasting.
    l_shape = feat_localizations.shape

    # feat_localizations = np.reshape(feat_localizations,
    #                                 (-1, l_shape[-1]))
    feat_localizations = np.reshape(feat_localizations,
                                    (-1, l_shape[-2], l_shape[-1]))
    anchor_xmin, anchor_ymin, anchor_xmax, anchor_ymax = anchor_bboxes

    xref = (anchor_xmin + anchor_xmax) /2.
    yref = (anchor_ymin + anchor_ymax) /2.

    href = anchor_ymax - anchor_ymin
    wref = anchor_xmax - anchor_xmin

    decode_bbox_center_x = prior_scaling[0] * feat_localizations[:, :, 0] * wref + xref
    decode_bbox_center_y = prior_scaling[1] * feat_localizations[:, :, 1] * href + yref
    decode_bbox_width = np.exp(feat_localizations[:, :, 2] * prior_scaling[2]) * wref
    decode_bbox_height = np.exp(feat_localizations[: , :, 3] * prior_scaling[3]) * href

    xmin = decode_bbox_center_x - decode_bbox_width / 2.
    ymin = decode_bbox_center_y - decode_bbox_height / 2.
    xmax = decode_bbox_center_x + decode_bbox_width / 2.
    ymax = decode_bbox_center_y + decode_bbox_height / 2.

    x1 = prior_scaling[0] * feat_localizations[:, :, 4] * wref + anchor_xmin
    y1 = prior_scaling[1] * feat_localizations[:, :, 5] * href + anchor_ymin
    x2 = prior_scaling[0] * feat_localizations[:, :, 6] * wref + anchor_xmax
    y2 = prior_scaling[1] * feat_localizations[:, :, 7] * href + anchor_ymin

    x3 = prior_scaling[0] * feat_localizations[:, :, 8] * wref + anchor_xmax
    y3 = prior_scaling[1] * feat_localizations[:, :, 9] * href + anchor_ymax
    x4 = prior_scaling[0] * feat_localizations[:, :, 10] * wref + anchor_xmin
    y4 = prior_scaling[1] * feat_localizations[:, :, 11] * href + anchor_ymax

    # bboxes: ymin, xmin, xmax, ymax.
    bboxes = np.zeros_like(feat_localizations)

    bboxes[:,:,0] = xmin
    bboxes[:,:,1] = ymin
    bboxes[:,:,2] = xmax
    bboxes[:,:,3] = ymax
    bboxes[:, :, 4] = x1
    bboxes[:, :, 5] = x2
    bboxes[:, :, 6] = x3
    bboxes[:, :, 7] = x4
    bboxes[:, :, 8] = y1
    bboxes[:, :, 9] = y2
    bboxes[:, :, 10] = y3
    bboxes[:, :, 11] = y4

    # Back to original shape.
    bboxes = np.reshape(bboxes, l_shape)
    return bboxes

def ssd_bboxes_select(predictions_net,
                      localizations_net,
                      anchors_net,
                      select_threshold=0.5,
                      img_shape=(384, 384),
                      num_classes=2,
                      decode=True):
    """Extract classes, scores and bounding boxes from network output layers.

    Return:
      classes, scores, bboxes: Numpy arrays...
    """
    l_classes = []
    l_scores = []
    l_bboxes = []
    for i in range(len(predictions_net)):
        classes, scores, bboxes = ssd_bboxes_select_layer(
            predictions_net[i], localizations_net[i], anchors_net[i],
            select_threshold, img_shape, num_classes, decode)
        l_classes.append(classes)
        l_scores.append(scores)
        l_bboxes.append(bboxes)
        # Debug information.
        # l_layers.append(i)
        # l_idxes.append((i, idxes))

    classes = np.concatenate(l_classes, 0)
    scores = np.concatenate(l_scores, 0)
    bboxes = np.concatenate(l_bboxes, 0)
    return classes, scores, bboxes

def ssd_bboxes_select_layer(predictions_layer,
                            localizations_layer,
                            anchors_layer,
                            select_threshold=0.2,
                            img_shape=(384, 384),
                            num_classes=2,
                            decode=True):
    """Extract classes, scores and bounding boxes from features in one layer.

    Return:
      classes, scores, bboxes: Numpy arrays...
    """
    # First decode localizations features if necessary.
    if decode:
        localizations_layer = ssd_bboxes_decode(localizations_layer, anchors_layer)

    # Reshape features to: Batches x N x N_labels | 4.
    p_shape = predictions_layer.shape
    batch_size = p_shape[0] if len(p_shape) == 5 else 1
    predictions_layer = np.reshape(predictions_layer,
                                   (batch_size, -1, p_shape[-1]))
    l_shape = localizations_layer.shape
    localizations_layer = np.reshape(localizations_layer,
                                     (batch_size, -1, l_shape[-1]))

    # Boxes selection: use threshold or score > no-label criteria.
    if select_threshold is None or select_threshold == 0:
        # Class prediction and scores: assign 0. to 0-class
        classes = np.argmax(predictions_layer, axis=2)
        scores = np.argmax(predictions_layer, axis=2)
        mask = (classes > 0)
        classes = classes[mask]
        scores = scores[mask]
        bboxes = localizations_layer[mask]
    else:#two preditcions
        sub_predictions = predictions_layer[:, :, 1:]
        idxes = np.where(sub_predictions > select_threshold)
        classes = idxes[-1]+1
        scores = sub_predictions[idxes]
        bboxes = localizations_layer[idxes[:-1]]

    return classes, scores, bboxes

# =========================================================================== #
# Common functions for bboxes handling and selection.
# =========================================================================== #
def bboxes_sort(classes, scores, bboxes, top_k=400):
    """Sort bounding boxes by decreasing order and keep only the top_k
    """
    # if priority_inside:
    #     inside = (bboxes[:, 0] > margin) & (bboxes[:, 1] > margin) & \
    #         (bboxes[:, 2] < 1-margin) & (bboxes[:, 3] < 1-margin)
    #     idxes = np.argsort(-scores)
    #     inside = inside[idxes]
    #     idxes = np.concatenate([idxes[inside], idxes[~inside]])
    idxes = np.argsort(-scores)
    classes = classes[idxes][:top_k]
    scores = scores[idxes][:top_k]
    bboxes = bboxes[idxes][:top_k]
    return classes, scores, bboxes

def bboxes_clip(bbox_ref, bboxes):
    """Clip bounding boxes with respect to reference bbox.
    """
    bboxes = np.copy(bboxes)
    bboxes = np.transpose(bboxes)
    bbox_ref = np.transpose(bbox_ref)
    bboxes[0] = np.maximum(bboxes[0], bbox_ref[0])
    bboxes[1] = np.maximum(bboxes[1], bbox_ref[1])
    bboxes[2] = np.minimum(bboxes[2], bbox_ref[2])
    bboxes[3] = np.minimum(bboxes[3], bbox_ref[3])

    bboxes[4] = np.maximum(bboxes[4], bbox_ref[4])
    bboxes[5] = np.maximum(bboxes[5], bbox_ref[5])
    bboxes[6] = np.maximum(bboxes[6], bbox_ref[6])
    bboxes[7] = np.maximum(bboxes[7], bbox_ref[7])
    bboxes[8] = np.maximum(bboxes[8], bbox_ref[8])
    bboxes[9] = np.maximum(bboxes[9], bbox_ref[9])
    bboxes[10] = np.maximum(bboxes[10], bbox_ref[10])
    bboxes[11] = np.maximum(bboxes[11], bbox_ref[11])

    bboxes = np.transpose(bboxes)
    return bboxes

def bboxes_resize(bbox_ref, bboxes):
    """Resize bounding boxes based on a reference bounding box,
    assuming that the latter is [0, 0, 1, 1] after transform.
    """
    bboxes = np.copy(bboxes)
    # Translate.
    bboxes[:, 0] -= bbox_ref[0]
    bboxes[:, 1] -= bbox_ref[1]
    bboxes[:, 2] -= bbox_ref[0]
    bboxes[:, 3] -= bbox_ref[1]
    # Resize.
    resize = [bbox_ref[2] - bbox_ref[0], bbox_ref[3] - bbox_ref[1]]
    bboxes[:, 0] /= resize[0]
    bboxes[:, 1] /= resize[1]
    bboxes[:, 2] /= resize[0]
    bboxes[:, 3] /= resize[1]
    return bboxes

def bboxes_jaccard(bboxes1, bboxes2):
    """Computing jaccard index between bboxes1 and bboxes2.
    Note: bboxes1 and bboxes2 can be multi-dimensional, but should broacastable.
    """
    bboxes1 = np.transpose(bboxes1)
    bboxes2 = np.transpose(bboxes2)
    # Intersection bbox and volume.
    int_ymin = np.maximum(bboxes1[0], bboxes2[0])
    int_xmin = np.maximum(bboxes1[1], bboxes2[1])
    int_ymax = np.minimum(bboxes1[2], bboxes2[2])
    int_xmax = np.minimum(bboxes1[3], bboxes2[3])

    int_h = np.maximum(int_ymax - int_ymin, 0.)
    int_w = np.maximum(int_xmax - int_xmin, 0.)
    int_vol = int_h * int_w
    # Union volume.
    vol1 = (bboxes1[2] - bboxes1[0]) * (bboxes1[3] - bboxes1[1])
    vol2 = (bboxes2[2] - bboxes2[0]) * (bboxes2[3] - bboxes2[1])
    jaccard = int_vol / (vol1 + vol2 - int_vol)
    return jaccard

def bboxes_intersection(bboxes_ref, bboxes2):
    """Computing jaccard index between bboxes1 and bboxes2.
    Note: bboxes1 and bboxes2 can be multi-dimensional, but should broacastable.
    """
    bboxes_ref = np.transpose(bboxes_ref)
    bboxes2 = np.transpose(bboxes2)
    # Intersection bbox and volume.
    int_ymin = np.maximum(bboxes_ref[0], bboxes2[0])
    int_xmin = np.maximum(bboxes_ref[1], bboxes2[1])
    int_ymax = np.minimum(bboxes_ref[2], bboxes2[2])
    int_xmax = np.minimum(bboxes_ref[3], bboxes2[3])

    int_h = np.maximum(int_ymax - int_ymin, 0.)
    int_w = np.maximum(int_xmax - int_xmin, 0.)
    int_vol = int_h * int_w
    # Union volume.
    vol = (bboxes_ref[2] - bboxes_ref[0]) * (bboxes_ref[3] - bboxes_ref[1])
    score = int_vol / vol
    return score

def bboxes_nms(classes, scores, bboxes, nms_threshold=0.45):
    """Apply non-maximum selection to bounding boxes.
    """
    keep_bboxes = np.ones(scores.shape, dtype=np.bool)
    for i in range(scores.size-1):
        if keep_bboxes[i]:
            # Computer overlap with bboxes which are following.
            overlap = bboxes_jaccard(bboxes[i], bboxes[(i+1):])
            # Overlap threshold for keeping + checking part of the same class
            keep_overlap = np.logical_or(overlap < nms_threshold, classes[(i+1):] != classes[i])
            keep_bboxes[(i+1):] = np.logical_and(keep_bboxes[(i+1):], keep_overlap)

    idxes = np.where(keep_bboxes)
    return classes[idxes], scores[idxes], bboxes[idxes]

def bboxes_nms_fast(classes, scores, bboxes, threshold=0.45):
    """Apply non-maximum selection to bounding boxes.
    """
    pass


text:
from dataclasses import dataclass
from typing import Dict, List, Optional, Set

import numpy as np
from nuplan.common.actor_state.vehicle_parameters import get_pacifica_parameters
from nuplan.planning.metrics.abstract_metric import AbstractMetricBuilder
from nuplan.planning.metrics.metric_result import MetricStatistics, MetricStatisticsType, Statistic
from nuplan.planning.metrics.utils.geometry import signed_lateral_distance, signed_longitudinal_distance
from nuplan.planning.simulation.history.simulation_history import SimulationHistory
from nuplan.planning.simulation.observation.observation_type import Detections
from nuplan.planning.simulation.observation.smart_agents.idm_agents.utils import box3d_to_polygon

@dataclass
class EgoToAgentDistances:
    """
    Class to keep track of the history of projected distances from ego to an agent.
    It also contains the length of the agent.
    """
    agent_length: float  # Length of agent [m]
    longitudinal_distances: List[float]  # Longitudinal distance from ego to the agent [m]
    lateral_distances: List[float]  # Lateral distance from ego to the agent [m]

class ClearanceFromStaticAgentsStatistics(AbstractMetricBuilder):

    def __init__(self, name: str, category: str, lateral_distance_threshold: float) -> None:
        """
        Metric on clearance while passing static vehicles.
        :param name: Metric name.
        :param category: Metric category.
        :param lateral_distance_threshold: Agents laterally further away than this threshold are not considered
        """

        self._name = name
        self._category = category
        self._lateral_distance_threshold = lateral_distance_threshold
        self._ego_half_length = get_pacifica_parameters().half_length

    @property
    def name(self) -> str:
        """
        Returns the metric name.
        :return: the metric name.
        """

        return self._name

    @property
    def category(self) -> str:
        """
        Returns the metric category.
        :return: the metric category.
        """

        return self._category

    def compute(self, history: SimulationHistory) -> List[MetricStatistics]:
        """
        Returns the estimated metric.
        :param history: History from a simulation engine.
        :return: the estimated metric.
        """

        # Compute projected distances
        agents_distances = self._extract_agent_projected_distances(history)

        clearances_during_passing = self._extract_passing_clearances(agents_distances)

        if not clearances_during_passing:
            return []

        statistics = {MetricStatisticsType.MAX: Statistic(name="max_clearance_overtaking_static_agent", unit="meters",
                                                          value=np.amax(clearances_during_passing)),
                      MetricStatisticsType.MIN: Statistic(name="min_clearance_overtaking_static_agent", unit="meters",
                                                          value=np.amin(clearances_during_passing)),
                      MetricStatisticsType.P90: Statistic(name="p90_clearance_overtaking_static_agent", unit="meters",
                                                          value=np.percentile(np.abs(clearances_during_passing), 90)),
                      }

        result = MetricStatistics(metric_computator=self.name,
                                  name="clearance_from_static_agents_statistics",
                                  statistics=statistics,
                                  time_series=None,
                                  metric_category=self.category)

        return [result]

    def get_overtake_start_idx(self, longitudinal_dist: List[float], idx_overtake: int,
                               critical_dist_abs: float) -> int:
        """
        Finds the index of the element which represents the start of the overtake
        :param longitudinal_dist: longitudinal distances
        :param idx_overtake: index of the distance closest to zero
        :param critical_dist_abs: critical distance which represent start of overtake
        :return: index of the start of overtake
        """
        offset = self._get_overtake_edge(longitudinal_dist[idx_overtake::-1], critical_dist_abs)
        return idx_overtake - offset if offset is not None else 0

    def get_overtake_end_idx(self, longitudinal_dist: List[float], idx_overtake: int, critical_dist_abs: float) -> int:
        """
        Finds the index of the element which represents the end of the overtake
        :param longitudinal_dist: longitudinal distances
        :param idx_overtake: index of the distance closest to zero
        :param critical_dist_abs: critical distance which represent end of overtake
        :return: index of the end of overtake
        """
        offset = self._get_overtake_edge(longitudinal_dist[idx_overtake:], critical_dist_abs)
        return idx_overtake + offset if offset is not None else -1

    @staticmethod
    def _get_overtake_edge(distances: List[float], critical_distance: float) -> Optional[int]:
        """
        Finds the index of the first element which exceeds the given amount in a list
        :param distances: list of distances
        :param critical_distance: threshold distance
        :return: index of the first element exceeding the given amount, None if it doesn't happen
        """
        for idx_start, d in enumerate(distances):
            if abs(d) > critical_distance:
                return idx_start
        return None

    def _extract_agent_projected_distances(self, history: SimulationHistory) -> Dict[str, EgoToAgentDistances]:
        """
        Computes the projected distances, for inactive agents only.

        :param history: The history of the scenario
        :return: A dict containing the projected distances to each inactive track in the entire scenario"""
        agents_distances: Dict[str, EgoToAgentDistances] = {}
        inactive_agents_scenario = self._get_inactive_agents_scenario(history)

        for sample in history.data:
            ego_state = sample.ego_state

            assert isinstance(sample.observation, Detections)
            # TODO: refactor the call below to use agent attributes once available
            inactive_agents = self._get_inactive_agents_sample(sample.observation, inactive_agents_scenario)

            lateral_dist = [signed_lateral_distance(ego_state.center, box3d_to_polygon(inactive_agent)) for
                            inactive_agent in inactive_agents.boxes]
            longitudinal_dist = [signed_longitudinal_distance(ego_state.center, box3d_to_polygon(inactive_agent)) for
                                 inactive_agent in inactive_agents.boxes]

            for agent, lon, lat in zip(inactive_agents.boxes, longitudinal_dist, lateral_dist):
                try:
                    agents_distances[agent.token].longitudinal_distances.append(lon)
                    agents_distances[agent.token].lateral_distances.append(lat)
                except KeyError:
                    agents_distances[agent.token] = EgoToAgentDistances(agent.length, [lon], [lat])

        return agents_distances

    def _extract_passing_clearances(self, agents_distances: Dict[str, EgoToAgentDistances]) -> List[float]:
        """
        Extracts the portion of projected distances relative to the passing of every agent and saves them to a list.

        :param agents_distances: The projected distances to each inactive agent
        :return: A list containing the lateral clearance of all inactive agents while ego is passing them.
        """
        clearances_during_overtake = []
        for distances in agents_distances.values():
            max_longitudinal_dist = max(distances.longitudinal_distances)
            idx_max = distances.longitudinal_distances.index(max_longitudinal_dist)
            min_longitudinal_dist = min(distances.longitudinal_distances)
            idx_min = distances.longitudinal_distances.index(min_longitudinal_dist)

            if max_longitudinal_dist > 0 and min_longitudinal_dist < 0 and idx_max < idx_min:
                overtake_idx = np.argmin(np.abs(distances.longitudinal_distances))
                if abs(distances.lateral_distances[overtake_idx]) < self._lateral_distance_threshold:
                    threshold = self._ego_half_length + distances.agent_length / 2.0
                    start_idx = self.get_overtake_start_idx(distances.longitudinal_distances, overtake_idx,
                                                            threshold)
                    end_idx = self.get_overtake_end_idx(distances.longitudinal_distances, overtake_idx,
                                                        threshold)
                    clearances_during_overtake.extend(np.abs(distances.lateral_distances[start_idx:end_idx + 1]))

        return clearances_during_overtake

    @staticmethod
    def _get_inactive_agents_scenario(history: SimulationHistory) -> Set[str]:
        """
        Get a set of agents which are inactive for the full length of the scenario.
        An inactive agents in this context is an agent that for the entire scenario never moves

        :param history: The history from the scenario
        :return: The set of inactive agents"""
        inactive_agents = set()
        active_agents = set()
        for sample in history.data:
            assert isinstance(sample.observation, Detections)
            for box in sample.observation.boxes:
                inactive_agents.add(box.token) if np.isclose(np.linalg.norm(box.velocity), 0.0) else active_agents.add(
                    box.token)

        return inactive_agents.difference(active_agents)

    @staticmethod
    def _get_inactive_agents_sample(agents: Detections, inactive_agents_set: Set[str]) -> Detections:
        """
        Get a set of agents which are inactive for the full length of the scenario.
        An inactive agents in this context is an agent that for the entire scenario never moves

        :param agents: The agents detected in the current sample
        :param inactive_agents_set: The tokens of the agents known to be inactive for the full scenario
        :return: The detections of inactive agents
        """
        return Detections([agent for agent in agents.boxes if agent.token in inactive_agents_set])

max_stars_repo_name: MCZhi/nuplan-devkit

max_issues_repo_name: MCZhi/nuplan-devkit

max_forks_repo_name: MCZhi/nuplan-devkit


text:
import time
import datetime

import matplotlib.pyplot as plt
import pandas as pd
import numpy as np

from dynamic_profile import DynamicProfile

joe = DynamicProfile('https://www.instagram.com/joebiden')
posts = joe.get_all_posts()

#Scrape just the first 500 posts
print(len(posts))
for i, post in enumerate(posts[:500]):
    #sleep between requests so Insta doesn't get mad
    time.sleep(4)
    print(i)
    try:
        post.scrape()
    except Exception as e:
        print(e)

max_stars_repo_path: examples/Joe Biden/profile_scrape.py

max_stars_repo_name: The-Cool-Coders/Instascrape

max_issues_repo_path: examples/Joe Biden/profile_scrape.py

max_issues_repo_name: The-Cool-Coders/Instascrape

max_forks_repo_path: examples/Joe Biden/profile_scrape.py

max_forks_repo_name: The-Cool-Coders/Instascrape


text:
import numpy as np

import arba.space
from .reg import Region

class RegionDiscriminate(Region):
    """ merges regions to minimize population pooled covariance for discrim

    Class Attributes:
        split (dict): keys are groups, values are lists of sbj in that grp

    Instance Attributes:
        t2 (float): hotelling's t^2 (delta.T @ (cov_pool)^-1 @ delta) where
                    delta is difference in grp means and cov_pool is the grp
                    pooled covariance
    """
    split = None

    @classmethod
    def set_split(cls, split):
        cls.split = split

    @staticmethod
    def from_data_img(data_img, ijk=None, pc_ijk=None):
        assert RegionDiscriminate.split is not None, 'call set_split'

        fs_dict = dict()
        for grp, sbj_list in RegionDiscriminate.split.items():
            fs_dict[grp] = data_img.get_fs(ijk=ijk, pc_ijk=pc_ijk,
                                           sbj_list=sbj_list)

        if pc_ijk is None:
            pc_ijk = arba.space.PointCloud([ijk], ref=data_img.ref)
        return RegionDiscriminate(pc_ijk=pc_ijk, fs_dict=fs_dict)

    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)

        self.delta, self.cov_pool, self.t2, self.f = \
            get_t2_stats(fs_dict=self.fs_dict)

        self.cov_pool_det = np.linalg.det(self.cov_pool)

    @staticmethod
    def get_error(reg_1, reg_2, reg_u=None):
        if reg_u is None:
            reg_u = reg_1 + reg_2

        # min mean cov_pool_det of resultant segmentation
        err = reg_u.cov_pool_det * len(reg_u) - \
              reg_1.cov_pool_det * len(reg_1) - \
              reg_2.cov_pool_det * len(reg_2)

        # assert err >= 0, 'err is negative'

        return err

def get_t2_stats(fs_dict, grp_tuple=None):
    if grp_tuple is None:
        grp0, grp1 = tuple(RegionDiscriminate.split.keys())
    else:
        grp0, grp1 = grp_tuple

    delta = fs_dict[grp1].mu - fs_dict[grp0].mu

    # get pool cov
    n0, n1 = fs_dict[grp0].n, fs_dict[grp1].n
    lam0, lam1 = n0 / (n0 + n1), n1 / (n0 + n1)
    cov_pool = fs_dict[grp0].cov * lam0 + \
               fs_dict[grp1].cov * lam1

    # compute t2
    t2 = delta @ np.linalg.inv(cov_pool) @ delta * (n0 * n1) / (n0 + n1)

    # compute f
    p = fs_dict[grp0].d
    f = t2 * (n0 + n1 - p - 1) / ((n0 + n1 - 2) * p)

    return delta, cov_pool, t2, f

max_stars_repo_name: matthigger/arba

max_issues_repo_name: matthigger/arba

max_forks_repo_name: matthigger/arba


text:
import warnings
import os
from PIL import Image

from torch.utils.data import Dataset

import torch

import numpy as np

from skimage import io
from skimage.color import rgb2lab, gray2rgb, rgba2rgb

try:
    import pyift.pyift as ift
except:
    warnings.warn("PyIFT is not installed.", ImportWarning)

__all__ = ["LIDSDataset", "ToTensor", "ToLAB"]

class LIDSDataset(Dataset):
    def __init__(self, root_dir, split_dir=None, transform=None, return_name=False):
        self.root_dir = root_dir
        self.split_dir = split_dir
        self.transform = transform

        self.return_name = return_name
        self.images_names = None
        self.opf_data = None
        self.opf_labels = None

        if self.root_dir.endswith('.zip'):
            self.opf_data, self.opf_labels = self._get_data_from_opfdataset()
        else:
            self.images_names = self._list_images_files()

    def __len__(self):
        if self.opf_data is not None:
            return self.opf_data.shape[0]
        return len(self.images_names)

    def __getitem__(self, idx):
        if self.opf_data is None:
            image_path = os.path.join(self.root_dir, self.images_names[idx])

            if image_path.endswith('mimg'):
                image = ift.ReadMImage(image_path).AsNumPy().squeeze()

            else:
                image = io.imread(image_path)

            if image.ndim == 2:
                image = gray2rgb(image)

            elif image.shape[2] == 4:
                image = rgba2rgb(image)

            image = rgb2lab(image)
            #image = image/(np.array([[116], [500], [200]])).reshape(1, 1, 3)

            image = image.astype(np.float32)

            label = self._label_of_image(self.images_names[idx])

        else:

            image = self.opf_data[idx]
            label = self.opf_labels[idx]

        if(self.transform):
            image = self.transform(image)

        if self.return_name:
            sample = (image, label, self.images_names[idx])
        else:
            sample = (image, label)

        return sample

    def _label_of_image(self, image_name):
        if not isinstance(image_name, str):
            raise TypeError("Parameter image_name must be a string.")
        i = image_name.index("_")
        label = int(image_name[0:i]) - 1

        return label

    def _list_images_files(self):
        if self.split_dir is not None:
            with open(self.split_dir, 'r') as f:
                _filenames = f.read()
                filenames = [filename for filename in _filenames.split('\n') if len(filename) > 0]
        else:

            filenames = os.listdir(self.root_dir)

        return filenames

    def _get_data_from_opfdataset(self):
        opfdataset = ift.ReadDataSet(self.root_dir)

        data = opfdataset.GetData()
        labels = opfdataset.GetTrueLabels().astype(np.int64) - 1

        return data, labels

    def weights_for_balance(self, nclasses):
        weights_dir = os.path.join(self.root_dir, '.weights-for-balance-{}.npy'.format(self.image_list))

        if(os.path.exists(weights_dir)):
            weight = np.load(weights_dir)
            return weight

        count = [0] * nclasses
        for image_name in self.images_names:
            label = self._label_of_image(image_name)
            count[label] += 1
        weight_per_class = [0.] * nclasses
        N = float(sum(count))
        for i in range(nclasses):
            weight_per_class[i] = N/float(count[i])
        weight = [0] * self.__len__()
        for idx, image_name in enumerate(self.images_names):
            label = self._label_of_image(image_name)
            weight[idx] = weight_per_class[label]
        weight = np.array(weight, dtype=np.float)

        np.save(weights_dir, weight)

        return weight

class ToTensor(object):
    def __call__(self, sample):
        image = np.array(sample)

        min = image.min()
        max = image.max()
        image = (image)/(max)

        if image.ndim > 2:
            image = image.transpose((2, 0, 1))

        return torch.from_numpy(image.copy())

class ToLAB(object):
    def __call__(self, sample):
        image = np.array(sample)

        image = rgb2lab(image)
        return image

max_stars_repo_name: LIDS-UNICAMP/FLIM

max_issues_repo_name: LIDS-UNICAMP/FLIM

max_forks_repo_name: LIDS-UNICAMP/FLIM


text:
import codecs
import datetime
#from itertools import chain
#import string
import math
import numpy
import io
import scipy.sparse as sp
from numpy.linalg import inv

class Person:

    userID = ""
    userLat = ""
    userLong = ""
    #Feature_Freq = {}
    total_words = 0
    Feature_Prob = {}
    outside_word_prob = 0
    fileFrom = ""

    def __init__(self, ID, latit, longit, F_Freq, file_from, F_All_Len):
        #self.userID = ID[ID.rfind("_")+1:]
        self.userID = ID
        self.userLat = latit
        self.userLong = longit
        #self.Feature_Freq = F_Freq
        self.fileFrom = file_from
        tw = 0
        for f in F_Freq:
            tw += int(F_Freq[f])
        self.total_words = tw
        self.CalcProb(F_Freq, 'smoothed', 30, False, F_All_Len)

    def CalcProb(self, F_Freq, method="likelihood", res_prob_space = 0, logNormal = False, F_All_Len=1):
        F_Prob = {}
        if method == "smoothed" and res_prob_space > 0:
            perct1 = float(self.total_words) / float(100.0 - (res_prob_space))
            other_words_space = perct1 * res_prob_space
            for word in F_Freq:
                F_Prob[word] = (float(F_Freq[word]) / (float(self.total_words) + float(other_words_space)))
                #print "####Actual Word#####"
                #print F_Prob[word]
            self.outside_word_prob = (float(other_words_space)/(F_All_Len - len(F_Freq))) / (float(self.total_words) + float(other_words_space))
            #print "#####Other Word#####"
            #print F_Prob[word]
        self.Feature_Prob = F_Prob

def BuildXMatrix(people, User_Lookup, std_dev_ref, noise_value=.000001):
    #print "Generating Noise Matrix with values from 0  to ", noise_value
    X_start_time = datetime.datetime.now()
    print "######Building X Matrix ##########"
    XMatrix = numpy.empty( (len(User_Lookup), len(std_dev_ref)) )
    g = 0
    for w in std_dev_ref:
        #noise = rand.random() * noise_value
        for p in people:
            if w in p.Feature_Prob:
                XMatrix[User_Lookup[p.userID]][std_dev_ref[w][0]] = ((p.Feature_Prob[w] - std_dev_ref[w][1]) / std_dev_ref[w][2])
            elif w in std_dev_ref:
                XMatrix[User_Lookup[p.userID]][std_dev_ref[w][0]] = ((p.outside_word_prob - std_dev_ref[w][1]) / std_dev_ref[w][2])
    X_end_time = datetime.datetime.now()
    print "######Finished Building X Matrix ##########"
    print X_end_time - X_start_time
    return XMatrix

def LoadY(wydirect, u):
    yf = open(wydirect + u + '_YVECTOR.bin', 'rb')
    yv = numpy.load(yf)
    yf.close()
    return yv

def LoadWu(wydirect, u, kerntype, numobs):

    wf = open(wydirect + u + '_' + kerntype + "_WUVECTOR.bin", 'rb')
    wv = sp.spdiags(numpy.load(wf).transpose(), [0], numobs, numobs).tocsr()
    wf.close()
    return wv

def StoreScore(u, B, scoredirect, kerntype):
    bf = open(scoredirect + 'Bu_' + u + '_' + kerntype + '.bin', 'wb')
    numpy.save(bf, B)
    bf.close()

def UpToInv(X, W, lam, wordsize):
    Ident = numpy.identity(wordsize)
    print X.shape
    print W.shape
    print Ident.shape
    firstpart = numpy.dot((X.transpose() * W), X)
    print "First part done"
    secondpart = lam * Ident
    print "Second part done"
    i = inv(numpy.add(secondpart, firstpart))
    return i

def train(f, rf_obs_in, rf_std_in, wu_y_direct, ulist, kerntype, lam, scoredirect):
    print "Import Train Worked"

    filename = f[f.rfind('/')+1:]
    print filename
    begin_time = datetime.datetime.now()
    personList = []
    x = 0
    y = 0

    #The length given here was observed from a smaller corpus of about 6000 people
    F_All_Len = 114616

    read_time_begin = datetime.datetime.now()
    z = 0

    with io.open(f, 'r', encoding='utf-8') as f:
        for person in f:
            x += 1
            #print "####NEW Person####"
            #print userID, latit, longit
            try:
                row = person.split('\t')
                #print row[0]
                userID = row[0]
                latit = row[1].split(',')[0]
                longit = row[1].split(',')[1]
                F_Freq = dict(f.split(':') for f in row[2].split(" "))
                #F_All = dict(chain(F_Freq.iteritems(), F_All.iteritems()))
                #print len(F_Freq)
                #F_Freq = {}
                newPerson = Person(userID, latit, longit, F_Freq, filename, F_All_Len)
                #newPerson.CalcProb('smoothed', 30, False, F_All_Len)
                personList.append(newPerson)
            except:
                print "@@@@@error reading user@@@@@@"
                print row
                print z
                break
            z += 1
            if z % 5000 == 0:
                print z
                print datetime.datetime.now()

    URef = {}
    with codecs.open(rf_obs_in, 'r', encoding='utf-8') as w:
        for line in w:
            row = line.strip().split('\t')
            URef[row[0]] = int(row[1])

    WordRef = {}
    i = 0
    with io.open(rf_std_in, 'r', encoding='utf-8') as w:
        for line in w:
            row = line.strip().split('\t')
            if len(row) == 3:
                WordRef[row[0]] = (i, float(row[1]), float(row[2]))
                i += 1
            else: print "Problem Reading row from std ref:", row

    print "Num Obs: ", len(URef)
    print "Num Vars: ", len(WordRef)
    X_M = BuildXMatrix(personList, URef, WordRef)
    print "######X Matrix Build Complete#######"

    for u in ulist:
        print "----------",u,"-------------"
        ustart = datetime.datetime.now()
        Y_V = LoadY(wu_y_direct, u)
        Wu = LoadWu(wu_y_direct, u, kerntype, len(URef))
        Inver = UpToInv(X_M, Wu, lam, len(WordRef))
        print "Invert Finished"
        B = numpy.dot((numpy.dot(Inver, X_M.transpose()) * Wu), Y_V)
        StoreScore(u, B, scoredirect, kerntype)
        print "Total time for u"
        print datetime.datetime.now() - ustart

    print "Total Train Time: "
    print datetime.datetime.now() - read_time_begin

max_stars_repo_name: grantdelozier/TextGWR

max_issues_repo_name: grantdelozier/TextGWR

max_forks_repo_name: grantdelozier/TextGWR

reason: import numpy,from numpy,import scipy