text:
# *************************************************************************************************
# Modified from rl/example/dqn_atari.py in keras-rl (https://github.com/matthiasplappert/keras-rl)
# by Heejin Chloe Jeong ( University of Pennsylvania, chloe.hjeong@gmail.com)
# @ Jan 2018
# *************************************************************************************************

from __future__ import division
import argparse

from PIL import Image
import numpy as np
import gym
from gym.wrappers import Monitor

import os, pdb, pickle, datetime, json
import tensorflow as tf

from keras.models import Sequential
from keras.layers import Dense, Activation, Flatten, Convolution2D, Permute
from keras.optimizers import Adam, RMSprop
import keras.backend as K
import keras.initializers as KInit
from keras.utils.generic_utils import get_custom_objects
from rl.policy import LinearAnnealedPolicy, BoltzmannQPolicy, EpsGreedyQPolicy
from rl.memory import SequentialMemory
from rl.core import Processor
from rl.callbacks import FileLogger, ModelIntervalCheckpoint

INPUT_SHAPE = (84, 84)
WINDOW_LENGTH = 4

class AtariProcessor(Processor):
    def process_observation(self, observation):
        assert observation.ndim == 3  # (height, width, channel)
        img = Image.fromarray(observation)
        img = img.resize(INPUT_SHAPE).convert('L')  # resize and convert to grayscale
        processed_observation = np.array(img)
        assert processed_observation.shape == INPUT_SHAPE
        return processed_observation.astype('uint8')  # saves storage in experience memory

    def process_state_batch(self, batch):
        # We could perform this processing step in `process_observation`. In this case, however,
        # we would need to store a `float32` array instead, which is 4x more memory intensive than
        # an `uint8` array. This matters if we store 1M observations.
        processed_batch = batch.astype('float32') / 255.
        return processed_batch

    def process_reward(self, reward):
        return np.clip(reward, -1., 1.)

# Default are from Hasselt's Double DQN paper.
parser = argparse.ArgumentParser()
parser.add_argument('--mode', choices=['train', 'test'], default='train')
parser.add_argument('--env_name', type=str, default='Breakout-v0')
parser.add_argument('--weights', type=str, default=None)
parser.add_argument('--double_dqn', type=int, default=0)
parser.add_argument('--dueling', type=int, default=0)
parser.add_argument('--num_layers', type=int, default=1)
parser.add_argument('--num_units', type=int, default=512)
parser.add_argument('--batch_size', type=int, default=32)
parser.add_argument('--init_step', type=int, default=0) # When you want to start from learned weights.
parser.add_argument('--nb_train_steps', type=int, default = 50000000)
parser.add_argument('--nb_steps_warmup', type=int, default = 50000)
parser.add_argument('--no_ops', type=int, default=0)
parser.add_argument('--nb_max_start_steps', type = int, default=0)
parser.add_argument('--epoch_steps', type=int, default = 50000)
parser.add_argument('--target_model_update', type=int, default=10000)
parser.add_argument('--learning_rate', type=float, default=0.00025)
parser.add_argument('--gamma', type=float, default=.99)
parser.add_argument('--log_dir', type=str, default='.')
parser.add_argument('--memory_size', type=int, default=500000)
parser.add_argument('--eps_max', type=float, default=1.)
parser.add_argument('--eps_min', type=float, default=.1)
parser.add_argument('--init_mean', type =float, default=1.)
parser.add_argument('--init_sd', type=float, default=50.)
parser.add_argument('--sd_max', type=float, default=10.)
parser.add_argument('--sd_min', type=float, default=0.01)
parser.add_argument('--sd_steps', type=int, default=50000000)
parser.add_argument('--gpu_memory', type=float, default=0.1)
parser.add_argument('--loss_type', type=str, default='kl')
parser.add_argument('--device', type=str, default='/cpu:0')
parser.add_argument('--alg', choices=['dqn','adfq'], default='dqn')
parser.add_argument('--record',type=int, default=0)

args = parser.parse_args()

# Get the environment and extract the number of actions.
env = gym.make(args.env_name)
if args.record == 1:
    env = Monitor(env, directory=args.log_dir)
np.random.seed(123)
env.seed(123)
nb_actions = env.action_space.n

# Next, we build our model. We use the same model that was described by Mnih et al. (2015).
input_shape = (WINDOW_LENGTH,) + INPUT_SHAPE

with tf.device(args.device):
    model = Sequential()
    if K.image_dim_ordering() == 'tf':
        # (width, height, channels)
        model.add(Permute((2, 3, 1), input_shape=input_shape))
    elif K.image_dim_ordering() == 'th':
        # (channels, width, height)
        model.add(Permute((1, 2, 3), input_shape=input_shape))
    else:
        raise RuntimeError('Unknown image_dim_ordering.')
    model.add(Convolution2D(32, 8, 8, subsample=(4, 4)))
    model.add(Activation('relu'))
    model.add(Convolution2D(64, 4, 4, subsample=(2, 2)))
    model.add(Activation('relu'))
    model.add(Convolution2D(64, 3, 3, subsample=(1, 1)))
    model.add(Activation('relu'))
    model.add(Flatten())
    for _ in range(args.num_layers):
        model.add(Dense(args.num_units))
        model.add(Activation('relu'))

    if args.alg=='dqn':
        from BRL.DRL.keras.dqn import DQNAgent
        agent_class = DQNAgent
        model.add(Dense(nb_actions))

    elif args.alg== 'adfq':
        from BRL.DRL.keras.adfq import ADFQAgent
        agent_class = ADFQAgent
        def custom_initializer(shape, dtype=None):
            v = [args.init_mean]*int(shape[0]/2) + [args.init_sd]*int(shape[0]/2)
            return K.constant(v, dtype=dtype)
        class CustomInitializer:
            def __call__(self, shape, dtype=None):
                return custom_initializer(shape, dtype=dtype)
        get_custom_objects().update({'custom_initializer': CustomInitializer})
        model.add(Dense(nb_actions*2, bias_initializer =  custom_initializer)) # mean and SD

    model.add(Activation('linear'))
    print(model.summary())

    memory = SequentialMemory(limit=args.memory_size, window_length=WINDOW_LENGTH)
    processor = AtariProcessor()
    policy = LinearAnnealedPolicy(EpsGreedyQPolicy(), attr='eps', value_max=1., value_min=.1,
                            value_test=.05, nb_steps=1000000)
    test_policy = EpsGreedyQPolicy(eps=0.05)

    if bool(args.double_dqn):
        print("DOUBLE DQN")
    if bool(args.dueling):
        print("DUELING NETWORK")

    if 'GPU' in args.device.split(':'):
        config = tf.ConfigProto(allow_soft_placement=True)
        config.gpu_options.per_process_gpu_memory_fraction = args.gpu_memory
        config.gpu_options.polling_inactive_delay_msecs = 25
        session = tf.Session(config=config)
        K.set_session(session)

    if args.alg== 'adfq':
        agent = agent_class(model=model, nb_actions=nb_actions, policy=policy, test_policy=test_policy,
                    memory=memory, processor=processor, nb_steps_warmup=args.nb_steps_warmup,
                    sd_min = args.sd_min, sd_max = args.sd_max, sd_steps = args.sd_steps,
                    gamma=args.gamma, batch_size=args.batch_size, target_model_update=args.target_model_update,
                    enable_double_dqn = bool(args.double_dqn), enable_dueling_network = bool(args.dueling),
                    train_interval=4, delta_clip=1., no_ops = args.no_ops)
        agent.compile(Adam(lr=args.learning_rate, clipnorm=10.0), loss_type = args.loss_type, metrics=[])
    elif args.alg == 'dqn':
        agent = agent_class(model=model, nb_actions=nb_actions, policy=policy, test_policy=test_policy,
                    memory=memory, processor=processor, nb_steps_warmup=args.nb_steps_warmup,
                    gamma=args.gamma, batch_size=args.batch_size, target_model_update=args.target_model_update,
                    enable_double_dqn = bool(args.double_dqn), enable_dueling_network = bool(args.dueling),
                    train_interval=4, delta_clip=1., no_ops = args.no_ops)
        #agent.compile(Adam(lr=args.learning_rate, clipnorm = 10.0), metrics=['mae'])
        agent.compile(RMSprop(lr=args.learning_rate, clipnorm = 10.0, rho = 0.95), metrics=['mae'])
    directory = os.path.join(args.log_dir, datetime.datetime.now().strftime("%m%d%H%M"))
    if args.mode == 'train':
        if not os.path.exists(directory):
            os.makedirs(directory)
        else:
            ValueError("The directory already exists...", directory)
        json.dump(vars(args), open(os.path.join(directory, 'learning_prop.json'), 'wb'))

        if args.weights:
            agent.load_weights(args.weights)
            agent.nb_steps_warmup += args.init_step

        weights_filename = os.path.join(directory, (args.alg + '_{}_weights.h5f').format(args.env_name))
        checkpoint_weights_filename = os.path.join(directory, args.alg + '_' + args.env_name + '_weights_{step}.h5f')
        log_filename = os.path.join(directory, (args.alg + '_{}_log.json').format(args.env_name))
        callbacks = [ModelIntervalCheckpoint(checkpoint_weights_filename, interval=1000000)]
        callbacks += [FileLogger(log_filename, interval=args.epoch_steps/2)] # interval is for episode intervals. # of points nb_train_steps/(interval*1000)
        _ = agent.fit(env, callbacks=callbacks, nb_steps=args.nb_train_steps,
            log_interval=args.epoch_steps, nb_steps_epoch = args.epoch_steps,
            nb_max_start_steps=args.nb_max_start_steps, init_step = args.init_step, epoch_file_path = os.path.join(directory, args.alg + '_epoch_test.pkl'))

        # After training is done, we save the final weights one more time.
        agent.save_weights(weights_filename, overwrite=True)

        # Finally, evaluate our algorithm for 10 episodes.
        agent.test(env, nb_episodes=10, visualize=False)

    elif args.mode == 'test':
        weights_filename = os.path.join(directory,(args.alg + '_{}_weights.h5f').format(args.env_name))
        if args.weights:
            weights_filename = args.weights
        agent.load_weights(weights_filename)
        agent.test(env, nb_episodes=10, visualize=True)

max_stars_repo_name: christopher-hsu/ray

max_issues_repo_name: christopher-hsu/ray

max_forks_repo_name: christopher-hsu/ray


text:
import pandas as pd
import re
import numpy as np
from collections import namedtuple
from typing import List

MergeFrame = namedtuple('MergeFrame', ['frame', 'left_on', 'right_on'])

def field_format(field: str, options: List[str], index: int) -> str:
    regex = '.{' + str(index) + '}.*'
    match = re.match(regex, field)
    if match is None:
        return field
    choice = np.random.choice(options)
    return field.format(choice)

def merge_frames(df: pd.DataFrame, to_merge: List[MergeFrame], options: List[List[str]] = []):
    result = df
    for frame in to_merge:
        left_on = frame.left_on
        for i, o in enumerate(options):
            left_on = field_format(left_on, o, i)

        right_on = frame.right_on
        for i, o in enumerate(options):
            right_on = field_format(right_on, o, i)

        result = result.merge(right=frame.frame, left_on=left_on, right_on=right_on, how='inner')
    return result

ext: py

reason: import numpy


text:
import numpy as np
import torch

from envs import ode_env
import disturb_models as dm
from constants import *

'''
Adapted from:
Quang Linh Lam, Antoneta Iuliana Bratcu, Delphine Riu. Frequency robust control in stand-alone
microgrids with PV sources : design and sensitivity analysis. Symposium de Génie Electrique, Jun
2016, Grenoble, France. ffhal-01361556
'''

