text:
from ._stopping_criterion import StoppingCriterion
from ..accumulate_data import MLMCData
from ..discrete_distribution import IIDStdUniform
from ..true_measure import Gaussian
from ..integrand import MLCallOptions
from ..util import MaxSamplesWarning, ParameterError, MaxLevelsWarning, ParameterWarning
from numpy import *
from scipy.stats import norm
from time import time
import warnings

class CubMCML(StoppingCriterion):
    """
    Stopping criterion based on multi-level monte carlo.

    >>> mlco = MLCallOptions(IIDStdUniform(seed=7))
    >>> sc = CubMCML(mlco,abs_tol=.05)
    >>> solution,data = sc.integrate()
    >>> solution
    10.440...
    >>> data
    Solution: 10.4406
    MLCallOptions (Integrand Object)
        option          european
        sigma           0.200
        k               100
        r               0.050
        t               1
        b               85
    IIDStdUniform (DiscreteDistribution Object)
        d               2^(6)
        seed            7
        mimics          StdUniform
    Gaussian (TrueMeasure Object)
        mean            0
        covariance      1
        decomp_type     pca
    CubMCML (StoppingCriterion Object)
        rmse_tol        0.019
        n_init          2^(8)
        levels_min      2^(1)
        levels_max      10
        theta           2^(-2)
    MLMCData (AccumulateData Object)
        levels          7
        dimensions      [ 1.  2.  4.  8. 16. 32. 64.]
        n_level         [7.795e+05 1.496e+04 5.916e+03 2.244e+03 7.560e+02 2.770e+02 1.060e+02]
        mean_level      [1.005e+01 1.777e-01 1.071e-01 5.340e-02 2.990e-02 1.167e-02 7.812e-03]
        var_level       [1.959e+02 1.441e-01 4.433e-02 1.093e-02 2.940e-03 7.304e-04 2.261e-04]
        cost_per_sample [ 1.  2.  4.  8. 16. 32. 64.]
        n_total         804118
        alpha           0.942
        beta            1.893
        gamma           1.000
        time_integrate  ...

    Original Implementation:

        http://people.maths.ox.ac.uk/~gilesm/mlmc/#MATLAB

    References:

        [1] M.B. Giles. 'Multi-level Monte Carlo path simulation'.
        Operations Research, 56(3):607-617, 2008.
        http://people.maths.ox.ac.uk/~gilesm/files/OPRE_2008.pdf.
    """

    parameters = ['rmse_tol','n_init','levels_min','levels_max','theta']

    def __init__(self, integrand, abs_tol=.05, alpha=.01, rmse_tol=None, n_init=256., n_max=1e10,
        levels_min=2., levels_max=10., alpha0=-1., beta0=-1., gamma0=-1.):
        """
        Args:
            integrand (Integrand): integrand with multi-level g method
            abs_tol (float): absolute tolerance. Reset if supplied, ignored if not.
            alpha (float): uncertaintly level.
                If rmse_tol not supplied, then rmse_tol = abs_tol/norm.ppf(1-alpha/2)
            rel_tol (float): relative tolerance. Reset if supplied, ignored if not.
                Takes priority over aboluste tolerance and alpha if supplied.
            n_init (int): initial number of samples
            n_max (int): maximum number of samples
            levels_min (int): minimum level of refinement >= 2
            levels_max (int): maximum level of refinement >= Lmin
            alpha0 (float): weak error is O(2^{-alpha0*level})
            beta0 (float): variance is O(2^{-bet0a*level})
            gamma0 (float): sample cost is O(2^{gamma0*level})

        Note:
            if alpha, beta, gamma are not positive, then they will be estimated
        """
        if levels_min < 2:
            raise ParameterError('needs levels_min >= 2')
        if levels_max < levels_min:
            raise ParameterError('needs levels_max >= levels_min')
        if n_init <= 0:
            raise ParameterError('needs n_init > 0')
        # initialization
        if rmse_tol:
            self.rmse_tol = float(rmse_tol)
        else: # use absolute tolerance
            self.rmse_tol =  float(abs_tol) / norm.ppf(1-alpha/2)
        self.n_init = float(n_init)
        self.n_max = float(n_max)
        self.levels_min = float(levels_min)
        self.levels_max = float(levels_max)
        self.theta = 0.25
        self.alpha0 = alpha0
        self.beta0 = beta0
        self.gamma0 = gamma0
        # QMCPy Objs
        self.integrand = integrand
        self.true_measure = self.integrand.true_measure
        self.discrete_distrib = self.integrand.discrete_distrib
        # Verify Compliant Construction
        allowed_levels = ['adaptive-multi']
        allowed_distribs = ["IIDStdUniform"]
        super(CubMCML,self).__init__(allowed_levels, allowed_distribs)

    def integrate(self):
        """ See abstract method. """
         # Construct AccumulateData Object to House Integration Data
        self.data = MLMCData(self, self.integrand, self.true_measure, self.discrete_distrib,
            self.levels_min, self.n_init, self.alpha0, self.beta0, self.gamma0)
        t_start = time()
        while self.data.diff_n_level.sum() > 0:
            self.data.update_data()
            self.data.n_total += self.data.diff_n_level.sum()
            # set optimal number of additional samples
            n_samples = self._get_next_samples()
            self.data.diff_n_level = maximum(0, n_samples-self.data.n_level)
            # check if over sample budget
            if (self.data.n_total + self.data.diff_n_level.sum()) > self.n_max:
                warning_s = """
                Alread generated %d samples.
                Trying to generate %d new samples, which would exceed n_max = %d.
                Stopping integration process.
                Note that error tolerances may no longer be satisfied""" \
                % (int(self.data.n_total), int(self.data.diff_n_level.sum()), int(self.n_max))
                warnings.warn(warning_s, MaxSamplesWarning)
                break
            # if (almost) converged, estimate remaining error and decide
            # whether a new level is required
            if (self.data.diff_n_level > 0.01*self.data.n_level).sum() == 0:
                range_ = arange(minimum(2.,self.data.levels-1)+1)
                idx = (self.data.levels-range_).astype(int)
                rem = (self.data.mean_level[idx] /
                        2.**(range_*self.data.alpha)).max() / (2.**self.data.alpha - 1)
                # rem = ml(l+1) / (2^alpha - 1)
                if rem > sqrt(self.theta)*self.rmse_tol:
                    if self.data.levels == self.levels_max:
                        warnings.warn(
                            'Failed to achieve weak convergence. levels == levels_max.',
                            MaxLevelsWarning)
                    else:
                        self.data.levels += 1
                        self.data.dimensions = hstack((self.data.dimensions,0))
                        self.data.var_level = hstack((self.data.var_level,
                            self.data.var_level[-1] / 2**self.data.beta))
                        self.data.cost_per_sample = hstack((self.data.cost_per_sample,
                            self.data.cost_per_sample[-1] * 2**self.data.gamma))
                        self.data.n_level = hstack((self.data.n_level, 0.))
                        self.data.sum_level = hstack((self.data.sum_level,
                            zeros((2,1))))
                        self.data.cost_level = hstack((self.data.cost_level, 0.))
                        n_samples = self._get_next_samples()
                        self.data.diff_n_level = maximum(0., n_samples-self.data.n_level)
        # finally, evaluate multilevel estimator
        self.data.solution = (self.data.sum_level[0,:]/self.data.n_level).sum()
        self.data.levels += 1
        self.data.time_integrate = time() - t_start
        return self.data.solution,self.data

    def _get_next_samples(self):
        ns = ceil( sqrt(self.data.var_level/self.data.cost_per_sample) *
                sqrt(self.data.var_level*self.data.cost_per_sample).sum() /
                ((1-self.theta)*self.rmse_tol**2) )
        return ns

    def set_tolerance(self, abs_tol=None, alpha=.01, rmse_tol=None):
        """
        See abstract method.

        Args:
            abs_tol (float): absolute tolerance. Reset if supplied, ignored if not.
            alpha (float): uncertaintly level.
                If rmse_tol not supplied, then rmse_tol = abs_tol/norm.ppf(1-alpha/2)
            rel_tol (float): relative tolerance. Reset if supplied, ignored if not.
                Takes priority over aboluste tolerance and alpha if supplied.
        """
        if rmse_tol != None:
            self.rmse_tol = float(rmse_tol)
        elif abs_tol != None:
            self.rmse_tol = (float(abs_tol) / norm.ppf(1-alpha/2.))

ext: py

max_stars_repo_name: jungtaekkim/QMCSoftware

max_issues_repo_name: jungtaekkim/QMCSoftware

max_forks_repo_name: jungtaekkim/QMCSoftware

reason: from numpy,from scipy


text:
from abc import ABCMeta, abstractmethod
import numpy as np
import time
import logging

from src.parameters import Parameters
from src.solvers import Records

""" The class Solver() can be declined in many subclasses, each defining a specific solver
    Each solver can be different from the others, but here they share:

    - Storing some essential information such as the Problem(), the Parameters()
    - Running one epoch via .run_epoch (it is used later in run and run_repetition)
    - Recording some value of iterest to be plotted later
        Those values are stored in Solver.records which is a dict of Records
        See the class Records() for more details
"""

class Solver(metaclass=ABCMeta):

    def __init__(self, problem):
        if not hasattr(self, 'name'): # sometimes the name is already defined
            self.name = ""
        self.problem = problem
        self.loss = problem.loss
        self.regularizer = problem.regularizer
        self.param = Parameters(problem.data) # can create problem with *reload*
        self.x = None
        self.epoch_running_time = 0.0
        self.total_running_time = 0.0
        self.records = {} # { record_name : Record() }
        if self.param.measure_time:
            self.append_records("time_epoch")
    """
    def append_records(self, *tuple_of_class):
        # tag the Records with the name of the Solver they'll belong to
        for cls in tuple_of_class:
            cls.solver_name = self.name
            cls.data_name = self.problem.data.name
        update = { cls.name : cls for cls in tuple_of_class }
        self.records = { **self.records, **update }
    """
    def append_records(self, *tuple_of_class_name):
        # all the existing records:
        list_records = [ record for record in Records.__subclasses__() if record.__name__[0] != '_' ]
        dict_of_records = {}
        # tag the Records with the name of the Solver they'll belong to
        for cls in list_records:
            if cls.name in tuple_of_class_name: # we want this Record for this Solver
                record_instance = cls(self)
                record_instance.solver_name = self.name
                record_instance.data_name = self.problem.data.name
                dict_of_records = { **dict_of_records, cls.name : record_instance }
        self.records = { **self.records, **dict_of_records }

    # Runs the Solver. It calls a handful of subroutines, defined after
    def run(self):
        # startup
        self.initialization_variables() # put all states at "zero"
        self.store_records() # updates self.records[all].value at iteration 0
        # main loop
        for cnt in range(1, self.param.nb_epoch + 1):
            # one epoch
            if self.param.measure_time:
                start_time = time.time()
                self.run_epoch()
                self.epoch_running_time = time.time() - start_time
                self.total_running_time += self.epoch_running_time
            else: # we just run the epoch
                self.run_epoch()
            # storing, verbosing, and stop-or-keep-looping
            self.store_records()
            self.log_information(cnt)
            if self.we_have_to_break():
                break # out of the for loop
        return

    # Some of the blocks needed to .run() the Solver
    def initialization_variables(self):
        # standard, can be overwritten if needed in child Solver
        # the initial iterate
        init = self.param.initialization
        if isinstance(init, str):
            if init == "zeros":
                self.x = np.zeros(self.problem.dim)
            # we could do "random" as well but need to be careful with seed
        else: # we assume it is a np.array
            self.x = init
        # the measure of time
        self.epoch_running_time = 0.0
        self.total_running_time = 0.0

    def store_records(self):
        # All the work is done in Records().store()
        for record in self.records.values():
            record.store(self)

    @abstractmethod
    def run_epoch(self):
        pass # This is to be defined for each child Solver

    def log_information(self, cnt):
        text = f"| Solver : {self.name} "
        text += f"| End of epoch #{cnt:d} "
        if self.param.measure_time:
            text += f"| ran in {self.epoch_running_time:f}s "
        for record in self.records.values():
            if record.name != "time_epoch": # to avoid doublons
                text += f"| {record.name} {record.value[-1]:f} "
        text += "|"
        logging.info(text)

    def we_have_to_break(self):
        return self.param.tol is not None and self.records["gradient_norm"].value[-1] <= self.param.tol

    # An upper function running a solver multiple times
    def run_repetition(self):
        np.random.seed(0)  # random seed to reproduce the experiments
        logging.info(f"------START {self.name}------")
        for i in range(self.param.n_repetition):
            logging.info(f"{i+1}-th repetition:")
            self.run() # one run of the solver, see above
            self.archive_records_repetition() # archive records and prepare for an another run
        logging.info(f"------END {self.name}------")
        return

    def archive_records_repetition(self):
        for record in self.records.values(): # record is a Records() class
            record.value_repetition.append(record.value) # carved in stone
            record.value = [] # run is over so we can wipe it
        return

    # Other needed methods
    def save_records(self):
        # save the list of lists contained in each Records.value_repetition
        # or maybe self.records itself?? TBC
        for record in self.records.values():
            record.save()

