text:
import numpy as np
import math
from random import sample
import os
import pandas as pd

from scripts.utils.fit_poly import fit_poly
from scripts.utils.plotting import plot_profiles, plot_along_strike_profile

pd.options.mode.chained_assignment = None

def load_data(filename, dist_between_profiles):
    """
    Loads csv file as a DataFrame
    :param filename: <string> The filename to be loaded
    :param dist_between_profiles: <int> The distance along strike between profiles
    :return: <DataFrame> The DataFrame
    """

    df = pd.read_csv(f"../data/{filename}")
    df['dist_along_fault'] = (df['profile'] - 1) * dist_between_profiles
    return df

def add_row(df, ls):
    """
    Given a DataFrame and a list, append the list as a new row to the DataFrame.
    :param df: <DataFrame> The original DataFrame
    :param ls: <list> The new row to be added
    :return: <DataFrame> The DataFrame with the newly appended row
    """

    num_el = len(ls)
    new_row = pd.DataFrame(np.array(ls).reshape(1, num_el), columns = list(df.columns))
    df = df.append(new_row, ignore_index = True)
    return df

def save_to_csv(df):
    """
    Saves DataFrame to csv file.
    :param df <pandas.DataFrame> A pandas DataFrame to save
    :return: pandas.DataFrame.to_csv
    """

    filename = 'manual.csv'
    save_path = os.path.join(os.path.dirname(os.path.abspath(__file__)), '..', 'outputs', filename)
    return df.to_csv(save_path, index = False)

def manual_analysis(df, num):
    """
    Perform manual analysis on profiles and save height, width and slope features to file
    :param df: <DataFrame> The original DataFrame
    :param num: <int> The number of profiles to perform manual analysis on
    """

    samples_profiles = sample(df['profile'].to_list(), num)
    samples_profiles.sort()

    df_out = pd.DataFrame(columns = ['profile', 'dist_along_fault', 'height', 'width', 'slope'])

    for n, profile_number in enumerate(samples_profiles):

        sample_profile = df[df['profile'] == profile_number]

        dist_along_fault = int(sample_profile['dist_along_fault'].values[0])

        # Calculates the slope and derivative of slope based on the elevation profiles
        dzdx = list(np.insert(np.diff(sample_profile.loc[:, 'elevation']), 0, 0))
        sample_profile['dzdx'] = dzdx

        d2zdx2 = list(np.insert(np.diff(sample_profile.loc[:, 'dzdx']), 0, 0))
        sample_profile['d2zdx2'] = d2zdx2

        base, crest = plot_profiles(sample_profile)

        rmse_upper, upper, upper_poly_height, p_upper = fit_poly(sample_profile, crest[0], base[0], where = "upper")
        rmse_lower, lower, lower_poly_height, p_lower = fit_poly(sample_profile, crest[0], base[0], where = "lower")
        rmse_scarp, scarp, scarp_poly_height, p_scarp = fit_poly(sample_profile, crest[0], base[0], where = "scarp")

        # Finds the point of maximum slope on the scarp, return scarp heights from
        # upper and lower slopes
        ind = np.nanargmin(scarp['elevation'])
        h = upper_poly_height[ind] - lower_poly_height[ind]
        w = base[0] - crest[0]
        slope = math.degrees(math.atan(p_scarp[0]))
        if h > 100:  # stops very large scarps being picked due to error
            h = None

        df_out = add_row(df_out, [profile_number, dist_along_fault, h, w, slope])
        save_to_csv(df_out)

    plot_along_strike_profile(df_out)

if __name__ == '__main__':
    df = load_data(filename = "sample.csv", dist_between_profiles = 100)
    manual_analysis(df, num = 10)

ext: py

max_stars_repo_name: mshodge/sparta

max_issues_repo_name: mshodge/sparta

max_forks_repo_name: mshodge/sparta

reason: import numpy


text:
import resonances
import numpy as np
import pytest

def test_periodogram():
    t = np.linspace(0, 10, 101)
    y = np.sin(np.pi * t)
    (frequency, power) = resonances.libration.periodogram(t, y, minimum_frequency=0.0, maximum_frequency=2.0)

    key = 0
    for i, elem in enumerate(frequency):
        if elem == 0.5:  # the frequency (1 period in 2s => 1/2 = 0.5 Hz)
            key = i

    assert 1.0 == pytest.approx(power[key], 0.1)

def test_find_peaks_with_position():
    t = np.linspace(0, 10, 101)
    y = np.sin(np.pi * t) + np.sin(2 * np.pi * t)  # 0.5 Hz and 1 Hz
    (frequency, power) = resonances.libration.periodogram(t, y, minimum_frequency=0.0, maximum_frequency=2.0)

    peaks = resonances.libration.find_peaks_with_position(frequency, power)

    keys = []
    for i, elem in enumerate(frequency):
        if (elem == 0.5) or (elem == 1.0):  # 0.5 and 1 Hz
            keys.append(i)

    # find_peaks returns times, not frequencies
    assert True == (peaks['position'][0][0] <= 1.0 / 1.0 <= peaks['position'][0][1])
    assert True == (peaks['position'][1][0] <= 1.0 / 0.5 <= peaks['position'][1][1])

    for key in keys:
        assert key in peaks['peaks']

def test_filter():
    t = np.linspace(0, 10, 101)
    y = np.sin(np.pi * t) + np.sin(2 * np.pi * t)  # 0.5 Hz and 1 Hz
    (frequency, power) = resonances.libration.periodogram(t, y, minimum_frequency=0.0, maximum_frequency=2.0)
    y_filtered = resonances.libration.butter_lowpass_filter(y, 0.6, 1, 2, 5)
    frequency, power = resonances.libration.periodogram(t, y_filtered, minimum_frequency=0.0, maximum_frequency=2.0)
    peaks = resonances.libration.find_peaks_with_position(frequency, power, height=0.2)

    assert True == (peaks['position'][0][0] <= 1.0 / 0.5 <= peaks['position'][0][1])
    assert 1 == len(peaks['position'])

max_stars_repo_name: apetrov/resonances

max_issues_repo_name: apetrov/resonances

max_forks_repo_name: apetrov/resonances


text:
from __future__ import division

import cv2
import numpy as np

def resize(img, ignore_mask, poses, shape):
    """resize img, mask and annotations"""
    height, width, _ = img.shape

    resized_img = cv2.resize(img, shape)
    ignore_mask = cv2.resize(ignore_mask.astype(
        np.uint8), shape).astype('bool')
    poses[:, :, :2] = (poses[:, :, :2] * np.array(shape) /
                       np.array((width, height)))
    return resized_img, ignore_mask, poses

max_stars_repo_name: iory/chainer-openpose

max_issues_repo_name: iory/chainer-openpose

max_forks_repo_name: iory/chainer-openpose


text:
import numpy as np
k = np.arange(1,7)
print(k)
k = k.reshape((2,3))
print(k)
b = k[np.newaxis,:]
k = k[:,np.newaxis]
print(k)
print(b)

max_stars_repo_path: Numpy practice/Numpy7.py

max_stars_repo_name: Macorov/Numpy-Python-Practice

max_issues_repo_path: Numpy practice/Numpy7.py

max_issues_repo_name: Macorov/Numpy-Python-Practice

max_forks_repo_path: Numpy practice/Numpy7.py

max_forks_repo_name: Macorov/Numpy-Python-Practice


text:
from numpy import loadtxt, where, zeros, e, array, log, ones, mean, where
from pylab import scatter, show, legend, xlabel, ylabel, plot
from scipy.optimize import fmin_bfgs

def sigmoid(X):
    '''Compute the sigmoid function '''
    #d = zeros(shape=(X.shape))
    den = 1.0 + e ** (-1.0 * X)
    d = 1.0 / den
    return d

def compute_cost(theta, X, y):
    '''
    Comput cost for logistic regression
    '''
    #Number of training samples

    theta.shape = (1, 3)

    m = y.size

    h = sigmoid(X.dot(theta.T))

    J = (1.0 / m) * ((-y.T.dot(log(h))) - ((1.0 - y.T).dot(log(1.0 - h))))

    return - 1 * J.sum()

def compute_grad(theta, X, y):

    #print theta.shape

    theta.shape = (1, 3)

    grad = zeros(3)

    h = sigmoid(X.dot(theta.T))

    delta = h - y

    l = grad.size

    for i in range(l):
        sumdelta = delta.T.dot(X[:, i])
        grad[i] = (1.0 / m) * sumdelta * - 1

    theta.shape = (3,)

    return  grad

#load the dataset
data = loadtxt('ex2data1.txt', delimiter=',')

X = data[:, 0:2]
y = data[:, 2]

pos = where(y == 1)
neg = where(y == 0)
scatter(X[pos, 0], X[pos, 1], marker='o', c='b')
scatter(X[neg, 0], X[neg, 1], marker='x', c='r')
xlabel('Exam 1 score')
ylabel('Exam 2 score')
legend(['Not Admitted', 'Admitted'])
#show()

m, n = X.shape

y.shape = (m, 1)

#Add intercept term to x and X_test
it = ones(shape=(m, 3))
it[:, 1:3] = X

def decorated_cost(it, y):
    def f(theta):
        return compute_cost(theta, it, y)

    def fprime(theta):
        return compute_grad(theta, it, y)

    #Initialize theta parameters
    theta = zeros(3)

    return fmin_bfgs(f, theta, fprime, disp=True, maxiter=400)

decorated_cost(it, y)
theta = [-25.161272, 0.206233, 0.201470]

#Plotting the decision boundary
plot_x = array([min(it[:, 1]) - 2, max(it[:, 2]) + 2])
plot_y = (- 1.0 / theta[2]) * (theta[1] * plot_x + theta[0])
plot(plot_x, plot_y)
legend(['Decision Boundary', 'Not admitted', 'Admitted'])
#show()

prob = sigmoid(array([1.0, 45.0, 85.0]).dot(array(theta).T))
print 'For a student with scores 45 and 85, we predict and admission ' + \
    'probability of %f' % prob

def predict(theta, X):
    '''Predict whether the label
    is 0 or 1 using learned logistic
    regression parameters '''
    m, n = X.shape
    p = zeros(shape=(m, 1))

    h = sigmoid(X.dot(theta.T))

    for it in range(0, h.shape[0]):
        if h[it] > 0.5:
            p[it, 0] = 1
        else:
            p[it, 0] = 0

    return p

#Compute accuracy on our training set
p = predict(array(theta), it)
print 'Train Accuracy: %f' % ((y[where(p == y)].size / float(y.size)) * 100.0)

reason: from numpy,from scipy


text:
import os
import sys
import torch
import torchvision
import torch.nn as nn
from torch.autograd import Variable

import time
import yaml
import argparse
import numpy as np
from printer import Printer
from dataset import get_data_loader
from model import Model
import datetime
import copy
from util import make_dir, get_optimizer, AverageMeter, save_train_info, norm_flow
from gyro import torch_QuaternionProduct, torch_QuaternionReciprocal, torch_norm_quat

os.environ["CUDA_VISIBLE_DEVICES"] = "0"

