text:
"""
Response from a Homogeneous Layer for Different Waveforms
=========================================================

Here we use the module *SimPEG.electromagnetics.viscous_remanent_magnetization*
to predict the characteristic VRM response over magnetically viscous layer.
We consider a small-loop, ground-based survey which uses a coincident loop
geometry. For this tutorial, we focus on the following:

    - How to define the transmitters and receivers
    - How to define the survey
    - How to define a diagnostic physical property
    - How to define the physics for the linear potential fields formulation
    - How the VRM response depends on the transmitter waveform

Note that for this tutorial, we are only modeling the VRM response. A separate
tutorial have been developed for modeling both the inductive and VRM responses.

"""

#########################################################################
# Import Modules
# --------------
#

import SimPEG.electromagnetics.viscous_remanent_magnetization as vrm

from discretize import TensorMesh
from discretize.utils import mkvc

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl

# sphinx_gallery_thumbnail_number = 2

##########################################################################
# Define Waveforms
# ----------------
#
# Under *SimPEG.electromagnetic.viscous_remanent_magnetization.waveform*
# there are a multitude of waveforms that can be defined (Step-off, square-pulse,
# piecewise linear, ...). Here we define a specific waveform for each transmitter.
# Each waveform is defined with a diferent set of parameters.
#

waveform_list = []

# Step-off waveform
waveform_list.append(vrm.waveforms.StepOff(t0=0))

# 20 ms square pulse with off-time at t = 0 s.
waveform_list.append(vrm.waveforms.SquarePulse(t0=0, delt=0.02))

# 30 ms trapezoidal waveform with off-time at t = 0 s.
t_wave = np.r_[-0.03, -0.02, -0.01, 0]
I_wave = np.r_[0.0, 1.0, 1.0, 0]
waveform_list.append(vrm.waveforms.ArbitraryPiecewise(t_wave=t_wave, I_wave=I_wave))

# 40 ms triangular waveform with off-time at t = 0 s.
t_wave = np.r_[-0.04, -0.02, 0]
I_wave = np.r_[0.0, 1.0, 0]
waveform_list.append(vrm.waveforms.ArbitraryPiecewise(t_wave=t_wave, I_wave=I_wave))

# Plot waveforms
fig = plt.figure(figsize=(8, 4))
mpl.rcParams.update({"font.size": 12})
ax1 = fig.add_axes([0.1, 0.1, 0.85, 0.85])
ax1.plot(np.r_[-0.04, 0.0, 0.0, 0.02], np.r_[1, 1, 0, 0], "b", lw=2)
ax1.plot(np.r_[-0.04, -0.02, -0.02, 0.0, 0.0, 0.04], np.r_[0, 0, 1, 1, 0, 0], "r", lw=2)
ax1.plot(np.r_[-0.04, -0.03, -0.02, -0.01, 0, 0.04], np.r_[0, 0, 1, 1, 0, 0], "k", lw=2)
ax1.plot(np.r_[-0.04, -0.02, 0, 0.04], np.r_[0, 1, 0, 0], "g", lw=2)
ax1.set_xlim((-0.04, 0.04))
ax1.set_ylim((-0.01, 1.1))
ax1.set_xlabel("time [s]")
ax1.set_ylabel("current [A]")
ax1.set_title("Waveforms")
ax1.legend(
    ["step-off", "20 ms square pulse", "30 ms trapezoidal", "40 ms triangular"],
    loc="upper right",
)

##########################################################################
# Survey
# ------
#
# Here we define the sources, the receivers and the survey. For this exercise,
# we are modeling the response for single transmitter-receiver pair with
# different transmitter waveforms.
#

# Define the observation times for the receivers. It is VERY important to
# define the first time channel AFTER the off-time.
time_channels = np.logspace(-4, -1, 31)

# Define the location of the coincident loop transmitter and receiver.
# In general, you can define the receiver locations as an (N, 3) numpy array.
xyz = np.c_[0.0, 0.0, 0.5]

# There are 4 parameters needed to define a receiver.
dbdt_receivers = [
    vrm.receivers.Point(xyz, times=time_channels, fieldType="dbdt", orientation="z")
]

# Define sources
source_list = []
dipole_moment = [0.0, 0.0, 1]
for pp in range(0, len(waveform_list)):

    # Define the transmitter-receiver pair for each waveform
    source_list.append(
        vrm.sources.MagDipole(
            dbdt_receivers, mkvc(xyz), dipole_moment, waveform_list[pp]
        )
    )

# Define the survey
survey = vrm.Survey(source_list)

##########################################################################
# Defining the Mesh
# -----------------
#
# Here we create the tensor mesh that will be used for this tutorial example.
# We are modeling the response from a magnetically viscous layer. As a result,
# we do not need to model the Earth at depth. For this example the layer is
# 10 m thick.
#

cs, ncx, ncy, ncz, npad = 2.0, 35, 35, 5, 5
hx = [(cs, npad, -1.3), (cs, ncx), (cs, npad, 1.3)]
hy = [(cs, npad, -1.3), (cs, ncy), (cs, npad, 1.3)]
hz = [(cs, ncz)]
mesh = TensorMesh([hx, hy, hz], "CCN")

##########################################################################
# Defining an Amalgamated Magnetic Property Model
# -----------------------------------------------
#
# For the linear potential field formulation, the magnetic viscosity
# characterizing each cell can be defined by an "amalgamated magnetic property"
# (see Cowan, 2016). Here we define an amalgamated magnetic property model.
# For other formulations of the forward simulation, you may define the parameters
# assuming a log-uniform or log-normal distribution of time-relaxation constants.
#

# Amalgamated magnetic property for the layer
model_value = 0.0001
model = model_value * np.ones(mesh.nC)

# Define the active cells. These are the cells that exhibit magnetic viscosity
# and/or lie below the surface topography.
ind_active = np.ones(mesh.nC, dtype="bool")

##########################################################################
# Define the Simulation
# ---------------------
#
# Here we define the formulation for solving Maxwell's equations. We have chosen
# to model the off-time VRM response. There are two important keyword arguments,
# *refinement_factor* and *refinement_distance*. These are used to refine the
# sensitivities of the cells near the transmitters. This improves the accuracy
# of the forward simulation without having to refine the mesh near transmitters.
#

# For this example, cells lying within 2 m of a transmitter will be modeled
# as if they are comprised of 4^3 equal smaller cells. Cells within 4 m of a
# transmitter will be modeled as if they are comprised of 2^3 equal smaller
# cells.
simulation = vrm.Simulation3DLinear(
    mesh,
    survey=survey,
    indActive=ind_active,
    refinement_factor=2,
    refinement_distance=[2.0, 4.0],
)

#######################################
# Predict Data and Plot
# ---------------------
#

# Predict VRM response
dpred = simulation.dpred(model)

# Reshape for plotting
n_times = len(time_channels)
n_waveforms = len(waveform_list)
dpred = np.reshape(dpred, (n_waveforms, n_times)).T

# Characteristic VRM decay for several waveforms.
fig = plt.figure(figsize=(6, 7))
ax1 = fig.add_axes([0.15, 0.1, 0.8, 0.85])
ax1.loglog(time_channels, -dpred[:, 0], "b", lw=2)
ax1.loglog(time_channels, -dpred[:, 1], "r", lw=2)
ax1.loglog(time_channels, -dpred[:, 2], "k", lw=2)
ax1.loglog(time_channels, -dpred[:, 3], "g", lw=2)
ax1.set_xlim((np.min(time_channels), np.max(time_channels)))
ax1.set_xlabel("time [s]")
ax1.set_ylabel("-dBz/dt [T/s]")
ax1.set_title("Characteristic Decay")
ax1.legend(
    ["step-off", "20 ms square pulse", "30 ms trapezoidal", "40 ms triangular"],
    loc="upper right",
)

max_stars_repo_name: Prithwijit-Chak/simpeg

max_issues_repo_name: thast/simpeg

max_forks_repo_name: thast/simpeg


text:
import rospy
import numpy
from gym import spaces
from openai_ros.robot_envs import sumitxl_env
from gym.envs.registration import register
from geometry_msgs.msg import Vector3
from geometry_msgs.msg import Point
from tf.transformations import euler_from_quaternion

# The path is __init__.py of openai_ros, where we import the SumitXlMazeEnv directly
timestep_limit_per_episode = 10000 # Can be any Value

register(
        id='SumitXlRoom-v0',
        entry_point='openai_ros:task_envs.sumit_xl.sumit_xl_room.SumitXlRoom',
        timestep_limit=timestep_limit_per_episode,
    )