reason: import numpy


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Credits: Benjamin Dartigues, Emmanuel Bouilhol, Hayssam Soueidan, Macha Nikolski

import pathlib
import time

import numpy as np
from loguru import logger

import constants
import helpers
import plot
from image_set import ImageSet
# this should be called as soon as possible
from path import global_root_dir
from plot import plot_heatmap

def compute_relative_densities(analysis_repo, molecule_type, quadrants_num=4, peripheral_flag=False):
    densities = {}
    stripes = constants.analysis_config['STRIPE_NUM']
    if peripheral_flag:
        stripes_flag = False
        stripes = 1
    else:
        stripes_flag = True

    if molecule_type == 'mrna':
        timepoints = constants.dataset_config['TIMEPOINTS_MRNA']
    else:
        timepoints = constants.dataset_config['TIMEPOINTS_PROTEIN']

    for gene in constants.analysis_config['PROTEINS']:
        median_densities, gene_clusters = [], []
        for timepoint in timepoints:
            image_set = ImageSet(analysis_repo, [molecule_type + "/{0}/{1}/".format(gene, timepoint)])
            arr = image_set.compute_normalised_quadrant_densities(quadrants_num=quadrants_num,
                                                                  peripheral_flag=peripheral_flag,
                                                                  stripes=stripes, stripes_flag=stripes_flag)
            num_images = arr.shape[0] // (quadrants_num * stripes)
            aligned_densities = arr[:, 0].reshape(num_images, quadrants_num * stripes)
            median_densities_per_slice = np.nanmedian(aligned_densities, axis=0)
            median_densities.append(median_densities_per_slice)
        densities[gene] = median_densities

    return densities

# configurations contain the order in which the degree of clustering is plotted
configurations = [
    ["src/analysis/colocalization/config_original.json"],
    ["src/analysis/colocalization/config_original_periph.json"],
    ["src/analysis/colocalization/config_nocodazole_arhgdia.json"],
    ["src/analysis/colocalization/config_nocodazole_pard3.json"]
]
# Figure 6C top Peripheral Colocalization Score (CS) for original data
# Figure 6C top Cytoplasmic Colocalization Score (CS) for original data
# Figure 6D Colocalization Score (CS) for original data (4 figures)
# Figure 7F top left Colocalization Score (CS) for arhgdia data
# Figure 7F top right Colocalization Score (CS) for nocodazole arhgdia data
# Figure 7F bottom left Colocalization Score (CS) for pard3 data
# Figure 7F bottom right Colocalization Score (CS) for nocodazole pard3 data

if __name__ == '__main__':
    np.random.seed(int(round(time.time())))
    for conf in configurations:
        logger.info("Colocalization Score")
        conf_full_path = pathlib.Path(global_root_dir, conf[0])
        constants.init_config(analysis_config_js_path=conf_full_path)
        repo = helpers.open_repo()
        peripheral_flag = "periph" in conf[0]
        stripes = constants.analysis_config["STRIPE_NUM"]
        mrna_densities = compute_relative_densities(repo, 'mrna', quadrants_num=8, peripheral_flag=peripheral_flag)
        prot_densities = compute_relative_densities(repo, 'protein', quadrants_num=8, peripheral_flag=peripheral_flag)

        css, results = [], {}
        for gene in constants.analysis_config['PROTEINS']:
            cs, p, ranking = helpers.calculate_colocalization_score(mrna_densities[gene], prot_densities[gene],
                                                                    constants.dataset_config['TIMEPOINTS_NUM_MRNA'],
                                                                    constants.dataset_config['TIMEPOINTS_NUM_PROTEIN'],
                                                                    peripheral_flag, stripes=stripes, quadrants=8)
            css.append(cs)
            results[gene] = [cs, p]
            tgt_image_name = constants.analysis_config['FIGURE_NAME_FORMAT_CS'].format(gene=gene)
            tgt_fp = pathlib.Path(constants.analysis_config['FIGURE_OUTPUT_PATH'].format(root_dir=global_root_dir),
                                  tgt_image_name)
            if "original" in conf[0]:
                plot_heatmap(ranking, gene, tgt_fp)
                logger.info("Generated image at {}", str(tgt_fp).split("analysis/")[1])
            else:
                plot_heatmap(ranking, gene, tgt_fp, size=2, xtickslabel=['3h', '5h'], ytickslabel=['3h', '5h'])
                logger.info("Generated image at {}", str(tgt_fp).split("analysis/")[1])

        tgt_image_name = constants.analysis_config['FIGURE_NAME_FORMAT_CS_HISTOGRAM']
        tgt_fp = pathlib.Path(constants.analysis_config['FIGURE_OUTPUT_PATH'].format(root_dir=global_root_dir), tgt_image_name)
        plot.bar_profile(css, tgt_fp, constants.analysis_config['PLOT_COLORS'])
        logger.info("Generated image at {}", str(tgt_fp).split("analysis/")[1])
        logger.info("Colocalization scores and p-values are: {}", results)

max_stars_repo_name: ebouilhol/dypfish

max_issues_repo_name: ebouilhol/dypfish

max_forks_repo_name: ebouilhol/dypfish


text:
import torch
import torch.utils.data as data
from PIL import Image
from PIL import ImageChops
import os
import math
import functools
import json
import copy

from utils import load_value_file
import random
import cv2
import numpy as np
import pdb

def pil_loader(path):
    # open path as file to avoid ResourceWarning (https://github.com/python-pillow/Pillow/issues/835)
    with open(path, 'rb') as f:
        with Image.open(f) as img:
            return img.convert('RGB')

def accimage_loader(path):
    try:
        import accimage
        return accimage.Image(path)
    except IOError:
        # Potentially a decoding problem, fall back to PIL.Image
        return pil_loader(path)

def get_default_image_loader():
    from torchvision import get_image_backend
    if get_image_backend() == 'accimage':
        return accimage_loader
    else:
        return pil_loader

def video_loader(video_dir_path, frame_indices, image_loader):
    video = []
    for i in frame_indices:
        image_path = os.path.join(video_dir_path, 'image_{:05d}.jpg'.format(i))
        image_path2 = os.path.join(video_dir_path, 'image_{:05d}.png'.format(i))
        if os.path.exists(image_path):
            video.append(image_loader(image_path))
        elif os.path.exists(image_path2):
            video.append(image_loader(image_path2))
        else:
            return video

    return video

def get_default_video_loader():
    image_loader = get_default_image_loader()
    return functools.partial(video_loader, image_loader=image_loader)

def make_dataset(root_path, subset, n_samples_for_each_video, sample_duration):

    video_input_path = os.path.join(root_path, 'X')
    if subset == 'training':
        video_target_path = os.path.join(root_path, 'Y')
    dataset = []
    for video_name in os.listdir(video_input_path):
        input_path = os.path.join(video_input_path, video_name)
        if subset == 'training':
            target_name = 'Y' + video_name[1:]
            target_path = os.path.join(video_target_path, target_name)
            if not (os.path.exists(input_path) and os.path.exists(target_path)):
                continue
        else:
            target_path = None
        n_frames = 125
        begin_t = 1
        end_t = n_frames
        sample = {
            'video':input_path,
            'segment': [begin_t, end_t],
            'n_frames': 125,
            'video_id':video_name[1:],
            'target_video' :target_path
        }
        if n_samples_for_each_video == 1:
            sample['frame_indices'] = list(range(1, n_frames + 1))
            dataset.append(sample)
        else: # TEST
            if n_samples_for_each_video > 1:
                step = max(1, math.ceil((n_frames - 1 - sample_duration) /
                                     (n_samples_for_each_video - 1)))
            elif n_samples_for_each_video ==0:
                step = sample_duration

            if n_samples_for_each_video ==-1: # SHRINK
                for j in range(1, n_frames+1): # 1 ~ 125
                    sample_j = copy.deepcopy(sample)
                    sample_j['frame_indices'] = [j] # j
                    dataset.append(sample_j)
            else:
                for j in range(1, n_frames, step):
                    sample_j = copy.deepcopy(sample)
                    sample_j['frame_indices'] = list( range(j, min(n_frames + 1, j + sample_duration)))
                    dataset.append(sample_j)

    return dataset

