text:
import torch
import torch.nn as nn
from torch.autograd import Variable
import torch.nn.functional as F
from medpy.metric import jc
import logging
import nibabel as nib

import numpy as np
import os

try:
    import cv2
except:
    logging.warning('Could not import opencv. Augmentation functions will be unavailable.')
else:
    def rotate_image(img, angle, interp=cv2.INTER_LINEAR):

        rows, cols = img.shape[:2]
        rotation_matrix = cv2.getRotationMatrix2D((cols / 2, rows / 2), angle, 1)
        return cv2.warpAffine(img, rotation_matrix, (cols, rows), flags=interp)

    def rotate_image_as_onehot(img, angle, nlabels, interp=cv2.INTER_LINEAR):

        onehot_output = rotate_image(convert_to_onehot(img, nlabels=nlabels), angle, interp)
        return np.argmax(onehot_output, axis=-1)

    def resize_image(im, size, interp=cv2.INTER_LINEAR):

        im_resized = cv2.resize(im, (size[1], size[0]), interpolation=interp)  # swap sizes to account for weird OCV API
        return im_resized

    def resize_image_as_onehot(im, size, nlabels, interp=cv2.INTER_LINEAR):

        onehot_output = resize_image(convert_to_onehot(im, nlabels), size, interp=interp)
        return np.argmax(onehot_output, axis=-1)

    def deformation_to_transformation(dx, dy):

        nx, ny = dx.shape

        # grid_x, grid_y = np.meshgrid(np.arange(nx), np.arange(ny))
        grid_y, grid_x = np.meshgrid(np.arange(nx), np.arange(ny), indexing="ij")  # Robin's change to make it work with non-square images

        map_x = (grid_x + dx).astype(np.float32)
        map_y = (grid_y + dy).astype(np.float32)

        return map_x, map_y

    def dense_image_warp(im, dx, dy, interp=cv2.INTER_LINEAR, do_optimisation=True):

        map_x, map_y = deformation_to_transformation(dx, dy)

        # The following command converts the maps to compact fixed point representation
        # this leads to a ~20% increase in speed but could lead to accuracy losses
        # Can be uncommented
        if do_optimisation:
            map_x, map_y = cv2.convertMaps(map_x, map_y, dstmap1type=cv2.CV_16SC2)
        return cv2.remap(im, map_x, map_y, interpolation=interp, borderMode=cv2.BORDER_REFLECT) #borderValue=float(np.min(im)))

    def dense_image_warp_as_onehot(im, dx, dy, nlabels, interp=cv2.INTER_LINEAR, do_optimisation=True):

        onehot_output = dense_image_warp(convert_to_onehot(im, nlabels), dx, dy, interp, do_optimisation=do_optimisation)
        return np.argmax(onehot_output, axis=-1)

def truncated_normal_(tensor, mean=0, std=1):
    size = tensor.shape
    tmp = tensor.new_empty(size + (4,)).normal_()
    valid = (tmp < 2) & (tmp > -2)
    ind = valid.max(-1, keepdim=True)[1]
    tensor.data.copy_(tmp.gather(-1, ind).squeeze(-1))
    tensor.data.mul_(std).add_(mean)

def init_weights(m):
    if type(m) == nn.Conv2d or type(m) == nn.ConvTranspose2d:
        nn.init.kaiming_normal_(m.weight, mode='fan_in', nonlinearity='relu')
        #nn.init.normal_(m.weight, std=0.001)
        #nn.init.normal_(m.bias, std=0.001)
        truncated_normal_(m.bias, mean=0, std=0.001)

def init_weights_orthogonal_normal(m):
    if type(m) == nn.Conv2d or type(m) == nn.ConvTranspose2d:
        nn.init.orthogonal_(m.weight)
        truncated_normal_(m.bias, mean=0, std=0.001)
        #nn.init.normal_(m.bias, std=0.001)

def l2_regularisation(m):
    l2_reg = None

    for W in m.parameters():
        if l2_reg is None:
            l2_reg = W.norm(2)
        else:
            l2_reg = l2_reg + W.norm(2)
    return l2_reg

def normalise_image(image):
    '''
    make image zero mean and unit standard deviation
    '''

    img_o = np.float32(image.copy())
    m = np.mean(img_o)
    s = np.std(img_o)
    return np.divide((img_o - m), s + 1e-6)

def normalise_images(X):
    '''
    Helper for making the images zero mean and unit standard deviation i.e. `white`
    '''

    X_white = np.zeros(X.shape, dtype=np.float32)

    for ii in range(X.shape[0]):

        Xc = X[ii,...]
        X_white[ii,...] = normalise_image(Xc)

    return X_white.astype(np.float32)

def ncc(a,v, zero_norm=True):

    a = a.flatten()
    v = v.flatten()

    if zero_norm:

        a = (a - np.mean(a)) / (np.std(a) * len(a))
        v = (v - np.mean(v)) / np.std(v)

    else:

        a = (a) / (np.std(a) * len(a))
        v = (v) / np.std(v)

    return np.correlate(a, v)

def generalised_energy_distance(sample_arr, gt_arr, nlabels=1, **kwargs):

    def dist_fct(m1, m2):

        label_range = kwargs.get('label_range', range(nlabels))

        per_label_iou = []
        for lbl in label_range:

            # assert not lbl == 0  # tmp check
            m1_bin = (m1 == lbl)*1
            m2_bin = (m2 == lbl)*1

            if torch.sum(m1_bin) == 0 and torch.sum(m2_bin) == 0:
                per_label_iou.append(1)
            elif torch.sum(m1_bin) > 0 and torch.sum(m2_bin) == 0 or torch.sum(m1_bin) == 0 and torch.sum(m2_bin) > 0:
                per_label_iou.append(0)
            else:
                per_label_iou.append(jc(m1_bin.detach().cpu().numpy(), m2_bin.detach().cpu().numpy()))

        # print(1-(sum(per_label_iou) / nlabels))

        return 1-(sum(per_label_iou) / nlabels)

    """
    :param sample_arr: expected shape N x X x Y
    :param gt_arr: M x X x Y
    :return:
    """

    N = sample_arr.shape[0]
    M = gt_arr.shape[0]

    d_sy = []
    d_ss = []
    d_yy = []

    for i in range(N):
        for j in range(M):
            # print(dist_fct(sample_arr[i,...], gt_arr[j,...]))
            d_sy.append(dist_fct(sample_arr[i,...], gt_arr[j,...]))

    for i in range(N):
        for j in range(N):
            # print(dist_fct(sample_arr[i,...], sample_arr[j,...]))
            d_ss.append(dist_fct(sample_arr[i,...], sample_arr[j,...]))

    for i in range(M):
        for j in range(M):
            # print(dist_fct(gt_arr[i,...], gt_arr[j,...]))
            d_yy.append(dist_fct(gt_arr[i,...], gt_arr[j,...]))

    return (2./(N*M))*sum(d_sy) - (1./N**2)*sum(d_ss) - (1./M**2)*sum(d_yy)

def variance_ncc_dist(sample_arr, gt_arr):

    def pixel_wise_xent(m_samp, m_gt, eps=1e-8):

        log_samples = np.log(m_samp + eps)

        return -1.0*np.sum(m_gt*log_samples, axis=0)

    """
    :param sample_arr: expected shape N x X x Y
    :param gt_arr: M x X x Y
    :return:
    """
    sample_arr = sample_arr.detach().cpu().numpy()
    gt_arr = gt_arr.detach().cpu().numpy()

    mean_seg = np.mean(sample_arr, axis=0)

    N = sample_arr.shape[0]
    M = gt_arr.shape[0]

    sX = sample_arr.shape[2]
    sY = sample_arr.shape[3]

    E_ss_arr = np.zeros((N,sX,sY))
    for i in range(N):
        E_ss_arr[i,...] = pixel_wise_xent(sample_arr[i,...], mean_seg)
        # print('pixel wise xent')
        # plt.imshow( E_ss_arr[i,...])
        # plt.show()

    E_ss = np.mean(E_ss_arr, axis=0)

    E_sy_arr = np.zeros((M,N, sX, sY))
    for j in range(M):
        for i in range(N):
            E_sy_arr[j,i, ...] = pixel_wise_xent(sample_arr[i,...], gt_arr[j,...])

    E_sy = np.mean(E_sy_arr, axis=1)

    ncc_list = []
    for j in range(M):

        ncc_list.append(ncc(E_ss, E_sy[j,...]))

    return (1/M)*sum(ncc_list)

def show_tensor(tensor):
    """Show images with matplotlib for debugging, only for 128,128"""
    with torch.no_grad():
        import matplotlib.pyplot as plt
        try:
            tensor = tensor.detach()
        except:
            pass

        height = tensor.shape[2]
        width = tensor.shape[3]

        batch_size = tensor.shape[0]

        result = tensor[0].view(height, width)
        for i in range(1, batch_size):
            result = torch.cat([result, tensor[i].view(height, width)], dim=1)

        plt.imshow(result, cmap='Greys_r')

# def convert_to_onehot(lblmap, nlabels):
#
#     output = torch.zeros((lblmap.shape[0], nlabels, lblmap.shape[2], lblmap.shape[3]))
#     for ii in range(nlabels):
#         output[:, ii, :, :] = (lblmap == ii).view(-1, lblmap.shape[2], lblmap.shape[3]).long()
#
#     assert output.shape == (lblmap.shape[0], nlabels, lblmap.shape[2], lblmap.shape[3])
#
#     return output
def convert_to_onehot(lblmap, nlabels):

    output = np.zeros((lblmap.shape[0], lblmap.shape[1], nlabels))
    for ii in range(nlabels):
        output[:, :, ii] = (lblmap == ii).astype(np.uint8)
    return output

# needs a torch tensor as input instead of numpy array
# accepts format HW and CHW
def convert_to_onehot_torch(lblmap, nlabels):
    if len(lblmap.shape) == 3:
        # 2D image
        output = torch.zeros((nlabels, lblmap.shape[-2], lblmap.shape[-1]))
        for ii in range(nlabels):
            lbl = (lblmap == ii).view(lblmap.shape[-2], lblmap.shape[-1])
            output[ii, :, :] = lbl
    elif len(lblmap.shape) == 4:
        # 3D images from brats are already one hot encoded
        output = lblmap
    return output.long()

def convert_batch_to_onehot(lblbatch, nlabels):
    out = []
    for ii in range(lblbatch.shape[0]):
        lbl = convert_to_onehot_torch(lblbatch[ii,...], nlabels)
        # TODO: check change
        out.append(lbl.unsqueeze(dim=0))

    result = torch.cat(out, dim=0)
    return result

def makefolder(folder):
    '''
    Helper function to make a new folder if doesn't exist
    :param folder: path to new folder
    :return: True if folder created, False if folder already exists
    '''
    if not os.path.exists(folder):
        os.makedirs(folder)
        return True
    return False

def convert_nhwc_to_nchw(tensor):
    result = tensor.transpose(1, 3).transpose(2, 3)
    return result

def convert_nchw_to_nhwc(tensor):
    result = tensor.transpose(1, 3).transpose(1, 2)
    assert torch.equal(tensor, convert_nhwc_to_nchw(result))
    return result

def setup_logger(name, log_file, level=logging.INFO):
    """To setup as many loggers as you want"""
    formatter = logging.Formatter('%(asctime)s %(levelname)s %(message)s')

    handler = logging.FileHandler(log_file, mode='w')
    handler.setFormatter(formatter)

    logger = logging.getLogger(name)
    logger.setLevel(level)
    logger.addHandler(handler)

    return logger

def load_nii(img_path):

    '''
    Shortcut to load a nifti file
    '''

    nimg = nib.load(img_path)
    return nimg.get_data(), nimg.affine, nimg.header

def save_nii(img_path, data, affine, header):
    '''
    Shortcut to save a nifty file
    '''

    nimg = nib.Nifti1Image(data, affine=affine, header=header)
    nimg.to_filename(img_path)

def create_and_save_nii(data, img_path):

    img = nib.Nifti1Image(data, np.eye(4))
    nib.save(img, img_path)

max_stars_repo_name: gigantenbein/UNet-Zoo


text:
import numpy as np
from PIL import Image
from recordclass import recordclass

FrameInfo = recordclass("FrameInfo", "cur_time index_hitobj info_index osr_index index_fp obj_endtime x_end y_end, break_index")
CursorEvent = recordclass("CursorEvent", "event old_x old_y")

def get_buffer(img, settings):
	np_img = img.reshape((settings.height, settings.width, 4))
	pbuffer = Image.fromarray(np_img, "RGBA")
	pbuffer.readonly = False
	return np_img, pbuffer

ext: py

reason: import numpy


