text:
#-------------------------------------------------------------------------------
# Name:        Image processing with TensorFlow example
# Purpose:     Experiments with processing images in TF 2.0
#
# Author:      kol
#
# Created:     10.01.2020
# Copyright:   (c) kol 2020
# Licence:     MIT
#-------------------------------------------------------------------------------

import os
import tensorflow as tf
import matplotlib.pyplot as plt
from pathlib import Path
import numpy as np
from PIL import Image
from functools import partial
from random import randrange
from time import time

AUTOTUNE = tf.data.experimental.AUTOTUNE
IMG_DIR = str(Path(__file__).absolute().parent.parent.joinpath('cc', 'p'))
MODEL_DIR =  str(Path(__file__).absolute().parent.joinpath('gbr'))
LOGS_DIR = str(Path(__file__).absolute().parent.joinpath('_logs'))
NUM_EPOCHS = 15
BATCH_SIZE = 32
IMG_HEIGHT = 20
IMG_WIDTH = 20
DISPLAY_COLS = 6

def map_fn(path, label):
    image = tf.image.decode_png(tf.io.read_file(path))

    image = tf.image.convert_image_dtype(image, tf.float32)
    image = tf.image.resize(image, [IMG_WIDTH, IMG_HEIGHT])

    return image, label

def plot_image(img, predicted_label, true_label, class_names):
    plt.grid(False)
    plt.xticks([])
    plt.yticks([])

    plt.imshow(img, cmap=plt.cm.binary)

    plt.title('{} ({})'.format(
        class_names[int(np.round(predicted_label, 0))],
        class_names[int(true_label)]))

def show_batch(images, labels, class_names):
    plt.figure(figsize=(5,5))
    for n in range(min(25, images.shape[0])):
        ax = plt.subplot(5,5,n+1)
        plt.imshow(images[n])
        if len(labels.shape) == 1:
            plt.title(class_names[int(labels[n])].title())
        else:
            m = np.argmax(labels[n])
            plt.title(class_names[int(labels[n, m])].title())
        plt.axis('off')

    plt.tight_layout()
    plt.show()

def plot_history(history):
    acc = history.history['accuracy']
    val_acc = history.history['val_accuracy']

    loss = history.history['loss']
    val_loss = history.history['val_loss']

    epochs = len(acc)
    epochs_range = range(epochs)

    plt.figure(figsize=(8, 8))
    plt.subplot(1, 2, 1)
    plt.plot(epochs_range, acc, label='Training Accuracy')
    plt.plot(epochs_range, val_acc, label='Validation Accuracy')
    plt.legend(loc='lower right')
    plt.title('Training and Validation Accuracy')

    plt.subplot(1, 2, 2)
    plt.plot(epochs_range, loss, label='Training Loss')
    plt.plot(epochs_range, val_loss, label='Validation Loss')
    plt.legend(loc='upper right')
    plt.title('Training and Validation Loss')
    plt.show()

def plot_predictions(predictions, images, labels, class_names):
    fig = plt.figure(figsize=(8, 8))
    num_rows = int(np.ceil(BATCH_SIZE / DISPLAY_COLS))

    for i in range(BATCH_SIZE-1):
        n = randrange(0, predictions.shape[0]-1)
        fig.add_subplot(num_rows, DISPLAY_COLS, i+1)
        plot_image(images[n], predictions[n], labels[n], class_names)

    plt.tight_layout()
    plt.show()

# Load images
print("==> Loading images")
class_names = np.array([item.name for item in Path(IMG_DIR).glob('*') if item.is_dir()])
print("Classes found: ", class_names)

class_names = []
file_names = []
file_labels = []
for n, d in enumerate(Path(IMG_DIR).glob('*')):
    class_names.extend([d.name])
    names = [str(f) for f in Path(IMG_DIR).joinpath(d).glob('*.png')]
    file_names.extend(names)
    file_labels.extend([1.0 * n] * len(names))

image_count = len(file_names)
print("Images found: ", image_count)

file_names = tf.convert_to_tensor(file_names, dtype=tf.string)
file_labels = tf.convert_to_tensor(file_labels)
file_labels = tf.expand_dims(file_labels, axis=-1)

# Check model exists and load/train it
if Path(MODEL_DIR).exists():
    print("==> Loading pre-trained model")
    model = tf.keras.models.load_model(MODEL_DIR)
else:
    # Generate dataset
    print("==> Generate datasets")
    dataset = tf.data.Dataset.from_tensor_slices((file_names, file_labels))
    dataset = dataset.repeat().shuffle(image_count, reshuffle_each_iteration=False)

    dataset = dataset.map(map_fn, num_parallel_calls=AUTOTUNE)
    dataset = dataset.batch(BATCH_SIZE)
    dataset = dataset.prefetch(buffer_size=AUTOTUNE)

    images, labels = next(iter(dataset))
    show_batch(images.numpy(), labels.numpy(), class_names)

    train_size = int(0.7 * image_count)
    val_size = int(0.15 * image_count)
    test_size = int(0.15 * image_count)

    train_dataset = dataset.take(train_size)
    test_dataset = dataset.skip(train_size)
    val_dataset = test_dataset.skip(test_size)
    test_dataset = test_dataset.take(test_size)

    print("==> Building model")
    model = tf.keras.models.Sequential([
        tf.keras.layers.Conv2D(16, 3, padding='same', activation='relu', input_shape=(IMG_HEIGHT, IMG_WIDTH, 3)),
        tf.keras.layers.MaxPooling2D((2, 2)),
        tf.keras.layers.Conv2D(32, 3, padding='same', activation='relu'),
        tf.keras.layers.MaxPooling2D((2, 2)),
        tf.keras.layers.Conv2D(64, 3, padding='same', activation='relu'),
        tf.keras.layers.MaxPooling2D((2, 2)),
        tf.keras.layers.Flatten(),
        tf.keras.layers.Dense(512, activation='relu'),
        tf.keras.layers.Dense(1, activation='sigmoid')
    ])
    model.summary()

    model.compile(optimizer='adam',
                  loss='binary_crossentropy',
                  metrics=['accuracy'])

    print("==> Training model")
    history = model.fit(train_dataset,
              epochs=NUM_EPOCHS,
              steps_per_epoch=train_size // BATCH_SIZE,
              validation_data=val_dataset,
              validation_steps=val_size // BATCH_SIZE,
              callbacks = [
                #tf.keras.callbacks.EarlyStopping(monitor='val_loss', min_delta=1e-5, patience=2),
                tf.keras.callbacks.TensorBoard(log_dir=LOGS_DIR,
                                               histogram_freq=1,
                                               write_graph=True,
                                               write_grads=True,
                                               write_images=True,
                                               update_freq="batch") ],
              use_multiprocessing=True)

    print("==> Evaluating results")
    model.evaluate(test_dataset, verbose=2)

    print("==> Saving model")
    model.save(MODEL_DIR)

    print("==> Plotting history")
    plot_history(history)

print("==> Generating prediction dataset")
predict_dataset = tf.data.Dataset.from_tensor_slices((file_names, file_labels))
predict_dataset = predict_dataset.shuffle(image_count, reshuffle_each_iteration=False)
predict_dataset = predict_dataset.map(map_fn, num_parallel_calls=AUTOTUNE)
predict_dataset = predict_dataset.batch(BATCH_SIZE)
predict_dataset = predict_dataset.take(1)

print("==> Predicting on {} images".format(BATCH_SIZE))
predictions = model.predict(predict_dataset)

print("==> Displaying results")
images, labels = next(iter(predict_dataset))
plot_predictions(predictions, images, labels, class_names)

max_stars_repo_name: skolchin/gbr

max_issues_repo_name: skolchin/gbr

max_forks_repo_name: skolchin/gbr


text:
from collections import namedtuple
from itertools import izip

import numpy as np
import torch
from torch.autograd import Variable

from gtd.ml.torch.utils import GPUVariable, conditional, is_binary, assert_tensor_equal
from gtd.ml.torch.utils import expand_dims_for_broadcast, NamedTupleLike
from gtd.ml.vocab import Vocab

