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

# * Copyright (c) 2021-2022. Authors: Ba Thien LE.
# *
# * Licensed under the Apache License, Version 2.0 (the "License");
# * you may not use this file except in compliance with the License.
# * You may obtain a copy of the License at
# *
# *      http://www.apache.org/licenses/LICENSE-2.0
# *
# * Unless required by applicable law or agreed to in writing, software
# * distributed under the License is distributed on an "AS IS" BASIS,
# * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# * See the License for the specific language governing permissions and
# * limitations under the License.

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

import cv2
import numpy as np
import torch
import os

from pathlib import Path
from PIL import Image

from shapely import wkt
from shapely.geometry import Polygon

from torchvision.transforms.functional import to_tensor, resize

from cytomine import CytomineJob
from cytomine.models import Annotation, AnnotationCollection

from src import NuClick, post_process

__author__ = "LE Ba Thien <ba.le@uliege.be>"

def to_uint8(array):
    if isinstance(array, torch.Tensor):
        return (array * 255).type(torch.uint8).numpy()

    if isinstance(array, np.ndarray):
        return (array * 255).astype(np.uint8)

def load_model(filename):
    model = NuClick()

    weight_path = os.path.join('/weights/', filename)

    model.load_state_dict(torch.load(weight_path, map_location='cpu'))
    model.eval()

    return model

def predict(model, image):
    shape = (image.shape[0], 2, image.shape[2], image.shape[3])
    signal = torch.zeros(shape, device=image.device)
    inputs = torch.cat([image, signal], dim=1)

    with torch.no_grad():
        preds = model(inputs)

    return preds

def main(argv):
    with CytomineJob.from_cli(argv) as cj:
        cj.job.update(progress=0, statusComment="Initialisation")

        # Create the image directory
        working_path = str(Path.home())
        images_path = os.path.join(working_path, 'images')
        os.makedirs(images_path, exist_ok=True)

        # Load the NuClick model
        model = load_model(f'weight-{cj.parameters.type}.pth')

        cj.job.update(progress=15, statusComment="Load the model")

        # Get the images of the project
        images = [
            int(id) for id in cj.parameters.cytomine_id_images.split(',')
        ]

        # Get the ROIs
        paths = []
        offsets = []
        for image in images:
            roi_annotations = AnnotationCollection(
                project=cj.parameters.cytomine_id_project,
                term=cj.parameters.roi_id_term,
                image=image,
                showWKT=True
            )
            roi_annotations.fetch()

            # Download the ROIs of the targeted project
            for roi in roi_annotations:
                filename = f'{roi.image}-{roi.id}.png'
                filepath = os.path.join(images_path, filename)
                roi.dump(dest_pattern=filepath)

                paths.append((roi.image, filepath))

                bbox = wkt.loads(roi.location).bounds
                offsets.append((bbox[0], bbox[3]))

        cj.job.update(progress=40, statusComment="Fetch the ROIs")

        # Predict the object inside the ROI
        annotations = AnnotationCollection()
        for (image_id, path), (min_x, min_y) in zip(paths, offsets):
            # Open the ROI image
            image = to_tensor(resize(Image.open(path), (512, 512)))

            # Predictions
            preds = predict(model, image.unsqueeze(0))

            masks = post_process(
                preds,
                min_size=cj.parameters.min_size,
                area_threshold=cj.parameters.area_threshold
            )
            masks = to_uint8(masks.squeeze(0).permute(1, 2, 0))

            contours, _ = cv2.findContours(
                masks,
                cv2.RETR_EXTERNAL,
                cv2.CHAIN_APPROX_SIMPLE
            )

            contours = map(np.squeeze, contours)
            for c in contours:
                # Add offset to have the correct coords
                c[:, 0] += int(min_x)
                c[:, 1] = int(min_y) - c[:, 1]

                mask = Polygon(c)
                annotations.append(
                    Annotation(
                        location=mask.wkt,
                        id_image=image_id,
                        id_project=cj.parameters.cytomine_id_project
                    )
                )

        cj.job.update(progress=90, statusComment="Predictions done")

        # Save all the annotations
        annotations.save()

        cj.job.update(progress=100, statusComment="Job terminated")

if __name__ == "__main__":
    import sys

    main(sys.argv[1:])