text:
from strategy.astrategy import AStrategy
from datetime import datetime, timedelta, timezone
from sklearn.impute import SimpleImputer
import numpy as np
import math
import pandas as pd
import pickle

class WeeklyPrice(AStrategy):
    def __init__(self):
        """
        Quarterly Financials Strat
        Attrs:
        year = year of model
        quarter = quarter of model
        ticker = ticker of model
        yearly_gap = number of years to project out to
        """
        super().__init__("weekly_price")
        self.suffix = "weekly_price"

    def transform(self,datasets,year,quarter,ticker,week_gap,training_years):
        """
        Required Datasets:
            1) funds: financial data in format ["year","quarter","financial metrics"]
            2) price: price data in format ["year","quarter","price]
        """
        drop_columns = ["_id"]
        ticker_regression = datasets["regression"]
        ticker_classification = datasets["classification"]
        ticker_classification.rename(columns={ticker:"y_class"},inplace=True)
        ticker_regression["y"] = ticker_regression[ticker]
        ticker_regression = ticker_regression.merge(ticker_classification[["year","quarter","week","y_class"]],how="left",on=["year","quarter","week"])
        ticker_regression["y"] = ticker_regression["y"].shift(-week_gap)
        ticker_regression["y_class"] = ticker_regression["y_class"].shift(-week_gap)
        ticker_regression = ticker_regression[:-week_gap]
        ticker_regression.dropna(inplace=True)
        ticker_regression.reset_index(inplace=True,drop=True)
        first = ticker_regression[(ticker_regression["year"] == year - training_years) & (ticker_regression["quarter"] == quarter)].index.values.tolist()[0]
        last = ticker_regression[(ticker_regression["year"] == year) & (ticker_regression["quarter"] == quarter)].index.values.tolist()[0]
        ticker_regression.drop(drop_columns,axis=1,inplace=True,errors="ignore")
        training_data = ticker_regression.iloc[first:last-1]
        prediction_data = ticker_regression[(ticker_regression["year"] == year) & (ticker_regression["quarter"] == quarter)]
        factors = list(prediction_data.columns)
        for column in prediction_data.columns:
            contains = True in [x == -99999 for x in prediction_data[column]]
            if contains:
                factors.remove(column)
            else:
                continue
        return training_data[factors], prediction_data[factors]

    def create_sim(self,ticker,year,quarter,models,factors,data):
            """
            Required Datasets:
                1) funds: financial data in format ["year","quarter","financial metrics"]
                2) price: price data in format ["year","quarter","price]
            """
            ## Creating Simulation Data
            current_sets = []
            prediction_set = data["prediction_set"]
            price_dict = []
            for week in prediction_set["week"]:
                price_dict.append({"year":year,"quarter":quarter,"week":week})
            price = pd.DataFrame(price_dict)
            for j in range(len(models)):
                model = models.iloc[j]
                price_model = model["model"]
                prediction = price_model.predict(prediction_set[factors])
                current_set = price.copy()
                col_name = model["model_type"]
                current_set["{}_{}_prediction".format(self.suffix,col_name)] = prediction
                current_set["{}_{}_score".format(self.suffix,col_name)] = model["score"]
                current_sets.append(current_set)
            sim = current_sets[0]
            for cs in current_sets[1:]:
                sim = sim.merge(cs,on=["year","quarter","week"],how="left")
            sim["ticker"] = ticker
            return sim

    ### This needs debugging
    def create_rec(self,models,factors,data):
        """
        Required Datasets:
            1) funds: financial data in format ["year","quarter","financial metrics"]
            2) price: price data in format ["year","quarter","price]
        """
        ## Creating Simulation Data
        # funds = data["funds"]
        # price = {"ticker":self.ticker,"year":self.year+self.yearly_gap,"quarter":self.quarter}
        # for i in range(len(models)):
        #     model = models.iloc[i]
        #     final = funds[(funds["year"] == self.year) & (funds["quarter"] == self.quarter)]
        #     prediction = model["model"].predict(final)[0]
        #     if model["api"] == "tf":
        #         p = []
        #         [p.extend(x) for x in prediction]
        #         prediction = p
        #     column_name = "quarterly_fundamental_{}".format(model["model_type"])
        #     price["{}_prediction".format(column_name)] = prediction
        #     price["{}_score".format(column_name)] = model["score"]
        # sim = pd.DataFrame([price])
        # self.db.connect()
        # self.db.store_data("rec",sim)
        # self.db.close()
        return "lel"


text:
import gym
from gym import (error, spaces)
from gym.utils import seeding
import numpy as np
from pydrake.all import (
    DiagramBuilder,
    Simulator,
)

class DrakeEnv(gym.Env):
    metadata = {'render.modes': ['human']}

    def __init__(self):
        '''
        Sets up the System diagram and creates a drake visualizer object to
        send LCM messages to during the render method.

        Subclasses must implement the methods below that throw a
        NotImplementedError
        '''

        # Create the Diagram.
        self.system = self.plant_system()
        self.context = self.system.CreateDefaultContext()

        # Create the simulator.
        self.simulator = Simulator(self.system, self.context)
        self.simulator.set_publish_every_time_step(False)

    def seed(self, seed=None):
        self.np_random, seed = seeding.np_random(seed)
        return [seed]

    def plant_system(self):
        '''
        Returns the fully constructed MDP diagram which is connected to the
        vector source for simulation. Each subclass should implement this
        method.
        '''
        raise NotImplementedError

    def get_input_port_action(self):
        '''
        Returns the system input port that corresponds to the action
        '''
        raise NotImplementedError

    def get_observation(self):
        '''
        Returns the system output port that corresponds to the observation
        '''
        raise NotImplementedError

    def get_state(self):
        '''
        Returns the system output port that corresponds to the observation
        '''
        raise NotImplementedError

    @property
    def action_space(self):
        '''
        Specifies the limits of tha action space. This must be overridden by subclasses.
        '''
        raise NotImplementedError

    @property
    def observation_space(self):
        '''
        Specifies the limits of tha action space. This must be overridden by subclasses.
        '''
        raise NotImplementedError

    def get_reward(self, state, action):
        '''
        Computes the reward from the state and an action
        '''
        raise NotImplementedError

    def reset_state(self):
        '''
        Computes the reward from the state and an action
        '''
        raise NotImplementedError

    def step(self, action):
        '''
        Simulates the system diagram for a short period of time
        '''
        # Only allow valid actions from the action space
        action = validate_action(action, self.action_space)

        # Fix input port with action and simulate
        action_fixed_input_port_value = self.context.FixInputPort(
            self.get_input_port_action().get_index(), action)
        self.simulator.StepTo(self.context.get_time() + self.dt)

        # Get the observation, reward, and terminal conditions
        state = self.get_state()
        observation = self.get_observation()
        reward = self.get_reward(state, action)
        done = self.is_done()
        info = {}

        return observation, reward, done, info

    def reset(self):
        '''
        Resets the state in the system diagram
        '''
        self.context.set_time(0)
        self.reset_state()
        return self.get_observation()

    def render(self, mode='human', close=False):
        '''
        Notifies the visualizer to draw the current state. Different implementations based on MultiBodyPlant and RigidBodyTree
        '''
        raise NotImplementedError

def validate_action(action, action_space):
    '''
    Ensure the action is within the action space
    '''
    if type(action_space) is spaces.Box:
        action = np.minimum(action, action_space.high)
        action = np.maximum(action, action_space.low)

    return action

max_stars_repo_name: HarvardAgileRoboticsLab/gym-drake

max_issues_repo_name: HarvardAgileRoboticsLab/gym-drake

max_forks_repo_name: HarvardAgileRoboticsLab/gym-drake


text:
#coding: utf-8

#######################################################
#Description:
# load julius dictation-kit-v4.4 dnn model data
# and calculate dnn by deep learning framework chainer
# This is based on dnnclient.py in dictation-kit-v4.4.zip of
#
# Copyright (c) 1991-2016 Kawahara Lab., Kyoto University
# Copyright (c) 2005-2016 Lee Lab., Nagoya Institute of Technology
#
# License: see LICENSE-Julius Dictation Kit.txt
#
#Date 2018-03-20
#By Shun
#########################################################

import os
import numpy as np
import chainer
from chainer import cuda, Variable
from chainer import Link, Chain, ChainList
import chainer.links as L
import chainer.functions as F

# Check version
#  Python 2.7.12 on win32 (Windows version)
#  numpy (1.14.0)
#  chainer (1.22.0)

class Class_chainer_DNN(Chain):
	def __init__(self, in_size, mid_size, out_size, net=None, state_prior=None):
		super(Class_chainer_DNN, self).__init__(

		l1 = L.Linear(in_size,mid_size,initialW=net.fc1_W if net else None ,initial_bias=net.fc1_b if net else None),
		l2 = L.Linear(mid_size,mid_size,initialW=net.fc2_W if net else None ,initial_bias=net.fc2_b if net else None),
		l3 = L.Linear(mid_size,mid_size,initialW=net.fc3_W if net else None ,initial_bias=net.fc3_b if net else None),
		l4 = L.Linear(mid_size,mid_size,initialW=net.fc4_W if net else None ,initial_bias=net.fc4_b if net else None),
		l5 = L.Linear(mid_size,mid_size,initialW=net.fc5_W if net else None ,initial_bias=net.fc5_b if net else None),
		l6 = L.Linear(mid_size,mid_size,initialW=net.fc6_W if net else None ,initial_bias=net.fc6_b if net else None),
		l7 = L.Linear(mid_size,mid_size,initialW=net.fc7_W if net else None ,initial_bias=net.fc7_b if net else None),
		l8 = L.Linear(mid_size,out_size,initialW=net.fc8_W if net else None ,initial_bias=net.fc8_b if net else None),
		)
		self.state_prior=state_prior

	def __call__(self,x):
		# dnn
		x1 = F.sigmoid( self.l1(x))
		x2 = F.sigmoid( self.l2(x1))
		x3 = F.sigmoid( self.l3(x2))
		x4 = F.sigmoid( self.l4(x3))
		x5 = F.sigmoid( self.l5(x4))
		x6 = F.sigmoid( self.l6(x5))
		x7 = F.sigmoid( self.l7(x6))
		# discriminator
		tmp = F.softmax(self.l8(x7))

		# apply state prior probability and log10
		if self.state_prior is not None:
			tmp /= self.state_prior
			tmp= F.log10(tmp)

		return tmp

class Class_net(object):
	def __init__(self, IN_DIR='model/dnn'):
		self.IN_DIR= IN_DIR
		self.W, self.b = self.load_init_DNN_Wb()
		self.num_states= self.W[7].shape[1]  # output dimension
		self.num_input= self.W[0].shape[0]   # input dimension
		self.num_mid= self.W[0].shape[1]     # middle dimension
		self.prior = self.get_prior()
	@property
	def fc1_W(self):
		return self.W[0].T
	@property
	def fc1_b(self):
		return self.b[0].reshape(self.b[0].size)
	@property
	def fc2_W(self):
		return self.W[1].T
	@property
	def fc2_b(self):
		return self.b[1].reshape(self.b[1].size)
	@property
	def fc3_W(self):
		return self.W[2] .T
	@property
	def fc3_b(self):
		return self.b[2].reshape(self.b[2].size)
	@property
	def fc4_W(self):
		return self.W[3].T
	@property
	def fc4_b(self):
		return self.b[3].reshape(self.b[3].size)
	@property
	def fc5_W(self):
		return self.W[4].T
	@property
	def fc5_b(self):
		return self.b[4].reshape(self.b[4].size)
	@property
	def fc6_W(self):
		return self.W[5].T
	@property
	def fc6_b(self):
		return self.b[5].reshape(self.b[5].size)
	@property
	def fc7_W(self):
		return self.W[6].T
	@property
	def fc7_b(self):
		return self.b[6].reshape(self.b[6].size)
	@property
	def fc8_W(self):
		return self.W[7].T
	@property
	def fc8_b(self):
		return self.b[7].reshape(self.b[7].size)

	def load_init_DNN_Wb(self,):
		w_filename = ["W_l1_f4.npy", "W_l2_f4.npy", "W_l3_f4.npy", "W_l4_f4.npy", "W_l5_f4.npy", "W_l6_f4.npy", "W_l7_f4.npy", "W_output_f4.npy"]
		b_filename = ["bias_l1_f4.npy", "bias_l2_f4.npy", "bias_l3_f4.npy", "bias_l4_f4.npy", "bias_l5_f4.npy", "bias_l6_f4.npy", "bias_l7_f4.npy", "bias_output_f4.npy"]
		listw=['fc1_W','fc2_W','fc3_W','fc4_W','fc5_W','fc6_W','fc7_W','fc8_W']
		listb=['fc1_b','fc2_b','fc3_b','fc4_b','fc5_b','fc6_b','fc7_b','fc8_b']
		for i, listx in enumerate(w_filename):
			f=os.path.join(self.IN_DIR,listx)
			if os.path.exists(f):
				listw[i] = np.load(f)
				print (' load of ', f)
			else:
				listw[i] = None
				print (' no file of ', f)
		for i, listx in enumerate(b_filename):
			f=os.path.join(self.IN_DIR,listx)
			if os.path.exists(f):
				listb[i] = np.load(f)
				print (' load of ', f)
			else:
				listb[i] = None
				print (' no file of ', f)
		return listw,listb

	def get_prior(self,):
		prior_filename = 'prior'
		f=os.path.join(self.IN_DIR,prior_filename)
		assert os.path.exists(f), 'error, no file of ' + f
		print (' open of ', f)
		state_prior = np.zeros(self.num_states)
		prior_factor = 1.0
		for line in open(f):
			state_id, state_p = line[:-1].split(' ')
			state_id = int(state_id)
			state_p = float(state_p) * prior_factor
			state_prior[state_id] = state_p

		return state_prior