class SequenceBatch(namedtuple('SequenceBatch', ['values', 'mask']), NamedTupleLike):
    """
    Attributes:
        values (Variable): of shape (batch_size, max_seq_length, X1, X2, ...)
        mask (Variable[FloatTensor]): of shape (batch_size, max_seq_length)
    """
    __slots__ = ()
    def __new__(cls, values, mask, left_justify=True):
        if not isinstance(values, Variable) or not isinstance(mask, Variable):
            raise ValueError('values and mask must both be of type Variable.')

        m = mask.data

        if len(m.size()) == 0:
            raise ValueError('Mask must not be 0-dimensional')

        # check that mask is binary
        if not is_binary(m):
            raise ValueError('Mask must be binary:\n{}'.format(mask))

        # check that mask is left-justified
        # since mask is binary, we just need to check that it is monotonically non-increasing from left to right
        batch_size, seq_len = m.size()
        if seq_len > 1 and left_justify:
            diffs = m[:, 1:] - m[:, :-1]  # (batch_size, max_seq_length - 1)
            non_increasing = diffs <= 0
            all_non_increasing = (torch.prod(non_increasing) == 1)
            if not all_non_increasing:
                raise ValueError('Mask must be left-justified:\n{}'.format(mask))

        self = super(SequenceBatch, cls).__new__(cls, values, mask)
        return self

    @classmethod
    def from_sequences(cls, sequences, vocab_or_vocabs, min_seq_length=0, volatile=False):
        """Convert a batch of sequences into a SequenceBatch.

        Args:
            sequences (list[list[unicode]])
            vocab_or_vocabs (WordVocab|list[WordVocab]): either a single vocab, or a list of vocabs, one per sequence
            min_seq_length (int): enforce that the Tensor representing the SequenceBatch have at least
                this many columns.
            volatile (bool): whether to make Variables volatile (don't track grads)

        Returns:
            SequenceBatch
        """
        # determine dimensions
        batch_size = len(sequences)
        if batch_size == 0:
            seq_length = 0
        else:
            seq_length = max(len(seq) for seq in sequences)  # max seq length in batch
        seq_length = max(seq_length, min_seq_length)  # make sure it is at least min_seq_length
        shape = (batch_size, seq_length)

        # set up vocabs
        if isinstance(vocab_or_vocabs, list):
            vocabs = vocab_or_vocabs
            assert len(vocabs) == batch_size
        else:
            # duplicate a single vocab
            assert isinstance(vocab_or_vocabs, Vocab)
            vocabs = [vocab_or_vocabs] * batch_size

        # build arrays
        values = np.zeros(shape, dtype=np.int64)  # pad with zeros
        mask = np.zeros(shape, dtype=np.float32)
        for i, (seq, vocab) in enumerate(izip(sequences, vocabs)):
            for j, word in enumerate(seq):
                values[i, j] = vocab.word2index(word)
                mask[i, j] = 1.0

        return SequenceBatch(GPUVariable(torch.from_numpy(values), volatile=volatile),
                             GPUVariable(torch.from_numpy(mask), volatile=volatile))

    def split(self):
        """Convert SequenceBatch into a list of Variables, where each element represents one time step.

        Returns:
            list[SequenceBatchElement]: a list of SequenceBatchElements, where for each list element:
                element.values has shape (batch_size, X1, X2, ...)
                element.mask has shape (batch_size, 1)
        """
        values_list = [v.squeeze(dim=1) for v in self.values.split(1, dim=1)]
        mask_list = self.mask.split(1, dim=1)
        return [SequenceBatchElement(v, m) for v, m in izip(values_list, mask_list)]

    @classmethod
    def cat(cls, elements):
        """Concatenate SequenceBatchElements to form a SequenceBatch.

        Args:
            elements (list[SequenceBatchElement])

        Returns:
            SequenceBatch
        """
        values = torch.cat([e.values.unsqueeze(1) for e in elements], 1)
        mask = torch.cat([e.mask for e in elements], 1)
        return SequenceBatch(values, mask)

    @classmethod
    def weighted_sum(cls, seq_batch, weights):
        """Compute weighted sum of elements in a SequenceBatch.

        Args:
            seq_batch (SequenceBatch): with values of shape (batch_size, seq_length, X1, X2, ...)
            weights (Variable): of shape (batch_size, seq_length)

        Returns:
            Variable: of shape (batch_size, X1, X2, ...)
        """
        values = seq_batch.values
        mask = seq_batch.mask
        weights = weights * mask  # ignore weights outside mask
        weights = expand_dims_for_broadcast(weights, values).expand(values.size())
        weighted = values * weights
        return torch.sum(weighted, dim=1).squeeze(dim=1)

    @classmethod
    def reduce_sum(cls, seq_batch):
        weights = GPUVariable(torch.ones(*seq_batch.mask.size()))
        return cls.weighted_sum(seq_batch, weights)

    @classmethod
    def reduce_prod(cls, seq_batch):
        """Compute the product of each sequence in a SequenceBatch.

        If a sequence is empty, we return a product of 1.

        Args:
            seq_batch (SequenceBatch): of shape (batch_size, seq_length, X1, X2, ...)

        Returns:
            Tensor: of shape (batch_size, X1, X2, ...)
        """
        mask = seq_batch.mask
        values = seq_batch.values

        # We set all pad values = 1, so that taking the log will not produce -inf
        mask_bcast = expand_dims_for_broadcast(mask, values).expand(values.size())  # (batch_size, seq_length, X1, X2, ...)
        values = conditional(mask_bcast, values, 1 - mask_bcast)

        logged = SequenceBatch(torch.log(values), seq_batch.mask)  # (batch_size, seq_length, X1, X2, ...)

        log_sum = SequenceBatch.reduce_sum(logged)  # (batch_size, X1, X2, ...)
        prod = torch.exp(log_sum)
        return prod

    @classmethod
    def reduce_mean(cls, seq_batch, allow_empty=False):
        """Compute the mean of each sequence in a SequenceBatch.

        Args:
            seq_batch (SequenceBatch): a SequenceBatch with the following attributes:
                values (Tensor): a Tensor of shape (batch_size, seq_length, X1, X2, ...)
                mask (Tensor): if the mask values are arbitrary floats (rather than binary), the mean will be
                a weighted average.
            allow_empty (bool): allow computing the average of an empty sequence. In this case, we assume 0/0 == 0, rather
                than NaN. Default is False, causing an error to be thrown.

        Returns:
            Tensor: of shape (batch_size, X1, X2, ...)
        """
        values, mask = seq_batch.values, seq_batch.mask
        # compute weights for the average
        sums = torch.sum(mask, dim=1)  # (batch_size, 1)

        if allow_empty:
            sums[sums == 0.0] = 1.0  # Modify in-place: replace zeros with ones
        else:
            if (sums.data == 0).any():
                raise ValueError("Averaging zero elements.")

        weights = mask / sums.expand(*mask.size())
        return cls.weighted_sum(seq_batch, weights)

    @classmethod
    def _empty_seqs(cls, seq_batch):
        return (torch.sum(seq_batch.mask, 1).data == 0).any()

    @classmethod
    def reduce_max(cls, seq_batch):
        if cls._empty_seqs(seq_batch):
            raise ValueError("Taking max over zero elements.")
        values, mask = seq_batch.values, seq_batch.mask

        inf_mask = mask.clone()  # (batch_size, seq_length)
        inf_mask[mask == 0] = float('inf')
        inf_mask[mask == 1] = 0
        # masked elements will never win the max, because we subtract infinity from them

        inf_mask_bcast = expand_dims_for_broadcast(inf_mask, values).expand_as(values)  # (batch_size, seq_length, X1, X2, ...)

        max_values, _ = torch.max(values - inf_mask_bcast, 1)  # (batch_size, 1, X1, X2, ...)
        max_values = torch.squeeze(max_values, 1)  # (batch_size, X1, X2, ...)

        return max_values

    @classmethod
    def log_sum_exp(cls, seq_batch):
        """Numerically stable computation of log-sum-exp.

        Mask must be left-justified.
        Does not allow empty rows.

        x = seq_batch.values
        lse[i] = log(exp(x[i, 0]) + exp(x[i, 1]) + ... + exp(x[i, n]))

        Uses the log-sum-exp stability trick:
        https://hips.seas.harvard.edu/blog/2013/01/09/computing-log-sum-exp/

        Args:
            seq_batch (SequenceBatch): of shape (batch_size, seq_length), where seq_batch.values of type FloatTensor

        Returns:
            lse (Variable): of shape (batch_size,)
        """
        values, mask = seq_batch.values, seq_batch.mask
        max_vals = cls.reduce_max(seq_batch)  # (batch_size,)
        shifted_values = values - max_vals.unsqueeze(1).expand(values.size())  # (batch_size, seq_length)
        exponentiated_values = torch.exp(shifted_values)
        sums = cls.reduce_sum(SequenceBatch(exponentiated_values, mask))  # (batch_size)
        log_sums = torch.log(sums)
        return log_sums + max_vals

    @classmethod
    def embed(cls, indices, embeds):
        """Embed a SequenceBatch of integers.

        Args:
            indices (SequenceBatch): of shape (batch_size, seq_length), with seq_batch.values of type LongTensor (ints)
            embeds (Variable): of shape (vocab_size, embed_dim)

        Returns:
            SequenceBatch: of shape (batch_size, seq_length, embed_dim)
        """
        values, mask = indices

        batch_size, seq_length = values.size()
        vocab_size, embed_dim = embeds.size()

        indices_flat = values.view(batch_size * seq_length)
        embedded_indices_flat = torch.index_select(embeds, 0, indices_flat)  # (batch_size * seq_length, embed_dim)
        embedded_indices = embedded_indices_flat.view(batch_size, seq_length, embed_dim)
        return SequenceBatch(embedded_indices, mask)

    @classmethod
    def multi_vocab_indices(self, sequences, vocabs, min_seq_length=0):
        """Convert a batch of sequences into indices, where each token is converted into a **tuple** of indices.

        Args:
            sequences (list[list[unicode]]): a batch of sequences
            vocabs (list[list[Vocab]]): vocabs[v] = a batch of vocabs (one per example in the batch), corresponding to vocab v
            min_seq_length (int): see SequenceBatch.from_sequences

        Returns:
            SequenceBatch:
                mask (Variable): has shape (batch_size, seq_length)
                values (Variable): has shape (batch_size, seq_length, num_vocabs)
        """
        seq_batches = [SequenceBatch.from_sequences(sequences, vocabs_v, min_seq_length) for vocabs_v in vocabs]
        values_list, mask_list = zip(*seq_batches)

        for mask in mask_list:
            assert_tensor_equal(mask, mask_list[0])  # all masks should be the same

        values = torch.stack(values_list, 2)

        return SequenceBatch(values, mask_list[0])

SequenceBatchElement = namedtuple('SequenceBatchElement', ['values', 'mask'])

ext: py

reason: import numpy


text:
import numpy as np
import pickle
import cv2

num_classes = 5

trainId_to_count = {}
for trainId in range(num_classes):
    trainId_to_count[trainId] = 0

train_label_img_paths = []
trains = open("../../data/signate/train_basefiles.txt").read().split("\n")[:-1]
train_label_img_paths += trains
vals = open("../../data/signate/val_basefiles.txt").read().split("\n")[:-1]
train_label_img_paths += vals
print(len(train_label_img_paths))

# get the total number of pixels in all train label_imgs that are of each object class:
for step, label_img_path in enumerate(train_label_img_paths):
    print(step)

    path = "../../data/signate/seg_train_dat/" + label_img_path + ".dat"
    fp = open(path,'rb')
    dat = np.fromfile(fp, np.uint8, -1)
    fp.close()
    label_img = dat.reshape((1216, 1936, 1))

    for trainId in range(num_classes):
        # count how many pixels in label_img which are of object class trainId:
        trainId_mask = np.equal(label_img, trainId)
        trainId_count = np.sum(trainId_mask)

        # add to the total count:
        trainId_to_count[trainId] += trainId_count

# compute the class weights according to the ENet paper:
class_weights = []
total_count = sum(trainId_to_count.values())
for trainId, count in trainId_to_count.items():
    trainId_prob = float(count)/float(total_count)
    trainId_weight = 1/np.log(1.02 + trainId_prob)
    class_weights.append(trainId_weight)

print (class_weights)

with open("./signate_class_weights.pkl", "wb") as file:
    pickle.dump(class_weights, file, protocol=2) # (protocol=2 is needed to be able to open this file with python2)


text:
import gym
import gym_line_follower

import numpy
import time

import RLAgents

import models.ppo_baseline_discrete.src.model            as Model
import models.ppo_baseline_discrete.src.config           as Config

class Wrapper(gym.Wrapper):
    def __init__(self, env):
        gym.Wrapper.__init__(self, env)

        self.action_space       = gym.spaces.Discrete(4)

        actions = []
        actions.append([0.0, 0.0])
        actions.append([0.0, 1.0])
        actions.append([1.0, 0.0])
        actions.append([1.0, 1.0])

        self.actions = numpy.array(actions, dtype=numpy.float32)

    def step(self, action):
        act = self.actions[action]
        return self.env.step(act)

    def reset(self):
        return self.env.reset()

path = "models/ppo_baseline_discrete/"

config  = Config.Config()

envs = RLAgents.MultiEnvSeq('LineFollower-v0', Wrapper, config.actors)

agent = RLAgents.AgentPPO(envs, Model, Config)

max_iterations = 1*(10**6)

trainig = RLAgents.TrainingIterations(envs, agent, max_iterations, path, 1000)
trainig.run()

'''
agent.load(path)
agent.disable_training()

while True:
    agent.main()
    envs.render(0)
'''


text:
#!/usr/bin/env python
# coding: utf-8
import os
from datetime import datetime, timedelta
import sys

import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
# from arviz import hpd

from rakott.mpl import savefig_bbox

sns.set_context('paper', font_scale=1.3)
red, blue, green, purple = sns.color_palette('Set1', 4)
colors = {'red':red, 'blue':blue, 'green':green, 'purple':purple}

def int_to_dt(t):
    return pd.to_datetime(start_date) + timedelta(days=t)

def calc_hpd(trace, mass_frac) :
    # Get sorted list
    d = np.sort(np.copy(trace))

    # Number of total samples taken
    n = len(trace)

    # Get number of samples that should be included in HPD
    n_samples = np.floor(mass_frac * n).astype(int)

    # Get width (in units of data) of all intervals with n_samples samples
    int_width = d[n_samples:] - d[:n-n_samples]

    # Pick out minimal interval
    min_int = np.argmin(int_width)

    # Return interval
    return np.array([d[min_int], d[min_int+n_samples]])