class SumitXlRoom(sumitxl_env.SumitXlEnv):
    def __init__(self):
        """
        This Task Env is designed for having the sumit_xl in the sumit_xl world
        closed room with columns.
        It will learn how to move around without crashing.
        """

        # Only variable needed to be set here
        number_actions = rospy.get_param('/sumit_xl/n_actions')
        self.action_space = spaces.Discrete(number_actions)

        # We set the reward range, which is not compulsory but here we do it.
        self.reward_range = (-numpy.inf, numpy.inf)

        # Actions and Observations
        self.linear_forward_speed = rospy.get_param('/sumit_xl/linear_forward_speed')
        self.linear_turn_speed = rospy.get_param('/sumit_xl/linear_turn_speed')
        self.angular_speed = rospy.get_param('/sumit_xl/angular_speed')
        self.init_linear_forward_speed = rospy.get_param('/sumit_xl/init_linear_forward_speed')
        self.init_linear_turn_speed = rospy.get_param('/sumit_xl/init_linear_turn_speed')

        self.new_ranges = rospy.get_param('/sumit_xl/new_ranges')
        self.min_range = rospy.get_param('/sumit_xl/min_range')
        self.max_laser_value = rospy.get_param('/sumit_xl/max_laser_value')
        self.min_laser_value = rospy.get_param('/sumit_xl/min_laser_value')
        self.max_linear_aceleration = rospy.get_param('/sumit_xl/max_linear_aceleration')

        self.max_distance = rospy.get_param('/sumit_xl/max_distance')

        # Get Desired Point to Get
        self.desired_point = Point()
        self.desired_point.x = rospy.get_param("/sumit_xl/desired_pose/x")
        self.desired_point.y = rospy.get_param("/sumit_xl/desired_pose/y")
        self.desired_point.z = rospy.get_param("/sumit_xl/desired_pose/z")

        # We create the arrays for the laser readings
        # We also create the arrays for the odometry readings
        # We join them toeguether.
        laser_scan = self._check_laser_scan_ready()
        num_laser_readings = len(laser_scan.ranges)/self.new_ranges
        high_laser = numpy.full((num_laser_readings), self.max_laser_value)
        low_laser = numpy.full((num_laser_readings), self.min_laser_value)

        # We place the Maximum and minimum values of the X,Y and YAW of the odometry
        # The odometry yaw can be any value in the circunference.
        high_odometry = numpy.array([self.max_distance, self.max_distance, 3.14])
        low_odometry = numpy.array([-1*self.max_distance, -1*self.max_distance, -1*3.14])

        # We join both arrays
        high = numpy.concatenate([high_laser, high_odometry])
        low = numpy.concatenate([low_laser, low_odometry])

        self.observation_space = spaces.Box(low, high)

        rospy.logdebug("ACTION SPACES TYPE===>"+str(self.action_space))
        rospy.logdebug("OBSERVATION SPACES TYPE===>"+str(self.observation_space))

        # Rewards
        self.closer_to_point_reward = rospy.get_param("/sumit_xl/closer_to_point_reward")
        self.not_ending_point_reward = rospy.get_param("/sumit_xl/not_ending_point_reward")

        self.end_episode_points = rospy.get_param("/sumit_xl/end_episode_points")

        self.cumulated_steps = 0.0

        # Here we will add any init functions prior to starting the MyRobotEnv
        super(SumitXlRoom, self).__init__()

    def _set_init_pose(self):
        """Sets the Robot in its init pose
        """
        self.move_base( self.init_linear_forward_speed,
                        self.init_linear_turn_speed,
                        epsilon=0.05,
                        update_rate=10)

        return True

    def _init_env_variables(self):
        """
        Inits variables needed to be initialised each time we reset at the start
        of an episode.
        :return:
        """
        # For Info Purposes
        self.cumulated_reward = 0.0
        # Set to false Done, because its calculated asyncronously
        self._episode_done = False

        odometry = self.get_odom()
        self.previous_distance_from_des_point = self.get_distance_from_desired_point(odometry.pose.pose.position)

    def _set_action(self, action):
        """
        This set action will Set the linear and angular speed of the SumitXl
        based on the action number given.
        :param action: The action integer that set s what movement to do next.
        """

        rospy.logdebug("Start Set Action ==>"+str(action))
        # We convert the actions to speed movements to send to the parent class CubeSingleDiskEnv
        if action == 0: #FORWARD
            linear_speed = self.linear_forward_speed
            angular_speed = 0.0
            self.last_action = "FORWARDS"
        elif action == 1: #LEFT
            linear_speed = self.linear_turn_speed
            angular_speed = self.angular_speed
            self.last_action = "TURN_LEFT"
        elif action == 2: #RIGHT
            linear_speed = self.linear_turn_speed
            angular_speed = -1*self.angular_speed
            self.last_action = "TURN_RIGHT"
        """
        elif action == 3: #STOP
            linear_speed = 0.0
            angular_speed = 0.0
            self.last_action = "STOP"
        """

        # We tell SumitXL the linear and angular speed to set to execute
        self.move_base(linear_speed, angular_speed, epsilon=0.05, update_rate=10)

        rospy.logdebug("END Set Action ==>"+str(action))

    def _get_obs(self):
        """
        Here we define what sensor data defines our robots observations
        To know which Variables we have acces to, we need to read the
        SumitXlEnv API DOCS

        WALL CLOSE LEFT [1, 1, 9, 0, 0, 0, -1.8, 0.46, 0.01]
        WALL CLOSE RIGHT [0, 0, 0, 10, 1, 2, -1.8, -0.61, 0.01]
        WALL BACK [0, 9, 1, 1, 6, 0, -1.8, -0.54, 1.59]
        WALL FRONT [2, 9, 0, 0, 2, 2, -1.83, 0.51, 1.58]

        0 in reality is arround front 0.4, back 0.5, sides 0.3

        :return:
        """
        rospy.logdebug("Start Get Observation ==>")
        # We get the laser scan data
        laser_scan = self.get_laser_scan()

        discretized_laser_scan = self.discretize_scan_observation(    laser_scan,
                                                                        self.new_ranges
                                                                        )
        # We get the odometry so that SumitXL knows where it is.
        odometry = self.get_odom()
        x_position = odometry.pose.pose.position.x
        y_position = odometry.pose.pose.position.y

        # We get the orientation of the cube in RPY
        roll, pitch, yaw = self.get_orientation_euler()
        # We round to only two decimals to avoid very big Observation space
        odometry_array = [  round(x_position, 2),
                            round(y_position, 2),
                            round(yaw, 2)]

        # We only want the X and Y position and the Yaw

        observations = discretized_laser_scan + odometry_array

        rospy.logdebug("Observations==>"+str(observations))
        rospy.logdebug("END Get Observation ==>")

        return observations

    def _is_done(self, observations):

        if self._episode_done:
            rospy.logerr("SumitXl is Too Close to wall==>")
        else:
            rospy.logdebug("SumitXl is NOT close to a wall ==>")

        # Now we check if it has crashed based on the imu
        imu_data = self.get_imu()
        linear_acceleration_magnitude = self.get_vector_magnitude(imu_data.linear_acceleration)
        if linear_acceleration_magnitude > self.max_linear_aceleration:
            rospy.logerr("SumitXl Crashed==>"+str(linear_acceleration_magnitude)+">"+str(self.max_linear_aceleration))
            self._episode_done = True
        else:
            rospy.logerr("DIDNT crash SumitXl ==>"+str(linear_acceleration_magnitude)+">"+str(self.max_linear_aceleration))

        current_position = Point()
        current_position.x = observations[-3]
        current_position.y = observations[-2]
        current_position.z = 0.0

        if abs(current_position.x) <= self.max_distance:
            if abs(current_position.y) <= self.max_distance:
                rospy.logdebug("SummitXL Position is OK ==>["+str(current_position.x)+","+str(current_position.y)+"]")
            else:
                rospy.logerr("SummitXL to Far in Y Pos ==>"+str(current_position.x))
                self._episode_done = True
        else:
            rospy.logerr("SummitXL to Far in X Pos ==>"+str(current_position.x))
            self._episode_done = True

        if self.is_in_desired_position(current_position):
            self._episode_done = True

        return self._episode_done

    def _compute_reward(self, observations, done):
        """
        We give reward to the robot when it gets closer to the desired point.
        We Dont give it contsnatly, but only if there is an improvement
        """

        # We get the current Position from the obervations
        current_position = Point()
        current_position.x = observations[-3]
        current_position.y = observations[-2]
        current_position.z = 0.0

        distance_from_des_point = self.get_distance_from_desired_point(current_position)

        distance_difference =  distance_from_des_point - self.previous_distance_from_des_point

        rospy.logwarn("current_position=" + str(current_position))
        rospy.logwarn("desired_point=" + str(self.desired_point))

        rospy.logwarn("total_distance_from_des_point=" + str(self.previous_distance_from_des_point))
        rospy.logwarn("distance_from_des_point=" + str(distance_from_des_point))
        rospy.logwarn("distance_difference=" + str(distance_difference))

        if not done:
            # If there has been a decrease in the distance to the desired point, we reward it
            if distance_difference < 0.0:
                rospy.logwarn("DECREASE IN DISTANCE GOOD")
                reward = self.closer_to_point_reward
            else:
                # If it didnt get closer, we give much less points in theory
                # This should trigger the behaviour of moving towards the point
                rospy.logwarn("NO DECREASE IN DISTANCE, so much less points")
                reward = self.not_ending_point_reward
        else:
            if self.is_in_desired_position(current_position):
                reward = self.end_episode_points
                rospy.logwarn("GOT TO DESIRED POINT ; DONE, reward=" + str(reward))
            else:
                reward = -1*self.end_episode_points
                rospy.logerr("SOMETHING WENT WRONG ; DONE, reward=" + str(reward))

        self.previous_distance_from_des_point = distance_from_des_point

        rospy.logwarn("reward=" + str(reward))
        self.cumulated_reward += reward
        rospy.logdebug("Cumulated_reward=" + str(self.cumulated_reward))
        self.cumulated_steps += 1
        rospy.logdebug("Cumulated_steps=" + str(self.cumulated_steps))

        return reward

    # Internal TaskEnv Methods

    def discretize_scan_observation(self,data,new_ranges):
        """
        Discards all the laser readings that are not multiple in index of new_ranges
        value.
        """
        self._episode_done = False

        discretized_ranges = []
        mod = len(data.ranges)/new_ranges

        rospy.logdebug("data=" + str(data))
        rospy.logdebug("new_ranges=" + str(new_ranges))
        rospy.logdebug("mod=" + str(mod))

        for i, item in enumerate(data.ranges):
            if (i%mod==0):
                if item == float ('Inf') or numpy.isinf(item):
                    discretized_ranges.append(self.max_laser_value)
                elif numpy.isnan(item):
                    discretized_ranges.append(self.min_laser_value)
                else:
                    discretized_ranges.append(int(item))

                if (self.min_range > item > 0):
                    rospy.logerr("done Validation >>> item=" + str(item)+"< "+str(self.min_range))
                    self._episode_done = True
                else:
                    rospy.logdebug("NOT done Validation >>> item=" + str(item)+"< "+str(self.min_range))

        return discretized_ranges

    def get_vector_magnitude(self, vector):
        """
        It calculated the magnitude of the Vector3 given.
        This is usefull for reading imu accelerations and knowing if there has been
        a crash
        :return:
        """
        contact_force_np = numpy.array((vector.x, vector.y, vector.z))
        force_magnitude = numpy.linalg.norm(contact_force_np)

        return force_magnitude

    def get_orientation_euler(self):
        # We convert from quaternions to euler
        orientation_list = [self.odom.pose.pose.orientation.x,
                            self.odom.pose.pose.orientation.y,
                            self.odom.pose.pose.orientation.z,
                            self.odom.pose.pose.orientation.w]

        roll, pitch, yaw = euler_from_quaternion(orientation_list)
        return roll, pitch, yaw

    def get_distance_from_desired_point(self, current_position):
        """
        Calculates the distance from the current position to the desired point
        :param start_point:
        :return:
        """
        distance = self.get_distance_from_point(current_position,
                                                self.desired_point)

        return distance

    def get_distance_from_point(self, pstart, p_end):
        """
        Given a Vector3 Object, get distance from current position
        :param p_end:
        :return:
        """
        a = numpy.array((pstart.x, pstart.y, pstart.z))
        b = numpy.array((p_end.x, p_end.y, p_end.z))

        distance = numpy.linalg.norm(a - b)

        return distance

    def is_in_desired_position(self,current_position, epsilon=0.05):
        """
        It return True if the current position is similar to the desired poistion
        """

        is_in_desired_pos = False

        x_pos_plus = self.desired_point.x + epsilon
        x_pos_minus = self.desired_point.x - epsilon
        y_pos_plus = self.desired_point.y + epsilon
        y_pos_minus = self.desired_point.y - epsilon

        x_current = current_position.x
        y_current = current_position.y

        x_pos_are_close = (x_current <= x_pos_plus) and (x_current > x_pos_minus)
        y_pos_are_close = (y_current <= y_pos_plus) and (y_current > y_pos_minus)

        is_in_desired_pos = x_pos_are_close and y_pos_are_close

        return is_in_desired_pos

ext: py

reason: import numpy


text:
import numpy as np
import pandas as pd
from python.plot import plot_seir
import python.datamanager as dm
import matplotlib.pyplot as plt
from datetime import datetime, timedelta
# S'(t) = -betaS(t)I(t)/N
# E'(t) = betaS(t)I(t)/N - sigmaE(t)
# I'(t) = sigaE(t) - gammaI(t)
# R'(t) = gammaI(t)
# S(t) --> Susceptibles
# E(t) --> Expuestas
# I(t) --> Infectados
# R(t) --> Recuperados

beta = 1          # Tasa de transmisión, probabilidad de que un susceptible se infecte al entrar en contacto con un infectado.
gamma = 1 / 5     # tasa de recuperación, su inversa es el tiempo medio de recuperación.
sigma = 1 / 7     # tasa de incubación, su inversa es el tiempo promedio de incubación.
N = 48000000  # N Población total
It0 = 10           # Número de infectados iniciales
St0 = N - It0    # Número de susceptibles iniciales
Et0 = 0           # Número de expuestos iniciales
Rt0 = 0           # Número de recuperado iniciales
St = St0
It = It0
Et = Et0
Rt = Rt0
dt = 1/1
days_total = (int)(240 / dt)
days_before = (int)(30 / dt)

model_evolution = np.zeros((days_total, 4))

for day in range(days_before):
    dS = -((beta * St * It)/N) * dt
    dE = ((beta * St * It) / N - (sigma * Et)) * dt
    dI = ((sigma * Et) - (gamma * It)) * dt
    dR = ((gamma * It)) * dt
    St += dS
    Et += dE
    It += dI
    Rt += dR
    model_evolution[day] = [St, Et, It, Rt]
    print(St, Et, It, Rt)