text:
import numpy as np
import timedelta
from datetime import datetime as dt
from datetime import timedelta as td

class Boss():
    def __init__(self, original_name, name, chance, location, respawn_hours):
        self.original_name = str(original_name.lower())
        self.name = str(name)
        self.chance = int(chance)
        self.location = str(location)
        self.respawn_hours = int(respawn_hours)
        self._generate_sound_names()

    def _generate_sound_names(self):
        self.sound_5min = "_".join([self.original_name.lower().replace(" ","-"),"5min"])
        self.sound_appear = "_".join([self.original_name.lower().replace(" ","-"),"appear"])

    def __lt__(self, other):
        return self.respawn_time_seconds < other.respawn_time_seconds

    def calc_respawn_time(self, force=False, respawn_time=None):
        self.respawn_time = (respawn_time if force else dt.now()) + td(hours=self.respawn_hours) - td(seconds=1)
        self.respawn_time_seconds = self.respawn_time.timestamp()

    def last_time(self):
        duration =  timedelta.Timedelta(self.respawn_time - dt.now())
        is_valid = self.respawn_time > dt.now()
        days, seconds = duration.days, duration.seconds
        hours = days * 24 + seconds // 3600
        minutes = (seconds % 3600) // 60
        seconds = seconds % 60
        return (hours, minutes, seconds, is_valid)

    def boss_tagged_str(self):
        (hours, minutes, seconds, _) = self.last_time()
        return(f"Босс **{self.name}** убит!\nЯ напомню о его респе через {hours}ч {minutes}м {seconds}с")

    def boss_status_str(self):
        (hours, minutes, seconds, is_valid) = self.last_time()
        moscow_time = self.respawn_time.strftime('%H:%M:%S')
        if is_valid and hours == 0 and minutes == 0 and seconds == 0:
            return(f"{moscow_time} **{self.name}** в {self.location} появился! | шанс {self.chance}% | {moscow_time} мск")
        elif is_valid:
            return(f"**{self.name}** в {self.location} через {hours}ч {minutes}м {seconds}с | шанс {self.chance}% | {moscow_time} мск")
        else:
            return(f"**{self.name}** уже реснулся {59-minutes}м {59-seconds}с назад! | шанс {self.chance}% | {moscow_time} мск")

    def export_msg(self):
        date_str = self.respawn_time.strftime('%Y-%m-%d %H:%M:%S')
        return(f"+{self.name.upper()} {date_str}")

    @staticmethod
    def from_export_str(astring, is_manual=False):
        line = astring.split()
        boss_name = str(line[0][1:])
        resp_last = str(" ".join(line[1:]))
        boss = bosses_dict.get(boss_name.upper())
        resp_last = dt.fromisoformat(resp_last) if is_manual else dt.fromisoformat(resp_last) - td(hours=boss.respawn_hours)
        boss.calc_respawn_time(force=True, respawn_time=resp_last)
        return boss

    def not_tagged(self):
        return(f"**{self.name}** должен был реснуться.\nНапиши +{self.name.lower()} для сброса таймера.")

    def boss_short_str(self):
        return(f"**{self.name}**")

    def __str__(self):
        return (f"{self.name}|{self.chance}|{self.location}|{self.respawn_hours}\n{self.respawn_time.strftime('%Y-%m-%d %H:%M:%S')}")

data = np.genfromtxt('./assets/data/bosses.csv', names=True, delimiter=',', dtype=None, encoding='utf8')
bosses_dict = dict(list(map(lambda x: (x[4], Boss(x[1], x[4], x[2], x[5], x[3])), data)))


text:
import numpy as np
import functools
import scipy.optimize as op

from math_helper import MathHelper
from model import Model

class NeuralNetwork:
    def __init__(self, lambda_val, input_layer_size, output_layer_size, hidden_layer_sizes, initial_thetas=None):
        if lambda_val is None or (type(lambda_val) != float and type(lambda_val) != int):
            raise TypeError("Lambda must be a numeric type.")

        if lambda_val < 0.:
            raise ValueError("Lambda must be a non-negative number.")

        if type(input_layer_size) != int or type(output_layer_size) != int:
            raise TypeError("Input and output layer sizes must be in int.")

        if input_layer_size <= 0 or output_layer_size <= 0:
            raise ValueError("Input and output layer sizes must be greater than zero.")

        if len(hidden_layer_sizes) < 1:
            raise ValueError("There must be at least one hidden layer.")

        if not all(isinstance(x, int) and x > 0 for x in hidden_layer_sizes):
            raise ValueError("All hidden layer sizes must be positive integer.")

        self.lambda_val = lambda_val
        self.input_layer_size = input_layer_size
        self.output_layer_size = output_layer_size
        self.hidden_layer_sizes = list(hidden_layer_sizes)
        self.hidden_layer_count = len(self.hidden_layer_sizes)
        self._initial_thetas = self._random_initialize_theta() if initial_thetas is None else [np.asarray(t) for t in
                                                                                initial_thetas]
        self.helper = MathHelper()

    @classmethod
    def init(cls, lambda_val, input_layer_size, output_layer_size, hidden_layer_sizes):
        return cls(lambda_val=lambda_val, input_layer_size=input_layer_size, output_layer_size=output_layer_size,
                   hidden_layer_sizes=hidden_layer_sizes, initial_thetas=None)

    @classmethod
    def init_with_theta(cls, lambda_val, thetas):
        if len(thetas) < 2:
            raise ValueError("There must be at least one hidden layer and hence at least two weight matrices.")

        def _shape(item):
            if hasattr(item, "shape"):
                return item.shape
            r = len(item)
            c = 0
            if r > 0:
                c = len(item[0])
            return r, c

        ils = _shape(thetas[0])[1] - 1
        hls = [_shape(t)[0] for t in thetas]
        ols = hls.pop()

        return cls(lambda_val=lambda_val, input_layer_size=ils, output_layer_size=ols, hidden_layer_sizes=hls,
                   initial_thetas=thetas)

    def _random_initialize_theta(self):
        delta = 0.12

        def random_matrix(row, col):
            return np.random.rand(row, col) * 2 * delta - delta

        first_size_array = [self.input_layer_size]
        [first_size_array.append(x) for x in self.hidden_layer_sizes]
        paired_copy = first_size_array[1:]
        paired_copy.append(self.output_layer_size)
        sizes = zip(paired_copy, first_size_array)

        return [random_matrix(r, c + 1) for (r, c) in sizes]

    @staticmethod
    def _unroll_matrices(matrices):
        def mapper(x):
            return np.ravel(x).tolist()

        def reducer(x, y):
            return x + y

        combined_list = functools.reduce(reducer, map(mapper, matrices))
        return np.array(combined_list)

    def _roll_into_matrices(self, unrolled_vector):
        taken, prev_layer_size, current_layer_size = 0, self.input_layer_size, self.hidden_layer_sizes[0]
        matrices = []
        for i in range(self.hidden_layer_count):
            current_layer_size = self.hidden_layer_sizes[i]
            matrices.append(np.reshape(unrolled_vector[taken: taken + current_layer_size * (prev_layer_size + 1)],
                                       (current_layer_size, prev_layer_size + 1)))
            taken += current_layer_size * (prev_layer_size + 1)
            prev_layer_size = current_layer_size

        matrices.append(np.reshape(unrolled_vector[taken:], (self.output_layer_size, prev_layer_size + 1)))
        return matrices

    def cost_regularization(self, current_thetas, train_data_size):
        if self.lambda_val == 0:
            return 0

        def mapper(x):
            matrix = np.zeros(x.shape)
            matrix[:, 1:] = x[:, 1:]
            return np.multiply(matrix, matrix).sum()

        def reducer(x, y):
            return x + y

        return self.lambda_val * functools.reduce(reducer, map(mapper, current_thetas), 0) / (2. * train_data_size)

    def _theta_regularization(self, theta, train_data_size):
        result = np.zeros(theta.shape)

        if self.lambda_val == 0:
            return result

        result[:, 1:] = (self.lambda_val / train_data_size) * theta[:, 1:]
        return result

    def _calculate_cost_gradient(self, unrolled_theta_vector, X_in, Y_in):
        thetas = self._roll_into_matrices(unrolled_theta_vector)
        theta_count = len(thetas)
        X = np.asmatrix(X_in)
        Y = np.asmatrix(Y_in)
        m, n = X.shape
        eps = 0.#np.finfo(float).eps

        def single_sample_mapper(pair):
            x, y = pair
            if np.any((y != 0)  & (y != 1)):
                raise ValueError(
                    "Output value cannot be anything other than 0 and 1. If you want more than two level of output, try converting the out values into a vector of 0 and 1 only.")

            activations = []
            prev_activation = np.ones((1, n+1))
            prev_activation[:,1:] = np.asmatrix(x)[:,:]
            for i in range(theta_count):
                theta = np.asmatrix(thetas[i])
                activations.append(prev_activation)
                zt = np.dot(prev_activation, np.transpose(theta))
                at = np.ones((zt.shape[0], zt.shape[1] + 1))
                at[:, 1:] = self.helper.sigmoid(zt)
                prev_activation = at

            ht = prev_activation[:, 1:]

            # for y==1 and y==0 items separately
            j_partial = sum(-np.log(ht[y == 1] + eps))
            j_partial += sum(-np.log(1 + eps - ht[y == 0]))

            deltas = []

            prev_delta = np.transpose(ht - y)
            deltas.append(prev_delta)

            for i in range(1, theta_count):
                a = np.transpose(activations[-i])

                delta = np.multiply(np.dot(np.transpose(thetas[-i]), prev_delta), np.multiply(a, 1 - a))[1:,:]
                deltas.append(delta)
                prev_delta = delta

            DELs = [np.dot(d, a) for (d,a) in zip(deltas[::-1], activations)]

            return {'cost': j_partial, 'deltas' : DELs}

        def sample_pair_reducer(x, y):
            return {'cost': x['cost'] + y['cost'], 'deltas': [dx + dy for (dx,dy) in zip(x['deltas'], y['deltas'])]}

        result = functools.reduce(
            sample_pair_reducer,
            map(single_sample_mapper, zip(X, Y)),
            {'cost': 0, 'deltas': [np.zeros(t.shape) for t in thetas]});

        J = result["cost"]/m + self.cost_regularization(thetas, m)

        gradients = [(d/m)+self._theta_regularization(t, m) for (d,t) in zip(result["deltas"], thetas)]

        return (J, self._unroll_matrices(gradients))

    def train(self, X, Y, maxiter = 200, tolerance = 1e-5):
        if X.shape[0] != Y.shape[0]:
            raise(ValueError("X and Y must have same number of rows."))
        t = self._unroll_matrices(self._initial_thetas)
        res = op.minimize(
            fun = self._calculate_cost_gradient,
            x0 = t,
            jac = True,
            tol = tolerance,
            method = 'CG',
            options = {'maxiter': maxiter},
            args = (X, Y)
        );

        return Model(self._roll_into_matrices(res.x))

max_stars_repo_name: sasoana/multiclassANN

max_issues_repo_name: sasoana/multiclassANN

max_forks_repo_name: sasoana/multiclassANN

reason: import numpy,import scipy


text:
import numpy as np
import numpy.random as random
import matplotlib.pyplot as pyplot

#Urbanscape v1.4

#Changes:
# randomized the income and rent distribution functions with a random multiplier
# added an element of randomness to the FAFH and FFR functions using a random float multiplier
# changed the FFR function such that the domain of incomes that yield a percent < 0.5% and beyond follows the function f(x) = 0.005
# added a GroceryStore Agent
# added an attribute to the UrbanScape class for income spent on groceries