def run_epoch(model, loader, cf, epoch, lr, optimizer=None, is_training=True, USE_CUDA=True, clip_norm=0):
    no_flo = False
    number_virtual, number_real = cf['data']["number_virtual"], cf['data']["number_real"]
    avg_loss = AverageMeter()
    if is_training:
        model.net.train()
        model.unet.train()
    else:
        model.net.eval()
        model.unet.eval()
    for i, data in enumerate(loader, 0):
        # get the inputs; data is a list of [inputs, labels]
        real_inputs, times, flo, flo_back, real_projections, real_postion, ois, real_queue_idx = data
        print("Fininsh Load data")

        real_inputs = real_inputs.type(torch.float) #[b,60,84=21*4]
        real_projections = real_projections.type(torch.float)
        flo = flo.type(torch.float)
        flo_back = flo_back.type(torch.float)
        ois = ois.type(torch.float)

        batch_size, step, dim = real_inputs.size()
        times = times.numpy()
        real_queue_idx = real_queue_idx.numpy()
        virtual_queue = loader.dataset.random_init_virtual_queue(batch_size, real_postion[:,0,:].numpy(), times[:,1]) # TODO
        # virtual_queue = [None] * batch_size
        loss = 0
        model.net.init_hidden(batch_size)
        for j in range(step):
            virtual_inputs, vt_1 = loader.dataset.get_virtual_data(
                virtual_queue, real_queue_idx, times[:, j], times[:, j+1], times[:, 0], batch_size, number_virtual, real_postion[:,j])

            real_inputs_step = real_inputs[:,j,:]
            inputs = torch.cat((real_inputs_step,virtual_inputs), dim = 1)

            # inputs = Variable(real_inputs_step)
            if USE_CUDA:
                real_inputs_step = real_inputs_step.cuda()
                virtual_inputs = virtual_inputs.cuda()
                inputs = inputs.cuda()
                if no_flo is False:
                    flo_step = flo[:,j].cuda()
                    flo_back_step = flo_back[:,j].cuda()
                else:
                    flo_step = None
                    flo_back_step = None
                vt_1 = vt_1.cuda()
                real_projections_t = real_projections[:,j+1].cuda()
                real_projections_t_1 = real_projections[:,j].cuda()
                real_postion_anchor = real_postion[:,j].cuda()
                ois_step = ois[:,j].cuda()

            if no_flo is False:
                b, h, w, _ = flo_step.size()
                flo_step = norm_flow(flo_step, h, w)
                flo_back_step = norm_flow(flo_back_step, h, w)

            if is_training:
                if no_flo is False:
                    flo_out = model.unet(flo_step, flo_back_step)
                else:
                    flo_out = None

                if j < 1:
                    for i in range(2):
                        out = model.net(inputs, flo_out, ois_step)
                else:
                    out = model.net(inputs, flo_out, ois_step)
            else:
                with torch.no_grad():
                    if no_flo is False:
                        flo_out = model.unet(flo_step, flo_back_step)
                    else:
                        flo_out = None

                    if j < 1:
                        for i in range(2):
                            out = model.net(inputs, flo_out, ois_step)
                    else:
                        out = model.net(inputs, flo_out, ois_step)

            if epoch <= 30:
                follow = True
            else:
                follow = False

            if epoch > 30:
                undefine = True
            else:
                undefine = False

            if epoch > 40:
                optical = True
            else:
                optical = False

            loss_step = model.loss(out, vt_1, virtual_inputs, real_inputs_step, \
                flo_step, flo_back_step, real_projections_t, real_projections_t_1, real_postion_anchor, \
                follow = follow, undefine = undefine, optical = optical, stay = optical)

            loss = loss_step

            virtual_position = virtual_inputs[:, -4:]
            pos = torch_QuaternionProduct(virtual_position, real_postion_anchor)
            out = torch_QuaternionProduct(out, pos)

            if USE_CUDA:
                out = out.cpu().detach().numpy()

            virtual_queue = loader.dataset.update_virtual_queue(batch_size, virtual_queue, out, times[:,j+1])

            if (j+1) % 10 == 0:
                print("Step: "+str(j+1)+"/"+str(step))
                print(loss)
            loss = torch.sum(loss)
            if is_training:
                optimizer.zero_grad()
                loss.backward(retain_graph=True)
                if clip_norm:
                    nn.utils.clip_grad_norm_(model.net.parameters(), max_norm=clip_norm)
                    nn.utils.clip_grad_norm_(model.unet.parameters(), max_norm=clip_norm)
                optimizer.step()

            avg_loss.update(loss.item(), batch_size)

    return avg_loss.avg

def train(args = None):
    torch.autograd.set_detect_anomaly(True)
    config_file = args.config
    cf = yaml.load(open(config_file, 'r'), Loader=yaml.Loader)

    USE_CUDA = cf['data']["use_cuda"]
    seed = cf['train']["seed"]

    torch.manual_seed(seed)
    if USE_CUDA:
        torch.cuda.manual_seed(seed)

    checkpoints_dir = cf['data']['checkpoints_dir']
    epochs = cf["train"]["epoch"]
    snapshot = cf["train"]["snapshot"]
    decay_epoch = cf['train']['decay_epoch']
    init_lr = cf["train"]["init_lr"]
    lr_decay = cf["train"]["lr_decay"]
    lr_step = cf["train"]["lr_step"]
    clip_norm = cf["train"]["clip_norm"]
    load_model = cf["model"]["load_model"]

    checkpoints_dir = make_dir(checkpoints_dir, cf)

    if load_model is None:
        log_file = open(os.path.join(cf["data"]["log"], cf['data']['exp']+'.log'), 'w+')
    else:
        log_file = open(os.path.join(cf["data"]["log"], cf['data']['exp']+'.log'), 'a')
    printer = Printer(sys.stdout, log_file).open()

    print('----Print Arguments Setting------')
    for key in cf:
        print('{}:'.format(key))
        for para in cf[key]:
            print('{:50}:{}'.format(para,cf[key][para]))
        print('\n')

    # Define the model
    model = Model(cf)
    optimizer = get_optimizer(cf["train"]["optimizer"], model, init_lr, cf)

    for idx, m in enumerate(model.net.children()):
        print('{}:{}'.format(idx,m))
    for idx, m in enumerate(model.unet.children()):
        print('{}:{}'.format(idx,m))

    if load_model is not None:
        print("------Load Pretrined Model--------")
        checkpoint = torch.load(load_model)
        model.net.load_state_dict(checkpoint['state_dict'])
        model.unet.load_state_dict(checkpoint['unet'])
        print("------Resume Training Process-----")
        optimizer.load_state_dict(checkpoint['optim_dict'])
        epoch_load = checkpoint['epoch']
        print("Epoch load: ", epoch_load)
    else:
        epoch_load = 0

    if USE_CUDA:
        model.net.cuda()
        model.unet.cuda()
        if load_model is not None:
            for state in optimizer.state.values():
                for k, v in state.items():
                    if isinstance(v, torch.Tensor):
                        state[k] = v.cuda()
            for param in optimizer.param_groups:
                init_lr = param['lr']

    print("-----------Load Dataset----------")
    train_loader, test_loader = get_data_loader(cf, no_flo = False)

    print("----------Start Training----------")
    currentDT = datetime.datetime.now()
    print(currentDT.strftime(" %Y-%m-%d %H:%M:%S"))

    start_time = time.time()

    if lr_step:
        decay_epoch = list(range(1+lr_step, epochs+1, lr_step))

    lr = init_lr

    for count in range(epoch_load+1, epochs+1):
        if decay_epoch != None and count in decay_epoch:
            lr *= lr_decay
            for param in optimizer.param_groups:
                param['lr'] *= lr_decay

        print("Epoch: %d, learning_rate: %.5f" % (count, lr))

        train_loss = run_epoch(model, train_loader, cf, count, lr, optimizer=optimizer, clip_norm=clip_norm, is_training=True, USE_CUDA=USE_CUDA)

        test_loss = run_epoch(model, test_loader, cf, count, lr, is_training=False, USE_CUDA=USE_CUDA)

        time_used = (time.time() - start_time) / 60
        print("Epoch %d done | TrLoss: %.4f | TestLoss: %.4f | Time_used: %.4f minutes" % (
            count, train_loss,  test_loss, time_used))

        if count % snapshot == 0:
            save_train_info("epoch", checkpoints_dir, cf, model, count, optimizer)
            save_train_info("last", checkpoints_dir, cf, model, count, optimizer)
            print("Model stored at epoch %d"%count)

    currentDT = datetime.datetime.now()
    print(currentDT.strftime(" %Y-%m-%d %H:%M:%S"))
    print("------------End Training----------")
    return

if __name__ == '__main__':
    parser = argparse.ArgumentParser("Training model")
    parser.add_argument("--config", default="./conf/stabilzation_train.yaml", help="Config file.")
    args = parser.parse_args()
    train(args = args)


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

from PIL import Image
from absl import app
from options import FLAGS as opts
import data
import datetime
import functools
import glob
import losses
import networks
import numpy as np
import options
import os.path as osp
import random
import skimage.measure
import staged_model
import tensorflow as tf
import time
import utils
import os

def build_model_fn(use_exponential_moving_average=True):
  """Builds and returns the model function for an estimator.

  Args:
    use_exponential_moving_average: bool. If true, the exponential moving
    average will be used.

  Returns:
    function, the model_fn function typically required by an estimator.
  """
  arch_type = opts.arch_type
  use_appearance = opts.use_appearance
  def model_fn(features, labels, mode, params):
    """An estimator build_fn."""
    del labels, params
    if mode == tf.compat.v1.estimator.ModeKeys.TRAIN:
      step = tf.compat.v1.train.get_global_step()

      x_in = features['conditional_input']
      x_gt = features['expected_output']  # ground truth output
      x_app = features['peek_input']

      if opts.training_pipeline == 'staged':
        ops = staged_model.create_computation_graph(x_in, x_gt, x_app=x_app,
                                                    arch_type=opts.arch_type)
        op_increment_step = tf.compat.v1.assign_add(step, 1)
        train_disc_op = ops['train_disc_op']
        train_renderer_op = ops['train_renderer_op']
        train_op = tf.compat.v1.group(train_disc_op, train_renderer_op, op_increment_step)

        utils.HookReport.log_tensor(ops['total_loss_d'], 'total_loss_d')
        utils.HookReport.log_tensor(ops['loss_d_real'], 'loss_d_real')
        utils.HookReport.log_tensor(ops['loss_d_fake'], 'loss_d_fake')
        utils.HookReport.log_tensor(ops['total_loss_g'], 'total_loss_g')
        utils.HookReport.log_tensor(ops['loss_g_gan'], 'loss_g_gan')
        utils.HookReport.log_tensor(ops['loss_g_recon'], 'loss_g_recon')
        utils.HookReport.log_tensor(step, 'global_step')

        return tf.compat.v1.estimator.EstimatorSpec(
            mode=mode, loss=ops['total_loss_d'] + ops['total_loss_g'],
            train_op=train_op)
      else:
        raise NotImplementedError('%s training is not implemented.' %
                                  opts.training_pipeline)
    elif mode == tf.compat.v1.estimator.ModeKeys.EVAL:
      raise NotImplementedError('Eval is not implemented.')
    else:  # all below modes are for difference inference tasks.
      # Build network and initialize inference variables.
      g_func = networks.RenderingModel(arch_type, use_appearance)
      if use_appearance:
        app_func = g_func.get_appearance_encoder()
      if use_exponential_moving_average:
        ema = tf.compat.v1.train.ExponentialMovingAverage(decay=0.999)
        var_dict = ema.variables_to_restore()
        tf.compat.v1.train.init_from_checkpoint(osp.join(opts.train_dir), var_dict)

      if mode == tf.compat.v1.estimator.ModeKeys.PREDICT:
        x_in = features['conditional_input']
        if use_appearance:
          x_app = features['peek_input']
          x_app_embedding, _, _ = app_func(x_app)
        else:
          x_app_embedding = None
        y = g_func(x_in, x_app_embedding)
        tf.compat.v1.logging.info('DBG: shape of y during prediction %s.' % str(y.shape))
        return tf.compat.v1.estimator.EstimatorSpec(mode=mode, predictions=y)

      # 'eval_subset' mode is same as PREDICT but it concatenates the output to
      # the input render, semantic map and ground truth for easy comparison.
      elif mode == 'eval_subset':
        x_in = features['conditional_input']
        x_gt = features['expected_output']
        if use_appearance:
          x_app = features['peek_input']
          x_app_embedding, _, _ = app_func(x_app)
        else:
          x_app_embedding = None
        y = g_func(x_in, x_app_embedding)
        tf.compat.v1.logging.info('DBG: shape of y during prediction %s.' % str(y.shape))
        x_in_rgb = tf.compat.v1.slice(x_in, [0, 0, 0, 0], [-1, -1, -1, 3])
        if opts.use_semantic:
          x_in_semantic = tf.compat.v1.slice(x_in, [0, 0, 0, 4], [-1, -1, -1, 3])
          output_tuple = tf.compat.v1.concat([x_in_rgb, x_in_semantic, y, x_gt], axis=2)
        else:
          output_tuple = tf.compat.v1.concat([x_in_rgb, y, x_gt], axis=2)
        return tf.compat.v1.estimator.EstimatorSpec(mode=mode, predictions=output_tuple)

      # 'compute_appearance' mode computes and returns the latent z vector.
      elif mode == 'compute_appearance':
        assert use_appearance, 'use_appearance is set to False!'
        x_app_in = features['peek_input']
        # NOTE the following line is a temporary hack (which is
        # specially bad for inputs smaller than 512x512).
        x_app_in = tf.image.resize_image_with_crop_or_pad(x_app_in, 512, 512)
        app_embedding, _, _ = app_func(x_app_in)
        return tf.estimator.EstimatorSpec(mode=mode, predictions=app_embedding)

      # 'interpolate_appearance' mode expects an already computed latent z
      # vector as input passed a value to the dict key 'appearance_embedding'.
      elif mode == 'interpolate_appearance':
        assert use_appearance, 'use_appearance is set to False!'
        x_in = features['conditional_input']
        x_app_embedding = features['appearance_embedding']
        y = g_func(x_in, x_app_embedding)
        tf.logging.info('DBG: shape of y during prediction %s.' % str(y.shape))
        return tf.estimator.EstimatorSpec(mode=mode, predictions=y)
      else:
        raise ValueError('Unsupported mode: ' + mode)

  return model_fn