beta = 0.5
for day in range(days_before, days_total):
    dS = -((beta * St * It)/N) * dt
    dE = ((beta * St * It) / N - (sigma * Et)) * dt
    dI = ((sigma * Et) - (gamma * It)) * dt
    dR = ((gamma * It)) * dt
    St += dS
    Et += dE
    It += dI
    Rt += dR
    model_evolution[day] = [St, Et, It, Rt]
    print(St, Et, It, Rt)

df = pd.DataFrame(model_evolution, columns=['S(t)', 'E(t)', 'I(t)', 'R(t)'])
plot_seir(data=df, features=['S(t)', 'E(t)', 'I(t)', 'R(t)'], title='')


text:
#!/usr/bin/env python
#
#
#
# $Header: /opt/cvs/python/packages/share1.5/AutoDockTools/Utilities24/prepare_covalent_flexres.py,v 1.10 2014/08/20 21:10:31 rhuey Exp $
#
import os
from PyBabel.atomTypes import AtomHybridization
from mglutil.math.rigidFit import RigidfitBodyAligner
from MolKit import Read
from MolKit.molecule import AtomSet
from MolKit.bondSelector import RotatableBondSelector, AmideBondSelector
from MolKit.bondSelector import GuanidiniumBondSelector, LeafBondSelector
from mglutil.math.rotax import rotax
from math import sqrt

from AutoDockTools.MoleculePreparation import AD4FlexibleReceptorPreparation, AD4LigandPreparation

try:
    import Numeric
except:
    import numpy.oldnumeric as Numeric

def dist(a1, a2):
    c1 = Numeric.array(a1.coords)
    c2 = Numeric.array(a2.coords)
    d = c2-c1
    return sqrt(sum(d*d))

def getBondedAtoms(atom):
    l = AtomSet()
    for b in atom.bonds:
        l.append(b.neighborAtom(atom))
    return l

def detect_clashes(atomset1, atomset2, cutoff=1.2, verbose=False):  #check this v. autodock/mdist.h
    #exclude bonded pairs
    clashes = []
    for at1 in atomset1:
        if verbose: print "processing ", at1.name,
        bonded = getBondedAtoms(at1)
        if verbose: print "len(bonded)=", len(bonded),
        for at2 in atomset2:
            if verbose: print "vs ", at2.name
            if at2 in bonded:
                continue
            d = dist(at1,at2)
            BondOrderSum = at1.bondOrderRadius + at2.bondOrderRadius
            if verbose: print 'distance = % 6.4f vs BondOrderSum*cutoff = % 6.4f' %(d, BondOrderSum*cutoff)
            if d <= BondOrderSum * cutoff:
                if verbose: print "%s - %s : % 6.4f - % 6.4f" %(at1.name, at2.name, d, BondOrderSum)
                clashes.append(((at1, at2),d))
    return clashes

def setAutoFlexFields(res):
    #process residues
    if hasattr(res, 'setup'):
        return
    res.setup = 1
    res.atoms.used = 0
    res.atoms.bonds[0].possibleTors = 0
    res.atoms.bonds[0].activeTors = 0
    backbone_names = ['C','N','O','HN','HN1','HN2', 'HA',
                'H1','H2','H3','HO', 'H']
    #includes CA
    sidechain = res.atoms.get(lambda x: x.name not in backbone_names)
    res.sideChain = sidechain
    bondlist = res.bondlist = sidechain.bonds[0]
    bondlist.leaf = 0
    bondlist.possibleTors = 0
    bondlist.activeTors = 0
    rbs = RotatableBondSelector()
    rotatables = rbs.select(bondlist)
    for b in rotatables:
        b.possibleTors = 1
        b.activeTors = 1
    amides = AmideBondSelector().select(bondlist)
    for b in amides:
        b.activeTors = 0
        b.possibleTors = 1
    guanidiniums = GuanidiniumBondSelector().select(bondlist)
    for b in guanidiniums:
        b.activeTors = 0
        b.possibleTors = 1
    leaves = LeafBondSelector().select(bondlist)
    for b in leaves:
        b.activeTors = 0
        b.possibleTors = 0
    res.torscount = len(bondlist.get(lambda x: x.activeTors==1))
    #this field is not used in AutoDock4
    res.torsdof = res.torscount
    res.torscount = len(bondlist.get(lambda x: x.activeTors==1))
    res.torsdof = res.torscount

    caAtoms = res.atoms.get(lambda x: x.name=='CA')
    #get returns an AtomSet
    if caAtoms:  #this checks for len(caAtoms)
        res.rootlist = caAtoms
    else:
        res.rootlist = AtomSet([res.atoms.get(lambda x: x._uniqIndex == 0)[0]])
        res.sideChain = res.atoms

def setRelativeTorsion(atom1, atom2, angle): #mov_atoms=None):
    #rat1, rat2, ??30degrees??
    mol = atom1.top
    if mov_atoms is None:
        mov_atoms = mol.subTree(atom1, atom2, mol.allAtoms)
    assert len(mov_atoms)
    mov_coords = Numeric.array(mov_atoms.coords)
    lenCoords = len(mov_coords)
    x = Numeric.array(atom1.coords)
    y = Numeric.array(atom2.coords)
    rot = (angle * 3.14159/180.)%(2 * Numeric.pi)
    matrix = rotax(x, y, rot)
    _ones = Numeric.ones(lenCoords, 'f')
    _ones.shape = (lenCoords,1)
    mov_coords = Numeric.concatenate((mov_coords, _ones),1)
    newcoords = Numeric.dot(mov_coords, matrix)
    nc = newcoords[:,:3].astype('f')
    for i in range(lenCoords):
        at = mov_atoms[i]
        at._coords[at.conformation] = nc[i].tolist()