#-------------------------------------
# Basic Classes: UrbanScape and Agent |
#-------------------------------------

class UrbanScape(object):
	# The urbanscape defines the context within which agents act
	population_per_block = 20

	#size = size of the urban grid
	#rent = yearly cost of commercial and residential space
	def __init__(self, size, rent, create_rule=None, gradient = None, randomize = False):
		blocks = np.ones((size,size))

		self.create_rule = create_rule
		self.size = size
		self.agents = []
		self.agent_locations = np.zeros((self.size,self.size))
		self.agent_coords = {}
		self.time = 0
		self.rent_ceiling = rent
		self.rent_floor = 5000
		self.rent = np.ones_like(blocks) * self.rent_ceiling	#distribution of rent on urbanscape

		if gradient == 'vertical':
			self.vertical_distribution()
		elif gradient == 'diagonal':
			self.diagonal_distribution()
		elif gradient == 'random':
			self.random_distribution()
		elif gradient == 'BDquadrants':
			self.businessdistricts_quadrant_distribution()
		elif gradient == 'CBD':
			self.centralbusinessdistrict_distribution()

		self.income = np.ones_like(blocks) * self.rent * 4

		if randomize == True:
			self.randomize_distribution()

		self.externalities = np.zeros_like(blocks)		#total count of (negative) externality effect
		self.ffcapture_number = np.zeros_like(blocks)	#count of FastFoodAgent effect in UrbanScape
		self.gscapture_number = np.zeros_like(blocks)   #count of GroceryStoreAgent effect in UrbanScape
		self.mobility = np.ones_like(blocks)			#a measure of disability
		self.food_away = np.ones_like(blocks)			#distribution of food-related expenditures per block
		self.fast_food = np.ones_like(blocks)
		self.grocery = np.ones_like(blocks)
		self.update_expenditures()

	def add_agent(self, agent):
		self.agents.append(agent)

	def remove_agent(self, agent):
		self.agents.remove(agent)

	# defines effect radius within particular radius and location
	def effect_radius(self, loc, radius):
		r = radius
		effect_coordinates = []
		x,y = loc
		for i in range(x-r,x+r+1):
			for j in range(y-r,y+r+1):
				effect_coordinates.append((i,j))
		return effect_coordinates

	# defines the $ amount spent on fast food and groceries in blocks of given coordinates
	def capture_expenditures(self, coords):
		coords = coords
		ffexposure_total = 0
		gsexposure_total = 0

		for k in range(len(coords)):
			x,y = coords[k]

			if 0 <= x <= (self.size-1) and 0 <= y <= (self.size-1):
				ffexposure = (self.fast_food[x,y] / (self.ffcapture_number[x,y] + 1e-10)) * self.population_per_block
				ffexposure_total += ffexposure

				gsexposure = (self.grocery[x,y] / (self.gscapture_number[x,y] + 1e-10)) * self.population_per_block
				gsexposure_total += gsexposure

		return ffexposure_total, gsexposure_total

	def update_agent_locations(self):
		#agent locations are specified as -1 if FastFoodAgent, 1 if GroceryStoreAgent, and -0.5 if block has both
		#these values correspond to red, green, and orange on the 'RdYlGn' colormap on matplotlib.pyplot
		self.agent_locations = np.zeros((self.size,self.size))
		self.agent_coords = {}
		for aa in self.agents:
			self.agent_coords[aa] = aa.loc     #keeps track of type coord and type of agent

		#this finds redundant and single-occurring coordinates.
		#redundant coordinates are assumed to contain both FastFoodAgents and GroceryStoreAgents
		coords = self.agent_coords.values()
		set_coords = sorted(set(coords))

		for i in set_coords:
			coords.remove(i)
		redundant_locs = coords

		for i in redundant_locs:
			try:
				set_coords.remove(i)
			except:
				pass
		single_locs = set_coords

		for agent, loc in self.agent_coords.items():
			string = str(agent)
			x,y = loc

			if loc in redundant_locs:
				self.agent_locations[x,y] = -0.5

			if loc in single_locs:
				if 'FastFoodAgent' in string:
					self.agent_locations[x,y] = -1

				if 'GroceryStoreAgent' in string:
					self.agent_locations[x,y] = 1

	def update_income(self):
		self.income = (self.rent * 4) * self.mobility

	def update_capture_number(self):
		self.ffcapture_number = np.zeros((self.size,self.size))
		self.gscapture_number = np.zeros((self.size,self.size))
		for agent in self.agents:
			string = str(agent)

			if 'FastFoodAgent' in string:
				for coords in agent.effect_coordinates:
					x,y = coords
					if 0 <= x <= (self.size-1) and 0<= y <= (self.size-1):
						self.ffcapture_number[x,y] += 1

			if 'GroceryStoreAgent' in string:
				for coords in agent.effect_coordinates:
					x,y = coords
					if 0 <= x <= (self.size-1) and 0<= y <= (self.size-1):
						self.gscapture_number[x,y] += 1

	def update_externalities(self):
		#externalities here are 'negative externalities':
		#ffscapture_number adds to externalities, and gscapture_number substracts from it.
		#diminishes externalities by heal_rate if capture_number = 0
		heal_rate = float(0.95)

		for i in range(self.size):
			for j in range(self.size):
				self.externalities[i,j] = self.externalities[i,j] + self.ffcapture_number[i,j] - self.gscapture_number[i,j]

				if self.externalities[i,j] <= 0:
					self.externalities[i,j] = 0

				if (self.ffcapture_number[i,j]-self.gscapture_number[i,j]) == 0:
					self.externalities[i,j] = int(self.externalities[i,j] * heal_rate)

	def update_mobility(self):
		#a logistic decay function
		#decreases mobility to limit of 'a3' as externality count approaches infinity

		a3 = float(0.75)
		b3 = float(1.5)
		c3 = float(50)

		for i in range(self.size):
			for j in range(self.size):
				E = float(self.externalities[i,j])
				M = a3 + ((1 - a3) / (1 + b3**(E - c3)))
				self.mobility[i,j] = M

	def update_expenditures(self):
		#Expenditures on 'Food Away From Home'
		#Defining constants a1, b1
		a1 = 1.5802
		b1 = -0.32

		for i in range(self.size):
			for j in range(self.size):
				rand_multiplier = self.random_float_range(0.75,1.25)
				FA = a1 * (self.income[i,j]**b1)
				self.food_away[i,j] = FA * self.income[i,j] * rand_multiplier

		#Expenditures on 'Fast Food Restaurants'
		#Defining constants a2, b2, c2
		a2 = 5.528e-6
		b2 = -0.332
		c2 = 0.40

		for i in range(self.size):
			for j in range(self.size):
				rand_multiplier = self.random_float_range(0.75,1.25)
				FF = -(a2 * self.income[i,j] + b2)**2 + c2
				if FF < 0.075:
					FF = 0.075
				self.fast_food[i,j] = FF * self.food_away[i,j] * rand_multiplier

		#Expenditures on 'Food At Home', which as assumed to be all grocery expenditures
		a1 = 73.969
		b1 = -0.632

		for i in range(self.size):
			for j in range(self.size):
				rand_multiplier = self.random_float_range(0.75,1.25)
				FH = a1 * (self.income[i,j]**b1)
				self.grocery[i,j] = FH * self.income[i,j] * rand_multiplier

#-----------------------------------------
# Rent and Income Distribution Functions |
#-----------------------------------------

	#Defining types of block rent distributions, which also determines income level
	def vertical_distribution(self):

		gradient = float((self.rent_ceiling - self.rent_floor)/(self.size-1))
		g_multiplier = float(gradient/self.rent_ceiling)
		fr_ratio = float(self.rent_floor)/float(self.rent_ceiling)

		for i in range(self.size):
			gg = float(fr_ratio + g_multiplier*i)
			self.rent[i] = self.rent[i] * gg

	def diagonal_distribution(self):

		gradient = float((self.rent_ceiling - self.rent_floor)/(self.size-1))
		g_multiplier = float(gradient/self.rent_ceiling)
		fr_ratio = float(self.rent_floor/self.rent_ceiling)

		for i in range(self.size):
			gg = float(fr_ratio + g_multiplier*i)
			self.rent[i] = self.rent[i] * gg
		for j in range(self.size):
			gg = float(fr_ratio + g_multiplier*j)
			self.rent[:,j] = self.rent[:,j] * gg

	def random_distribution(self):
		dist = np.random.randint(self.rent_floor, self.rent_ceiling,(self.size,self.size))
		self.rent = dist

        #calculates the rent based on distance from the business district block
        def distance_rent_function(self,quadrant_coords,bd_coords):
                q = quadrant_coords
                bd = bd_coords
                x_dist = float(q[0] - bd[0])
                y_dist = float(q[1] - bd[1])
                dist = np.sqrt(x_dist**2 + y_dist**2)
                return dist

        def exponential_rent_function(self,distance,coord):
                i,j = coord
                Lambda = np.log(float(self.rent_floor)/float(self.rent_ceiling)) * (-1/self.dmax)
                rent_value = self.rent_ceiling * np.e**(-Lambda*distance)
                self.rent[i][j] = rent_value

        def linear_rent_function(self,distance,coord):
                i,j = coord
                Alpha = (float(self.rent_floor)-float(self.rent_ceiling))/self.dmax
                rent_value = self.rent_ceiling + (Alpha * distance)
                self.rent[i][j] = rent_value

        #the approximate center of each quadrant of the urbanscape is a 'business district'
        #the assumption being that rent and income decreases as one moves further away
        def businessdistricts_quadrant_distribution(self):
                district1 = [(i,j)
                    for i in range(int(np.ceil(float(self.size)*0.5)))
                    for j in range(int(np.ceil(float(self.size)*0.5)))]
                district2 = [(i,j)
                    for i in range(int(np.ceil(float(self.size)*0.5)),self.size)
                    for j in range(int(np.ceil(float(self.size)*0.5)))]
                district3 = [(i,j)
                    for i in range(int(np.ceil(float(self.size)*0.5)))
                    for j in range(int(np.ceil(float(self.size)*0.5)),self.size)]
                district4 = [(i,j)
                    for i in range(int(np.ceil(float(self.size)*0.5)),self.size)
                    for j in range(int(np.ceil(float(self.size)*0.5)),self.size)]

                d1_coords = (np.ceil(float(self.size) * 0.25), np.ceil(float(self.size)* 0.25))
                d2_coords = (np.ceil(float(self.size) * 0.75), np.ceil(float(self.size)* 0.25))
                d3_coords = (np.ceil(float(self.size) * 0.25), np.ceil(float(self.size)* 0.75))
                d4_coords = (np.ceil(float(self.size) * 0.75), np.ceil(float(self.size)* 0.75))

                d1_distances = [self.distance_rent_function(q,d1_coords) for q in district1]
                d2_distances = [self.distance_rent_function(q,d2_coords) for q in district2]
                d3_distances = [self.distance_rent_function(q,d3_coords) for q in district3]
                d4_distances = [self.distance_rent_function(q,d4_coords) for q in district4]

                self.dmax = max(d1_distances)

                for i in range(len(district1)):
                    self.exponential_rent_function(d1_distances[i], district1[i])
                    self.exponential_rent_function(d2_distances[i], district2[i])
                    self.exponential_rent_function(d3_distances[i], district3[i])
                    self.exponential_rent_function(d4_distances[i], district4[i])

        def centralbusinessdistrict_distribution(self):
                block_coords = [(i,j)
                    for i in range(self.size)
                    for j in range(self.size)]
                district_coords = (np.ceil(float(self.size) *0.5), np.ceil(float(self.size)*0.5))
                block_distances = [self.distance_rent_function(q,district_coords) for q in block_coords]
                self.dmax = max(block_distances)

                for i in range(len(block_coords)):
                    #self.linear_rent_function(block_distances[i], block_coords[i])
                    self.exponential_rent_function(block_distances[i],block_coords[i])

        #a function that randomizes the rent and income values (independently) within an error of 20%
        def randomize_distribution(self):
                for i in range(self.size):
                    for j in range(self.size):
                        randfloat = self.random_float_range(0.75,1.25)

                        if self.rent[i][j] * randfloat > self.rent_ceiling:
                            pass
                        else:
                            self.rent[i][j] *= randfloat

                for i in range(self.size):
                    for j in range(self.size):
                        randfloat = self.random_float_range(0.75,1.25)

                        if self.income[i][j] * randfloat > (self.rent_ceiling * 4):
                            pass
                        else:
                            self.income[i][j] *= randfloat

        #a function that generates a random float between a specified range
        def random_float_range(self, low, high):
                randfloat = np.random.random() * (high - low) + low
                return randfloat

	def step(self):
		self.create_rule(self)

		#Update externality-related attributes
		self.update_agent_locations()
		self.update_capture_number()
		self.update_externalities()
		self.update_mobility()
		self.update_income()
		self.update_expenditures()

		#Agent actions
		for agent in self.agents:
			agent.step(self)

		#Time step
		self.time += 1