class VideoDecaptionData(data.Dataset):
    """
    Args:
        root (string): Root directory path.
        spatial_transform (callable, optional): A function/transform that  takes in an PIL image
            and returns a transformed version. E.g, ``transforms.RandomCrop``
        temporal_transform (callable, optional): A function/transform that  takes in a list of frame indices
            and returns a transformed version
        target_transform (callable, optional): A function/transform that takes in the
            target and transforms it.
        loader (callable, optional): A function to load an video given its path and frame indices.
     Attributes:
        classes (list): List of the class names.
        class_to_idx (dict): Dict with items (class_name, class_index).
        imgs (list): List of (image path, class_index) tuples
    """

    def __init__(self,
                 root_path,
                 subset,
                 n_samples_for_each_video=1,
                 spatial_transform=None,
                 temporal_transform=None,
                 sample_duration=16,
                 get_loader=get_default_video_loader,
                 opt=None):
        self.subset = subset
        if subset == 'testing':
            n_samples_for_each_video = -1
        self.data = make_dataset(root_path, subset, n_samples_for_each_video, sample_duration)

        self.spatial_transform = spatial_transform
        self.target_transform = spatial_transform
        self.temporal_transform = temporal_transform

        self.loader = get_loader()

        self.opt = opt
        self.lr_flip = opt.lr_flip
        self.tb_flip = opt.tb_flip
        self.t_stride = opt.t_stride
        self.t_shrink = opt.t_shrink
    def __getitem__(self, index):
        """
        Args:
            index (int): Index
        Returns:
            tuple: (image, target) where target is class_index of the target class.
        """
        if self.subset == 'training':
            path = self.data[index]['video']
            target_path = self.data[index]['target_video']
            frame_indices = self.data[index]['frame_indices']

            if self.temporal_transform is not None:
                frame_indices = self.temporal_transform(frame_indices)
                if self.t_stride > 1:
                    frame_indices = frame_indices[0::self.t_stride]

            clip = self.loader(path, frame_indices)
            if self.t_shrink:
                mid_idx = len(frame_indices)//2
                target_indices = [frame_indices[mid_idx]]
                target_clip = self.loader(target_path, target_indices)
            else:
                target_clip = self.loader(target_path, frame_indices)

            if self.spatial_transform is not None:
                self.spatial_transform.randomize_parameters()
                clip = [self.spatial_transform(img) for img in clip]
                target_clip = [self.spatial_transform(img) for img in target_clip]

            clip = torch.stack(clip, 0).permute(1, 0, 2, 3)
            target_clip = torch.stack(target_clip, 0).permute(1, 0, 2, 3)

            """ Random Horizontal & Vertical Flip """
            if self.lr_flip and random.random() < 0.5:
                idx = [i for i in range(clip.size(3)-1, -1, -1)]
                idx = torch.LongTensor(idx)
                clip = clip.index_select(3, idx)
                target_clip = target_clip.index_select(3, idx)
            if self.tb_flip and random.random() < 0.5:
                idx = [i for i in range(clip.size(2)-1, -1, -1)]
                idx = torch.LongTensor(idx)
                clip = clip.index_select(2, idx)
                target_clip = target_clip.index_select(2, idx)

            return clip, target_clip

        else:
            path = self.data[index]['video']
            frame_indices = self.data[index]['frame_indices']
            if self.temporal_transform is not None:
                if self.t_shrink:
                    frame_indices_ = []
                    for i in range(self.opt.sample_duration):
                        frame_indices_.append(frame_indices[0]+self.t_stride*i)
                    offset = (len(frame_indices_)//2)*self.t_stride
                    frame_indices = [x-offset for x in frame_indices_]
                    frame_indices = [-x+2 if x<=0 else x for x in frame_indices]
                    frame_indices = [125*2-x if x>125 else x for x in frame_indices]
                else:
                    frame_indices = self.temporal_transform(frame_indices)
                    if self.t_stride > 1:
                        frame_indices = frame_indices[0::self.t_stride]

            clip = self.loader(path, frame_indices)

            if self.spatial_transform is not None:
                self.spatial_transform.randomize_parameters()
                clip = [self.spatial_transform(img) for img in clip]

            clip = torch.stack(clip, 0).permute(1, 0, 2, 3)

            return clip, path

    def __len__(self):
        return len(self.data)


text:
import numpy as np
import cv2

image = cv2.imread('scan.png')

img = cv2.resize(image, (960, 540))
imCanny = cv2.Canny(img, 150, 200)

contours, hierarchy = cv2.findContours(imCanny, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)
imgsh = img.copy()

final_con = []
for cnt in contours:
    area = cv2.contourArea(cnt)
    if area < 7900 and area > 1300:
        final_con.append(cnt)
        peri = cv2.arcLength(cnt, True)
        approx = cv2.approxPolyDP(cnt, 0.02 * peri, True)
        cv2.drawContours(imgsh, [approx], -1, (0, 255, 0), 2)

idx = 0
for c in final_con:
    x,y,w,h = cv2.boundingRect(c)
    if w<100:
        idx+=1
        new_img=img[y:y+h,x:x+w]
        cv2.imwrite(str(idx) + '.png', new_img)

cv2.imshow('finalImg2', imCanny)
cv2.imshow('finalImg', imgsh)
cv2.imshow("im",img)
cv2.waitKey(0)

max_stars_repo_name: FatinShadab/opencv-projects-python

max_issues_repo_name: FatinShadab/opencv-projects-python

max_forks_repo_name: FatinShadab/opencv-projects-python


text:
# --------------------------------------------------------
# Tensorflow Faster R-CNN
# Licensed under The MIT License [see LICENSE for details]
# Written by Zheqi He, Xinlei Chen, based on code from Ross Girshick
# --------------------------------------------------------
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import _init_paths
from model.train_val import get_training_roidb, train_net
from model.config import cfg, cfg_from_file, cfg_from_list, get_output_dir, get_output_tb_dir
from datasets.factory import get_imdb, get_vg_imdb
import datasets.imdb
import argparse
import pprint
import numpy as np
import sys

import tensorflow as tf
from nets.vgg16 import vgg16
from nets.resnet_v1 import resnetv1
from nets.mobilenet_v1 import mobilenetv1

def parse_args():
    """
    Parse input arguments
    """
    parser = argparse.ArgumentParser(description='Train a Fast R-CNN network')
    parser.add_argument('--cfg', dest='cfg_file',
                        help='optional config file',
                        default='/home/wangwenbin/Method/scene_graph/tf-faster-rcnn/experiments/cfgs/vgg16.yml',
                        type=str)
    parser.add_argument('--weight', dest='weight',
                        help='initialize with pretrained model weights',
                        default='/home/wangwenbin/Method/scene_graph/tf-faster-rcnn-relPN/data/imagenet_weights/vgg16.ckpt',
                        type=str)
    parser.add_argument('--imdb', dest='imdb_name',
                        help='dataset to train on',
                        default='visual_genome_train_det', type=str)
    parser.add_argument('--imdbval', dest='imdbval_name',
                        help='dataset to validate on',
                        default='visual_genome_val_det', type=str)
    parser.add_argument('--iters', dest='max_iters',
                        help='number of iterations to train',
                        default=70000, type=int)
    parser.add_argument('--tag', dest='tag',
                        help='tag of the model',
                        default=None, type=str)
    parser.add_argument('--net', dest='net',
                        help='vgg16, res50, res101, res152, mobile',
                        default='vgg16', type=str)
    parser.add_argument('--set', dest='set_cfgs',
                        help='set config keys', default=None,
                        nargs=argparse.REMAINDER)
    parser.add_argument('--train_size', dest='train_size',
                        default=-1, type=int)
    parser.add_argument('--val_size', dest='val_size',
                        default=1, type=int)

    args = parser.parse_args()
    return args

def combined_roidb(split, task, num_im, num_val_im):
    """
    Combine multiple roidbs
    """
    imdb = get_vg_imdb(split, task, num_im, num_val_im)
    print('Loaded dataset `{:s}` for training'.format(imdb.name))
    roidb = get_training_roidb(imdb)
    return imdb, roidb

if __name__ == '__main__':
    args = parse_args()

    print('Called with args:')
    print(args)

    if args.cfg_file is not None:
        cfg_from_file(args.cfg_file)
    if args.set_cfgs is not None:
        cfg_from_list(args.set_cfgs)

    print('Using config:')
    pprint.pprint(cfg)

    np.random.seed(cfg.RNG_SEED)

    # train set
    split, task = args.imdb_name.split('_')[2:]
    imdb, roidb = combined_roidb(split, task, num_im=args.train_size, num_val_im=args.val_size)
    print('Loaded imdb `{:s}` for training'.format(args.imdb_name))
    print('{:d} roidb entries'.format(len(roidb)))

    # also add the validation set, but with no flipping images
    orgflip = cfg.TRAIN.USE_FLIPPED
    cfg.TRAIN.USE_FLIPPED = False
    split, task = args.imdbval_name.split('_')[2:]
    valimdb, valroidb = combined_roidb(split, task, num_im=args.val_size, num_val_im=-1)
    cfg.TRAIN.USE_FLIPPED = orgflip
    print('Loaded imdb `{:s}` for validating'.format(args.imdbval_name))
    print('{:d} validation roidb entries'.format(len(valroidb)))

    # output directory where the models are saved
    output_dir = get_output_dir(imdb, args.tag)
    print('Output will be saved to `{:s}`'.format(output_dir))

    # tensorboard directory where the summaries are saved during training
    tb_dir = get_output_tb_dir(imdb, args.tag)
    print('TensorFlow summaries will be saved to `{:s}`'.format(tb_dir))

    # load network
    if args.net == 'vgg16':
        net = vgg16()
    elif args.net == 'res50':
        net = resnetv1(num_layers=50)
    elif args.net == 'res101':
        net = resnetv1(num_layers=101)
    elif args.net == 'res152':
        net = resnetv1(num_layers=152)
    elif args.net == 'mobile':
        net = mobilenetv1()
    else:
        raise NotImplementedError

    train_net(net, imdb, roidb, valimdb, valroidb, output_dir, tb_dir,
              pretrained_model=args.weight,
              max_iters=args.max_iters)


text:
import keras
from keras.models import Sequential
from keras.layers import Dense, Activation, Dropout
from keras.layers import LSTM
from keras.optimizers import RMSprop
import numpy as np
import os

import datetime
import cPickle

import argparse
import sys

TRAIN_FILENAME= 'data.pkl'
RNN_SIZE = 256
N_EPOCHS = 20
N_LAYERS = 1
LEARNING_RATE= 0.01
DROPOUT = 0.0
RECURRENT_DROPOUT = 0.0

TEMPERATURE = 1.0
PRIMETEXT = 'Diagnosis'
LENGTH = 1000

def get_parser():
    parser = argparse.ArgumentParser(prog='p3b1_baseline',
                                     formatter_class=argparse.ArgumentDefaultsHelpFormatter)

    parser.add_argument("-v", "--verbose", action="store_true",
                        default= True,
                        help="increase output verbosity")

    parser.add_argument("-e", "--n_epochs", action="store",
                        default=N_EPOCHS, type=int,
                        help="number of training epochs")
    parser.add_argument('-l', '--learning_rate', action='store',
                        default= LEARNING_RATE, type=float,
                        help='learning rate')
    parser.add_argument("--dropout", action="store",
                        default=DROPOUT, type=float,
                        help="ratio of dropout")
    parser.add_argument("--recurrent_dropout", action="store",
                        default=RECURRENT_DROPOUT, type=float,
                        help="ratio of recurrent dropout")

    parser.add_argument("--train_data", action="store",
                        default= TRAIN_FILENAME,
                        help='training data filename')

    parser.add_argument("--output_dir", action="store",
                        default='out',
                        help="output directory")

    parser.add_argument("--rnn_size", action="store",
                        default= RNN_SIZE,
                        help='size of LSTM internal state')
    parser.add_argument("--n_layers", action="store",
                        default= N_LAYERS,
                        help='number of layers in the LSTM')

    parser.add_argument("--do_sample", action="store_true",
                        default= True,
                        help="generate synthesized text")
    parser.add_argument("--temperature", action="store",
                        default= TEMPERATURE, type= float,
                        help="variability of text synthesis")
    parser.add_argument("--primetext", action="store",
                        default= PRIMETEXT,
                        help= 'seed string of text synthesis' )
    parser.add_argument("--length", action="store",
                        default= LENGTH, type= int,
                        help= 'length of synthesized text')

    return parser

class LossHistory( keras.callbacks.Callback ):
    def on_train_begin( self, logs= {} ):
        self.losses = []

    def on_batch_end( self, batch, logs= {} ):
        self.losses.append( logs.get( 'loss' ) )

def sample( preds, temperature= 1.0 ):
    # helper function to sample an index from a probability array
    preds = np.asarray( preds ).astype( 'float64' )
    preds = np.log( preds ) / temperature
    exp_preds = np.exp( preds )
    preds = exp_preds / np.sum( exp_preds )
    probas = np.random.multinomial( 1, preds, 1 )
    return np.argmax( probas )

def char_rnn(
             rnn_size= 256,
             n_layers= 1,
             learning_rate= 0.01,
             dropout= 0.0,
             recurrent_dropout= 0.0,
             n_epochs= 100,
             data_train= 'data.pkl',
             verbose= 1,
             savedir= '',
             do_sample= True,
             temperature= 1.0,
             primetext= 'Diagnosis',
             length= 1000
             ):

        # load data from pickle
        f = open( data_train, 'r' )

        classes = cPickle.load( f )
        chars = cPickle.load( f )
        char_indices = cPickle.load( f )
        indices_char = cPickle.load( f )

        maxlen = cPickle.load( f )
        step = cPickle.load( f )

        X_ind = cPickle.load( f )
        y_ind = cPickle.load( f )

        f.close()

        [ s1, s2 ] = X_ind.shape
        print( X_ind.shape )
        print( y_ind.shape )
        print( maxlen )
        print( len( chars ) )

        X = np.zeros( ( s1, s2, len( chars ) ), dtype=np.bool )
        y = np.zeros( ( s1, len( chars ) ), dtype=np.bool )

        for i in range( s1 ):
            for t in range( s2 ):
                X[ i, t, X_ind[ i, t ] ] = 1
            y[ i, y_ind[ i ] ] = 1

        # build the model: a single LSTM
        if verbose:
            print( 'Build model...' )

        model = Sequential()

        # for rnn_size in rnn_sizes:
        for k in range( n_layers ):
            if k < n_layers - 1:
                ret_seq = True
            else:
                ret_seq = False

            if k == 0:
                model.add( LSTM( rnn_size, input_shape= ( maxlen, len( chars ) ), return_sequences= ret_seq,
                                 dropout= dropout, recurrent_dropout= recurrent_dropout ) )
            else:
                model.add( LSTM( rnn_size, dropout= dropout, recurrent_dropout= recurrent_dropout, return_sequences= ret_seq ) )

        model.add( Dense( len( chars ) ) )
        model.add( Activation( 'softmax' ) )

        optimizer = RMSprop( lr= learning_rate )
        model.compile( loss= 'categorical_crossentropy', optimizer= optimizer )

        if verbose:
            model.summary()

        for iteration in range( 1, n_epochs ):
            if verbose:
                print()
                print('-' * 50)
                print('Iteration', iteration)

            history = LossHistory()
            model.fit( X, y, batch_size= 100, epochs= 1, callbacks= [ history ] )

            loss = history.losses[ -1 ]
            if verbose:
                print( loss )

            dirname = savedir
            if len( dirname ) > 0 and not dirname.endswith( '/' ):
                dirname = dirname + '/'

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

            # serialize model to JSON
            model_json = model.to_json()
            with open( dirname + "/model_" + str( iteration ) + "_" + str( round( loss, 6 ) ) + ".json", "w" ) as json_file:
                json_file.write( model_json )

            # serialize weights to HDF5
            model.save_weights( dirname + "/model_" + str( iteration ) + "_" + str( round( loss, 6 ) ) + ".h5" )

            if verbose:
                print( "Checkpoint saved." )

            if do_sample:
                outtext = open( dirname + "/example_" + str( iteration ) + "_" + str( round( loss, 6 ) ) + ".txt", "w" )

                diversity = temperature

                outtext.write('----- diversity:' + str( diversity ) + "\n" )

                generated = ''
                seedstr = primetext

                outtext.write('----- Generating with seed: "' + seedstr + '"' + "\n" )

                sentence = " " * maxlen

                # class_index = 0
                generated += sentence
                outtext.write( generated )

                for c in seedstr:
                    sentence = sentence[1:] + c
                    x = np.zeros( ( 1, maxlen, len( chars ) ) )
                    for t, char in enumerate(sentence):
                        x[ 0, t, char_indices[ char ] ] = 1.

                    preds = model.predict( x, verbose= verbose )[ 0 ]
                    next_index = sample( preds, diversity )
                    next_char = indices_char[ next_index ]

                    generated += c

                    outtext.write( c )

                for i in range( length ):
                    x = np.zeros( ( 1, maxlen, len( chars ) ) )
                    for t, char in enumerate( sentence ):
                        x[ 0, t, char_indices[ char ] ] = 1.

                    preds = model.predict( x, verbose= verbose )[ 0 ]
                    next_index = sample( preds, diversity )
                    next_char = indices_char[ next_index ]

                    generated += next_char
                    sentence = sentence[ 1 : ] + next_char

                    outtext.write( next_char )

                outtext.write( "\n" )

            outtext.close()

if __name__  == "__main__":
    parser = get_parser()
    args = parser.parse_args()
    print( 'Args:', args )

    if args.verbose:
        verbose = 1
    else:
        verbose = 0

    ## Read files
    file_path = os.path.dirname(os.path.realpath(__file__))
    print file_path
    lib_path = os.path.abspath(os.path.join(file_path, '..', '..', 'common'))
    sys.path.append(lib_path)

    from data_utils import get_file
    origin = 'http://ftp.mcs.anl.gov/pub/candle/public/benchmarks/P3B2/P3B2_data.tgz'
    data_loc = get_file('P3B2_data.tgz', origin, untar=True, md5_hash=None, cache_subdir='P3B2')

    print 'Data downloaded and stored at: ' + data_loc
    data_path = os.path.dirname(data_loc)
    print data_path

    char_rnn(
        rnn_size= args.rnn_size,
        n_layers= args.n_layers,
        learning_rate= args.learning_rate,
        dropout= args.dropout,
        recurrent_dropout= args.recurrent_dropout,
        n_epochs= args.n_epochs,
        data_train= data_path+'/data.pkl',
        verbose= args.verbose,
        savedir= args.output_dir,
        do_sample= args.do_sample,
        temperature= args.temperature,
        primetext= args.primetext,
        length= args.length
    )


text:
#!/usr/bin/env python2
# -*- coding: utf-8 -*-

import csv
import cv2
import numpy as np
import os
import random
from keras.models import Sequential
from keras.layers import Flatten, Dense, Lambda, Conv2D, AveragePooling2D, Cropping2D
from sklearn.model_selection import train_test_split

os.environ['KMP_DUPLICATE_LIB_OK']='True'

samples = []
images = []
measurements = []

def LoadMeasurements(dirs):
    '''
        Load the measurement from all listed directories
    '''
    for dir in dirs:
        LoadMeasurement(dir)

def LoadMeasurement(dir):
    '''
        Load measurement from one directory
    '''
    with open(dir + 'driving_log.csv') as csvfile:
        reader = csv.reader(csvfile)
        for line in reader:
            samples.append((line, dir))

LoadMeasurements(['data/track_1_data/', 'data/recoveries2/'])

train_samples, validation_samples = train_test_split(samples, test_size=0.2)

def generator(samples, batch_size=32):
    '''
        Setup generator to generation input data
    '''
    num_samples = len(samples)
    while 1:
        random.shuffle(samples)
        for offset in range(0, num_samples, batch_size):
            batch_samples = samples[offset:offset+batch_size]

            images = []
            measurements = []
            for batch_sample in batch_samples:
                line, dir = batch_sample
                image_dir = dir + 'IMG/'

                filename = line[0].split('/')[-1]
                image = cv2.imread(image_dir + filename)
                image = cv2.cvtColor(image, cv2.COLOR_BGR2RGB)

                # Use left and write pictures with a correction factor
                filename_left = line[1].split('/')[-1]
                image_left = cv2.imread(image_dir + filename_left)
                #Convert from BGR to RGB
                image_left = cv2.cvtColor(image_left, cv2.COLOR_BGR2RGB)

                filename_right = line[2].split('/')[-1]
                image_right = cv2.imread(image_dir + filename_right)
                #Convert from BGR to RGB
                image_right = cv2.cvtColor(image_right, cv2.COLOR_BGR2RGB)

                images.append(image)
                images.append(image_left)
                images.append(image_right)

                correction = 0.2

                measurement = float(line[3])
                measurement_left = measurement + correction
                measurement_right = measurement - correction

                measurements.append(measurement)
                measurements.append(measurement_left)
                measurements.append(measurement_right)

            X_train = np.array(images)
            y_train = np.array(measurements)
            yield sklearn.utils.shuffle(X_train, y_train)

# Set our batch size
batch_size=32

# compile and train the model using the generator function
train_generator = generator(train_samples, batch_size=batch_size)
validation_generator = generator(validation_samples, batch_size=batch_size)

X_train = np.array(images)
y_train = np.array(measurements)

def MorePowerfulNet():
    '''
        Setup NVidia conv net architecture
    '''
    model.add(Cropping2D(cropping=((50,20), (0,0)), input_shape=(160,320,3)))
    model.add(Lambda(lambda x: (x/255.0) - 0.5))
    model.add(Conv2D(filters=24, strides=(2,2), kernel_size=(5, 5), activation='relu'))
    model.add(Conv2D(filters=36,  strides=(2,2), kernel_size=(5, 5), activation='relu'))
    model.add(Conv2D(filters=48, kernel_size=(5, 5), activation='relu'))
    model.add(Conv2D(filters=64, kernel_size=(3, 3), activation='relu'))
    model.add(Conv2D(filters=64, kernel_size=(3, 3), activation='relu'))
    model.add(Flatten())
    model.add(Dense(100))
    model.add(Dense(50))
    model.add(Dense(10))
    model.add(Dense(1))

model = Sequential()
MorePowerfulNet()
model.compile(loss='mse', optimizer='adam', metrics=['accuracy'])

model.fit_generator(train_generator, steps_per_epoch=math.ceil(len(train_samples)/batch_size),
                    validation_data=validation_generator,
                    validation_steps=math.ceil(len(validation_samples)/batch_size), epochs=3, verbose=1)
model.save('model.h5')


text:
"""
Expose each GPU device directly
"""
from __future__ import print_function, absolute_import, division
import functools
from numba import servicelib
from .driver import hsa as driver, Context as _Context

class _culist(object):
    """A thread local list of GPU instances
    """

    def __init__(self):
        self._lst = None

    @property
    def _gpus(self):
        if not self._lst:
            self._lst = self._init_gpus()
        return self._lst

    def _init_gpus(self):
        gpus = []
        for com in driver.components:
            gpus.append(CU(com))
        return gpus

    def __getitem__(self, item):
        return self._gpus[item]

    def append(self, item):
        return self._gpus.append(item)

    def __len__(self):
        return len(self._gpus)

    def __nonzero__(self):
        return bool(self._gpus)

    def __iter__(self):
        return iter(self._gpus)

    __bool__ = __nonzero__

    def reset(self):
        for gpu in self:
            gpu.reset()

    @property
    def current(self):
        """Get the current GPU object associated with the thread
        """
        return _custack.top

cus = _culist()
del _culist

class CU(object):
    def __init__(self, cu):
        self._cu = cu
        self._context = None

    def __getattr__(self, key):
        """Redirect to self._gpu
        """
        if key.startswith('_'):
            raise AttributeError(key)
        return getattr(self._cu, key)

    def __repr__(self):
        return repr(self._cu)

    def associate_context(self):
        """Associate the context of this GPU to the running thread
        """
        # No context was created for this GPU
        if self._context is None:
            self._context = self._cu.create_context()

        return self._context

    def __enter__(self):
        self.associate_context()
        _custack.push(self)

    def __exit__(self, exc_type, exc_val, exc_tb):
        assert _get_device() is self
        self._context.pop()
        _custack.pop()

    def reset(self):
        if self._context:
            self._context.reset()
            self._context = None

_cpu_context = None

def get_cpu_context():
    global _cpu_context
    if _cpu_context is None:
        cpu_agent = [a for a in driver.agents if not a.is_component][0]
        _cpu_context = _Context(cpu_agent)
    return _cpu_context

def get_gpu(i):
    return cus[i]

def get_num_gpus():
    return len(cus)

_custack = servicelib.TLStack()

def _get_device(devnum=0):
    """Get the current device or use a device by device number.
    """
    if not _custack:
        _custack.push(get_gpu(devnum))
    return _custack.top

def get_context(devnum=0):
    """Get the current device or use a device by device number, and
    return the HSA context.
    """
    return _get_device(devnum=devnum).associate_context()

def get_all_contexts():
    return [get_context(i) for i in range(get_num_gpus())]

def require_context(fn):
    """
    A decorator to ensure a context for the HSA subsystem
    """

    @functools.wraps(fn)
    def _require_cu_context(*args, **kws):
        get_context()
        return fn(*args, **kws)

    return _require_cu_context

def reset():
    cus.reset()
    _custack.clear()

max_stars_repo_name: mawanda-jun/numba

max_issues_repo_name: olivier-be/lumberyard

max_forks_repo_name: olivier-be/lumberyard

reason: from numba


text:
#!/usr/bin/env python

import numpy as np
import pandas as pd
import slab2functions as s2f
import psutil

def loop1(lons, lats, testarea, slab, depgrid, strgrid, dipgrid,
        slab1query, eventlist, seismo_thick,
        alen, blen, mdist, sdr, ddr, mindip, maxID2, AA_data,
        TR_data, maxdist, maxthickness, minstk,
        tomo_sets, meanBA, slab1guide, spacing, slab1data, dipthresh, datainfo, nodeinfo, j):

    ''' Arguments:  lons - list of longitudes (list of floats)[deg]
                    lats - list of latitudes (list of floats)[deg]
                    nodeIDs - list of nodeIDs (list of ints)
                    testarea - list defined as [lonmin,lonmax,latmin,latmax]
                            nodes within this box will print troubleshooting
                            statements and PDFs and data used to create them
                            are saved (list of floats)[deg]
                    slab - three letter code to indicate slab of interest (str)
                    depgrid - clipped depth GMT grid file from Slab 1.0 (grd)
                    strgrid - grid of strike values calculated from
                            Slab 1.0 (grd)
                    dipgrid - grid of dip values calculated from Slab 1.0 (grd)
                    slab1data - dataframe with columns [lon,lat,depth,strike,
                            dip] representing Slab1.0 (dataframe)
                    eventlist - dataframe with columns [lat,lon,depth,unc,etype,
                            ID,mag,S1,D1,R1,S2,D2,R2,src] representing input
                            dataset (dataframe)
                    seismo_thick - depth defining where earthquakes at and
                            below are intra-slab as opposed to on the
                            slab surface. (float)[km]
                    alen - long radius of ellipsoid (float)[km]
                    blen - short radius of search ellipsoid (float)[km]
                    mdist - alen*ecentricity of ellipse (float)[km]
                    sdr - shallow search range of ellipsoid (float)[km]
                    ddr - deep search range of ellipsoid (float)[km]
                    mindip - upper limit of dip to search vertically. Nodes with
                            dips greater than mindip have a perpendicular search
                            ellipsoid (float)[deg]
                    maxID2 - maxID of original input dataset (int)
                    AA_data - dataframe with columns [dist,depth].
                            if slab == 'alu' or slab == 'him', columns =
                            [dist,depth,avlat,avlon] (dataframe)
                    TR_data - dataframe with columns [lon,lat,az,
                                bound,slab] (dataframe)
                    maxdist - useless variable, need to remove
                    maxthickness - maximum thickness of this slab as calculated
                            by a function relating thickness to age (float)[km]
                    minstk - maximum strike difference between node and farthest
                            data point constraining that node (float)[deg]
                    tomo_sets - list of tomography datasets for this slab.
                            (list of strings)
                    OGmulti - dataframe with columns [lon,lat,depth,nID] used to
                            record which nodes have a multi-modal distribution
                            and the depths to each peak (dataframe)
                    i - index representing this node. (int)

        Returns:    lon - lons[i], search longitude of this node (float)[deg]
                    lat - lats[i], search latitude of this node (float)[deg]
                    locdep - search depth of this node (float)[km]
                    locstr - search strike of this node (float)[deg]
                    locdip - search dip of this node (float)[deg]
                    nID - integer value representing this node (int)
                    farpoint - length of long axis of ellipse after
                            refiltering based on strike difference from
                            node center (float)[km]
                    used_TO - list of tomography data points used to constrain
                            nodes on this slab. 2D numpy array with columns
                            [mean(|tomo_depths-EQdepths|),index representing
                            TO source] (array of floats)
                    used_tmp - list of data points used to constrain the depth
                            at this node. 2D numpy array with columns
                            [nodeID,dataID] (array of ints)
                    trimmedAA - dataframe with columns [lat,lon,depth,unc,etype,
                            ID,mag,S1,D1,R1,S2,D2,R2,src] representing average
                            active source and receiver functions to
                            add to input dataset (dataframe) '''

    # Define lon,lat coordinate, integer ID of this node, and ID for added data
    lat = lats[j]
    lon = lons[j]
    if lat>0:
        nID = int('%i%i'%(lon*10,lat*10))
    else:
        nID = int('%i0%i'%(lon*10,lat*-10))
    maxID = maxID2+nID-1
    OGmulti = []

    if slab == 'sol' and lon > 149 and lon < 155:
        alen = 20.0
    if slab == 'puy' and lat < -48:
        alen = 100

    if slab == 'ryu' and lon > 138:
        alen = 20.0

    if slab == 'alu' and lat > 63:
        alen = 50.0

    if slab == 'man' and lat > 17.0:
        alen = 50.0

    if slab == 'izu' and lat < 27 and lat > 23:
        eventlist = eventlist[(eventlist.lon<lon+4.5)&(eventlist.lat<lat+4.5)&(eventlist.lon>lon-4.5)&(eventlist.lat>lat-4.5)]
    else:
        eventlist = eventlist[(eventlist.lon<lon+2.5)&(eventlist.lat<lat+2.5)&(eventlist.lon>lon-2.5)&(eventlist.lat>lat-2.5)]
    if len(eventlist)<2:
        #print ('no data here',lon,lat,nID)
        return lon, lat, np.nan, np.nan, np.nan, nID, np.nan, np.nan, np.nan, [], [], [], [], False, 0, 0

    # If in troubleshooting range, print ID and set testprint = True
    if lon>testarea[0] and lon<testarea[1]:
        if lat>testarea[2] and lat<testarea[3]:
            testprint = True
            print('___________________ %i ______________________'% nID)
        else:
            testprint = False
    else:
        testprint = False

    if testprint:
        f = open(nodeinfo, 'a')
        f.write('-%i- pre-anyfilters, len(eventlist) %i\n'%(nID, len(eventlist)))
        f.close()
    # Get depth, strike, and dip of Slab1.0 or calculate from trend at edge
    slab1, strtmp, diptmp, inside, extended, out, extlon, extlat = s2f.getslab12(slab1guide,slab1query,lon,lat,spacing,depgrid,strgrid,dipgrid,testprint,TR_data,meanBA,slab)

    if slab1 > 650 and slab == 'sam' and lat<-30:
        sdr = 60
    if testprint:
        f = open(nodeinfo, 'a')
        f.write('-%i- slab1, strtmp, diptmp, inside, extended, out %0.2f,%0.2f,%0.2f,%s,%s,%s \n'%(nID, slab1, strtmp, diptmp, inside, extended, out))
        f.close()

    if slab == 'izu' and lat < 27 and lat > 23:
        dipthresh = 45

    if diptmp > dipthresh:
        newlats = [lat,lat,lat+spacing/2.0,lat-spacing/2.0,lat+spacing/2.0,lat-spacing/2.0,lat+spacing/2.0,lat-spacing/2.0]
        newlons = [lon-spacing/2.0,lon+spacing/2.0,lon,lon,lon+spacing/2.0,lon-spacing/2.0,lon-spacing/2.0,lon+spacing/2.0]

        if slab == 'ker':
            newlats1 = [lat,lat,lat+spacing/4.0,lat-spacing/4.0,lat+spacing/4.0,lat-spacing/4.0,lat+spacing/4.0,lat-spacing/4.0]
            newlons1 = [lon-spacing/4.0,lon+spacing/4.0,lon,lon,lon+spacing/4.0,lon-spacing/4.0,lon-spacing/4.0,lon+spacing/4.0]

            newlats2 = [lat,lat,lat+3*spacing/4.0,lat-3*spacing/4.0,lat+3*spacing/4.0,lat-3*spacing/4.0,lat+3*spacing/4.0,lat-3*spacing/4.0]
            newlons2 = [lon-3*spacing/4.0,lon+3*spacing/4.0,lon,lon,lon+3*spacing/4.0,lon-3*spacing/4.0,lon-3*spacing/4.0,lon+3*spacing/4.0]

            newlats.extend(newlats1)
            newlats.extend(newlats2)
            newlons.extend(newlons1)
            newlons.extend(newlons2)

        newnodes = np.zeros((len(newlons),2))
        newnodes[:,0] = np.round(newlons,1)
        newnodes[:,1] = np.round(newlats,1)
        #diptmp = 90
        if slab == 'izu' and lat>23 and lat<27 and slab1 > 400:
            clen = 400
        elif slab != 'sol' and slab != 'hal' and slab != 'phi':
            clen = 200
        elif slab == 'sol' and lon>148:
            clen = 300
        else:
            clen = alen
    else:
        eventlist = eventlist[(eventlist.lon<lon+1.5)&(eventlist.lat<lat+1.5)&(eventlist.lon>lon-1.5)&(eventlist.lat>lat-1.5)]
        if testprint:
            f = open(nodeinfo, 'a')
            f.write('-%i- after reducing OGsearch range (dip<dipthresh), len(eventlist) %i\n'%(nID, len(eventlist)))
            f.close()
        clen = alen
        if len(eventlist)<2 and len(eventlist[(eventlist.etype == 'AS')|(eventlist.etype == 'RF')|(eventlist.etype == 'CP')])<1 and slab != 'ryu':
            #print ('no data here',lon,lat,nID)
            return lon, lat, slab1, strtmp, diptmp, nID, alen, blen, clen, [], [], [], [], False, 0, 0
        elif slab == 'ryu' and len(eventlist)<2 and len(eventlist[(eventlist.etype == 'AS')|(eventlist.etype == 'RF')|(eventlist.etype == 'TO')|(eventlist.etype == 'CP')])<1:
            return lon, lat, slab1, strtmp, diptmp, nID, alen, blen, clen, [], [], [], [], False, 0, 0
        else:
            newnodes = []

    if slab == 'phi' and lat > 11:
        minstk /= 2

    if extended:
        alen = s2f.refilter4(strtmp,minstk,extlon,extlat,eventlist,alen,blen,slab1query,slab,testprint)
    else:
        alen = s2f.refilter4(strtmp,minstk,lon,lat,eventlist,alen,blen,slab1query,slab,testprint)

    if lon > 170.5 and lon < 172 and lat > -22 and lat < -21:
        alen = 20
    if testprint:
        f = open(nodeinfo, 'a')
        f.write('-%i- alen, blen, clen, minstk, seismo_thick %0.2f,%0.2f,%0.2f,%0.2f,%0.2f \n'%(nID, alen, blen, clen, minstk, seismo_thick))
        f.close()

    if alen == 1:
    	return lon, lat,slab1, strtmp, diptmp, nID, alen, blen, clen, [], [], [], [], False, 0, 0
    else:
    	farpoint = alen

    if slab == 'cam' and lon < 260 and lat > 20.5 and lon > 257.5:
        alen = 20
        clen = 20

    if slab == 'man' and lat < 14:
        alen = 20
        if slab1 < 200:
            clen = 20

    if slab == 'mue':
        alen = 150
        clen = 150
        strtmp = 270

    # Filter data for present node
    trimmed, test, sdepth, ddepth, locstr, locdip, maxID, locdep = s2f.allFilters(
                                eventlist, lat, lon, inside, slab1, strtmp, diptmp, seismo_thick, alen, blen,
                                clen, mdist, sdr, ddr, mindip, maxID, out,
                                AA_data, TR_data, slab, maxdist, testprint, extended, datainfo, nodeinfo,nID)

    trimmed, test = s2f.doublecheckEREQ(trimmed, lon, lat)

    if not test:
        if testprint:
            print ('not test',lon,lat,nID)
        return lon, lat, locdep, locstr, locdip, nID, alen, blen, clen, [], [], [], newnodes, False, sdepth, ddepth

    trimmedAA = trimmed[trimmed.etype == 'AA']
    trimmedBA = trimmed[trimmed.etype == 'BA']
    trimmedAS = trimmed[trimmed.etype == 'AS']
    trimmedRF = trimmed[trimmed.etype == 'RF']
    trimmedTO = trimmed[trimmed.etype == 'TO']
    trimmedCP = trimmed[trimmed.etype == 'CP']
    # Save data used
    # Add this trimmed list of events to the list of used input data
    if len(trimmed) > 1 or (len(trimmed) < 2 and (len(trimmedAA)>0 or len(trimmedBA)>0 or len(trimmedAS)>0 or len(trimmedRF)>0 or len(trimmedCP)>0)):
        trimmed_index = np.array(trimmed.ID)
        node_number = np.ones(len(trimmed_index)) * nID
        used_tmp = np.array(list(zip(node_number, trimmed_index)))
    elif len(trimmed) > 1 or (len(trimmed) < 2 and len(trimmedTO)>0 and slab == 'ryu'):
        trimmed_index = np.array(trimmed.ID)
        node_number = np.ones(len(trimmed_index)) * nID
        used_tmp = np.array(list(zip(node_number, trimmed_index)))
    else:
        if testprint:
            print ('trimmed too small',lon,lat,nID)
        return lon, lat, locdep, locstr, locdip, nID, alen, blen, clen, [], [], [], newnodes, False, sdepth, ddepth

    # We only need to do the rest of this loop if there is tomography data
    if tomo_sets == 0 or slab != 'samzzz' or slab =='ryu':
        return lon, lat, locdep, locstr, locdip, nID, alen, blen, clen, [], used_tmp, trimmedAA, newnodes, True, sdepth, ddepth

    '''nothing beyond this point in this loop gets used '''

    # Trim dataset for PDF calculation and separate tomography data
    # We want to run PDF on non-TO, non-RF data
    TOdata = trimmed[trimmed.etype == 'TO']
    trimmed = trimmed[trimmed.etype != 'RF']
    trimmed = trimmed[trimmed.etype != 'TO']
    trimmed = trimmed[trimmed.etype != 'CP']

    # Calculate PDF
    if locdip <= mindip: #or len(trimmedAS) > 0 or len(trimmedAA) >0 or len(trimmedRF) >0:
        #print 'full',lon,lat
        peak_depth, stdv, test2, npeaks, OGmulti, centsurf = s2f.fullPDFcalc(trimmed, sdepth, ddepth, testprint, nID, lat, lon, locdep, 'pre', slab, locstr, locdip)
        peak_lon = lon
        peak_lat = lat
    else:
        peak_lon, peak_lat, peak_depth, stdv, test2, npeaks, OGmulti, centsurf, rlist = s2f.perpPDFcalc(trimmed, sdepth, ddepth, testprint, nID, lat, lon, locdep, 'pre', slab, locstr, locdip, maxthickness)
        if len(rlist) > 0:
            removeIDs = np.array(rlist.ID)
            used_tmp = used_tmp[~np.in1d(used_tmp[:, 1], removeIDs)]

    # For each node with a peak depth and tomography data, calculate differences
    # Do separately for each tomography dataset
    if test2:
        if len(TOdata) > 0:
            for idx, src in enumerate(tomo_sets):
                data = TOdata[TOdata.src == src]
                if len(data) > 0:
                    TOmean = np.mean(data.depth)
                    diff = peak_depth - TOmean
                    temp = [diff, idx]
                    temp2 = np.array(temp)
                    if used_TO is not None:
                        used_TO = np.vstack((used_TO, temp2))
                    else:
                        used_TO = temp2.copy()
                else:
                    return lon, lat, locdep, locstr, locdip, nID, farpoint, used_TO, used_tmp, trimmedAA, newnodes

    return lon, lat, locdep, locstr, locdip, nID, alen, blen, clen, used_TO, used_tmp, trimmedAA, newnodes

def loop2(testarea, lons, lats, nIDs1, locdep, locstr, locdip, used_all, eventlist, sdr, ddr, seismo_thick, slab, maxthickness, rlist, mindip, aalen, ablen, aclen, j):
    these_parameters = []

    premulti = []
    bilats, bilons, binods, bistds = np.nan, np.nan, np.nan, np.nan
    biindx, bistrs, bidips, bideps = np.nan, np.nan, np.nan, np.nan

    lat = lats[j]
    lon = lons[j]
    nID = nIDs1[j]
    cstr = locstr[j]
    cdip = locdip[j]
    alen = aalen[j]
    blen = ablen[j]
    clen = aclen[j]

    if lon>testarea[0] and lon<testarea[1] and lat>testarea[2] and lat<testarea[3]:
        testprint = True
        print('____________________________________ %i ___________________________________________'% nID)
    else:
        testprint = False

    # Collect used data associated with this node
    loc_used_data = used_all[:, 1][used_all[:, 0] == nID]

    if len(loc_used_data) > 1:
        trimmed = eventlist[eventlist['ID'].isin(loc_used_data)]
    elif len(loc_used_data) > 0:
        trimmed = eventlist[eventlist['ID'].isin(loc_used_data)]
        trimmedAA = trimmed[trimmed.etype == 'AA']
        trimmedBA = trimmed[trimmed.etype == 'BA']
        trimmedAS = trimmed[trimmed.etype == 'AS']
        trimmedRF = trimmed[trimmed.etype == 'RF']
        trimmedCP = trimmed[trimmed.etype == 'CP']
        if len(trimmedAA)<1 and len(trimmedBA)<1 and len(trimmedAS)<1 and len(trimmedRF)<1 and len(trimmedCP)<1:
            return np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, rlist, premulti, alen, blen, clen, 0
    else:
        return np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, rlist, premulti, alen, blen, clen, 0

    # Remove RF data for PDF analysis
    if locdep[j] > seismo_thick:
        rft = trimmed[(trimmed.etype == 'RF') | (trimmed.etype == 'AA') | (trimmed.etype == 'AS') | (trimmed.etype == 'CP')]
        trimmed = trimmed[(trimmed.etype != 'RF') & (trimmed.etype != 'AA') & (trimmed.etype != 'AS') & (trimmed.etype != 'CP')]
    else:
        rft = []

    if len(trimmed)<2 and len(rft)>0:
        trimmed = pd.concat([rft, trimmed],sort=True)
        if abs(rft['depth'].mean()-trimmed['depth'].mean())>maxthickness:
            trimmed = trimmed[(trimmed.etype == 'RF') | (trimmed.etype == 'AA') | (trimmed.etype == 'AS') | (trimmed.etype == 'CP')]
    if len(rft)>0:
        loc_depth = np.mean(trimmed['depth'].values)
        if np.isnan(loc_depth):
            if testprint:
                print ('1 ** nan depth found for this node, returning nan output in this directory',nID)
                trimmed.to_csv('%i_faileddata.csv'%nID,header=True,index=False,na_rep=np.nan)
            return np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, rlist, premulti, alen, blen, clen, 0
    else:
        loc_depth = locdep[j]
        if np.isnan(loc_depth):
            if testprint:
                print ('2 ** nan depth found for this node, returning nan output in this directory',nID)
                trimmed.to_csv('%i_faileddata.csv'%nID,header=True,index=False,na_rep=np.nan)
            return np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, rlist, premulti, alen, blen, clen, 0

    place_holder, sdepth, ddepth, depthwritten = s2f.depthRange(loc_depth, sdr, ddr, seismo_thick, trimmed, slab, these_parameters, True)

    if cdip <= mindip:
        if cdip > mindip:
            cdip = mindip
        #try:
        peak_depth, stdv, test2, npeaks, premulti, centsurf = s2f.fullPDFcalc(trimmed, sdepth, ddepth, testprint, nID, lat, lon, loc_depth, 'pre2', slab, cstr, cdip)
        peak_lon = lon
        peak_lat = lat
    else:
        peak_lon, peak_lat, peak_depth, stdv, test2, npeaks, premulti, centsurf, rlist = s2f.perpPDFcalc(trimmed, sdepth, ddepth, testprint, nID, lat, lon, loc_depth, 'pre2', slab, cstr, cdip, maxthickness)

    if test2:
        #tmp_res[j,0],tmp_res[j,1],tmp_res[j,2],tmp_res[j,3],tmp_res[j,4],tmp_res[j,5],tmp_res[j,6],tmp_res[j,7],tmp_res[j,8],tmp_res[j,9] = peak_lat,peak_lon,peak_depth,std,nID,lat,lon,cstr,cdip,centsurf
        jpeak_lat, jpeak_lon, jpeak_depth, jstd, jnID = peak_lat, peak_lon, peak_depth, stdv, nID
        jlat, jlon, jcstr, jcdip, jcentsurf = lat, lon, cstr, cdip, centsurf
        # If there is more than one peak
        if npeaks > 1:
            bilats, bilons, binods, bistds= lat, lon, nID, stdv
            biindx, bistrs, bidips, bideps = j, cstr, cdip, loc_depth
    else:
        jpeak_lon, jpeak_lat, jpeak_depth, jstd, jnID = np.nan, np.nan, np.nan, np.nan, np.nan
        jlat, jlon, jcstr, jcdip, jcentsurf = np.nan, np.nan, np.nan, np.nan, np.nan

    eventtypes = list(set(trimmed.etype))
    if len(eventtypes) == 1 and eventtypes[0] == 'TO':
        onlyTO = 1
    else:
        onlyTO = 0

    return jpeak_lon, jpeak_lat, jpeak_depth, jstd, jnID, jlat, jlon, jcstr, jcdip, jcentsurf, bilats, bilons, binods, bistds, biindx, bistrs, bidips, bideps, rlist, premulti, alen, blen, clen, onlyTO

def loop3(shift_out, testarea, used_all, eventlist, sdr, ddr, seismo_thick, these_parameters, slab, maxthickness, mindip, taper, nID):

    postmulti = []
    bilats, bilons, binods, bistds = np.nan, np.nan, np.nan, np.nan
    biindx, bistrs, bidips, bideps = np.nan, np.nan, np.nan, np.nan

    if len(shift_out)>0: # used to be if slab == 'sam' bc of variable taper
        taper = 0.0

    loc_shift_out = shift_out[shift_out.nID == nID]
    lat = loc_shift_out['lat'].values[0]
    lon = loc_shift_out['lon'].values[0]
    cstr = loc_shift_out['ogstr'].values[0]
    cdip = loc_shift_out['ogdip'].values[0]
    bflat = loc_shift_out['bzlat'].values[0]
    bflon = loc_shift_out['bzlon'].values[0]
    dep = loc_shift_out['psdepth'].values[0]

    if len(loc_shift_out) == 0:
        #print 'misnamed ID',nID
        return np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, nID, postmulti, np.nan, np.nan
    else:
        testprint = False
    if lon>testarea[0] and lon<testarea[1] and lat>testarea[2] and lat<testarea[3]:
        testprint = True
    else:
        testprint = False
    if testprint:
        print('____________________________________ %i ___________________________________________'% nID)

    # Collect used data associated with this node
    loc_used_data = used_all[:, 1][used_all[:, 0] == nID]
    if len(loc_used_data) > 1:
        trimmed = eventlist[eventlist['ID'].isin(loc_used_data)]
    elif len(loc_used_data) > 0:
        trimmed = eventlist[eventlist['ID'].isin(loc_used_data)]
        trimmedAA = trimmed[trimmed.etype == 'AA']
        trimmedBA = trimmed[trimmed.etype == 'BA']
        trimmedAS = trimmed[trimmed.etype == 'AS']
        trimmedRF = trimmed[trimmed.etype == 'RF']
        trimmedCP = trimmed[trimmed.etype == 'CP']
        if len(trimmedAA)<1 and len(trimmedBA)<1 and len(trimmedAS)<1 and len(trimmedRF)<1 and len(trimmedCP)<1:
            #if testprint:
            print('skipping but returning something wonky', lon, lat, dep)
            return np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, nID, postmulti, np.nan, np.nan
    else:
        return np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, nID, postmulti, np.nan, np.nan
    # Calculate average strike, rake, dip of CMTs for final output
    avg_strike, avg_dip, avg_rake = s2f.avStrDipRak(trimmed)

    if dep < seismo_thick or (len(trimmed[trimmed.etype == 'RF']) < 1 and \
                    len(trimmed[trimmed.etype == 'AA']) < 1 and \
                    len(trimmed[trimmed.etype == 'CP']) < 1 and \
                    len(trimmed[trimmed.etype == 'AS']) < 1):
        #print 'shallow',lon,lat,nID
        rstd = loc_shift_out['stdv'].values[0]
        depth = loc_shift_out['depth'].values[0]
        if dep<seismo_thick:
            return depth, rstd, avg_strike, avg_dip, avg_rake, bilats, bilons, binods, bistds, biindx, bistrs, bidips, bideps, nID, postmulti, bflon, bflat
        else:
            return depth, rstd, avg_strike, avg_dip, avg_rake, bilats, bilons, binods, bistds, biindx, bistrs, bidips, bideps, nID, postmulti, lon, lat

    if len(trimmed[(trimmed.etype != 'RF') & (trimmed.etype != 'CP')])<1 and len(trimmed[(trimmed.etype == 'RF') | (trimmed.etype == 'CP')])>0:
        rstd = loc_shift_out['stdv'].values[0]
        depth = loc_shift_out['psdepth'].values[0]
        return depth, rstd, avg_strike, avg_dip, avg_rake, bilats, bilons, binods, bistds, biindx, bistrs, bidips, bideps, nID, postmulti, bflon, bflat

    if len(trimmed[(trimmed.etype != 'RF') & (trimmed.etype != 'AS') & (trimmed.etype != 'AA') & (trimmed.etype != 'CP')]) > 0:
        # Calculate local depth difference between shifted and non-shifted results
        dep_dif = dep - loc_shift_out['depth'].values[0]
        lon_dif = bflon - lon
        lat_dif = bflat - lat

        # Apply shift to non-RF data
        trimmed['depth'][(trimmed.etype != 'RF') & (trimmed.etype != 'AS') & (trimmed.etype != 'AA') & (trimmed.etype != 'CP')] -= dep_dif
        trimmed['lon'][(trimmed.etype != 'RF') & (trimmed.etype != 'AS') & (trimmed.etype != 'AA') & (trimmed.etype != 'CP')] -= lon_dif
        trimmed['lat'][(trimmed.etype != 'RF') & (trimmed.etype != 'AS') & (trimmed.etype != 'AA') & (trimmed.etype != 'CP')] -= lat_dif

    # Calculate new PDF
    loc_depth = trimmed['depth'].mean()
    place_holder, sdepth, ddepth, depthwritten = s2f.depthRange(loc_depth, sdr, ddr, seismo_thick, trimmed, slab, these_parameters, True)

    if cdip <= mindip:
        if cdip > mindip:
            cdip = mindip
        try:
            new_peak_depth, new_std, test3, npeaks, postmulti, centsurf = s2f.fullPDFcalc(trimmed, sdepth, ddepth, testprint, nID, lat, lon, loc_depth, 'post', slab, cstr, cdip)
            peak_lon = lon
            peak_lat = lat
        except:
            print('post full pdf didnt work!!', lon, lat, nID)
            new_peak_depth, new_std, test3, npeaks, postmulti, centsurf = loc_shift_out['depth'].values[0], loc_shift_out['stdv'].values[0], False, 1, pd.DataFrame(), 2000
            peak_lon = lon
            peak_lat = lat
    else:
        try:
            peak_lon, peak_lat, new_peak_depth, new_std, test3, npeaks, postmulti, centsurf, rlist = s2f.perpPDFcalc(trimmed, sdepth, ddepth, testprint, nID, lat, lon, loc_depth, 'post', slab, cstr, cdip, maxthickness)
        except:
            print('post perp pdf didnt work!!', lon, lat, nID)
            peak_lon, peak_lat, new_peak_depth, new_std, test3, npeaks, postmulti, centsurf, rlist = lon, lat, loc_shift_out['depth'].values[0], loc_shift_out['stdv'].values[0], False, 1, pd.DataFrame(), 2000, pd.DataFrame()
    # Populate results array
    if test3:
        rdepth, rstd, rstrike, rdip, rrake = new_peak_depth, new_std, avg_strike, avg_dip, avg_rake
        if npeaks > 1:
            bilats, bilons, binods, bistds= lat, lon, nID, new_std
            biindx, bistrs, bidips, bideps = 1, cstr, cdip, loc_depth
    else:
        rdepth, rstd, rstrike, rdip, rrake = np.nan, np.nan, np.nan, np.nan, np.nan

    return rdepth, rstd, rstrike, rdip, rrake, bilats, bilons, binods, bistds, biindx, bistrs, bidips, bideps, nID, postmulti, peak_lon, peak_lat


text:
from collections import Iterable, deque
from slm_lab.agent.memory.base import Memory
from slm_lab.lib import logger, util
from slm_lab.lib.decorator import lab_api
import copy
import numpy as np
import pydash as ps

logger = logger.get_logger(__name__)

class OnPolicyReplay(Memory):
    '''
    Stores agent experiences and returns them in a batch for agent training.

    An experience consists of
        - state: representation of a state
        - action: action taken
        - reward: scalar value
        - next state: representation of next state (should be same as state)
        - done: 0 / 1 representing if the current state is the last in an episode
        - priority (optional): scalar value, unnormalized

    The memory does not have a fixed size. Instead the memory stores data from N episodes, where N is determined by the user. After N episodes, all of the examples are returned to the agent to learn from.

    When the examples are returned to the agent, the memory is cleared to prevent the agent from learning from off policy experiences. This memory is intended for on policy algorithms.

    Differences vs. Replay memory:
        - Experiences are nested into episodes. In Replay experiences are flat, and episode is not tracked
        - The entire memory constitues a batch. In Replay batches are sampled from memory.
        - The memory is cleared automatically when a batch is given to the agent.

    e.g. memory_spec
    "memory": {
        "name": "OnPolicyReplay"
    }
    '''

    def __init__(self, memory_spec, algorithm, body):
        super(OnPolicyReplay, self).__init__(memory_spec, algorithm, body)
        # NOTE for OnPolicy replay, frequency = episode; for other classes below frequency = frames
        util.set_attr(self, self.agent_spec['algorithm'], ['training_frequency'])
        self.state_buffer = deque(maxlen=0)  # for API consistency
        # Don't want total experiences reset when memory is
        self.is_episodic = True
        self.total_experiences = 0
        self.warn_size_once = ps.once(lambda msg: logger.warn(msg))
        self.reset()

    @lab_api
    def reset(self):
        '''Resets the memory. Also used to initialize memory vars'''
        self.data_keys = ['states', 'actions', 'rewards', 'next_states', 'dones', 'priorities']
        for k in self.data_keys:
            setattr(self, k, [])
        self.cur_epi_data = {k: [] for k in self.data_keys}
        self.most_recent = [None] * len(self.data_keys)
        self.true_size = 0  # Size of the current memory

        self.state_buffer.clear()
        for _ in range(self.state_buffer.maxlen):
            self.state_buffer.append(np.zeros(self.body.state_dim))

    @lab_api
    def update(self, action, reward, state, done):
        '''Interface method to update memory'''
        self.base_update(action, reward, state, done)
        if not np.isnan(reward):  # not the start of episode
            self.add_experience(self.last_state, action, reward, state, done)
        self.last_state = state

    def add_experience(self, state, action, reward, next_state, done, priority=1):
        '''Interface helper method for update() to add experience to memory'''
        self.most_recent = [state, action, reward, next_state, done, priority]
        for idx, k in enumerate(self.data_keys):
            self.cur_epi_data[k].append(self.most_recent[idx])
        # If episode ended, add to memory and clear cur_epi_data
        if done:
            for k in self.data_keys:
                getattr(self, k).append(self.cur_epi_data[k])
            self.cur_epi_data = {k: [] for k in self.data_keys}
            # If agent has collected the desired number of episodes, it is ready to train
            # length is num of epis due to nested structure
            if len(self.states) == self.body.agent.algorithm.training_frequency:
                self.body.agent.algorithm.to_train = 1
        # Track memory size and num experiences
        self.true_size += 1
        if self.true_size > 1000:
            self.warn_size_once('Large memory size: {}'.format(self.true_size))
        self.total_experiences += 1

    def get_most_recent_experience(self):
        '''Returns the most recent experience'''
        return self.most_recent

    def sample(self):
        '''
        Returns all the examples from memory in a single batch
        Batch is stored as a dict.
        Keys are the names of the different elements of an experience. Values are nested lists of the corresponding sampled elements. Elements are nested into episodes
        e.g.
        batch = {
            'states'     : [[s_epi1], [s_epi2], ...],
            'actions'    : [[a_epi1], [a_epi2], ...],
            'rewards'    : [[r_epi1], [r_epi2], ...],
            'next_states': [[ns_epi1], [ns_epi2], ...],
            'dones'      : [[d_epi1], [d_epi2], ...],
            'priorities' : [[p_epi1], [p_epi2], ...]}
        '''
        batch = {k: getattr(self, k) for k in self.data_keys}
        self.reset()
        return batch

class OnPolicySeqReplay(OnPolicyReplay):
    '''
    Same as OnPolicyReplay Memory but returns the last `seq_len` states and next_states for input to a recurrent network.
    Experiences with less than `seq_len` previous examples are padded with a 0 valued state and action vector.

    e.g. memory_spec
    "memory": {
        "name": "OnPolicySeqReplay"
    }
    * seq_len provided by net_spec
    '''

    def __init__(self, memory_spec, algorithm, body):
        self.seq_len = algorithm.net_spec['seq_len']
        super(OnPolicySeqReplay, self).__init__(memory_spec, algorithm, body)
        self.state_buffer = deque(maxlen=self.seq_len)

    def preprocess_state(self, state, append=True):
        '''Transforms the raw state into format that is fed into the network'''
        if append:
            assert id(state) != id(self.state_buffer[-1]), 'Do not append to buffer other than during action'
            self.state_buffer.append(state)
        processed_state = np.stack(self.state_buffer)
        return processed_state

    def sample(self):
        '''
        Returns all the examples from memory in a single batch
        Batch is stored as a dict.
        Keys are the names of the different elements of an experience. Values are nested lists of the corresponding sampled elements. Elements are nested into episodes
        states and next_states have are further nested into sequences containing the previous `seq_len` - 1 relevant states
        e.g.
        let s_seq_0 be [0, ..., s0] (zero-padded), s_seq_k be [s_{k-seq_len}, ..., s_k], so the states are nested for passing into RNN.
        batch = {
            'states'    : [
                [s_seq_0, s_seq_1, ..., s_seq_k]_epi_1,
                [s_seq_0, s_seq_1, ..., s_seq_k]_epi_2,
                ...]
            'actions'   : [[a_epi1], [a_epi2], ...],
            'rewards'   : [[r_epi1], [r_epi2], ...],
            'next_states: [
                [ns_seq_0, ns_seq_1, ..., ns_seq_k]_epi_1,
                [ns_seq_0, ns_seq_1, ..., ns_seq_k]_epi_2,
                ...]
            'dones'     : [[d_epi1], [d_epi2], ...],
            'priorities': [[p_epi1], [p_epi2], ...]}
        '''
        batch = {}
        batch['states'] = self.build_seqs(self.states)
        batch['actions'] = self.actions
        batch['rewards'] = self.rewards
        batch['next_states'] = self.build_seqs(self.next_states)
        batch['dones'] = self.dones
        batch['priorities'] = self.priorities
        self.reset()
        return batch

    def build_seqs(self, data):
        '''Construct the epi-nested-seq data for sampling'''
        all_epi_data_seq = []
        for epi_data in data:
            data_seq = []
            # make [0, ..., *epi_data]
            padded_epi_data = copy.deepcopy(epi_data)
            padding = np.zeros_like(epi_data[0])
            for i in range(self.seq_len - 1):
                padded_epi_data.insert(0, padding)
            # slide seqs and build for one epi
            for i in range(len(epi_data)):
                data_seq.append(padded_epi_data[i:i + self.seq_len])
            all_epi_data_seq.append(data_seq)
        return all_epi_data_seq

class OnPolicyBatchReplay(OnPolicyReplay):
    '''
    Same as OnPolicyReplay Memory with the following difference.

    The memory does not have a fixed size. Instead the memory stores data from N experiences, where N is determined by the user. After N experiences or if an episode has ended, all of the examples are returned to the agent to learn from.

    In contrast, OnPolicyReplay stores entire episodes and stores them in a nested structure. OnPolicyBatchReplay stores experiences in a flat structure.

    e.g. memory_spec
    "memory": {
        "name": "OnPolicyBatchReplay"
    }
    * batch_size is training_frequency provided by algorithm_spec
    '''

    def __init__(self, memory_spec, algorithm, body):
        super(OnPolicyBatchReplay, self).__init__(memory_spec, algorithm, body)
        self.is_episodic = False

    def add_experience(self, state, action, reward, next_state, done, priority=1):
        '''Interface helper method for update() to add experience to memory'''
        self.most_recent = [state, action, reward, next_state, done, priority]
        for idx, k in enumerate(self.data_keys):
            getattr(self, k).append(self.most_recent[idx])
        # Track memory size and num experiences
        self.true_size += 1
        if self.true_size > 1000:
            self.warn_size_once('Large memory size: {}'.format(self.true_size))
        self.total_experiences += 1
        # Decide if agent is to train
        if done or len(self.states) == self.body.agent.algorithm.training_frequency:
            self.body.agent.algorithm.to_train = 1

    def sample(self):
        '''
        Returns all the examples from memory in a single batch
        Batch is stored as a dict.
        Keys are the names of the different elements of an experience. Values are a list of the corresponding sampled elements
        e.g.
        batch = {
            'states'     : states,
            'actions'    : actions,
            'rewards'    : rewards,
            'next_states': next_states,
            'dones'      : dones,
            'priorities' : priorities}
        '''
        return super(OnPolicyBatchReplay, self).sample()

class OnPolicySeqBatchReplay(OnPolicyBatchReplay):
    '''
    Same as OnPolicyBatchReplay Memory but returns the last `seq_len` states and next_states for input to a recurrent network.
    Experiences with less than `seq_len` previous examples are padded with a 0 valued state and action vector.

    e.g. memory_spec
    "memory": {
        "name": "OnPolicySeqBatchReplay"
    }
    * seq_len provided by net_spec
    * batch_size is training_frequency provided by algorithm_spec
    '''

    def __init__(self, memory_spec, algorithm, body):
        self.is_episodic = False
        self.seq_len = algorithm.net_spec['seq_len']
        super(OnPolicySeqBatchReplay, self).__init__(memory_spec, algorithm, body)
        self.state_buffer = deque(maxlen=self.seq_len)

    def reset(self):
        '''Initializes the memory arrays, size and head pointer'''
        super(OnPolicySeqBatchReplay, self).reset()
        self.state_buffer.clear()
        for _ in range(self.state_buffer.maxlen):
            self.state_buffer.append(np.zeros(self.body.state_dim))

    def preprocess_state(self, state, append=True):
        '''Transforms the raw state into format that is fed into the network'''
        if append:
            assert id(state) != id(self.state_buffer[-1]), 'Do not append to buffer other than during action'
            self.state_buffer.append(state)
        processed_state = np.stack(self.state_buffer)
        return processed_state

    def sample(self):
        '''
        Returns all the examples from memory in a single batch
        Batch is stored as a dict.
        Keys are the names of the different elements of an experience. Values are a list of the corresponding sampled elements.
        states and next_states have are further nested into sequences containing the previous `seq_len` - 1 relevant states
        e.g.
        let s_seq_0 be [0, ..., s0] (zero-padded), s_seq_k be [s_{k-seq_len}, ..., s_k], so the states are nested for passing into RNN.
        batch = {
            'states'     : [[s_seq_0, s_seq_1, ..., s_seq_k]],
            'actions'    : actions,
            'rewards'    : rewards,
            'next_states': [[ns_seq_0, ns_seq_1, ..., ns_seq_k]],
            'dones'      : dones,
            'priorities' : priorities}
        '''
        batch = {}
        batch['states'] = self.build_seqs(self.states)
        batch['actions'] = self.actions
        batch['rewards'] = self.rewards
        batch['next_states'] = self.build_seqs(self.next_states)
        batch['dones'] = self.dones
        batch['priorities'] = self.priorities
        self.reset()
        return batch

    def build_seqs(self, data):
        '''Construct the seq data for sampling'''
        data_seq = []
        # make [0, ..., *data]
        padded_data = copy.deepcopy(data)
        padding = np.zeros_like(data[0])
        for i in range(self.seq_len - 1):
            padded_data.insert(0, padding)
        # slide seqs and build for one epi
        for i in range(len(data)):
            data_seq.append(padded_data[i:i + self.seq_len])
        return data_seq

max_stars_repo_name: cclauss/SLM-Lab

max_issues_repo_name: cclauss/SLM-Lab

max_forks_repo_name: cclauss/SLM-Lab


text:
import abc
try:
    import anndata
except ImportError:
    anndata = None
import numpy as np

from .external import _ModelGLM

class Model(_ModelGLM, metaclass=abc.ABCMeta):
    """
    Generalized Linear Model (GLM) with negative binomial noise.
    """

    def link_loc(self, data):
        return np.log(data)

    def inverse_link_loc(self, data):
        return np.exp(data)

    def link_scale(self, data):
        return np.log(data)

    def inverse_link_scale(self, data):
        return np.exp(data)

    @property
    def eta_loc(self) -> np.ndarray:
        eta = np.matmul(self.design_loc, self.a)
        if self.size_factors is not None:
            eta += self.size_factors
        eta = self.np_clip_param(eta, "eta_loc")
        return eta

    def eta_loc_j(self, j) -> np.ndarray:
        # Make sure that dimensionality of sliced array is kept:
        if isinstance(j, int) or isinstance(j, np.int32) or isinstance(j, np.int64):
            j = [j]
        eta = np.matmul(self.design_loc, self.a[:, j])
        if self.size_factors is not None:
            eta += self.size_factors
        eta = self.np_clip_param(eta, "eta_loc")
        return eta

    # Re-parameterizations:

    @property
    def mu(self) -> np.ndarray:
        return self.location

    @property
    def phi(self) -> np.ndarray:
        return self.scale

max_stars_repo_name: le-ander/batchglm

max_issues_repo_name: le-ander/batchglm

max_forks_repo_name: le-ander/batchglm