if __name__ == '__main__':
    import sys
    import getopt

    def usage():
        "Print helpful, accurate usage statement to stdout."
        print "Usage: prepare_covalent_flexres.py -r receptor_filename -l ligand_to_superimpose_filename -b atom1_atom2_atom3 -R residue_name"
        print "    Description of command..."
        print "         -r   receptor_filename"
        print "         -l   ligand_filename"
        print "         -b   ligand atom1_atom2_atom3 to superimpose (eg SD_CG_CB)"
        print "         -R   name of receptor residue onto which the bond atoms "
        print "                   in the ligand will be superimposed"
        print "    Optional parameters:"
        print "        [-B]  alternate atom1_atom2_atom3 in receptor residue to superimpose "
        print "            (default uses atom1_atom2_atom3 names specified with -b flag)"
        print "        [-v]  verbose output"
        print "        [-g pdbqt_filename] (rigid output filename)"
        print "        [-x pdbqt_filename] (flexible output filename)"
        print "        [-n pdbqt_filename] (flexible output filename)"
        print "                                              "
        print "  NOTE the order of the atoms specified for the bond is crucial:"
        print "        '-b atom1_atom2_atom3' "
        print "   specifies this superposition:"
        print "       LIGAND-atom1-atom2-atom3"
        print "              atom1-atom2-atom3-RECEPTOR"
        print "  Process removes atom1-atom2-atom3-RECEPTOR bond and adds bond LIGAND-atom1-atom2-atom3-RECEPTOR"
        print "                                   ^                                                    ^"

    # process command arguments
    try:
        opt_list, args = getopt.getopt(sys.argv[1:], 'r:l:b:R:B:g:x:n:vh')
    except getopt.GetoptError, msg:
        print 'prepare_covalent_flexres.py: %s' %msg
        usage()
        sys.exit(2)

    # initialize required parameters
    #-r: receptor
    receptor_filename =  None
    #-l: ligand
    ligand_filename =  None
    #-b: atom1_atom2_atom3
    atom1_atom2_atom3 =  None
    #-R: res_in_receptor
    res_in_receptor =  None
    # optional parameters
    ratom1_ratom2_ratom3 = None
    # optional parameters
    verbose = None
    #-g  : rigid output filename
    rigid_filename = None
    #-x  : flexible output filename
    flexres_filename = None
    #-n :  newflexible output filename
    new_flexres_filename = None
    #'r:l:b:R:B:g:x:n:vh'
    for o, a in opt_list:
        #print "o=", o, " a=", a
        if o in ('-r', '--r'):
            receptor_filename = a
            #CHECK THAT IT IS ALREADY IN PDBQT FORMAT
            ext = os.path.splitext(receptor_filename)[1]
            if ext!=".pdbqt":
                error = receptor_filename + " not in .pdbqt format"
                raise error
            if verbose: print 'set receptor_filename to ', a
        if o in ('-l', '--l'):
            ligand_filename = a
            ext = os.path.splitext(ligand_filename)[1]
            if ext!=".pdbqt":
                error = ligand_filename + " not in .pdbqt format"
                raise error
            if verbose: print 'set ligand_filename to ', a
        if o in ('-b', '--b'):
            atom1_atom2_atom3 = a
            if verbose: print 'set atom1_atom2_atom3 to ', a
        if o in ('-R', '--R'):
            res_in_receptor = a
            if verbose: print 'set res_in_receptor to ', a
        if o in ('-B', '--B'):
            ratom1_ratom2_ratom3 = a
            if verbose: print 'set ratom1_ratom2_ratom3 to ', a
        if o in ('-v', '--v'):
            verbose = True
            if verbose: print 'set verbose to ', True
        if o in ('-g', '--g'):
            rigid_filename = a
            if verbose: print 'set rigid_filename to ', a
        if o in ('-x', '--x'):
            flexres_filename = a
            if verbose: print 'set flexres_filename to ', a
        if o in ('-n', '--n'):
            new_flexres_filename = a
            if verbose: print 'set new_flexres_filename to ', a
        if o in ('-h', '--h'):
            usage()
            sys.exit()

    if not  receptor_filename:
        print 'prepare_covalent_flexres: receptor filename must be specified!\n'
        usage()
        sys.exit()

    if not  ligand_filename:
        print 'prepare_covalent_flexres: ligand_filename to move must be specified!\n'
        usage()
        sys.exit()

    if not  atom1_atom2_atom3:
        print 'prepare_covalent_flexres: atom1_atom2_atom3 to superimpose must be specified!\n'
        usage()
        sys.exit()

    if not  res_in_receptor:
        print 'prepare_covalent_flexres: res_in_receptor onto which the ligand will superimposed must be specified!\n'
        usage()
        sys.exit()

    if not ratom1_ratom2_ratom3:
        ratom1_ratom2_ratom3 = atom1_atom2_atom3

    #initialize the receptor
    extension = os.path.splitext(receptor_filename)[1]
    if extension!=".pdbqt":
        print 'prepare_covalent_flexres: receptor file must be in .pdbqt format\n'
        usage()
        sys.exit()
    rM = Read(receptor_filename)[0]
    rM.buildBondsByDistance()
    if verbose: print 'read ', receptor_filename

    #initialize the flexible-residue-to-be in the receptor A:MET36
    nsc = res_in_receptor.count(':')
    if nsc==0:
        resP=res_in_receptor
        chINDEX = 0
    elif nsc==1:
        chID,resP = res_in_receptor.split(':')
        chINDEX = rM.chains.id.index(chID)
    elif nsc==2:
        molN,chID,resP = res_in_receptor.split(':')
        chINDEX = rM.chains.id.index(chID)
    else:
        print res_in_receptor + " not properly formatted"
        raise "line 274 ERROR"
    resPSet = rM.chains[chINDEX].residues.get(resP)
    if not len(resPSet):
        msg= resP.name + " is not a residue in " + receptor_filename
        raise msg
    if len(resPSet)>1:
        msg= resP.name + " specifies more than one residue:" + resP.full_name()
        raise msg
    resP = resPSet[0]
    if verbose: print "receptor residue is %s which has %d atoms" %(resP.full_name(), len(resP.atoms))

    #initialize the ligand
    lM = Read(ligand_filename)[0]
    if verbose: print 'read ', ligand_filename
    lM.buildBondsByDistance()            # try to assign bond order;set to 1 if exception
    #lM.allAtoms.bonds[0].bondOrder = 1   #@@@@@@@@CHECK THIS!!!!@@@@@@@@@
    ah = AtomHybridization()
    ah.assignHybridization(lM.allAtoms)

    #process the ligand:
    #(1)locate specified atom1, atom2 and atom3  SD-CG-CB
    #NB: it is assumed that atom3 has exactly 2 bonds
    at_names = atom1_atom2_atom3.split('_')
    assert len(at_names)==3
    at1Set = lM.allAtoms.get(at_names[0])
    if not len(at1Set):
        msg= at_names[0] + " is not in " + ligand_filename
        raise msg
    at1 = at1Set[0] #SD, closest to ligand, farthest from receptor
    if verbose: print "first ligand atom to superimpose is %s" %(at1.full_name())
    at2Set = lM.allAtoms.get(at_names[1])
    if not len(at2Set):
        msg =  at_names[1] + " is not in " + ligand_filename
        raise msg
    at2 = at2Set[0] #CG
    if verbose: print "second ligand atom to superimpose is %s" %(at2.full_name())
    at3Set = lM.allAtoms.get(at_names[2])
    if not len(at3Set):
        msg =  at_names[2] + " is not in " + ligand_filename
        raise msg
    at3 = at3Set[0] #CB, farthest from ligand, closest to receptor
    if verbose: print "third ligand atom to superimpose is %s" %(at3.full_name())

    #(2)setup shortcut to these two atoms !IN ORDER!
    lig_atoms = AtomSet([at1,at2,at3])

    #process the receptor:
    #(1) locate specified receptor residue
    rN,chainID,resN = res_in_receptor.split(':')
    #resSet = rM.chains.residues.get(res_in_receptor.split(':')[-1])
    chS = rM.chains.get(chainID)
    chain = chS[0]
    resSet = chain.residues.get(resN)
    resP = resSet[0]
    if not len(resSet):
        msg = res_in_receptor + " is not in " + receptor_filename
        raise msg
    if len(resSet)!=1:
        msg = res_in_receptor + " selected more than one residue:" + resSet.name
        raise msg
    # SD-CG-CB
    #(2) locate specified ratom1 in that residue
    rat_names = ratom1_ratom2_ratom3.split('_')
    assert len(rat_names)==3
    rat1Set = resP.atoms.get(rat_names[0])
    if not len(rat1Set):
        msg= rat_names[0] + " is not in " + resP.full_name()
        raise msg
    rat1 = rat1Set[0] #SD
    if verbose: print "first receptor atom for superposition is  %s" %(rat1.full_name())
    #(3) locate specified ratom2 in that residue
    rat2Set = resP.atoms.get(rat_names[1])
    if not len(rat2Set):
        msg =  rat_names[1] + " is not in " + resP.full_name()
        raise msg
    rat2 = rat2Set[0] #CG
    if verbose: print "second receptor atom for superposition is  %s" %(rat2.full_name())
    #(4) locate specified ratom3 in that residue
    rat3Set = resP.atoms.get(rat_names[2])
    if not len(rat3Set):
        msg =  rat_names[2] + " is not in " + resP.full_name()
        raise msg
    rat3 = rat3Set[0] #CB, closest to receptor
    if verbose: print "third receptor atom for superposition is  %s" %(rat3.full_name())
    #(5)setup shortcut to these two atoms !IN ORDER!
    ref_atoms = AtomSet([rat1,rat2,rat3])
    if verbose:
        print "for superposition: compute matrix to fit: "
        print "lig_atoms =[at1 = %s, at2 = %s, at3 = %s]" %(at1.full_name(), at2.full_name(), at3.full_name())
        print "onto "
        print "ref_atoms =[rat1 = %s, rat2 = %s, rat3 = %s]" %(rat1.full_name(), rat2.full_name(), rat3.full_name())
    if verbose: print "all input is valid!"
    #(6)locate 'bonds to break' which is the bond from rat3 not to rat2 , here CB: -SD--CG--CB==OtherAtom
    bond_to_break = None
    for b in rat3.bonds:
        otherAt = b.neighborAtom(rat3)
        if otherAt!=rat2:
            bond_to_break = b
            atom_to_attach = otherAt
            if verbose: print "found bond_to_break in receptor => %s" %(b.atom1.name +':'+ b.atom2.name)
            # break the bond
            if verbose: print "receptor atom for ligand attachment in residue %s is %s" %(resP.full_name(),otherAt.full_name())
            otherAt.bonds.remove(b)
            rat3.bonds.remove(b)
            if verbose: print "removed %s-%s bond from %s and from %s"%(b.atom1.name, b.atom2.name, otherAt.full_name(), rat3.full_name())
    if bond_to_break is None:
        print "unable to locate the bond to break in receptor residue %s: expected a bond from %s to an atom other than %s"%(resP.full_name(), rat3.name, rat2.name)
        raise 'invalid input'

    #setup for the superposition
    ref_coords = ref_atoms.coords
    lig_coords = lig_atoms.coords
    #(1) create the tool
    rigidfitAligner = RigidfitBodyAligner()
    rigidfitAligner.setRefCoords(ref_coords)
    rigidfitAligner.rigidFit(lig_coords)
    #(2) get new coordinates for the moving atoms
    #    here all the ligand atoms move
    new_coords = rigidfitAligner.transformCoords(lM.allAtoms.coords)
    # get the correct shape for the coords (natoms,3)
    #new_coords.shape
    #(25, 4)
    new_coords=new_coords[:,:3]
    #(3) add the new to ligand
    lM.allAtoms.addConformation(new_coords)
    confNum = len(lM.allAtoms[0]._coords)-1
    #(4) set the coordinates to use to the new ones
    lM.allAtoms.setConformation(confNum)
    if verbose: print "new coordinates added to ligand"

    #add the ligand atoms other than at1 + at2 to the residue
    #(1) get the atoms to add
    atoms_to_adopt = lM.allAtoms[:]
    #atoms_to_adopt = lM.allAtoms.get(lambda x: not x in lig_atoms)
    if verbose: print "there are %d atoms_to_adopt" %(len(atoms_to_adopt))
    #(3) remove the bonds to the three excluded atoms
    ###attachment_atom.bonds.remove(bond_to_break)
    ###print "removed %s-%s bond from %s"%(bond_to_break.atom1.name, bond_to_break.atom2.name, attachment_atom.full_name())
    ###bond_to_break.neighborAtom(attachment_atom).bonds.remove(bond_to_break)
    if verbose: print 'before removing rat1,rat2 and rat3, residue %s has %d atoms.' %( resP.name, len(resP.atoms))
    #remove rat1, rat2 and rat3
    resP.atoms.remove(rat1)
    if verbose: print 'removed rat1 %s'%rat1.full_name()
    resP.atoms.remove(rat2)
    if verbose: print 'removed rat2 %s'%rat2.full_name()
    resP.atoms.remove(rat3)
    if verbose: print 'removed rat3 %s'%rat3.full_name()
    rM.allAtoms = rM.chains.residues.atoms
    if verbose: print 'after removing rat1,rat2 and rat3 but before adding ligand atoms, residue %s has %d atoms.' %( resP.name, len(resP.atoms))
    if verbose: print 'and the receptor has %d atoms'%len(rM.allAtoms)
    #keep handle to existing atoms for possible clashes
    res_atoms_to_check = resP.atoms[:]
    #(4) add the ligand atoms to the residue
    #for h in atoms_to_adopt:
    for h in lM.allAtoms:
        resP.adopt(h)
        if verbose: print "%s added to %s"%(h.full_name(), resP.full_name())
    if verbose: print 'after adding %d ligand atoms, residue %s has %d atoms.' %( len(lM.allAtoms), resP.name, len(resP.atoms))
    rM.allAtoms = rM.chains.residues.atoms
    if verbose: print 'and the receptor has %d atoms'%len(rM.allAtoms)
    #(5) add the bonds between at3 'CB' in the ligand and the attachment atom in the receptor
    rM.addBond(at3, atom_to_attach)
    if verbose:
        print 'built bond between %s from ligand and  atom %s  in residue %s' %( at3.name,atom_to_attach.name, resP.name)
        print "now the bonds of atom %s are:"%(at3.name)
        for b in at3.bonds:
            print " %s -%s"%(b.atom1.full_name(), b.atom2.full_name())

    #(6) do some necessary updates of the protein data structure
    rM.allAtoms = rM.chains.residues.atoms
    num_ats = len(rM.allAtoms)
    rM.allAtoms.number = range(1, num_ats+1)
    # LOOK FOR CLOSE CONTACTS BETWEEN atoms_to_adopt and res_atoms_to_check
    #if find close contacts, iterate rotation around bond
    #in steps of 30degrees using setAngle stuff...
    clashes = detect_clashes(res_atoms_to_check, atoms_to_adopt, verbose=verbose)
    print "@@@ detected %d clashes! @@@"%len(clashes)
    for cl in clashes:
        print "%s - %s = %6.4f"%(cl[0][0].name, cl[0][1].name, cl[1])
    if len(clashes)>0:
        print 'detected %d clash(s)'%(len(clashes))

    # Setup to use AutoDockTools code to write formatted flexres file
    #(1) detect rotatable bonds etc
    setAutoFlexFields(resP)
    if verbose: print "set autodock fields in atoms in %s" %resP.name
    #to write output
    #(2) setup default output filename if necessary
    #-g  : rigid output filename
    if rigid_filename==None:
        rigid_filename = rM.name + '_rigid.pdbqt'
    if verbose: print "rigid filename is %s" %rigid_filename
    #-x  : flexible output filename
    if flexres_filename==None:
        flexres_filename = rM.name + '_flex.pdbqt'
    if verbose: print "flexible residue filename is %s" %flexres_filename
    #(3) use this class in AutoDockTools for writing formatted outputfiles...
    frp = AD4FlexibleReceptorPreparation(rM, mode='interactive',
                rigid_filename=rigid_filename,
                residues=[resP], flexres_filename=flexres_filename)

    frp.write_flex([resP], flexres_filename)
    frp.write_rigid(rM, rigid_filename)
    if new_flexres_filename is None:
        new_flexres_filename = "new_" + flexres_filename
    flexresMol = Read(flexres_filename)[0]
    LPO = AD4LigandPreparation(flexresMol, root=0, outputfilename = new_flexres_filename)

# To execute this command type:
# prepare_covalent_flexres.py -l filename -r receptor_filename -l ligand_to_superimpose_filename -b atom1_atom2 -R residue
#specific example:
#python -i prepare_covalent_flexres.py -r 1pwc-protein.pdbqt -l 1pwc-ligand.pdbqt -b OG_CB -R SER62  -g 1pwc_rigid.pdbqt -x 1pwc_flex.pdbqt

max_stars_repo_name: bio-hpc/metascreener

max_issues_repo_name: bio-hpc/metascreener

max_forks_repo_name: bio-hpc/metascreener


text:
import time
import numpy as np
import ibmiotf.gateway
from sense_hat import SenseHat
import math

R = [255, 0, 0] # Red
W = [255, 255, 255]  # White
G = [0, 255, 0]  # Green
B = [0, 0, 0] # Black