class MicrogridEnv(ode_env.NLDIEnv):

    def __init__(self, random_seed=None, device=None):
        if random_seed is not None:
            np.random.seed(random_seed)
            torch.manual_seed(random_seed+1)

        self.n = 3
        self.m = 2
        self.wp = 1
        self.wq = 1
        self.isD0 = True   # TODO?

        # for per-unit normalization
        v_base = 585
        p_base = 1.0015
        i_base = p_base / v_base
        r_base = v_base / i_base

        omega_b = 3.5  # between 2.61 and 5.22 given in paper
        R_dc = 50 / r_base   # based on 5-200 ohms range at https://doc.ingeniamc.com/venus/product-manual/installation-and-configuration/motor-output-wiring/shunt-braking-resistor
        C_dc = 500e-6 * r_base     # based on 200-800 micro-farads given here: https://www.elmomc.com/wp-content/uploads/2019/08/Simple-Capacitor-white-paper.pdf
                        # and also approx here: https://doc.ingeniamc.com/titan/manuals/titan-go-product-manual/wiring-and-connections/dc-bus-bulk-capacitance
        T_diesel = 0.01  # based on plausible value from https://core.ac.uk/download/pdf/52115363.pdf (TODO switch to p.u.?)
        s_diesel = 0.04  # arbitrary, based on plausible value from Wikipedia: https://en.wikipedia.org/wiki/Droop_speed_control
        H = 0.7305  # based on plausible value from https://core.ac.uk/download/pdf/52115363.pdf (TODO switch to p.u.?)
        D_load = 0.9   # arbitrary
        alpha_ce = 0.585
        beta_de = 0.4
        v_sce = 585 / v_base
        R_sc = 30e-3 / r_base    # based on plausible internal supercapacitor resistence of 30 milli-ohms from https://en.wikipedia.org/wiki/Supercapacitor#Internal_resistance
        i_se = -2.5 / i_base

        self.A = 0.001 * torch.tensor((
            (-omega_b/(R_dc*C_dc),      0,      0),
            (0,     -1/T_diesel,    -1/(T_diesel * s_diesel)),
            (0,     1/(2*H),        -D_load/(2*H))
        ), dtype=TORCH_DTYPE, device=device)

        self.B = torch.tensor((
            (omega_b * alpha_ce / C_dc, -omega_b * beta_de / C_dc),
            (0,     0),
            ( (v_sce - 2*R_sc*i_se)/(2*H),  0)
        ), dtype=TORCH_DTYPE, device=device)

        self.G = torch.tensor((
                (0, ),
                (0, ),
                (-1/(2*H), )
        ), dtype=TORCH_DTYPE, device=device)

        # Capture some (arbitrary) dependence between voltage/freq variation and disturb
        #    Note: Doesn't always solve, but solves for random seed 10
        # self.C = torch.tensor((-0.05, 0, 0.05), dtype=TORCH_DTYPE, device=device).unsqueeze(0)
        self.C = 5.0 * torch.tensor(np.random.randn(3), dtype=TORCH_DTYPE, device=device).unsqueeze(0)
        self.D = torch.zeros(self.wq, self.m, dtype=TORCH_DTYPE, device=device)

        # Objective: Assign weight 1 to entries in output vector y,
        #    and small weight to other values
        self.Q = torch.tensor((
            (1,  0,    0),
            (0,  0.1,  0),
            (0,  0,    1)
        ), dtype=TORCH_DTYPE, device=device)

        self.R = torch.tensor((
            (0.1,   0),
            (0,     0.1)
        ), dtype=TORCH_DTYPE, device=device)

        self.disturb_f = dm.NLDIDisturbModel(self.C, self.D, self.n, self.m, self.wp)
        if device is not None:
            self.disturb_f.to(device=device, dtype=TORCH_DTYPE)

        self.adversarial_disturb_f = None

    def xdot_f(self, x, u, t):
        w = self.disturb_f(x, u, t)
        return x @ self.A.T + u @ self.B.T + w @ self.G.T

    def xdot_adversarial_f(self, x, u, t):
        if self.adversarial_disturb_f is None:
            raise ValueError('You must initialize adversarial_disturb_f before running in adversarial mode')
        w = self.adversarial_disturb_f(x, u, t)
        return x @ self.A.T + u @ self.B.T + w @ self.G.T

    def cost_f(self, x, u, t):
        return ((x @ self.Q) * x).sum(-1) + ((u @ self.R) * u).sum(-1)

    def get_nldi_linearization(self):
        return self.A, self.B, self.G, self.C, self.D, self.Q, self.R

    def gen_states(self, num_states, device=None):
        return torch.tensor(np.random.rand(num_states, self.n), device=device, dtype=TORCH_DTYPE)

    def __copy__(self):
        new_env = MicrogridEnv.__new__(MicrogridEnv)
        new_env.__dict__.update(self.__dict__)
        return new_env

max_stars_repo_name: locuslab/robust-nn-control

max_issues_repo_name: locuslab/robust-nn-control

max_forks_repo_name: locuslab/robust-nn-control


text:
from base64 import b64decode, b64encode

import numpy as np
import cv2

def base64_to_bw_img(blob):
    """Convert base64 to black and white image."""
    data = b64decode(blob)
    data = np.frombuffer(data, np.uint8)
    img = cv2.imdecode(data, 0)
    return img

def img_to_base64(img):
    retval, buffer = cv2.imencode('.png', img)
    return b64encode(buffer)


text:
# coding: utf-8
#
#def ceti_exp(
#    GHZ_inner, GHZ_outer, tau_awakening, tau_survive, D_max, t_max
#):
##{{{
#
#    import ceti_tools as ct
#
#    random.seed(420)
#    np.random.seed(420)
#
#    # lista de CETIs alguna vez activas
#    CETIs = dict()
#
#    # lista de CETIs actualmente activas
#    CHATs = []
#    CHATs_idx = []
#
#    # inicializacion del tiempo: scalar
#    t_now = 0
#    # inicializacion del ID: index
#    ID = 0
#
#    # lista de tiempos de eventos futuros: ordered list
#    # [time, ID_emit, ID_receive, case]
#    t_forthcoming = ct.OrderedList()
#
#    # estructura de arbol para buscar vecinos
#    try:
#        del tree
#    except NameError:
#        pass
#
#    # INITIALIZATION
#    # Simulation starts when the first CETI appears:
#    next_event = [0., 0, None, 1]
#    t_forthcoming.add(next_event)
#
#    # SIMULATION LOOP OVER TIME
#    while (t_now < t_max):
#
#        t_now, ID_emit, ID_hear, case = t_forthcoming.head.getData()
#
#        if case == 1:
#
#            # print('desaparece CETI con id:%d' % ID_emit)
#            ID_new = ID_emit
#            ID_next = ID_new + 1
#            t_new_hola = t_now
#
#            # sortear el lugar donde aparece dentro de la GHZ
#            r = np.sqrt(random.random()*GHZ_outer**2 + GHZ_inner**2)
#            o = random.random()*2.*np.pi
#            x = r * np.cos(o)  # X position on the galactic plane
#            y = r * np.sin(o)  # Y position on the galactic plane
#
#            # sortear el tiempo de actividad
#            t_new_active = np.random.exponential(tau_survive, 1)[0]
#            t_new_chau = t_new_hola + t_new_active
#
#            # agregar el tiempo de desaparición a la lista de tiempos
#            next_event = [t_new_chau, ID_new, None, 2]
#            t_forthcoming.add(next_event)
#
#            # agregar la CETI a la lista histórica
#            CETIs[ID_new] = list()
#            CETIs[ID_new].append(
#                (ID_new, ID_new, x, y, t_new_hola, t_new_chau))
#
#            # sortear el tiempo de aparición de la próxima CETI
#            t_next_awakening = np.random.exponential(tau_awakening, 1)[0]
#            t_next_awakening = t_new_hola + t_next_awakening
#            if t_next_awakening < t_max:
#                next_event = [t_next_awakening, ID_next, None, 1]
#                t_forthcoming.add(next_event)
#
#            if len(CHATs_idx) > 0:
#
#                # if there are other CETIs, compute contacts:
#                # encontrar todas (los IDs de) las CETIs dentro del radio D_max
#                query_point = [x, y]
#                idx = tree.query_ball_point(query_point, r=D_max)
#
#                # traverse all CETIs within reach
#                for k in idx:
#
#                    ID_old = CHATs_idx[k]
#
#                    Dx = np.sqrt(
#                                ((np.array(query_point) -
#                                 np.array(CETIs[ID_old][0][2:4]))**2).sum()
#                                )
#
#                    t_old_hola, t_old_chau = CETIs[ID_old][0][4:6]
#
#                    # check if contact is possible
#                    new_can_see_old = (
#                                    (Dx < D_max) &
#                                    (t_new_hola < t_old_chau + Dx) &
#                                    (t_new_chau > t_old_hola + Dx))
#
#                    if new_can_see_old:  # (·) new sees old
#
#                        # :start (type 3 event)
#                        t_new_see_old_start = max(t_old_hola + Dx, t_new_hola)
#                        next_event = [t_new_see_old_start, ID_new, ID_old, 3]
#                        t_forthcoming.add(next_event)
#
#                        # :end (type 4 event)
#                        t_new_see_old_end = min(t_old_chau + Dx, t_new_chau)
#                        next_event = [t_new_see_old_end, ID_new, ID_old, 4]
#                        t_forthcoming.add(next_event)
#
#                        contact = (ID_new, ID_old,
#                                   CETIs[ID_old][0][2], CETIs[ID_old][0][3],
#                                   t_new_see_old_start, t_new_see_old_end)
#                        CETIs[ID_new].append(contact)
#
#                    # check if contact is possible
#                    old_can_see_new = (
#                        (Dx < D_max) & (t_new_hola+Dx > t_old_hola) &
#                        (t_new_hola+Dx < t_old_chau))
#
#                    if old_can_see_new:  # (·) old sees new
#
#                        # :start (type 3 event)
#                        t_old_see_new_start = t_new_hola + Dx
#                        next_event = [t_old_see_new_start, ID_old, ID_new, 3]
#                        t_forthcoming.add(next_event)
#
#                        # :end (type 4 event)
#                        t_old_see_new_end = min(t_new_chau+Dx, t_old_chau)
#                        next_event = [t_old_see_new_end, ID_old, ID_new, 4]
#                        t_forthcoming.add(next_event)
#
#                        contact = (ID_old, ID_new,
#                                   CETIs[ID_new][0][2], CETIs[ID_new][0][3],
#                                   t_old_see_new_start, t_old_see_new_end)
#                        CETIs[ID_old].append(contact)
#
#            # agregar la CETI a la lista de posiciones de CETIs activas (CHATs)
#            CHATs.append([x, y])
#            CHATs_idx.append(ID_new)
#
#            # rehacer el árbol
#            tree = sp.cKDTree(data=CHATs)
#
#        if case == 2:
#            ID_bye = ID_emit
#
#            # eliminar la CETI a la lista de CETIs activas
#            # [ID, x, y, t_new_hola, t_new_chau]
#            try:
#                id_loc = CHATs_idx.index(ID_bye)
#                del CHATs[id_loc]
#                del CHATs_idx[id_loc]
#
#            except TypeError:
#                pass
#
#            # rehacer el árbol
#            if len(CHATs) > 0:
#                tree = sp.cKDTree(data=CHATs)
#
#        if case == 3: pass
#        if case == 4: pass
#
#        # eliminar el tiempo actual
#        t_forthcoming.remove_first()
#        # salir si no queda nada para hacer:
#        if t_forthcoming.size() < 1:
#            break
#        t_forthcoming.show()
#
#    return(CETIs)
##}}}
#

import ceti_exp

import numpy as np
import random
import pandas as pd
from scipy import spatial as sp
import time
import sys
import pickle
import ceti_tools as ct

# para recargar el modulo si se hacen modificaciones:
from importlib import reload
ct = reload(ct)

# para guardar la salida
from io import StringIO

# DEBUGING
# import pdb; pdb.set_trace()

# CASES:
# 1 - awakening
# 2 - doomsday
# 3 - contact
# 4 - blackout

random.seed(420)
np.random.seed(420)

message = "PRESS%ENTER%TO%CONTINUE"
message = "~~~~~~~~~~~~~~~~~~~~~~~"
message = "\n"+message*3+"\n"
message = ""
state = ['awakening','doomsday','contact ','blackout']

# FIXED AND SIMULATION VARIABLES
# {{{

# PARAMETERS :::

# radio interno de la zona galactica habitable, años luz
GHZ_inner = 1000.0

# radio interno de la zona galactica habitable, años luz
GHZ_outer = 5000.0

# tiempo medio, en años, que hay que esperar para que
# aparezca otra CETI en la galaxia
tau_awakening = 500.

# Tiempo medio, en años, durante el cual una CETI esta activa
tau_survive = 1000.

# Maxima distancia, en años luz, a la cual una CETI puede
# enviar o recibir mensajes
D_max = 3000.

# maximo tiempo para simular
t_max = 5000.

# flag to set interactive session
interactive = True

# }}}

# {{{

# VARIABLES :::

# lista de CETIs alguna vez activas: dictionary
CETIs = dict()