def make_sample_grid_and_save(est, dataset_name, dataset_parent_dir, grid_dims,
                              output_dir, cur_nimg):
  """Evaluate a fixed set of validation images and save output.

  Args:
    est: tf,estimator.Estimator, TF estimator to run the predictions.
    dataset_name: basename for the validation tfrecord from which to load
      validation images.
    dataset_parent_dir: path to a directory containing the validation tfrecord.
    grid_dims: 2-tuple int for the grid size (1 unit = 1 image).
    output_dir: string, where to save image samples.
    cur_nimg: int, current number of images seen by training.

  Returns:
    None.
  """
  num_examples = grid_dims[0] * grid_dims[1]
  def input_val_fn():
    dict_inp = data.provide_data(
        dataset_name=dataset_name, parent_dir=dataset_parent_dir, subset='val',
        batch_size=1, crop_flag=True, crop_size=opts.train_resolution,
        seeds=[0], max_examples=num_examples,
        use_appearance=opts.use_appearance, shuffle=0)
    x_in = dict_inp['conditional_input']
    x_gt = dict_inp['expected_output']  # ground truth output
    x_app = dict_inp['peek_input']
    return x_in, x_gt, x_app

  def est_input_val_fn():
    x_in, _, x_app = input_val_fn()
    features = {'conditional_input': x_in, 'peek_input': x_app}
    return features

  images = [x for x in est.predict(est_input_val_fn)]
  images = np.array(images, 'f')
  images = images.reshape(grid_dims + images.shape[1:])
  utils.save_images(utils.to_png(utils.images_to_grid(images)), output_dir,
                    cur_nimg)

def visualize_image_sequence(est, dataset_name, dataset_parent_dir,
                             input_sequence_name, app_base_path, output_dir):
  """Generates an image sequence as a video and stores it to disk."""
  batch_sz = opts.batch_size
  def input_seq_fn():
    dict_inp = data.provide_data(
        dataset_name=dataset_name, parent_dir=dataset_parent_dir,
        subset=input_sequence_name, batch_size=batch_sz, crop_flag=False,
        seeds=None, use_appearance=False, shuffle=0)
    x_in = dict_inp['conditional_input']
    return x_in

  # Compute appearance embedding only once and use it for all input frames.
  app_rgb_path = app_base_path + '_reference.png'
  app_rendered_path = app_base_path + '_color.png'
  app_depth_path = app_base_path + '_depth.png'
  app_sem_path = app_base_path + '_seg_rgb.png'
  x_app = _load_and_concatenate_image_channels(
      app_rgb_path, app_rendered_path, app_depth_path, app_sem_path)
  def seq_with_single_appearance_inp_fn():
    """input frames with a fixed latent appearance vector."""
    x_in_op = input_seq_fn()
    x_app_op = tf.convert_to_tensor(x_app)
    x_app_tiled_op = tf.tile(x_app_op, [tf.shape(x_in_op)[0], 1, 1, 1])
    return {'conditional_input': x_in_op,
            'peek_input': x_app_tiled_op}

  images = [x for x in est.predict(seq_with_single_appearance_inp_fn)]
  for i, gen_img in enumerate(images):
    output_file_path = osp.join(output_dir, 'out_%04d.png' % i)
    print('Saving frame #%d to %s' % (i, output_file_path))
    with tf.gfile.Open(output_file_path, 'wb') as f:
      f.write(utils.to_png(gen_img))

def train(dataset_name, dataset_parent_dir, load_pretrained_app_encoder,
          load_trained_fixed_app, save_samples_kimg=50):
  """Main training procedure.

  The trained model is saved in opts.train_dir, the function itself does not
   return anything.

  Args:
    save_samples_kimg: int, period (in KiB) to save sample images.

  Returns:
    None.
  """
  image_dir = osp.join(opts.train_dir, 'images')  # to save validation images.
  tf.compat.v1.gfile.MakeDirs(image_dir)
  session_config = tf.ConfigProto()
  session_config.gpu_options.allow_growth = True
  config = tf.compat.v1.estimator.RunConfig(
      save_summary_steps=(1 << 10) // opts.batch_size,
      save_checkpoints_steps=(save_samples_kimg << 10) // opts.batch_size,
      keep_checkpoint_max=5,
      log_step_count_steps=1 << 30,
      session_config=session_config)
  model_dir = opts.train_dir
  if (opts.use_appearance and load_trained_fixed_app and
      not tf.compat.v1.train.latest_checkpoint(model_dir)):
    tf.compat.v1.logging.warning('***** Loading resume_step from %s!' %
                       opts.fixed_appearance_train_dir)
    resume_step = utils.load_global_step_from_checkpoint_dir(
        opts.fixed_appearance_train_dir)
  else:
    tf.compat.v1.logging.warning('***** Loading resume_step (if any) from %s!' %
                       model_dir)
    resume_step = utils.load_global_step_from_checkpoint_dir(model_dir)
  if resume_step != 0:
    tf.compat.v1.logging.warning('****** Resuming training at %d!' % resume_step)

  model_fn = build_model_fn()  # model function for TFEstimator.

  hooks = [utils.HookReport(1 << 12, opts.batch_size)]

  if opts.use_appearance and load_pretrained_app_encoder:
    tf.compat.v1.logging.warning('***** Will warm-start from %s!' %
                       opts.appearance_pretrain_dir)
    ws = tf.compat.v1.estimator.WarmStartSettings(
        ckpt_to_initialize_from=opts.appearance_pretrain_dir,
        vars_to_warm_start='appearance_net/.*')
  elif opts.use_appearance and load_trained_fixed_app:
    tf.compat.v1.logging.warning('****** Finetuning will warm-start from %s!' %
                       opts.fixed_appearance_train_dir)
    ws = tf.compat.v1.estimator.WarmStartSettings(
        ckpt_to_initialize_from=opts.fixed_appearance_train_dir,
        vars_to_warm_start='.*')
  else:
    ws = None
    tf.compat.v1.logging.warning('****** No warm-starting; using random initialization!')

  est = tf.compat.v1.estimator.Estimator(model_fn, model_dir, config, params={},
                               warm_start_from=ws)

  for next_kimg in range(opts.save_samples_kimg, opts.total_kimg + 1,
                         opts.save_samples_kimg):
    next_step = (next_kimg << 10) // opts.batch_size
    if opts.num_crops == -1:  # use random crops
      crop_seeds = None
    else:
      crop_seeds = list(100 * np.arange(opts.num_crops))
    input_train_fn = functools.partial(
        data.provide_data, dataset_name=dataset_name,
        parent_dir=dataset_parent_dir, subset='train',
        batch_size=opts.batch_size, crop_flag=True,
        crop_size=opts.train_resolution, seeds=crop_seeds,
        use_appearance=opts.use_appearance)
    est.train(input_train_fn, max_steps=next_step, hooks=hooks)
    tf.compat.v1.logging.info('DBG: kimg=%d, cur_step=%d' % (next_kimg, next_step))
    tf.compat.v1.logging.info('DBG: Saving a validation grid image %06d to %s' % (
        next_kimg, image_dir))
    make_sample_grid_and_save(est, dataset_name, dataset_parent_dir, (3, 3),
                              image_dir, next_kimg << 10)

def _build_inference_estimator(model_dir):
  model_fn = build_model_fn()
  est = tf.estimator.Estimator(model_fn, model_dir)
  return est

def evaluate_sequence(dataset_name, dataset_parent_dir, virtual_seq_name,
                      app_base_path):
  output_dir = osp.join(opts.train_dir, 'seq_output_%s' % virtual_seq_name)
  tf.gfile.MakeDirs(output_dir)
  est = _build_inference_estimator(opts.train_dir)
  visualize_image_sequence(est, dataset_name, dataset_parent_dir,
                           virtual_seq_name, app_base_path, output_dir)

def evaluate_image_set(dataset_name, dataset_parent_dir, subset_suffix,
                       output_dir=None, batch_size=6):
  if output_dir is None:
    output_dir = osp.join(opts.train_dir, 'validation_output_%s' % subset_suffix)
  tf.compat.v1.gfile.MakeDirs(output_dir)
  model_fn_old = build_model_fn()
  def model_fn_wrapper(features, labels, mode, params):
    del mode
    return model_fn_old(features, labels, 'eval_subset', params)
  model_dir = opts.train_dir
  est = tf.compat.v1.estimator.Estimator(model_fn_wrapper, model_dir)
  est_inp_fn = functools.partial(
      data.provide_data, dataset_name=dataset_name,
      parent_dir=dataset_parent_dir, subset=subset_suffix,
      batch_size=batch_size, use_appearance=opts.use_appearance, shuffle=0)

  print('Evaluating images for subset %s' % subset_suffix)
  images = [x for x in est.predict(est_inp_fn)]
  print('Evaluated %d images' % len(images))
  for i, img in enumerate(images):
    output_file_path = osp.join(output_dir, 'out_%04d.png' % i)
    print('Saving file #%d: %s' % (i, output_file_path))
    with tf.compat.v1.gfile.Open(output_file_path, 'wb') as f:
      f.write(utils.to_png(img))

def _load_and_concatenate_image_channels(rgb_path=None, rendered_path=None,
                                         depth_path=None, seg_path=None,
                                         size_multiple=64):
  """Prepares a single input for the network."""
  if (rgb_path is None and rendered_path is None and depth_path is None and
      seg_path is None):
    raise ValueError('At least one of the inputs has to be not None')

  channels = ()
  if rgb_path is not None:
    rgb_img = np.array(Image.open(rgb_path)).astype(np.float32)
    rgb_img = utils.crop_to_multiple(rgb_img, size_multiple)
    channels = channels + (rgb_img,)
  if rendered_path is not None:
    rendered_img = np.array(Image.open(rendered_path)).astype(np.float32)
    if not opts.use_alpha:
      rendered_img = rendered_img[:, :, :3]  # drop the alpha channel
    rendered_img = utils.crop_to_multiple(rendered_img, size_multiple)
    channels = channels + (rendered_img,)
  if depth_path is not None:
    depth_img = np.array(Image.open(depth_path))
    depth_img = depth_img.astype(np.float32)
    depth_img = utils.crop_to_multiple(depth_img[:, :, np.newaxis],
                                       size_multiple)
    channels = channels + (depth_img,)
    # depth_img = depth_img * (2.0 / 255) - 1.0
  if seg_path is not None:
    seg_img = np.array(Image.open(seg_path)).astype(np.float32)
    seg_img = utils.crop_to_multiple(seg_img, size_multiple)
    channels = channels + (seg_img,)
  # Concatenate and normalize channels
  img = np.dstack(channels)
  img = np.expand_dims(img, axis=0)
  img = img * (2.0 / 255) - 1.0
  return img