maskLeft = [[],[10],[10,18],[10,18,26],[10,18,26,34],[10,18,26,34,42],[10,18,26,34,42,50],[10,18,26,34,42,50,9],[10,18,26,34,42,50,9,17],[10,18,26,34,42,50,9,17,25]]
maskRight = [[],[14],[14,22],[14,22,30],[14,22,30,38],[14,22,30,38,46],[14,22,30,38,46,54],[14,22,30,38,46,54,13],[14,22,30,38,46,54,13,21],[14,22,30,38,46,54,13,21,29]]

sense = SenseHat()
def getPixelMap(color, number): #not working for negative numbers
    redPlane = np.tile(R,64)
    redPlane.shape = (64,3)

    greenPlane = np.tile(G,64)
    greenPlane.shape = (64,3)

    blackPlane = np.tile(B, 64)
    blackPlane.shape = (64, 3)

    if color == "red":
        plane = redPlane
    elif color == "green":
        plane = greenPlane
    elif color == "black":
        plane = blackPlane

    number = math.floor(number)
    if number >= 10:
        leftDigit = int(str(number)[0])
        rightDigit = int(str(number)[1])

        plane[maskLeft[leftDigit]] = W
        plane[maskRight[rightDigit]] = W

    elif number <=9:
        plane[maskRight[number]] = W

    plane.shape = (64,3)
    return plane.tolist()

client = None

# Keep track of current led bg color
currBgColor = "black"

def myCommandCallback(cmd):
    if cmd.command == "display":
        command = cmd.data['screen']
        if command == "on":
            # TODO insert your code here
	    # print('command on')
            global currBgColor
	    currBgColor="red"
	    sense.set_pixels(getPixelMap("red", sense.get_temperature()))
            pass
        elif command == "off":
            # TODO insert your code herei
	    # print('command off')
            global currBgColor
	    currBgColor="green"
	    sense.set_pixels(getPixelMap("green", sense.get_temperature()))
            pass

try:
    gatewayOptions = {"org": "xnf2sg", "type": "raspiGateway", "id": "rpi123", "auth-method": "token", "auth-token": "password"}
    gatewayCli = ibmiotf.gateway.Client(gatewayOptions)

    gatewayCli.connect()
    gatewayCli.deviceCommandCallback = myCommandCallback
    gatewayCli.subscribeToDeviceCommands(deviceType='raspiGateway', deviceId='rpi123', command='display',format='json',qos=2)

    while True:
        temp = sense.get_temperature() # TODO insert your code here
        myData = {"d": { "temperature" : temp} } # TODO insert your code here
	sense.set_pixels(getPixelMap(currBgColor, sense.get_temperature()))
        gatewayCli.publishDeviceEvent("senseHAT", "senrpi123", "Event", "json", myData, qos=1 )
        time.sleep(2)

except ibmiotf.ConnectionException  as e:
    print(e)


text:
from typing import Iterable, Union

import numpy as np

import tensorflow as tf

from .base import ModelBase

__all__ = ["MCBN"]

class MCBN(ModelBase):
    """
    Monte Carlo Batch Normalization

    References
    ----------
    M. Teye et al.,
    "Bayesian Uncertainty Estimation for Batch Normalized Deep Networks",
    arXiv 1802.06455, 2018
    """
    BN_class = tf.keras.layers.BatchNormalization

    def __init__(self,units: Iterable[int],
                 noise_variance: float,*,
                 input_shape: Iterable[int]=(1,),
                 dtype: Union[tf.dtypes.DType,np.dtype,str]=tf.float32):
        """
        Initialize MCBN

        Parameters
        ----------
        units : Iterable[int]
            Number of units at hidden layers and output layer.
        noise_variance : float
            Variance of observation noise. (Hyper parameter)

        Raises
        ------
        ValueError
            If `noise_variance` is negative value
        """
        super().__init__(dtype, input_shape,units[-1])

        S = tf.keras.Sequential
        D = tf.keras.layers.Dense
        B = tf.keras.layers.BatchNormalization
        A = tf.keras.layers.Activation

        # Input Layer
        layers = [D(units[0], dtype=self.dtype,
                    kernel_regularizer='l2', use_bias=False,
                    input_shape = self.input_shape),
                  B(),
                  A("relu")]

        # Hidden Layers
        for u in units[1:-1]:
            layers.extend([D(u, dtype=self.dtype,
                             kernel_regularizer='l2', use_bias=False),
                           B(),
                           A("relu")])

        # Output Layer
        layers.append(D(units[-1], dtype=self.dtype))

        self.network = S(layers)

        self.train_data = None
        self.batch_size = None

        if noise_variance < 0.0:
            raise ValueError("`noise_variance` must be positive.")
        self.noise_variance = noise_variance * tf.eye(self.network.layers[-1].units)

    def compile(self,*args,**kwargs):
        """
        Compile network

        Any arguments to be passed to tf.keras.Model.compile
        """
        return self.network.compile(*args,**kwargs)

    def fit(self,x,y,batch_size: int=32,*args,**kwargs):
        """
        Fit network

        Parameters
        ----------
        x : array-like
            Network input values
        y : array-like
            Network output values
        batch_size : int, optional
            Mini batch size. The default size is 32.
        *args
            Arguments to be passed to tf.keras.Model.fit
        **kwargs
            Keyword arguments to be passed to tf.keras.Model.fit

        Raises
        ------
        ValueError
            If batch_size is different from that of the previous fit.
        """
        x = self._ensure_input(x)
        y = self._ensure_output(y)
        if self.train_data is None:
            self.train_data = tf.data.Dataset.from_tensor_slices(x)
        else:
            self.train_data.concatenate(tf.data.Dataset.from_tensor_slices(x))

        if self.batch_size and self.batch_size != batch_size:
            raise ValueError("Batch Size is inconsistent with the previous training")

        self.batch_size = batch_size
        kwargs["batch_size"] = batch_size
        return self.network.fit(x,y,*args,**kwargs)

    def predict(self, x, n_batches: int=100):
        """
        Predict Mean and Covariance

        Parameters
        ----------
        x : array-like
            Input Values
        n_batches : int, optional
            Number of batches to calculate mean and variance. The default is 100

        Returns
        -------
        m : tf.Tensor
            Mean of prediction. [batch size, output units]
        cov : tf.Tensor
            Covariance of prediction. [batch size, output units, output units]
        """
        x = self._ensure_input(x)
        batch_size = tf.constant(self.batch_size,dtype=tf.int64)
        n_batches = tf.constant(n_batches,dtype=tf.int64)

        return self._predict(x,batch_size,n_batches)

    @tf.function
    def _predict(self, x: tf.Tensor, batch_size: tf.Tensor, n_batches: tf.Tensor):
        sum  = 0.0
        sum2 = 0.0

        for B in self._mini_bathes(batch_size, n_batches):
            x_ = tf.identity(x)
            for L in self.network.layers:
                if not isinstance(L, self.BN_class):
                    x_ = L(x_,training=False)
                    B  = L(B,training=False)
                else:
                    mL = tf.reduce_mean(B, axis=0, keepdims=True)
                    vL = (tf.reduce_mean(tf.square(B), axis=0, keepdims=True)
                          - tf.square(mL))

                    norm = tf.sqrt(1.0/(vL + 1e-12))

                    x_ = (x_ - mL) * norm
                    B  = (B  - mL) * norm

                    if L.scale:
                        x_ = L.gamma * x_
                        B  = L.gamma * B

                    if L.center:
                        x_ += L.beta
                        B  += L.beta

            sum  = sum  + x_
            sum2 = sum2 + tf.expand_dims(x_,1) * tf.expand_dims(x_,2)

        inv_batch = 1.0 / tf.cast(n_batches,sum.dtype)
        m = sum  * inv_batch
        v = sum2 * inv_batch - tf.expand_dims(m,1) * tf.expand_dims(m,2)

        return m, v + tf.expand_dims(self.noise_variance,0)

    @tf.function
    def _mini_bathes(self, batch_size: tf.Tensor, n_batches: tf.Tensor):
        buffer_size = self.train_data.cardinality()

        if buffer_size < 0: # -1: INFINITE_CARDINALITY, -2: UNKNOWN_CARDINALITY
            buffer_size = tf.constant(1024,dtype=buffer_size.dtype)

        _dataset = self.train_data.shuffle(buffer_size, reshuffle_each_iteration=True)
        _dataset = _dataset.repeat().batch(batch_size, drop_remainder=True)

        return _dataset.take(n_batches)


text:
#!/usr/bin/env python3

import pandas as pd
import numpy as np
import scipy
from sklearn.cluster import DBSCAN
from sklearn.metrics import accuracy_score

def find_permutation(real_labels, labels):
    permutation = []
    m = max(labels) + 1
    for i in range(m):
        # Initialize np-arrays with no negative labels and real labels
        no_negatives = np.array(labels)
        reals = np.array(real_labels)

        # Get the indices of the labels that are not outliers
        indices = no_negatives != -1

        # Get rid off outliers
        reals = reals[indices]
        no_negatives = no_negatives[indices]

        # Get the true-false values for the labels that match the current label i
        idx = no_negatives == i
        # idx = labels == 1

        # Choose the most common label among data points in the cluster
        new_label = scipy.stats.mode(reals[idx])[0][0]
        # new_label = scipy.stats.mode(real_labels[idx])[0][0]

        # Append the label to the permutation array
        permutation.append(new_label)

    return permutation

def nonconvex_clusters():
    # Dataframe
    df = pd.read_csv("src/data.tsv", sep="\t")
    X = df[["X1", "X2"]]
    y = df.y
    print(X, X.values)

    # Dataframe arrays
    eps = []
    acc_scores = []
    no_clusters = []
    no_outliers = []

    # Loop over eps values and collect results
    for ep in np.arange(0.05, 0.2, 0.05):
        # Initialize model and predict the labels
        model = DBSCAN(eps=ep)
        labels = model.fit_predict(X)

        # Number of labels
        clusters = max(labels) + 1

        # Find the permutation vector and get the new labels
        permutation = find_permutation(y, labels)

        # Get the indices where labels == -1
        idx = labels == -1

        # New labels
        new_labels = [permutation[label] for label in labels[~idx]]

        # Get the accuracy score
        acc = (accuracy_score(y[~idx], new_labels)
               if clusters == 2 else np.nan)

        # Outliers
        outliers = np.sum(idx)

        # Append the values to the arrays
        eps.append(ep)
        acc_scores.append(acc)
        no_clusters.append(clusters)
        no_outliers.append(outliers)

    return pd.DataFrame(list(zip(eps, acc_scores, no_clusters, no_outliers)), columns=["eps", "Score", "Clusters", "Outliers"]).astype(float)

def main():
    print(nonconvex_clusters())

if __name__ == "__main__":
    main()

reason: import numpy,import scipy


text:
"""BLPose Train: Dataset

Author: Bo Lin (@linbo0518)
Date: 2021-01-01
"""
import os
import copy
import cv2
import numpy as np
import matplotlib.pyplot as plt
from pycocotools import mask as _mask
from pycocotools.coco import COCO

from .dataset_utils import COCOUtils

__all__ = ["COCOKeypoints", "COCOUtils"]