#---------------------------------
# Base class for all Food Agents |
#---------------------------------

class Agent(object):
    def __init__(self, loc, wealth):
		self.loc = loc		# x,y coordinates of location in urbanscape
		self.wealth = wealth	# $ amount

	# Called every time-step
    def step(self, urbanscape):
		self.check_bankruptcy(urbanscape)

	# Remove the agent if broke
    def check_bankruptcy(self, urbanscape):
		if self.wealth < 0:
			self.wealth = 0
			urbanscape.remove_agent(self)

class FastFoodAgent(Agent):
	operations = 50000	# $ per year
	initial_wealth = 10000	# $ amount
	radius = 2	# number of block effected around location
	effect_coordinates = []

	def __init__(self, loc, urbanscape):
		super(FastFoodAgent, self).__init__(loc, self.initial_wealth)
		self.define_radius(urbanscape)
		self.operating_costs = (urbanscape.rent[loc]) + (FastFoodAgent.operations)
	def step(self, urbanscape):
		# 1. Gather revenues for the year
		self.capture_revenue(urbanscape)

		# 2. Perform parent class step function
		super(FastFoodAgent, self).step(urbanscape)

	#this function defines how much of fast food expenditures are
	#captured by the FoodAgent. If there are no competing FastFoodAgents,
	#assumes that FoodAgent captures all of block expenditures.

	def define_radius(self, urbanscape):
		self.effect_coordinates = urbanscape.effect_radius(self.loc, self.radius)

	def capture_revenue(self, urbanscape):
		x,y = self.loc
		self.wealth += (urbanscape.capture_expenditures(self.effect_coordinates)[0]) - self.operating_costs
		#capture_expenditures returns a list of two values, the 0th index being the ff_revenues

class GroceryStoreAgent(Agent):
	operations = 200000	# $ per year
	initial_wealth = 25000	# $ amount
	radius = 2	# number of block effected around location
	effect_coordinates = []

	def __init__(self, loc, urbanscape):
		super(GroceryStoreAgent, self).__init__(loc, self.initial_wealth)
		self.define_radius(urbanscape)
		self.operating_costs = (urbanscape.rent[loc]) + (GroceryStoreAgent.operations)
	def step(self, urbanscape):
		# 1. Gather revenues for the year
		self.capture_revenue(urbanscape)

		# 2. Perform parent class step function
		super(GroceryStoreAgent, self).step(urbanscape)

	#this function defines how much of grocery store expenditures are
	#captured by the FoodAgent. If there are no competing FoodAgents,
	#assumes that FoodAgent captures all of block expenditures.

	def define_radius(self, urbanscape):
		self.effect_coordinates = urbanscape.effect_radius(self.loc, self.radius)

	#assumes that the rest of income spent on food away from home is spent on groceries
	def capture_revenue(self, urbanscape):
		x,y = self.loc
		self.wealth += (urbanscape.capture_expenditures(self.effect_coordinates)[1]) - self.operating_costs
		#capture_expenditures returns a list of two values, the 1st index being the ff_revenues

#------------------------
# Create Rule Functions |
#------------------------

def no_create_rule(urbanscape):
	pass

#creates FastFoodAgents randomly on the grid
def random_create_rule(urbanscape):
	x = random.randint(0,urbanscape.size)
	y = random.randint(0,urbanscape.size)

	if urbanscape.time % 5 == 0:
		agent_type_decision = random.rand()
		if agent_type_decision < 0.5:
			urbanscape.add_agent(FastFoodAgent((x,y),urbanscape))
		else:
			urbanscape.add_agent(GroceryStoreAgent((x,y),urbanscape))

#adds new FastFoodAgent in location with highest potential profit to cost ratio
def profit_probability_create_rule(urbanscape):
	size = urbanscape.size
	potential_fflocations = []
	potential_gslocations = []

	#makes a list of available coordinates in UrbanScape
	#only one FastFoodAgent and one GroceryStoreAgent can occupy a block at any given time
	for i in range(size):
		for j in range(size):
		        ff_locs = [loc for agent,loc in urbanscape.agent_coords.items() if 'FastFoodAgent' in str(agent)]
		        gs_locs = [loc for agent,loc in urbanscape.agent_coords.items() if 'GroceryStoreAgent' in str(agent)]
			if (i,j) not in ff_locs:
				potential_fflocations.append((i,j))
			if (i,j) not in gs_locs:
			        potential_gslocations.append((i,j))

	#makes a list of profit:cost ratios at potential_location coordinates
	potential_fflist = potential_profits(urbanscape, potential_fflocations, agent = 'FF')
	potential_gslist = potential_profits(urbanscape, potential_gslocations, agent = 'GS')

	if not potential_fflist and not potential_gslist:
		print 'no suitable locations at time ' + str(urbanscape.time)

	if potential_fflist:
	       potential_creations(urbanscape, potential_fflocations, potential_fflist,agent='FF')

	if potential_gslist:
	       potential_creations(urbanscape, potential_gslocations, potential_gslist,agent='GS')

#a function that creates agents
def potential_creations(urbanscape, potential_locations, potential_list,agent = ''):
        #generates a random float to make the 'create agent' decision
       	prob_max = 0.15
        rand = random.random()
        maximum = max(potential_list)
	create_locations = []

	#generates a 'create' probability based on profit:cost ratios in potential_list
	#modifies prob_max value by a multiplier based on profit ranges
	for i in range(len(potential_list)):
		prob_multiplier = float(potential_list[i])/float(maximum)
		if prob_multiplier < 0:
			prob_multiplier = 0
		probability = prob_max * prob_multiplier

		#adds the location to create_locations if the 'create agent' decision
		#is less that the create probability
		if rand < probability:
			create_locations.append(potential_locations[i])

	#picks a random location from list of locations that made it through the
	#random 'create agent' decision.
	if create_locations:
		choice = random.randint(0,len(create_locations))
		create_coord = create_locations[choice]

		if agent == 'FF':
		      urbanscape.add_agent(FastFoodAgent(create_coord, urbanscape))

		if agent == 'GS':
		      urbanscape.add_agent(GroceryStoreAgent(create_coord, urbanscape))

def potential_profits(urbanscape, potential_locations, agent = ''):
	urbanscape = urbanscape
	locations = potential_locations
	location_profit_list = []
	for loc in locations:
	        if agent == 'FF':
	               startup = (urbanscape.rent[loc]) + (FastFoodAgent.operations)*2
	               potential_radius = urbanscape.effect_radius(loc,FastFoodAgent.radius)
	               revenue = potential_revenue(urbanscape,potential_radius,agent)
	               profits_short = revenue - startup
	               profits_long = (revenue - FastFoodAgent.operations) * 5
	               profit_margin = (profits_short/startup) + (profits_long / (FastFoodAgent.operations*5))
	               location_profit_list.append(profit_margin)

	        if agent == 'GS':
	               startup = (urbanscape.rent[loc]) + (GroceryStoreAgent.operations)*2
	               potential_radius = urbanscape.effect_radius(loc,GroceryStoreAgent.radius)
	               revenue = potential_revenue(urbanscape,potential_radius,agent)
	               profits_short = revenue - startup
	               profits_long = (revenue - GroceryStoreAgent.operations) * 5
	               profit_margin = (profits_short/startup) + (profits_long / (GroceryStoreAgent.operations*5))
	               location_profit_list.append(profit_margin)

	return location_profit_list

#define create rule variables
def potential_revenue(urbanscape, coords, agent = ''):
		coords = coords
		s = urbanscape.size
		ff = urbanscape.fast_food
		gs = urbanscape.grocery
		ff_cn = urbanscape.ffcapture_number
		gs_cn = urbanscape.gscapture_number
		capture_total = 0

		if agent == 'FF':
		      for k in range(len(coords)):
		          x,y = coords[k]
		          if 0 <= x <= (s-1) and 0 <= y <= (s-1):
		              cc = (ff[x,y] / (ff_cn[x,y]+1)) * urbanscape.population_per_block
		              capture_total += cc

		if agent == 'GS':
		      for k in range(len(coords)):
		          x,y = coords[k]
		          if 0 <= x <= (s-1) and 0 <= y <= (s-1):
		              cc = (gs[x,y] / (gs_cn[x,y]+1)) * urbanscape.population_per_block
		              capture_total += cc

		return capture_total

#---------------------------------------
# Functions for Visualizing UrbanScape |
#---------------------------------------

def plot_urbanscape(urbanscape):
	fig = pyplot.figure()

	pyplot.clf()
	pyplot.subplot(221)
	ii = urbanscape.income
	pyplot.imshow(ii, interpolation='nearest')
	pyplot.title("Income Distribution")
	pyplot.subplot(221).axes.get_xaxis().set_ticks([])
	pyplot.subplot(221).axes.get_yaxis().set_ticks([])

	pyplot.subplot(222)
	ff = urbanscape.ffcapture_number
	pyplot.imshow(ff, cmap = 'YlOrRd', interpolation='nearest')
	pyplot.title("FFR Effect Radii")
	pyplot.subplot(222).axes.get_xaxis().set_ticks([])
	pyplot.subplot(222).axes.get_yaxis().set_ticks([])

	pyplot.subplot(223)
	mm = urbanscape.externalities
	pyplot.imshow(mm, cmap = 'Purples', interpolation='nearest')
	pyplot.title("Negative Externalities")
	pyplot.subplot(223).axes.get_xaxis().set_ticks([])
	pyplot.subplot(223).axes.get_yaxis().set_ticks([])

	pyplot.subplot(224)
	ll = urbanscape.agent_locations
	pyplot.imshow(ll, cmap = 'RdYlGn', interpolation='nearest',vmin=-1,vmax=1)
	pyplot.title("Food Agent Locations")
	pyplot.subplot(224).axes.get_xaxis().set_ticks([])
	pyplot.subplot(224).axes.get_yaxis().set_ticks([])

	#return fig

#------------------------------------------------------------------------------
# Running Simulations that Returns Externalities Exposures by Income Quintile |
#------------------------------------------------------------------------------

def run_experiment(grid_size, rent_ceiling, create_rule, distribution, randomize = True, steps = 100):

	n = grid_size
	ceiling = rent_ceiling
	create_rule = create_rule
	steps = steps

	#defining the urbanscape
	#u = UrbanScape(20,10000,profit_probability_create_rule,'random')
	u = UrbanScape(n, ceiling, create_rule, distribution, randomize)

	externality_quintiles = np.array(([],[],[],[],[]))

	for i in range(steps):
		u.step()
		income_distribution = set()
		for uu in u.income:
			for ii in uu:
				income_distribution.add(ii)

		idist_sorted = sorted(income_distribution)
		num = len(idist_sorted)
		quintile_cutoffs = sorted([int((float(num)/float(5))*i - 1) for i in range(1,6)])

		incomes_q1 = [idist_sorted[i] for i in range(0,quintile_cutoffs[0])]
		incomes_q2 = [idist_sorted[i] for i in range(quintile_cutoffs[0],quintile_cutoffs[1])]
		incomes_q3 = [idist_sorted[i] for i in range(quintile_cutoffs[1],quintile_cutoffs[2])]
		incomes_q4 = [idist_sorted[i] for i in range(quintile_cutoffs[2],quintile_cutoffs[3])]
		incomes_q5 = [idist_sorted[i] for i in range(quintile_cutoffs[3],quintile_cutoffs[4])]

		externalities_q1 = []
		externalities_q2 = []
		externalities_q3 = []
		externalities_q4 = []
		externalities_q5 = []

		for i in range(n):
			for j in range(n):

				if u.income[i,j] in incomes_q1:
					externalities_q1.append(u.externalities[i,j])
				elif u.income[i,j] in incomes_q2:
					externalities_q2.append(u.externalities[i,j])
				elif u.income[i,j] in incomes_q3:
					externalities_q3.append(u.externalities[i,j])
				elif u.income[i,j] in incomes_q4:
					externalities_q4.append(u.externalities[i,j])
				elif u.income[i,j] in incomes_q5:
					externalities_q5.append(u.externalities[i,j])

		avg_ext_q1 = sum(externalities_q1)/len(externalities_q1)
		avg_ext_q2 = sum(externalities_q2)/len(externalities_q2)
		avg_ext_q3 = sum(externalities_q3)/len(externalities_q3)
		avg_ext_q4 = sum(externalities_q4)/len(externalities_q4)
		avg_ext_q5 = sum(externalities_q5)/len(externalities_q5)

		externality_quintiles = np.append(externality_quintiles, [[avg_ext_q1],[avg_ext_q2],[avg_ext_q3],[avg_ext_q4],[avg_ext_q5]],axis=1)

	terminal_mobility_dist = np.ones_like(u.mobility) * u.mobility
	terminal_income_dist = np.ones_like(u.income) * u.income

	return externality_quintiles, terminal_mobility_dist, terminal_income_dist, u.agent_coords