def infer_dir(model_dir, input_dir, output_dir):
  tf.gfile.MakeDirs(output_dir)
  est = _build_inference_estimator(opts.train_dir)

  def read_image(base_path, is_appearance=False):
    if is_appearance:
      ref_img_path = base_path + '_reference.png'
    else:
      ref_img_path = None
    rendered_img_path = base_path + '_color.png'
    depth_img_path = base_path + '_depth.png'
    seg_img_path = base_path + '_seg_rgb.png'
    img = _load_and_concatenate_image_channels(
        rgb_path=ref_img_path, rendered_path=rendered_img_path,
        depth_path=depth_img_path, seg_path=seg_img_path)
    return img

  def get_inference_input_fn(base_path, app_base_path):
    x_in = read_image(base_path, False)
    x_app_in = read_image(app_base_path, True)
    def infer_input_fn():
      return {'conditional_input': x_in, 'peek_input': x_app_in}
    return infer_input_fn

  file_paths = sorted(glob.glob(osp.join(input_dir, '*_depth.png')))
  base_paths = [x[:-10] for x in file_paths]  # remove the '_depth.png' suffix
  for inp_base_path in base_paths:
    est_inp_fn = get_inference_input_fn(inp_base_path, inp_base_path)
    img = next(est.predict(est_inp_fn))
    basename = osp.basename(inp_base_path)
    output_img_path = osp.join(output_dir, basename + '_out.png')
    print('Saving generated image to %s' % output_img_path)
    with tf.gfile.Open(output_img_path, 'wb') as f:
      f.write(utils.to_png(img))

def joint_interpolation(model_dir, app_input_dir, st_app_basename,
                        end_app_basename, camera_path_dir):
  """
  Interpolates both viewpoint and appearance between two input images.
  """
  # Create output direcotry
  output_dir = osp.join(model_dir, 'joint_interpolation_out')
  tf.gfile.MakeDirs(output_dir)

  # Build estimator
  model_fn_old = build_model_fn()
  def model_fn_wrapper(features, labels, mode, params):
    del mode
    return model_fn_old(features, labels, 'interpolate_appearance', params)
  def appearance_model_fn(features, labels, mode, params):
    del mode
    return model_fn_old(features, labels, 'compute_appearance', params)
  config = tf.estimator.RunConfig(
      save_summary_steps=1000, save_checkpoints_steps=50000,
      keep_checkpoint_max=50, log_step_count_steps=1 << 30)
  model_dir = model_dir
  est = tf.estimator.Estimator(model_fn_wrapper, model_dir, config, params={})
  est_app = tf.estimator.Estimator(appearance_model_fn, model_dir, config,
                                   params={})

  # Compute appearance embeddings for the two input appearance images.
  app_inputs = []
  for app_basename in [st_app_basename, end_app_basename]:
    app_rgb_path = osp.join(app_input_dir, app_basename + '_reference.png')
    app_rendered_path = osp.join(app_input_dir, app_basename + '_color.png')
    app_depth_path = osp.join(app_input_dir, app_basename + '_depth.png')
    app_seg_path = osp.join(app_input_dir, app_basename + '_seg_rgb.png')
    app_in = _load_and_concatenate_image_channels(
        rgb_path=app_rgb_path, rendered_path=app_rendered_path,
        depth_path=app_depth_path, seg_path=app_seg_path)
    # app_inputs.append(tf.convert_to_tensor(app_in))
    app_inputs.append(app_in)

  embedding1 = next(est_app.predict(
      lambda: {'peek_input': app_inputs[0]}))
  embedding1 = np.expand_dims(embedding1, axis=0)
  embedding2 = next(est_app.predict(
      lambda: {'peek_input': app_inputs[1]}))
  embedding2 = np.expand_dims(embedding2, axis=0)

  file_paths = sorted(glob.glob(osp.join(camera_path_dir, '*_depth.png')))
  base_paths = [x[:-10] for x in file_paths]  # remove the '_depth.png' suffix

  # Compute interpolated appearance embeddings
  num_interpolations = len(base_paths)
  interpolated_embeddings = []
  delta_vec = (embedding2 - embedding1) / (num_interpolations - 1)
  for delta_iter in range(num_interpolations):
    x_app_embedding = embedding1 + delta_iter * delta_vec
    interpolated_embeddings.append(x_app_embedding)

  # Generate and save interpolated images
  for frame_idx, embedding in enumerate(interpolated_embeddings):
    # Read in input frame
    frame_render_path = osp.join(base_paths[frame_idx] + '_color.png')
    frame_depth_path = osp.join(base_paths[frame_idx] + '_depth.png')
    frame_seg_path = osp.join(base_paths[frame_idx] + '_seg_rgb.png')
    x_in = _load_and_concatenate_image_channels(
        rgb_path=None, rendered_path=frame_render_path,
        depth_path=frame_depth_path, seg_path=frame_seg_path)

    img = next(est.predict(
        lambda: {'conditional_input': tf.convert_to_tensor(x_in),
                 'appearance_embedding': tf.convert_to_tensor(embedding)}))
    output_img_name = '%s_%s_%03d.png' % (st_app_basename, end_app_basename,
                                          frame_idx)
    output_img_path = osp.join(output_dir, output_img_name)
    print('Saving interpolated image to %s' % output_img_path)
    with tf.gfile.Open(output_img_path, 'wb') as f:
      f.write(utils.to_png(img))

def interpolate_appearance(model_dir, input_dir, target_img_basename,
                           appearance_img1_basename, appearance_img2_basename):
  # Create output direcotry
  output_dir = osp.join(model_dir, 'interpolate_appearance_out')
  tf.gfile.MakeDirs(output_dir)

  # Build estimator
  model_fn_old = build_model_fn()
  def model_fn_wrapper(features, labels, mode, params):
    del mode
    return model_fn_old(features, labels, 'interpolate_appearance', params)
  def appearance_model_fn(features, labels, mode, params):
    del mode
    return model_fn_old(features, labels, 'compute_appearance', params)
  config = tf.estimator.RunConfig(
      save_summary_steps=1000, save_checkpoints_steps=50000,
      keep_checkpoint_max=50, log_step_count_steps=1 << 30)
  model_dir = model_dir
  est = tf.estimator.Estimator(model_fn_wrapper, model_dir, config, params={})
  est_app = tf.estimator.Estimator(appearance_model_fn, model_dir, config,
                                   params={})

  # Compute appearance embeddings for the two input appearance images.
  app_inputs = []
  for app_basename in [appearance_img1_basename, appearance_img2_basename]:
    app_rgb_path = osp.join(input_dir, app_basename + '_reference.png')
    app_rendered_path = osp.join(input_dir, app_basename + '_color.png')
    app_depth_path = osp.join(input_dir, app_basename + '_depth.png')
    app_seg_path = osp.join(input_dir, app_basename + '_seg_rgb.png')
    app_in = _load_and_concatenate_image_channels(
        rgb_path=app_rgb_path, rendered_path=app_rendered_path,
        depth_path=app_depth_path, seg_path=app_seg_path)
    # app_inputs.append(tf.convert_to_tensor(app_in))
    app_inputs.append(app_in)

  embedding1 = next(est_app.predict(
      lambda: {'peek_input': app_inputs[0]}))
  embedding2 = next(est_app.predict(
      lambda: {'peek_input': app_inputs[1]}))
  embedding1 = np.expand_dims(embedding1, axis=0)
  embedding2 = np.expand_dims(embedding2, axis=0)

  # Compute interpolated appearance embeddings
  num_interpolations = 10
  interpolated_embeddings = []
  delta_vec = (embedding2 - embedding1) / num_interpolations
  for delta_iter in range(num_interpolations + 1):
    x_app_embedding = embedding1 + delta_iter * delta_vec
    interpolated_embeddings.append(x_app_embedding)

  # Read in the generator input for the target image to render
  rendered_img_path = osp.join(input_dir, target_img_basename + '_color.png')
  depth_img_path = osp.join(input_dir, target_img_basename + '_depth.png')
  seg_img_path = osp.join(input_dir, target_img_basename + '_seg_rgb.png')
  x_in = _load_and_concatenate_image_channels(
      rgb_path=None, rendered_path=rendered_img_path,
      depth_path=depth_img_path, seg_path=seg_img_path)

  # Generate and save interpolated images
  for interpolate_iter, embedding in enumerate(interpolated_embeddings):
    img = next(est.predict(
        lambda: {'conditional_input': tf.convert_to_tensor(x_in),
                 'appearance_embedding': tf.convert_to_tensor(embedding)}))
    output_img_name = 'interpolate_%s_%s_%s_%03d.png' % (
        target_img_basename, appearance_img1_basename, appearance_img2_basename,
        interpolate_iter)
    output_img_path = osp.join(output_dir, output_img_name)
    print('Saving interpolated image to %s' % output_img_path)
    with tf.gfile.Open(output_img_path, 'wb') as f:
      f.write(utils.to_png(img))

def main(argv):
  del argv
  configs_str = options.list_options()
  tf.compat.v1.gfile.MakeDirs(opts.train_dir)
  with tf.compat.v1.gfile.Open(osp.join(opts.train_dir, 'configs.txt'), 'wb') as f:
    f.write(configs_str)
  tf.compat.v1.logging.info('Local configs\n%s' % configs_str)

  if opts.run_mode == 'train':
    dataset_name = opts.dataset_name
    dataset_parent_dir = opts.dataset_parent_dir
    load_pretrained_app_encoder = opts.load_pretrained_app_encoder
    load_trained_fixed_app = opts.load_from_another_ckpt
    batch_size = opts.batch_size
    train(dataset_name, dataset_parent_dir, load_pretrained_app_encoder,
          load_trained_fixed_app)
  elif opts.run_mode == 'eval':  # generate a camera path output sequence from TFRecord inputs.
    dataset_name = opts.dataset_name
    dataset_parent_dir = opts.dataset_parent_dir
    virtual_seq_name = opts.virtual_seq_name
    inp_app_img_base_path = opts.inp_app_img_base_path
    evaluate_sequence(dataset_name, dataset_parent_dir, virtual_seq_name,
                      inp_app_img_base_path)
  elif opts.run_mode == 'eval_subset':  # generate output for validation set (encoded as TFRecords)
    dataset_name = opts.dataset_name
    dataset_parent_dir = opts.dataset_parent_dir
    virtual_seq_name = opts.virtual_seq_name
    evaluate_image_set(dataset_name, dataset_parent_dir, virtual_seq_name,
                       opts.output_validation_dir, opts.batch_size)
  elif opts.run_mode == 'eval_dir':  # evaluate output for a directory with input images
    input_dir = opts.inference_input_path
    output_dir = opts.inference_output_dir
    model_dir = opts.train_dir
    infer_dir(model_dir, input_dir, output_dir)
  elif opts.run_mode == 'interpolate_appearance':  # interpolate appearance only between two images.
    model_dir = opts.train_dir
    input_dir = opts.inference_input_path
    target_img_basename = opts.target_img_basename
    app_img1_basename = opts.appearance_img1_basename
    app_img2_basename = opts.appearance_img2_basename
    interpolate_appearance(model_dir, input_dir, target_img_basename,
                           app_img1_basename, app_img2_basename)
  elif opts.run_mode == 'joint_interpolation':  # interpolate viewpoint and appearance between two images
    model_dir = opts.train_dir
    app_input_dir = opts.inference_input_path
    st_app_basename = opts.appearance_img1_basename
    end_app_basename = opts.appearance_img2_basename
    frames_dir = opts.frames_dir
    joint_interpolation(model_dir, app_input_dir, st_app_basename,
                        end_app_basename, frames_dir)
  else:
    raise ValueError('Unsupported --run_mode %s' % opts.run_mode)

if __name__ == '__main__':
  os.environ["TF_MIN_GPU_MULTIPROCESSOR_COUNT"]="2"
  os.environ["CUDA_VISIBLE_DEVICES"]="0,1"
  app.run(main)


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

"""
Handle annotated whole-slide images (WSI; also known as virtual slides).

"""

import warnings
from distutils.version import StrictVersion
import random

import numpy as np
import openslide
import cv2

import annotation as reader
import tissue