# lista de CETIs actualmente activas: ndarray
CHATs = []
CHATs_idx = []

# inicializacion del tiempo: scalar
t_now = 0

# inicializacion del ID: index
ID = 0

# lista de tiempos de eventos futuros: ordered list
# [time, ID_emit, ID_receive, case]
t_forthcoming = ct.OrderedList()

# estructura de arbol para buscar vecinos
try:
    del tree
except NameError:
    pass

# INITIALIZATION
# Simulation starts when the first CETI appears:
next_event = [0., 0, None, 1]
t_forthcoming.add(next_event)

if not interactive:
    old_stdout = sys.stdout
    result = StringIO()
    sys.stdout = result

# SIMULATION LOOP OVER TIME
while (t_now < t_max):

    t_now, ID_emit, ID_hear, case = t_forthcoming.head.getData()

    if interactive: wait = input(message)

    print('<<< t = %f >>> %10s of %3d' % (t_now, state[case-1], ID_emit))

    #print('[ case:%d | id:%d ]      <<< t = %f >>>' % (case, ID_emit, t_now))
    #print('   active CETIs', CHATs_idx)
    print(CHATs_idx)
    #t_forthcoming.show()

    # sys.stdout.write("\rTime: %f  (max=%f)\n" % (t_now, t_max))
    # sys.stdout.flush()
    # print('tiempo: %9.3f | ID_emit: %5d|CASE %1d'%(t_now,ID_emit,case))

    if case == 1:
        # -------------------------------------
        # NEW CETI
        # -------------------------------------
        # {{{

        # print('desaparece CETI con id:%d' % ID_emit)
        ID_new = ID_emit
        ID_next = ID_new + 1
        t_new_hola = t_now

        # sortear el lugar donde aparece dentro de la GHZ
        r = np.sqrt(random.random()*GHZ_outer**2 + GHZ_inner**2)
        o = random.random()*2.*np.pi
        x = r * np.cos(o)  # X position on the galactic plane
        y = r * np.sin(o)  # Y position on the galactic plane
        #print('posiciones para el ID %5d   :::: %8.3f %8.3f' % (ID_new, x, y))

        # sortear el tiempo de actividad
        t_new_active = np.random.exponential(tau_survive, 1)[0]
        t_new_chau = t_new_hola + t_new_active
        #print('tiempos para el ID %2d   :::: %8.3f %8.3f' % (ID_new, t_new_active, t_new_chau))

        # agregar el tiempo de desaparición a la lista de tiempos
        next_event = [t_new_chau, ID_new, None, 2]
        t_forthcoming.add(next_event)

        # agregar la CETI a la lista histórica
        CETIs[ID_new] = list()
        CETIs[ID_new].append(
            (ID_new, ID_new, x, y, t_new_hola, t_new_chau))

        # sortear el tiempo de aparición de la próxima CETI
        t_next_awakening = np.random.exponential(tau_awakening, 1)[0]
        t_next_awakening = t_new_hola + t_next_awakening
        if t_next_awakening < t_max:
            next_event = [t_next_awakening, ID_next, None, 1]
            t_forthcoming.add(next_event)

        if len(CHATs_idx) > 0:

            #print('entering loop to compute contacts')

            # if there are other CETIs, compute contacts:
            # encontrar todas (los IDs de) las CETIs dentro del radio D_max
            query_point = [x, y]
            idx = tree.query_ball_point(query_point, r=D_max)

            # traverse all CETIs within reach
            for k in idx:

                #print(k)

                ID_old = CHATs_idx[k]

                Dx = np.sqrt(
                            ((np.array(query_point) -
                             np.array(CETIs[ID_old][0][2:4]))**2).sum()
                            )

                t_old_hola, t_old_chau = CETIs[ID_old][0][4:6]

                # check if contact is possible
                new_can_see_old = (
                                (Dx < D_max) &
                                (t_new_hola < t_old_chau + Dx) &
                                (t_new_chau > t_old_hola + Dx))

                if new_can_see_old:  # (·) new sees old

                    # :start (type 3 event)
                    t_new_see_old_start = max(t_old_hola + Dx, t_new_hola)
                    next_event = [t_new_see_old_start, ID_new, ID_old, 3]
                    t_forthcoming.add(next_event)

                    # :end (type 4 event)
                    t_new_see_old_end = min(t_old_chau + Dx, t_new_chau)
                    next_event = [t_new_see_old_end, ID_new, ID_old, 4]
                    t_forthcoming.add(next_event)

                    contact = (ID_new, ID_old,
                               CETIs[ID_old][0][2], CETIs[ID_old][0][3],
                               t_new_see_old_start, t_new_see_old_end)
                    CETIs[ID_new].append(contact)

                # check if contact is possible
                old_can_see_new = (
                    (Dx < D_max) & (t_new_hola+Dx > t_old_hola) &
                    (t_new_hola+Dx < t_old_chau))

                if old_can_see_new:  # (·) old sees new

                    # :start (type 3 event)
                    t_old_see_new_start = t_new_hola + Dx
                    next_event = [t_old_see_new_start, ID_old, ID_new, 3]
                    t_forthcoming.add(next_event)

                    # :end (type 4 event)
                    t_old_see_new_end = min(t_new_chau+Dx, t_old_chau)
                    next_event = [t_old_see_new_end, ID_old, ID_new, 4]
                    t_forthcoming.add(next_event)

                    contact = (ID_old, ID_new,
                               CETIs[ID_new][0][2], CETIs[ID_new][0][3],
                               t_old_see_new_start, t_old_see_new_end)
                    CETIs[ID_old].append(contact)

        # agregar la CETI a la lista de posiciones de CETIs activas (CHATs)
        CHATs.append([x, y])
        CHATs_idx.append(ID_new)

        # rehacer el árbol
        tree = sp.cKDTree(data=CHATs)

        # }}}

    if case == 2:
        # -------------------------------------
        # END CETI
        # -------------------------------------
        # {{{
        # print('desaparece CETI con id:%d' % ID_emit)

        ID_bye = ID_emit

        # eliminar la CETI a la lista de CETIs activas
        # [ID, x, y, t_new_hola, t_new_chau]
        try:
            id_loc = CHATs_idx.index(ID_bye)
            del CHATs[id_loc]
            del CHATs_idx[id_loc]

        except TypeError:
            pass

        # rehacer el árbol
        if len(CHATs) > 0:
            tree = sp.cKDTree(data=CHATs)

        # }}}

    if case == 3:
        # -------------------------------------
        # NEW CONTACT
        # -------------------------------------
        # {{{
        pass
        # print('primer contacto entre %d y %d' % (ID_emit, ID_hear))

        # }}}

    if case == 4:
        # -------------------------------------
        # END CONTACT
        # -------------------------------------
        # {{{
        pass
        # print('ultimo contacto entre %d y %d' % (ID_emit, ID_hear))
        # }}}

    # eliminar el tiempo actual
    t_forthcoming.remove_first()

    # salir si no queda nada para hacer:
    if t_forthcoming.size() < 1:
        break

    t_forthcoming.show()
    print(CHATs_idx)

print('\n\n')
ct.ShowCETIs(CETIs)

if not interactive:
    sys.stdout = old_stdout
    result_string = result.getvalue()

    i = input('experiment code: ').zfill(3)

    fname = 'experiment_' + i + '.dat'
    savef = open(fname,'w')
    savef.write(result_string)
    savef.close()

# El experimento consiste es correr interactivamente hasta que se
# corrobora que anda bien.  Luego puse la funcion en un modulo, la
# llamo con la misma semilla y anda igual
#
#fCETIs = ceti_exp.ceti_exp(
#    GHZ_inner, GHZ_outer, tau_awakening, tau_survive, D_max, t_max
#)
#ct.ShowCETIs(fCETIs)

max_stars_repo_name: mlares/hearsay

max_issues_repo_name: mlares/hearsay

max_forks_repo_name: mlares/hearsay

reason: import numpy,from scipy


text:
# -*- coding: utf-8 -*-
# test_sinc_decimate.py
# This module provides the tests for the sinc_decimate function.
# Copyright 2014 Giuseppe Venturini
# This file is part of python-deltasigma.
#
# python-deltasigma is a 1:1 Python replacement of Richard Schreier's
# MATLAB delta sigma toolbox (aka "delsigma"), upon which it is heavily based.
# The delta sigma toolbox is (c) 2009, Richard Schreier.
#
# python-deltasigma is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# LICENSE file for the licensing terms.

"""This module provides the test class for the sinc_decimate() function.
"""

import unittest
import numpy as np
import deltasigma as ds

from nose.tools import raises

class TestSincDecimate(unittest.TestCase):
    """Test class for sinc_decimate()"""

    def test_sinc_decimate_1(self):
        """Test function for sinc_decimate() 1/4"""
        x = [1]*10
        self.assertTrue(np.allclose(ds.sinc_decimate(x, 1, 10), [1.]))

    def test_sinc_decimate_2(self):
        """Test function for sinc_decimate() 2/4"""
        x = [1]*10
        self.assertTrue(np.allclose(ds.sinc_decimate(x, 2, 5), [.6, 1.]))

    def test_sinc_decimate_3(self):
        """Test function for sinc_decimate() 3/4"""
        x = [1]*10
        x = np.cumsum(np.asarray(x))
        self.assertTrue(np.allclose(ds.sinc_decimate(x, 1, 5), [3., 8.]))
        self.assertTrue(np.allclose(ds.sinc_decimate(x, 2, 5), [1.4, 6.]))
        self.assertTrue(np.allclose(ds.sinc_decimate(x, 3, 3), [0.55555556, 3., 6.]))

    @raises(ValueError)
    def test_sinc_decimate_4(self):
        """Test function for sinc_decimate() 4/4"""
        x = [1]*10
        ds.sinc_decimate(x, 2, 15)

max_issues_repo_name: thomasrussellmurphy/python-deltasigma

max_forks_repo_name: thomasrussellmurphy/python-deltasigma


text:
# -*- coding: utf-8 -*-
"""
Created on Fri Jan  7 12:22:59 2022

@author: kkrao
"""

import os

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
import sklearn.metrics

import init

years = range(2016, 2022)

# df = pd.DataFrame()
# for year in years:
#     df = df.append(pd.read_csv(os.path.join(init.dir_root, "data",\
#             "lightning", "nldn",f"nldn-tiles-{year:4d}.csv"), header=2,\
#             names = ["date","longitude","latitude","count"],
#             dtype = {"date":str,"longitude":float,"latitude":float,"count":int}))
# # df.date = pd.to_datetime(df.date, format = "%Y%m%d")
# df.date = df.date.str[:4]+"-"+df.date.str[4:6]+"-"+df.date.str[6:]
# df.date = df.date.astype(str)
# df.index.name = "fid"
# df.to_csv(os.path.join(init.dir_root, "data",\
#             "lightning", "nldn",f"nldn-tiles-15-feb-2022.csv"))
# # df.head()
# # df.head().to_csv(r"C:\Users\kkrao\Desktop\trial.csv")

# """
# While uploading as asset to GEE leave the Date fiel blank, but in lat and lon
# type in "longitude" and "latitude"
# """

df = pd.read_csv(os.path.join(init.dir_root, "data","gee",\
                               "lightnings_22_feb_2022_2016_2021_California.csv"))
# df = pd.read_csv(os.path.join(init.dir_root, "data","gee",\
                              # "lightnings_11_jan_2022_California.csv"))

    #     ## has lagged vpd 4 months, ppt 1 year, erc 15 days
df = df.rename(columns = {"vpd_4":"vpd_4m","ppt_1":"ppt_1y","erc_15":"erc_15d"})
# df = df.rename(columns = {"vpd":"vpd_4m","ppt":"ppt_1y","erc_15":"erc_15d"})
# rdf = pd.read_csv(os.path.join(init.dir_root, "data","gee",\
#                   "lightnings_18_jan_2022_climatology_ndvi_California.csv"))
# merge_cols = ["fid","system:index","date",".geo"]
# cols = ["mean_t","mean_vpd","mean_ppt","mean_ndvi"]
# df = df.merge(rdf[merge_cols + cols], \
#           on =merge_cols)

# # df = df.join(rdf[])
# rdf = pd.read_csv(os.path.join(init.dir_root, "data","gee",\
#                   "lightnings_18_jan_2022_fire_size_California.csv"))

# cols = ["fire_size"]
# df = df.merge(rdf[merge_cols + cols], \
#           on =merge_cols)