class COCOKeypoints:
    def __init__(
        self,
        coco_dir,
        year: int = 2017,
        mode: str = "train",
        to_func=COCOUtils.to_blpose,
        from_func=COCOUtils.from_blpose,
    ):
        if not os.path.isdir(coco_dir):
            raise ValueError(f"{coco_dir} is not a valid dir for coco dataset")

        self.image_dir = os.path.join(coco_dir, f"{mode}{year}")
        annotation_file = os.path.join(
            coco_dir, f"annotations{os.sep}person_keypoints_{mode}{year}.json"
        )

        self.coco = COCO(annotation_file=annotation_file)
        self.image_ids = self.coco.getImgIds(
            catIds=self.coco.getCatIds(catNms="person")
        )
        self.to_func = to_func
        self.from_func = from_func

    def __repr__(self):
        self.coco.info()
        return ""

    def __len__(self):
        return len(self.image_ids)

    def __getitem__(self, item):
        image_id = self.image_ids[item]
        image = cv2.imread(os.path.join(self.image_dir, f"{image_id:012}.jpg"))
        anno = self.coco.loadAnns(self.coco.getAnnIds(imgIds=image_id))
        keypoints_all, mask_all, mask_loss = self._parse_annotations(image.shape, anno)
        return image, keypoints_all, mask_all, mask_loss, anno

    def generate_annotations(self, annotations, keypoints_all, mask_all, mask_loss):
        annotations = copy.deepcopy(annotations)
        new_annotations = []
        offset = 0
        # keypoints_all and mask_all
        for idx, (keypoints, mask) in enumerate(zip(keypoints_all, mask_all)):
            anno = annotations[idx + offset]
            while anno["num_keypoints"] == 0 or anno["iscrowd"]:
                offset += 1
                anno = annotations[idx + offset]

            anno["keypoints"] = self.from_func(keypoints)

            rle = _mask.encode(np.asfortranarray(mask))
            anno["segmentation"] = rle
            anno["area"] = _mask.area(rle)
            anno["bbox"] = _mask.toBbox(rle).tolist()
            new_annotations.append(anno)
        # mask_loss
        rle = _mask.encode(np.asfortranarray(1 - mask_loss))
        new_annotations.append(
            {
                "segmentation": rle,
                "num_keypoints": 0,
                "area": _mask.area(rle),
                "iscrowd": 1,
                "keypoints": [0 for _ in range(17 * 3)],
                "image_id": new_annotations[0]["image_id"],
                "bbox": _mask.toBbox(rle),
                "category_id": new_annotations[0]["category_id"],
                "id": new_annotations[0]["id"],
            }
        )
        return new_annotations

    def show_annotations(self, image, annotations):
        plt.imshow(cv2.cvtColor(image, cv2.COLOR_BGR2RGB))
        self.coco.showAnns(annotations)
        plt.axis("off")
        plt.show()

    def _parse_annotations(self, shape, annotations):
        h, w, _ = shape
        keypoints_all = []
        mask_all = []
        mask_loss = np.full((h, w), True)

        for anno in annotations:
            if anno["num_keypoints"] == 0 or anno["iscrowd"] == 1:
                mask_loss &= self.coco.annToMask(anno) == 0
            else:
                keypoints_all.append(self.to_func(anno["keypoints"]))
                mask_body = self.coco.annToMask(anno) != 0
                mask_all.append(mask_body.astype(np.uint8))
        mask_loss = mask_loss.astype(np.uint8)
        return keypoints_all, mask_all, mask_loss


text:
import streamlit as st

import numpy as np
import pandas as pd
import json

import matplotlib.pyplot as plt
import matplotlib.patches as patches
from matplotlib.offsetbox import OffsetImage, AnnotationBbox
from highlight_text import fig_text

from team_optimization import Team_Optimization

@st.cache
def get_data():
    """ Get player names, xpts for first GW and next GW number

    Returns:
        tuple: Series, int, series
    """
    to = Team_Optimization(
        team_id=35868,
        horizon=5,
        noise=False,
        premium=True)

    return to.data.Name, to.start, to.data[[f'{to.start}_Pts']]