class _AnnotatedOpenSlide(openslide.OpenSlide):
    """An open *annotated* WSI.

    Wrapper of `openslide.OpenSlide` to handle annotation polygons together
    with their corresponding OpenSlide Python WSI instances.

    Parameters
    ----------
    filename: str
        Path to slide file.
    annotation_filename: str
        Path to XML annotation file.
    xml_style: str
        The annotation XML style. Typically associated with a computational
        histology challenge releasing the dataset.
        * 'asap': CAMELYON grand challenges in pathology.
        * 'bach': BACH Grand Challenge on Breast Cancer Histology Images.

    Attributes
    ----------
    filename: str
        The WSI file name to read.
    annotation_filename: str
        Corresponding XML annotation file.
    polygons: list of lists of 2-tuples
        Polygon region annotations.
    labels: list of str
        Polygon region annotation labels.
    label_map: dict {str: int}
        Correspondence between polygon region annotation labels and integer
        values used in numpy masks.

    """

    def __init__(self, filename, annotation_filename, xml_style):
        openslide.OpenSlide.__init__(self, filename)
        self.filename = filename  # Useful to name predicted annotations
        self.annotation_filename = annotation_filename
        self.xml_style = xml_style
        self.polygons = None
        self.label_map = None

        # Label annotations as:
        #   1: Normal tissue
        #   2: Tumor tissue
        # 0 is automatically assigned to all unannotated pixels
        if self.annotation_filename is not None:
            if self.xml_style == 'asap':
                self.polygons, self.labels = reader.asap_annotations(
                    self.annotation_filename)
                # CAMELYON17 data
                if 'metastases' in self.labels:
                    # Avoid value 0 (used by default for unlabeled regions)
                    self.label_map = {'metastases': 2, 'normal': 1}
                # CAMELYON16 data
                elif '_0' in self.labels or '_1' in self.labels:
                    self.label_map = {'_0': 2, '_1': 2, '_2': 1}
                elif 'Tumor' in self.labels or 'Exclusion' in self.labels:
                    self.label_map = {'Tumor': 2, 'Exclusion': 1}
                else:  # Predicted annotations
                    self.label_map = {'predicted_tumor': 1}
            elif self.xml_style == 'bach':
                # Value 1 is reserved for 'normal' tissue annotations
                self.label_map = {'Benign': 2, 'Carcinoma in situ': 3,
                                  'Invasive carcinoma': 4}
                self.polygons, self.labels = reader.bach_annotations(
                    self.annotation_filename)
            else:
                raise ValueError(
                    '"xml_style" value must be either "asap" or "bach".')
        else:
            if self.xml_style is not None:
                warnings.warn('"xml_style" is only used if an ' +
                              '"annotation_filename" is provided.')

class Slide(_AnnotatedOpenSlide):
    """An open WSI, with or without annotations.

    Added functionality to OpenSlide WSI instances, optionally with
    annotations.

    Parameters
    ----------
    filename: str
        Path to slide file. Namesake argument to `openslide.OpenSlide` and
        `AnnotatedOpenSlide`.
    annotation_filename: str
        Path to XML annotation file. Namesake argument to `AnnotatedOpenSlide`.
    xml_style: str
        The computational histology challenge releasing the dataset. Namesake
        argument to `AnnotatedOpenSlide`.

    Attributes
    ----------
    tissue_mask: Numpy 2D array
        Tissue RoI ``annotation``.
    tissue_label_map: dict {str: int}
        Correspondence between polygon region annotation labels and integer
        values used in numpy masks.
    downsampling_factor: float
        The scaling factor relative to level 0 dimensions.
    downsampled_slide: PIL image
        The scaled down slide generated when running the tissue detector.

    Examples
    --------
    >>> from wsipre import slide
    >>> wsi = slide.Slide('tumor_001.tif', 'tumor_001.xml', 'asap')
    >>> wsi.label_map
    {'_0': 2, '_1': 2, '_2': 1}
    >>> wsi.level_count
    10
    >>> len(wsi.polygons)
    2

    >>> thumbnail, mask, dwnspl_factor = wsi.get_thumbnail_with_annotation(
    ...     size=(3000, 3000), polygon_type='line', line_thickness=5)

    >>> slide_region = wsi.read_region_with_annotation(
    ...     location=(65000, 110000), level=4, size=(1000, 1000))

    """

    def __init__(self, filename, annotation_filename=None, xml_style=None):
        self.tissue_mask = None
        self.tissue_label_map = None
        self.downsampling_factor = None
        self.downsampled_slide = None

        # OpenSlide version < 3.4.1 lacks support for some recent file formats
        version = openslide.__library_version__
        if StrictVersion(version) < StrictVersion('3.4.1'):
            warnings.warn(
                f'OpenSlide version {version} lacks support for some ' +
                'whole-slide image file formats. It is highly recommended ' +
                'to update to a more recent version.')

        if annotation_filename is None:
            if xml_style is not None:
                warnings.warn(
                    '"xml_style" is not used (no annotation was provided).')
            openslide.OpenSlide.__init__(self, filename)
        else:
            _AnnotatedOpenSlide.__init__(
                self, filename, annotation_filename, xml_style)

    def _draw_polygons(self, mask, polygons, polygon_type,
                       line_thickness=None):
        """Convert polygon vertex coordinates to polygon drawing."""
        for poly, label in zip(polygons, self.labels):
            if polygon_type == 'line':
                mask = cv2.polylines(mask, [poly], True,
                                     int(self.label_map[label]),
                                     line_thickness)
            elif polygon_type == 'area':
                if line_thickness is not None:
                    warnings.warn('"line_thickness" is only used if ' +
                                  '"polygon_type" is "line".')

                mask = cv2.fillPoly(mask, [poly], int(self.label_map[label]))
            else:
                raise ValueError(
                    'Accepted "polygon_type" values are "line" or "area".')

        return mask

    def _draw_tissue_polygons(self, mask, polygons, polygon_type,
                              line_thickness=None):
        """Convert tissue polygon vertex coordinates to polygon drawing."""
        tissue_label = 1

        for poly in polygons:
            if polygon_type == 'line':
                mask = cv2.polylines(
                    mask, [poly], True, tissue_label, line_thickness)
            elif polygon_type == 'area':
                if line_thickness is not None:
                    warnings.warn('"line_thickness" is only used if ' +
                                  '"polygon_type" is "line".')

                mask = cv2.fillPoly(mask, [poly], tissue_label)
            else:
                raise ValueError(
                    'Accepted "polygon_type" values are "line" or "area".')

        return mask

    def get_thumbnail_with_annotation(self, size, polygon_type='area',
                                      line_thickness=None):
        """Convert *annotated* WSI to a scaled-down thumbnail.

        Parameters
        ----------
        size: 2-tuple
            Maximum size of compressed image thumbnail as (width, height).
        polygon_type: str
            Type of polygon drawing:
            * 'line'
            * 'area'
        line_thickness: int
            Polygon edge line thickness. Required if ``polygon_type`` is
            ``line``.

        Returns
        -------
        Scaled down slide (as PIL RGBA image) and annotation mask (as numpy 2D
        array), plus the scaling factor.

        """
        if self.polygons is None:
                raise ValueError(
                    'No annotation is available. Please load an annotation ' +
                    'or consider using "get_thumbnail" instead.')

        thumb = self.get_thumbnail(size)

        width_factor = self.dimensions[0] / size[0]
        height_factor = self.dimensions[1] / size[1]
        downsampling_factor = max(width_factor, height_factor)

        polygons = np.array([np.array(poly) for poly in self.polygons])
        polygons = polygons / downsampling_factor
        polygons = [poly.astype('int32') for poly in polygons]

        mask = np.zeros(thumb.size[::-1])  # Reverse width and height
        mask = self._draw_polygons(
            mask, polygons, polygon_type, line_thickness)

        return thumb, mask, downsampling_factor

    def read_region_with_annotation(self, location, level, size,
                                    polygon_type='area', line_thickness=None):
        """Crop a smaller region from an *annotated* WSI.

        Get a defined region from a WSI, together with its annotation mask.

        Parameters
        ----------
        location: 2-tuple
            X and y coordinates of the top left pixel. Namesake argument to
            OpenSlide's "read_region" function.
        level: int
            Slide level. Namesake argument to OpenSlide's "read_region"
            function.
        size: 2-tuple
            Crop size (width, height). Namesake argument to OpenSlide's
            "read_region" function.
        polygon_type: str
            Type of polygon drawing:
            * 'line'
            * 'area'
        line_thickness: int
            Polygon edge line thickness. Required if 'polygon_type' is 'line'.

        Returns
        -------
        Image region (as PIL RGBA image) and the annotation mask (as numpy 2d
        array).

        """
        if self.polygons is None:
                raise ValueError(
                    'No annotation is available. Please load an annotation ' +
                    'or consider using "read_region" instead.')
        downsampling_factor = self.level_downsamples[level]
        slide_region = self.read_region(location, level, size)

        polygons = np.array([np.array(poly) for poly in self.polygons])
        polygons = polygons / downsampling_factor
        polygons = [poly.astype('int32') for poly in polygons]

        location = tuple(int(coord / downsampling_factor)
                         for coord in location)

        # Convert polygon coordinates to patch coordinates
        polygons = [poly - np.array(location) for poly in polygons]

        mask_region = np.zeros(size[::-1])
        mask_region = self._draw_polygons(
           mask_region, polygons, polygon_type, line_thickness)

        return slide_region, mask_region

    def get_tissue_mask(self, downsampling_factor=64, polygon_type='area',
                        line_thickness=None):
        """Generate tissue region annotation.

        Make a binary mask annotating tissue regions of interest (RoI) on the
        WSI using an automatic threshold-based segmentation method inspired by
        the one used by Wang *et al.* [1]_. Briefly, the method consists on the
        following steps, starting from a WSI:
            * Select downsampling level (typically a factor of 64)
            * Transfer from the RGB to the HSV color space
            * Determine optimal threshold value in the saturation channel using
              the Otsu algorithm [2]_
            * Threshold image to generate a binary mask
            * Fill in small holes and remove small objects

        .. [1] Dayong Wang, Aditya Khosla, Rishab Gargeya, Humayun Irshad,
           Andrew H. Beck, "Deep Learning for Identifying Metastatic Breast
           Cancer", arXiv:1606.05718.

        .. [2] Nobuyuki Otsu, "A Threshold Selection Method from Gray-level
           Histograms", IEEE Trans Syst Man Cybern., 9(1):62–66, 1979.

        Parameters
        ----------
        downsampling_factor: int
            The desired factor to downsample the image by, since full WSIs will
            not fit in memory. The image's closest level downsample is found
            and used.
        polygon_type: str
            Type of polygon drawing:
            * 'line'
            * 'area'
        line_thickness: int
            Polygon edge line thickness. Required if ``polygon_type`` is
            ``line``.

        Returns
        -------
        Binary mask as numpy 2D array, RGB slide image (in the used
        downsampling level, to allow visualization) and downsampling factor.

        """
        (mask_contours,
         self.downsampled_slide,
         self.downsampling_factor) = tissue.detect_tissue(
             self, downsampling_factor)

        mask = np.zeros(self.downsampled_slide.shape[:2])
        self.tissue_mask = self._draw_tissue_polygons(
           mask, mask_contours, polygon_type, line_thickness)
        self.tissue_label_map = {'background': 0, 'tissue': 1}

        return self

    def _get_random_coordinates(self, pixels_in_roi, level,
                                mask_downsampling_factor, size):
        """Get region (patch) coordinates at random."""
        coordinates = random.choice(pixels_in_roi)

        # Scale coordinates up to level-0 dimensions
        coordinates = [int(x * mask_downsampling_factor) for x in coordinates]

        # OpenSlide Python's read_region takes top left corner position, which
        # effectively excludes tissue above and to the left of annotations
        # Fix it by offsetting coordinates: convert to center pixel position
        # (subtract half of patch size from x and y coordinates)
        half_size = [x * self.level_downsamples[level] / 2 for x in size]
        row_location = coordinates[0] - half_size[0]
        col_location = coordinates[1] - half_size[1]

        # Make sure coordinates are still within slide margins
        cols_edge = self.dimensions[0] - size[0]
        rows_edge = self.dimensions[1] - size[1]
        row_location = max(0, min(row_location, rows_edge))
        col_location = max(0, min(col_location, cols_edge))

        return int(col_location), int(row_location)  # OpenSlide: width, height

    def _get_pixels_avoiding_labels(self, avoid_labels):
        # Can only do this if an annotation is available
        if self.polygons is None:
            raise ValueError('No annotation is available.')

        # Generate annotated thumbnail with exact dimensions of tissue mask
        rows, cols = self.tissue_mask.shape
        _, mask, _ = self.get_thumbnail_with_annotation(
            size=(cols, rows), polygon_type='area')

        # Locate pixels in tissue and not in "avoid_labels"
        new_mask = np.where(self.tissue_mask == 1, 1, mask)
        new_mask = np.where(np.isin(mask, avoid_labels), 2, new_mask)

        return tuple(zip(*np.where(new_mask == 1)))

    def _max_repeated_pixel_ratio(self, image):
        """Compute ratio of count of most common pixel value to total count."""
        image = np.array(image)
        _, counts = np.unique(image, return_counts=True)

        return np.max(counts) / image.size

    def _get_area_ratio(self, mask, target_class):
        """Compute ratio of pixels labeled as 'target_class' to all pixels."""
        class_pixels = mask[np.where(mask == target_class)].size

        return class_pixels / mask.size

    def read_random_tissue_patch(self, level, size, avoid_labels=None):
        """Crop random patch from detected tissue on WSI.

        Parameters
        ----------
        level: int
            Slide level. Namesake argument to OpenSlide's "read_region"
            function.
        size: 2-tuple
            Crop size (width, height). Namesake argument to OpenSlide's
            "read_region" function.
        avoid_labels: iterable
            Labels from provided annotation to avoid when targeting tissue.
            An annotation is must be available.

        Returns
        -------
        Image region or patch as PIL RGBA image (plus annotation mask as 2D
        Numpy array if labels to avoid are provided).

        """
        if self.tissue_mask is None:
            self.get_tissue_mask()

        if avoid_labels is None:
            # Select coordinates matching tissue (labeled as 1)
            pixels_in_roi = tuple(zip(*np.where(self.tissue_mask == 1)))
        else:
            # Select coordinates matching tissue (labeled as 1) and avoiding
            # indicated labels
            # Can only do this if an annotation is available
            if self.polygons is None:
                raise ValueError('No annotation is available.')

            # Generate annotated thumbnail with exact dimensions of tissue mask
            rows, cols = self.tissue_mask.shape
            _, mask, _ = self.get_thumbnail_with_annotation(
                size=(cols, rows), polygon_type='area')

            # Locate pixels in tissue and not in "avoid_labels"
            new_mask = np.where(self.tissue_mask == 1, 1, mask)
            new_mask = np.where(np.isin(mask, avoid_labels), 2, new_mask)

            pixels_in_roi = tuple(zip(*np.where(new_mask == 1)))

        coordinates = self._get_random_coordinates(
            pixels_in_roi, level, self.downsampling_factor, size)
        patch = self.read_region(coordinates, level, size)

        # Avoid false positive regions where a large proportion of pixels is
        # exactly the same
        while self._max_repeated_pixel_ratio(patch) > 0.5:
            coordinates = self._get_random_coordinates(
                pixels_in_roi, level, self.downsampling_factor, size)
            patch = self.read_region(coordinates, level, size)

        # Also, make sure the patch is mostly normal tissue
        # (or backgound, not actually excluding it here...)
        if avoid_labels is not None:
            patch, mask = self.read_region_with_annotation(
                coordinates, level, size, polygon_type='area')
            avoid_ratios = [self._get_area_ratio(mask, x)
                            for x in avoid_labels]
            i = 0
            while any(x > 0.1 for x in avoid_ratios):
                i += 1
                coordinates = self._get_random_coordinates(
                    pixels_in_roi, level, self.downsampling_factor, size)
                patch, mask = self.read_region_with_annotation(
                    coordinates, level, size, polygon_type='area')
                avoid_ratios = [self._get_area_ratio(mask, x)
                                for x in avoid_labels]
                if i > 15:
                    raise ValueError(
                        'Cannot seem to find patch matching requirements.')

            return patch, mask

        return patch

    def _pick_random_coordinates(self, level, size, target_class):
        """Pick random patch to crop from WSI."""
        # Get thumbnail and select coordinates at random
        _, mask, downsampling_factor = self.get_thumbnail_with_annotation(
             size=(5000, 5000), polygon_type='area')

        # Select coordinates matching target class
        pixels_in_roi = tuple(zip(*np.where(mask == target_class)))

        coordinates = self._get_random_coordinates(
            pixels_in_roi, level, downsampling_factor, size)

        # Get ratio (polygon type must be 'area')
        _, mask = self.read_region_with_annotation(
            coordinates, level, size, polygon_type='area')
        area_ratio = self._get_area_ratio(mask, target_class)

        return coordinates, area_ratio

    def read_random_patch(self, level, size, target_class,
                          min_class_area_ratio, polygon_type='area',
                          line_thickness=None):
        """Crop random patch from WSI.

        Select random location within target class according to provided
        annotation.

        Parameters
        ----------
        level: int
            Slide level. Namesake argument to OpenSlide's "read_region"
            function.
        size: 2-tuple
            Crop size (width, height). Namesake argument to OpenSlide's
            "read_region" function.
        target_class: int
            The class annotation of the central pixel of the patch.
            function.
        min_class_area_ratio: float [0, 1]
            Minimum ratio of target class pixels to total pixels.
        polygon_type: str
            Type of polygon drawing:
            * 'line'
            * 'area'
        line_thickness: int
            Polygon edge line thickness. Required if ``polygon_type`` is
            ``line``.

        Returns
        -------
        Image region or patch as PIL RGBA image plus annotation mask as 2D
        Numpy array.

        """
        if not 0 <= min_class_area_ratio <= 1:
            raise ValueError(
                '"min_class_area_ratio" must be in the interval [0, 1].')

        coordinates, area_ratio = self._pick_random_coordinates(
            level, size, target_class)

        # When the slide level is high, the following while loop may be
        # infinite, since the tumor area will never be high enough.
        # Try to sample a few times and throw error if not successful
        i = 0
        while area_ratio < min_class_area_ratio:
            i += 1
            coordinates, area_ratio = self._pick_random_coordinates(
                level, size, target_class)
            if i > 15:
                raise ValueError(
                    'Cannot seem to find patch with a minimum of ' +
                    f'{min_class_area_ratio * 100}% of class {target_class}.' +
                    ' The chosen slide level may be too high.')

        # Get actual data (chosen polygon type)
        patch, mask = self.read_region_with_annotation(
            coordinates, level, size, polygon_type, line_thickness)

        return patch, mask