if __name__ == '__main__':
	job_id = sys.argv[1]
	output_folder = r'../../output-tmp/{}/'.format(job_id)
	country = sys.argv[2]
	quiet = len(sys.argv) > 3 and sys.argv[3] == '-q'
	nburn = 2_000_000

	NPI_dates = pd.read_csv('../data/NPI_dates.csv')

	if country == 'all':
		countries = [country.replace(' ', '_') for country in NPI_dates['Country']]
	else:
		countries = [country]

	if country=='Spain':
		delete_chain_less_than = 15 #TODO add as input parameter
	else:
		delete_chain_less_than = None

	for country in countries:
		npz_path = os.path.join(output_folder, 'inference', '{}.npz').format(country)
		print("Loading inference data from", npz_path)
		data = np.load(npz_path)

		var_names = list(data['var_names'])
		start_date = data['start_date']
		start_date = pd.to_datetime(start_date)
		sample = data['chain']
		log_posterior = data['lnprobability']

		if delete_chain_less_than:
			if len((sample[:,1_000_000, var_names.index('τ')]<delete_chain_less_than).nonzero())>1:
				raise AssertionError('too many bad chains')
			bad_chain_ind = (sample[:,1_000_000, var_names.index('τ')]<delete_chain_less_than).nonzero()[0][0]
			sample = np.delete(sample, bad_chain_ind, axis=0)
			log_posterior = np.delete(log_posterior, bad_chain_ind, axis=0)

		first_date = pd.to_datetime(NPI_dates.loc[NPI_dates['Country'] == country.replace('_', ' '), 'First'].values[0])
		last_date = pd.to_datetime(NPI_dates.loc[NPI_dates['Country'] == country.replace('_', ' '), 'Last'].values[0])
		first_date_days = (first_date - start_date).days
		last_date_days = (last_date - start_date).days

		τ_sample = sample[:, :, -1]

		thin = 100
		plt.plot(τ_sample[:, ::thin].T)
		plt.ylabel('τ')
		plt.xlabel('Iteration')
		plt.title('Trace plot: {}'.format(country));
		plt.axvline(nburn/thin, color='k');
		if not quiet: plt.show()

		plt.plot(τ_sample[:,nburn::thin].T)
		plt.ylabel('τ')
		plt.xlabel('Iteration')
		plt.title('Trace plot (burnin): {}'.format(country));
		if not quiet: plt.show()

		fig, ax = plt.subplots()
		for i in range(τ_sample.shape[0]):
		    ax.hist(τ_sample[i, nburn:], density=False, alpha=0.4)
		ax.set(xlabel='τ', ylabel='Posterior probability')
		ax.set_title('Posterior: {}'.format(country))
		if not quiet: plt.show()

		idx_max = log_posterior[:, nburn:].ravel().argmax()
		τ_sample = τ_sample[:, nburn:].ravel()
		τ_map = τ_sample[idx_max]
		print("τ MAP = {:.2f}".format(τ_map), int_to_dt(τ_map))
		τ_mean = np.mean(τ_sample)
		print("τ mean = {:.2f}".format(τ_mean), int_to_dt(τ_mean))
		τ_med = np.median(τ_sample)
		print("τ median = {:.2f}".format(τ_med), int_to_dt(τ_med))
		quantile = 0.95
		τ_ci = np.quantile(abs(τ_sample - τ_med), quantile)
		print("+-{:.2f} days for {:.2f} quantile".format(τ_ci, quantile))
		# hpi_ = hpd(τ_sample, quantile)#, multimodal=True)
		# print("HPI {:.2f}-{:.2f} for {:.2f} confidence".format(*hpi_, quantile))
		τ_hpd_from, τ_hpd_to = calc_hpd(τ_sample, quantile)
		# print(τ_hpd_from,τ_hpd_to)
		# print(hpd(τ_sample, quantile)) #arviz hpd is the same

		fig, ax = plt.subplots(figsize=(6, 3))
		xmin = τ_sample.min()-2.5
		xmax = max(last_date_days, np.ceil(τ_sample.max())) + 2.5
		density, bins, _ = ax.hist(τ_sample, bins=np.arange(0, xmax, 1), density=True, align='mid')
		ymax = density.max() * 1.2
		ax.axvline(τ_med, color='k')
		# ax.fill_between([τ_med - τ_ci, τ_med + τ_ci], 0, ymax, alpha=0.4, color='k')
		ax.fill_between([τ_hpd_from, τ_hpd_to], 0, ymax, alpha=0.4, color='k')
		# for hpi_ in hpis:
		# ax.fill_between([hpi_[0], hpi_[1]], 0, ymax, alpha=0.4, color='k')
		# ax.axvline(first_date_days, color=red)
		ax.axvline(last_date_days, color=red)
		# ax.fill_between([first_date_days, last_date_days], 0, ymax, alpha=0.4, color=red)
		ax.set(xlabel=r'Effective start of NPI, $\tau$',
		       ylabel=r'Posterior density, $P(\tau \mid \vec{X})$',
		       xlim=(τ_sample.min()-1, τ_sample.max()+1),
		       ylim=(0, ymax)
		)
		days = np.arange(xmin, xmax, 3)
		deltas = [timedelta(int(x)) for x in days]
		ax.set_xticks(days);
		txt = ax.set_xticklabels([(start_date + d).strftime('%b %d') for d in deltas], rotation=45)

		ax.annotate(r'$\hat{\tau}$', (τ_med+0.1, ymax*0.875))
		# ax.annotate(r'$\tau^0$', (first_date_days+0.1, ymax*0.875))
		ax.annotate(r'$\tau^*$', (last_date_days+0.1, ymax*0.875))

		ax.set_title(country.replace('_', ' '))
		sns.despine()
		if not quiet: plt.show()
		fig_filename = os.path.join(output_folder, 'figures', '{}_τ_posterior_hpd.pdf'.format(country))
		print("Saving to", fig_filename)
		fig.savefig(fig_filename, dpi=100, **savefig_bbox(*txt))

max_stars_repo_name: yoavram-lab/EffectiveNPI

max_issues_repo_name: yoavram-lab/EffectiveNPI

max_forks_repo_name: yoavram-lab/EffectiveNPI


text:
import numpy as np
import pandas as pd
import utils
import copy

import config

class RectangularWorld():

    def __init__(self, noise_location, round_length = 6, full_noise_file = True, center_radius = None, edge_goal = False):

        if noise_location is not None:
            if full_noise_file:
                self.noise_location = noise_location + '/'
                with open(self.noise_location + 't0.csv') as f:
                    self.noise_line_width = len(f.readline())
                assert not edge_goal
            else:
                self.centers = np.array(pd.read_csv(noise_location))

        self.center_radius = center_radius
        self.edge_goal = edge_goal

        self.full_noise_file = full_noise_file

        self.world = config.WORLD
        self.shape = 'rectangle'

        self.tick_frequency = 125
        self.ticks_per_sec = 1000/125
        self.us_min_wage_per_tick = 7.25 / (60*60*(1000 / self.tick_frequency))
        self.round_length = round_length
        self.max_bonus = 1.25*6/self.round_length
        self.game_length = int(self.round_length * 60 * self.ticks_per_sec)

        self.min_speed = 17 / float(self.ticks_per_sec)
        self.max_speed = 57 / float(self.ticks_per_sec)

        self.size = config.SIZE

        self.pos_limits = config.POS_LIMITS

    def get_random_position(self):
        """
        >>> np.random.seed(1)
        >>> w = RectangularWorld('../test/')
        >>> np.round(np.max([w.get_random_position()[0] for i in range(10000)]))
        482.0
        >>> round(np.min([w.get_random_position()[0] for i in range(10000)]))
        3.0
        >>> round(np.min([w.get_random_position()[1] for i in range(10000)]))
        3.0
        >>> round(np.max([w.get_random_position()[1] for i in range(10000)]))
        277.0
        """

        pos = None

        if self.edge_goal:
            while True:
                pos = np.array([np.random.uniform(self.pos_limits["x_min"], self.pos_limits["x_max"]),
                                np.random.uniform(self.pos_limits["y_min"], self.pos_limits["y_max"])])
                if utils.check_collision(pos, self.pos_limits, self.shape, update = False, extended = True):
                    break
        else:
            pos = np.array([np.random.uniform(self.pos_limits["x_min"], self.pos_limits["x_max"]),
                            np.random.uniform(self.pos_limits["y_min"], self.pos_limits["y_max"])])

        return pos

    def get_random_angle(self):
        """
        >>> np.random.seed(1)
        >>> w = RectangularWorld('../test/')
        >>> np.max([w.get_random_angle() for i in range(10000)])
        359.0
        >>> np.min([w.get_random_angle() for i in range(10000)])
        0.0
        """
        return np.floor(np.random.random() * 360)

    def get_score(self, pos, time):
        """
        >>> w = RectangularWorld('../test/')
        >>> w.get_score(np.array([2.1,2.9]), 0)
        0.0
        """

        if self.full_noise_file:
            return utils.get_score(pos, time, self.noise_location, self.pos_limits, self.noise_line_width)
        else:
            if self.edge_goal:
                return utils.wall_score(pos, self.centers[time], self.center_radius, self.pos_limits, self.shape)
            else:
                return utils.calculate_score(pos, self.centers[time], self.center_radius, self.pos_limits, self.shape)

if __name__ == "__main__":
    import doctest
    doctest.testmod()

max_stars_repo_name: hawkrobe/fish

max_issues_repo_name: hawkrobe/fish


text:
# -*- coding: utf-8 -*-
# This file is part of the pyMOR project (http://www.pymor.org).
# Copyright Holders: Rene Milk, Stephan Rave, Felix Schindler
# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)

import numpy as np

from pymor.operators.basic import OperatorBase
from pymor.operators.interfaces import OperatorInterface
from pymor.vectorarrays.block import BlockVectorArray
from pymor.vectorarrays.interfaces import VectorSpace

class BlockOperator(OperatorBase):
    """A sparse matrix of arbitrary |Operators|.

    This operator can be :meth:`applied <pymor.operators.interfaces.OperatorInterface.apply>`
    to :class:`BlockVectorArrays <pymor.vectorarrays.block.BlockVectorArray>` of an
    appropriate :attr:`~pymor.vectorarrays.interfaces.VectorArrayInterface.subtype`.

    Parameters
    ----------
    blocks
        Two-dimensional |NumPy| array where each entry is an |Operator| or `None`.
    """

    def _operators(self):
        """Iterate over operators (not None)"""
        for row in self._blocks:
            for entry in row:
                if entry is not None:
                    yield entry

    def __init__(self, blocks):
        blocks = np.array(blocks)
        assert isinstance(blocks, np.ndarray) and blocks.ndim == 2
        self._blocks = blocks

        assert all(isinstance(op, OperatorInterface) for op in self._operators())
        assert all(any(self._blocks[i, j] is not None for j in range(self._blocks.shape[1]))
                   for i in range(self._blocks.shape[0]))
        assert all(any(self._blocks[i, j] is not None for i in range(self._blocks.shape[0]))
                   for j in range(self._blocks.shape[1]))

        source_types = [None for j in range(self._blocks.shape[1])]
        range_types = [None for i in range(self._blocks.shape[0])]
        for (i, j), op in np.ndenumerate(self._blocks):
            if op is not None:
                assert source_types[j] is None or op.source == source_types[j]
                source_types[j] = op.source
                assert range_types[i] is None or op.range == range_types[i]
                range_types[i] = op.range

        self.source = VectorSpace(BlockVectorArray, tuple(source_types))
        self.range = VectorSpace(BlockVectorArray, tuple(range_types))
        self._source_dims = tuple(space.dim for space in self.source.subtype)
        self._range_dims = tuple(space.dim for space in self.range.subtype)
        self.num_source_blocks = len(source_types)
        self.num_range_blocks = len(range_types)
        self.linear = all(op.linear for op in self._operators())
        self.build_parameter_type(inherits=list(self._operators()))

    @classmethod
    def hstack(cls, operators):
        """Horizontal stacking of |Operators|.

        Parameters
        ----------
        operators
            A tuple, list, array, or iterable of |Operators|.
        """
        blocks = np.array([[op for op in operators]])
        return cls(blocks)

    @classmethod
    def vstack(cls, operators):
        """Vertical stacking of |Operators|.

        Parameters
        ----------
        operators
            A tuple, list, array, or iterable of |Operators|.
        """
        blocks = np.array([[op] for op in operators])
        return cls(blocks)

    def apply(self, U, ind=None, mu=None):
        assert U in self.source
        assert U.check_ind(ind)

        blocks = [None for i in range(self.num_range_blocks)]
        for i in range(self.num_range_blocks):
            for j in range(self.num_source_blocks):
                op = self._blocks[i, j]
                if op is not None:
                    V = op.apply(U.block(j), ind=ind, mu=mu)
                    if blocks[i] is None:
                        blocks[i] = V
                    else:
                        blocks[i] += V

        return BlockVectorArray(blocks)

    def apply_adjoint(self, U, ind=None, mu=None, source_product=None, range_product=None):
        assert U in self.range
        assert U.check_ind(ind)
        assert source_product is None or source_product.source == source_product.range == self.source
        assert range_product is None or range_product.source == range_product.range == self.range

        if range_product is not None:
            U = range_product.apply(U)

        blocks = [None for j in range(self.num_source_blocks)]
        for j in range(self.num_source_blocks):
            for i in range(self.num_range_blocks):
                op = self._blocks[i, j]
                if op is not None:
                    V = op.apply_adjoint(U.block(i), ind=ind, mu=mu)
                    if blocks[j] is None:
                        blocks[j] = V
                    else:
                        blocks[j] += V

        V = BlockVectorArray(blocks)
        if source_product is not None:
            V = source_product.apply_inverse(V)

        return V