# rdf = pd.read_csv(os.path.join(init.dir_root, "data","gee",\
#               "lightnings_25_jan_2022_prism_instantenous_California.csv"))

# cols = ["vpd","ppt","erc"]
# df = df.merge(rdf[merge_cols + cols], \
#           on =merge_cols)

# rdf = pd.read_csv(os.path.join(init.dir_root, "data","gee",\
#               "lightnings_25_jan_2022_biomass_canopy_height_California.csv"))
# cols = ["agb","canopy_height"]
# df = df.merge(rdf[merge_cols + cols], \
#           on =merge_cols)

# rdf = pd.read_csv(os.path.join(init.dir_root, "data","gee",\
#               "lightnings_31_jan_2022_elevation_California.csv"))
# cols = ["elevation"]
# df = df.merge(rdf[merge_cols + cols], \
#           on =merge_cols)

# rdf = pd.read_csv(os.path.join(init.dir_root, "data","gee",\
#               "lightnings_31_jan_2022_ndvi_1m_California.csv"))
# cols = ["ndvi_1m"]
# df = df.merge(rdf[merge_cols + cols], \
#           on =merge_cols)

# rdf = pd.read_csv(os.path.join(init.dir_root, "data","gee",\
#               "lightnings_31_jan_2022_wind_1d_California.csv"))
# cols = ["wind"]
# df = df.merge(rdf[merge_cols + cols], \
#           on =merge_cols)

df.lc = df.lc.round()
df.fire = 1-df.fire.isnull()

df["longitude"] = df[".geo"].str.split("[").str[1].str.split(",").str[0]
df["latitude"] = df[".geo"].str.split("[").str[1].str.split(",").str[1].str[:-2]

df.longitude = df.longitude.astype(float).round(1)
df.latitude = df.latitude.astype(float).round(1)

df["system:time_start"] = df["system:time_start"]/1000/60/60/24
df["system:time_start"] = df["system:time_start"].astype(int)

df = df.dropna()
df = df.loc[df.lc.isin(init.thresh["extreme"].keys())].copy()
# df.shape

df["z_extreme"] = (df.lfmc <= list( map(init.thresh["extreme"].get, df.lc) )).astype(int)
df["z_high"] = (df.lfmc <= list( map(init.thresh["high"].get, df.lc) )).astype(int)
df["z_moderate"]= (df.lfmc <= list( map(init.thresh["moderate"].get, df.lc) )).astype(int)
df["z_ndvi"] = (df.ndvi_1m >= df.mean_ndvi).astype(int)
df["z_wind"] = (df.wind >= df.wind.quantile(0.50)).astype(int)
df["z_vpd"] = (df.vpd >= df.mean_vpd).astype(int)
# df.z_moderate.mean()
df.to_csv(os.path.join(init.dir_root, "data","r","rct_22_apr_2022.csv"))
# df.head()['.geo'].apply(get_lon)
fig, ax = plt.subplots()
sns.kdeplot(x = "lfmc", data =  df.loc[(df.fire==0)],alpha = 0.5,\
            ax = ax, label = "no fire")
sns.kdeplot(x = "lfmc", data =  df.loc[(df.fire==1)],alpha = 0.5,\
            ax = ax, label = "fire")
ax.legend()

"""

41.0        5
42.0     9334
43.0       79
52.0    24823
71.0     2886
81.0      177
Name: lc, dtype: int64"""

cm = sklearn.metrics.confusion_matrix(df.fire, df.z_ndvi)
cm[1,1]*cm[0,0]/cm[1,0]/cm[0,1] #= 4.59
print("Extreme\n",sklearn.metrics.confusion_matrix(df.fire, df.z_extreme))
print("High\n",sklearn.metrics.confusion_matrix(df.fire, df.z_high))
print("Moderate\n",sklearn.metrics.confusion_matrix(df.fire, df.z_moderate))


text:
import io
import logging
import os

import cv2 as cv
import en_core_web_lg
import fitz
import numpy as np
import streamlit as st
from PIL import Image
from presidio_analyzer import AnalyzerEngine
from presidio_image_redactor import ImageAnalyzerEngine

SCORE_THRESHOLD = 0.5

def local_css(file_name):
    with open(file_name) as f:
        st.markdown(''.format(f.read()), unsafe_allow_html=True)

# Load models - Because this is cached it will only happen once.
@st.cache(allow_output_mutation=True)
def load_nlp():
    logging.info("loading spacy nlp")
    return en_core_web_lg.load()

@st.cache(allow_output_mutation=True)
def load_analyzer():
    logging.info("loading analyzer")
    return AnalyzerEngine()

@st.cache(allow_output_mutation=True)
def load_image_analyzer():
    logging.info("loading image analyzer")
    return ImageAnalyzerEngine()

@st.cache(allow_output_mutation=True)
def load_face_detector():
    logging.info("loading face cascade detector")
    face_cascade = cv.CascadeClassifier()
    face_cascade.load("require/haarcascade_frontalface_default.xml")
    return face_cascade

def process_face(image):
    gray = cv.cvtColor(image, cv.COLOR_BGR2GRAY)
    gray = cv.equalizeHist(gray)
    face_detector = load_face_detector()
    faces = face_detector.detectMultiScale(gray)
    return faces

def process_image_ocr(image):
    ocr_analyzer = load_image_analyzer()
    res = ocr_analyzer.analyze(image)
    return np.array([[spot.left, spot.top, spot.width, spot.height] for spot in res])

def process_image(image):
    faces = process_face(image)
    spots = process_image_ocr(image)
    if (len(faces) > 0) & (len(spots) > 0):
        return np.concatenate((faces, spots), axis=0)
    elif len(faces) > 0:
        return faces
    elif len(spots) > 0:
        return spots
    else:
        return np.array([])

def process_text(text):
    analyzer = load_analyzer()
    # Call analyzer to get results
    results = analyzer.analyze(text=text, language='en', )
    return filter(lambda pii: pii.score > SCORE_THRESHOLD, results)

def process_page(page):
    ## scan for images
    data = page.get_text('dict')
    img_data = list(filter(lambda block: block['type'] == 1, data['blocks']))

    for i, img in enumerate(img_data):
        img_bytes = img['image']
        if img['bpc'] == 8:
            decoded = cv.imdecode(np.frombuffer(img_bytes, np.uint8), -1)
        elif img['bpc'] == 16:
            decoded = cv.imdecode(np.frombuffer(img_bytes, np.uint16), -1)
        else:
            logging.error(f"Unknown bitrate - img.bpc:{img['bpc']}")
            raise (f"Unknown bitrate - img.bpc:{img['bpc']}")
        faces_spots = process_image(decoded)
        img_data[i]['faces_spots'] = faces_spots

    ## scan text
    text = page.get_text('text')
    # text_data = list(filter(lambda block: block['type'] == 0, data['blocks']))
    # text_lines = [span for block in text_data for line in block['lines'] for span in line['spans']]
    piis = process_text(text)
    return img_data, piis

def process_pdf(doc):
    anot_texts = []
    page_imgs = []
    for j in range(doc.page_count):
        page = doc.load_page(j)
        pdf_img_data, pdf_page_piis = process_page(page)

        ## convert page to numpy image
        pix = page.get_pixmap()
        with Image.frombytes("RGB", [pix.width, pix.height], pix.samples) as pil_:
            page_img = np.array(pil_)
        ## draw faces on page
        for block in pdf_img_data:
            x0, y0, x1, y1 = block['bbox']
            for face in block['faces_spots']:
                x, y, w, h = face
                rel_x0 = (x / block['width']) * (x1 - x0)
                rel_x1 = ((x + w) / block['width']) * (x1 - x0)
                rel_y0 = (y / block['height']) * (y1 - y0)
                rel_y1 = ((y + h) / block['height']) * (y1 - y0)
                cv.rectangle(page_img, (int(x0 + rel_x0), int(y0 + rel_y0)),
                             (int(x0 + rel_x1), int(y0 + rel_y1)), (255, 0, 0), 2)
        # construct colored text
        text = page.get_text('text')
        i = 0
        pretty_text = "
"
        for pii in pdf_page_piis:
            pretty_text += text[i:pii.start]
            pretty_text += " " + \
                           text[pii.start:pii.end] + \
                           f"{str(pii.entity_type)} "
            i = pii.end
        pretty_text += "
"

        page_imgs.append(page_img)
        anot_texts.append(pretty_text)
    return anot_texts, page_imgs

def process_image_final(image):
    faces_spots = process_image(image)
    output_img = image.copy()
    for spot in faces_spots:
        x, y, w, h = spot
        cv.rectangle(output_img, (int(x), int(y)), (int(x + w), int(y + h)), (255, 0, 0), 2)

    anot_texts = []  # text from ocr is not used now, only detected PII
    page_imgs = [output_img]
    return anot_texts, page_imgs

def process_file(buff, ext):
    anot_texts, page_imgs = [], []
    if ext == 'pdf':
        with fitz.open(stream=buff.read(), filetype="pdf") as doc:
            anot_texts, page_imgs = process_pdf(doc)
    elif ext in ["jpg", "jpeg", "png", "PNG", 'JPG']:
        with Image.open(buff) as buf:
            image = np.array(buf)
        if image.shape[-1] == 4:
            image = cv.cvtColor(image, cv.COLOR_RGBA2RGB)
        anot_texts, page_imgs = process_image_final(image)
    with st.beta_expander("Analyzed text"):
        # svg = page.get_svg_image()
        # b64 = base64.b64encode(svg.encode("utf-8")).decode("utf-8")
        # html = r''.format(b64)
        # st.write(html, unsafe_allow_html=True)
        for pretty_text in anot_texts:
            st.markdown(pretty_text, unsafe_allow_html=True)

    for page_img in page_imgs:
        st.image(page_img)

def main():
    st.title('Detect Personal Information')
    img_file_buffer = st.file_uploader("Upload an pdf/image", type=["png", "jpg", "jpeg", "pdf", "PNG", 'JPG'])

    if img_file_buffer is not None:
        ext = os.path.splitext(img_file_buffer.name)[-1].replace(".", "")
        file_buff = img_file_buffer
    else:
        logging.info("No input file - using sample image.")
        ext = 'png'
        with open("data/tds_jessica.png", 'rb') as buf:
            file_buff = io.BytesIO(buf.read())
    process_file(file_buff, ext)

if __name__ == "__main__":
    logging.basicConfig(level=logging.INFO)
    # logging.info(fitz.__doc__)
    local_css("style.css")
    main()

max_stars_repo_name: konvica/read-presonal-data

max_issues_repo_name: konvica/read-presonal-data

max_forks_repo_name: konvica/read-presonal-data


text:
#!/usr/bin/env python3

"""
Benchmark representation of unique side-chain contacts in the PDB for a fragment

"""

import os
import pickle
from pprint import pprint
import random
from multiprocessing import Pool

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
from matplotlib.offsetbox import OffsetImage, AnnotationBbox
import pandas as pd
import seaborn as sns

from rdkit import Chem
from rdkit.Chem import Draw
from rdkit.Chem.Draw.MolDrawing import MolDrawing, DrawingOptions

from ..utils import *

def fragment_representation_bm(args):
    """
    Bootstrap AUC plot for unique side-chain fragment contact recovery for specific fragments

    Usage:
      bsff benchmark fragment_representation <user_defined_dir> <clusters_dir> [options]
      bsff benchmark fragment_representation figures <user_defined_dir> <clusters_dir> <bm_csv_dir> [options]
      bsff benchmark fragment_representation draw_fragments <fragment_inputs_csv>

    Arguments:
      <user_defined_dir>      Path to project root directory
      <clusters_dir>          Path to directory containing clusters generated by clusters command
      <bm_csv_dir>            Path to directory containing fragment recovery benchmark csvs
      <fragment_inputs_csv>   Path to directory containing defined fragment PDBs

    Options:
      --csv=<csv_path>        Previously generated benchmark data
      -e, --error_bands       Use error bands instead of bars for standard deviations
      -m, --maintext          Generate Figures for main text
      -p, --processes         Number of processes to run in parallel
    """

    """
    Before this function gets called:
      * Defined fragment (Use full amino acids to standardize atom names, need to figure out params for fragments)
      * Assemble fuzzball (SAVE REDUNDANT CONTACTS!!!)

    What happens here:
      * Assemble Cluster/PDB mapping dataframe from ProDy clusters
      * Remove singletons (keep track of this)
      * Find source of all structure codes in PDB
      * Bootstrap by taking 5% increments of PDB and intersect with dataframe
      * Save proportion of clusters recovered, plot

    Bam, done.

    All steps up to clustering to not rely on params files, so I might be able to get away with just generating fragments
    from ideal ligands from the PDB/LigandExpo.
    """

    user_defined_dir = args['<user_defined_dir>']
    clusters_dir = args['<clusters_dir>']
    processes = args['--processes'] if args['--processes'] else 12

    # Why.
    if args['draw_fragments']:
        draw_fragments(args['<fragment_inputs_csv>'])

    if args['figures']:
        png_dir = generate_fragment_pngs(user_defined_dir)
        # generate_cluster_counts_hist(clusters_dir, format_maintext=args['--maintext'])
        generate_auc_figures(args['<bm_csv_dir>'], png_dir, error_bands=args['--error_bands'], format_maintext=args['--maintext'])

    else:
        pool = Pool(processes=processes)
        cluster_dir_list = [cluster_dir for cluster_dir in directory_check(clusters_dir, base_only=False)]
        pool.map(perform_fragment_recovery_bm, cluster_dir_list)
        pool.close()
        pool.join()