max_stars_repo_name: JoostvDoorn/wsipre

max_issues_repo_name: JoostvDoorn/wsipre

max_forks_repo_name: JoostvDoorn/wsipre


text:
import numpy as np
from numpy.linalg import eigvalsh, eigh, matrix_power
import matplotlib.pyplot as plt
from scipy.special import gamma
from scipy.optimize import curve_fit
from qca import multipage
from figure3 import select, moving_average

# plotting defaults
import matplotlib as mpl

mpl.rcParams["text.latex.preamble"] = [r"\usepackage{amsmath}"]
font = {"size": 9, "weight": "normal"}
mpl.rcParams["mathtext.fontset"] = "stix"
mpl.rcParams["font.family"] = "STIXGeneral"
mpl.rcParams["pdf.fonttype"] = 42
mpl.rc("font", **font)

def brody_fit(x, n):
    def brody_func(x, eta, A):
        b = (gamma((eta + 2) / (eta + 1))) ** (eta + 1.0)
        return A * b * (eta + 1.0) * x ** eta * np.exp(-b * x ** (eta + 1.0))

    popt, pcov = curve_fit(brody_func, x, n, p0=[0.0, 1.0], bounds=[0, 1])

    def func(x):
        return brody_func(x, *popt)

    return func, popt, pcov

from measures import renyi_entropy
L = 18
IC = "c1_f0"
Skey = [13, 14, 1, 6]
cs = ["darkturquoise", "darkorange", "limegreen", "crimson"]
for j, (c, S) in enumerate(zip(cs, Skey)):
    sim = select(L, S, IC, V="H", BC="0")
    h5file = sim["h5file"]
    d = h5file["cut_half"][:]
    for ti, rho in enumerate(d[100:101]):
        spec, vecs = eigh(rho)
        fig, axs = plt.subplots(1, 2, figsize=(6, 3))
        print(renyi_entropy(rho))
        axs[0].set_title("spectrum")
        axs[0].semilogy(spec, color=c, marker="o")

        axs[1].set_title("column vector magnitude")
        axs[1].imshow(np.abs(vecs), cmap="gist_gray_r")
        fig.suptitle("$T_{%d}$"%S)