class BlockDiagonalOperator(BlockOperator):
    """Block diagonal |Operator| of arbitrary |Operators|.

    This is a specialization of :class:`BlockOperator` for the
    block diagonal case.
    """

    def __init__(self, blocks):
        n = len(blocks)
        blocks2 = np.array([[None for j in range(n)] for i in range(n)])
        for i, op in enumerate(blocks):
            blocks2[i, i] = op
        super().__init__(blocks2)

    def apply_inverse(self, V, ind=None, mu=None, least_squares=False):
        assert V in self.range
        assert V.check_ind(ind)

        U = [None for i in range(self.num_source_blocks)]
        for i in range(self.num_source_blocks):
            U[i] = self._blocks[i, i].apply_inverse(V.block(i), ind=ind, mu=mu, least_squares=least_squares)

        return BlockVectorArray(U)

    def apply_inverse_adjoint(self, U, ind=None, mu=None, source_product=None, range_product=None,
                              least_squares=False):
        assert U in self.source
        assert U.check_ind(ind)
        assert source_product is None or source_product.source == source_product.range == self.source
        assert range_product is None or range_product.source == range_product.range == self.range

        if source_product or range_product:
            return super().apply_inverse_adjoint(
                U, ind=ind, mu=mu, source_product=source_product, range_product=range_product)

        V = [None for i in range(self.num_source_blocks)]
        for i in range(self.num_source_blocks):
            V[i] = self._blocks[i, i].apply_inverse_adjoint(U.block(i), ind=ind, mu=mu, least_squares=least_squares)

        return BlockVectorArray(V)

max_stars_repo_name: JuliaBru/pymor

max_issues_repo_name: JuliaBru/pymor

max_forks_repo_name: JuliaBru/pymor


text:
''' Functions

Summary:
    This file contains

Example:
    Usage of

Todo:
    *
'''

''' Imports '''

# data array handling
import numpy as np

# data table handling
import pandas as pd

''' Current-Voltage File Format Parse Functions '''

def parse_format_halm(_file_path):

    ''' Parse HALM Current-Voltage Format

        Import current-voltage measurement settings and data from a HALM format file

    # inputs
        _file_path (str): full filepath

    Returns:
        dict: extracted data and parameters
    '''

    # file extract storage
    _headers = []; _data = []

    # import data segments and headers from txt file, store in lists
    with open(_file_path, 'r', encoding = 'iso-8859-1') as file:

        lines = file.readlines()
        # row numbers [header, data start, n rows] for data segments, iterate
        segs = [
            [12, 13, 199], # 1.0 suns IV
            [12, 217, 199], # 0.5 suns IV
            [434, 435, 100], # dark reverse bias IV
            [548, 549, 200], # shunt res
            [762, 763, 200], # series res
        ]
        for s in segs:
            head, start, rows = s
            _headers.append(lines[head])
            _data.append([ lines[i] for i in range(len(lines))[start:start+rows] ])

    # parse header string, discard Nr, split for (name, units) by column, store array
    for i in range(len(_headers)):
        _headers[i] = np.array([ [head[str.index(head, ']')+1:],
                                 head[str.index(head, '[')+1:str.index(head, ']')]]
                                for head in str.split(str.strip(_headers[i]), '\t')[1:] ])

    # parse data rows, discard Nr, split each column, parse to float else zero, store array
    chk = lambda s: str.isdigit(str.join('',[ c for c in s if c not in ['.', '-', '+'] ]))
    for i in range(len(_data)):
        _data[i] = np.stack([ [ float(v) if chk(v) else 0.
                               for v in str.split(str.strip(row), '\t')[1:] ]
                             for row in _data[i] ])

    names = ['full', 'half', 'dark', 'shunt', 'series']

    # store each data array in dict with header as key
    #data_entry = { _headers[i][0][0]:_data[i][:,:] for i in list(range(len(_headers))) }
    data_entry = { names[i]: { 'voltage': _data[i][:,6][::-1],
                               'current': -_data[i][:,7][::-1],
                               'intensity': _data[i][:,8][::-1], }
        for i in range(len(names)) }

    # return data dict
    return data_entry

def parse_format_wavelabs(_file_path):

    ''' Parse WaveLabs Current-Voltage Format

        Import current-voltage measurement settings and data from a WaveLabs IV format file

    # inputs
        _file_path (str): full filepath

    Returns:
        dict: extracted data and parameters
    '''

    # read in 4 header rows, convert to numpy array, discard null 3rd row, discard first 3 columns
    jv_head = pd.read_excel(io = _file_path, sheet_name = 'IV-Raw', header = None, nrows = 4).values[[0, 1, 3], 3:]

    # read raw jv data from file, convert to numpy array, skip 4 header rows, discard first 3 columns
    jv_raw = pd.read_excel(io = _file_path, sheet_name = 'IV-Raw', header = None, skiprows = 4,
        index_col = None).values[:, 3:]

    # get index of data columns (discard null cols)
    jv_datas_mask = [ i-1 for i in range(jv_head.shape[1]+1) if i%3 != 0 ]

    # get data columns (discard null cols)
    jv_datas = jv_raw[:, jv_datas_mask].astype(np.float32)

    # get name of each data column pair
    jv_data_names = jv_head[1, jv_datas_mask[1::2]].astype(np.str)

    # get step name of each data column pair
    jv_data_steps = jv_head[0, jv_datas_mask[1::2]].astype(np.str)

    # get units of each data column
    jv_data_units = jv_head[2, jv_datas_mask].astype(np.str)

    # set new unit labels
    relabel = {'Current [A]':'current', 'E':'intensity', 'Voltage [V]':'voltage'}

    # iterate and update unit labels
    for i in range(len(jv_data_units)):
        if jv_data_units[i] in relabel.keys():
            jv_data_units[i] = relabel[ jv_data_units[i] ]

    # set new step labels
    relabel = {'1 sun': 'full', '1/2 sun': 'half', 'SunsVoc': 'suns_voc', 'UV light': 'uv', 'IR light': 'ir',
        'dark': 'dark'}

    # iterate and update unit labels
    for i in range(len(jv_data_steps)):
        if jv_data_steps[i] in relabel.keys():
            jv_data_steps[i] = relabel[ jv_data_steps[i] ]

    # initialise data storage dict
    data = {}

    # iterate steps
    for i in range(len(jv_datas[1])):

        # unpack data and headers
        _data = jv_datas[:,i]
        step = jv_data_steps[int(i/2)]
        name = jv_data_names[int(i/2)]
        unit = jv_data_units[i]

        # only keep raw current, voltage, intensity data for select steps
        if step in ['full', 'half', 'suns_voc', 'dark']:

            # only keep raw current, voltage, intensity data for select steps
            if (name in ['Raw IV Data', 'Intensity']) and (unit in ['current', 'voltage', 'intensity']):

                # add step to dict
                if not step in data.keys():
                    data[step] = {}

                # store data array in dict by units, strip tail nan
                data[step][unit] = _data[np.where( np.invert(np.isnan(_data)) )]

    # return data dict
    return data

''' File Import and Processing Functions '''

def process_raw_data(_raw_data, _filter = False):

    ''' Process Raw Data

        Process raw data parsed from file, file format independent; perform calculations, filter for return only
        required data and parameters

    Args:
        _raw_data (dict): raw data parsed from file
        _filter (bool): filter processed raw data or return complete dataset

    Returns:
        dict: processed raw data
    '''

    # filter for only required data and parameters to return
    if _filter:

        # define list of required
        required = []

        # return processed data, filtered by required
        #return { key: value for key, value in _raw_data.items() if key in required }
        pass

    else:
        # return complete raw data
        return _raw_data

def parse_current_voltage(_file_path, _format):

    ''' Parse File

        Parse source data file of given format, run post parse processing, and return imported data

    Args:
        _file_path (str): full file path including name with extension
        _format (dict): data file format information

    Returns:
        dict: parsed, processed data and parameters from file
    '''

    # HALM IV, SIRF
    if _format['equipment'] == 'halm':

        # process raw data, format independent, return result
        return process_raw_data( parse_format_halm(_file_path) )

    # WaveLabs Solar Simulator, 140 TETB
    if _format['equipment'] == 'wavelabs':

        # process raw data, format independent, return result
        return process_raw_data( parse_format_wavelabs(_file_path) )


text:
import os

import numpy as np
import torch
import torch.nn.functional as F
from torch import nn
from collections import OrderedDict
from torch.autograd import Variable

def gram_matrix(y):
    (b, ch, h, w) = y.size()
    features = y.view(b, ch, w * h)
    features_t = features.transpose(1, 2)
    gram = features.bmm(features_t) / (ch * h * w)
    return gram

def poly_lr_scheduler(optimizer, init_lr, iter, lr_decay_iter=1, max_iter=30000, power=0.9,):
    """Polynomial decay of learning rate
        :param init_lr is base learning rate
        :param iter is a current iteration
        :param lr_decay_iter how frequently decay occurs, default is 1
        :param max_iter is number of maximum iterations
        :param power is a polymomial power
    """
    if iter % lr_decay_iter or iter > max_iter:
        return optimizer

    for param_group in optimizer.param_groups:
        tmp = (1 - iter/max_iter)**power
        param_group['lr'] = init_lr*tmp

def wct(content_feat, style_feat):
    content_feat = content_feat.data.cpu()
    content_feat = content_feat.squeeze(0).double()

    style_feat = style_feat.data.cpu()
    style_feat = style_feat.squeeze(0).double()

    C, W, H = content_feat.size()

    transfered = whiten_and_color(content_feat.view(C, -1),
            style_feat.view(C, -1))

    transfered = transfered.view_as(content_feat).float().unsqueeze(0)

    return Variable(transfered).cuda()