def perform_fragment_recovery_bm(cluster_dir):
    """
    Perform fragment recovery benchmark
    :return:
    """

    # Import all PDBIDs from file
    pdb_entries_path = os.path.join(os.path.dirname(__file__), '..', '..', 'Additional_Files', 'entries.idx')
    pdb_df = pd.read_csv(pdb_entries_path,
                         delimiter='\t',
                         names=['IDCODE', 'HEADER', 'ACCESSION DATE', 'COMPOUND', 'SOURCE', 'AUTHOR LIST', 'RESOLUTION',
                                'EXPERIMENT TYPE'],
                         skiprows=2
                         )

    pdb_set = set([idcode.lower() for idcode in pdb_df['IDCODE']])
    all_pdbids_count = len(pdb_set)
    pdb_coverage_fractions = list(np.arange(0, 1.05, 0.05))

    recovery_list = list()

    # For loops all the way down...
    fragment_name = os.path.basename(cluster_dir)

    with open(os.path.join(cluster_dir, f'{fragment_name}-PDB_Sources.pickle'), 'rb') as pdb_source_pickle:
        cluster_source_dict = pickle.load(pdb_source_pickle)

    # todo: flatten cluster_source_dict!

    flattened_cluster_source_dict = dict()
    for cat_cluster_idx, contact_cluster_dict in cluster_source_dict.items():
        for contact_cluster_idx, contact_clusters in contact_cluster_dict.items():
            flattened_cluster_source_dict[(cat_cluster_idx, contact_cluster_idx)] = contact_clusters

    cluster_source_dict = flattened_cluster_source_dict
    pprint(cluster_source_dict)

    # Remove singletons
    singleton_cluster_ids = list()
    total_clusters = 0
    total_cluster_members = 0

    for cluster_id in cluster_source_dict:
        if len(cluster_source_dict[cluster_id]['list']) <= 1:
            singleton_cluster_ids.append(cluster_id)
        else:
            total_clusters += 1
            total_cluster_members += len(cluster_source_dict[cluster_id]['list'])

    for cluster_id in singleton_cluster_ids:
        del cluster_source_dict[cluster_id]

    print(total_clusters, total_cluster_members)

    for fraction in pdb_coverage_fractions:
        samples = int(all_pdbids_count * fraction)

        for i in range(1000):
            cluster_hits = 0
            selected_pdbids = set(random.sample(pdb_set, samples))

            for cluster_id, cluster_members in cluster_source_dict.items():
                cluster_pop = len(set([a.lower() for a in cluster_members['list']]) & selected_pdbids)
                if cluster_pop > 0:
                    cluster_hits += 1

            trial_dict = {'fragment': fragment_name,
                          'fraction': fraction,
                          'total_residues': total_cluster_members,
                          'total_clusters': total_clusters,
                          'cluster_hits': cluster_hits,
                          'recovered': cluster_hits / total_clusters
                          }
            recovery_list.append(trial_dict)

            print(fragment_name, fraction, total_clusters, cluster_hits)

    df = pd.DataFrame(recovery_list)
    df.to_csv(f'{os.path.basename(fragment_name)}-cluster_recovery.csv')

    return df

def generate_fragment_pngs(user_defined_dir):
    """
    Create pngs of fragments for embedding into plots
    """
    fragment_inputs = os.path.join(user_defined_dir, 'Inputs', 'Fragment_Inputs')
    fragment_inputs_csv = os.path.join(fragment_inputs, 'Fragment_inputs.csv')
    df = pd.read_csv(fragment_inputs_csv)

    frag_smarts_dict = {f"Fragment_{row['Fragment']}": re.sub(r'[^\[\]]+(?=\])', lambda x: f"{x.group().split(';')[0]}",
                                              row['SMILES_fragment']).upper() for index, row in df.iterrows()}

    png_dir = os.path.join(user_defined_dir, 'Fragment_PNGs')
    os.makedirs(png_dir, exist_ok=True)

    for fragment_name, fragment_smarts in frag_smarts_dict.items():
        fragment_mol_H = Chem.MolFromSmarts(fragment_smarts)
        if fragment_mol_H is None:
            fragment_mol_H = Chem.MolFromPDBFile(os.path.join(fragment_inputs, f'{fragment_name}.pdb'), removeHs=True)

        fragment_mol = Chem.RemoveHs(fragment_mol_H)
        fragment_mol.UpdatePropertyCache(strict=False)  # https://github.com/rdkit/rdkit/issues/146
        Chem.rdDepictor.Compute2DCoords(fragment_mol)
        img = Draw.MolToMPL(fragment_mol, fitImage=True)
        plt.axis('off')
        img.savefig(os.path.join(png_dir, f'{fragment_name}.png'), bbox_inches='tight', transparent=True)

    return png_dir

def generate_cluster_counts_hist(clusters_dir, format_maintext=False):
    """Histogram to visualize cluster count/occupancy"""
    for cluster_dir in directory_check(clusters_dir, base_only=False):
        fragment_name = os.path.basename(cluster_dir)
        df = pd.read_csv(os.path.join(clusters_dir, fragment_name, f'{fragment_name}_report.csv'))
        df.sort_values(by=['cluster_members'], ascending=False, inplace=True)

        fig = plt.figure()
        ax = sns.barplot(x=list(range(1, len(df['cluster_members']) + 1)), y=df['cluster_members'], palette=['#007CBE']*len(df['cluster_members']))

        ax.set(ylim=(0, df['cluster_members'].max()))
        ax.set(xlim=(0, len(df['cluster_members'])))

        # 5% Vertical line
        plt.axvline(x=len(df['cluster_members'])/20, color='#EB093C')

        x_clusters = len(df['cluster_members'])
        x_cluster_ticks = len(df['cluster_members']) // 10

        plt.xticks(range(0, x_clusters, x_cluster_ticks), range(0, x_clusters, x_cluster_ticks))
        ax.set_xlabel('Clusters')
        ax.set_ylabel('Cluster Members')

        # Format axes for main text/suppmat
        if format_maintext:
            ax.tick_params(labelsize=16, labelrotation=45)
            plt.rcParams.update({'font.size': 16})

        fig.add_axes(ax)
        figure_name = f'{fragment_name}-cluster_occupancy.png'
        fig.savefig(figure_name, dpi=300, bbox_inches='tight', transparent=True)
        plt.close()

def generate_auc_figures(bm_csv_dir, png_dir, error_bands=False, format_maintext=False):
    """
    Generate figures for fragment recovery benchmark
    :param df:
    :return:
    """
    global fraction
    fraction_recovery_list = []

    annotation_x_positions = (0.3, 0.3, 0.3, 0.3, 0.3) if format_maintext else (0.3, 0.3, 0.3, 0.2, 0.2)
    annotation_y_positions = (0.33, 0.26, 0.19, 0.12, 0.05) if format_maintext else (0.3, 0.25, 0.2, 0.15, 0.1, 0.05)
    annotation_fontsize = 16 if format_maintext else 10

    for csv in os.listdir(bm_csv_dir):
        if csv.endswith('.csv'):
            df = pd.read_csv(os.path.join(bm_csv_dir, csv))

            for fragment, fragment_df in df.groupby('fragment'):
                print(fragment)

                fig = plt.figure()
                plt.rcParams.update({'font.size': annotation_fontsize})

                if error_bands:
                    muh_plot = sns.lineplot(x='fraction', y='recovered', data=fragment_df, ci='sd', err_style='band', estimator=np.mean)
                    sns.regplot(x='fraction', y='recovered', data=fragment_df, x_estimator=np.mean, x_ci=None, fit_reg=False, ax=muh_plot, scatter_kws={'color': list(sns.color_palette("Blues_r")[1])})
                else:
                    muh_plot = sns.regplot(x='fraction', y='recovered', data=fragment_df, x_ci='sd', x_estimator=np.mean, fit_reg=False)

                # Format axes for main text/suppmat
                if format_maintext:
                    muh_plot.tick_params(labelsize=16, labelrotation=45)

                muh_plot.set(ylim=(0, 1.1))
                muh_plot.set(xlim=(0, 1.1))

                muh_plot.set_xlabel('PDB Fraction')
                muh_plot.set_ylabel('Clusters Recovered')

                # Fraction of PDB for 80% recovery
                for fraction, fraction_df in fragment_df.groupby('fraction'):
                    avg_recovered = fraction_df['recovered'].mean()
                    if avg_recovered > 0.8: break
                fraction_recovery_list.append({'fraction': fraction, 'fragment': fragment})

                # AUC calculation
                # auc_info = [(fraction, recovered['recovered'].mean()) for fraction, recovered in fragment_df.groupby('fraction')]
                # x_coords = [a[0] for a in auc_info]
                # y_coords = [a[1] for a in auc_info]
                # muh_plot.text(0.7, 0.05, f'AUC: {np.trapz(y_coords, x_coords):.4f}')

                # embed fragment RDKit image
                if not format_maintext:
                    fragment_png = mpimg.imread(os.path.join(png_dir, f'{fragment}.png'))
                    fragment_inlay = OffsetImage(fragment_png, zoom=0.15)
                    ab = AnnotationBbox(fragment_inlay, (0.915, 0.25), pad=0)
                    muh_plot.add_artist(ab)

                # Add annotations
                # if not format_maintext:
                #     muh_plot.text(annotation_x_positions[0], annotation_y_positions[0], f'{fragment}')

                muh_plot.text(annotation_x_positions[1], annotation_y_positions[1], f'Total Residues: {fragment_df["total_residues"].median():.0f}')
                muh_plot.text(annotation_x_positions[2], annotation_y_positions[2], f'Total Clusters: {fragment_df["total_clusters"].median():.0f}')
                muh_plot.text(annotation_x_positions[3], annotation_y_positions[3], f'Fraction PDB sampled for')
                muh_plot.text(annotation_x_positions[4], annotation_y_positions[4], f'>80% contact recovery: {fraction:.2f}')

                if format_maintext:
                    muh_plot.tick_params(labelsize=16)

                fig.add_axes(muh_plot)
                figure_name = f'{fragment}-cluster_recovery-{"band" if error_bands else "bars"}.png'
                fig.savefig(figure_name, dpi=600, bbox_inches='tight', transparent=True)
                plt.close()

    # Fraction recovered Histogram
    df = pd.DataFrame(fraction_recovery_list)

    sns.set(style="whitegrid", palette="hls")
    sns.set_context("paper", font_scale=1.5)

    fig, ax = plt.subplots(figsize=(6,6))
    sns.countplot(df['fraction'].round(2), color='#007CBE', ax=ax)
    fig.suptitle('Fraction of PDB Required for >80%\nContact Recovery Across All Tested Fragments')
    ax.set_xlabel('Fraction', fontsize=12)
    ax.set_ylabel('Count', fontsize=12)

    fig.savefig('RecoveryHist.tiff', dpi=300)
    df.to_csv('FractionRecovered.csv')