multipage("figures/figure4/eigenRDM.pdf")
plt.close("all")
print("done")
if __name__ == "__main__":
    L = 18
    IC = "c1_f0"
    Skey = [13, 14, 1, 6]
    cs = ["darkturquoise", "darkorange", "limegreen", "crimson"]
    fig, axs = plt.subplots(1, 1, figsize=(4, 3), sharex=False)
    fig2, axs2 = plt.subplots(1, 1, figsize=(3, 2), sharex=True)
    fig2.subplots_adjust(left=0.2, bottom=0.2, hspace=0.1)

    fig3s = []

    for j, (c, S) in enumerate(zip(cs, Skey)):
        sim = select(L, S, IC, V="H", BC="0")
        h5file = sim["h5file"]
        try:
            espec = h5file["espec"]
        except:
            d = h5file["cut_half"]
            espec = np.zeros((d.shape[0], d.shape[1]))
            for ti, rho in enumerate(d):
                espec[ti, :] = eigvalsh(rho)
            h5file["espec"] = espec

        etas = []
        detas = []
        svns = []
        ii = 0
        t0 = 10

        fig3, axs3 = plt.subplots(3, 3, figsize=(4, 4), sharex=True, sharey=True)
        for ti, es in enumerate(espec[t0:1000]):
            es = es[es > 1e-6]
            NN = len(es)
            es = np.sort(es)
            es = es[NN // 3 : 2 * NN // 3]
            ns = range(len(es))
            s = es[1:] - es[:-1]
            s /= np.mean(s)
            n, bin, _ = axs.hist(
                s, density=True, alpha=1, histtype="step", bins=10, log=False
            )

            x = (bin[1:] + bin[:-1]) / 2.0
            xs = np.linspace(x[0], x[-1], 100)
            xs = xs[xs > 0]

            func, popt, pcov = brody_fit(x, n)
            detas.append(np.sqrt(np.diag(pcov)[0]))
            etas.append(popt[0])

            if (ti+t0) % 100 == 0:
                row, col = ii // 3, ii % 3
                ax3 = axs3[row, col]
                dx = x[1] - x[0]
                n = np.insert(n, 0, 0)
                n = np.insert(n, len(n), 0)
                x = np.insert(x, 0, x[0] - dx / 2)
                x = np.insert(x, len(x), x[-1] + dx / 2)
                ax3.step(x, n, where="mid")
                ax3.plot(xs, func(xs))
                ax3.set_title(f"t={t0+ti}", pad=-13)
                fig3.suptitle("$R = %d$" % S)
                ii += 1
                if col == 1 and row == 2:
                    ax3.set_xlabel("$\delta E/\overline{\delta E}$")
                if col == 0 and row == 1:
                    ax3.set_ylabel("density")
                ax3.tick_params(direction="inout")

        fig3.subplots_adjust(hspace=0, wspace=0)
        fig3s.append(fig3)
        etas = np.array(etas)
        detas = np.array(detas)
        ts = np.arange(2, len(etas) + 2)
        mask = detas < 1
        etas = etas[mask]
        detas = detas[mask]
        ts = ts[mask]

        if S == 6:
            pass
        else:
            if S == 13:
                label = r"$R = %s$" % S
            else:

                label = str(S)
            aetas = moving_average(etas, n=L)
            axs2.plot(aetas, marker=None, color=c, label=label, lw=1)
            avgerr = np.mean(detas)
            axs2.plot(ts, etas, c=c, alpha=0.3)
            #axs2.errorbar(ts, etas, yerr=detas, color=c)
        axs2.set_xticks([0, 250, 500, 750, 1000])

    axs2.legend(loc="lower right")
    axs.set_xlabel("$\delta E / \overline{\delta E}$")
    axs2.set_xlabel("$t$")
    axs2.set_ylabel("$\eta$")
    fig.tight_layout()
    fig2.tight_layout()

    multipage(
        "figures/figure4/spectrum_statistics_fixed-10bins.pdf",
        figs=[fig2] + fig3s,
        clip=True,
        dpi=10 * fig.dpi,
    )

max_stars_repo_name: lhillber/qops

max_issues_repo_name: lhillber/qops

max_forks_repo_name: lhillber/qops

reason: import numpy,from numpy,from scipy


text:
import csv
import numpy as np
import os
import re
import scipy.stats
import sys
import time

# Given a regular expression, list the files that match it, and ask for user input
def selectFile(regex, subdirs = False):
	files = []
	if subdirs:
		for (dirpath, dirnames, filenames) in os.walk('.'):
			for file in filenames:
				path = os.path.join(dirpath, file)
				if path[:2] == '.\\': path = path[2:]
				if bool(re.match(regex, path)):
					files.append(path)
	else:
		for file in os.listdir(os.curdir):
			if os.path.isfile(file) and bool(re.match(regex, file)):
				files.append(file)

	print()
	if len(files) == 0:
		print(f'No files were found that match "{regex}"')
		print()
		return ''

	print('List of files:')
	for i, file in enumerate(files):
		print(f'  File {i + 1}  -  {file}')
	print()

	selection = None
	while selection is None:
		try:
			i = int(input(f'Please select a file (1 to {len(files)}): '))
		except KeyboardInterrupt:
			sys.exit()
		except:
			pass
		if i > 0 and i <= len(files):
			selection = files[i - 1]
	print()
	return selection

# Given a regular expression, list the directories that match it, and ask for user input
def selectDir(regex, subdirs = False):
	dirs = []
	if subdirs:
		for (dirpath, dirnames, filenames) in os.walk('.'):
			if dirpath[:2] == '.\\': dirpath = dirpath[2:]
			if bool(re.match(regex, dirpath)):
				dirs.append(dirpath)
	else:
		for obj in os.listdir(os.curdir):
			if os.path.isdir(obj) and bool(re.match(regex, obj)):
				dirs.append(obj)

	print()
	if len(dirs) == 0:
		print(f'No directories were found that match "{regex}"')
		print()
		return ''

	print('List of directories:')
	for i, directory in enumerate(dirs):
		print(f'  Directory {i + 1}  -  {directory}')
	print()

	selection = None
	while selection is None:
		try:
			i = int(input(f'Please select a directory (1 to {len(dirs)}): '))
		except KeyboardInterrupt:
			sys.exit()
		except:
			pass
		if i > 0 and i <= len(dirs):
			selection = dirs[i - 1]
	print()
	return selection

# List the files with a regular expression
def listFiles(regex, directory = ''):
	path = os.path.join(os.curdir, directory)
	return [os.path.join(path, file) for file in os.listdir(path) if os.path.isfile(os.path.join(path, file)) and bool(re.match(regex, file))]

# Change the color of text in the terminal
# Leaving the forground or background blank will reset the color to its default
# Providing a message will return the colored message (reset to default afterwards)
# If it's not working for you, be sure to call os.system('cls') before modifying colors
# Usage:
# - print(color('black', 'white', 'Inverted') + ' Regular')
# - print(color('black', 'white') + 'Inverted' + color() + ' Regular')
def color(foreground = '', background = '', message = ''):
	fg = {
		'red': '1;31',
		'green': '1;32',
		'yellow': '1;33',
		'blue': '1;34',
		'purple': '1;35',
		'cyan': '1;36',
		'white': '1;37',
		'black': '0;30',
		'gray': '1;30'
	}
	bg = {
		'red': ';41m',
		'green': ';42m',
		'yellow': ';43m',
		'blue': ';44m',
		'purple': ';45m',
		'cyan': ';46m',
		'white': ';47m',
		'black': ';48m'
	}
	if foreground not in fg or background not in bg: return '\033[0m' # Reset
	color = f'\033[0m\033[{ fg[foreground.lower()] }{ bg[background.lower()] }'

	if message == '': return color
	else: return f'{ color }{ str(message) }\033[0m'

# Return the +- value after a set of data
def percentile(data, p=95):
	a = 1.0 * np.array(data)
	return np.percentile(a, p)

directory = selectDir('.*Logs.*', True)

csvFiles = listFiles(r'.*\_parsed_handshake.csv', directory)

for inputFileName in csvFiles:
	print(f'Processing "{inputFileName}"...')

	readerFile = open(inputFileName, 'r')
	reader = csv.reader(x.replace('\0', '') for x in readerFile)

	temp_header = next(reader)
	handshake_durations = []

	for row in reader:
		handshake_duration = float(row[3])
		handshake_durations.append(handshake_duration)

	top_95_quantile = percentile(handshake_durations, 95)
	print(f'    Top 95% percentile: {top_95_quantile}')
	readerFile.close()

	# 95 quantile has now been computed
	outputFilePath = os.path.join(directory, os.path.basename(inputFileName[:-4] + '_remove95percentile.csv'))
	outputFile = open(outputFilePath, 'w', newline='')
	outputWriter = csv.writer(outputFile)

	readerFile = open(inputFileName, 'r')
	reader = csv.reader(x.replace('\0', '') for x in readerFile)

	temp_header = next(reader)
	outputWriter.writerow(temp_header)

	removedRows = 0
	totalRows = 0

	for row in reader:
		totalRows += 1
		handshake_duration = float(row[3])
		if handshake_duration >= top_95_quantile:
			removedRows += 1
			continue

		outputWriter.writerow(row)

	readerFile.close()
	outputFile.close()
	print(f'    {removedRows} / {totalRows} rows were removed.')
	print(f'    Saved to {outputFilePath}')

print('Goodbye.')
sys.exit()

max_stars_repo_name: prchander/lsquic

max_issues_repo_name: prchander/lsquic

max_forks_repo_name: prchander/lsquic

reason: import numpy,import scipy


text:
import numpy as np
import csv
import os

import util

# Reproducible random numbers
np.random.seed(0)

# Prefix
path = 'results/hw1/'

if not os.path.exists(path):
    os.mkdir(path)

reportLines = []

# 100 children
res = util.measurement(path, taskId = 'hw1_1', childrenPerLevel = (100, ), initSecPrice = [1, 100, 100, 100], logNormalMean = 2, logNormalSigma = 1)
reportLines.append(res)

# 1000 children
res = util.measurement(path, taskId = 'hw1_2', childrenPerLevel = (1000, ), initSecPrice = [1, 100, 100, 100], logNormalMean = 2, logNormalSigma = 1)
reportLines.append(res)

# 10'000 children
res = util.measurement(path, taskId = 'hw1_3', childrenPerLevel = (10000, ), initSecPrice = [1, 100, 100, 100], logNormalMean = 2, logNormalSigma = 1)
reportLines.append(res)

# 100'000 children
res = util.measurement(path, taskId = 'hw1_4', childrenPerLevel = (100000, ), initSecPrice = [1, 100, 100, 100], logNormalMean = 2, logNormalSigma = 1)
reportLines.append(res)

# 100 children
res = util.measurement(path, taskId = 'hw1_5', childrenPerLevel = (100, ), initSecPrice = [1, 100, 100, 100], logNormalMean = 0, logNormalSigma = 2)
reportLines.append(res)

# 1000 children
res = util.measurement(path, taskId = 'hw1_6', childrenPerLevel = (1000, ), initSecPrice = [1, 100, 100, 100], logNormalMean = 0, logNormalSigma = 2)
reportLines.append(res)

# 10'000 children
res = util.measurement(path, taskId = 'hw1_7', childrenPerLevel = (10000, ), initSecPrice = [1, 100, 100, 100], logNormalMean = 0, logNormalSigma = 2)
reportLines.append(res)

# 100'000 children
res = util.measurement(path, taskId = 'hw1_8', childrenPerLevel = (100000, ), initSecPrice = [1, 100, 100, 100], logNormalMean = 0, logNormalSigma = 2)
reportLines.append(res)

# csv report
with open(path + 'report.csv','w') as report:
    csvWriter = csv.writer(report)
    csvWriter.writerow(['Option price', 'Wall time', 'CPU Time', 'Lognormal Mean', 'Lognormal Sigma', 'Files'])
    for row in reportLines:
        csvWriter.writerow(row)


text:
## Saman Farahmand
## This python script seelct random tiles from each slide

# To add a new cell, type '# %%'
# To add a new markdown cell, type '# %% [markdown]'
# %%
import os
import cv2
import numpy as np
import random
import shutil

# %%
mydir="/project/umb_kourosh_zarringhalam/saman/HER2/images/TCGA_raw_output_filtered"
target_dir="/project/umb_kourosh_zarringhalam/saman/HER2/scripts/"
os.chdir(mydir)
# get a number of tiles you want to randomly select for each slide
rand_num=100
# number of slides based on subdirectories
tile_size=512

# %%
subdirs = [x[0] for x in os.walk(mydir)]
target_dirs=[]
for subdir in subdirs:
    if('20.0' in subdir):
        target_dirs.append(subdir)

# %%
#print(target_dirs)
vert=[]
for subdir in target_dirs:
       print(subdir)
       if not os.path.exists(str(subdir.strip('20.0') + '22.0')):
           os.mkdir(str(subdir.strip('20.0') + '22.0'))
       files = os.listdir(subdir)
       print("lenght is: " + str(len(files)))
       if(len(files)!=0):
          for j in range(rand_num):
             index = random.randrange(0, len(files))
             src=str(subdir +'/' +files[index])
             trg=str(subdir.strip('20.0') + '22.0')
             shutil.copy(src, trg)


text:
import numpy as np
import cv2

from skimage.feature import peak_local_max

NOCS_CAMERA_MAT = np.array([[591.0125 ,   0.     , 322.525  ,   0.     ],
       						[  0.     , 590.16775, 244.11084,   0.     ],
       						[  0.     ,   0.     ,   1.     ,   0.     ],
       						[  0.     ,   0.     ,   0.     ,   1.     ]])

open_3d_lines = [
        [0, 1],
        [7,3],
        [1, 3],
        [2, 0],
        [3, 2],
        [0, 4],
        [1, 5],
        [2, 6],
        # [4, 7],
        [7, 6],
        [6, 4],
        [4, 5],
        [5, 7],
    ]

# Ref: https://github.com/zubair-irshad/CenterSnap/blob/5422258475c30c37807566c60996f4d8b3a810e7/utils/nocs_utils.py#L7
def load_depth(depth_path):
	""" Load depth image from img_path. """
	# depth_path = depth_path + '_depth.png'
	# print("depth_path", depth_path)
	depth = cv2.imread(depth_path, -1)
	if len(depth.shape) == 3:
		# This is encoded depth image, let's convert
		# NOTE: RGB is actually BGR in opencv
		depth16 = depth[:, :, 1]*256 + depth[:, :, 2]
		depth16 = np.where(depth16==32001, 0, depth16)
		depth16 = depth16.astype(np.uint16)
	elif len(depth.shape) == 2 and depth.dtype == 'uint16':
		depth16 = depth
	else:
		assert False, '[ Error ]: Unsupported depth type.'
	return depth16

# Ref: https://github.com/zubair-irshad/CenterSnap/blob/5422258475c30c37807566c60996f4d8b3a810e7/utils/nocs_utils.py#L24
def load_img_NOCS(rgm_img_path, depth_path):

	left_img = cv2.imread(rgm_img_path)
	depth = load_depth(depth_path)
	depth_norm = np.array(depth, dtype=np.float32)/255.0

	return left_img, depth_norm, depth

# Ref: https://github.com/zubair-irshad/CenterSnap/blob/5422258475c30c37807566c60996f4d8b3a810e7/simnet/lib/net/post_processing/pose_outputs.py#L121
def find_nearest(peaks,value):
	newList = np.linalg.norm(peaks-value, axis=1)
	return peaks[np.argsort(newList)]

# Ref: https://github.com/zubair-irshad/CenterSnap/blob/5422258475c30c37807566c60996f4d8b3a810e7/simnet/lib/net/post_processing/pose_outputs.py#L126
def extract_peaks_from_centroid_sorted(centroid_heatmap,min_confidence=0.15, min_distance=10):
	peaks = peak_local_max(centroid_heatmap, min_distance=min_distance, threshold_abs=min_confidence)
	peaks = find_nearest(peaks,[0,0])
	return peaks

# Ref: https://github.com/zubair-irshad/CenterSnap/blob/5422258475c30c37807566c60996f4d8b3a810e7/simnet/lib/net/post_processing/pose_outputs.py#L131
def extract_latent_emb_from_peaks(heatmap_output, peaks, latent_emb_output, scale_factor=8):
	assert peaks.shape[1] == 2
	latent_embeddings = []
	indices = []
	scores = []
	for ii in range(peaks.shape[0]):
		index = np.zeros([2])
		index[0] = int(peaks[ii, 0] / scale_factor)
		index[1] = int(peaks[ii, 1] / scale_factor)
		index = index.astype(np.int)
		latent_emb = latent_emb_output[index[0], index[1], :]
		latent_embeddings.append(latent_emb)
		indices.append(index*scale_factor)
		scores.append(heatmap_output[peaks[ii, 0], peaks[ii, 1]])
	return latent_embeddings, indices, scores

# Ref: https://github.com/zubair-irshad/CenterSnap/blob/5422258475c30c37807566c60996f4d8b3a810e7/simnet/lib/net/post_processing/pose_outputs.py#L147
def extract_abs_pose_from_peaks(peaks, abs_pose_output, scale_factor=8):
	assert peaks.shape[1] == 2
	abs_poses = []
	scales = []
	for ii in range(peaks.shape[0]):
		index = np.zeros([2])
		index[0] = int(peaks[ii, 0] / scale_factor)
		index[1] = int(peaks[ii, 1] / scale_factor)
		index = index.astype(np.int)

		abs_pose_values = abs_pose_output[index[0], index[1],:]
		rotation_matrix = np.array([[abs_pose_values[0], abs_pose_values[1], abs_pose_values[2]],
																[abs_pose_values[3], abs_pose_values[4], abs_pose_values[5]],
																[abs_pose_values[6], abs_pose_values[7], abs_pose_values[8]]])
		translation_vector = np.array([abs_pose_values[9], abs_pose_values[10], abs_pose_values[11]])

		transformation_mat = np.eye(4)
		transformation_mat[:3,:3] = rotation_matrix
		transformation_mat[:3,3] = translation_vector

		scale = abs_pose_values[12]
		scale_matrix = np.eye(4)
		scale_mat = scale*np.eye(3, dtype=float)
		scale_matrix[0:3, 0:3] = scale_mat
		scales.append(scale_matrix)

		abs_poses.append((transformation_mat, scale_matrix))
	return abs_poses

# Ref: https://github.com/zubair-irshad/CenterSnap/blob/5422258475c30c37807566c60996f4d8b3a810e7/simnet/lib/net/post_processing/pose_outputs.py#L177
def draw_peaks(centroid_target, peaks):
	centroid_target = np.clip(centroid_target, 0.0, 1.0) * 255.0
	heatmap_img = cv2.applyColorMap(centroid_target.astype(np.uint8), cv2.COLORMAP_JET)
	for ii in range(peaks.shape[0]):
		point = (int(peaks[ii, 1]), int(peaks[ii, 0]))
		heatmap_img = cv2.putText(heatmap_img,str(ii),
		point,
		cv2.FONT_HERSHEY_SIMPLEX,
		1,
		(255,255,255),
		2)
		cv2.line(heatmap_img, point, (0,0), (0, 255, 0), thickness=3, lineType=8)
	return heatmap_img

# Ref: https://github.com/zubair-irshad/CenterSnap/blob/5422258475c30c37807566c60996f4d8b3a810e7/simnet/lib/net/post_processing/abs_pose_outputs.py#L120
def compute_point_cloud_embeddings(heatmap_output, latent_emb_output, min_confidence):
	peaks = extract_peaks_from_centroid_sorted(np.copy(heatmap_output), min_confidence)
	#peaks_image = None
	peaks_image = draw_peaks(np.copy(heatmap_output), np.copy(peaks))
	latent_embs, indices, scores = extract_latent_emb_from_peaks(
			np.copy(heatmap_output), np.copy(peaks), np.copy(latent_emb_output)
	)
	return latent_embs, peaks, peaks_image,scores, indices

# Ref: https://github.com/zubair-irshad/CenterSnap/blob/5422258475c30c37807566c60996f4d8b3a810e7/simnet/lib/net/post_processing/abs_pose_outputs.py#L129
def compute_pointclouds_and_poses(
		heatmap_output,
		latent_emb_output,
		abs_pose_output,
		min_confidence
):
	latent_embeddings , peaks, img,scores, indices = compute_point_cloud_embeddings(np.copy(heatmap_output), np.copy(latent_emb_output), min_confidence)

	abs_pose_outputs = extract_abs_pose_from_peaks(np.copy(peaks), abs_pose_output)
	return latent_embeddings, abs_pose_outputs, img, scores, indices

# Ref: https://github.com/zubair-irshad/CenterSnap/blob/5422258475c30c37807566c60996f4d8b3a810e7/simnet/lib/camera.py#L106
def convert_points_to_homopoints(points):
	"""Project 3d points (3xN) to 4d homogenous points (4xN)"""
	assert len(points.shape) == 2
	assert points.shape[0] == 3
	points_4d = np.concatenate([
	  points,
	  np.ones((1, points.shape[1])),
	], axis=0)
	assert points_4d.shape[1] == points.shape[1]
	assert points_4d.shape[0] == 4
	return points_4d

# Ref: https://github.com/zubair-irshad/CenterSnap/blob/5422258475c30c37807566c60996f4d8b3a810e7/simnet/lib/camera.py#L119
def convert_homopoints_to_points(points_4d):
	"""Project 4d homogenous points (4xN) to 3d points (3xN)"""
	assert len(points_4d.shape) == 2
	assert points_4d.shape[0] == 4
	points_3d = points_4d[:3, :] / points_4d[3:4, :]
	assert points_3d.shape[1] == points_3d.shape[1]
	assert points_3d.shape[0] == 3
	return points_3d

# Ref: https://github.com/zubair-irshad/CenterSnap/blob/5422258475c30c37807566c60996f4d8b3a810e7/utils/transform_utils.py#L29
def project(K, p_3d):
	projections_2d = np.zeros((2, p_3d.shape[1]), dtype='float32')
	p_2d = np.dot(K, p_3d)
	projections_2d[0, :] = p_2d[0, :]/p_2d[2, :]
	projections_2d[1, :] = p_2d[1, :]/p_2d[2, :]
	return projections_2d

# Ref: https://github.com/zubair-irshad/CenterSnap/blob/5422258475c30c37807566c60996f4d8b3a810e7/utils/transform_utils.py#L118
def calculate_2d_projections(coordinates_3d, intrinsics):
	"""
	Input:
		coordinates: [3, N]
		intrinsics: [3, 3]
	Return
		projected_coordinates: [N, 2]
	"""
	projected_coordinates = intrinsics @ coordinates_3d
	projected_coordinates = projected_coordinates[:2, :] / projected_coordinates[2, :]
	projected_coordinates = projected_coordinates.transpose()
	projected_coordinates = np.array(projected_coordinates, dtype=np.int32)

	return projected_coordinates

# Ref: https://github.com/zubair-irshad/CenterSnap/blob/5422258475c30c37807566c60996f4d8b3a810e7/utils/transform_utils.py#L65
def get_3d_bbox(size, shift=0):
	"""
	Args:
		size: [3] or scalar
		shift: [3] or scalar
	Returns:
		bbox_3d: [3, N]
	"""
	bbox_3d = np.array([[+size[0] / 2, +size[1] / 2, +size[2] / 2],
					[+size[0] / 2, +size[1] / 2, -size[2] / 2],
					[-size[0] / 2, +size[1] / 2, +size[2] / 2],
					[-size[0] / 2, +size[1] / 2, -size[2] / 2],
					[+size[0] / 2, -size[1] / 2, +size[2] / 2],
					[+size[0] / 2, -size[1] / 2, -size[2] / 2],
					[-size[0] / 2, -size[1] / 2, +size[2] / 2],
					[-size[0] / 2, -size[1] / 2, -size[2] / 2]]) + shift
	return bbox_3d

# Ref: https://github.com/zubair-irshad/CenterSnap/blob/5422258475c30c37807566c60996f4d8b3a810e7/utils/transform_utils.py#L102
def transform_coordinates_3d(coordinates, RT):
	"""
	Input:
		coordinates: [3, N]
		RT: [4, 4]
	Return
		new_coordinates: [3, N]
	"""
	assert coordinates.shape[0] == 3
	coordinates = np.vstack([coordinates, np.ones((1, coordinates.shape[1]), dtype=np.float32)])
	new_coordinates = RT @ coordinates
	new_coordinates = new_coordinates[:3, :]/new_coordinates[3, :]
	return new_coordinates

# Ref: https://github.com/zubair-irshad/CenterSnap/blob/5422258475c30c37807566c60996f4d8b3a810e7/utils/transform_utils.py#L37
def get_pointclouds_3d(pose, pc):

	pc_homopoints = convert_points_to_homopoints(pc.T)
	morphed_pc_homopoints = pose[0] @ (pose[1] @ pc_homopoints)
	morphed_pc_homopoints = convert_homopoints_to_points(morphed_pc_homopoints).T

	pc_hp = convert_points_to_homopoints(pc.T)
	scaled_homopoints = (pose[1] @ pc_hp)
	scaled_homopoints = convert_homopoints_to_points(scaled_homopoints).T
	size = 2 * np.amax(np.abs(scaled_homopoints), axis=0)
	box = get_3d_bbox(size)
	unit_box_homopoints = convert_points_to_homopoints(box.T)
	morphed_box_homopoints = pose[0] @ unit_box_homopoints
	morphed_box_points = convert_homopoints_to_points(morphed_box_homopoints).T

	xyz_axis = 0.3*np.array([[0, 0, 0], [0, 0, 1], [0, 1, 0], [1, 0, 0]]).transpose()
	sRT = pose[0] @ pose[1]
	transformed_axes = transform_coordinates_3d(xyz_axis, sRT)

	return morphed_pc_homopoints, morphed_box_points, transformed_axes, size

def get_pointclouds_2d(pointcloud_3d, box_3d, axes_3d, camera_mat):

	points_3d_hp = convert_points_to_homopoints(pointcloud_3d.T)
	points_2d_mesh = project(camera_mat, points_3d_hp)
	points_2d_mesh = points_2d_mesh.T

	box_3d_hp = convert_points_to_homopoints(np.array(box_3d).T)
	box_2d = project(camera_mat, box_3d_hp)
	box_2d = box_2d.T

	projected_axes = calculate_2d_projections(axes_3d, camera_mat[:3,:3])

	return points_2d_mesh, box_2d, projected_axes

# Ref: https://github.com/zubair-irshad/CenterSnap/blob/5422258475c30c37807566c60996f4d8b3a810e7/utils/viz_utils.py#L238
def draw_bboxes(img, img_pts, axes, color):

	img_pts = np.int32(img_pts).reshape(-1, 2)
	# draw ground layer in darker color

	# color_ground = (int(color[0]*0.3), int(color[1]*0.3), int(color[2]*0.3))
	color_ground = (int(color[0]), int(color[1]), int(color[2]))

	for i, j in zip([4, 5, 6, 7], [5, 7, 4, 6]):
		img = cv2.line(img, tuple(img_pts[i]), tuple(img_pts[j]), color_ground, 2)
	# draw pillars in minor darker color
	# color_pillar = (int(color[0]*0.6), int(color[1]*0.6), int(color[2]*0.6))
	color_pillar = (int(color[0]), int(color[1]), int(color[2]))
	for i, j in zip(range(4), range(4, 8)):
		img = cv2.line(img, tuple(img_pts[i]), tuple(img_pts[j]), color_pillar, 2)
	# draw top layer in original color
	for i, j in zip([0, 1, 2, 3], [1, 3, 0, 2]):
		img = cv2.line(img, tuple(img_pts[i]), tuple(img_pts[j]), color, 2)

	# draw axes
	img = cv2.arrowedLine(img, tuple(axes[0]), tuple(axes[1]), (0, 0, 255), 2)
	img = cv2.arrowedLine(img, tuple(axes[0]), tuple(axes[3]), (255, 0, 0), 2)
	img = cv2.arrowedLine(img, tuple(axes[0]), tuple(axes[2]), (0, 255, 0), 4) ## y last

	return img