def write():
    st.title('FPL - Sensitivity Analysis Model')
    st.header(
        """
        Sensitivity Analysis FPL Optimization.
        """)

    plt.style.use(".streamlit/style.mplstyle")
    player_names, start, xpts = get_data()

    with st.expander('Basics'):

        col1, col2 = st.columns(2)
        with col1:
            horizon = st.slider("Horizon", min_value=1, max_value=min(39-start, 8), value=min(39-start, 5), step=1)
        with col2:
            premium = st.selectbox("Data type", ['Premium', 'Free'], 0)

        col1, col2, col3, col4 = st.columns(4)
        with col1:
            gk_weight = st.slider("GK Weight", min_value=0.01, max_value=1., value=0.03, step=0.02)
        with col2:
            first_bench_weight = st.slider("1st Weight", min_value=0.01, max_value=1., value=0.21, step=0.02)
        with col3:
            second_bench_weight = st.slider("2nd Weight", min_value=0.01, max_value=1., value=0.06, step=0.02)
        with col4:
            third_bench_weight = st.slider("3rd Weight", min_value=0.01, max_value=1., value=0.01, step=0.02)

        col1, col2, col3, col4, col5 = st.columns(5)
        with col1:
            decay = st.slider("Decay rate", min_value=0., max_value=1., value=0.9, step=0.02)
        with col2:
            vicecap_decay = st.slider("Vicecap rate", min_value=0., max_value=1., value=0.1, step=0.02)
        with col3:
            ft_val = st.slider("FT value", min_value=0., max_value=5., value=1.5, step=0.2)
        with col4:
            hit_val = st.slider("Hit value", min_value=2., max_value=8., value=6., step=0.5)
        with col5:
            itb_val = st.slider("ITB value", min_value=0., max_value=1., value=0.008, step=0.02)

    with st.expander('Chip selection'):

        col1, col2, col3, col4 = st.columns(4)
        with col1:
            wc_gw = st.selectbox("Wildcard", [None] + [start + gw for gw in np.arange(horizon)], 0)
        with col2:
            fh_gw = st.selectbox("Freehit", [None] + [start + gw for gw in np.arange(horizon)], 0)
        with col3:
            tc_gw = st.selectbox("Triple Captain", [None] + [start + gw for gw in np.arange(horizon)], 0)
        with col4:
            bb_gw = st.selectbox("Bench Boost", [None] + [start + gw for gw in np.arange(horizon)], 0)

    with st.expander('Parameters', expanded=True):
        col1, col2 = st.columns(2)
        with col1:
            repeats = st.slider("Number of Experiments", min_value=1, max_value=25, value=5)
        with col2:
            iterations = st.slider("Iterations per exp.", min_value=1, max_value=25, value=7)

    if st.button('Run Optimization'):

        with st.spinner("Running Optimization ..."):
            if fh_gw == 0 and (horizon > 1 or iterations > 1):
                st.warning('Should not have more than 1 iteration or longer than 1 gw horizon')

            else:
                to = Team_Optimization(
                    team_id=35868,
                    horizon=horizon,
                    noise=False,
                    premium=True if premium=='Premium' else False)

                my_bar = st.progress(0)
                progress = 0

                for tmp in to.sensitivity_analysis(
                        repeats=repeats,
                        iterations=iterations,
                        parameters={
                            'model_name':'sensitivity_analysis',
                            'freehit_gw':fh_gw-start if fh_gw is not None else -1,
                            'wildcard_gw':wc_gw-start if wc_gw is not None else -1,
                            'bboost_gw':bb_gw-start if bb_gw is not None else -1,
                            'threexc_gw':tc_gw-start if tc_gw is not None else -1,
                            'objective_type':'decay' if decay != 0 else 'linear',
                            'decay_gameweek':decay,
                            'vicecap_decay':vicecap_decay,
                            'decay_bench':[gk_weight, first_bench_weight, second_bench_weight, third_bench_weight],
                            'ft_val':ft_val,
                            'itb_val':itb_val,
                            'hit_val':hit_val
                        }):

                    progress += 1
                    my_bar.progress(progress / repeats)

                player_names = (
                    pd
                    .read_csv('data/fpl_official/vaastav/data/2021-22/player_idlist.csv')
                    .set_index('id'))

                # Get first GW optimal teams
                with open('optimization/tmp/hashes.json', 'r') as f:
                    hashes = json.load(f)

                df = pd.read_csv("optimization/tmp/podium.csv")

                df['Total'] = df.apply(
                    lambda x: sum(x[str(col)] for col in np.arange(1, iterations+1)),
                    axis=1)

                (
                    df['Transfer in'],
                    df['Transfer out']) = zip(
                        *df['Unnamed: 0'].apply(
                            lambda x: (hashes[str(x)][0], hashes[str(x)][1])))

                max_cols = [col for col in df.columns if 'EV_' in col]
                evs = np.unique(df[max_cols].values)

                df['Mean'] = (
                    df.apply(
                        lambda x: - sum(x[col] if pd.notnull(x[col]) else 0 for col in max_cols) / x["Total"],
                        axis=1))
                df['Std'] = (
                    df.apply(
                        lambda x: np.sqrt(
                            np.sum(
                                np.power(
                                    [x[col] + x['Mean'] if pd.notnull(x[col]) else 0 for col in max_cols],
                                    2)
                                ) / x["Total"]),
                        axis=1))

                df[['Total', '1', '2', '3']] = df[['Total', '1', '2', '3']].astype('int32')
                # df = df.sort_values(['Total', '1', '2', '3'], ascending=False)
                # st.dataframe(df[['Transfer', 'Total', '1', '2', '3', 'Mean', 'Std']])

                if fh_gw == 0 or wc_gw == 0:
                    # Freehit graph

                    # Freehit lineups
                    freehit_teams = pd.DataFrame(
                        zip(*df['Transfer in'].apply(lambda x: write_freehit(x, player_names)))
                        ).T

                    # Gather data on the players in the team
                    percent = pd.DataFrame(columns=['Player', 'Pos', 'Appearences'])
                    percent['Player'] = np.unique(freehit_teams)
                    player_pos = (
                        pd
                        .read_csv('data/fpl_official/vaastav/data/2021-22/cleaned_players.csv')
                        [['first_name', 'second_name', 'element_type']])
                    percent['Pos'] = (
                        percent['Player']
                        .apply(
                            lambda x: player_pos.loc[
                                player_pos.first_name + ' ' + player_pos.second_name == x
                                ]['element_type'].values[0]))
                    percent['Appearences'] = percent['Player'].apply(
                        lambda x: np.sum(np.sum(freehit_teams == x)))

                    percent['Mean'] = percent.apply(
                        lambda x: np.sum(np.sum(freehit_teams == x['Player'], axis=1)*df['Mean'])/x['Appearences'],
                        axis=1)
                    percent['Std'] = percent.apply(
                        lambda x: np.sum(np.sum(freehit_teams == x['Player'], axis=1)*df['Std'])/x['Appearences'],
                        axis=1)

                    with st.expander("Percentages", expanded=False):
                        st.write('Goalkeepers')
                        st.dataframe(
                            percent
                            .loc[percent.Pos == 'GK']
                            .sort_values(
                                by=['Appearences', 'Mean'],
                                ascending=[False, False])
                            [['Player', 'Appearences', 'Mean', 'Std']]
                            .reset_index(drop=True)
                            )

                        st.write('Defenders')
                        st.dataframe(
                            percent
                            .loc[percent.Pos == 'DEF']
                            .sort_values(
                                by=['Appearences', 'Mean'],
                                ascending=[False, False])
                            [['Player', 'Appearences', 'Mean', 'Std']]
                            .reset_index(drop=True)
                            )

                        st.write('Midfielders')
                        st.dataframe(
                            percent
                            .loc[percent.Pos == 'MID']
                            .sort_values(
                                by=['Appearences', 'Mean'],
                                ascending=[False, False])
                            [['Player', 'Appearences', 'Mean', 'Std']]
                            .reset_index(drop=True)
                            )

                        st.write('Forwards')
                        st.dataframe(
                            percent
                            .loc[percent.Pos == 'FWD']
                            .sort_values(
                                by=['Appearences', 'Mean'],
                                ascending=[False, False])
                            [['Player', 'Appearences', 'Mean', 'Std']]
                            .reset_index(drop=True)
                            )

                    # Get team names
                    photos = (
                        pd.read_csv('data/fpl_official/vaastav/data/2021-22/players_raw.csv')
                        [['first_name', 'second_name', 'photo', 'team']])
                    url = "https://resources.premierleague.com/premierleague/photos/players/110x140/p{index}.png"

                    percent['Team'] = percent.apply(
                        lambda x: photos.loc[
                            photos.first_name + ' ' + photos.second_name == x.Player
                            ]['team'].values[0],
                        axis=1)

                    percent = pd.merge(
                        percent,
                        pd.read_csv('data/fpl_official/vaastav/data/2021-22/teams.csv')[['id', 'short_name']],
                        left_on='Team',
                        right_on='id',
                    ).drop(['Team', 'id'], axis=1)

                    # Add upcoming fixtures
                    fixtures = pd.read_csv('data/fpl_official/vaastav/data/2021-22/fixtures.csv')
                    fixtures = fixtures.loc[fixtures.event == start][['team_h', 'team_a']]
                    fixtures = pd.merge(
                        fixtures,
                        pd.read_csv('data/fpl_official/vaastav/data/2021-22/teams.csv')[['id', 'short_name']],
                        left_on='team_h',
                        right_on='id',
                    ).drop(['team_h', 'id'], axis=1).rename(columns={'short_name': 'team_h'})

                    fixtures = pd.merge(
                        fixtures,
                        pd.read_csv('data/fpl_official/vaastav/data/2021-22/teams.csv')[['id', 'short_name']],
                        left_on='team_a',
                        right_on='id',
                    ).drop(['team_a', 'id'], axis=1).rename(columns={'short_name': 'team_a'})

                    def get_fixtures(x, y):
                        if y.team_h == x:
                            return y.team_a

                        if y.team_a == x:
                            return y.team_h.lower()

                    percent['GW'] = percent['short_name'].apply(
                        lambda x: ' + '.join(
                            [
                                game for game in fixtures.apply(
                                    lambda y: get_fixtures(x, y),
                                    axis=1
                                ).values if game is not None]
                            )
                        )

                    # Add xpts data
                    xpts = pd.merge(
                        xpts,
                        player_names,
                        left_index=True,
                        right_index=True,
                    )
                    xpts['Player'] = xpts.apply(lambda x: x.first_name + ' ' + x.second_name, axis=1)

                    percent = pd.merge(
                        percent,
                        xpts,
                        left_on='Player',
                        right_on='Player'
                    )

                    fig, ax = plt.subplots(figsize=(12, 14))

                    # Goalkeeper
                    df = (
                        percent
                        .loc[percent.Pos == 'GK']
                        .sort_values(
                            by=['Appearences', 'Mean'],
                            ascending=[False, False])
                        .reset_index(drop=True)
                        .head(2))

                    for j, row in df.iterrows():
                        rectangle = patches.Rectangle(
                            (25+j*35, 3),
                            17, .25,
                            facecolor='#656B73')
                        ax.add_patch(rectangle)
                        rx, ry = rectangle.get_xy()
                        cx = rx + rectangle.get_width()/2.0
                        cy = ry + rectangle.get_height()/2.0
                        name = row['Player'] if '-' not in row['Player'] else row['Player'].split(' ')[-1]
                        ax.annotate(
                            f"{name}\n{row['GW']}\n{np.round(row['Appearences']/repeats*100)} | xPts: {row[f'{start}_Pts']}",
                            (cx, cy),
                            color='w',
                            weight='bold',
                            fontsize=11,
                            ha='center',
                            va='center')
                        # Plot the portrait
                        imscatter(
                            x=cx, y=3.25+.31,
                            image=url.format(index=(
                                photos.loc[
                                    photos.first_name + ' ' +
                                    photos.second_name == row['Player']
                                ]['photo'].values[0][:-4])),
                            ax=ax, zoom=.4)

                    # Defender
                    df = (
                        percent
                        .loc[percent.Pos == 'DEF']
                        .sort_values(
                            by=['Appearences', 'Mean'],
                            ascending=[False, False])
                        .reset_index(drop=True)
                        .head(5))

                    for j, row in df.iterrows():
                        rectangle = patches.Rectangle(
                            (2.5+j*20, 2),
                            17, .25,
                            facecolor='#656B73')
                        ax.add_patch(rectangle)
                        rx, ry = rectangle.get_xy()
                        cx = rx + rectangle.get_width()/2.0
                        cy = ry + rectangle.get_height()/2.0
                        name = row['Player'] if '-' not in row['Player'] else row['Player'].split(' ')[-1]
                        ax.annotate(
                            # name + '\n' + row['GW'] + '\n' + str(row['Appearences']),
                            f"{name}\n{row['GW']}\n{np.round(row['Appearences']/repeats*100)} | xPts: {row[f'{start}_Pts']}",
                            (cx, cy),
                            color='w',
                            weight='bold',
                            fontsize=11,
                            ha='center',
                            va='center')
                        # Plot the portrait
                        imscatter(
                            x=cx, y=2.25+.31,
                            image=url.format(index=(
                                photos.loc[
                                    photos.first_name + ' ' +
                                    photos.second_name == row['Player']
                                ]['photo'].values[0][:-4])),
                            ax=ax, zoom=.4)

                    # Midfielder
                    df = (
                        percent.loc[percent.Pos == 'MID']
                        .sort_values(
                            by=['Appearences', 'Mean'],
                            ascending=[False, False])
                        .reset_index(drop=True)
                        .head(5))

                    for j, row in df.iterrows():
                        rectangle = patches.Rectangle(
                            (2.5+j*20, 1),
                            17, .25,
                            facecolor='#656B73')
                        ax.add_patch(rectangle)
                        rx, ry = rectangle.get_xy()
                        cx = rx + rectangle.get_width()/2.0
                        cy = ry + rectangle.get_height()/2.0
                        name = row['Player'] if row['Player'] != "Bruno Miguel Borges Fernandes" else "Bruno Fernandes"
                        name = name if '-' not in name else name.split(' ')[-1]
                        ax.annotate(
                            # name + '\n' + row['GW'] + '\n' + str(row['Appearences']),
                            f"{name}\n{row['GW']}\n{np.round(row['Appearences']/repeats*100)} | xPts: {row[f'{start}_Pts']}",
                            (cx, cy),
                            color='w',
                            weight='bold',
                            fontsize=11,
                            ha='center',
                            va='center')
                        # Plot the portrait
                        imscatter(
                            x=cx, y=1.25+.31,
                            image=url.format(index=(
                                photos.loc[
                                    photos.first_name + ' ' +
                                    photos.second_name == row['Player']
                                ]['photo'].values[0][:-4])),
                            ax=ax, zoom=.4)

                    # Forward
                    df = (
                        percent.loc[percent.Pos == 'FWD']
                        .sort_values(
                            by=['Appearences', 'Mean'],
                            ascending=[False, False])
                        .reset_index(drop=True)
                        .head(3))

                    for j, row in df.iterrows():
                        rectangle = patches.Rectangle(
                            (17.5+j*25, 0),
                            17, .25,
                            facecolor='#656B73')
                        ax.add_patch(rectangle)
                        rx, ry = rectangle.get_xy()
                        cx = rx + rectangle.get_width()/2.0
                        cy = ry + rectangle.get_height()/2.0
                        name = row['Player'] if row['Player'] != "Cristiano Ronaldo dos Santos Aveiro" else "Cristiano Ronaldo"
                        name = name if '-' not in name else name.split(' ')[-1]
                        ax.annotate(
                            # name + '\n' + row['GW'] + '\n' + str(row['Appearences']),
                            f"{name}\n{row['GW']}\n{np.round(row['Appearences']/repeats*100)} | xPts: {row[f'{start}_Pts']}",
                            (cx, cy),
                            color='w',
                            weight='bold',
                            fontsize=11,
                            ha='center',
                            va='center')
                        # Plot the portrait
                        imscatter(
                            x=cx, y=.25+.31,
                            image=url.format(index=(
                                photos.loc[
                                    photos.first_name + ' ' +
                                    photos.second_name == row['Player']
                                ]['photo'].values[0][:-4])),
                            ax=ax, zoom=.4)

                    ax.set_ylim(0, 4.25)
                    ax.set_xlim(0, 100)
                    ax.axis('off')

                    fig_text(
                        x=0.14, y=.855,
                        s="<Players appearing most on randomized FH Team>",
                        highlight_textprops=[{"fontweight": "bold"}],
                        fontsize=24, fontfamily="DejaVu Sans", color='w')

                    st.pyplot(fig, ax)
                    plt.close(fig)

                else:
                    df['Transfer'] = (
                        df[['Transfer in', 'Transfer out']]
                        .apply(lambda x: write_transfer(x, player_names), axis=1))

                    fig, ax = plt.subplots(figsize=(16, 12))
                    # Set up the axis limits with a bit of padding
                    ax.set_ylim(0, 5.5)
                    ax.set_xlim(0, 25+(4+1)*16+2.5)
                    ax.axis('off')
                    header_pos = 5.25

                    newaxes = []
                    for j, row in df.head(7)[::-1].reset_index().iterrows():
                        rectangle = patches.Rectangle(
                            (0, j*.75),
                            25, .5,
                            facecolor='#656B73')
                        ax.add_patch(rectangle)
                        rx, ry = rectangle.get_xy()
                        cx = rx + rectangle.get_width()/2.0
                        cy = ry + rectangle.get_height()/2.0
                        ax.annotate(
                            row['Transfer'],
                            (cx, cy),
                            color='w',
                            weight='bold',
                            fontsize=14,
                            ha='center',
                            va='center')

                        newax = fig.add_axes([.63, 0.11+j*.105, .25 , .07])
                        newax.boxplot(
                            [-row[col] for col in max_cols if pd.notnull(row[col]) and row[col] != 0],
                            vert=False,
                            patch_artist=True,
                            capprops=dict(color='#656B73'),
                            boxprops=dict(facecolor='#656B73', color='#656B73'),
                            whiskerprops=dict(color='#656B73'),
                            flierprops=dict(markerfacecolor='r'),
                            medianprops=dict(color='w'),
                            widths=.6)

                        newax.set_xlim(
                            -np.max(evs[evs != 0])-1,
                            -np.min(evs[evs != 0])+1)
                        newaxes.append(newax)

                    # Column headers
                    ax.text(
                        cx, header_pos,
                        'Transfer',
                        fontsize=16, weight='bold', ha='center')

                    # Header separator
                    ax.plot(
                        [0, 25],
                        [5.125, 5.125],
                        ls='-', lw='2.5', c='grey')

                    for i, col in enumerate(['Total', '1', '2', '3']):
                        for j, row in df.head(7)[::-1].reset_index().iterrows():
                            rectangle = patches.Rectangle(
                                (30+i*10, j*.75),
                                5, .5,
                                facecolor='#656B73')
                            ax.add_patch(rectangle)
                            rx, ry = rectangle.get_xy()
                            cx = rx + rectangle.get_width()/2.0
                            cy = ry + rectangle.get_height()/2.0
                            ax.annotate(
                                row[col],
                                (cx, cy),
                                color='w',
                                weight='bold',
                                fontsize=14,
                                ha='center',
                                va='center')

                        # Column headers
                        ax.text(
                            cx, header_pos,
                            col,
                            fontsize=16, weight='bold', ha='center')

                        # Header separator
                        ax.plot(
                            [30+i*10, 30+i*10+5],
                            [5.125, 5.125],
                            ls='-', lw='2.5', c='grey')

                    st.pyplot(fig, ax)
                    plt.close(fig)