# todo: how do I embed these into the plots?
def draw_fragments(fragment_inputs_csv):
    """
    Use RDKit to draw fragments and output to .png
    :return:
    """

    df = pd.read_csv(fragment_inputs_csv)
    frag_smarts_dict = {f"Fragment_{row['Fragment']}": re.sub(r'[^\[\]]+(?=\])', lambda x: f"{x.group().split(';')[0]}", row['SMILES_fragment']).upper() for index, row in df.iterrows()}

    name_list = list()
    mol_list = list()

    for fragment_name, fragment_smarts in frag_smarts_dict.items():
        print(fragment_smarts)
        fragment_mol = Chem.RemoveHs(Chem.MolFromSmarts(fragment_smarts))
        fragment_mol.UpdatePropertyCache(strict=False)  # https://github.com/rdkit/rdkit/issues/146
        Chem.rdDepictor.Compute2DCoords(fragment_mol)
        Draw.MolToFile(fragment_mol, f'{fragment_name}.png')

        name_list.append(fragment_name)
        mol_list.append(fragment_mol)

    img = Draw.MolsToGridImage(mol_list, molsPerRow=4, subImgSize=(200, 200), legends=name_list)
    img.save('AUC_Fragments.png')

max_stars_repo_name: jaaamessszzz/BindingSitesFromFragments

max_issues_repo_name: jaaamessszzz/BindingSitesFromFragments

max_forks_repo_name: jaaamessszzz/BindingSitesFromFragments


text:
## -------------------------------------------------------------------------------------------------
## -- Project : FH-SWF Automation Technology - Common Code Base (CCB)
## -- Package : mlpro
## -- Module  : test_pool_policies
## -------------------------------------------------------------------------------------------------
## -- History :
## -- yyyy-mm-dd  Ver.      Auth.    Description
## -- 2021-12-08  1.0.0     MRD      Creation
## -------------------------------------------------------------------------------------------------

"""
Ver. 1.0.0 (2021-12-08)

Unit test classes for policies.
"""

import pytest
import random
import gym
import numpy as np
from mlpro.rl.models import *
from mlpro.wrappers.openai_gym import WrEnvGYM2MLPro
from mlpro.rl.pool.policies.dummy import MyDummyPolicy

## Instructions
# Please Include your own policy implementation class on the test list parameter

## -------------------------------------------------------------------------------------------------
@pytest.mark.parametrize("policy_cls", [MyDummyPolicy])
def test_pool_policies(policy_cls):
    class MyScenario (RLScenario):

        C_NAME      = 'Matrix'

        def _setup(self, p_mode, p_ada, p_logging):
            # 1 Setup environment
            gym_env     = gym.make('CartPole-v1')
            self._env   = WrEnvGYM2MLPro(gym_env, p_logging=False)

            # 2 Setup and return standard single-agent with own policy
            return Agent(
                p_policy=policy_cls(
                    p_observation_space=self._env.get_state_space(),
                    p_action_space=self._env.get_action_space(),
                    p_buffer_size=10,
                    p_ada=p_ada,
                    p_logging=p_logging
                ),
                p_envmodel=None,
                p_name='Smith',
                p_ada=p_ada,
                p_logging=p_logging
            )

    # 2.4 Create and run training object
    training = RLTraining(
            p_scenario_cls=MyScenario,
            p_cycle_limit=100,
            p_max_adaptations=0,
            p_max_stagnations=0,
            p_logging=False
    )

    training.run()

max_stars_repo_name: fhswf/MLPro

max_issues_repo_name: fhswf/MLPro

max_forks_repo_name: fhswf/MLPro


text:
import streamlit as st
import numpy as np
import pandas as pd
from sklearn import datasets
from sklearn.model_selection import train_test_split
import matplotlib.pyplot as plt
import seaborn as sns
plt.style.use('seaborn')  # change the default style

def app():
    st.set_option('deprecation.showPyplotGlobalUse', False)

    st.title('EDA')

    st.subheader("Portuguese Bank Marketing Dataset")
    df = pd.read_csv('projectdataset-1.csv')
    df.Class.replace((1, 2), ('no', 'yes'), inplace=True)
# the table
    st.write(df)

#Pie Chart so show % of sub/not sub
    st.subheader("Subscribed vs Not Subscribed")
    labels = ["Not \nsubscribed", "Subscribed"]
    explode = (0, 0.1)
    fig = plt.figure()
    ax = fig.add_axes([0,0,1,1])

    ax.pie(df['Class'].value_counts(),
           labels = labels,
           explode = explode,
           autopct ='%1.2f%%',
           frame = True,
           textprops = dict(color ="black", size=12))

    ax.axis('equal')
    plt.title('Subcription to the term deposit\n% of Total Clients',
         loc='left',
         color = 'black',
         fontsize = '18')

    plt.show()
    st.pyplot()

#correlation matrix for 9 features
    X = df.drop(['Class', 'job', 'marital', 'education','default','loan', 'campaign', 'previous'], axis=1)
    plt.subplots(figsize=(15,10))
    ax = plt.axes()
    ax.set_title("Marketing Characteristic Heatmap")
    corr = X.corr()
    sns.heatmap(corr,
                xticklabels=corr.columns.values,
                yticklabels=corr.columns.values,
               cmap="Blues")
    plt.show()
    st.subheader("Correlation Matrix")

    st.pyplot()

#Chart for Age
    st.subheader("Age Feature")
    plt.figure(figsize=(19, 9))
    sns.countplot(data=df, x='age', hue='Class')
    st.pyplot()

    st.write("")
# Chart for month
    df_age = df
    st.subheader("Month Feature")
    fig, ax = plt.subplots(figsize = (15, 5))
    sns.countplot(x = 'month', data = df_age, order = ['jan','feb','mar', 'apr', 'may', 'jun', 'jul', 'aug', 'sep', 'oct', 'nov', 'dec'])
    ax.set_xlabel("Months")
    ax.set_ylabel("Count")
    ax.set_title("Count of contacts made in each month")
    st.pyplot()

    df_age.loc[df_age['month']=='jan','month']=1
    df_age.loc[df_age['month']=='feb','month']=2
    df_age.loc[df_age['month']=='mar','month']=3
    df_age.loc[df_age['month']=='apr','month']=4
    df_age.loc[df_age['month']=='may','month']=5
    df_age.loc[df_age['month']=='jun','month']=6
    df_age.loc[df_age['month']=='jul','month']=7
    df_age.loc[df_age['month']=='aug','month']=8
    df_age.loc[df_age['month']=='sep','month']=9
    df_age.loc[df_age['month']=='oct','month']=10
    df_age.loc[df_age['month']=='nov','month']=11
    df_age.loc[df_age['month']=='dec','month']=12

    dict1=dict(list(df_age.groupby(['month','Class'])))
    list1=[1,2,3,4,5,6,7,8,9,10,11,12]
    no=[]
    yes=[]
    months=[]
    for i in list1:
      months.append(i)
      for j in ['no','yes']:
        if(j=='no'):
          no.append(dict1[i,j].count()['Class'])
        else:
          yes.append(dict1[i,j].count()['Class'])

    total_count_per_month=[]
    dict2=dict(list(df.groupby(['month'])))
    for i in list1:
      total_count_per_month.append(dict2[i].count()['Class'])

    month_wise=pd.DataFrame()
    month_wise['Months']=months
    month_wise['Total Enteries per month']=total_count_per_month
    month_wise['Count of Subscribed']=yes
    month_wise['Count of Not Sub']=no
    month_wise['Subscription Rate %']=round((month_wise['Count of Subscribed']/month_wise['Total Enteries per month'])*100)
    month_wise['Not Sub Rate %']=round((month_wise['Count of Not Sub']/month_wise['Total Enteries per month'])*100)
    month_wise=month_wise.sort_values("Months",ascending=True)

    st.write("Based off the chart above, May has the most contact made to the customer (",13766,") but the least amount of subscription rate (",7,"%) compared to March which has the second least contact rate (",477,") but the highest subscription rate (",52,"%)")
    st.dataframe(month_wise)

##Housing
    st.subheader("Housing Feature")
    housing_n=df[df['housing']=='no']
    housing_y=df[df['housing']=='yes']

    total=df.shape[0]
    no=housing_n.count()['Class']
    yes=housing_y.count()['Class']
    sns.countplot(df['housing'])
    st.pyplot()

    st.write(round((yes/total)*100,2),"% of customers, which mean",yes," has a house loan and" ,round((no/total)*100,2),"% do not have a house loan with a count of",no)

    no=housing_y[housing_y['Class']=='no'].count()['Class']
    yes=housing_y[housing_y['Class']=='yes'].count()['Class']
    total=housing_y.count()['Class']
    st.write("Out of the total amount of customers that have a house loan",round((yes/total)*100,2)," % subscribed to term deposit and",round((no/total)*100,2),"% did not subscribe")

##Contact
    st.subheader("Contact Feature")
    sns.countplot(df['contact'],hue=df['Class'])
    st.pyplot()

##duration
    st.subheader("Duration Feature")
    plt.figure(figsize=(20,5))
    plt.xticks(np.arange(0,5000,150))
    sns.scatterplot(df['duration'],df['campaign'],hue=df['Class'])
    st.pyplot()
    st.write("For duration of the calls, if the call had a shorter duration the customers least likely subscribed to the term deposit while when calls lasted longer you can see more customers subscribing.")


text:
# coding: utf-8
#
# This code is part of qclib.
#
# Copyright (c) 2021, Dylan Jones

import numpy as np
from scipy.linalg import expm
from itertools import product
from ..math import kron, P0, P1

# ======================== SINGLE QUBIT GATES =======================

X_GATE = np.array([[0, 1], [1, 0]])
Y_GATE = np.array([[0, -1j], [1j, 0]])
Z_GATE = np.array([[1, 0], [0, -1]])
HADAMARD_GATE = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
PHASE_GATE = np.array([[1, 0], [0, 1j]])
T_GATE = np.array([[1, 0], [0, np.exp(1j*np.pi/4)]])

def id_gate(*args):
    """Functional for the single-qubit identity gate."""
    return np.eye(2)

def x_gate(*args):
    """Functional for the single-qubit Pauli-X gate."""
    return X_GATE

def y_gate(*args):
    """Functional for the single-qubit Pauli-Y gate."""
    return Y_GATE

def z_gate(*args):
    """Functional for the single-qubit Pauli-Z gate."""
    return Z_GATE

def h_gate(*args):
    """Functional for the single-qubit Hadamard (H) gate."""
    return HADAMARD_GATE

def s_gate(*args):
    """Functional for the single-qubit phase (S) gate."""
    return PHASE_GATE

def t_gate(*args):
    """Functional for the single-qubit T gate."""
    return T_GATE

def rx_gate(phi: float = 0):
    """Functional for the single-qubit Pauli-X rotation-gate.

    Parameters
    ----------
    phi : float
        Rotation angle (in radians)

    Returns
    -------
    rx : (2, 2) np.ndarray
    """
    arg = phi / 2
    return np.array([[np.cos(arg), -1j*np.sin(arg)], [-1j*np.sin(arg), np.cos(arg)]])

def ry_gate(phi: float = 0):
    """Functional for the single-qubit Pauli-Y rotation-gate.

    Parameters
    ----------
    phi : float
        Rotation angle (in radians)

    Returns
    -------
    ry : (2, 2) np.ndarray
    """
    arg = phi / 2
    return np.array([[np.cos(arg), -np.sin(arg)], [+np.sin(arg), np.cos(arg)]])

def rz_gate(phi: float = 0):
    """Functional for the single-qubit Pauli-Z rotation-gate.

    Parameters
    ----------
    phi : float
        Rotation angle (in radians)

    Returns
    -------
    rz : (2, 2) np.ndarray
    """
    arg = 1j * phi / 2
    return np.array([[np.exp(-arg), 0], [0, np.exp(arg)]])

def single_gate(qbits, gates, num_qubits=None):
    """Builds the matrix of a single-qubit gate.

    Parameters
    ----------
    qbits : int or Sequence of int
        Index of the qubits.
    gates : list of ndarray or ndarray
        Matrix-representation of the single qubit gate.
        A matrix for each qubit index in 'qbits' must be given.
    num_qubits : int, optional
        Total number of qubits. If not specified number of involved qubits is used.

    Returns
    -------
    gate : np.ndarray
        The full matrix of the single-qubit gate.
    """
    if not hasattr(qbits, "__len__"):
        qbits = [qbits]
        gates = [gates]
    if num_qubits is None:
        num_qubits = max(qbits) + 1

    arrays = np.zeros((num_qubits, 2, 2), dtype=np.complex64)
    idx = [0, 1]
    arrays[:, idx, idx] = 1

    for idx, gate in zip(qbits, gates):
        arrays[idx] = gate
    return kron(arrays)

def _get_projection(*items, num_qubits):
    parts = [np.eye(2)] * num_qubits
    for idx, proj in items:
        parts[idx] = proj
    return parts