def run_batch_experiments(batches, grid_size, rent_ceiling, create_rule, distribution, randomize=True, steps = 100):
	batches = batches
	total_externalities = np.zeros((5,steps))
	total_terminal_mobility = np.zeros((grid_size,grid_size))
	total_terminal_income = np.zeros((grid_size,grid_size))
	total_ffagent_locations = np.zeros((grid_size,grid_size))
	total_gsagent_locations = np.zeros((grid_size,grid_size))

	for i in range(batches):
		experiment = run_experiment(grid_size, rent_ceiling, create_rule, distribution, randomize, steps)
		total_externalities += experiment[0]
		total_terminal_mobility += experiment[1]
		total_terminal_income += experiment[2]
		agent_coords = experiment[3]

		ffagent_locations = np.zeros((grid_size,grid_size))
		gsagent_locations = np.zeros((grid_size,grid_size))

		for agent, loc in agent_coords.items():
		        string = str(agent)
		        x,y = loc

		        if 'FastFoodAgent' in string:
		            ffagent_locations[x,y] = 1

		        if 'GroceryStoreAgent' in string:
		            gsagent_locations[x,y] = 1

		total_ffagent_locations += ffagent_locations
		total_gsagent_locations += gsagent_locations

	avg_externalities = total_externalities/batches
	avg_terminal_mobility = total_terminal_mobility/batches
	avg_terminal_income = total_terminal_income/batches
	avg_ffagent_locations = total_ffagent_locations/batches
	avg_gsagent_locations = total_gsagent_locations/batches

	ax1 = pyplot.subplot(321)
	pyplot.plot(avg_externalities[0], color = 'm', label = 'low-income',linewidth=2.0)
	pyplot.plot(avg_externalities[1], color = 'r', label = 'low-middle income',linewidth=2.0)
	pyplot.plot(avg_externalities[2], color = 'y', label = 'middle income',linewidth=2.0)
	pyplot.plot(avg_externalities[3], color = 'g', label = 'middle-high income',linewidth=2.0)
	pyplot.plot(avg_externalities[4], color = 'b', label = 'high income',linewidth=2.0)
	pyplot.legend(bbox_to_anchor=(1.05, 1), loc=2, prop={'size':7}, borderaxespad=0.)

	pyplot.xlabel('time',fontsize=10,fontname='serif')
	pyplot.ylabel('total exposure',fontsize=10,fontname='serif')
	pyplot.title('Accumulation of Negative Externalities',fontsize=10,fontname='serif')

	ax2 = pyplot.subplot(323)
	pyplot.imshow(avg_terminal_income, interpolation='nearest')
	pyplot.title('Terminal Income Distribution',fontsize=10,fontname='serif')

	ax3 = pyplot.subplot(324)
	pyplot.imshow(avg_terminal_mobility, cmap = 'Purples', interpolation='nearest')
	pyplot.title('Terminal Mobility Distribution',fontsize=10,fontname='serif')

	ax4 = pyplot.subplot(325)
	pyplot.imshow(avg_ffagent_locations, cmap = 'Reds', interpolation='nearest')
	pyplot.title('Fast Food Agent Locations',fontsize=10,fontname='serif')

	ax5 = pyplot.subplot(326)
	pyplot.imshow(avg_gsagent_locations, cmap = 'Greens', interpolation='nearest')
	pyplot.title('Grocery Store Agent Locations',fontsize=10,fontname='serif')

	ax1.tick_params(axis='x',labelsize=9)
        ax1.tick_params(axis='y',labelsize=9)
	ax2.tick_params(axis='x',labelsize=9)
        ax2.tick_params(axis='y',labelsize=9)
        ax3.tick_params(axis='x',labelsize=9)
        ax3.tick_params(axis='y',labelsize=9)
        ax4.tick_params(axis='x',labelsize=9)
        ax4.tick_params(axis='y',labelsize=9)
        ax5.tick_params(axis='x',labelsize=9)
        ax5.tick_params(axis='y',labelsize=9)

	pyplot.show()

def plot_experiment(externality_quintiles, urbanscape):
	pyplot.subplot(221)
	pyplot.plot(externality_quintiles[0], color = 'm', label = 'low-income',linewidth=2.0)
	pyplot.plot(externality_quintiles[1], color = 'r', label = 'low-middle income',linewidth=2.0)
	pyplot.plot(externality_quintiles[2], color = 'y', label = 'middle income',linewidth=2.0)
	pyplot.plot(externality_quintiles[3], color = 'g', label = 'middle-high income',linewidth=2.0)
	pyplot.plot(externality_quintiles[4], color = 'b', label = 'high income',linewidth=2.0)
	pyplot.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)

	pyplot.xlabel('time')
	pyplot.ylabel('total exposure')
	pyplot.title('Accumulation of Negative Externalities')

	pyplot.subplot(223)
	pyplot.imshow(urbanscape.income, interpolation='nearest')
	pyplot.title('Income Distribution')

	pyplot.subplot(224)
	pyplot.imshow(urbanscape.agent_locations, cmap = 'YlOrRd', interpolation='nearest')
	pyplot.title('Fast Food Locations')
	pyplot.show()

#EXPERIMENTS TO RUN ON SHELL

#u = UrbanScape(5, 100000,no_create_rule,'vertical', randomize = False)

#u = UrbanScape(20,200000,profit_probability_create_rule,'random')
#u = UrbanScape(20,200000,profit_probability_create_rule,'BDquadrants')
#u = UrbanScape(20,200000,profit_probability_create_rule,'CBD')

#run_experiment(20,250000,profit_probability_create_rule,'random',steps=10)
#run_experiment(20,250000,profit_probability_create_rule,'BDquadrants',steps=10)
#run_experiment(20,250000,profit_probability_create_rule,'CBD',steps=100)

#run_batch_experiments(3,20,200000,random_create_rule,'random',steps=50)
#run_batch_experiments(3,20,200000,profit_probability_create_rule,'random',steps=50)
#run_batch_experiments(3,20,200000,profit_probability_create_rule,'CBD',steps=200)

#run_batch_experiments(50,20,250000,profit_probability_create_rule,'random',steps=200)
#run_batch_experiments(50,20,250000,profit_probability_create_rule,'vertical',steps=200)
#run_batch_experiments(50,20,250000,profit_probability_create_rule,'CBD',steps=200)

max_stars_repo_name: cosmicBboy/Urbanscape

max_issues_repo_name: cosmicBboy/Urbanscape

max_forks_repo_name: cosmicBboy/Urbanscape


text:
"""
General utility functions.
"""

import glob
import warnings
from pathlib import Path

import numpy as np
from osgeo import gdal
from PIL import Image

def nodata(image):
    """
    Return no-data flag, and mask of no-data pixels, if any.

    Parameters
    ----------
    image: GDAL image
        Input image

    Returns
    -------
    has_nodata: bool
        If the image contains no-data pixels
    """
    band = image.GetRasterBand(1)
    value = band.GetNoDataValue()
    array = band.ReadAsArray()
    has_nodata = value in array
    mask_nodata = np.zeros(array.shape, dtype=bool)
    if has_nodata:
        mask_nodata = value != array
    return has_nodata, mask_nodata

def convert_image(input_path, output_path, suffix='.jpg', backend='gdal', n_channels=1,
                  min_value=None, max_value=None, exclude_nodata=False, mode='I;16'):
    """
    Convert images to a specified format.

    Parameters
    ----------
    mode: str
        Image mode if as greyscale, support backend pillow only
    input_path: str
        Path to input image
    output_path: str
        Path to output image
    suffix: str
        Suffix of filetype to convert to, 'jpg' supported
    backend: str
        Backend of conversion, 'pillow' or 'gdal'
    n_channels: int
        Number of channels to output
    min_value: None or int
        Min value to scale to
    max_value: None or int
        Max value to scale to
    exclude_nodata: bool
        exclude images with no-data pixels if set True
    """
    assert suffix == ".jpg" or suffix == ".jpeg" or suffix == ".png"
    if exclude_nodata and backend != 'gdal':
        backend = 'gdal'
        warnings.warn('excluded_nodata enabled. backend switched to gdal')

    Path(output_path).parent.mkdir(exist_ok=True, parents=True)

    if backend == 'gdal':
        image = gdal.Open(input_path)
        has_nodata, _ = nodata(image)
        if exclude_nodata and has_nodata:
            # print(f'{input_path} contains no-data pixel(s). excluded from processing.')
            return input_path

        options_list = [
            f'-ot Byte',
            f'-of JPEG',
            f'-b {n_channels}',
            f'-scale {min_value} {max_value}' if min_value is not None and max_value is not None else f'-scale'
        ]
        options_string = " ".join(options_list)

        gdal.Translate(
            output_path,
            input_path,
            options=options_string
        )

    elif backend == 'pillow':
        image = Image.open(input_path)
        if n_channels == 1:
            image = image.convert(mode)
        elif n_channels == 3:
            image = image.convert('RGB')
        else:
            raise ValueError(f'n_channels has to be 1 or 3. got {n_channels} instead.')
        image.save(output_path, quality=100)
    else:
        raise ValueError(f'backend has to be gdal or pillow. got {backend} instead.')

def glob_tile_ids(pattern):
    """
    Get available tile folders.

    Parameters
    ----------
    pattern: str
        Pattern path to search for

    Returns
    -------
    as_list: list[Path]
        list of available file paths
    """
    tile_ids = glob.glob(pattern)
    return [Path(i).stem for i in tile_ids]

max_stars_repo_name: chenzhaiyu/ahninpainter

max_issues_repo_name: chenzhaiyu/ahninpainter

max_forks_repo_name: chenzhaiyu/ahninpainter


text:
"""
Wrappers for drawing plots for reports.
"""

from matplotlib import pyplot as plt
from matplotlib import patheffects
from typing import List, Tuple, Optional
import numpy as np
import itertools
from pathlib import Path
import matplotlib.colors as colors
from matplotlib import gridspec
from matplotlib.collections import LineCollection
from math import floor

def time_series_plot(
        outpath: Optional[Path],
        title: str,
        xtitle: str,
        xunit: str,
        ytitle: str,
        yunit: str,
        xdata: List,
        ydata: List,
        trimxvalues: bool = True,
        skipfirst: bool = False,
        figsize: Tuple = (15, 8.5),
        bins: int = 20):
    """
    Draws time series plot.

    Used i.e. for timeline of resource usage.

    It also draws the histogram of values that appeared throughout the
    experiment.

    Parameters
    ----------
    outpath : Optional[Path]
        Output path for the plot image. If None, the plot will be displayed.
    title : str
        Title of the plot
    xtitle : str
        Name of the X axis
    xuint : str
        Unit for the X axis
    ytitle : str
        Name of the Y axis
    yunit : str
        Unit for the Y axis
    xdata : List
        The values for X dimension
    ydata : List
        The values for Y dimension
    trimxvalues : bool
        True if all values for the X dimension should be subtracted by
        the minimal value on this dimension
    skipfirst: bool
        True if the first entry should be removed from plotting.
    figsize: Tuple
        The size of the figure
    bins: int
        Number of bins for value histograms
    """
    start = 1 if skipfirst else 0
    xdata = np.array(xdata[start:], copy=True)
    ydata = np.array(ydata[start:], copy=True)
    if trimxvalues:
        minx = min(xdata)
        xdata = [x - minx for x in xdata]
    fig, (axplot, axhist) = plt.subplots(
        ncols=2,
        tight_layout=True,
        figsize=figsize,
        sharey=True,
        gridspec_kw={'width_ratios': (8, 3)}
    )
    fig.suptitle(title, fontsize='x-large')
    axplot.scatter(xdata, ydata, c='purple', alpha=0.5)
    xlabel = xtitle
    if xunit is not None:
        xlabel += f' [{xunit}]'
    ylabel = ytitle
    if yunit is not None:
        ylabel += f' [{yunit}]'
    axplot.set_xlabel(xlabel, fontsize='large')
    axplot.set_ylabel(ylabel, fontsize='large')
    axplot.grid()

    axhist.hist(ydata, bins=bins, orientation='horizontal', color='purple')
    axhist.set_xscale('log')
    axhist.set_xlabel('Value histogram', fontsize='large')
    axhist.grid(which='both')
    plt.setp(axhist.get_yticklabels(), visible=False)

    if outpath is None:
        plt.show()
    else:
        plt.savefig(outpath)