def whiten_and_color(cF,sF):
    cFSize = cF.size()
    c_mean = torch.mean(cF,1) # c x (h x w)
    c_mean = c_mean.unsqueeze(1).expand_as(cF)
    cF = cF - c_mean

    contentConv = torch.mm(cF,cF.t()).div(cFSize[1]-1) + torch.eye(cFSize[0]).double()
    c_u,c_e,c_v = torch.svd(contentConv,some=False)

    k_c = cFSize[0]
    for i in range(cFSize[0]):
        if c_e[i] < 0.00001:
            k_c = i
            break

    sFSize = sF.size()
    s_mean = torch.mean(sF,1)
    sF = sF - s_mean.unsqueeze(1).expand_as(sF)
    styleConv = torch.mm(sF,sF.t()).div(sFSize[1]-1)
    s_u,s_e,s_v = torch.svd(styleConv,some=False)

    k_s = sFSize[0]
    for i in range(sFSize[0]):
        if s_e[i] < 0.00001:
            k_s = i
            break

    c_d = (c_e[0:k_c]).pow(-0.5)
    step1 = torch.mm(c_v[:,0:k_c],torch.diag(c_d))
    step2 = torch.mm(step1,(c_v[:,0:k_c].t()))
    whiten_cF = torch.mm(step2,cF)

    s_d = (s_e[0:k_s]).pow(0.5)
    targetFeature = torch.mm(torch.mm(torch.mm(s_v[:,0:k_s],torch.diag(s_d)),(s_v[:,0:k_s].t())),whiten_cF)
    targetFeature = targetFeature + s_mean.unsqueeze(1).expand_as(targetFeature)

    return targetFeature

def calc_mean_std(feat, eps=1e-5):
    # eps is a small value added to the variance to avoid divide-by-zero.
    size = feat.data.size()
    assert (len(size) == 4)
    N, C = size[:2]
    feat_var = feat.view(N, C, -1).var(dim=2) + eps
    feat_std = feat_var.sqrt().view(N, C, 1, 1)
    feat_mean = feat.view(N, C, -1).mean(dim=2).view(N, C, 1, 1)
    return feat_mean, feat_std

def adaptive_instance_normalization(content_feat, style_feat):
    assert (content_feat.data.size()[:2] == style_feat.data.size()[:2])
    size = content_feat.data.size()
    style_mean, style_std = calc_mean_std(style_feat)
    content_mean, content_std = calc_mean_std(content_feat)

    normalized_feat = (content_feat - content_mean.expand(
        size)) / content_std.expand(size)
    return normalized_feat * style_std.expand(size) + style_mean.expand(size)

def load_model_filter(model, snapshot, prefix=False):
    pretrained_dict = torch.load(snapshot)
    if prefix:
        new_state_dict = OrderedDict()
        for k, v in pretrained_dict.items():
            name = k[7:] # remove `enc.` or `dec.`
            new_state_dict[name] = v
        pretrained_dict = new_state_dict

    model_dict = model.state_dict()

    # 1. filter out unnecessary keys
    pretrained_dict = {k: v for k, v in pretrained_dict.items() if k in model_dict}
    # 2. overwrite entries in the existing state dict
    model_dict.update(pretrained_dict)
    # 3. load the new state dict
    model.load_state_dict(pretrained_dict)

    return model

def check_mkdir(dir_name):
    if not os.path.exists(dir_name):
        os.mkdir(dir_name)

def initialize_weights(*models):
    for model in models:
        for module in model.modules():
            if isinstance(module, nn.Conv2d) or isinstance(module, nn.Linear):
                nn.init.kaiming_normal(module.weight)
                if module.bias is not None:
                    module.bias.data.zero_()
            elif isinstance(module, nn.BatchNorm2d):
                module.weight.data.fill_(1)
                module.bias.data.zero_()

def get_upsampling_weight(in_channels, out_channels, kernel_size):
    factor = (kernel_size + 1) // 2
    if kernel_size % 2 == 1:
        center = factor - 1
    else:
        center = factor - 0.5
    og = np.ogrid[:kernel_size, :kernel_size]
    filt = (1 - abs(og[0] - center) / factor) * (1 - abs(og[1] - center) / factor)
    weight = np.zeros((in_channels, out_channels, kernel_size, kernel_size), dtype=np.float64)
    weight[range(in_channels), range(out_channels), :, :] = filt
    return torch.from_numpy(weight).float()

class CrossEntropyLoss2d(nn.Module):
    def __init__(self, weight=None, size_average=True, ignore_index=255):
        super(CrossEntropyLoss2d, self).__init__()
        self.nll_loss = nn.NLLLoss2d(weight, size_average, ignore_index)

    def forward(self, inputs, targets):
        return self.nll_loss(F.log_softmax(inputs), targets)

def cross_entropy2d(input, target, weight=None, size_average=True, ignore_index=255):
    n, c, h, w = input.size()
    log_p = F.log_softmax(input, dim=1)
    log_p = log_p.transpose(1, 2).transpose(2, 3).contiguous().view(-1, c)
    log_p = log_p[target.view(n * h * w, 1).repeat(1, c) >= 0]
    log_p = log_p.view(-1, c)

    mask = target >= 0
    target = target[mask]
    loss = F.nll_loss(log_p, target, ignore_index=ignore_index,
                      weight=weight, size_average=False)
    if size_average:
        loss /= mask.data.sum()
    return loss

def bootstrapped_cross_entropy2d(input, target, K, weight=None, size_average=True):

    batch_size = input.size()[0]

    def _bootstrap_xentropy_single(input, target, K, weight=None, size_average=True):
        n, c, h, w = input.size()
        log_p = F.log_softmax(input, dim=1)
        log_p = log_p.transpose(1, 2).transpose(2, 3).contiguous().view(-1, c)
        log_p = log_p[target.view(n * h * w, 1).repeat(1, c) >= 0]
        log_p = log_p.view(-1, c)

        mask = target >= 0
        target = target[mask]
        loss = F.nll_loss(log_p, target, weight=weight, ignore_index=255,
                          reduce=False, size_average=False)
        topk_loss, _ = loss.topk(K)
        reduced_topk_loss = topk_loss.sum() / K

        return reduced_topk_loss

    loss = 0.0
    # Bootstrap from each image not entire batch
    for i in range(batch_size):
        loss += _bootstrap_xentropy_single(input=torch.unsqueeze(input[i], 0),
                                           target=torch.unsqueeze(target[i], 0),
                                           K=K,
                                           weight=weight,
                                           size_average=size_average)
    return loss / float(batch_size)

def _fast_hist(label_pred, label_true, num_classes):
    mask = (label_true >= 0) & (label_true < num_classes)
    hist = np.bincount(
        num_classes * label_true[mask].astype(int) +
        label_pred[mask], minlength=num_classes ** 2).reshape(num_classes, num_classes)
    return hist

### T-SNE
def Hbeta(D=np.array([]), beta=1.0):
    """
        Compute the perplexity and the P-row for a specific value of the
        precision of a Gaussian distribution.
    """

    # Compute P-row and corresponding perplexity
    P = np.exp(-D.copy() * beta)
    sumP = sum(P)
    H = np.log(sumP) + beta * np.sum(D * P) / sumP
    P = P / sumP
    return H, P

def x2p(X=np.array([]), tol=1e-5, perplexity=30.0):
    """
        Performs a binary search to get P-values in such a way that each
        conditional Gaussian has the same perplexity.
    """

    # Initialize some variables
    print("Computing pairwise distances...")
    (n, d) = X.shape
    sum_X = np.sum(np.square(X), 1)
    D = np.add(np.add(-2 * np.dot(X, X.T), sum_X).T, sum_X)
    P = np.zeros((n, n))
    beta = np.ones((n, 1))
    logU = np.log(perplexity)

    # Loop over all datapoints
    for i in range(n):

        # Print progress
        if i % 500 == 0:
            print("Computing P-values for point %d of %d..." % (i, n))

        # Compute the Gaussian kernel and entropy for the current precision
        betamin = -np.inf
        betamax = np.inf
        Di = D[i, np.concatenate((np.r_[0:i], np.r_[i+1:n]))]
        (H, thisP) = Hbeta(Di, beta[i])

        # Evaluate whether the perplexity is within tolerance
        Hdiff = H - logU
        tries = 0
        while np.abs(Hdiff) > tol and tries < 50:

            # If not, increase or decrease precision
            if Hdiff > 0:
                betamin = beta[i].copy()
                if betamax == np.inf or betamax == -np.inf:
                    beta[i] = beta[i] * 2.
                else:
                    beta[i] = (beta[i] + betamax) / 2.
            else:
                betamax = beta[i].copy()
                if betamin == np.inf or betamin == -np.inf:
                    beta[i] = beta[i] / 2.
                else:
                    beta[i] = (beta[i] + betamin) / 2.

            # Recompute the values
            (H, thisP) = Hbeta(Di, beta[i])
            Hdiff = H - logU
            tries += 1

        # Set the final row of P
        P[i, np.concatenate((np.r_[0:i], np.r_[i+1:n]))] = thisP

    # Return final P-matrix
    print("Mean value of sigma: %f" % np.mean(np.sqrt(1 / beta)))
    return P

def pca(X=np.array([]), no_dims=50):
    """
        Runs PCA on the NxD array X in order to reduce its dimensionality to
        no_dims dimensions.
    """

    print("Preprocessing the data using PCA...")
    (n, d) = X.shape
    X = X - np.tile(np.mean(X, 0), (n, 1))
    (l, M) = np.linalg.eig(np.dot(X.T, X))
    Y = np.dot(X, M[:, 0:no_dims])
    return Y

def tsne(X=np.array([]), no_dims=2, initial_dims=50, perplexity=30.0, max_iter=1000):
    """
        Runs t-SNE on the dataset in the NxD array X to reduce its
        dimensionality to no_dims dimensions. The syntaxis of the function is
        `Y = tsne.tsne(X, no_dims, perplexity), where X is an NxD NumPy array.
    """

    # Check inputs
    if isinstance(no_dims, float):
        print("Error: array X should have type float.")
        return -1
    if round(no_dims) != no_dims:
        print("Error: number of dimensions should be an integer.")
        return -1

    # Initialize variables
    X = pca(X, initial_dims).real
    (n, d) = X.shape
    max_iter = max_iter
    initial_momentum = 0.5
    final_momentum = 0.8
    eta = 500
    min_gain = 0.01
    Y = np.random.randn(n, no_dims)
    dY = np.zeros((n, no_dims))
    iY = np.zeros((n, no_dims))
    gains = np.ones((n, no_dims))

    # Compute P-values
    P = x2p(X, 1e-5, perplexity)
    P = P + np.transpose(P)
    P = P / np.sum(P)
    P = P * 4.									# early exaggeration
    P = np.maximum(P, 1e-12)

    # Run iterations
    for iter in range(max_iter):

        # Compute pairwise affinities
        sum_Y = np.sum(np.square(Y), 1)
        num = -2. * np.dot(Y, Y.T)
        num = 1. / (1. + np.add(np.add(num, sum_Y).T, sum_Y))
        num[range(n), range(n)] = 0.
        Q = num / np.sum(num)
        Q = np.maximum(Q, 1e-12)

        # Compute gradient
        PQ = P - Q
        for i in range(n):
            dY[i, :] = np.sum(np.tile(PQ[:, i] * num[:, i], (no_dims, 1)).T * (Y[i, :] - Y), 0)

        # Perform the update
        if iter < 20:
            momentum = initial_momentum
        else:
            momentum = final_momentum
        gains = (gains + 0.2) * ((dY > 0.) != (iY > 0.)) + \
                (gains * 0.8) * ((dY > 0.) == (iY > 0.))
        gains[gains < min_gain] = min_gain
        iY = momentum * iY - eta * (gains * dY)
        Y = Y + iY
        Y = Y - np.tile(np.mean(Y, 0), (n, 1))

        # Compute current value of cost function
        if (iter + 1) % 10 == 0:
            C = np.sum(P * np.log(P / Q))
            print("Iteration %d: error is %f" % (iter + 1, C))

        # Stop lying about P-values
        if iter == 100:
            P = P / 4.

    # Return solution
    return Y