def cgate(control, target, gate, num_qubits=None, trigger=1):
    """Builds matrix of a n-qubit control gate.

    Parameters
    ----------
    control : int or list of int
        Index of control-qubit(s)
    target : int
        Index of target qubit
    gate : array_like
        Matrix of the gate that is controlled
    num_qubits : int, optional
        Total number of qubits. If not specified
        number of involved qubits is used
    trigger : int, optional
        Value of control qubits that triggers gate. The default is 1.

    Returns
    -------
    gate : np.ndarray
    """
    if not hasattr(control, "__len__"):
        control = [control]
    num_qubits = num_qubits or max(*control, target) + 1
    gate = np.asarray(gate)
    array = 0
    for vals in product([0, 1], repeat=len(control)):
        # Projections without applying gate
        projections = [P1 if x else P0 for x in vals]
        items = list(zip(control, projections))
        if np.all(vals) == trigger:
            # Projection with applying gate
            items.append((target, gate))
        # Build projection array and add to matrix
        array = array + kron(_get_projection(*items, num_qubits=num_qubits))
    return array

max_stars_repo_name: dylanljones/qclib

max_issues_repo_name: dylanljones/qclib

max_forks_repo_name: dylanljones/qclib


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
The geneticoptimizer module provides optimizers and wrappers for GA.

As for now I didn't tried a convincing Genetic Algorithm library. This is why
you it is built-in here. Feel free to propose a copyleft library on GitHub!

Created on Thu Jun 10 21:20:47 2021.