def draw_confusion_matrix(
        confusion_matrix: np.ndarray,
        outpath: Optional[Path],
        title: str,
        class_names: List[str],
        cmap=None,
        figsize: Optional[Tuple] = None,
        dpi: Optional[int] = None):
    """
    Creates a confusion matrix plot.

    Parameters
    ----------
    confusion_matrix : ArrayLike
        Square numpy matrix containing the confusion matrix.
        0-th axis stands for ground truth, 1-st axis stands for predictions
    outpath : Optional[Path]
        Path where the plot will be saved. If None, the plot will be displayed.
    title : str
        Title of the plot
    class_names : List[str]
        List of the class names
    cmap : Any
        Color map for the plot
    figsize : Optional[Tuple]
        The size of the plot
    dpi : Optional[int]
        The dpi of the plot
    """
    if cmap is None:
        cmap = plt.get_cmap('BuPu')

    confusion_matrix = np.array(confusion_matrix, dtype=np.float32, copy=True)

    # normalize confusion matrix
    confusion_matrix /= confusion_matrix.sum(axis=0)

    # compute overall accuracy
    accuracy = np.trace(confusion_matrix) / np.sum(confusion_matrix)

    # compute sensitivity
    correctactual = confusion_matrix.diagonal() / confusion_matrix.sum(axis=1)
    correctactual = correctactual.reshape(1, len(class_names))

    # compute precision
    correctpredicted = \
        confusion_matrix.diagonal() / confusion_matrix.sum(axis=0)
    correctpredicted = correctpredicted.reshape(len(class_names), 1)

    if figsize is None:
        figsize = confusion_matrix.shape

    if dpi is None:
        dpi = 216

    # create axes
    fig = plt.figure(figsize=figsize, dpi=dpi)
    gs = gridspec.GridSpec(len(class_names) + 1, len(class_names) + 1)
    axConfMatrix = fig.add_subplot(gs[0:len(class_names), 0:len(class_names)])
    axPredicted = fig.add_subplot(
        gs[len(class_names), 0:len(class_names)],
        sharex=axConfMatrix
    )
    axActual = fig.add_subplot(
        gs[0:len(class_names), len(class_names)],
        sharey=axConfMatrix
    )
    axTotal = fig.add_subplot(
        gs[len(class_names), len(class_names)],
        sharex=axActual,
        sharey=axPredicted
    )

    # define ticks for classes
    ticks = np.arange(len(class_names))

    # configure and draw confusion matrix
    axConfMatrix.set_xticks(ticks)
    axConfMatrix.set_xticklabels(class_names, fontsize='large', rotation=90)
    axConfMatrix.set_yticks(ticks)
    axConfMatrix.set_yticklabels(class_names, fontsize='large')
    axConfMatrix.set_xlabel('Actual class', fontsize='x-large')
    axConfMatrix.set_ylabel('Predicted class', fontsize='x-large')
    img = axConfMatrix.imshow(
        confusion_matrix,
        norm=colors.LogNorm(),
        interpolation='nearest',
        cmap=cmap,
        aspect='auto',
    )
    axConfMatrix.xaxis.set_ticks_position('top')
    axConfMatrix.xaxis.set_label_position('top')

    # add percentages for confusion matrix
    for i, j in itertools.product(
            range(len(class_names)),
            range(len(class_names))):
        txt = axConfMatrix.text(
            i, j,
            ('100' if confusion_matrix[i, j] == 1.0
                else f'{100.0 * confusion_matrix[i,j]:3.1f}'),
            ha='center',
            va='center',
            color='black',
            fontsize='medium')
        txt.set_path_effects([
            patheffects.withStroke(linewidth=5, foreground='w')
        ])

    # configure and draw sensitivity percentages
    axPredicted.set_xticks(ticks)
    axPredicted.set_yticks([0])
    axPredicted.set_xlabel('Sensitivity', fontsize='large')
    axPredicted.imshow(
        correctactual,
        interpolation='nearest',
        cmap='RdYlGn',
        aspect='auto',
        vmin=0.0,
        vmax=1.0
    )
    for i in range(len(class_names)):
        txt = axPredicted.text(
            i, 0,
            ('100' if correctactual[0, i] == 1.0
                else f'{100.0 * correctactual[0, i]:3.1f}'),
            ha='center',
            va='center',
            color='black',
            fontsize='medium')
        txt.set_path_effects([
            patheffects.withStroke(linewidth=5, foreground='w')
        ])

    # configure and draw precision percentages
    axActual.set_xticks([0])
    axActual.set_yticks(ticks)
    axActual.set_ylabel('Precision', fontsize='large')
    axActual.yaxis.set_label_position('right')
    axActual.imshow(
        correctpredicted,
        interpolation='nearest',
        cmap='RdYlGn',
        aspect='auto',
        vmin=0.0,
        vmax=1.0
    )
    for i in range(len(class_names)):
        txt = axActual.text(
            0, i,
            ('100' if correctpredicted[i, 0] == 1.0
                else f'{100.0 * correctpredicted[i, 0]:3.1f}'),
            ha='center',
            va='center',
            color='black',
            fontsize='medium')
        txt.set_path_effects([
            patheffects.withStroke(linewidth=5, foreground='w')
        ])

    # configure and draw total accuracy
    axTotal.set_xticks([0])
    axTotal.set_yticks([0])
    axTotal.set_xlabel('Accuracy', fontsize='large')
    axTotal.imshow(
        np.array([[accuracy]]),
        interpolation='nearest',
        cmap='RdYlGn',
        aspect='auto',
        vmin=0.0,
        vmax=1.0
    )
    txt = axTotal.text(
        0, 0,
        f'{100 * accuracy:3.1f}',
        ha='center',
        va='center',
        color='black',
        fontsize='medium'
    )
    txt.set_path_effects([
        patheffects.withStroke(linewidth=5, foreground='w')
    ])

    # disable axes for other matrices than confusion matrix
    for a in (axPredicted, axActual, axTotal):
        plt.setp(a.get_yticklabels(), visible=False)
        plt.setp(a.get_xticklabels(), visible=False)

    # draw colorbar for confusion matrix
    cbar = fig.colorbar(
        img,
        ax=[axPredicted, axConfMatrix, axActual, axTotal],
        shrink=0.5,
        pad=0.1
    )
    for t in cbar.ax.get_yticklabels():
        t.set_fontsize('medium')
    suptitlehandle = fig.suptitle(
        f'{title} (ACC={accuracy:.5f})',
        fontsize='xx-large'
    )
    if outpath is None:
        plt.show()
    else:
        plt.savefig(
            outpath,
            dpi=dpi,
            bbox_inches='tight',
            bbox_extra_artists=[suptitlehandle],
            pad_inches=0.1
        )

def recall_precision_curves(
        outpath: Optional[Path],
        title: str,
        lines: List[Tuple[List, List]],
        class_names: List[str],
        figsize: Tuple = (15, 15)):
    """
    Draws Recall-Precision curves for AP measurements.

    Parameters
    ----------
    outpath : Optional[Path]
        Output path for the plot image. If None, the plot will be displayed.
    title : str
        Title of the plot
    lines : List[List[List]]
        Per-class list of tuples with list of recall values and precision
        values
    class_names : List[str]
        List of the class names
    figsize: Tuple
        The size of the figure
    """
    fig = plt.figure(figsize=figsize)
    ax = fig.add_subplot(111)
    colormap = plt.cm.nipy_spectral
    colors = [colormap(i) for i in np.linspace(0, 1, len(class_names))]
    linestyles = ['-', '--', '-.', ':']
    for i, (cls, line) in enumerate(zip(class_names, lines)):
        ax.plot(
            line[0], line[1],
            label=cls, c=colors[i], linewidth=3,
            linestyle=linestyles[i % len(linestyles)],
            alpha=0.8
        )
    legendhandle = ax.legend(
        bbox_to_anchor=(0.5, -0.3),
        loc='lower center',
        ncol=10
    )
    ax.set_xlabel('recall')
    ax.set_ylabel('precision')
    ax.set_xlim((0.0, 1.01))
    ax.set_ylim((0.0, 1.01))
    ax.grid('on')
    ax.set_xticks(np.arange(0, 1.1, 0.1))
    ax.set_yticks(np.arange(0, 1.1, 0.1))
    ax.set_title(title)

    if outpath is None:
        plt.show()
    else:
        fig.savefig(
            outpath,
            bbox_extra_artists=[legendhandle],
            bbox_inches='tight'
        )

def true_positive_iou_histogram(
        outpath: Optional[Path],
        title: str,
        lines: List[float],
        class_names: List[str],
        figsize: Tuple = (10, 25)):
    """
    Draws per-class True Positive IoU precision plot

    Parameters
    ----------
    outpath : Optional[Path]
        Output path for the plot image. If None, the plot will be displayed.
    title : str
        Title of the plot
    lines : List[float]
        Per-class list of floats with IoU values
    class_names : List[str]
        List of the class names
    figsize: Tuple
        The size of the figure
    """
    plt.figure(figsize=figsize)
    plt.barh(
        class_names,
        np.array(lines),
        orientation='horizontal',
        color='purple'
    )
    plt.ylim((-1, len(class_names)))
    plt.yticks(np.arange(0, len(class_names)))
    plt.xticks(np.arange(0, 1.1, 0.1))
    plt.xlabel('IoU precision')
    plt.ylabel('classes')
    plt.title(f'{title}')
    plt.tight_layout()

    if outpath is None:
        plt.show()
    else:
        plt.savefig(outpath)

def true_positives_per_iou_range_histogram(
        outpath: Optional[Path],
        title: str,
        lines: List[float],
        range_fraction: float = 0.05,
        figsize: Tuple = (10, 10)):
    """
    Draws histogram of True Positive IoU values

    Parameters
    ----------
    outpath : Optional[Path]
        Output path for the plot image. If None, the plot will be displayed.
    title : str
        Title of the plot
    lines : List[float]
        All True Positive IoU values
    range_fraction : float
        Fraction by which the range should be divided (1/number_of_segments)
    figsize: Tuple
        The size of the figure
    """
    lower_bound = floor(min(lines)*10)/10
    x_range = np.arange(lower_bound, 1.01, (1-lower_bound)*range_fraction)
    plt.figure(figsize=figsize)
    plt.hist(
        lines,
        x_range,
        color='purple'
    )
    plt.xlabel('IoU ranges')
    plt.xticks(x_range, rotation=45)
    plt.ylabel('Number of masks in IoU range')
    plt.title(f'{title}')
    plt.tight_layout()

    if outpath is None:
        plt.show()
    else:
        plt.savefig(outpath)

def recall_precision_gradients(
        outpath: Optional[Path],
        title: str,
        lines: List[Tuple[List, List]],
        class_names: List[str],
        aps: List[float],
        figsize: Tuple = (10, 25)):
    """
    Draws per-class gradients of precision dependent to recall.

    Provides per-class AP and mAP values.

    Parameters
    ----------
    outpath : Optional[Path]
        Output path for the plot image. If None, the plot will be displayed.
    title : str
        Title of the plot
    lines : List[Tuple[List, List]]
        Per-class list of tuples with list of recall values and precision
        values
    class_names : List[str]
        List of the class names
    aps: List[float]
        Per-class AP values
    figsize: Tuple
        The size of the figure
    """
    plt.figure(figsize=figsize)
    clsticks = []
    for i, (cls, line, averageprecision) \
            in enumerate(zip(class_names, lines, aps)):
        clscoords = np.ones(len(line[0])) * i
        points = np.array([line[0], clscoords]).T.reshape(-1, 1, 2)
        segments = np.concatenate([points[:-1], points[1:]], axis=1)
        lc = LineCollection(
            segments,
            cmap='RdYlGn',
            norm=plt.Normalize(0, 1.0)
        )
        lc.set_array(line[1])
        lc.set_linewidth(10)
        plt.gca().add_collection(lc)
        clsticks.append(f'{cls} (AP={averageprecision})')
    plt.ylim((-1, len(class_names)))
    plt.yticks(np.arange(0, len(clsticks)), labels=clsticks)
    plt.xticks(np.arange(0, 1.1, 0.1))
    plt.xlabel('recall')
    plt.ylabel('classes')
    plt.colorbar(
        plt.cm.ScalarMappable(norm=plt.Normalize(0, 1.0), cmap='RdYlGn'),
        orientation='horizontal',
        label='precision',
        fraction=0.1,
        pad=0.05
    )
    plt.title(f'{title} (mAP={np.mean(aps)})')
    plt.tight_layout()

    if outpath is None:
        plt.show()
    else:
        plt.savefig(outpath)