def evaluate(predictions, gts, num_classes):
    hist = np.zeros((num_classes, num_classes))
    for lp, lt in zip(predictions, gts):
        hist += _fast_hist(lp.flatten(), lt.flatten(), num_classes)
    # axis 0: gt, axis 1: prediction
    acc = np.diag(hist).sum() / hist.sum()
    acc_cls = np.diag(hist) / hist.sum(axis=1)
    acc_cls = np.nanmean(acc_cls)
    iu = np.diag(hist) / (hist.sum(axis=1) + hist.sum(axis=0) - np.diag(hist))
    mean_iu = np.nanmean(iu)
    freq = hist.sum(axis=1) / hist.sum()
    fwavacc = (freq[freq > 0] * iu[freq > 0]).sum()
    return acc, acc_cls, mean_iu, fwavacc, iu

class AverageMeter(object):
    def __init__(self):
        self.reset()

    def reset(self):
        self.val = 0
        self.avg = 0
        self.sum = 0
        self.count = 0

    def update(self, val, n=1):
        self.val = val
        self.sum += val * n
        self.count += n
        self.avg = self.sum / self.count

class PolyLR(object):
    def __init__(self, optimizer, curr_iter, max_iter, lr_decay):
        self.max_iter = float(max_iter)
        self.init_lr_groups = []
        for p in optimizer.param_groups:
            self.init_lr_groups.append(p['lr'])
        self.param_groups = optimizer.param_groups
        self.curr_iter = curr_iter
        self.lr_decay = lr_decay

    def step(self):
        for idx, p in enumerate(self.param_groups):
            p['lr'] = self.init_lr_groups[idx] * (1 - self.curr_iter / self.max_iter) ** self.lr_decay

class LogFile:
    def __init__(self, fl):
        open(fl,'w').close()
        self.fl = fl

    def log(self, log_str):
        with open(self.fl, 'a') as f:
            f.write(log_str+'\n')


text:
import os
from shutil import rmtree
import random
import string
import atexit
import warnings
from typing import Union, Sequence
from collections import OrderedDict
from tempfile import NamedTemporaryFile
# disable for now
# from ecogdata.parallel.mproc import Process
import numpy as np
import h5py
from scipy.signal import lfilter_zi, lfilter
from numpy.linalg import LinAlgError
# from tqdm import tqdm

from ecogdata.expconfig import load_params
from ecogdata.util import ToggleState
import ecogdata.parallel.sharedmem as shm

from .basic import ElectrodeDataSource, calc_new_samples, PlainArraySource
from .array_abstractions import HDF5Buffer, BufferBinder, slice_to_range, slice_data_buffer, BackgroundRead

__all__ = ['TempFilePool', 'MappedSource', 'MemoryBlowOutError', 'downsample_and_load', 'bfilter']

# Make one unique temp path per process
_rand_temp_path = 'MAPPED_TEMPFILES_' + ''.join([random.choice(string.ascii_letters) for _ in range(8)])

class TempFilePool:
    """
    Create a temporary file in the temp "pool" that can be closed and reopened without going away. The entire pool of
    temporary files will be cleared when the process exits.

    Example usage:

    >>> with TempFilePool() as f:
    >>>    name = str(f)
    >>> with h5py.File(name, 'w') as hdf_write:
    >>>    hdf_write.create_dataset(...)
    >>> hdf_read = h5py.File(name, 'r')

    When the resulting tempfile is meant to remain open (i.e. not in a file-open context), then
    use "register_to_close" to attempt closing before deleting

    >>> f.register_to_close(hdf_read)

    """

    pool_dir = os.path.join(os.path.abspath(os.path.curdir), _rand_temp_path)

    def __init__(self, *args, **kwargs):
        kwargs['dir'] = TempFilePool.pool_dir
        kwargs['delete'] = False
        if not os.path.exists(TempFilePool.pool_dir):
            os.makedirs(TempFilePool.pool_dir)
            with open(os.path.join(TempFilePool.pool_dir, 'README.txt'), 'w') as fw:
                fw.write('This directory contains temporary files for memory mapped arrays. '
                         'Under normal circumstances, this will be deleted when the process ends. '
                         'If you are reading this after a process has stopped, please delete!\n')
        self.tf = NamedTemporaryFile(*args, **kwargs)

    def __str__(self):
        return self.tf.name

    def __repr__(self):
        return repr(self.tf.name)

    def __enter__(self):
        self.tf.__enter__()
        return self

    def __exit__(self, exc, value, tb):
        return self.tf.__exit__(exc, value, tb)

    def register_to_close(self, file_handle):
        def closer():
            try:
                file_handle.close()
                print('closed file {}'.format(self))
            except Exception:
                print('Could not close file {}'.format(self))
                pass
        atexit.register(closer)

def _remove_pool():
    if os.path.exists(TempFilePool.pool_dir):
        print('Removing temp directory', TempFilePool.pool_dir)
        rmtree(TempFilePool.pool_dir)

atexit.register(_remove_pool)

class MemoryBlowOutError(Exception):
    """Raise this if a memory-mapped slice will be ginormous."""
    pass