def write_transfer(x, player_names):
    """Write into a string the transfers made for a gw

    Args:
        x (array): _description_
        player_names (pd.DataFrame): player names

    Returns:
        str: transfered players
    """
    if len(x['Transfer in']) == len(x['Transfer out']) == 0:
        return 'Roll Transfer'

    else:
        in_list, out_list = [], []
        for transfer_out in x['Transfer out']:
            out_list.append(player_names.loc[transfer_out]['second_name'])

        for transfer_in in x['Transfer in']:
            in_list.append(player_names.loc[transfer_in]['second_name'])

        return " + ".join(out_list) + " -> " + " + ".join(in_list)

def write_freehit(x, player_names):
    """Write into a list the freehit team

    Args:
        x (array): _description_
        player_names (pd.DataFrame): player names

    Returns:
        list: transfered players
    """
    fh_list = []

    for fh_player in x:
        fh_list.append(
            player_names.loc[fh_player]['first_name'] +
            ' ' +
            player_names.loc[fh_player]['second_name'])

    return fh_list

def imscatter(x, y, image, ax=None, zoom=1):
    """stackoverflow.com/questions/35651932/plotting-img-with-matplotlib/35651933"""
    if ax is None:
        ax = plt.gca()

    try:
        image = plt.imread(image)
    except TypeError:
        pass

    im = OffsetImage(image, zoom=zoom)
    x, y = np.atleast_1d(x, y)
    artists = []
    for x0, y0 in zip(x, y):
        ab = AnnotationBbox(im, (x0, y0), xycoords='data', frameon=False)
        artists.append(ax.add_artist(ab))

    ax.update_datalim(np.column_stack([x, y]))
    ax.autoscale()
    return artists

max_stars_repo_name: Fournierp/FPL

max_issues_repo_name: Fournierp/FPL

max_forks_repo_name: Fournierp/FPL


text:
# coding: utf-8

import variational_inference_DP_ngram_mix_shared_base_no_label as vin
import numpy as np
import os,datetime,sys,argparse
import pandas as pd

def get_result_dir():
	with open('result_dir_info.txt') as f:
		dir=f.readlines()[0].strip()
	return dir

def main_loop(
		training_data,
		n,
		sl_T,
		T_base,
		concent_priors,
		max_iters,
		min_increase,
		decoder,
		data_prefix,
		test_data,
		):
# 	data_fileroot=os.path.splitext(datapath.split('/')[-1])[0]
	now = datetime.datetime.now().strftime('%y-%m-%d-%H-%M-%S-%f')
	result_path=os.path.join(
						get_result_dir(),
						data_prefix+('_train-and-test_with_%i-gram_shared_base' % n),
						now
						)
	os.makedirs(result_path)
	vi = vin.VariationalInference(
					sl_T,
					training_data,
					n,
					T_base,
					concent_priors,
					result_path
					)
	vi.train(max_iters, min_increase)

	df_post_pred=pd.DataFrame()
	df_post_pred['log_pred_prob']=vi.get_log_posterior_pred(test_data)
	df_post_pred['test_data_id']=df_post_pred.index
	df_post_pred.to_csv(os.path.join(result_path,'log_pred_probs.csv'), index=False)

	vi.save_results(decoder)

if __name__=='__main__':
	parser = argparse.ArgumentParser()
	parser.add_argument(
		"-n",
		"--ngram",
		default=3,
		type=int,
		help="Context length of ngram"
		)
	parser.add_argument(
		"-d",
		"--data",
		default="../data/BCCWJ_frequencylist_suw_ver1_0_core-nouns_170712180028",
		type=str,
		help="Path to data"
		)
	parser.add_argument(
		"-i",
		"--iterations",
		default=10000,
		type=int,
		help="Maxmum # of iterations"
		)
	parser.add_argument(
		"-T",
		"--threshold",
		default=0.01,
		type=np.float64,
		help="Threshold to detect convergence"
		)
	parser.add_argument(
		"-s",
		"--sublex",
		default=10,
		type=int,
		help="Max # of sublexica"
		)
	options=parser.parse_args()

	datapath = options.data
	data_prefix=datapath.split('/')[-1]
	df_train = pd.read_csv(os.path.join(datapath,data_prefix+'_training_data.tsv'),
						sep='\t',
						encoding='utf-8'
						)
	df_test = pd.read_csv(os.path.join(datapath,data_prefix+'_test_data.tsv'),
						sep='\t',
						encoding='utf-8'
						)
	training_data,test_data,encoder,decoder = vin.code_data(df_train.IPA_csv,test_data=df_test.IPA_csv)

	n = options.ngram
	# T_scalar = 2 # (# of tables) = 2*(# of symbols).
	sl_T=options.sublex
	T_base=len(decoder)*2
	concent_priors = np.array((10.0,10.0)) # Gamma parameters (shape, INVERSE of scale) for prior on concentration.
	max_iters = options.iterations
	min_increase = options.threshold
	start = datetime.datetime.now()
	main_loop(
		training_data,
		n,
		sl_T,
		T_base,
		concent_priors,
		max_iters,
		min_increase,
		decoder,
		data_prefix,
		test_data
		)
	print 'Time spent',str(datetime.datetime.now()-start)


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on May, 2020
@author: Nathan de Lara <ndelara@enst.fr>
"""
from abc import ABC
from typing import Union, Iterable, Tuple

import numpy as np

from sknetwork.utils.base import Algorithm

class BaseLinkPred(Algorithm, ABC):
    """Base class for link prediction algorithms."""

    def _predict_base(self, source: int, targets: Iterable):
        """Prediction for a single node and multiple targets"""
        raise NotImplementedError

    def _predict_node(self, node: int):
        """Prediction for a single node."""
        raise NotImplementedError

    def _predict_nodes(self, nodes: np.ndarray):
        """Prediction for multiple nodes."""
        preds = []
        for node in nodes:
            preds.append(self._predict_node(node))
        return np.array(preds)

    def _predict_edge(self, source: int, target: int):
        """Prediction for a single edge."""
        return self._predict_base(source, [target])[0]

    def _predict_edges(self, edges: np.ndarray):
        """Prediction for a list of edges."""
        preds = []
        for edge in edges:
            i, j = edge[0], edge[1]
            preds.append(self._predict_edge(i, j))
        return np.array(preds)

    def predict(self, query: Union[int, Iterable, Tuple]):
        """Compute similarity scores.

        Parameters
        ----------
        query : int, list, array or Tuple
            * If int i, return the similarities s(i, j) for all j.
            * If list or array integers, return s(i, j) for i in query, for all j as array.
            * If tuple (i, j), return the similarity s(i, j).
            * If list of tuples or array of shape (n_queries, 2), return s(i, j) for (i, j) in query as array.

        Returns
        -------
        predictions : int, float or array
            The prediction scores.
        """
        if np.issubdtype(type(query), np.integer):
            return self._predict_node(query)
        if isinstance(query, Tuple):
            return self._predict_edge(query[0], query[1])
        if isinstance(query, list):
            query = np.array(query)
        if isinstance(query, np.ndarray):
            if query.ndim == 1:
                return self._predict_nodes(query)
            elif query.ndim == 2 and query.shape[1] == 2:
                return self._predict_edges(query)
            else:
                raise ValueError("Query not understood.")
        else:
            raise ValueError("Query not understood.")

    def fit_predict(self, adjacency, query):
        """Fit algorithm to data and compute scores for requested edges."""
        self.fit(adjacency)
        return self.predict(query)

max_stars_repo_name: altana-tech/scikit-network

max_issues_repo_name: altana-tech/scikit-network

max_forks_repo_name: altana-tech/scikit-network


text:
import os
import gensim
import numpy as np
import pytest
from matchms import Spectrum
from spec2vec import SpectrumDocument
from spec2vec.model_building import (set_learning_rate_decay,
                                     train_new_word2vec_model)

def test_set_learning_rate_decay():
    """Test if correct alpha and min_alpha are calculated."""
    alpha, min_alpha = set_learning_rate_decay(0.5, 0.05, 8)
    assert alpha == 0.5, "Expected different alpha."
    assert min_alpha == 0.5 - 8 * 0.05, "Expected different min_alpha"

def test_set_learning_rate_decay_rate_too_high():
    """Test if correct alpha and min_alpha are calculated if rate is too high."""
    alpha, min_alpha = set_learning_rate_decay(0.5, 0.05, 20)
    assert alpha == 0.5, "Expected different alpha."
    assert min_alpha == 0.0, "Expected different min_alpha"

def test_train_new_word2vec_model():
    """Test training of a dummy model."""
    # Create fake corpus
    documents = []
    for i in range(100):
        spectrum = Spectrum(mz=np.linspace(i, 9+i, 10),
                            intensities=np.ones((10)).astype("float"),
                            metadata={})
        documents.append(SpectrumDocument(spectrum, n_decimals=1))
    model = train_new_word2vec_model(documents, iterations=20, vector_size=20,
                                     progress_logger=False)
    assert model.sg == 0, "Expected different default value."
    assert model.negative == 5, "Expected different default value."
    assert model.window == 500, "Expected different default value."
    assert model.alpha == 0.025, "Expected different default value."
    assert model.min_alpha == 0.02, "Expected different default value."
    assert model.epochs == 20, "Expected differnt number of epochs."
    assert model.wv.vector_size == 20, "Expected differnt vector size."
    assert len(model.wv) == 109, "Expected different number of words in vocab."
    assert model.wv.get_vector(documents[0].words[1]).shape[0] == 20, "Expected differnt vector size."

def test_train_new_word2vec_model_with_logger_and_saving(tmp_path):
    """Test training of a dummy model and save it."""
    # Create fake corpus
    documents = []
    for i in range(100):
        spectrum = Spectrum(mz=np.linspace(i, 9+i, 10),
                            intensities=np.ones((10)).astype("float"),
                            metadata={})
        documents.append(SpectrumDocument(spectrum, n_decimals=1))
    # Train model and write to file
    filename = os.path.join(tmp_path, "test.model")
    model = train_new_word2vec_model(documents, iterations=20, filename=filename,
                                     vector_size=20, progress_logger=True)

    # Test if file exists
    assert os.path.isfile(filename), "Could not find saved model file."

    # Test if saved model seems to be correct
    model = gensim.models.Word2Vec.load(filename)
    assert model.sg == 0, "Expected different default value."
    assert model.negative == 5, "Expected different default value."
    assert model.window == 500, "Expected different default value."
    assert model.alpha == 0.025, "Expected different default value."
    assert model.min_alpha == 0.02, "Expected different default value."
    assert model.epochs == 20, "Expected differnt number of epochs."
    assert model.wv.vector_size == 20, "Expected differnt vector size."
    assert len(model.wv) == 109, "Expected different number of words in vocab."
    assert model.wv.get_vector(documents[0].words[1]).shape[0] == 20, "Expected differnt vector size."

def test_train_new_word2vec_model_wrong_entry():
    """Test training of a dummy model with not-accepted gensim argument entry."""
    # Create fake corpus
    documents = []
    for i in range(10):
        spectrum = Spectrum(mz=np.linspace(i, 9+i, 10),
                            intensities=np.ones((10)).astype("float"),
                            metadata={})
        documents.append(SpectrumDocument(spectrum, n_decimals=1))

    with pytest.raises(AssertionError) as msg:
        _ = train_new_word2vec_model(documents, iterations=20, alpha=0.01,
                                     progress_logger=False)

    expected_message_part = "Expect 'learning_rate_initial' instead of 'alpha'."
    assert expected_message_part in str(msg.value), "Expected particular error message."