@author: HugoFara
"""
import os.path
import json
import multiprocessing as mp
import numpy as np
from numpy.random import rand, normal, randint
# Progress bar
import tqdm
from pylinkage.geometry import dist

def tqdm_verbosity(iterable, verbose=True, *args, **kwargs):
    """Wrapper for tqdm, that let you specify if you want verbosity."""
    if verbose:
        for i in tqdm.tqdm(iterable, *args, **kwargs):
            yield i
    else:
        for i in iterable:
            yield i

def kwargs_switcher(arg_name, kwargs, default=None):
    """Simple function to return the good element from a kwargs dict."""
    out = default
    if arg_name in kwargs:
        out = kwargs.pop(arg_name) or out
    return out

def load_population(file_path):
    """Return a population from a given file."""
    with open(file_path, 'r') as file:
        data = json.load(file)
        pop = data[-1]['population']
    return pop

def save_population(file_path, population, verbose=False, data_descriptors=None):
    """
    Save the population to a json file.

    Parameters
    ----------
    file_path : str
    population : list of dna
    verbose : bool
    data_descriptors : dict
        Any additional value you want to save for the current generation.
    """
    if not os.path.exists(file_path):
        data = []
        turn = 0
    else:
        with open(file_path, 'r') as file:
            data = json.load(file)
        turn = data[-1]['turn'] + 1
    data.append({
        'turn': turn,
        'population': population
    })
    if data_descriptors is not None:
        data[-1].update(data_descriptors)
    with open(file_path, 'w') as file:
        json.dump(data, file)
    if verbose:
        print('Data saved.')

def birth(par1, par2, prob):
    """
    Return a new individual with par1 and par2 as parents (two sequences).

    Child are generated by a uniform crossover followed by a random "resetting"
    mutation of each gene. The resetting is a normal law.

    Initial positions come from one of the two parents randomly.

    Parameters
    ----------
    par1 : list[float, tuple of float, tuple of tuple of float]
        Dna of first parent.
    par2 : list[float, tuple of float, tuple of tuple of float]
        Dna of second parent.
    prob : list[float]
        Probability for each gene to mutate, width of a normal law.

    Returns
    -------
    child : list[float, tuple of float, tuple of tuple of float]
        Dna of the child.

    """
    child = [0, [], []]
    for gene1, gene2 in zip(par1[1], par2[1]):
        child[1].append(normal((gene1 if rand() < .5 else gene2), prob))
    for pos1, pos2 in zip(par1[2], par2[2]):
        child[2].append(pos1 if rand() < .5 else pos2)
    return child

def evaluate_individual(dna, fitness, fitness_args):
    """Simple evaluation for a single individual.

    Parameters
    ---------
    dna : list[float, tuple of float, tuple of tuple of float]
        List of the individuals' DNAs
    fitness : callable
        Fitness function of signature fitness(dna, fitness_args) → float.
    fitness_args : tuple
        Additional arguments to pass to the fitness function. Usually the
        initial positions of the joints.

    Returns
    -------
    tuple
        Score then initial coordinates.

    See Also
    --------
    evaluate_population : counterpart for an entire population.
    """
    if fitness_args is None:
        fit = fitness(dna)
    else:
        fit = fitness(dna, *fitness_args)
    if len(fit[1]):
        return fit[0], fit[1]
    # Don't change initial positions for unbuildable individuals.
    return fit[0], dna[2]

def evaluate_population(pop, fitness, fitness_args, verbose=True, processes=1):
    """
    Evaluate the whole population, attribute scores.

    Parameters
    ---------
    pop : list of list[float, tuple of float, tuple of tuple of float]
        List of the individuals' DNAs
    fitness : callable
        Fitness function of signature fitness(dna, fitness_args) → float.
    fitness_args : tuple
        Additional arguments to pass to the fitness function. Usually the
        initial positions of the joints.
    verbose : bool, default=True
        To display informations about population evaluation.
    processes : int, default=1
        Number of processes involved for a multiprocessors evaluation.

    See Also
    --------
    evaluate_individual : same function but on a single DNA.
    """
    # For multiprocessing we load the processes
    if processes > 1:
        res = [None] * len(pop)
        with mp.Pool(processes=processes) as pool:
            # Load the processes
            for i, dna in enumerate(pop):
                res[i] = pool.apply_async(
                    evaluate_individual, (dna, fitness, fitness_args)
                )
            # Then get data
            for result, dna in zip(res, pop):
                dna[0], dna[2] = result.get()
    else:
        for dna in pop:
            dna[0], dna[2] = evaluate_individual(dna, fitness, fitness_args)

    if verbose:
        diversity = np.linalg.norm(
            np.var([dna[1] for dna in pop], axis=0)
        )
        print("Scores:", [dna[1] for dna in pop])
        print("Genetic diversity: ", diversity)

def select_parents(pop, verbose=True):
    """Selection 1/4 of the population as parents."""
    median = np.median([dna[0] for dna in pop])
    # Index of best individual
    best_index, best_dna = max(enumerate(pop), key=lambda x: x[1][0])
    # Parents selection, 1/4 of population
    parents = []
    indexes = []
    for j, individual in enumerate(pop):
        # Parents whose score is above median.
        # Individuals with best fitness are more likely to be selected
        if (
                .5 * (individual[0] - best_dna[0]) / (best_dna[0] - median)
        ) + 1 > max(rand(), .5):
            parents.append(individual)
            indexes.append(j)
    # Add best individual if needed
    if best_index not in indexes:
        parents.insert(0, best_dna)
        indexes.append(best_index)
    # Add a random parent if odd number
    if len(parents) % 2:
        for j, individual in enumerate(pop):
            if j not in indexes:
                parents.append(individual)
                indexes.append(j)
                break
    if verbose:
        print(f"Median score: {median}, {len(parents)} parents")
    return parents

def make_children(parents, prob, max_genetic_dist=float('inf')):
    children = []
    j = 0
    while len(parents) > 1 and j < 100:
        par1 = parents.pop(randint(len(parents) - 1))
        if len(parents) > 1:
            par2 = parents.pop(randint(len(parents) - 1))
        else:
            par2 = parents.pop()
        if dist(par1[1], par2[1]) < max_genetic_dist:
            children.append(birth(par1, par2, prob))
        elif parents:
            parents.append(par1)
            parents.append(par2)
        j += 1
    return children

def evolutionary_optimization_builtin(
        dna,
        prob,
        fitness,
        iters,
        **kwargs
):
    """
    Optimization by genetic algorithm (GA).

    Parameters
    ----------
    dna : list[float, tuple of float, tuple of tuple of float]
        DNA of individuals.
    prob : list of float
        List of probabilities of the good transmission of one
        characteristics.
    fitness : callable
        Evaluation function.
        Return float.
        The
    iters : int
        Number of iterations.
    **kwargs : dict
        Other useful parameters for the optimization.

        max_pop : int, default=11
            Maximum number of individuals. The default is 11.
        init_pop : sequence of object, default=None
            Initial population, for wider INITIAL genetic diversity.
            The default is None.
        max_genetic_dist : float, default=.7
            Maximum genetic distance, before individuals
            cannot reproduce (separated species). The default is .7.
        startnstop : bool, default=False
            Ability to close program without loosing population.
            If True, we verify at initialization the existence of a data file.
            Population is save every int(250 / max_pop) iterations.
            The default is False.
        fitness_args : tuple
            Keyword arguments to send to the fitness function.
            The default is None (no argument sent).
        verbose : int
            Level of verbosity.
            0 : no verbose, do not print anything.
            1 : show a progress bar.
            2 : complete report for each turn.
            The default is 1.
        processes : int, default=1
            Number of processes that will evaluate the linkages.

    Returns
    -------
    list[float, tuple of float, tuple of tuple of float]
        List of 3-tuples: best dimensions, best score and initial positions.
        The list is sorted by score order.
    """
    file_path = 'Population evolution.json'
    startnstop = kwargs_switcher('startnstop', kwargs, False)
    if startnstop and os.path.exists(file_path):
        pop = load_population(file_path)
    else:
        # At least two parents to begin with
        pop = [[dna[0], list(dna[1]), list(dna[2])] for _ in range(2)]

    max_pop = kwargs_switcher('max_pop', kwargs, 11)
    max_genetic_dist = kwargs_switcher('max_genetic_dist', kwargs, .7)
    verbose = kwargs_switcher('verbose', kwargs, 1)
    fitness_args = kwargs_switcher('fitness_args', kwargs, None)
    # Number of evaluations to run in parallel
    processes = kwargs_switcher('processes', kwargs, default=1)
    # "Garden of Eden" phase, add enough children to get as much individuals as
    # required
    init_pop = kwargs_switcher('init_pop', kwargs, default=max_pop)
    for i in range(len(pop), init_pop):
        pop.append(
            birth(
                pop[randint(len(pop) - 1)],
                pop[randint(len(pop) - 1)],
                prob
            )
        )
    postfix = [
        "best_score", max(x[0] for x in pop),
        "best_dimensions", max(pop, key=lambda x: x[0])[1]
    ]
    iterations = tqdm_verbosity(
        range(iters),
        verbose=verbose == 1,
        total=iters,
        desc='Evolutionary optimization',
        postfix=postfix
    )
    for i in iterations:
        if verbose > 1:
            print(f"Turn: {i}, {len(pop)} individuals.")
        # Individuals evaluation
        evaluate_population(
            pop, fitness, fitness_args,
            verbose=verbose > 1,
            processes=processes
        )
        # Population selection
        # Minimal score before death
        death_score = np.quantile([j[0] for j in pop], 1 - max_pop / len(pop))
        if np.isnan(death_score):
            death_score = - float('inf')
        # We only keep max_pop individuals
        pop = list(filter(lambda x: x[0] >= death_score, pop))
        parents = select_parents(pop, verbose=verbose > 1)
        # We select the best fit individual to show off, we now he is a parent
        best_id = max(enumerate(parents), key=lambda x: x[1][0])[0]
        postfix[1] = parents[best_id][0]
        postfix[3] = parents[best_id][1]
        if startnstop:
            save_population(
                file_path, pop, verbose > 1,
                {
                    'best_score': parents[best_id][0],
                    'best_individual_id': best_id
                }
            )
        # Children generation
        children = make_children(parents, prob, max_genetic_dist)
        # Add to population
        pop.extend(children)

    out = []
    for dna in pop:
        fit = fitness(dna)
        if isinstance(fit, tuple):
            out.append((fit[0], fit[1], dna[2]))
        else:
            out.append((fit, dna[1], dna[2]))
    out.sort(key=lambda x: x[0], reverse=True)
    return out

def evolutionary_optimization(
        dna,
        fitness,
        iters,
        prob=.07,
        **kwargs
):
    """
    Run the Genetic Optimizer.

    Genetic Optimization is a procedural algorithm based on Darwinian evolution
    models.

    As of today, this function will only use the built-in algorithm, but you
    should still use it because we may implement another GA library.

    Parameters
    ----------
    dna : list[float, tuple of float, tuple of tuple of float]
        DNA of a linkage in format (score, dimensions, initial coordinates).
    prob : float or list of float, default=.07
        Mutation probability for each gene.
    fitness : callable
        Evaluation function for an MAXIMISATION problem.
        Must return a float.
    iters : int
        Number of iterations.
    **kwargs : dict
        Other useful parameters for the optimization.

        max_pop : int, optional
            Maximum number of individuals. The default is 11.
        init_pop : sequence of object, optional
            Initial population, for wider INITIAL genetic diversity.
            The default is None.
        max_genetic_dist : float, optional
            Maximum genetic distance, before individuals
            cannot reproduce (separated species). The default is .7.
        startnstop : bool, optional
            Ability to close program without loosing population.
            If True, we verify at initialization the existence of a data file.
            Population is save every int(250 / max_pop) iterations.
            The default is False.
        fitness_args : tuple
            Positional arguments to send to the fitness function.
            The default is None (no argument sent).
        verbose : int, optional
            Level of verbosity.
            0 : no verbose, do not print anything.
            1 : show a progress bar.
            2 : complete report for each turn.
            The default is 1.
        processes : int, default=1
            Number of processes that will evaluate the linkages.

    Returns
    -------
    list[dna]
        An iterable of the best fit individuals, in format
        (score, dimensions, initial coordinates).

    See Also
    --------
    evolutionary_optimization_builtin : built-in genetic algorithm.

    """
    # Legacy fallback
    return evolutionary_optimization_builtin(
        dna, prob, fitness,
        iters=iters,
        **kwargs
    )

max_stars_repo_name: HugoFara/leggedsnake

max_issues_repo_name: HugoFara/leggedsnake

max_forks_repo_name: HugoFara/leggedsnake

reason: import numpy,from numpy


text:
#!/usr/bin/env python3

import sys
import math
import pprint
import os
import re
import doctest
import itertools
import types
import logging
from collections import deque
from collections import defaultdict
#import networkx as nx
from copy import deepcopy
try:
   import matplotlib.pyplot as plt
except ImportError:
   plt = None

# create absolute mydir
mydir = os.path.abspath(os.path.dirname(__file__))
sys.path.append(os.path.join(mydir, '../lib'))

from advent import *
from intcode import *

def print_seats(seats):
    for s in seats:
        print(''.join(s))
    print()
    print()

def get_seat(seats, pos):
    # print(pos)
    if pos[0] < 0 or pos[0] >= len(seats):
        return '.'
    if pos[1] < 0 or pos[1] >= len(seats[0]):
        return '.'
    return seats[pos[0]][pos[1]]

def alldirs():
    for x in [-1, 0, 1]:
        for y in [-1, 0, 1]:
            p = (x, y)
            if p != (0, 0):
                yield p

def process_seats(seats):
    newseats = deepcopy(seats)
    occupied = 0
    for y, row in enumerate(seats):
        for x, s in enumerate(row):
            mypos = (y, x)
            if s == '.':
                continue
            adj = sum(1 if '#' == get_seat(seats, addpos(mypos, p)) else 0 for p in alldirs())
            sval = s
            if s == 'L':
                if adj == 0:
                    sval = '#'
            else:
                # print(adj, mypos)
                sval = 'L' if adj >= 4 else '#'
            occupied += 1 if sval == '#' else 0
            newseats[y][x] = sval
    return occupied, newseats

if __name__ == '__main__':
    if len(sys.argv) == 2 and sys.argv[1] == "TEST":
        import doctest
        doctest.testmod()
        sys.exit(0)

    logging.basicConfig(level=logging.INFO)
    path = "day11.txt"
    if len(sys.argv) > 1:
        path = sys.argv[1]

    with open(path) as f:
        seats = [list(s.strip()) for s in f.readlines()]

    c = 0
    while True:
        print(c)
        print_seats(seats)

        c2, seats = process_seats(seats)
        if c == c2:
            break
        c = c2

reason: import networkx


text:
import torch
import os
import datasets
import models
import pickle
import time
import numpy
import utils
import loss
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d

device = 'cuda' if torch.cuda.is_available() else 'cpu'

# load parameters of a learned model
loaddir = sys.argv[1]
savedir = sys.argv[2]

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

with open(loaddir + 'args.pkl', 'rb') as f:
    variables = pickle.load(f)
args = variables[0]

params = torch.load(loaddir + 'cnn_microscopy.pth', map_location=device)

model = models.Network().to(device)
model.load_state_dict(params)
model.eval()

# load real data as validataset
validateset = datasets.Tubulin_191216()
# validateset = datasets.Tubulin_191216_masked()
validate_loader = torch.utils.data.DataLoader(validateset, batch_size=50, shuffle=False)

# size of a real image
image_size = torch.Tensor([256, 256, 4]).to(device)
# minimum coordinates of each voxel of input image
defaults = utils.gen_defaults(1/image_size, device, centering=False)
lossclass = loss.LocationAndConfidenceLoss(defaults, 1/image_size, default_centering=False, negapos_ratio=3).to(device)
# thresholding value for binarizing classification results
thresh = 0.5

start_time = time.time()
# output = numpy.array([])
particles = numpy.empty(shape=[0, 3])
num_particles = []

# size of the target space
target_space_size = torch.Tensor([
    256 * 192, 256 * 192, 4 * 400]).to(device)
# the minimum coordinate of the target space
target_space_min = torch.Tensor([0, 0, -200]).to(device)

with torch.no_grad():
    for i, data in enumerate(validate_loader):
        if (i+1) % 10 == 0:
            print(i+1)
            end_time = time.time()
            tdiff = end_time - start_time
            print('elapsed time:', tdiff)
            print('fps={0:.1f}'.format(i*validate_loader.batch_size/tdiff))
        data = data.to(device)
        predictions = model(data).to(float)

        # rescale predictions
        predictions[:, :, :3] = (predictions[:, :, :3] * (1/image_size) + defaults) *target_space_size + target_space_min

        # save predicted coordinates as a text file
        f = open(savedir + 'coord_{0:.2f}.txt'.format(thresh), 'a')
        for pred in predictions:
            coordinates = pred[pred[:, 3] > thresh].cpu().numpy()
            numpy.savetxt(f, coordinates, delimiter=',')
        f.close()

end_time = time.time()
tdiff = end_time - start_time
print('elapsed time:', tdiff)
print('fps={0:.1f}'.format(len(validate_loader.dataset)/tdiff))

max_stars_repo_name: t-aritake/voxel-wise-regression-SMLM

max_issues_repo_name: t-aritake/voxel-wise-regression-SMLM

max_forks_repo_name: t-aritake/voxel-wise-regression-SMLM


text:
import csv
import matplotlib.animation as animation
'''
timestamo: data[11]
linearacceleration x,y,z[12,13,14]
'''

import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

def randrange(n, vmin, vmax):
    return (vmax-vmin)*np.random.rand(n) + vmin

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
# n = 100
# for c, m, zl, zh in [('r', 'o', -50, -25), ('b', '^', -30, -5)]:
#     xs = randrange(n, 23, 32)
#     ys = randrange(n, 0, 100)
#     zs = randrange(n, zl, zh)
#     ax.scatter(xs, ys, zs, c=c, marker=m)

# ax.set_xlabel('X Label')
# ax.set_ylabel('Y Label')
# ax.set_zlabel('Z Label')
#
# plt.show()

with open('signal1_data.csv','rb') as f_in:
    reader = csv.reader(f_in)
    raw_data = list(reader)

# print data[0]

for index,data in enumerate(raw_data):
    if index==0:
        continue
    for index2,k in enumerate(data):
        raw_data[index][index2]=float(raw_data[index][index2])

v_x=0
v_y=0
v_z=0
s_x=0
s_y=0
s_z=0

sports_result=list()
sports_x=[]
sports_y=[]
sports_z=[]
moving_time=[]
a_x=[]
a_y=[]
a_z=[]
v_x_list=[]
v_y_list=[]
v_z_list=[]

for index,data in enumerate(raw_data):
    if index==0 or index==1:
        continue
    time=data[11] if index==1 else data[11]-raw_data[index-1][11]
    time_x=s_x+v_x*time+(0.5)*raw_data[index-1][12]*(time**2)
    time_y=s_y+v_y*time+(0.5)*raw_data[index-1][13]*(time**2)
    time_z=s_z+v_z*time+(0.5)*raw_data[index-1][14]*(time**2)
    sports_result.append([data[11],time_x,time_y,time_z])

    a_x.append(raw_data[index-1][12])
    a_y.append(raw_data[index-1][13])
    a_z.append(raw_data[index-1][14])

    v_x_list.append(v_x)
    v_y_list.append(v_y)
    v_z_list.append(v_z)

    sports_x.append(time_x*10)
    sports_y.append(time_y*10)
    sports_z.append(time_z*10)
    moving_time.append(data[11])
    # if time_x<s_x or time_y<s_y or time_z<s_z:
    #     print "yes"
    s_x=time_x
    s_y=time_y
    s_z=time_z
    v_x+=raw_data[index-1][12]*time
    v_y+=raw_data[index-1][13]*time
    v_z+=raw_data[index-1][14]*time
    # v_x=0
    # v_y=0
    # v_z=0
# print moving_time

def data_gen():
    t = data_gen.t
    cnt = 0
    while cnt < len(sports_x):
        cnt+=1
        t += 0.05
        print "postion of x, y, z: ",sports_x[cnt],sports_y[cnt],sports_z[cnt]
        print "time: ",moving_time[cnt]
        print "x accrleration: ",a_x[cnt]," y accrleration: ",a_y[cnt]," z accrleration: ",a_z[cnt]
        print "x speed: ",v_x_list[cnt]," y speed: ",v_y_list[cnt]," z speed: ",v_z_list[cnt]
        yield sports_x[cnt],sports_y[cnt],sports_z[cnt]
    print "end"
data_gen.t = 0

xdata=[]
ydata=[]
zdata=[]

def run(data):
    # update the data
    x,y,z = data
    # xdata.append(x)
    # ydata.append(y)
    # zdata.append(z)
    # line.set_data(xdata, ydata)
    ax.scatter([x],[y],[z],color='blue')

# ani = animation.FuncAnimation(fig, run, data_gen, blit=False, interval=1,
#     repeat=False)
# ax.set_xlabel('X Label')
# ax.set_ylabel('Y Label')
# ax.set_zlabel('Z Label')
# plt.show()

ax.scatter(sports_x, sports_y, sports_z,color='blue')
ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')

plt.show()
# print "end"
# print sports_result
# import json
# out_put={'data':sports_result}
#
# with open('moving_track_1.json','wb') as f_out:
#     f_out.write(json.dumps(out_put))

max_stars_repo_name: luckyharryji/physical-modeling-band

max_issues_repo_name: luckyharryji/physical-modeling-band

max_forks_repo_name: luckyharryji/physical-modeling-band


text:
import os
import numpy as np
import sys
from PIL import Image
import shutil
base_dir = sys.argv[1]
im_dir = os.path.join(base_dir, 'Images')
label_dir = os.path.join(base_dir, 'Labels')
list_file = os.path.join(base_dir, 'all.txt')
fds = [fd for fd in os.listdir(base_dir) if os.path.isdir(os.path.join(base_dir, fd))]
list_objs=[]
for fd in fds:
    print fd
    sub_fds = os.listdir(os.path.join(base_dir, fd, 'data'))
    sub_fds = sorted(sub_fds)
    for sub_fd in sub_fds:
        list_objs.append(fd+'_'+sub_fd)
        im_list = os.listdir(os.path.join(base_dir, fd,'data',sub_fd,'shots','001','images'))
        new_im_fd = os.path.join(im_dir, fd+'_'+sub_fd)
        new_label_fd = os.path.join(label_dir, fd+'_'+sub_fd)
        if not os.path.exists(new_im_fd):
            os.makedirs(new_im_fd)
            os.makedirs(new_label_fd)
        for im_name in im_list:
            shutil.copy(os.path.join(base_dir, fd,'data',sub_fd,'shots','001','images',im_name),
                    os.path.join(new_im_fd, im_name))
            label = Image.open(os.path.join(base_dir, fd,'data',sub_fd,'shots','001','labels',im_name[:-4]+'.jpg'))
            label = Image.fromarray(((np.array(label)>127) * 255).astype(np.uint8))
            label.save(os.path.join(new_label_fd, im_name))
with open(list_file,'w') as f:
    f.write('\n'.join(list_objs))