class MappedSource(ElectrodeDataSource):
    # TODO:
    #  2. a mapped source should have a start/stop index that only exposes data inside a range

    def __init__(self, data_buffer: Union[HDF5Buffer, BufferBinder], electrode_field='data', electrode_channels=None,
                 channel_mask=None, aligned_arrays=None, transpose=False, raise_on_big_slice=True):
        """
        Provides a memory-mapped data source. This object should rarely be constructed by-hand, but it should not be
        overwhelmingly difficult to do so.

        Parameters
        ----------
        data_buffer: HDF5Buffer or BufferBinder
            A single- or multi-source buffer that has a MappedBuffer interface. Mapped array(s) may be Channel x Time
            or Time x Channel. Output slices from this object are always Channel x Time.
        electrode_field: str
            The dataset field name for electrode recording channels. Only needed if mirroring to a new file with
            identical dataset fields.
        electrode_channels: None or sequence
            The subset of channels in the electrode data array that contain electrode signals. A value of None
            indicates all channels contain electrode signals.
        channel_mask: boolean ndarray or None
            A binary mask the same length as electrode_channels or the number of rows in the signal array is
            electrode channels is None. The set of active electrodes is where the binary mask is True. The number of
            rows returned in a slice from this data source is sum(channel_mask). A value of None indicates all
            electrode channels are active.
        aligned_arrays: dict
            Any other datasets in source_file that should be aligned with the electrode signal array. These fields
            will be kept at the same sampling rate and index alignment as the electrode signal. They will also be
            preserved if this source is mirrored or copied to another source. Element of aigned_arrays are
            name:buffer pairs.
        units_scale: float or 2-tuple
            Either the scaling value or (offset, scaling) values such that signal = (memmap + offset) * scaling
        transpose: bool
            Is the mapped array stored in transpose (Time x Channels)?
        raise_on_big_slice: bool
            If True (default), then raise an exception if the proposed read slice will be larger than the memory

        """

        self._electrode_field = electrode_field
        # The data buffer object is an array slicing cache with possible units conversion.
        # This class uses mapping logic to expose only (active) electrode channels, depending on the
        # electrode_channels list and the current channel mask
        self.data_buffer = data_buffer
        # electrode_channels is a list of of data channel indices where electrode data are (immutable)
        if electrode_channels is None:
            n_channels = self.data_buffer.shape[1] if transpose else self.data_buffer.shape[0]
            self._electrode_channels = list(range(n_channels))
        else:
            self._electrode_channels = electrode_channels

        self._transpose = transpose
        self.dtype = self.data_buffer.dtype

        # channel_mask is the list of data indices for the active electrode set (mutable)
        self._active_channels = self._electrode_channels
        self.set_channel_mask(channel_mask)

        # this is toggle than when called creates a context manager with the supplied state (ToggleState)
        self._allow_big_slicing = ToggleState(init_state=False)
        self._raise_on_big_slice = raise_on_big_slice
        # Normally self[:10] will return an array. With this state toggled, it will return a new MappedSource
        self._slice_channels_as_maps = ToggleState(init_state=False)
        # Allow for a subprocess that will cache buffer reads in the background
        # (disabled for now)
        self._caching_process = None
        self._can_cache = False  # can't find out SWMR status from the dataset (gh #1580)
        # if data_buffer.subprocess_readable:
        #     self._can_cache = True
        # else:
        #     self._can_cache = False

        # The other aux fields will be setattr'd like
        # setattr(self, field, source_file[field]) ??
        if aligned_arrays is None:
            aligned_arrays = dict()
        for key in aligned_arrays:
            setattr(self, key, aligned_arrays[key])
        self.aligned_arrays = tuple(aligned_arrays.keys())

    @classmethod
    def from_hdf_sources(cls, source_files: Sequence[h5py.File], electrode_field, aligned_arrays=(), units_scale=None,
                         transpose=False, real_slices=True, **kwargs):
        """
        Constructs a MappedSource from one or more h5py.File objects. This builder creates the appropriate
        MappedBuffer-like object first.

        Parameters
        ----------
        source_files: Sequence[h5py.File]
            An opened HDF5 file: the file access mode will be respected.
        electrode_field: str
            The electrode recording channels are in source_file[electrode_field]. Mapped array may be Channel x Time
            or Time x Channel. Output slices from this object are always Channel x Time.
        aligned_arrays: sequence
            Any other datasets in source_files that should be aligned with the electrode signal array. These fields
            will be kept at the same sampling rate and index alignment as the electrode signal. They will also be
            preserved if this source is mirrored or copied to another source. Elements of the sequence can be dataset
            names in the mapped file, or (name, channel-list) pairs in the case that the aligned array is particular
            subset of channels in the same array (e.g. same as the electrode data).
        units_scale: float or 2-tuple
            Either the scaling value or (offset, scaling) values such that signal = (memmap + offset) * scaling
        transpose: bool
            Is the mapped array stored in transpose (Time x Channels)?
        kwargs:
            Further arguments detailing channel subsets, etc., for MappedSource

        Returns
        -------
        mapdata: MappedSource

        """

        if isinstance(source_files, h5py.File):
            source_files = (source_files,)

        # Set up underlying data buffer(s). Use a BufferBinder if there are multiple sources
        main_buffers = [HDF5Buffer(hdf[electrode_field], units_scale=units_scale, real_slices=real_slices)
                        for hdf in source_files]
        # Skip check for now -- SWMR issue pending
        # for b, hdf in zip(main_buffers, source_files):
        #     # if b.writeable and not hdf.swmr_mode:
        #     if not hdf.swmr_mode:
        #         try:
        #             hdf.swmr_mode = True
        #         except ValueError as e:
        #             if str(e).startswith('Unable to convert file format'):
        #                 print('Single-writer multiple-reader mode was not supported for this data file. '
        #                       'Block iteration on this MappedSource may be crashy.')
        #                 print('Original error:', str(e))
        #             else:
        #                 # unknown error -- reraise it
        #                 raise e

        concat_axis = 0 if transpose else 1
        if len(main_buffers) > 1:
            data_buffer = BufferBinder(main_buffers, axis=concat_axis)
        else:
            data_buffer = main_buffers[0]

        # Handle aligned arrays similarly to the "electrode_field" except no unit conversion
        array_names = OrderedDict()
        for field in aligned_arrays:
            if isinstance(field, str):
                aux_buffers = [HDF5Buffer(hdf[field]) for hdf in source_files]
                if len(aux_buffers) > 1:
                    array_names[field] = BufferBinder(aux_buffers, axis=concat_axis)
                else:
                    array_names[field] = aux_buffers[0]
            else:
                # this actually loads the channels to main memory??
                channels = np.s_[:, field[1]] if transpose else np.s_[field[1], :]
                field = field[0]
                arrays = [mb[channels] for mb in main_buffers]
                if len(arrays) > 1:
                    array_names[field] = np.concatenate(arrays, axis=concat_axis)
                else:
                    array_names[field] = arrays[0]
        return MappedSource(data_buffer, electrode_field=electrode_field, transpose=transpose,
                            aligned_arrays=array_names, **kwargs)

    def __len__(self):
        return self.shape[0]

    @property
    def shape(self):
        if self._transpose:
            return len(self._active_channels), self.data_buffer.shape[0]
        else:
            return len(self._active_channels), self.data_buffer.shape[1]

    @property
    def writeable(self):
        return self.data_buffer.writeable

    @property
    def big_slices(self):
        return self._allow_big_slicing

    @property
    def channels_are_maps(self):
        return self._slice_channels_as_maps

    @property
    def is_direct_map(self):
        all_active = self.binary_channel_mask.all()
        # quick short-circuit: False if there are masked channels
        if not all_active:
            return False
        num_disk_channels = self.data_buffer.shape[1] if self._transpose else self.data_buffer.shape[0]
        mapped_channels = np.array(self._electrode_channels)
        all_mapped = np.array_equal(mapped_channels, np.arange(num_disk_channels))
        return all_mapped

    @property
    def binary_channel_mask(self):
        if not hasattr(self, '_bmask'):
            self._bmask = np.ones(len(self._electrode_channels), '?')
        if self._bmask.sum() == len(self._active_channels):
            return self._bmask
        active_channels = set(self._active_channels)
        for n, c in enumerate(self._electrode_channels):
            self._bmask[n] = c in active_channels
        return self._bmask

    @property
    def _auto_block_length(self):
        # a bit hacky now...
        if isinstance(self.data_buffer, HDF5Buffer):
            chunks = self.data_buffer.chunks
        else:
            chunks = self.data_buffer.chunks[0]
        # This is the block size recommended by chunking
        block_size = chunk_size = chunks[0] if self._transpose else chunks[1]
        # If memory supports, then increase the block size to a multiple of the chunk size
        mem = load_params().memory_limit
        block_rows = self.shape[1] if self._transpose else self.shape[0]
        num_blocks = mem // (block_rows * block_size * self.dtype.itemsize)
        num_blocks = min(num_blocks, 100)
        if num_blocks < 1:
            warnings.warn('Memory limit is too low to support minimum block size', RuntimeWarning)
        else:
            block_size = num_blocks * chunk_size
        # Don't return a block size that's longer than the long axis
        block_size = min(block_size, (self.shape[0] if self._transpose else self.shape[1]))
        return block_size

    def join(self, other_source: 'MappedSource') -> 'MappedSource':
        if not isinstance(other_source, MappedSource):
            raise ValueError('Cannot append source type {} to this MappedSource'.format(type(other_source)))
        if set(self.aligned_arrays) != set(other_source.aligned_arrays):
            raise ValueError('Mismatch in aligned arrays')
        if self._transpose is not other_source._transpose:
            raise ValueError('Mismatch in transposed sources')
        # data buffer and aligned arrays can be a plain buffer or a binder, and transpose can be T/F
        keys = ['data_buffer'] + list(self.aligned_arrays)
        new_sources = OrderedDict()
        # For all sources (data_buffer plus others):
        # 1) promote this buffer to a binder if needed
        # 2) tell binder to extend to other source
        # This will raise error for various inconsistencies:
        # * dimension mismatch, concat axis mismatch, read/write mode mismatch
        # * missing aligned array keys
        # TODO: check consistency for electrode channels.. if the mapped channels are not consistent, it would be
        #  possible to join them after mirroring each source to direct maps with copy='all'
        if self._electrode_channels != other_source._electrode_channels:
            raise ValueError('Joining sources with different channel maps is not yet supported.')
        for k in keys:
            this_buf = getattr(self, k)
            if isinstance(this_buf, HDF5Buffer):
                this_buf = BufferBinder([this_buf], axis=0 if self._transpose else 1)
            new_buf = this_buf + getattr(other_source, k)
            new_sources[k] = new_buf
        data_buffer = new_sources.pop('data_buffer')
        binary_mask = self.binary_channel_mask & other_source.binary_channel_mask
        return MappedSource(data_buffer, electrode_field=self._electrode_field,
                            electrode_channels=self._electrode_channels, channel_mask=binary_mask,
                            aligned_arrays=new_sources, transpose=self._transpose,
                            raise_on_big_slice=self._raise_on_big_slice)

    def set_channel_mask(self, channel_mask):
        """
        Update the internal electrode channel mapping to exclude channel
        masked from an electrode vector binary mask.

        Parameters
        ----------
        channel_mask: binary ndarray or None
            A binary mask vector which is 1) the same length of the number of electrode channels
            or 2) the length of the number of currently active channels. If None, then the set
            of active channels resets to the original channels.

        """

        n_electrodes = len(self._electrode_channels)
        n_active_electrodes = len(self._active_channels)
        if channel_mask is None or not len(channel_mask):
            # (re)set to the full set of electrode channels
            self._active_channels = self._electrode_channels
        elif len(channel_mask) == n_electrodes:
            # reset the entire active electrode set
            self._active_channels = list()
            for n, c in enumerate(self._electrode_channels):
                if channel_mask[n]:
                    self._active_channels.append(c)
        elif len(channel_mask) == n_active_electrodes:
            # remove channels at the indices marked by the mask
            removed_idx = np.where(~channel_mask)[0]
            removed = [self._active_channels[i] for i in removed_idx]
            active_set = set(self._active_channels) - set(removed)
            self._active_channels = sorted(list(active_set))
        else:
            raise ValueError('channel_mask must be length {} (full map) '
                             'or {} (active set)'.format(n_electrodes, n_active_electrodes))

    def _slice_logic(self, slicer):
        """Translate the slicing object to point to the correct data channels on disk"""
        if not np.iterable(slicer):
            slicer = (slicer,)
        if self.is_direct_map:
            return slicer[::-1] if self._transpose else slicer
        chan_range = slicer[0]
        if len(slicer[1:]):
            time_range = slicer[1]
        else:
            # If the slice is for channels only and we're in map-subset mode, then simply return the slice without
            # modifying the range
            if self.channels_are_maps:
                return slicer
            time_range = slice(None)
        if isinstance(chan_range, (int, np.integer)):
            get_chan = self._active_channels[chan_range]
            return (time_range, get_chan) if self._transpose else (get_chan, time_range)
        # Map the sliced channels to corresponding active channels
        if isinstance(chan_range, slice):
            r_max = len(self._active_channels)
            chan_range = list(slice_to_range(chan_range, r_max))
        get_chans = [self._active_channels[c] for c in chan_range]
        return (time_range, get_chans) if self._transpose else (get_chans, time_range)

    def _check_slice_size(self, slicer):
        if not self._raise_on_big_slice:
            return
        shape = self.data_buffer.get_output_array(slicer, only_shape=True)
        size = np.prod(shape) * self.data_buffer.dtype.itemsize
        state = self._allow_big_slicing
        if size > load_params().memory_limit and not state:
            raise MemoryBlowOutError('A read with shape {} will be *very* large. Use the big_slices context to '
                                     'proceed'.format(shape))

    def slice_subset(self, slicer):
        new_electrode_channels = np.array(self._active_channels)[slicer].tolist()
        aligned = OrderedDict([(k, getattr(self, k)) for k in self.aligned_arrays])
        return MappedSource(self.data_buffer, electrode_field=self._electrode_field,
                            electrode_channels=new_electrode_channels,
                            aligned_arrays=aligned, transpose=self._transpose,
                            raise_on_big_slice=self._raise_on_big_slice)

    def cache_slice(self, slicer, sharedmem=False, **kwargs):
        slicer = self._slice_logic(slicer)
        self._check_slice_size(slicer)
        # TODO: not sure how plain memory will play with actual background reads?
        with self.data_buffer.transpose_reads(self._transpose), self.data_buffer.shared_memory(sharedmem):
            output = self.data_buffer.get_output_array(slicer)
        # Sub-process issues on Windows and unaccountable lag on MacOS, so
        # make the slice in this main process
        # print('can subprocess slice:', self._can_cache)
        if self._can_cache:
            p = BackgroundRead(self.data_buffer, slicer, transpose=self._transpose, output=output)
            p.start()
            self._caching_process = p
        else:
            # non-dispatched slicing
            slice_data_buffer(self.data_buffer, slicer, transpose=self._transpose, output=output)
        self._cache_output = output

    def get_cached_slice(self):
        # Remove sub-process slice retrieval for now
        if self._can_cache:
            if self._caching_process == None:
                print('There is no cached read pending')
                return
            self._caching_process.join()
            self._caching_process = None
        # try to release memory
        out = self._cache_output
        self._cache_output = None
        return out

        # if self._cache_output is None:
        #     print('There is no slice available')
        #     return
        # out = self._cache_output
        # self._cache_output = None
        # return out

    def __getitem__(self, slicer):
        """Return the sub-series of samples selected by slicer on (possibly a subset of) all channels"""
        # data goes out as [subset]channels x time
        slicer = self._slice_logic(slicer)
        if len(slicer) == 1 and self.channels_are_maps:
            return self.slice_subset(slicer)
        self._check_slice_size(slicer)
        return slice_data_buffer(self.data_buffer, slicer, self._transpose)

    def __setitem__(self, slicer, data):
        """Write the sub-series of samples selected by slicer (from possibly a subset of channels)"""
        # data comes in as [subset]channels x time
        if not self.writeable:
            print('Cannot write to this file')
            return
        # if both these conditions are try, something is probably wrong
        suspect_call = self.writeable and self._transpose
        # make sure subset mode is off
        with self.channels_are_maps(False):
            slicer = self._slice_logic(slicer)
        try:
            self.data_buffer[slicer] = data
        except IndexError as e:
            if suspect_call:
                tb = e.__traceback__
                msg = 'This mapped array is both writeable and in transpose mode, which probably is not what you mean '
                'to be doing and probably is the source of this error.'
                raise Exception(msg).with_traceback(tb)
            else:
                raise e

    def iter_channels(self, chans_per_block=None, use_max_memory=True, return_slice=False):
        # TODO: argument that supplies a channel order permutation
        # Just change the signature to use maximum memory by default
        return super(MappedSource, self).iter_channels(chans_per_block=chans_per_block,
                                                       use_max_memory=use_max_memory, return_slice=return_slice)

    def filter_array(self, **kwargs):
        from ecogdata.filt.time import filter_array
        kwargs['block_filter'] = bfilter
        # kwargs['inplace'] = True
        filter_array(self, **kwargs)
        if kwargs.get('out', None) is not None:
            return kwargs['out']
        return self

    def notch_filter(self, *args, **kwargs):
        from ecogdata.filt.time import notch_all
        kwargs['block_filter'] = bfilter
        # kwargs['inplace'] = True
        notch_all(self, *args, **kwargs)
        if kwargs.get('out', None) is not None:
            return kwargs['out']
        return self

    def mirror(self, new_rate_ratio=None, writeable=True, mapped=True, channel_compatible=False, filename='',
               copy='', new_sources=dict(), **map_args):
        # TODO: channel order permutation in mirrored source
        """
        Create an empty ElectrodeDataSource based on the current source, possibly with a new sampling rate and new
        access permissions.

        Parameters
        ----------
        new_rate_ratio: int or None
            Ratio of old to new sample rate for the mirrored array (> 1).
        writeable: bool
            Make any new MappedSource arrays writeable. This implies 1) datatype casting to floats, and 2) there is
            no more units conversion on the primary array.
        mapped: bool
            If False, mirror to a PlainArraySource (in memory). Else mirror into a new MappedSource.
        channel_compatible: bool
            If True, preserve the same number of raw data channels in a MappedSource. Otherwise, reduce the channels
            to just the set of active channels. Currently not supported for non-mapped mirrors.
        filename: str
            Name of the new MappedSource. If empty, use a self-destroying temporary file.
        copy: str
            Code whether to copy any arrays, which is only valid when new_rate_ratio is None or 1. 'aligned' copies
            aligned arrays. 'electrode' copies electrode data: only valid if channel_compatible is False.
            'all' copies all arrays. By default, nothing is copied.
        new_sources: dict
            If mapped=False, then pre-allocated arrays can be provided for each source (i.e. 'data_buffer' and any
            aligned arrays).
        map_args: dict
            Any other MappedSource arguments

        Returns
        -------
        datasource: ElectrodeDataSource subtype

        """

        T = self.shape[1]
        if new_rate_ratio:
            T = calc_new_samples(T, new_rate_ratio)

        # if casting to floating point, check for preferred precision
        fp_precision = load_params().floating_point.lower()
        fp_dtype = 'f' if fp_precision == 'single' else 'd'
        fp_dtype = np.dtype(fp_dtype)

        # unpack copy mode
        copy_electrodes_coded = copy.lower() in ('all', 'electrode')
        copy_aligned_coded = copy.lower() in ('all', 'aligned')
        diff_rate = T != self.shape[1]
        if copy_electrodes_coded:
            if diff_rate or channel_compatible:
                copy_electrodes = False
                print('Not copying electrode channels. Diff rate '
                      '({}) or indirect channel map ({})'.format(diff_rate, channel_compatible))
            else:
                copy_electrodes = True
        else:
            copy_electrodes = False
        if copy_aligned_coded:
            if diff_rate:
                copy_aligned = False
                print('Not copying aligned arrays: different sample rate')
            else:
                copy_aligned = True
        else:
            copy_aligned = False

        if mapped:
            if channel_compatible:
                electrode_channels = self._electrode_channels
                C = self.data_buffer.shape[1] if self._transpose else self.data_buffer.shape[0]
                channel_mask = self.binary_channel_mask
            else:
                C = self.shape[0]
                electrode_channels = None
                channel_mask = None
            if writeable:
                new_dtype = fp_dtype
                reopen_mode = 'r+'
                units_scale = None
            else:
                new_dtype = self.data_buffer.map_dtype
                reopen_mode = 'r'
                units_scale = self.data_buffer.units_scale
            tempfile = not filename
            if tempfile:
                with TempFilePool(mode='ab') as f:
                    # punt on the unlink-on-close issue for now with "delete=False"
                    # f.file.close()
                    filename = f
            # open as string just in case its a TempFilePool
            # TODO: (need to figure out how to make it work seamless as a string)
            with h5py.File(str(filename), 'w', libver='latest') as fw:
                # Create all new datasets as non-transposed
                fw.create_dataset(self._electrode_field, shape=(C, T), dtype=new_dtype, chunks=True)
                if copy_electrodes:
                    for block, sl in self.iter_blocks(return_slice=True, sharedmem=False):
                        fw[self._electrode_field][sl] = block
                for name in self.aligned_arrays:
                    arr = getattr(self, name)
                    if len(arr.shape) > 1:
                        dims = (arr.shape[1], T) if self._transpose else (arr.shape[0], T)
                    # is this correct ???
                    dtype = fp_dtype if writeable else arr.dtype
                    fw.create_dataset(name, shape=dims, dtype=dtype, chunks=True)
                    if copy_aligned:
                        aligned = getattr(self, name)[:]
                        fw[name][:] = aligned.T if self._transpose else aligned
                # set single write, multiple read mode AFTER datasets are created
                # Skip for now -- SWMR issue pending
                # fw.swmr_mode = True
            print(str(filename), 'reopen mode', reopen_mode)
            f_mapped = h5py.File(str(filename), reopen_mode, libver='latest')
            # Skip for now -- SWMR issue pending
            # if writeable:
            #     f_mapped.swmr_mode = True
            if isinstance(filename, TempFilePool):
                filename.register_to_close(f_mapped)
            # If the original buffer's unit scaling is 1.0 then we should keep real slices on
            real_slices = units_scale == 1.0
            return MappedSource.from_hdf_sources(f_mapped, self._electrode_field, units_scale=units_scale,
                                                 real_slices=real_slices, aligned_arrays=self.aligned_arrays,
                                                 transpose=False, electrode_channels=electrode_channels,
                                                 channel_mask=channel_mask, **map_args)
            # return MappedSource(f_mapped, self._electrode_field, electrode_channels=electrode_channels,
            #                     channel_mask=channel_mask, aligned_arrays=self.aligned_arrays,
            #                     transpose=False, **map_args)
        else:
            if channel_compatible:
                print('RAM mirrors are not channel compatible--use mapped=True')
            self._check_slice_size(np.s_[:, :T])
            C = self.shape[0]
            if 'data_buffer' in new_sources:
                new_source = new_sources['data_buffer']
            else:
                new_source = shm.shared_ndarray((C, T), fp_dtype.char)
            if copy_electrodes:
                for block, sl in self.iter_blocks(return_slice=True, sharedmem=False):
                    new_source[sl] = block
            # Kind of tricky with aligned fields -- assume that transpose means the same thing for them?
            # But also un-transpose them on this mirroring step
            aligned_arrays = dict()
            for name in self.aligned_arrays:
                arr = getattr(self, name)
                if len(arr.shape) > 1:
                    dims = (arr.shape[1], T) if self._transpose else (arr.shape[0], T)
                else:
                    dims = (T,)
                if name in new_sources:
                    aligned_arrays[name] = new_sources[name]
                else:
                    aligned_arrays[name] = shm.shared_ndarray(dims, fp_dtype.char)
                if copy_aligned:
                    aligned = getattr(self, name)[:]
                    aligned_arrays[name][:] = aligned.T if self._transpose else aligned
            return PlainArraySource(new_source, **aligned_arrays)