def draw_plot(
        outpath: Optional[Path],
        title: str,
        xtitle: str,
        xunit: str,
        ytitle: str,
        yunit: str,
        line: Tuple[List, List],
        figsize: Tuple = (15, 15)):
    """
    Draws plot.

    Parameters
    ----------
    outpath : Optional[Path]
        Output path for the plot image. If None, the plot will be displayed.
    title : str
        Title of the plot
    xtitle : str
        Name of the X axis
    xuint : str
        Unit for the X axis
    ytitle : str
        Name of the Y axis
    yunit : str
        Unit for the Y axis
    line : Tuple[List, List]
        Per-class list of tuples with list of recall values and precision
        values
    figsize: Tuple
        The size of the figure
    """
    plt.figure(figsize=figsize)
    plt.plot(line[0], line[1], c='purple', linewidth=3)
    xlabel = xtitle
    if xunit is not None:
        xlabel += f' [{xunit}]'
    ylabel = ytitle
    if yunit is not None:
        ylabel += f' [{yunit}]'
    plt.xlabel(xlabel, fontsize='large')
    plt.ylabel(ylabel, fontsize='large')
    plt.grid()
    plt.title(title)

    if outpath is None:
        plt.show()
    else:
        plt.savefig(outpath)

max_stars_repo_name: antmicro/edge-ai-tester

max_issues_repo_name: antmicro/edge-ai-tester

max_forks_repo_name: antmicro/edge-ai-tester


text:
"""
Generate a PDF of a given quantity also accounting for boostrap errors

"""
from __future__ import print_function, absolute_import, division, unicode_literals

#some global staff first
import glob
import pickle
import matplotlib.pyplot as plt
import matplotlib
import numpy as np
matplotlib.rcParams.update({'font.size': 16})
import os.path
import h5py

from astropy.table import Table

from xastropy.xutils import xdebug as xdb

mblue='#1E90FF'
mred='#DC143C'

#now the work function
def work_generate_pdf(modelsdir, quantity, data_fil, outh5):

    print('Processing ', modelsdir, quantity)

    if(quantity == 'met'):
        minval=-5
        maxval=1.5
        binsize=0.15
        precomp=False
    elif(quantity == 'dens'):
        minval=-5
        maxval=1
        binsize=0.15
        precomp=False
    elif(quantity == 'jgal'):
        minval=-5.4
        maxval=2.4
        binsize=0.15
        precomp=False
    elif(quantity == 'jqso'):
        minval=-6.8
        maxval=1.0
        binsize=0.15
        precomp=False
    elif(quantity == 'fstar'):
        minval=-2
        maxval=2
        binsize=0.15
        precomp=False
    elif(quantity == 'temp'):
        minval=3.6
        maxval=6.4
        binsize=0.15
        precomp=False
    elif(quantity == 'upar'):
        precomp=True
        filein=modelsdir+'/'+modelsdir+'_ionisation.pkl'
        if(os.path.isfile(filein)):
            try:
                pdfcom=pickle.load(open(filein))
            except:
                pdfcom=pickle.load(open(filein,'rb'),encoding='latin1')
        else:
            print('Need to compute {} first'.format(filein))
            exit()
        #load precomputed histograms
        binsize=pdfcom['U_base'][1]-pdfcom['U_base'][0]
        minval=np.min(pdfcom['U_base'])
        maxval=np.max(pdfcom['U_base'])+binsize
        base_hist=pdfcom['U_base']
        pdf_hist=pdfcom['U_hist']
    elif(quantity == 'xhi'):
        precomp=True
        filein=modelsdir+'/'+modelsdir+'_ionisation.pkl'
        if(os.path.isfile(filein)):
            try:
                pdfcom=pickle.load(open(filein))
            except:
                pdfcom=pickle.load(open(filein,'rb'),encoding='latin1')
        else:
            print('Need to compute {} first'.format(filein))
            exit()
        #load precomputed histograms
        binsize=pdfcom['XHI_base'][1]-pdfcom['XHI_base'][0]
        minval=np.min(pdfcom['XHI_base'])
        maxval=np.max(pdfcom['XHI_base'])+binsize
        base_hist=pdfcom['XHI_base']
        pdf_hist=pdfcom['XHI_hist']
    elif(quantity == 'si2'):
        precomp=True
        filein=modelsdir+'/'+modelsdir+'_ionisation.pkl'
        if(os.path.isfile(filein)):
            try:
                pdfcom=pickle.load(open(filein))
            except:
                pdfcom=pickle.load(open(filein,'rb'),encoding='latin1')
        else:
            print('Need to compute {} first'.format(filein))
            exit()
        #load precomputed histograms
        binsize=pdfcom['XSiII_base'][1]-pdfcom['XSiII_base'][0]
        minval=np.min(pdfcom['XSiII_base'])
        maxval=np.max(pdfcom['XSiII_base'])+binsize
        base_hist=pdfcom['XSiII_base']
        pdf_hist=pdfcom['XSiII_hist']
    else:
        print('Option not found!')
        exit()

    #some util
    nsample=1000
    figname=modelsdir

    #load all data
    info = Table.read(data_fil, format='ascii')
    #info=np.loadtxt(data_fil,dtype={'names':('name','zabs','ezabs','ion','logN','elogN','flag','sample'),
    #                           'formats':('S30', 'f8', 'f8','S4','f6','f6','i4','S5')})

    #if quantity is not precomputed, compute it
    if(precomp == False):
        #load the models
        listin=glob.glob(modelsdir+'/*emcee*hd5')
        nsys=len(listin)
        print('Found {} fits'.format(nsys))

        #This defines the bin edges, including the rightmost edge
        base_hist=np.arange(minval,maxval,binsize)
        pdf_hist=[]
        sample_tag=[]
        outh5[quantity]['left_edge_bins'] = base_hist[0:-1]
        outh5[quantity]['right_edge_bins'] = base_hist[1:]
        outh5[quantity]['left_edge_bins'].attrs['BINSIZE'] = binsize

        #first, loop over the sample to construct PDFs for individual chains
        print('Generating PDFs for individual systems')

        for fit in listin:

            # Load
            fh5=h5py.File(fit, 'r')
            #xdb.set_trace()
            #fl=open(fit)
            #result=pickle.load(fl)
            #fl.close()

            #find sample tag
            #match=np.where(info['name'] == result['info']['name'])
            match=np.where(info['name'] == fh5['inputs'].attrs['name'])
            sample_tag.append(info['sample'][match])

            #populate histo
            #qnt=result['tags'].index(quantity)
            qnt= np.where(fh5['outputs']['tags'].value == quantity)[0][0]
            #hist,edge=np.histogram(result['pdfs'][:,qnt],bins=base_hist)
            hist,edge=np.histogram(fh5['outputs']['pdfs'][:,qnt],bins=base_hist)

            #normalise
            hist=hist/binsize/len(fh5['outputs']['pdfs'][:,qnt])

            #Append
            pdf_hist.append(hist)

            # Write to hdf5
            outh5[quantity][fh5['inputs'].attrs['name']] = hist
            # Inputs
            try:
                outh5['inputs'].create_group(fh5['inputs'].attrs['name'])
            except ValueError:  # No repeat
                pass
            else:
                for key in fh5['inputs']:
                    outh5['inputs'][fh5['inputs'].attrs['name']][key] = fh5['inputs'][key].value
    else:
        print("Work with precomputed quantity")
        nsys=len(pdfcom['name'])
        print('Found {} fits'.format(nsys))

        #Now need to find the sample tag
        sample_tag=[]

        for fit in pdfcom['name']:

            #find sample tag
            match=np.where(info['name'] == fit)
            sample_tag.append(info['sample'][match][0])

    #now loop over trials for entire sample
    pdf_trials=np.zeros((nsample,base_hist.size-1))

    print('Bootstrap all sample')

    for trial in range(nsample):

        #now generate random indexes with repetition
        thisindex=np.random.random_integers(0,high=nsys-1,size=nsys)

        #fill the pdf
        pdf_thisset=np.zeros(base_hist.size-1)

        for ii in thisindex:
            pdf_thisset=pdf_thisset+pdf_hist[ii]

        #normalise and stack
        pdf_thisset=pdf_thisset/nsys
        pdf_trials[trial,:]=pdf_thisset

    #now compute percentiles
    pdf_med=np.percentile(pdf_trials,50,axis=0)
    pdf_90=np.percentile(pdf_trials,90,axis=0)
    pdf_75=np.percentile(pdf_trials,75,axis=0)
    pdf_25=np.percentile(pdf_trials,25,axis=0)
    pdf_10=np.percentile(pdf_trials,10,axis=0)

    #plot
    plt.plot(base_hist[0:-1],pdf_med)
    plt.plot(base_hist[0:-1],pdf_90)
    plt.plot(base_hist[0:-1],pdf_10)
    plt.show()

    #now extract the statistical subset
    sample_tag=np.array(sample_tag)
    substat=np.where(sample_tag == 'P15')
    substat=substat[0]
    nstat=len(substat)
    pdf_hist_sub=[]

    for ii in substat:
        pdf_hist_sub.append(pdf_hist[ii])

    #now repeat sampling for the statistical subset
    pdf_trials=np.zeros((nsample,base_hist.size-1))

    return
    #####################################################################
    #####################################################################
    #####################################################################

    boot_stat = False
    if boot_stat:
        print('Bootstrap stat sample')

        for trial in range(nsample):

            #now generate random indexes with repetition
            try:
                thisindex=np.random.random_integers(0,high=nstat-1,size=nstat)
            except ValueError:
                xdb.set_trace()

            #fill the pdf
            pdf_thisset=np.zeros(base_hist.size-1)

            for ii in thisindex:
                pdf_thisset=pdf_thisset+pdf_hist_sub[ii]

            #normalise to nsys (not nstat to have fraction) and stack
            pdf_thisset=pdf_thisset/nsys
            pdf_trials[trial,:]=pdf_thisset

        #now compute percentiles
        pdf_med_stat=np.percentile(pdf_trials,50,axis=0)
        pdf_90_stat=np.percentile(pdf_trials,90,axis=0)
        pdf_75_stat=np.percentile(pdf_trials,75,axis=0)
        pdf_25_stat=np.percentile(pdf_trials,25,axis=0)
        pdf_10_stat=np.percentile(pdf_trials,10,axis=0)

        #write
        nameout="{0}/{1}_{2}_pdfdata.pkl".format(modelsdir,figname,quantity)
        save={'base':base_hist,'nsample':nsample,'minval':minval,'maxval':maxval,'binsize':binsize,'quantity':quantity,
              'pdf50':pdf_med,'pdf90':pdf_90,'pdf75':pdf_75,'pdf25':pdf_25,'pdf10':pdf_10,'nsys':nsys,'nstat':nstat,
              'pdf50stat':pdf_med_stat,'pdf90stat':pdf_90_stat,'pdf75stat':pdf_75_stat,'pdf25stat':pdf_25_stat,'pdf10stat':pdf_10_stat}
        pickle.dump(save,open(nameout,'w'))

        #plot
        plt.plot(base_hist[0:-1],pdf_med_stat)
        plt.plot(base_hist[0:-1],pdf_90_stat)
        plt.plot(base_hist[0:-1],pdf_10_stat)
        plt.show()

if __name__ == "__main__":

    #start with general quantities
    workers=[]
    #allmod=['minimalnhi','sourcenhi','cafgnhi','dustnhi']
    #allmod=['MCMC_TST']
    allmod=['./Lehner13-MCMC-CRC/MCMC_FULL/']
    alltype=['met', 'dens']
    mcmc_fil = './Lehner13-MCMC-CRC/MgIIsurvey_MTL.ascii'
    hdf5_outfil = './MgIIsurvey_MTL.hdf5'

    #allmod=['dustnhi']
    #alltype=['fstar']

    #allmod=['dustnhi']
    #alltype=['upar','xhi','si2']
    #alltype=['xhi']

    outh5 = h5py.File(hdf5_outfil, 'w')
    outh5.create_group('inputs')

    for mm in allmod:
        for tt in alltype:
            outh5.create_group(tt)
            work_generate_pdf(mm, tt, mcmc_fil, outh5)

    print("Writing {:s}".format(hdf5_outfil))
    outh5.close()

max_stars_repo_name: pyigm/pyigm

max_issues_repo_name: pyigm/pyigm

max_forks_repo_name: pyigm/pyigm

reason: import numpy,from astropy