text:
    def test_divide_on_different_shapes(self):
        x = MaskedArray(np.arange(6, dtype=float))
        x.shape = (2, 3)
        y = MaskedArray(np.arange(3, dtype=float))

        z = x / y
        assert_masked_equal(z, MaskedArray([[X, 1., 1.], [X, 4., 2.5]]))

        z = x / y[None,:]
        assert_masked_equal(z, MaskedArray([[X, 1., 1.], [X, 4., 2.5]]))

        y = MaskedArray(np.arange(2, dtype=float))
        z = x / y[:, None]
        assert_masked_equal(z, MaskedArray([[X, X, X], [3., 4., 5.]]))

    def test_mixed_arithmetic(self):
        # Tests mixed arithmetics.
        na = np.array([1])
        ma = MaskedArray([1])
        assert_(isinstance(na + ma, MaskedArray))
        assert_(isinstance(ma + na, MaskedArray))

    def test_limits_arithmetic(self):
        tiny = np.finfo(float).tiny
        a = MaskedArray([tiny, 1. / tiny, 0.])
        assert_equal((a / 2).mask, [0, 0, 0])
        assert_equal((2 / a).mask, [0, 0, 1]) # changed from numpy

    def test_masked_singleton_arithmetic(self):
        # Tests some scalar arithmetics on MaskedArrays.
        # Masked singleton should remain masked no matter what
        xm = MaskedArray(0, mask=1)
        assert_((1 / MaskedArray(0)).mask)
        assert_((1 + xm).mask)
        assert_((-xm).mask)
        assert_(np.maximum(xm, xm).mask)
        assert_(np.minimum(xm, xm).mask)

    def test_scalar_arithmetic(self):
        # Next lines from ogigianl tests, no longer valid. XXX so this means .filled used to return a view?
        #x = MaskedArray(0, mask=0)
        #assert_equal(x.filled().ctypes.data, x.ctypes.data)

        # Make sure we don't lose the shape in some circumstances
        xm = MaskedArray((0, 0)) / 0.
        assert_equal(xm.shape, (2,))
        assert_equal(xm.mask, [1, 1])

    def test_mod(self):
        # Tests mod
        (x, y, a10, m1, m2, xm, ym, z, zm, xf) = self.d
        mx = MaskedArray(x)
        my = MaskedArray(y)
        assert_equal(np.mod(x, y), np.mod(mx, my).filled(np.nan))

        test = np.mod(ym, xm)
        assert_equal(test.mask, xm.mask | ym.mask)
        test = np.mod(xm, ym)
        assert_equal(test.mask, xm.mask | ym.mask | (ym == 0).filled(True))

    def test_numpyarithmetics(self):
        # Check that the mask is not back-propagated when using numpy functions
        a = MaskedArray([-1, 0, 1, 2, X])
        assert_masked_equal(np.log(a), MaskedArray([X, X, 0, np.log(2), X]))
        assert_equal(a.mask, [0, 0, 0, 0, 1])

    def test_ndarray_mask(self):
        # Check that the mask of the result is a ndarray (not a MaskedArray...)
        a = MaskedArray([-1, 0, 1, 2, X])
        test = np.sqrt(a)
        control = MaskedArray([X, 0, 1, np.sqrt(2), X])
        assert_equal(test, control)
        assert_equal(test.mask, control.mask)
        assert_(not isinstance(test.mask, MaskedArray))

    def test_inplace_division_scalar_float(self):
        # Test of inplace division
        (x, y, xm) = self.floatdata
        x = x.copy()
        xm = xm.copy()

        x /= 2.0
        assert_equal(x, y / 2.0)
        xm /= np.arange(10)
        assert_masked_equal(xm, MaskedArray([X, 1, X, 1, 1, 1, 1, 1, 1, 1.]))

    def test_inplace_division_array_float(self):
        # Test of inplace division
        (x, y, xm) = self.floatdata
        x = x.copy()
        xm = xm.copy()

        m = xm.mask
        a = MaskedArray(np.arange(10, dtype=float))
        a[-1] = X
        x /= a
        xm /= a
        assert_equal(x, y / a)
        assert_equal(xm, y / a)
        assert_equal(xm.mask, m | (a == 0).filled(True))

    def test_inplace_division_misc(self):
        x = MaskedArray([X, 1., 1., -2., pi / 2., 4., X, -10., 10., 1., 2., 3.])
        y = MaskedArray([5., 0., X, 2., -1., X, X, -10., 10., 1., 0., X])
        control = MaskedArray([X, X, X, -1, -pi/2, X, X, 1, 1, 1, X, X])

        z = xm / ym
        assert_masked_equal(z, control)

        xm = xm.copy()
        xm /= ym
        assert_masked_equal(xm, control)

    def test_inplace_floor_division_array_type(self):
        # Test of inplace division
        for t in self.othertypes:
            with warnings.catch_warnings(record=True) as w:
                warnings.filterwarnings("always")
                (x, y, xm) = (_.astype(t) for _ in self.uint8data)
                m = xm.mask
                a = MaskedArray(np.arange(10, dtype=t))
                a[-1] = X
                x //= a
                xm //= a
                assert_equal(x, y // a)
                assert_equal(xm, y // a)
                assert_equal(xm.mask, m | (a == t(0)).filled(True))

                assert_equal(len(w), 0, "Failed on type=%s." % t)

    def test_inplace_floor_division_array_type(self):
        # Test of inplace division
        for t in self.othertypes:
            with warnings.catch_warnings(record=True) as w:
                warnings.filterwarnings("always")
                (x, y, xm) = (_.astype(t) for _ in self.uint8data)
                m = xm.mask
                a = MaskedArray(np.arange(10, dtype=t))
                a[-1] = X
                x //= a
                xm //= a
                assert_equal(x, y // a)
                assert_equal(xm, y // a)
                assert_equal(xm.mask, m)

                assert_equal(len(w), 0, "Failed on type=%s." % t)

    def test_inplace_division_array_type(self):
        # Test of inplace division
        for t in self.othertypes:
            with suppress_warnings() as sup:
                sup.record(UserWarning)
                (x, y, xm) = (_.astype(t) for _ in self.uint8data)
                m = xm.mask
                a = MaskedArray(np.arange(10, dtype=t))
                a[-1] = X

                # May get a DeprecationWarning or a TypeError.
                #
                # This is a consequence of the fact that this is true divide
                # and will require casting to float for calculation and
                # casting back to the original type. This will only be raised
                # with integers. Whether it is an error or warning is only
                # dependent on how stringent the casting rules are.
                #
                # Will handle the same way.
                try:
                    x /= a
                    assert_equal(x, y / a)
                except (DeprecationWarning, TypeError) as e:
                    warnings.warn(str(e), stacklevel=1)
                try:
                    xm /= a
                    assert_equal(xm, y / a)
                    assert_equal(xm.mask, m | (a == t(0)).filled(True))
                except (DeprecationWarning, TypeError) as e:
                    warnings.warn(str(e), stacklevel=1)

                if issubclass(t, np.integer):
                    assert_equal(len(sup.log), 2, "Failed on type=%s." % t)
                else:
                    assert_equal(len(sup.log), 0, "Failed on type=%s." % t)

ext: py

reason: from numpy


text:
import os
from os import listdir
from os.path import isfile, join, splitext, isdir
from scipy.signal import butter, lfilter, filtfilt
import numpy as np
from sklearn import preprocessing

def butter_bandpass_filter(data, lowcut, highcut, fs, order=9):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq

    b, a = butter(order, [low, high], btype='bandpass',analog=False)
    y = lfilter(b, a, data)
    return y

def Normalization(X):
    X_max, X_min = X.max(), X.min()
    X = (X-X_min)/(X_max-X_min)
    return X

max_stars_repo_name: Lukious/bisat

max_issues_repo_name: Lukious/bisat

max_forks_repo_name: Lukious/bisat

reason: import numpy,from scipy


text:
import comet_ml
comet_ml.config.save(api_key="WRmA8ms9A78K85fLxcv8Nsld9")
from comet_ml.api import API
import numpy as np
#import pandas as pd
import matplotlib.pyplot as plt

api = API()
exp_urls_lr1 = {"run1":"yasmeenvh/growspaceenv-ppo/91df8e59877e484f96022b6b3548fcf8",
                "run2":"yasmeenvh/growspaceenv-ppo/a633ea5fc5064259a013889cf46dfe64",
                "run3":"yasmeenvh/growspaceenv-ppo/b9d798bcf5d545c884162ac3373c8d8e"}

# exp_urls_lr2 = {"run1":"https://www.comet.ml/yasmeenvh/growspaceenv-ppo/91df8e59877e484f96022b6b3548fcf8",
#                 "run2":"https://www.comet.ml/yasmeenvh/growspaceenv-ppo/a633ea5fc5064259a013889cf46dfe64",
#                 "run3":"https://www.comet.ml/yasmeenvh/growspaceenv-ppo/b9d798bcf5d545c884162ac3373c8d8e"}

comet_metrics = ["Reward Mean", "Reward Max", "Reward Min", "Reward Max", "Number of Mean Branches"]
def get_data(url_link, comet_metrics):
    min_numframes = 0000000
    max_numframes = 100000000
    reward_means = []
    for run, link in url_link.items():
        experiment = api.get(link)

        reward_mean = experiment.get_metrics(comet_metrics[0])[:10000]
        reward_mean = np.array(reward_mean).transpose()
        print(reward_mean)
        reward_mean = [r_m['metricValue'] for r_m in reward_mean]# = reward_mean[:,np.where(np.logical_and(reward_mean[0] >= 000000, reward_mean[0] < 1000000))[0]]
        reward_mean = [float(x) for x in reward_mean]
        reward_means.append(reward_mean)
        print(reward_mean)

    return reward_means

def plot(value_comet):
    rewards = np.mean(value_comet, axis=0)
    # print("these are av rewards:",av_oracle)
    # err_oracle = np.std(oracle_rewards, axis=0)
    #
    # av_random = np.mean(random_rewards, axis=0)
    err_value = np.std(value_comet, axis=0)
    x = np.linspace(0,1000000,10000)
    #
    fig, ax = plt.subplots(1)
    ax.plot(x, rewards, color = 'green')
    #
    #ax.plot(x, av_random, color = 'magenta')
    ax.fill_between(x, rewards+err_value, rewards-err_value, alpha=0.5, color='green')
    # ax.fill_between(x, av_random + err_random, av_random - err_random, alpha=0.5, color='magenta')
    ax.set_title(r'Average rewards')
    ax.set_ylabel("rewards")
    ax.set_xlabel("steps")
    ax.grid()
    plt.show()

if __name__ == "__main__":
    test = get_data(exp_urls_lr1, comet_metrics)
    print(len(test[0]))
    print(len(test[1]))
    print(len(test[2]))
    plot(test)

reason: import numpy


text:
# -*- coding: utf-8 -*-
"""
Created on Thu Aug 19 19:35:05 2021

@author: Gabriel
"""
import numpy as np

def floor(x: float) -> int:
    return int(x)

def ceil(x: float) -> int:
    return int(x) + 1

def frac(x: float) -> int:
    return x - int(x)

def fact(n: int) -> int:
    '''
    Return the factorial of a number

    Parameters
    ----------
    n : int
        Integer.

    Returns
    -------
    int
        Factorial of n.

    '''

    return 1 if n == 0 else n * fact(n-1)

def C(n: int, r: int) -> int:
    '''
    Number of combinations of n things taken r at a time

    Parameters
    ----------
    n : int
        DESCRIPTION.
    r : int
        DESCRIPTION.

    Returns
    -------
    int
        DESCRIPTION.

    '''

    assert n >= r, "n must be greater or equal than r"

    return fact(n) // (fact(r) * fact(n-r))

def P(n: int, r: int) -> int:
    '''
    Number of combinations of n things taken r at a time

    Parameters
    ----------
    n : int
        DESCRIPTION.
    r : int
        DESCRIPTION.

    Returns
    -------
    int
        DESCRIPTION.

    '''
    assert n >= r, "n must be greater or equal than r"

    return fact(n) // fact(n-r)

def summation(f, a: int, b: int) -> float:
    '''
    Cummulative sum of f.

    Parameters
    ----------
    f : function
        A numerical function.
    a : int
        Summation lower bound.
    b : int
        Summation upper bound.

    Returns
    -------
    float
        Cummulative value of f in the interval specified.

    '''

    s = 0
    for i in range(a, b):
        s += np.round(f(i), 4)

    return s

def integral(f, a: float, b: float, n: int = 1000) -> float:
    '''
    Integral of f.

    Parameters
    ----------
    f : function
        A numerical function.
    a : float
        Integral lower bound.
    b : float
        Integral upper bound.
    n : int
        Number of subdivisions of the interval.

    Returns
    -------
    float
        Integral of f from a to b.

    '''

    dx = (b - a) / n
    F = 0
    x = a
    for i in range(n + 1):
        F += np.round(f(x) * dx, 4)
        x += dx

    return F

def generate_range(a, b, n):
    '''
    Generates a n+1 numbers betwen a and b with a constant step.
    Almost equivalent to numpy linspace

    Parameters
    ----------
    a : int
        Lower value.
    b : int
        Upper value.
    n : int
        Size of the list. Must be even more increased as |b-a| gets bigger.

    Returns
    -------
    list
        List of n+1 numbers inside the [a, b] interval.

    '''
    temp = [i/n for i in range(n+1)]
    diff = (b - a)
    return [i * diff + a for i in temp]


text:
from math import isclose

import pandas as pd
import numpy as np

import matplotlib.pyplot as plt

import metpy.calc as mpcalc
import metpy.constants as c
from metpy.plots import add_metpy_logo, add_timestamp, SkewT
from metpy.units import units

from datetime import datetime

import warnings
warnings.filterwarnings("ignore", category=RuntimeWarning)
warnings.filterwarnings("ignore", category=FutureWarning)

def plot_skewt(df,plot_stability=True,plot_cin_cape=True,plot_indices=True,output_pdf=False,output_display=True):
    """
    Plots annotated Skew-T log-P diagram using metpy module, given a meteorological profile

    Required inputs:
        df (Pandas dataframe) = sounding data, containing columns for:
                                'pressure' [hPa], 'temperature' [˚C], 'dewpoint' [˚C],
                                wind 'speed' [kts] and 'direction' [˚]
    Optional inputs:         Default:
        plot_stability (bool) [True]  = Include/exclude shading that shows stability in each interval of profile
        plot_cin_cape  (bool) [True]  = Include/exclude shading that shows CIN and CAPE
        plot_indices   (bool) [True]  = Include/exclude values for key indices estimated from profile
        output_pdf     (bool) [False] = Produce a PDF file of the plot in current directory
        output_display (bool) [True]  = Generate a plot interactively on the screen

    Outputs:
        skew         (object)         = The generated plot as an object
    """
    ###############################################################
    ### Set up interctive display
    plt.close()
    if output_display:
        plt.ion()
    else:
        plt.ioff()

    ###############################################################
    ### Parse meteorological variables

    # We will pull the data out of the sounding data into individual variables and assign units.
    p = df['pressure'].values * units.hPa
    T = df['temperature'].values * units.degC
    Td = df['dewpoint'].values * units.degC
    wind_speed = df['speed'].values * units.knots
    wind_dir = df['direction'].values * units.degrees
    u, v = mpcalc.wind_components(wind_speed, wind_dir)

    ###############################################################
    ### Initialize plot

    # Create a new figure. The dimensions here give a good aspect ratio.
    fig = plt.figure(figsize=(14, 10))
    skew = SkewT(fig,rect=(0.05,0.1,0.7,0.8),rotation=45) # `rotation` sets the slope of the isotherms
###    skew = SkewT(fig, rotation=45) # `rotation` sets the slope of the isotherms
    text_edge = 0.76

    ###############################################################
    ### Block for calculating and plotting stability

    if plot_stability:

        # For showing stability categories - color choices
        c_stab = ["indigo","silver","cadetblue","turquoise","yellowgreen","gold","deeppink"]

        # Other profile quantities needed to estimate stability in each interval
        abs_tol = 0.2    # The ± temperature range for defining isothermal, adiabatic profile intervals
        t_inv = T[:len(p)-1] - T[1:]  # Negative is inversion
        t_dry = [np.nan] * (len(p)-1) # Negative is dry stable, positive unstable
        t_moi = [np.nan] * (len(p)-1) # Negative is stable
                                      #   positive while t_dry negative is conditionally unstable
                                      #   both positive is moist adiabatic

        for i in range(len(p)-1):     # Every interval between measurements in the sounding...
            t_dry[i] = mpcalc.dry_lapse(p[i+1],T[i],reference_pressure=p[i]) - T[i+1]
            t_moi[i] = mpcalc.moist_lapse(p[i+1:i+2],T[i],reference_pressure=p[i]) - T[i+1]

            t_type = [2] * (len(p)-1) # Codes for stability; default is generic "stable"

        # Here we use codes for profiles in each interval:
        #  0 = inversion
        #  1 = isothermal (±abs_tol)
        #  2 = stable
        #  3 = moist adiabatic (±abs_tol)
        #  4 = conditionally unstable
        #  5 = neutral or adiabatic (±abs_tol)
        #  6 = absolutely unstable
        for i in range(len(p)-1):     # Every interval between measurements in the sounding...
            if t_inv[i] <= 0:
                t_type[i] = 0         # inversion
            if t_dry[i] > 0:
                t_type[i] = 6         # unstable
            if t_moi[i].size == 1:
                if t_dry[i] < 0 and t_moi[i] > 0:
                    t_type[i] = 4     # conditionally unstable
                if isclose(t_moi[i].magnitude, 0, abs_tol=abs_tol):
                    t_type[i] = 3     # moist adiabatic
            if isclose(t_dry[i].magnitude, 0, abs_tol=abs_tol):
                t_type[i] = 5         # neutral or adiabatic
            if isclose(t_inv[i].magnitude, 0, abs_tol=abs_tol):
                t_type[i] = 1         # isothermal

        for i in range(len(p)-1):     # Plot stability shading in each interval
            y = [p[i],p[i+1]]
            x1 = [T[i],T[i+1]]
            x2 = [(50.0 * units.degC),(50.0 * units.degC)]
            skew.shade_area(y, x1, x2, which='both', facecolor=c_stab[t_type[i]])

        plt.figtext(text_edge, 0.68,"Unstable", ha='left', va='center',fontsize=20,c=c_stab[6])
        plt.figtext(text_edge, 0.71,"Neutral or Adiabatic", ha='left', va='center',fontsize=20,c=c_stab[5])
        plt.figtext(text_edge, 0.74,"Conditionally Unstable", ha='left', va='center',fontsize=20,c=c_stab[4])
        plt.figtext(text_edge, 0.77,"Moist Adiabatic", ha='left', va='center',fontsize=20,c=c_stab[3])
        plt.figtext(text_edge, 0.80,"Stable", ha='left', va='center',fontsize=20,c=c_stab[2])
        plt.figtext(text_edge, 0.83,"Isothermal", ha='left', va='center',fontsize=20,c=c_stab[1])
        plt.figtext(text_edge, 0.86,"Inversion", ha='left', va='center',fontsize=20,c=c_stab[0])

    ###############################################################
    ### Standard sounding plot

    skew.ax.set_ylim(1050, 100) # set pressure level range
    skew.ax.set_xlim(-40, 50)   # set temperature range on bottom of plot
    skew.plot(p, T, 'crimson')
    skew.plot(p, Td, 'teal')
    p_top = np.where(p.magnitude>=100.0)[0][-1] # truncate at top of plot
    skew.plot_barbs(p[:p_top], u[:p_top], v[:p_top])

    # Calculate LCL height and plot as black dot
    lcl_pressure, lcl_temperature = mpcalc.lcl(p[0], T[0], Td[0])
    skew.plot(lcl_pressure, lcl_temperature, 'ko', markerfacecolor='black')

    # Place markers at surface values of T, Td
    skew.plot(p[0], T[0], color='maroon', marker='x')
    skew.plot(p[0], Td[0], color='darkslategrey', marker='x')

    # Calculate full lifted parcel profile and add to plot as black line
    prof = mpcalc.parcel_profile(p, T[0], Td[0]).to('degC')
    skew.plot(p, prof, 'k', linewidth=1)

    # Mark the LFC, EL
    t_parcel = mpcalc.parcel_profile(p, T[0], Td[0])
    lfc_pressure,lfc_temperature = mpcalc.lfc(p, T, Td, t_parcel)
    skew.plot(lfc_pressure, lfc_temperature, 'wo', markerfacecolor='red')
    el_pressure,el_temperature = mpcalc.el(p, T, Td, which='most_cape')
    skew.plot(el_pressure, el_temperature, 'wo', markerfacecolor='sienna')

    # An example of a slanted line at constant T -- in this case the 0
    # isotherm
    skew.ax.axvline(0, color='grey', linestyle='-.', linewidth=1)

    ###############################################################
    ### Shade areas of CAPE and CIN, plot stats
    if plot_cin_cape:
        skew.shade_cin(p, T, t_parcel)
        skew.shade_cape(p, T, t_parcel)

    ###############################################################
    ### Plot stability indices, other statistics
    if plot_indices:
        c_cape,c_cin = mpcalc.cape_cin(p, T, Td, t_parcel, which_lfc='bottom', which_el='top')
        c_lcl_p,c_lcl_t = mpcalc.lcl(p[0], T[0], Td[0], max_iters=50, eps=1e-05)
        c_lcl_z = ((T[0]-Td[0]) * (125 * units.m / units.degK)).to_base_units()
        c_el_p,c_el_t = mpcalc.el(p, T, Td, which='most_cape')
        c_pw = mpcalc.precipitable_water(p,Td)
        c_lfc_p,c_lfc_t = mpcalc.lfc(p, T, Td, t_parcel)
        c_lfc_t = c_lfc_t.to(units.degC)

        if not np.isnan(c_el_p.magnitude):
            plt.figtext(text_edge, 0.47,"EL: "+('%.0f' % c_el_p.magnitude)+" hPa", ha='left', va='center',fontsize=20,c='sienna')
            plt.figtext(text_edge, 0.44,"EL: "+('%.1f' % c_el_t.magnitude)+"˚C", ha='left', va='center',fontsize=20,c='sienna')
        if not np.isnan(c_lfc_p.magnitude):
            plt.figtext(text_edge, 0.39,"LFC: "+('%.1f' % c_lfc_p.magnitude)+" hPa", ha='left', va='center',fontsize=20,c='tab:red')
            plt.figtext(text_edge, 0.36,"LFC: "+('%.1f' % c_lfc_t.magnitude)+"˚C", ha='left', va='center',fontsize=20,c='tab:red')
        plt.figtext(text_edge, 0.305,"CAPE: "+('%.0f' % c_cape.magnitude)+" J/kg", ha='left', va='center',fontsize=20,c='firebrick')
        plt.figtext(text_edge, 0.265,"CIN: "+('%.0f' % c_cin.magnitude)+" J/kg", ha='left', va='center',fontsize=20,c='blue')
        plt.figtext(text_edge, 0.21,"PW: "+('%.0f' % c_pw.magnitude)+" mm", ha='left', va='center',fontsize=20,c='lightseagreen')
        plt.figtext(text_edge, 0.16,"LCL: "+('%.0f' % c_lcl_p.magnitude)+" hPa", ha='left', va='center',fontsize=20,c='black')
        plt.figtext(text_edge, 0.13,"LCL: "+('%.1f' % c_lcl_t.magnitude)+"˚C", ha='left', va='center',fontsize=20,c='black')

    ###############################################################
    ### Plot the relevant thermodynamic lines
    skew.plot_dry_adiabats(alpha=0.25, color='orangered')
    skew.plot_moist_adiabats(alpha=0.25, color='tab:green')
    skew.plot_mixing_lines(pressure=np.arange(1050, 80, -20) * units.hPa,linestyle='dotted', color='tab:blue',
            mixing_ratio=np.array([2.5e-6,1e-5,3e-5,8e-5,2e-4,5e-4,0.001,0.002,0.004,0.007,0.01,0.016,0.024,0.032]).reshape(-1, 1))

    plt.tick_params(axis = 'y', which = 'major', labelsize = 16)
    plt.tick_params(axis = 'x', which = 'major', labelsize = 16, labelrotation=45)

    # Add the timestamp for the data to the plot
    s_site = df['station'][0]
    s_time = '{dt:%Y%m%d_%H%M}'.format(dt=df['time'][0])
    add_timestamp(skew.ax, datetime.strptime(s_time, '%Y%m%d_%H%M'), pretext='Valid: ', y=1.02, x=-0.1, ha='left', fontsize=17)
    skew.ax.set_title(s_site,fontsize=28,x=0.66)

    ###############################################################
    ### Produce a PDF file of the plot
    if output_pdf:
        filename = s_site+"_"+s_time+".pdf"
        plt.savefig(filename)

    return skew


text:
import NumberToPattern
import numpy as np
import PatternToNumber

def computing_frequencies(text, k):
	m = 4**k
	frequency_array = np.zeros(shape = (1, m), dtype = int)

	for i in range(len(text) - k + 1):
		pattern = text[i:i+k]
		j = PatternToNumber.pattern_to_number(pattern)
		frequency_array[0 , j] = frequency_array[0, j] + 1
	return frequency_array.flatten()

def faster_frequent_word(text, k):
	frequent_words = set()
	frequency_array = computing_frequencies(text, k).tolist()
	max_count = max(frequency_array)
	for i in range(4**k):
		if frequency_array[i] == max_count:
			pattern = NumberToPattern.number_to_pattern(i, k)
			frequent_words.add(pattern)
	return frequent_words

def results(text, k):
	result = faster_frequent_word(text, k)
	print(' '.join(sorted(result)))


text:
import tensorflow as tf
import numpy as np
import time
import os
from tqdm import tqdm

# Import matplotlib libraries

import util as ut
import save_utils as su
import cropping as cr

def movenet(pic_paths, pic_dir_names, output_path, pose_type, model_type="li"):
    # get input folder name to create output file name
    # input_dir_name = os.path.basename(pic_path).split("./", 1)[0]
    # start_date = datetime.now().strftime("--%H-%M--%d-%m-%Y")
    # output_pic_dir_name = args.output + "/" + input_dir_name + start_date
    output_csv_dir_names = [os.path.join(
        output_path, pdn)+".csv" for i, pdn in enumerate(pic_dir_names)]
    # output_csv_dir_name = os.path.join(output_path, pic_dir_name)+".csv"

    movenet, input_size = choose_model(model_type)

    for pp, ocsv, pose in tqdm(zip(pic_paths, output_csv_dir_names, pose_type), desc="VID", ascii=True, total=len(pic_paths)):
        filenames = ut.get_filenames(0, pp)
        n_images = len(filenames)

        initial_image = tf.io.read_file(filenames[0])
        initial_image = tf.image.decode_jpeg(initial_image)
        image_height, image_width, _ = initial_image.shape
        crop_region = cr.init_crop_region(image_height, image_width)

        # output_images = []
        start = time.time()

        # PREPARE CSV FILE
        csv_file = open(ocsv, "ab")
        np.savetxt(csv_file, su.KEYPOINT_LABELS, delimiter=",", fmt="%s")

        for fname in tqdm(filenames, desc="PIC", ascii=True, total=n_images):
            # Load the input image.
            image = tf.io.read_file(fname)
            image = tf.image.decode_jpeg(image)
            keypoints_with_scores = cr.run_inference(
                movenet, image[:, :, :], crop_region,
                crop_size=[input_size, input_size])

            write = np.hstack(
                [fname, pose, keypoints_with_scores.squeeze().T.flatten()]).reshape([1, 53])
            np.savetxt(csv_file, write, delimiter=",", fmt='%s')

            # output_images.append(du.draw_prediction_on_image(
            #     image[:, :, :].numpy().astype(np.int32),
            #     keypoints_with_scores, crop_region=None,
            #     close_figure=True, output_image_height=image_height))
            crop_region = cr.determine_crop_region(
                keypoints_with_scores, image_height, image_width)

        print(': : : Average FPS:', n_images / (time.time() - start))
        csv_file.close()

    # output = np.stack(output_images, axis=0)
    # du.to_gif(output, fps=10)
    return output_csv_dir_names

def choose_model(model_name="movenet_lightning"):
    if (model_name == "movenet_lightning" or model_name == 'li'):
        module = tf.saved_model.load(
            "./models/lightning")
        input_size = 192
    elif (model_name == "movenet_thunder" or model_name == 'th'):
        module = tf.saved_model.load("./models/thunder")
        input_size = 256
    else:
        raise ValueError("Unsupported model name: %s" % model_name)

    return module.signatures['serving_default'], input_size

max_stars_repo_name: Shandelier/py-tf-movenet-fork

max_issues_repo_name: Shandelier/python-movenet

max_forks_repo_name: Shandelier/python-movenet


text:
import csv
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns

from matplotlib.ticker import FormatStrFormatter
from tensorflow.keras import backend as K

def load_data_from_csv(filename, shuffle=False, only_positive=False, multichannel=False, magRange=[20, 26]):
    """
    ['Label', 'Angle', 'Magnitude', 'Length', 'Pixel data']
    If only_positive = True, loads only positive training examples.
    """

    with open(filename) as file:
        csv_reader = csv.reader(file, delimiter=',', quotechar='"', quoting=csv.QUOTE_MINIMAL)
        first_line = True
        temp_images = []
        temp_labels = []
        width = 0
        skipped = 0
        for row in csv_reader:
            if first_line:
                first_line = False
                imageIndexStart = len(row) - 1
            elif ((only_positive is True and int(row[0]) == 1) or only_positive is False
                  ) and (magRange[0] <= float(row[6]) <= magRange[1] or float(row[6]) == 0):
                temp_labels.append(row[0:8])
                image_data = np.asarray(row[imageIndexStart:], dtype="uint32")
                if width == 0:
                    if multichannel is False:
                        width = int(np.sqrt(len(image_data)))
                    else:
                        width = int(np.sqrt(len(image_data)/4))
                if multichannel is False:
                    image_data_array = np.array_split(image_data, width)
                else:
                    image_data_array = image_data.reshape(width, width, 4)
                    # Converts image from channel first to channel last shape
                    #image_data_array = np.moveaxis(image_data_array, 0, 2)
                temp_images.append(image_data_array)
            else:
                #print("passing", only_positive, only_positive is False, ((only_positive is True and int(row[0]) == 1) or only_positive is False))
                skipped += 1
                pass
        labels = np.array(temp_labels).astype('float')
        images = np.array(temp_images).astype('uint32')
        if only_positive:
            print("  Only_positive mode:", only_positive, "  Magnitude range:", magRange, "  Skipped examples:", skipped)
    if shuffle:
        images, labels = unison_shuffled_copies(images, labels)

    return images, labels

def unison_shuffled_copies(a, b):
    """
    Shuffles two arrays the same way.
    """

    assert len(a) == len(b)
    p = np.random.permutation(len(a))

    return a[p], b[p]

def custom_YOLO_loss(y_true, y_pred):
    """
    Loss function for combined label and streak end coordinate prediction
    """

    # Weights for coordinates and labels
    params = [5.0, 0.5]

    lambda_coord = y_true[..., 0] * params[0]
    lambda_noobj = K.abs(y_true[..., 0] - y_pred[..., 0]) * params[1]

    lp = K.sum(K.square(y_true[..., 0] - y_pred[..., 0]) * lambda_noobj)
    lx1 = K.sum(K.square(y_pred[..., 1] - y_true[..., 1]) * lambda_coord)
    ly1 = K.sum(K.square(y_pred[..., 2] - y_true[..., 2]) * lambda_coord)
    lx2 = K.sum(K.square(y_pred[..., 3] - y_true[..., 3]) * lambda_coord)
    ly2 = K.sum(K.square(y_pred[..., 4] - y_true[..., 4]) * lambda_coord)

    return lp + lx1 + ly1 + lx2 + ly2

def f1_metric(y_true, y_pred):
    """
    F1 metric that combines precision and recall
    """

    true_positives = K.sum(K.round(K.clip(y_true * y_pred, 0, 1)))
    possible_positives = K.sum(K.round(K.clip(y_true, 0, 1)))
    predicted_positives = K.sum(K.round(K.clip(y_pred, 0, 1)))
    precision = true_positives / (predicted_positives + K.epsilon())
    recall = true_positives / (possible_positives + K.epsilon())
    f1_val = 2 * (precision * recall) / (precision + recall + K.epsilon())

    return f1_val

def analyze_5unit_errors(predictions, Y_test, image_size=20):
    """
    Compares predictions of labels and coordinates to ground truth
    """

    predictionsLabel = predictions[:, 0]
    predictionsPixel = predictions[:, 1:] * image_size
    Y_testLabel = Y_test[:, 0]
    Y_testPixel = Y_test[:, 1:] * image_size

    pixErrorList = []
    tp = 0
    tn = 0
    fp = 0
    fn = 0
    for i in range(len(predictionsPixel)):
        trueLabel = Y_testLabel[i]
        predictedLabel = round(predictionsLabel[i])
        if trueLabel == 1 and predictedLabel == 1:
            tp += 1
            pixError = abs(predictionsPixel[i] - Y_testPixel[i])
            pixErrorList.append(pixError)
        elif trueLabel == 0 and predictedLabel == 0:
            tn += 1
        elif trueLabel == 0 and predictedLabel == 1:
            fp += 1
        elif trueLabel == 1 and predictedLabel == 0:
            fn += 1
        else:
            pass
    accuracy = (tp + tn) / len(predictionsLabel) * 100
    precision = tp / (tp + fp) * 100
    recall = tp / (tp + fn) * 100
    pixErrorList = np.asarray(pixErrorList)
    medianError = np.median(pixErrorList, axis=0)
    meanError = np.mean(pixErrorList, axis=0)
    medianErrorTotal = np.median(pixErrorList)
    meanErrorTotal = np.mean(pixErrorList)

    print("\nClassification accuracy, precision, recall:", accuracy, precision, recall)
    print("TP, TN, FP, FN:", tp, tn, fp, fn)
    print("Median coordinate error (x1, y1, x2, y2) [pixels]:", medianError)
    print("Mean coordinate error (x1, y1, x2, y) [pixels]:", meanError)
    print("Median total coordinate error [pixels]:", medianErrorTotal)
    print("Mean total coordinate error [pixels]:", meanErrorTotal)

    return

def create_histogram_2d(predictions, Y, binsMag, binsLength):
    """
    Prepares results for plotting with plot_results_heatmap.
    """

    minMag = 20
    maxMag = 26

    gtMags = Y[:, 6]
    gtLengths = Y[:, 7]
    predLabels = predictions[:, 0]

    gtLengths = convert_pixels_to_arcsec_per_h(gtLengths)

    minLength = np.min(gtLengths)
    maxLength = np.max(gtLengths)

    range2d = [[minMag, maxMag], [minLength, maxLength]]
    bins = [binsMag, binsLength]

    tptnLengths = []
    tptnMags = []
    for i in range(len(Y)):
        if int(np.round(predLabels[i])) == int(Y[i][0]):
            tptnLengths.append(Y[i][7])
            tptnMags.append(Y[i][6])

    tptnLengths = np.asarray(tptnLengths)
    tptnMags = np.asarray(tptnMags)
    tptnLengths = convert_pixels_to_arcsec_per_h(tptnLengths)

    histGT, xEdgesGT, yEdgesGT = np.histogram2d(gtMags, gtLengths, bins=bins, range=range2d, normed=None, weights=None, density=None)
    histTPTN, xEdgesTPTN, yEdgesTPTN = np.histogram2d(tptnMags, tptnLengths, bins=bins, range=range2d, normed=None, weights=None, density=None)

    print("\n2D histogram counts")
    print("Mean count", np.mean(histGT))
    print("Min count", np.min(histGT))
    print("Length range:", minLength, maxLength)

    histResult = histTPTN / histGT * 100

    dataframe2D = pd.DataFrame(histResult, columns=yEdgesGT[:-1], index=xEdgesGT[:-1])

    return dataframe2D

def convert_pixels_to_arcsec_per_h(pixels):
    """
    Converts streak length in pixels to arcsec/h.
    Values are specific to ESA Euclid mission.
    """

    coefficient = 565/3600*10
    arcsecPerH = pixels/coefficient

    return arcsecPerH

def plot_results_heatmap(dataframe2D, binsMag, title, savepdf=True):
    """
    Plots the recall (completeness of predictions as a heatmap).
    Saves the heatmap to a pdf.
    """
    cmap = "PRGn"
    xticklabels = np.append(dataframe2D.columns.values.round(0).astype(int), 80)
    yticklabels = np.append(dataframe2D.index.values.round(2), 26.0)
    fig4, ax4 = plt.subplots(figsize=(10, 6), dpi=100)
    sns.heatmap(dataframe2D,
                ax=ax4,
                cmap=cmap,
                annot=True,
                fmt='3.0f',
                cbar_kws={'label': 'Completeness [%]'},
                xticklabels=xticklabels,
                yticklabels=yticklabels,
                vmin=0,
                vmax=100,
                linewidths=0.5)
    ax4.set_ylim(0, binsMag)
    ax4.set_title(title)
    ax4.set_xlabel(r'Sky motion [$\rm arcsec\,h^{-1}$]')
    ax4.set_ylabel('Magnitude')
    ax4.yaxis.set_major_formatter(FormatStrFormatter('%5.2f'))
    yticks = ax4.get_yticks() - 0.5
    ax4.set_yticks(yticks)
    xticks = ax4.get_xticks() - 0.5
    ax4.set_xticks(xticks)
    ax4.set_xticklabels(ax4.get_xticklabels(), rotation=0)  # , horizontalalignment='right')
    ax4.set_yticklabels(ax4.get_yticklabels(), rotation=0)
    #PRGn
    #PuBuGn
    fig4.tight_layout()
    if savepdf:
        fig4.savefig('Plot_2D_histogram_CNN.pdf', dpi=300, format="pdf")
    return

max_stars_repo_name: MikkoPontinen/EuclidCNN

max_issues_repo_name: MikkoPontinen/EuclidCNN

max_forks_repo_name: MikkoPontinen/EuclidCNN


text:
#Marcus Keil and Jon Holdship 13/03/2020
#Examples of a simple grid of models run in parallel
from __future__ import print_function
import uclchem
import numpy as np
import pandas as pd
from multiprocessing import Pool
import time

#uclchem general takes a dictionary of parameters where outSpecies=number of outspecies
#and a string outspeciesin which is a delimited list of species
#this wrapper helps Pool to work and also converts a delimited list in the dictionary into a number
#and send list separately as required.
def run_uclchem(param_dict):
    outSpecies = (param_dict['outSpecies'])
    param_dict['outSpecies'] = len(outSpecies.split())

    #this will run uclchem to final time step and return an array of abundances
    #arrays cannot be variably sized with f2py so you need to set length of output array in src/wrap.f90
    abunds=uclchem.wrap.run_model_for_abundances(dictionary=param_dict, outspeciesin=outSpecies)

    #altenatively, you can run uclchem to file as you normally would
    #just put columnfile and/or output file in the dictionary
    #uclchem.wrap.run_model_to_file(dictionary, outSpeciesIn)
    return abunds

#basic set of parameters we'll use for this grid.
ParameterDictionary = {"phase": 1, "switch": 0, "collapse": 0, "readAbunds": 0, "writeStep": 1,"fn": 6.1e-9,
                       "outSpecies": 'SO H3O+ NH3 HNC', "initialDens": 1e4,"finalTime":2.0e6,"baseAv":10}

# This part can be substituted with any choice of grid
# here we just combine various initial and final densities into an easily iterable array
initialTempGrid = np.linspace(50, 300, 6)
crgrid = np.logspace(1, 5, 5)
parameterSpace = np.asarray(np.meshgrid(initialTempGrid, crgrid)).reshape(2, -1)

print(parameterSpace.shape)

#then we loop through parameters, update parameter dictionary and run
#however, to use Pool we just store dictionaries for now
models=[]

#we'll number our models so store the parameters we used in a list
for k in range(np.shape(parameterSpace)[1]):
    paramDict=ParameterDictionary.copy()
    paramDict["initialTemp"] = parameterSpace[0,k]
    paramDict["zeta"] = parameterSpace[1,k]
    models.append(paramDict)

m=pd.DataFrame(models)
m.to_csv("models.csv",index=False)
#use pool.map to run each dictionary throuh our helper function
start=time.time()
pool=Pool(6)
result=pool.map(run_uclchem,models)
result=np.asarray(result)
pool.close()
pool.join()
df=pd.DataFrame({"Temperature":parameterSpace[0,:],"CR":parameterSpace[1,:],
                "SO":result[:,0],"H3O+":result[:,1],"NH3":result[:,2],"HNC":result[:,3]})

df.to_csv("grid.csv",index=False)
end=time.time()
end=(end-start)/60.0
print(f"grid in {end} minutes")

max_stars_repo_name: NormanRH/UCLCHEM

max_issues_repo_name: NormanRH/UCLCHEM

max_forks_repo_name: NormanRH/UCLCHEM


text:
"""
Implement one-hot encoder.
"""

from typing import Any, List
import dask.dataframe as dd
import numpy as np

class OneHotEncoder:
    """One-hot encoder for encoding categorical values
    Attributes:
        name
            Name of encoder
        unique_list
            Unique categorical values in provided data columns
        unique_num
            Number of unique categorical values in provided data columns
    """

    def __init__(self) -> None:
        """
        This function initiate numerical scaler.
        """

        self.name = "OneHotEncoder"
        self.unique_list = np.zeros(1)
        self.unique_num = 0

    def fit(self, col_df: dd.Series) -> Any:
        """
        Extract unique categorical values for one-hot encoder according to the provided column.

        Parameters
        ----------
        col_df
            Provided data column.
        """

        self.unique_list = col_df.drop_duplicates().values
        self.unique_num = col_df.drop_duplicates().count()
        return self

    def transform(self, col_df: dd.Series) -> dd.Series:
        """
        Transform the provided data column with the extracted unique values.

        Parameters
        ----------
        col_df
            Provided data column.
        """

        result = col_df.map(self.compute_val)
        return result

    def fit_transform(self, col_df: dd.Series) -> dd.Series:
        """
        Extract unique categorical values for one-hot encoder according to the data column.
        Transform the data column with the extracted unique values.

        Parameters
        ----------
        col_df
            Data column.
        """

        return self.fit(col_df).transform(col_df)

    def compute_val(self, val: str) -> List[float]:
        """
        Compute one-hot encoding of provided value.

        Parameters
        ----------
        val
            Value should be transferred to one-hot encoding.
        """
        temp_result = np.zeros(len(self.unique_list))
        idx = self.unique_list.tolist().index(val)
        temp_result[idx] = 1
        result: List[float] = temp_result.tolist()
        return result


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

import h5py
import numpy as np
import os
import sys
from openEPhys_DACQ.HelperFunctions import tetrode_channels, channels_tetrode, closest_argmin
from pprint import pprint
from copy import copy
import argparse
import importlib
from tqdm import tqdm

def OpenEphys_SamplingRate():
    return 30000

def bitVolts():
    return 0.195

def spike_waveform_leftwards_shift():
    """Returns the leftwards shift of waveforms from detection point in seconds.

    :return:
    :rtype: float
    """
    return 6 * (1.0 / OpenEphys_SamplingRate())

def get_filename(path):
    if not os.path.isfile(path):
        return os.path.join(path, 'experiment_1.nwb')
    else:
        return path

def delete_path_in_file(filename, path):
    with h5py.File(filename, 'r+') as h5file:
        del h5file[path]

def get_recordingKey(filename):
    with h5py.File(filename, 'r') as h5file:
        return list(h5file['acquisition']['timeseries'].keys())[0]

def get_all_processorKeys(filename):
    with h5py.File(filename, 'r') as h5file:
        return list(h5file['acquisition']['timeseries'][get_recordingKey(filename)]['continuous'].keys())

def get_processorKey(filename):
    return get_all_processorKeys(filename)[0]

def get_all_processor_paths(filename):
    return ['/acquisition/timeseries/' + get_recordingKey(filename)
            + '/continuous/' + processorKey
            for processorKey in get_all_processorKeys(filename)]

def get_processor_path(filename):
    return '/acquisition/timeseries/' + get_recordingKey(filename) \
           + '/continuous/' + get_processorKey(filename)

def check_if_open_ephys_nwb_file(filename):
    """
    Returns True if processor path can be identified
    in the file and False otherwise.
    """
    try:
        processor_path = get_processor_path(filename)
        with h5py.File(filename, 'r') as h5file:
            return processor_path in h5file
    except:
        return False

def get_downsampled_data_paths(filename):
    """
    Returns paths to downsampled data in NWB file.

    :param filename: path to NWB file
    :type filename: str
    :return: paths
    :rtype: dict
    """
    processor_path = get_processor_path(filename)
    return {'tetrode_data': processor_path + '/downsampled_tetrode_data/',
            'aux_data': processor_path + '/downsampled_AUX_data/',
            'timestamps': processor_path + '/downsampled_timestamps/',
            'info': processor_path + '/downsampling_info/'}

def check_if_downsampled_data_available(filename):
    """
    Checks if downsampled data is available in the NWB file.

    :param filename: path to NWB file
    :type filename: str
    :return: available
    :rtype: bool
    """
    paths = get_downsampled_data_paths(filename)
    with h5py.File(filename, 'r') as h5file:
        # START Workaround for older downsampled datasets
        if '/acquisition/timeseries/recording1/continuous/processor102_100/tetrode_lowpass' in h5file:
            return True
        # END Workaround for older downsampled datasets
        for path in [paths[key] for key in paths]:
            if not (path in h5file):
                return False
        if h5file[paths['tetrode_data']].shape[0] == 0:
            return False
        if h5file[paths['tetrode_data']].shape[0] \
                != h5file[paths['timestamps']].shape[0] \
                != h5file[paths['aux_data']].shape[0]:
            return False

    return True

def get_raw_data_paths(filename):
    """
    Returns paths to downsampled data in NWB file.

    :param filename: path to NWB file
    :type filename: str
    :return: paths
    :rtype: dict
    """
    processor_path = get_processor_path(filename)
    return {'continuous': processor_path + '/data',
            'timestamps': processor_path + '/timestamps'}

def check_if_raw_data_available(filename):
    """
    Returns paths to raw data in NWB file.

    :param filename:
    :type filename: str
    :return: paths
    :rtype: dict
    """
    paths = get_raw_data_paths(filename)
    if all([check_if_path_exists(filename, paths[key]) for key in paths]):
        return True
    else:
        return False

def save_downsampling_info_to_disk(filename, info):
    # Get paths to respective dataset locations
    paths = get_downsampled_data_paths(filename)
    # Ensure dictionary fields are in correct format
    info = {'original_sampling_rate': np.int64(info['original_sampling_rate']),
            'downsampled_sampling_rate': np.int64(info['downsampled_sampling_rate']),
            'downsampled_channels': np.array(info['downsampled_channels'], dtype=np.int64)}
    # Write data to disk
    with h5py.File(filename, 'r+') as h5file:
        recursively_save_dict_contents_to_group(h5file, paths['info'], info)

def save_downsampled_data_to_disk(filename, tetrode_data, timestamps, aux_data, info):
    # Get paths to respective dataset locations
    paths = get_downsampled_data_paths(filename)
    # Write data to disk
    save_downsampling_info_to_disk(filename, info)
    with h5py.File(filename, 'r+') as h5file:
        h5file[paths['tetrode_data']] = tetrode_data
        h5file[paths['timestamps']] = timestamps
        h5file[paths['aux_data']] = aux_data

def delete_raw_data(filename, only_if_downsampled_data_available=True):
    if only_if_downsampled_data_available:
        if not check_if_downsampled_data_available(filename):
            print('Warning', 'Downsampled data not available in NWB file. Raw data deletion aborted.')
            return None
    if not check_if_raw_data_available(filename):
        print('Warning', 'Raw data not available to be deleted in: ' + filename)
    else:
        raw_data_paths = get_raw_data_paths(filename)
        with h5py.File(filename,'r+') as h5file:
            for path in [raw_data_paths[key] for key in raw_data_paths]:
                del h5file[path]

def repack_NWB_file(filename, replace_original=True, check_validity_with_downsampled_data=True):
    # Create a repacked copy of the file
    os.system('h5repack ' + filename + ' ' + (filename + '.repacked'))
    # Check that the new file is not corrupted
    if check_validity_with_downsampled_data:
        if not check_if_downsampled_data_available(filename):
            raise Exception('Downsampled data cannot be found in repacked file. Original file not replaced.')
    # Replace original file with repacked file
    if replace_original:
        os.system('mv ' + (filename + '.repacked') + ' ' + filename)

def repack_all_nwb_files_in_directory_tree(folder_path, replace_original=True,
                                           check_validity_with_downsampled_data=True):
    # Commence directory walk
    for dir_name, subdirList, fileList in os.walk(folder_path):
        for fname in fileList:
            fpath = os.path.join(dir_name, fname)
            if fname == 'experiment_1.nwb':
                print('Repacking file {}'.format(fpath))
                repack_NWB_file(fpath, replace_original=replace_original,
                                check_validity_with_downsampled_data=check_validity_with_downsampled_data)

def list_AUX_channels(filename, n_tetrodes):
    data = load_continuous(filename)
    n_channels = data['continuous'].shape[1]
    data['file_handle'].close()
    aux_chan_list = range(n_tetrodes * 4 - 1, n_channels)

    return aux_chan_list

def load_continuous(filename):
    # Load data file
    h5file = h5py.File(filename, 'r')
    # Load timestamps and continuous data
    recordingKey = get_recordingKey(filename)
    processorKey = get_processorKey(filename)
    path = '/acquisition/timeseries/' + recordingKey + '/continuous/' + processorKey
    if check_if_path_exists(filename, path + '/data'):
        continuous = h5file[path + '/data'] # not converted to microvolts!!!! need to multiply by 0.195
        timestamps = h5file[path + '/timestamps']
        data = {'continuous': continuous, 'timestamps': timestamps, 'file_handle': h5file}
    else:
        data = None

    return data

def load_raw_data_timestamps_as_array(filename):
    data = load_continuous(filename)
    timestamps = np.array(data['timestamps']).squeeze()
    data['file_handle'].close()

    return timestamps

def load_data_columns_as_array(filename, data_path, first_column, last_column):
    """
    Loads a contiguous columns of dataset efficiently from HDF5 dataset.
    """
    with h5py.File(filename, 'r') as h5file:
        data = h5file[data_path]
        data = h5file[data_path][:, first_column:last_column]

    return data

def load_data_as_array(filename, data_path, columns):
    """
    Fast way of reading a single column or a set of columns.

    filename - str - full path to file
    columns  - list - column numbers to include (starting from 0).
               Single column can be given as a single list element or int.
               Columns in the list must be in sorted (ascending) order.
    """
    # Make columns variable into a list if int given
    if isinstance(columns, int):
        columns = [columns]
    # Check that all elements of columns are integers
    if isinstance(columns, list):
        for column in columns:
            if not isinstance(column, int):
                raise ValueError('columns argument must be a list of int values.')
    else:
        raise ValueError('columns argument must be list or int.')
    # Check that column number are sorted
    if sorted(columns) != columns:
        raise ValueError('columns was not in sorted (ascending) order.')
    # Check that data is available, otherwise return None
    if not check_if_path_exists(filename, data_path):
        raise ValueError('File ' + filename + '\n'
                         + 'Does not contain path ' + data_path)
    # Find contiguous column groups
    current_column = columns[0]
    column_groups = [current_column]
    for i in range(1, len(columns)):
        if (columns[i] - columns[i - 1]) == 1:
            column_groups.append(current_column)
        else:
            column_groups.append(columns[i])
            current_column = columns[i]
    # Find start and end column numbers for contiguous groups
    column_ranges = []
    for first_channel in sorted(set(column_groups)):
        last_channel = first_channel + column_groups.count(first_channel)
        column_ranges.append((first_channel, last_channel))
    # Get contiguous column segments for each group
    column_group_data = []
    for column_range in column_ranges:
        column_group_data.append(
            load_data_columns_as_array(filename, data_path, *column_range))
    # Concatenate column groups
    data = np.concatenate(column_group_data, axis=1)

    return data

def load_continuous_as_array(filename, channels):
    """
    Fast way of reading a single channel or a set of channels.

    filename - str - full path to file
    channels - list - channel numbers to include (starting from 0).
               Single channel can be given as a single list element or int.
               Channels in the list must be in sorted (ascending) order.
    """
    # Generate path to raw continuous data
    root_path = '/acquisition/timeseries/' + get_recordingKey(filename) \
                + '/continuous/' + get_processorKey(filename)
    data_path = root_path + '/data'
    timestamps_path = root_path + '/timestamps'
    # Check that data is available, otherwise return None
    if not check_if_path_exists(filename, data_path):
        return None
    if not check_if_path_exists(filename, timestamps_path):
        return None
    # Load continuous data
    continuous = load_data_as_array(filename, data_path, channels)
    # Load timestamps for data
    with h5py.File(filename, 'r') as h5file:
        timestamps = np.array(h5file[timestamps_path])
    # Arrange output into a dictionary
    data = {'continuous': continuous, 'timestamps': timestamps}

    return data

def remove_surrounding_binary_markers(text):
    if text.startswith("b'"):
        text = text[2:]
    if text.endswith("'"):
        text = text[:-1]
    return text

def get_downsampling_info_old(filename):
    # Generate path to downsampling data info
    root_path = '/acquisition/timeseries/' + get_recordingKey(filename) \
                + '/continuous/' + get_processorKey(filename)
    data_path = root_path + '/downsampling_info'
    # Load info from file
    with h5py.File(filename, 'r') as h5file:
        data = h5file[data_path]
        data = [str(i) for i in data]
    # Remove b'x' markers from strings if present. Python 3 change.
    data = list(map(remove_surrounding_binary_markers, data))
    # Parse elements in loaded data
    info_dict = {}
    for x in data:
        key, value = x.split(' ')
        if key == 'original_sampling_rate':
            info_dict[key] = np.int64(value)
        elif key == 'downsampled_sampling_rate':
            info_dict[key] = np.int64(value)
        elif key == 'downsampled_channels':
            info_dict[key] = np.array(list(map(int, value.split(',')))).astype(np.int64)

    return info_dict

def get_downsampling_info(filename):
    root_path = '/acquisition/timeseries/' + get_recordingKey(filename) \
                + '/continuous/' + get_processorKey(filename)
    data_path = root_path + '/downsampling_info/'
    with h5py.File(filename, 'r') as h5file:
        return recursively_load_dict_contents_from_group(h5file, data_path)

def load_downsampled_tetrode_data_as_array(filename, tetrode_nrs):
    """
    Returns a dict with downsampled continuous data for requested tetrodes

    filename                    - str - full path to file
    tetrode_nrs                 - list - tetrode numbers to include (starting from 0).
                                Single tetrode can be given as a single list element or int.
                                Tetrode numbers in the list must be in sorted (ascending) order.
                                If data is not available for a given tetrode number, error is raised.
    """
    # Generate path to raw continuous data
    root_path = '/acquisition/timeseries/' + get_recordingKey(filename) \
                + '/continuous/' + get_processorKey(filename)
    data_path = root_path + '/downsampled_tetrode_data'
    timestamps_path = root_path + '/downsampled_timestamps'
    # Check that data is available, otherwise return None
    if not check_if_path_exists(filename, data_path):
        return None
    if not check_if_path_exists(filename, timestamps_path):
        return None
    # Get info on downsampled data
    info = get_downsampling_info(filename)
    sampling_rate = int(info['downsampled_sampling_rate'])
    downsampled_channels = list(info['downsampled_channels'])
    # Map tetrode_nrs elements to columns in downsampled_tetrode_data
    columns = []
    channels_used = []
    tetrode_nrs_remaining = copy(tetrode_nrs)
    for tetrode_nr in tetrode_nrs:
        for chan in tetrode_channels(tetrode_nr):
            if chan in downsampled_channels:
                columns.append(downsampled_channels.index(chan))
                channels_used.append(chan)
                tetrode_nrs_remaining.pop(tetrode_nrs_remaining.index(tetrode_nr))
                break
    # Check that all tetrode numbers were mapped
    if len(tetrode_nrs_remaining) > 0:
        raise Exception('The following tetrodes were not represented in downsampled data\n' \
                        + ','.join(list(map(str, tetrode_nrs_remaining))))
    # Load continuous data
    continuous = load_data_as_array(filename, data_path, columns)
    # Load timestamps for data
    with h5py.File(filename, 'r') as h5file:
        timestamps = np.array(h5file[timestamps_path])
    # Arrange output into a dictionary
    data = {'continuous': continuous, 'timestamps': timestamps,
            'tetrode_nrs': tetrode_nrs, 'channels': channels_used,
            'sampling_rate': sampling_rate}

    return data

def empty_spike_data():
    """
    Creates a fake waveforms of 0 values and at timepoint 0
    """
    waveforms = np.zeros((1,4,40), dtype=np.int16)
    timestamps = np.array([0], dtype=np.float64)

    return {'waveforms': waveforms, 'timestamps': timestamps}

def get_tetrode_nrs_if_spikes_available(filename, spike_name='spikes'):
    """
    Returns a list of tetrode numbers if spikes available in NWB file.
    """
    spikes_path = '/acquisition/timeseries/' + get_recordingKey(filename) + '/' + spike_name + '/'
    # Get tetrode keys if available
    with h5py.File(filename, 'r') as h5file:
        if not (spikes_path in h5file):
            # Return empty list if spikes data not available
            return []
        tetrode_keys = list(h5file[spikes_path].keys())
    # Return empty list if spikes not available on any tetrode
    if len(tetrode_keys) == 0:
        return []
    # Extract tetrode numbers
    tetrode_nrs = []
    for tetrode_key in tetrode_keys:
        tetrode_nrs.append(int(tetrode_key[9:]) - 1)
    # Sort tetrode numbers in ascending order
    tetrode_nrs.sort()

    return tetrode_nrs

def construct_paths_to_tetrode_spike_data(filename, tetrode_nrs, spike_name='spikes'):
    spikes_path = '/acquisition/timeseries/' + get_recordingKey(filename) + '/' + spike_name + '/'
    return [(spikes_path + 'electrode' + str(tetrode_nr + 1) + '/') for tetrode_nr in tetrode_nrs]

def count_spikes(filename, tetrode_nrs, spike_name='spikes', use_idx_keep=False):
    """
    :param filename: full path to NWB file
    :type filename: str
    :param tetrode_nrs: tetrode numbers to count spikes for
    :type tetrode_nrs: list
    :param spike_name: type of spikes to look for (field in NWB file)
    :type spike_name: str
    :param use_idx_keep: If False (default) all spikes are counted, otherwise only filtered spikes are counted
    :type use_idx_keep: bool
    :return: total number of spikes on each tetrode
    :rtype: list
    """
    tetrode_paths = construct_paths_to_tetrode_spike_data(filename, tetrode_nrs, spike_name=spike_name)
    count = []
    with h5py.File(filename, 'r') as h5file:
        for tetrode_path in tetrode_paths:
            if use_idx_keep:
                count.append(sum(np.array(h5file[tetrode_path + 'idx_keep'][()]).squeeze()))
            else:
                count.append(h5file[tetrode_path + 'timestamps/'].shape[0])

    return count

def load_spikes(filename, spike_name='spikes', tetrode_nrs=None, use_idx_keep=False,
                use_badChan=False, no_waveforms=False, clustering_name=None, verbose=True):
    """
    Inputs:
        filename - pointer to NWB file to load
        tetrode_nrs [list]    - can be a list of tetrodes to load (from 0)
        use_idx_keep [bool]   - if True, only outputs spikes according to idx_keep of tetrode, if available
        use_badChan [bool]    - if True, sets all spikes on badChannels to 0
        no_waveforms [bool]   - if True, waveforms are not loaded
        clustering_name [str] - if specified, clusterID will be loaded from:
                              electrode[nr]/clustering/clustering_name
        verbose [bool]        - prints out loading progress bar if True (default)
    Output:
        List of dictionaries for each tetrode in correct order where:
        List is empty, if no spike data detected
        'waveforms' is a list of tetrode waveforms in the order of channels
        'timestamps' is a list of spike detection timestamps corresponding to 'waveforms'
        If available, two more variables will be in the dictionary
        'idx_keep' is boolan index for 'waveforms' and 'timestamps' indicating the spikes
            that are to be used for further processing (based on filtering for artifacts etc)
        'clusterIDs' is the cluster identities of spikes in 'waveforms'['idx_keep',:,:]
    """
    # If not provided, get tetrode_nrs
    if tetrode_nrs is None:
        tetrode_nrs = get_tetrode_nrs_if_spikes_available(filename, spike_name=spike_name)
    tetrode_paths = construct_paths_to_tetrode_spike_data(filename, tetrode_nrs, spike_name=spike_name)
    with h5py.File(filename, 'r') as h5file:
        # Put waveforms and timestamps into a list of dictionaries in correct order
        data = []
        if verbose:
            print('Loading tetrodes from {}'.format(filename))
        iterable = zip(tetrode_nrs, tetrode_paths)
        for nr_tetrode, tetrode_path in (tqdm(iterable, total=len(tetrode_nrs)) if verbose else iterable):
            # Load waveforms and timestamps
            if no_waveforms:
                waveforms = empty_spike_data()['waveforms']
            else:
                waveforms = h5file[tetrode_path + 'data/'][()]
            timestamps = h5file[tetrode_path + 'timestamps/'][()]
            if not isinstance(timestamps, np.ndarray):
                timestamps = np.array([timestamps])
            if waveforms.shape[0] == 0:
                # If no waveforms are available, enter one waveform of zeros at timepoint zero
                waveforms = empty_spike_data()['waveforms']
                timestamps = empty_spike_data()['timestamps']
            # Arrange waveforms, timestamps and nr_tetrode into a dictionary
            tet_data = {'waveforms': waveforms,
                        'timestamps': timestamps,
                        'nr_tetrode': nr_tetrode}
            # Include idx_keep if available
            idx_keep_path = tetrode_path + 'idx_keep'
            if idx_keep_path in h5file:
                tet_data['idx_keep'] = np.array(h5file[idx_keep_path][()])
                if use_idx_keep:
                    # If requested, filter wavefoms and timestamps based on idx_keep
                    if np.sum(tet_data['idx_keep']) == 0:
                        tet_data['waveforms'] = empty_spike_data()['waveforms']
                        tet_data['timestamps'] = empty_spike_data()['timestamps']
                    else:
                        if not no_waveforms:
                            tet_data['waveforms'] = tet_data['waveforms'][tet_data['idx_keep'], :, :]
                        tet_data['timestamps'] = tet_data['timestamps'][tet_data['idx_keep']]
            # Include clusterIDs if available
            if clustering_name is None:
                clusterIDs_path = tetrode_path + 'clusterIDs'
            else:
                clusterIDs_path = tetrode_path + '/clustering/' + clustering_name
            if clusterIDs_path in h5file:
                tet_data['clusterIDs'] = np.int16(h5file[clusterIDs_path][()]).squeeze()
            # Set spikes to zeros for channels in badChan list if requested
            if use_badChan and not no_waveforms:
                badChan = listBadChannels(filename)
                if len(badChan) > 0:
                    for nchan in tetrode_channels(nr_tetrode):
                        if nchan in badChan:
                            tet_data['waveforms'][:, np.mod(nchan, 4), :] = 0
            data.append(tet_data)

        return data

def save_spikes(filename, tetrode_nr, data, timestamps, spike_name='spikes', overwrite=False):
    """
    Stores spike data in NWB file in the same format as with OpenEphysGUI.
    tetrode_nr=0 for first tetrode.
    """
    if data.dtype != np.int16:
        raise ValueError('Waveforms are not int16.')
    if timestamps.dtype != np.float64:
        raise ValueError('Timestamps are not float64.')
    recordingKey = get_recordingKey(filename)
    path = '/acquisition/timeseries/' + get_recordingKey(filename) + '/' + spike_name + '/' + \
           'electrode' + str(tetrode_nr + 1) + '/'
    if check_if_path_exists(filename, path):
        if overwrite:
            # If overwrite is true, path is first cleared
            with h5py.File(filename, 'r+') as h5file:
                del h5file[path]
        else:
            raise Exception('Spikes already in file and overwrite not requested.\n' \
                            + 'File: ' + filename + '\n' \
                            + 'path: ' + path)
    with h5py.File(filename, 'r+') as h5file:
        h5file[path + 'data'] = data
        h5file[path + 'timestamps'] = np.float64(timestamps).squeeze()

def processing_method_and_spike_name_combinations():
    """
    Outputs a list of potential processing_method and spike_name combinations
    """
    processing_methods = ['klustakwik', 'klustakwik_raw', 'kilosort']
    spike_names = ['spikes', 'spikes_raw', 'spikes_kilosort']

    return processing_methods, spike_names

def get_spike_name_for_processing_method(processing_method):
    processing_methods, spike_names = processing_method_and_spike_name_combinations()
    spike_name = spike_names[processing_methods.index(processing_method)]

    return spike_name

def load_events(filename, internally_generated=False):
    # Outputs a dictionary timestamps and eventIDs for TTL signals received
    # timestamps are in seconds, aligned to timestamps of continuous recording
    # eventIDs indicate TTL channel number (starting from 1) and are positive for rising signals

    if internally_generated:
        ttl_type = 'ttl2'
    else:
        ttl_type = 'ttl1'

    # Load data file
    recordingKey = get_recordingKey(filename)
    with h5py.File(filename, 'r') as h5file:
        # Load timestamps and TLL signal info
        timestamps = h5file['acquisition']['timeseries'][recordingKey]['events'][ttl_type]['timestamps'][()]
        eventID = h5file['acquisition']['timeseries'][recordingKey]['events'][ttl_type]['data'][()]
        data = {'eventID': eventID, 'timestamps': timestamps}

    return data

def load_GlobalClock_timestamps(filename, GlobalClock_TTL_channel=1):
    """
    Returns timestamps of GlobalClock TTL pulses.
    """
    data = load_events(filename)
    return data['timestamps'][data['eventID'] == GlobalClock_TTL_channel]

def load_open_ephys_generated_ttl_events(filename):
    """Returns Open Ephys generated TTL pulse events with channel numbers and timestamps

    :param str filename: full path to NWB file
    :return: channel_event, timestamps
    """
    data = load_events(filename, internally_generated=True)
    return data['eventID'], data['timestamps']

def load_network_events(filename):
    """returns network_events_data

    Extracts the list of network messages from NWB file
    and returns it along with corresponding timestamps
    in dictionary with keys ['messages', 'timestamps']

    'messages' - list of str

    'timestamps' - list of float

    :param filename: full path to NWB file
    :type filename: str
    :return: network_events_data
    :rtype: dict
    """
    # Load data file
    recordingKey = get_recordingKey(filename)
    with h5py.File(filename, 'r') as h5file:
        # Load timestamps and messages
        timestamps = h5file['acquisition']['timeseries'][recordingKey]['events']['text1']['timestamps'][()]
        messages = h5file['acquisition']['timeseries'][recordingKey]['events']['text1']['data'][()]
    messages = [x.decode('utf-8') for x in messages]
    timestamps = [float(x) for x in timestamps]

    data = {'messages': messages, 'timestamps': timestamps}

    return data

def check_if_path_exists(filename, path):
    with h5py.File(filename, 'r') as h5file:
        return path in h5file

def save_list_of_dicts_to_group(h5file, path, dlist, overwrite=False, list_suffix='_NWBLIST'):
    # Check that all elements are dictionaries
    for dic in dlist:
        if not isinstance(dic, dict):
            raise Exception('List elements must be dictionaries')
    # Write elements to file
    for i, dic in enumerate(dlist):
        recursively_save_dict_contents_to_group(h5file, (path + str(i) + '/'), dic,
                                                overwrite=overwrite, list_suffix=list_suffix)

def recursively_save_dict_contents_to_group(h5file, path, dic, overwrite=False, list_suffix='_NWBLIST', verbose=False):
    """
    h5file - h5py.File
    path   - str       - path to group in h5file. Must end with '/'
    overwrite - bool   - any dictionary elements or lists that already exist are overwritten.
                         Default is False, if elements already exist in NWB file, error is raised.
    list_suffix - str  - suffix used to highlight paths created from lists of dictionaries.
                         Must be consistent when saving and loading data.
    verbose - bool     - If True (default is False), h5file path used is printed for each recursion

    Only works with: numpy arrays, numpy int64 or float64, strings, bytes, lists of strings and dictionaries these are contained in.
    Also works with lists dictionaries as part of the hierachy.
    Long lists of dictionaries are discouraged, as individual groups are created for each element.
    """
    if verbose:
        print(path)
    if len(dic) == 0:
        if path in h5file:
            del h5file[path]
        h5file.create_group(path)
    for key, item in dic.items():
        if isinstance(item, (int, float)):
            item = np.array(item)
        if isinstance(item, (np.ndarray, np.int64, np.float64, str, bytes)):
            if overwrite:
                if path + key in h5file:
                    del h5file[path + key]
            h5file[path + key] = item
        elif isinstance(item, dict):
            recursively_save_dict_contents_to_group(h5file, path + key + '/', item,
                                                    overwrite=overwrite, list_suffix=list_suffix,
                                                    verbose=verbose)
        elif isinstance(item, list):
            if all(isinstance(i, str) for i in item):
                if overwrite:
                    if path + key in h5file:
                        del h5file[path + key]
                asciiList = [n.encode("ascii", "ignore") for n in item]
                h5file[path + key] = h5file.create_dataset(None, (len(asciiList),),'S100', asciiList)
            else:
                if overwrite:
                    if path + key + list_suffix in h5file:
                        del h5file[path + key + list_suffix]
                save_list_of_dicts_to_group(h5file, path + key + list_suffix + '/', item,
                                            overwrite=overwrite, list_suffix=list_suffix)
        elif item is None:
            h5file.create_group(path + key)
        else:
            raise ValueError('Cannot save %s type'%type(item) + ' from ' + path + key)

def convert_bytes_to_string(b):
    """
    If input is bytes, returns str decoded with utf-8

    :param b:
    :type b: bytes
    :return: string decoded with utf-8 if input is bytes object, otherwise returns unchanged input
    :rtype: str
    """
    if isinstance(b, bytes):
        if sys.version_info >= (3, 0):
            return str(b, 'utf-8')
        else:
            return str(b.decode('utf-8'))
    else:
        return b

def load_list_of_dicts_from_group(h5file, path, list_suffix='_NWBLIST', ignore=()):
    # Load all elements on this path
    items = []
    for key in list(h5file[path].keys()):
        items.append(
            (int(key), recursively_load_dict_contents_from_group(h5file, path + key + '/',
                                                                 list_suffix=list_suffix,
                                                                 ignore=ignore))
        )
    # Create a list from items sorted by group keys
    ans = [item for _, item in sorted(items)]

    return ans

def recursively_load_dict_contents_from_group(h5file, path, list_suffix='_NWBLIST', ignore=()):
    """
    Returns value at path if it has no further items

    h5file - h5py.File
    path   - str       - path to group in h5file. Must end with '/'
    list_suffix - str  - suffix used to highlight paths created from lists of dictionaries.
                         Must be consistent when saving and loading data.
    ignore - tuple     - paths including elements matching any element in this tuple return None
    """
    if not path.endswith('/'):
        raise ValueError('Input path must end with "/"')

    if path.split('/')[-2] in ignore or path.split('/')[-2][:-len(list_suffix)] in ignore:
        ans = None

    elif path[:-1].endswith(list_suffix):
        ans = load_list_of_dicts_from_group(h5file, path, list_suffix=list_suffix,
                                            ignore=ignore)
    elif hasattr(h5file[path], 'items'):

        ans = {}
        for key, item in h5file[path].items():

            if key.endswith(list_suffix):
                ans[str(key)[:-len(list_suffix)]] = load_list_of_dicts_from_group(
                    h5file, path + key + '/', list_suffix=list_suffix,
                    ignore=ignore
                )

            elif isinstance(item, h5py._hl.dataset.Dataset):
                if 'S100' == item.dtype:
                    tmp = list(item[()])
                    ans[str(key)] = [convert_bytes_to_string(i) for i in tmp]
                elif item.dtype == 'bool' and item.ndim == 0:
                    ans[str(key)] = np.array(bool(item[()]))
                else:
                    ans[str(key)] = convert_bytes_to_string(item[()])

            elif isinstance(item, h5py._hl.group.Group):
                ans[str(key)] = recursively_load_dict_contents_from_group(h5file, path + key + '/',
                                                                          ignore=ignore)

    else:
        ans = convert_bytes_to_string(h5file[path][()])

    return ans

def save_settings(filename, Settings, path='/'):
    """
    Writes into an existing file if path is not yet used.
    Creates a new file if filename does not exist.
    Only works with: numpy arrays, numpy int64 or float64, strings, bytes, lists of strings and dictionaries these are contained in.
    To save specific subsetting, e.g. TaskSettings, use:
        Settings=TaskSetttings, path='/TaskSettings/'
    """
    full_path = '/general/data_collection/Settings' + path
    if os.path.isfile(filename):
        write_method = 'r+'
    else:
        write_method = 'w'
    with h5py.File(filename, write_method) as h5file:
        recursively_save_dict_contents_to_group(h5file, full_path, Settings)

def load_settings(filename, path='/', ignore=()):
    """
    By default loads all settings from path
        '/general/data_collection/Settings/'
    or for example to load animal ID, use:
        path='/General/animal/'

    ignore - tuple - any paths including any element of ignore are returned as None
    """
    full_path = '/general/data_collection/Settings' + path
    with h5py.File(filename, 'r') as h5file:
        data = recursively_load_dict_contents_from_group(h5file, full_path, ignore=ignore)

    return data

def check_if_settings_available(filename, path='/'):
    """
    Returns whether settings information exists in NWB file
    Specify path='/General/badChan/' to check for specific settings
    """
    full_path = '/general/data_collection/Settings' + path
    with h5py.File(filename, 'r') as h5file:
        return full_path in h5file

def save_analysis(filename, data, overwrite=False, complete_overwrite=False, verbose=False):
    """Stores analysis results from nested dictionary to /analysis path in NWB file.

    See :py:func:`NWBio.recursively_save_dict_contents_to_group` for details on supported data structures.

    :param str filename: path to NWB file
    :param dict data: analysis data to be stored in NWB file
    :param bool overwrite: if True, any existing data at same dictionary keys
                           as in previously saved data is overwritten.
                           Default is False.
    :param bool complete_overwrite: if True, all previous analysis data is discarded before writing.
                                    Default is False.
    :param bool verbose: if True (default is False), the path in file for each element is printed.
    """
    with h5py.File(filename, 'r+') as h5file:
        if complete_overwrite:
            del h5file['/analysis']
        recursively_save_dict_contents_to_group(h5file, '/analysis/', data, overwrite=overwrite, verbose=verbose)

def load_analysis(filename, ignore=()):
    """Loads analysis results from /analysis path in NWB file into a dictionary.

    :param str filename: path to NWB file
    :param tuple ignore: paths containing any element of ignore are terminated with None.
        In the output dictionary any elements downstream of a key matching any element of ignore
        is not loaded and dictionary tree is terminated at that point with value None.
    """
    with h5py.File(filename, 'r') as h5file:
        return recursively_load_dict_contents_from_group(h5file, '/analysis/', ignore=ignore)

def listBadChannels(filename):
    if check_if_settings_available(filename,'/General/badChan/'):
        badChanString = load_settings(filename,'/General/badChan/')
        # Separate input string into a list using ',' as deliminaters
        if badChanString.find(',') > -1: # If more than one channel specified
            # Find all values tetrode and channel values listed
            badChanStringList = badChanString.split(',')
        else:
            badChanStringList = [badChanString]
        # Identify any ranges specified with '-' and append these channels to the list
        for chanString in badChanStringList:
            if chanString.find('-') > -1:
                chan_from = chanString[:chanString.find('-')]
                chan_to = chanString[chanString.find('-') + 1:]
                for nchan in range(int(chan_to) - int(chan_from) + 1):
                    badChanStringList.append(str(nchan + int(chan_from)))
                badChanStringList.remove(chanString) # Remove the '-' containing list element
        # Reorder list of bad channels
        badChanStringList.sort(key=int)
        badChan = list(np.array(list(map(int, badChanStringList))) - 1)
    else:
        badChan = []

    return badChan

def save_tracking_data(filename, TrackingData, ProcessedPos=False, overwrite=False):
    """
    TrackingData is expected as dictionary with keys for each source ID
    If saving processed data, TrackingData is expected to be numpy array
        Use ProcessedPos=True to store processed data
        Use overwrite=True to force overwriting existing processed data
    """
    if os.path.isfile(filename):
        write_method = 'r+'
    else:
        write_method = 'w'
    recordingKey = get_recordingKey(filename)
    with h5py.File(filename, write_method) as h5file:
        full_path = '/acquisition/timeseries/' + recordingKey + '/tracking/'
        if not ProcessedPos:
            recursively_save_dict_contents_to_group(h5file, full_path, TrackingData)
        elif ProcessedPos:
            processed_pos_path = full_path + 'ProcessedPos/'
            # If overwrite is true, path is first cleared
            if overwrite:
                if full_path in h5file and 'ProcessedPos' in list(h5file[full_path].keys()):
                    del h5file[processed_pos_path]
            h5file[processed_pos_path] = TrackingData

def get_recording_cameraIDs(filename):
    path = '/general/data_collection/Settings/CameraSettings/CameraSpecific'
    with h5py.File(filename, 'r') as h5file:
        if path in h5file:
            return list(h5file[path].keys())

def load_raw_tracking_data(filename, cameraID, specific_path=None):
    path = '/acquisition/timeseries/' + get_recordingKey(filename) + '/tracking/' + cameraID + '/'
    if not (specific_path is None):
        path = path + '/' + specific_path + '/'
    with h5py.File(filename, 'r') as h5file:
        if path in h5file:
            return recursively_load_dict_contents_from_group(h5file, path)

def load_processed_tracking_data(filename, subset='ProcessedPos'):
    path = '/acquisition/timeseries/' + get_recordingKey(filename) + '/tracking/'
    path = path + subset
    with h5py.File(filename, 'r') as h5file:
        return np.array(h5file[path][()])

def get_processed_tracking_data_timestamp_edges(filename, subset='ProcessedPos'):
    if check_if_processed_position_data_available(filename):
        data = load_processed_tracking_data(filename, subset=subset)
        edges = [data[0, 0], data[-1, 0]]
    else:
        print('Warning! ProcessedPos not available. Using continuous data timestamps')
        h5file = h5py.File(filename, 'r')
        recordingKey = get_recordingKey(filename)
        processorKey = get_processorKey(filename)
        path = '/acquisition/timeseries/' + recordingKey + '/continuous/' + processorKey + '/timestamps'
        edges = [h5file[path][0], h5file[path][-1]]
        h5file.close()

    return edges

def check_if_tracking_data_available(filename):
    if check_if_settings_available(filename, path='/General/Tracking/'):
        return load_settings(filename, path='/General/Tracking/')
    else:
        return False

def check_if_processed_position_data_available(filename, subset='ProcessedPos'):
    path = '/acquisition/timeseries/' + get_recordingKey(filename) + '/tracking/'
    path = path + subset
    return check_if_path_exists(filename, path)

def check_if_binary_pos(filename):
    # Checks if binary position data exists in NWB file
    path = '/acquisition/timeseries/' + get_recordingKey(filename) + '/events/binary1/'
    return check_if_path_exists(filename, path)

def save_tetrode_idx_keep(filename, ntet, idx_keep, spike_name='spikes', overwrite=False):
    path = '/acquisition/timeseries/' + get_recordingKey(filename) + '/' + spike_name + '/' + \
           'electrode' + str(ntet + 1) + '/idx_keep'
    with h5py.File(filename, 'r+') as h5file:
        if path in h5file:
            if overwrite:
                del h5file[path]
            else:
                raise ValueError('Tetrode ' + str(ntet + 1) + ' idx_keep already exists in ' + filename)
        h5file[path] = idx_keep

def save_tetrode_clusterIDs(filename, ntet, clusterIDs, spike_name='spikes', overwrite=False):
    path = '/acquisition/timeseries/' + get_recordingKey(filename) + '/' + spike_name + '/' + \
           'electrode' + str(ntet + 1) + '/clusterIDs'
    with h5py.File(filename, 'r+') as h5file:
        if path in h5file:
            if overwrite:
                del h5file[path]
            else:
                raise ValueError('Tetrode ' + str(ntet + 1) + ' clusterIDs already exists in ' + filename)
        h5file[path] = np.int16(clusterIDs).squeeze()

def fill_empty_dictionary_from_source(selection, src_dict):
    """
    Populates a dictionary with None values with values from a source
    dictionary with identical structure.
    """
    dst_dict = copy(selection)
    for key, item in dst_dict.items():
        if isinstance(item, dict):
            dst_dict[key] = fill_empty_dictionary_from_source(item, src_dict[key])
        elif item is None:
            dst_dict[key] = src_dict[key]
        else:
            raise ValueError('Destination dictionary has incorrect.')

    return dst_dict

def get_recording_start_timestamp_offset(filename):
    """Returns the first timestamp of raw or downsampled continuous data.

    :param str filename: path to NWB file
    :return: first timestamp of continuous data
    :rtype: float
    """
    if check_if_raw_data_available(filename):
        path = get_raw_data_paths(filename)['timestamps']
    elif check_if_downsampled_data_available(filename):
        path = get_downsampled_data_paths(filename)['timestamps']
    else:
        raise Exception('NWB file does not contain raw or downsampled data ' + filename)
    with h5py.File(filename, 'r') as h5file:
        return float(h5file[path][0:1])

def get_recording_full_duration(filename):
    """Returns the total duration from first to last timestamp of
    raw or downsampled continuous data.

    :param str filename: path to NWB file
    :return: total duration from first to last timestamp of continuous data
    :rtype: float
    """
    if check_if_raw_data_available(filename):
        path = get_raw_data_paths(filename)['timestamps']
    elif check_if_downsampled_data_available(filename):
        path = get_downsampled_data_paths(filename)['timestamps']
    else:
        raise Exception('NWB file does not contain raw or downsampled data ' + filename)
    with h5py.File(filename, 'r') as h5file:
        return float(h5file[path][-1]) - float(h5file[path][0:1])

def import_task_specific_log_parser(task_name):
    """
    Returns LogParser module for the specific task.

    :param task_name: name of the task
    :type task_name: str
    :return: TaskLogParser
    :rtype: module
    """
    if task_name == 'Pellets_and_Rep_Milk_Task':  # Temporary workaround to function with older files
        task_name = 'Pellets_and_Rep_Milk'
    try:
        return importlib.import_module('.Tasks.' + task_name + '.LogParser', package='openEPhys_DACQ')
    except ModuleNotFoundError:
        print('Task {} LogParser not found. Returning None.'.format(task_name))
        return None

def load_task_name(filename):
    """
    Returns the name of the task active in the recording.

    :param filename: absolute path to NWB recording file
    :type filename: str
    :return: task_name
    :rtype: str
    """
    return load_settings(filename, path='/TaskSettings/name/')

def get_recording_log_parser(filename, final_timestamp=None):
    """Finds task specific LogParser class and returns it initialized
    with network events from that recording.

    :param str filename:
    :return: Task specific log parser initialized with network events
    :rtype: LogParser class
    """
    task_log_parser = import_task_specific_log_parser(load_task_name(filename))
    if task_log_parser is None:
        return None
    else:
        return task_log_parser.LogParser(task_settings=load_settings(filename, path='/TaskSettings/'),
                                         final_timestamp=final_timestamp,
                                         **load_network_events(filename))

def get_channel_map(filename):
    return load_settings(filename, '/General/channel_map/')

def list_tetrode_nrs_for_area_channel_map(area_channel_map):
    return list(set([channels_tetrode(chan) for chan in list(area_channel_map['list'])]))

def get_channel_map_with_tetrode_nrs(filename):
    channel_map = get_channel_map(filename)
    for area in channel_map:
        channel_map[area]['tetrode_nrs'] = list_tetrode_nrs_for_area_channel_map(channel_map[area])

    return channel_map

def check_if_channel_maps_are_same(channel_map_1, channel_map_2):
    """
    Determines if two channel maps are identical
    """
    # Check that there are same number of areas in the dictionary
    if len(channel_map_1) != len(channel_map_2):
        return False
    # Sort the area names because dictionary is not ordered
    channel_map_1_keys = sorted(list(channel_map_1.keys()))
    channel_map_2_keys = sorted(list(channel_map_2.keys()))
    # Check that the areas have the same name
    for n_area in range(len(channel_map_1_keys)):
        if channel_map_1_keys[n_area] != channel_map_2_keys[n_area]:
            return False
    # Check that the channel lists are the same
    for area in channel_map_1_keys:
        if not all(channel_map_1[area]['list'] == channel_map_2[area]['list']):
            return False

    return True

def estimate_open_ephys_timestamps_from_other_timestamps(open_ephys_global_clock_times, other_global_clock_times,
                                                         other_times, other_times_divider=None):
    """Returns Open Ephys timestamps for each timestamp from another device by synchronising with global clock.

    Note, other times must be in same units as open_ephys_global_clock_times. Most likely seconds.
    For example, Raspberry Pi camera timestamps would need to be divided by 10 ** 6

    :param numpy.ndarray open_ephys_global_clock_times: shape (N,)
    :param numpy.ndarray other_global_clock_times: shape (M,)
    :param numpy.ndarray other_times: shape (K,)
    :param int other_times_divider: if provided, timestamps from the other devices are divided by this value
        before matching to Open Ephys time. This allows inputting timestamps from other device in original units.
        In case of Raspberry Pi camera timestamps, this value should be 10 ** 6.
        If this value is not provided, all provided timestamps must be in same units.
    :return: open_ephys_times
    :rtype: numpy.ndarray
    """

    # Crop data if more timestamps recorded on either system.
    if open_ephys_global_clock_times.size > other_global_clock_times.size:
        open_ephys_global_clock_times = open_ephys_global_clock_times[:other_global_clock_times.size]
        print('[ Warning ] OpenEphys recorded more GlobalClock TTL pulses than other system.\n' +
              'Dumping extra OpenEphys timestamps from the end.')
    elif open_ephys_global_clock_times.size < other_global_clock_times.size:
        other_global_clock_times = other_global_clock_times[:open_ephys_global_clock_times.size]
        print('[ Warning ] Other system recorded more GlobalClock TTL pulses than Open Ephys.\n' +
              'Dumping extra other system timestamps from the end.')

    # Find closest other_global_clock_times indices to each other_times
    other_times_gc_indices = closest_argmin(other_times, other_global_clock_times)

    # Compute difference from the other_global_clock_times for each value in other_times
    other_times_nearest_global_clock_times = other_global_clock_times[other_times_gc_indices]
    other_times_global_clock_delta = other_times - other_times_nearest_global_clock_times

    # Convert difference values to Open Ephys timestamp units
    if not (other_times_divider is None):
        other_times_global_clock_delta = other_times_global_clock_delta / float(other_times_divider)

    # Use other_times_global_clock_delta to estimate timestamps in OpenEphys time
    other_times_nearest_open_ephys_global_clock_times = open_ephys_global_clock_times[other_times_gc_indices]
    open_ephys_times = other_times_nearest_open_ephys_global_clock_times + other_times_global_clock_delta

    return open_ephys_times

def extract_recording_info(filename, selection='default'):
    """
    Returns recording info for the recording file.

    selection - allows specifying which data return
        'default' - some hard-coded selection of data
        'all' - all of the recording settings
        dict - a dictionary with the same exact keys and structure
               as the recording settings, with None for item values
               and missing keys for unwanted elements. The dictionary
               will be returned with None values populated by values
               from recording settings.
    """
    recording_info = {}
    if isinstance(selection, str) and selection == 'default':
        recording_info.update(load_settings(filename, '/General/'))
        del recording_info['experimenter']
        del recording_info['rec_file_path']
        del recording_info['root_folder']
        if recording_info['TaskActive']:
            recording_info.update({'TaskName': load_settings(filename, '/TaskSettings/name/')})
        for key in list(recording_info['channel_map'].keys()):
            del recording_info['channel_map'][key]['list']
        pos_edges = get_processed_tracking_data_timestamp_edges(filename)
        recording_info['duration'] = pos_edges[1] - pos_edges[0]
        recording_info['duration (min)'] = int(round((pos_edges[1] - pos_edges[0]) / 60))
        recording_info['time'] = load_settings(filename, '/Time/')
    elif isinstance(selection, str) and selection == 'all':
        recording_info = load_settings(filename)
    elif isinstance(selection, dict):
        full_recording_info = load_settings(filename)
        recording_info = fill_empty_dictionary_from_source(selection, full_recording_info)

    return recording_info

def display_recording_data(root_path, selection='default'):
    """
    Prints recording info for the whole directory tree.
    """
    for dirName, subdirList, fileList in os.walk(root_path):
        for fname in fileList:
            if fname == 'experiment_1.nwb':
                filename = os.path.join(dirName, fname)
                recording_info = extract_recording_info(filename, selection)
                print('Data on path: ' + dirName)
                pprint(recording_info)

if __name__ == '__main__':
    # Input argument handling and help info
    parser = argparse.ArgumentParser(description='Extract info from Open Ephys.')
    parser.add_argument('root_path', type=str, nargs=1,
                        help='Root directory for recording(s)')
    args = parser.parse_args()
    # Get paths to recording files
    display_recording_data(args.root_path[0])

max_stars_repo_name: Barry-lab/SpatialAutoDACQ

max_issues_repo_name: Barry-lab/SpatialAutoDACQ

max_forks_repo_name: Barry-lab/SpatialAutoDACQ


text:
#!/usr/bin/python
import numpy as np
import matplotlib
import time

def matplotlib_numpy_test():
    # evenly sampled time at 200ms intervals
    t = np.arange(0., 5., 0.2)
    time.sleep(0.5)


text:
from theano import *
import numpy
import theano.tensor as T

x = T.dmatrix('x')
y = T.dmatrix('y')
z = x + y
f = function([x,y],z)

print(f(\
  [[2,3],[4,5]],\
  [[2,3],[4,5]]\
))

max_stars_repo_name: ibiscum/Learning-Concurrency-in-Python

reason: import numpy,import theano,from theano