def bfilter(b, a, x, out=None, filtfilt=False, verbose=False, **extra):
    """
    Apply linear filter inplace over array x streaming from disk.

    Parameters
    ----------
    b: ndarray
        Transfer function numerator polynomial coefficients.
    a: ndarray
        Transfer function denominator polynomial coefficients.
    x: MappedSource
        Mapped data source
    out: MappedSource
        Mapped source to store filter output, otherwise store in x.
    filtfilt: bool
        If True, make a second filtering in reverse time sequence to create zero phase.
    verbose: bool
        Make a progress bar.

    """
    try:
        zii = lfilter_zi(b, a)
    except LinAlgError:
        # the integrating filter doesn't have valid zi
        zii = np.array([0.0])

    zi_sl = [np.newaxis] * 2
    zi_sl[1] = slice(None)
    zi_sl = tuple(zi_sl)
    xc_sl = [slice(None)] * 2
    xc_sl[1] = slice(0, 1)
    xc_sl = tuple(xc_sl)
    # fir_size = len(b)
    zi = None
    itr = x.iter_blocks(return_slice=True, sharedmem=False)
    # Disable verbose & tqdm -- garbage collection issue
    # verbose = False
    # if verbose:
    #     itr = tqdm(itr, desc='Blockwise filtering', leave=True, total=len(itr))
    for xc, sl in itr:
        if zi is None:
            zi = zii[zi_sl] * xc[xc_sl]
        # treat xc as mutable, since it is sliced from a mapped source
        # zi = lfilter_void(b, a, xc, xc, zi, axis=1)
        xc, zi = lfilter(b, a, xc, zi=zi, axis=1)
        if out is None:
            x[sl] = xc
        else:
            out[sl] = xc

    if not filtfilt:
        del xc
        return

    # Nullify initial conditions for reverse run (not sure if it's correct, but it's compatible with bfilter from
    # the ecogdata.filter.time.blocked_filter module
    zi = None
    if out is None:
        itr = x.iter_blocks(return_slice=True, reverse=True, sharedmem=False)
    else:
        itr = out.iter_blocks(return_slice=True, reverse=True, sharedmem=False)
    # if verbose:
    #     itr = tqdm(itr, desc='Blockwise filtering (reverse)', leave=True, total=len(itr))
    for xc, sl in itr:
        if zi is None:
            zi = zii[zi_sl] * xc[xc_sl]
        # zi = lfilter_void(b, a, xc, xc, zi, axis=1)
        xc, zi = lfilter(b, a, xc, zi=zi, axis=1)
        # write out with negative step slices (buffer will correct the write order)
        if out is None:
            x[sl] = xc
        else:
            out[sl] = xc
    del xc

def downsample_and_load(mapped_source, downsample_ratio, **kwargs):
    """

    Parameters
    ----------
    mapped_source: MappedSource
    downsample_ratio: int

    Returns
    -------
    new_source: PlainArraySource

    """
    kwargs.setdefault('filter_inplace', True)
    kwargs.setdefault('aggregate_aligned', True)
    new_source = mapped_source.mirror(new_rate_ratio=downsample_ratio, mapped=False)
    mapped_source.batch_change_rate(downsample_ratio, new_source, **kwargs)
    return new_source

max_stars_repo_name: miketrumpis/ecogdata

max_issues_repo_name: miketrumpis/ecogdata

max_forks_repo_name: miketrumpis/ecogdata

reason: import numpy,from numpy,from scipy


text:
import numpy as np
from PIL import Image
from skimage.io import imread, imsave
from skimage.transform import resize as imresize

def combine_and_save(source_files, output_file, arrangement='horizontal'):
    images = [Image.open(x) for x in source_files]
    widths, heights = zip(*(i.size for i in images))
    total_width = sum(widths)
    max_height = max(heights)
    if arrangement == 'vertical':
        new_im = Image.new('RGB', (max_height, total_width))
    else:
        new_im = Image.new('RGB', (total_width, max_height))
    x_offset = 0
    for im in images:
        if arrangement == 'vertical':
            new_im.paste(im, (0, x_offset))
        else:
            new_im.paste(im, (x_offset, 0))
        x_offset += im.size[0]
    new_im.save(output_file)

def resize_and_save(source_file, output_file, target_size=(64, 64), grayscale=False):
    img = imread(source_file, as_gray=grayscale)
    img = imresize(img, target_size, anti_aliasing=True)
    img = 255 * img
    img = img.astype(np.uint8)
    imsave(output_file, img)

max_stars_repo_name: dalikewara/pangsitpy

max_issues_repo_name: dalikewara/pangsitpy

max_forks_repo_name: dalikewara/pangsitpy