text:
import cv2
import numpy as np
from math import sqrt
from collections import Counter
from .BrailleCharacter import BrailleCharacter

class SegmentationEngine(object):
    def __init__(self, image = None):
        self.image = image
        self.initialized = False
        self.dots = []
        self.diameter = 0.0
        self.radius = 0.0
        self.next_epoch = 0
        self.characters = []
        return;

    def __iter__(self):
        return self

    def __next__(self):
        return self.next()

    def next(self):
        if not self.initialized:
            self.initialized = True
            contours = self.__process_contours()
            if len(contours) == 0:
                # Since we have no dots.
                self.__clear()
                raise StopIteration()
            enclosingCircles = self.__get_min_enclosing_circles(contours)
            if len(enclosingCircles) == 0:
                self.__clear()
                raise StopIteration()

            diameter,dots,radius = self.__get_valid_dots(enclosingCircles)
            if len(dots) == 0:
                self.__clear()
                raise StopIteration()
            self.diameter = diameter
            self.dots = dots
            self.radius = radius
            self.next_epoch = 0
            self.characters = []

        if len(self.characters) > 0:
            r = self.characters[0]
            del self.characters[0]
            return r

        cor = self.__get_row_cor(self.dots, epoch=self.next_epoch) # do not respect breakpoints
        if cor is None:
            self.__clear()
            raise StopIteration()

        top = int(cor[1] - int(self.radius*1.5)) # y coordinate
        self.next_epoch = int(cor[1] + self.radius)

        cor = self.__get_row_cor(self.dots,self.next_epoch,self.diameter,True)
        if cor is None:
            # Assume next epoch
            self.next_epoch = int(self.next_epoch + (2*self.diameter))
        else:
            self.next_epoch = int(cor[1] + self.radius)

        cor = self.__get_row_cor(self.dots,self.next_epoch,self.diameter,True)
        if cor is None:
            self.next_epoch = int(self.next_epoch + (2*self.diameter))
        else:
            self.next_epoch = int(cor[1] + self.radius)

        bottom = self.next_epoch
        self.next_epoch += int(2*self.diameter)

        DOI = self.__get_dots_from_region(self.dots, top, bottom)
        xnextEpoch = 0
        while True:
            xcor = self.__get_col_cor(DOI, xnextEpoch)
            if xcor is None:
                break

            left = int(xcor[0] - self.radius) # x coordinate
            xnextEpoch = int(xcor[0] + self.radius)
            xcor = self.__get_col_cor(DOI,xnextEpoch,self.diameter,True)
            if xcor is None:
                # Assumed
                xnextEpoch += int(self.diameter*1.5)
            else:
                xnextEpoch = int(xcor[0]) + int(self.radius)
            right = xnextEpoch
            box = (left, right, top, bottom)
            dts = self.__get_dots_from_box(DOI, box)
            char = BrailleCharacter(dts, self.diameter, self.radius, self.image)
            char.left = left
            char.right = right
            char.top = top
            char.bottom = bottom
            self.characters.append(char)

        if len(self.characters) < 1:
            self.__clear()
            raise StopIteration()

        r = self.characters[0]
        del self.characters[0]
        return r

    def __clear(self):
        self.image = None
        self.initialized = False
        self.dots = []
        self.diameter = 0.0
        self.radius = 0.0
        self.next_epoch = 0
        self.characters = []

    def update(self, image):
        self.__clear()
        self.image = image
        return True

    def __get_row_cor(self, dots, epoch = 0, diameter = 0, respectBreakpoint = False):
        if len(dots) == 0:
            return None
        minDot = None
        for dot in dots:
            x,y = dot[0]
            if y < epoch:
                continue

            if minDot is None:
                minDot = dot
            else:
                v = int(y - epoch)
                minV = int(minDot[0][1] - epoch)
                if minV > v:
                    minDot = dot
                else:
                    continue
        if minDot is None:
            return None
        if respectBreakpoint:
            v = int(minDot[0][1] - epoch)
            if v > (2*diameter):
                return None # indicates that the entire row is not set
        return minDot[0] # (X,Y)

    def __get_col_cor(self, dots, epoch = 0, diameter = 0, respectBreakpoint = False):
        if len(dots) == 0:
            return None
        minDot = None
        for dot in dots:
            x,y = dot[0]
            if x < epoch:
                continue

            if minDot is None:
                minDot = dot
            else:
                v = int(x - epoch)
                minV = int(minDot[0][0] - epoch)
                if minV > v:
                    minDot = dot
                else:
                    continue
        if minDot is None:
            return None
        if respectBreakpoint:
            v = int(minDot[0][0] - epoch)
            if v > (2*diameter):
                return None # indicates that the entire row is not set
        return minDot[0] # (X,Y)

    def __get_dots_from_box(self, dots, box):
        left,right,top,bottom = box
        result = []
        for dot in dots:
            x,y = dot[0]
            if x >= left and x <= right and y >= top and y <= bottom:
                result.append(dot)
        return result

    def __get_dots_from_region(self, dots, y1, y2):
        D = []
        if y2 < y1:
            return D

        for dot in dots:
            x,y = dot[0]
            if y > y1 and y < y2:
                D.append(dot)
        return D

    def __get_valid_dots(self, circles):
        tolerance = 0.45
        radii = []
        consider = []
        bin_img = self.image.get_binary_image()
        for circle in circles:
            x,y = circle[0]
            rad = circle[1]
            # OpenCV uses row major
            # Since we do a bitwise not, white pixels belong to the dot.

            # Go through the x axis and check if all those are white
            # pixels till you reach the rad
            it = 0
            while it < int(rad):
                if bin_img[y,x+it] > 0 and bin_img[y+it,x] > 0:
                    it += 1
                else:
                    break
            else:
                if bin_img[y,x] > 0:
                    consider.append(circle)
                    radii.append(rad)

        baserad = Counter(radii).most_common(1)[0][0]
        dots = []
        for circle in consider:
            x,y = circle[0]
            rad = circle[1]
            if rad <= int(baserad * (1+tolerance)) and rad >= int(baserad * (1-tolerance)):
                dots.append(circle)

        # Remove duplicate enclosing circles
        # (i.e) Remove circle enclosed by another other circle.
        for dot in dots:
            X1,Y1 = dot[0]
            C1 = dot[1]
            for sdot in dots:
                if dot == sdot:
                    continue
                X2,Y2 = sdot[0]
                C2 = sdot[1]
                D = sqrt(((X2 - X1)**2) + ((Y2-Y1)**2))
                if C1 > (D + C2):
                    dots.remove(sdot)

        # Filtered base radius
        radii = []
        for dot in dots:
            rad = dot[1]
            radii.append(rad)
        baserad = Counter(radii).most_common(1)[0][0]
        return 2*(baserad), dots, baserad

    def __get_min_enclosing_circles(self, contours):
        circles = []
        radii = []
        for contour in contours:
            (x,y), radius = cv2.minEnclosingCircle(contour)
            center = (int(x), int(y))
            radius = int(radius)
            radii.append(radius)
            circles.append((center, radius))
        return circles

    def __process_contours(self):
        edg_bin_img = self.image.get_edged_binary_image()
        contours = cv2.findContours(edg_bin_img, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
        if len(contours) == 2:
            contours = contours[0]
        else:
            contours = contours[1]
        return contours

max_stars_repo_name: antony-jr/OpticalBrailleRecognition

max_issues_repo_name: antony-jr/OpticalBrailleRecognition

max_forks_repo_name: antony-jr/OpticalBrailleRecognition


text:
# ==============================================================================
# Copyright 2019 - Philip Paquette
#
# NOTICE:  Permission is hereby granted, free of charge, to any person obtaining
#   a copy of this software and associated documentation files (the "Software"),
#   to deal in the Software without restriction, including without limitation the
#   rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
#   sell copies of the Software, and to permit persons to whom the Software is
#   furnished to do so, subject to the following conditions:
#
#   The above copyright notice and this permission notice shall be included in all
#   copies or substantial portions of the Software.
# ==============================================================================
""" Graph Convolution
    - Implements the graph convolution algorithm (ArXiv 1609.02907)
    - Partially adapted from https://github.com/tkipf/gcn/ (MIT License)
"""
import sys
assert 'tensorflow' in sys.modules, 'You need to import TF before importing this module.'

import numpy as np
from diplomacy_research.models.layers.initializers import he, zeros
from diplomacy_research.models.state_space import NB_NODES
from diplomacy_research.utils.tensorflow import tf, batch_norm

# Method described in - SEMI-SUPERVISED CLASSIFICATION WITH GRAPH CONVOLUTIONAL NETWORKS
# ArXiv 1609.02907
def preprocess_adjacency(adjacency_matrix):
    """ Symmetrically normalize the adjacency matrix for graph convolutions.
        :param adjacency_matrix: A NxN adjacency matrix
        :return: A normalized NxN adjacency matrix
    """
    # Computing A^~ = A + I_N
    adj = adjacency_matrix
    adj_tilde = adj + np.eye(adj.shape[0])

    # Calculating the sum of each row
    sum_of_row = np.array(adj_tilde.sum(1))

    # Calculating the D tilde matrix ^ (-1/2)
    d_inv_sqrt = np.power(sum_of_row, -0.5).flatten()
    d_inv_sqrt[np.isinf(d_inv_sqrt)] = 0.
    d_mat_inv_sqrt = np.diag(d_inv_sqrt)

    # Calculating the normalized adjacency matrix
    norm_adj = adj_tilde.dot(d_mat_inv_sqrt).transpose().dot(d_mat_inv_sqrt)
    return np.array(norm_adj, dtype=np.float32)

class GraphConvolution():
    """ Performs a graph convolution (ArXiV 1609.02907) """
    # pylint: disable=too-few-public-methods, too-many-arguments

    def __init__(self, input_dim, output_dim, norm_adjacency, activation_fn=tf.nn.relu, residual=False, bias=False,
                 scope=None, reuse=None):
        """ Initializes the graph convolutional network
            :param input_dim: The number of features per node in the input
            :param output_dim: The number of features per node desired in the output
            :param norm_adjacency: [PLACEHOLDER] The sparse normalized adjacency matrix (NxN matrix)
            :param activation_fn: The activation function to use after the graph convolution
            :param residual: Use residual connection or not.
            :param bias: Boolean flag that indicates we also want to include a bias term
            :param scope: Optional. The scope to use for this layer
            :param reuse: Optional. Boolean. Whether or not the layer and its variables should be reused.
        """
        self.activation_fn = activation_fn if activation_fn is not None else lambda x: x
        self.norm_adjacency = norm_adjacency
        self.bias = bias
        self.var_w, self.var_b = None, None
        self.residual = residual

        # Initializing variables
        with tf.variable_scope(scope, 'GraphConv', reuse=reuse):
            self.var_w = he('W', [NB_NODES, input_dim, output_dim])
            if self.bias:
                self.var_b = zeros('b', [output_dim])

    def __call__(self, inputs):
        """ Actually performs the graph convolution
            :param inputs: The input feature matrix
            :return: The activated output matrix
        """
        # Performs the convolution
        pre_act = tf.transpose(inputs, perm=[1, 0, 2])              # (b, N, in )               => (N, b, in )
        pre_act = tf.matmul(pre_act, self.var_w)                    # (N, b, in) * (N, in, out) => (N, b, out)
        pre_act = tf.transpose(pre_act, perm=[1, 0, 2])             # (N, b, out)               => (b, N, out)
        pre_act = tf.matmul(self.norm_adjacency, pre_act)           # (b, N, N) * (b, N, out)   => (b, N, out)

        # Adds the bias
        if self.bias:
            pre_act += self.var_b                                   # (b, N, out) + (1,1,out) => (b, N, out)

        # Applying activation fn and residual connection
        post_act = self.activation_fn(pre_act)
        if self.residual:
            post_act += inputs
        return post_act                                             # (b, N, out)

def film_gcn_res_block(inputs, gamma, beta, gcn_out_dim, norm_adjacency, is_training, residual=True,
                       activation_fn=tf.nn.relu):
    """ Following the design here https://arxiv.org/pdf/1709.07871.pdf
        :param inputs: A tensor of [b, NB_NODES, gcn_in_dim]
        :param gamma: A tensor of [b, 1, gcn_out_dim]. Used for film
        :param beta: A tensor of [b, 1, gcn_out_dim]. Used for film
        :param gcn_out_dim: number of output channels of graph conv
        :param norm_adjacency: The adjacency matrix for graph conv
        :param is_training: a flag for train/test behavior.
        :param residual: Use residual connection or not.
        :param activation_fn: The activation function on the output. default: relu
        :return: A tensor of [b, NB_NODES, gcn_out_dim]
    """
    gcn_in_dim = inputs.shape.as_list()[-1]
    assert gcn_in_dim == gcn_out_dim or not residual, 'For residual blocks, the in and out dims must be equal'

    gcn_result = GraphConvolution(input_dim=gcn_in_dim,
                                  output_dim=gcn_out_dim,
                                  norm_adjacency=norm_adjacency,
                                  activation_fn=None,
                                  residual=False,
                                  bias=True)(inputs)
    gcn_bn_result = batch_norm(gcn_result, is_training=is_training, fused=True)
    film_result = gamma * gcn_bn_result + beta

    # Applying activation fn and residual connection
    if activation_fn is not None:
        film_result = activation_fn(film_result)
    if residual:
        film_result += inputs
    return film_result

ext: py

reason: import numpy


text:
from flask import Flask,make_response,session
import pandas as pd
import numpy as np
import uuid
import dash
import os

import webbrowser

def get_dynamic_values(id_type):
    res={}
    for inputs_list in dash.callback_context.inputs_list:
        if isinstance(inputs_list,list): #dynamic types
            for cur_input in inputs_list:
                if 'value' in cur_input and isinstance(cur_input['id'],dict) and cur_input['id']['type']==id_type:
                    res[cur_input['id']['index']]=cur_input['value']
    return res

def generate_unique_id():
    return uuid.uuid4().hex[:6].upper()

def open_browser():
    webbrowser.open_new('http://127.0.0.1:5000/cc/')

def load_environ_variables():
    with open('.env') as f:
        for line in f:
            if not line.startswith('#') and '=' in line:
                props=line.split('=')
                if len(props)==2 and not props[0] in os.environ:
                    os.environ[props[0]]=props[1]

max_stars_repo_name: sdswart/pysts

max_issues_repo_name: sdswart/pysts

max_forks_repo_name: sdswart/pysts


text:
import argparse
import sys
from libs.processing import processing
from libs.embedding import embedding
from libs.getModel import get_model
from matplotlib import pyplot as plt
from keract import get_activations
import pandas as pd
import numpy as np

parser = argparse.ArgumentParser(description='keras attention text generator')

parser.add_argument('--file', required=True, help='train data file')
parser.add_argument('--epoch', required=True, default=100, help='epoch')
parser.add_argument('--batch_size', required=False, default=64, help='batch_size')

args = parser.parse_args()

f = open(args.file,'r')
data = processing(f.read())
f.close()

x,y,chars,n_vocab,dataX,seq_length = embedding(data)

n = 300000

m = get_model(x,y)
m.compile(optimizer='adam', loss='categorical_crossentropy')

history  = m.fit(x, y, epochs=int(args.epoch), batch_size=int(args.batch_size), validation_split=0)
import sys
int_to_char = dict((i, c) for i, c in enumerate(chars))

start = np.random.randint(0, len(dataX)-1)
pattern = dataX[start]
print ("Seed:")
print ("\"", ''.join([int_to_char[value] for value in pattern]), "\"")
# generate characters
for i in range(100):
    out_x = np.reshape(pattern, (1, len(pattern), 1))
    out_x = out_x / float(n_vocab)
    prediction = m.predict(out_x, verbose=0)
    index = np.argmax(prediction)
    result = int_to_char[index]
    seq_in = [int_to_char[value] for value in pattern]
    sys.stdout.write(result)
    pattern.append(index)
    pattern = pattern[1:len(pattern)]
print ("\nDone.")

num_simulations = x.shape[2]
attention_vectors = np.zeros(shape=(num_simulations, seq_length))
for i in range(num_simulations):
    #testing_inputs_1, testing_outputs = get_data_recurrent(1, TIME_STEPS, INPUT_DIM)
    activations = get_activations(m, x, layer_name='attention_weight')
    #activations = K.function([m.layers[0].input], [m.layers[1].output])
    attention_vec = np.mean(activations['attention_weight'], axis=0).squeeze()
    assert np.abs(np.sum(attention_vec) - 1.0) < 1e-5
    attention_vectors[i] = attention_vec
print("attention vector: ",attention_vectors)

max_stars_repo_name: ytype/keras-attention-text-generator

max_issues_repo_name: ytype/keras-attention-text-generator

max_forks_repo_name: ytype/keras-attention-text-generator


text:
"""Utilities for graphs"""
from typing import Dict, List, Set, TypeVar

import numpy

Node = TypeVar('Node')
Nodes = Set[Node]
Graph = Dict[Node, Nodes]

def get_transpose_graph(graph: Graph):
    """Get the transpose graph"""
    transpose: Graph = {node: set() for node in graph.keys()}
    for node, target_nodes in graph.items():
        for target_node in target_nodes:
            transpose[target_node].add(node)
    return transpose

def get_indegrees(graph: Graph):
    """Get a dict with indegrees for all nodes"""
    transpose = get_transpose_graph(graph)
    return {node: len(target_nodes) for node, target_nodes in transpose.items()}

def get_induced_subgraph(graph: Graph, nodes: Nodes):
    """Get the nodes-induced subgraph G[S] for a graph G and a subset of nodes S"""
    return {node: graph[node].intersection(nodes) for node in nodes}

def get_dependencies(graph: Graph, node: Node):
    """Get the dependencies of a node"""
    dependencies: Set[Node] = set()
    def traverse_nodes(nodes):
        for candidate in nodes:
            if candidate not in dependencies:
                dependencies.add(candidate)
                traverse_nodes(graph[candidate])
    traverse_nodes(graph[node])
    dependencies.discard(node)
    return dependencies

def get_adjacency_matrix(node_list: List[Node], graph: Graph):
    """Get the adjacency matrix of a graph"""
    node_to_index = {node: index for index, node in enumerate(node_list)}
    adjacency_matrix = numpy.zeros((len(node_list), len(node_list)), dtype=int)
    for node in node_list:
        for dependency in graph[node]:
            if dependency != node:
                adjacency_matrix[
                    node_to_index[node],
                    node_to_index[dependency]
                    ] = 1
    return adjacency_matrix

max_stars_repo_name: andrenarchy/stellar-observatory

max_issues_repo_name: andrenarchy/stellar-observatory

max_forks_repo_name: andrenarchy/stellar-observatory


text:
#!/usr/bin/python

import numpy as np
import nibabel as nib
import cortex
import os

# a script to batch semi-automatic convertion of fsaverage altas Destrieux to pycortex script

def read_label(filename):
    with open(filename, 'r') as f:
        for line in f:
            if line[0] == '#':
                continue
            ss = line.split()
            if len(ss) <= 1:
                continue
            yield int(ss[0])

labels_dir = 'fsaverage'
labels_vertices = dict()
for f in os.listdir(labels_dir):
    labels_vertices[f] = list(read_label(os.path.join(labels_dir, f)))

m = cortex.get_mapper('fsaverage', 'grp_bold')
# print len(labels_vertices)
labels = set(map(lambda x: x.split('.')[1], labels_vertices.keys()))

labels = []
with open('./labels', 'r') as f:
    for line in f:
        labels.append(line.strip())
labels.sort()
print labels

for label in labels:
    lh = np.array(labels_vertices['lh.'+label+'.label'])
    rh = np.array(labels_vertices['rh.'+label+'.label'])
    print label

    vertices = (lh, rh)
    vols = m.backwards(vertices)
    prefix = ('lh', 'rh')
    # print (vols[0] > 0).sum(), lh.shape[0]

    for (hemi, vol) in zip(prefix, vols):
        name = hemi+'.'+label
        print '______', name, vol.sum()
        if vol.sum() > 10:
            volume = cortex.Volume(vol, 'fsaverage', 'grp_bold')
            cortex.add_roi(volume, name=name)
            raw_input()

max_stars_repo_name: kshmelkov/forrestgump

max_issues_repo_name: kshmelkov/forrestgump

max_forks_repo_name: kshmelkov/forrestgump


text:
"""Augmentations pipeline for self-supervised pre-training.

Two variations of augmentations are used. Jittering, where random uniform noise is added to the EEG signal depending on its peak-to-peak values,
along with masking, where signals are masked randomly. Flipping, where the EEG signal is horizontally flipped randomly, and scaling, where EEG
signal is scaled with Gaussian noise.

This file can also be imported as a module and contains the following functions:

    * add_noise - Adds low-frequency and high-frequency noise to the EEG signal.
    * jitter - Applies random uniform noise to the EEG signal.
    * scaling - Applies Gaussian noise to the EEG signal.
    * masking - Masks a single segment of the EEG signal randomly.
    * multi_masking - Masks multiple segments of the EEG signal randomly.
    * flip - Applies a horizontal flip to the EEG signal.
    * augment - Builds the augmentations pipeline.

"""
__author__ = "Likith Reddy, Vamsi Kumar"
__version__ = "1.0.0"
__email__ = "likith012@gmail.com, vamsi81523@gmail.com"

from typing import Optional, Tuple

import numpy as np
import torch
from scipy.interpolate import interp1d

def add_noise(x: torch.Tensor, degree: float) -> torch.Tensor:
    """Adds low-frequency and high-frequency noise to the EEG signal.

    Parameters
    ----------
    x: torch.Tensor
        Input EEG signals.
    degree: float
        Degree of low and high frequency noise to be added.

    Returns
    -------
    torch.Tensor
        Augmented EEG signals.
    """

    len_x = len(x)
    num_range = np.ptp(x) + 1e-4  # add a small number for flat signal

    noise_high_frequency = degree * num_range * (2.0 * np.random.rand(len_x) - 1)
    noise_low_frequency = degree * num_range * (2.0 * np.random.rand(len_x // 100) - 1)
    x_old = np.linspace(0, 1, num=len_x // 100, endpoint=True)
    x_new = np.linspace(0, 1, num=len_x, endpoint=True)
    interpolation = interp1d(x_old, noise_low_frequency, kind="linear")
    noise2 = interpolation(x_new)
    out_x = x + noise_low_frequency + noise_high_frequency
    return out_x

def jitter(x: torch.Tensor, config, degree: float = 1.0) -> torch.Tensor:
    """Applies random uniform noise to the EEG signal.

    Parameters
    ----------
    x: torch.Tensor
        Input EEG signals.
    config
        Configuration object.
    degree: float, optional
        Degree of uniform noise to be added.

    Returns
    -------
    torch.Tensor
        Augmented EEG signals.
    """

    ret = []
    for ch in range(x.shape[0]):
        ret.append(add_noise(x[ch], config.degree * degree))
    ret = np.vstack(ret)
    ret = torch.from_numpy(ret)
    return ret

def scaling(x: torch.Tensor, config, degree: float = 2.0) -> torch.Tensor:
    """Applies Gaussian noise to the EEG signal.

    Parameters
    ----------
     x: torch.Tensor
        Input EEG signals.
    config
        Configuration object.
    degree: float, optional
        Degree of Gaussian noise to be added.

    Returns
    -------
    torch.Tensor
        Augmented EEG signals.
    """

    ret = np.zeros_like(x)
    degree = config.degree * (degree + np.random.rand())
    factor = 2.0 * np.random.normal(size=x.shape[1]) - 1
    factor = 1.5 + (2.0 * np.random.rand()) + degree * factor
    for i in range(x.shape[0]):
        ret[i] = x[i] * factor
    ret = torch.from_numpy(ret)
    return ret

def masking(x: torch.Tensor, config) -> torch.Tensor:
    """Masks a single segment of the EEG signal randomly.

    Parameters
    ----------
    x: torch.Tensor
        Input EEG signals.
    config
        Configuration object.

    Returns
    -------
    torch.Tensor
        Augmented EEG signals.
    """

    segments = config.mask_min_points + int(
        np.random.rand() * (config.mask_max_points - config.mask_min_points)
    )
    points = np.random.randint(0, 3000 - segments)
    ret = x.detach().clone()
    for i, k in enumerate(x):
        ret[i, points : points + segments] = 0

    return ret

def multi_masking(
    x: np.ndarray,
    mask_min: int = 40,
    mask_max: int = 10,
    min_seg: int = 8,
    max_seg: int = 14,
) -> torch.Tensor:
    """Masks multiple segments of the EEG signal randomly.

    Parameters
    ----------
    x: torch.Tensor
        Input EEG signals.
    mask_min: int, optional
        Minimum number of points to be masked in a segment.
    mask_max: int, optional
        Maximum number of points to be masked in a segment.
    min_seg: int, optional
        Minimum number of points segments.
    max_seg: int, optional
        Maximum number of points segments.

    Returns
    -------
    torch.Tensor
        Augmented EEG signals.
    """

    fin_masks = []
    segments = min_seg + int(np.random.rand() * (max_seg - min_seg))
    for seg in range(segments):
        fin_masks.append(mask_min + int(np.random.rand() * (mask_max - mask_min)))
    points = np.random.randint(0, 3000 - segments, size=segments)
    ret = x.clone()
    for i, _ in enumerate(x):
        for seg in range(segments):
            ret[i, points[seg] : points[seg] + fin_masks[seg]] = 0
    return ret

def flip(x: torch.Tensor) -> torch.Tensor:
    """Flips the input EEG Signal.

    Parameters
    ----------
    inputs: torch.Tensor
        Input EEG signal.

    Returns
    -------
    torch.Tensor
        Flipped EEG signal.

    """

    if np.random.rand() > 0.5:
        return torch.tensor(np.flip(x.numpy(), 1).copy())
    else:
        return x

def augment(
    x: torch.Tensor, config, masking_type: Optional[str] = None, degree: float = 3.0
) -> Tuple[torch.Tensor, torch.Tensor]:
    """Applies strong and weak kind of augmentations on the EEG Signal.

    Parameters
    ----------
    inputs: torch.tensor
        Input EEG signal.
    config
        Configuration object.
    masking_type: str, optional
        Type of masking to be applied.
    degree: float, optional
        Degree of noise to be added.

    Returns
    -------
    Tuple[torch.Tensor, torch.Tensor]
       Two different augmented views of the input EEG Signal.

    """

    if masking_type is None:
        weak_augment = masking(jitter(x, config), config)
    else:
        weak_augment = multi_masking(jitter(x, config), config)
    strong_augment = scaling(flip(x, config), config, degree=degree)
    return weak_augment, strong_augment

reason: import numpy,from scipy


text:
#Karan Vombatkere
#Settlers of Catan, 2020

#Imports
from string import *
import numpy as np
from hexTile import *
from hexLib import *
from player import *
#import networkx as nx
#import matplotlib.pyplot as plt
import pygame

pygame.init()

#Class to implement Catan board logic
#Use a graph representation for the board
class catanBoard(hexTile, Vertex):
    'Class Definition for Catan Board Logic'
    #Object Creation - creates a random board configuration with hexTiles
    def __init__(self):
        self.hexTileDict = {} #Dict to store all hextiles, with hexIndex as key
        self.vertex_index_to_pixel_dict = {} #Dict to store the Vertices coordinates with vertex indices as keys
        self.boardGraph = {} #Dict to store the vertex objects with the pixelCoordinates as keys
        self.resourcesList = self.getRandomResourceList()

        self.edgeLength = 80 #Specify for hex size
        self.size = self.width, self.height = 1000, 800
        self.flat = Layout(layout_flat, Point(self.edgeLength, self.edgeLength), Point(self.width/2, self.height/2)) #specify Layout

        #Get a random permutation of indices 0-18 to use with the resource list
        randomIndices = np.random.permutation([i for i in range(len(self.resourcesList))])

        hexIndex_i = 0 #initialize hexIndex at 0
        #Neighbors are specified in adjacency matrix - hard coded

        print("Initializing Game Board...")
        #Generate the hexes and the graphs with the Index, Centers and Resources defined
        for rand_i in randomIndices:
            #Get the coordinates of the new hex, indexed by hexIndex_i
            hexCoords = self.getHexCoords(hexIndex_i)

            #Create the new hexTile with index and append + increment index
            newHexTile = hexTile(hexIndex_i, self.resourcesList[rand_i], hexCoords)
            if(newHexTile.resource.type == 'DESERT'): #Initialize robber on Desert
                newHexTile.robber = True

            self.hexTileDict[hexIndex_i] = newHexTile
            hexIndex_i += 1

        #Create the vertex graph
        self.vertexIndexCount = 0 #initialize vertex index count to 0
        self.generateVertexGraph()

        self.updatePorts() #Add the ports to the graph

        #Initialize DevCardStack
        self.devCardStack = {'KNIGHT':15, 'VP':5, 'MONOPOLY':2, 'ROADBUILDER':2, 'YEAROFPLENTY':2}

        return None

    def getHexCoords(self, hexInd):
        #Dictionary to store Axial Coordinates (q, r) by hexIndex
        coordDict = {0:Axial_Point(0,0), 1:Axial_Point(0,-1), 2:Axial_Point(1,-1), 3:Axial_Point(1,0), 4:Axial_Point(0,1), 5:Axial_Point(-1,1), 6:Axial_Point(-1,0), 7:Axial_Point(0,-2), 8:Axial_Point(1,-2), 9:Axial_Point(2,-2), 10:Axial_Point(2,-1),
                        11:Axial_Point(2,0), 12:Axial_Point(1,1), 13:Axial_Point(0,2), 14:Axial_Point(-1,2), 15:Axial_Point(-2,2), 16:Axial_Point(-2,1), 17:Axial_Point(-2,0), 18:Axial_Point(-1,-1)}
        return coordDict[hexInd]

    #Function to generate a random permutation of resources
    def getRandomResourceList(self):
        #Define Resources as a dict
        Resource_Dict = {'DESERT':1, 'ORE':3, 'BRICK':3, 'WHEAT':4, 'WOOD':4, 'SHEEP':4}
        #Get a random permutation of the numbers
        NumberList = np.random.permutation([2,3,3,4,4,5,5,6,6,8,8,9,9,10,10,11,11,12])
        numIndex = 0

        resourceList = []
        for r in Resource_Dict.keys():
            numberofResource = Resource_Dict[r]
            if(r != 'DESERT'):
                for n in range(numberofResource):
                    resourceList.append(Resource(r, NumberList[numIndex]))
                    numIndex += 1
            else:
                resourceList.append(Resource(r, None))

        return resourceList

    #Function to generate the entire board graph
    def generateVertexGraph(self):
        for hexTile in self.hexTileDict.values():
            hexTileCorners = polygon_corners(self.flat, hexTile.hex) #Get vertices of each hex
            #Create vertex graph with this list of corners
            self.updateVertexGraph(hexTileCorners, hexTile.index)

        #Once all hexTiles have been added  get edges
        self.updateGraphEdges()

    #Function to update a graph of the board with each vertex as a node
    def updateVertexGraph(self, vertexCoordList, hexIndx):
        for v in vertexCoordList:
            #Check if vertex already exists - update adjacentHexList if it does
            if v in self.vertex_index_to_pixel_dict.values():
                for existingVertex in self.boardGraph.keys():
                    if(existingVertex == v):
                        self.boardGraph[v].adjacentHexList.append(hexIndx)

            else:#Create new vertex if it doesn't exist
                #print('Adding Vertex:', v)
                newVertex = Vertex(v, hexIndx, self.vertexIndexCount)
                self.vertex_index_to_pixel_dict[self.vertexIndexCount] = v #Create the index-pixel key value pair
                self.boardGraph[v] = newVertex
                self.vertexIndexCount += 1 #Increment index for future

    #Function to add adges to graph given all vertices
    def updateGraphEdges(self):
        for v1 in self.boardGraph.keys():
            for v2 in self.boardGraph.keys():
                if(self.vertexDistance(v1, v2) == self.edgeLength):
                    self.boardGraph[v1].edgeList.append(v2)

    @staticmethod
    def vertexDistance(v1, v2):
        dist = ((v1.x - v2.x)**2 + (v1.y - v2.y)**2)**0.5
        return round(dist)

    #View the board graph info
    def printGraph(self):
        print(len(self.boardGraph))
        for node in self.boardGraph.keys():
            print("Pixel:{}, Index:{}, NeighborVertexCount:{}, AdjacentHexes:{}".format(node, self.boardGraph[node].vertexIndex, len(self.boardGraph[node].edgeList), self.boardGraph[node].adjacentHexList))

    #Update Board vertices with Port info
    def updatePorts(self):
        #Use this dictionary to map vertex indices to specific ports as per the game board - can add randomization later
        port_dict = {'2:1 BRICK':[43,44], '2:1 SHEEP':[33,34], '2:1 WOOD':[45,49], '2:1 WHEAT':[27,53], '2:1 ORE':[24,29], '3:1 ?':[30,31,36,39,41,42,51,52]}

        #Iterate thru each port and update vertex info
        for portType, portVertexIndex_list in port_dict.items():
            for v_index in portVertexIndex_list: #Each vertex
                vertexPixel = self.vertex_index_to_pixel_dict[v_index] #Get the pixel coordinates to update the boardgraph
                self.boardGraph[vertexPixel].port = portType #Update the port type

    #Function to Display Catan Board Info
    def displayBoardInfo(self):
        for tile in self.hexTileList.values():
            tile.displayHexInfo()
        return None

    #Function to get the list of potential roads a player can build.
    #Return these roads as a dictionary where key=vertex coordinates and values is the rect
    def get_potential_roads(self, player):
        colonisableRoads = {}
        #Check potential roads from each road the player already has
        for existingRoad in player.buildGraph['ROADS']:
            for vertex_i in existingRoad: #Iterate over both vertices of this road
                #Check neighbors from this vertex
                for indx, v_i in enumerate(self.boardGraph[vertex_i].edgeList):
                    if((self.boardGraph[vertex_i].edgeState[indx][1] == False) and (self.boardGraph[vertex_i].state['Player'] in [None, player])): #Edge currently does not have a road and vertex isn't colonised by another player
                        if((v_i, vertex_i) not in colonisableRoads.keys() and (vertex_i, v_i) not in colonisableRoads.keys()): #If the edge isn't already there in both its regular + opposite orientation
                            #Use boolean to keep track of potential roads
                            colonisableRoads[(vertex_i, v_i)] = True
                            #print(vertex_i, v_i)

        return colonisableRoads

    #Function to get available settlements for colonisation for a particular player
    #Return these settlements as a dict of vertices with their Rects
    def get_potential_settlements(self, player):
        colonisableVertices = {}
        #Check starting from each road the player already has
        for existingRoad in player.buildGraph['ROADS']:
            for vertex_i in existingRoad: #Iterate over both vertices of this road
                #Check if vertex isn't already in the potential settlements - to remove double checks
                if(vertex_i not in colonisableVertices.keys()):
                    if(self.boardGraph[vertex_i].isColonised): #Check if this vertex is already colonised
                        break

                    canColonise = True
                    for v_neighbor in self.boardGraph[vertex_i].edgeList: #Check each of the neighbors from this vertex
                        if(self.boardGraph[v_neighbor].isColonised):
                            canColonise = False
                            break

                #If all checks are good add this vertex and its rect as the value
                if(canColonise):
                    #colonisableVertices[vertex_i] = self.draw_possible_settlement(vertex_i, player.color)
                    colonisableVertices[vertex_i] = True

        return colonisableVertices

    #Function to get available cities for colonisation for a particular player
    #Return these cities as a dict of vertex-vertexRect key value pairs
    def get_potential_cities(self, player):
        colonisableVertices = {}
        #Check starting from each settlement the player already has
        for existingSettlement in player.buildGraph['SETTLEMENTS']:
            #colonisableVertices[existingSettlement] = self.draw_possible_city(existingSettlement, player.color)
            colonisableVertices[existingSettlement] = True

        return colonisableVertices

    #Special function to get potential first settlements during setup phase
    def get_setup_settlements(self, player):
        colonisableVertices = {}
        #Check every vertex and every neighbor of that vertex, amd if both are open then we can build a settlement there
        for vertexCoord in self.boardGraph.keys():
            canColonise = True
            potentialVertex = self.boardGraph[vertexCoord]
            if(potentialVertex.isColonised):  #First check if vertex is colonised
                canColonise = False

            #Check each neighbor
            for v_neighbor in potentialVertex.edgeList:
                if(self.boardGraph[v_neighbor].isColonised):  #Check if any of first neighbors are colonised
                    canColonise = False
                    break

            if(canColonise): #If the vertex is colonisable add it to the dict with its Rect
                #colonisableVertices[vertexCoord] = self.draw_possible_settlement(vertexCoord, player.color)
                colonisableVertices[vertexCoord] = True

        return colonisableVertices

    #Special function to get potential first roads during setup phase
    def get_setup_roads(self, player):
        colonisableRoads = {}
        #Can only build roads next to the latest existing player settlement
        latestSettlementCoords = player.buildGraph['SETTLEMENTS'][-1]
        for v_neighbor in self.boardGraph[latestSettlementCoords].edgeList:
            possibleRoad = (latestSettlementCoords, v_neighbor)
            #colonisableRoads[possibleRoad] = self.draw_possible_road(possibleRoad, player.color)
            colonisableRoads[possibleRoad] = True

        return colonisableRoads

    #Function to update boardGraph with Road by player
    def updateBoardGraph_road(self, v_coord1, v_coord2, player):
        #Update edge from first vertex v1
        for indx, v in enumerate(self.boardGraph[v_coord1].edgeList):
            if(v == v_coord2):
                self.boardGraph[v_coord1].edgeState[indx][0] = player
                self.boardGraph[v_coord1].edgeState[indx][1] = True

        #Update edge from second vertex v2
        for indx, v in enumerate(self.boardGraph[v_coord2].edgeList):
            if(v == v_coord1):
                self.boardGraph[v_coord2].edgeState[indx][0] = player
                self.boardGraph[v_coord2].edgeState[indx][1] = True

        #self.draw_road([v_coord1, v_coord2], player.color) #Draw the settlement

    #Function to update boardGraph with settlement on vertex v
    def updateBoardGraph_settlement(self, v_coord, player):
        self.boardGraph[v_coord].state['Player'] = player
        self.boardGraph[v_coord].state['Settlement'] = True
        self.boardGraph[v_coord].isColonised = True

        #self.draw_settlement(v_coord, player.color) #Draw the settlement

    #Function to update boardGraph with settlement on vertex v
    def updateBoardGraph_city(self, v_coord, player):
        self.boardGraph[v_coord].state['Player'] = player
        self.boardGraph[v_coord].state['Settlement'] = False
        self.boardGraph[v_coord].state['City'] = True

        #Remove settlement from player's buildGraph
        player.buildGraph['SETTLEMENTS'].remove(v_coord)

    #Function to update boardGraph with Robber on hexTile
    def updateBoardGraph_robber(self, hexIndex):
        #Set all flags to false
        for hex_tile in self.hexTileDict.values():
            hex_tile.robber = False

        self.hexTileDict[hexIndex].robber = True

    #Function to get possible robber hexTiles
    #Return robber hex spots with their hexIndex - rect representations as key-value pairs
    def get_robber_spots(self):
        robberHexDict = {}
        for indx, hex_tile in self.hexTileDict.items():
            if(hex_tile.robber == False):
                #robberHexDict[indx] = self.draw_possible_robber(hex_tile.pixelCenter)
                robberHexDict[indx] = hex_tile

        return robberHexDict

    #Get a Dict of players to rob based on the hexIndex of the robber, with the circle Rect as the value
    def get_players_to_rob(self, hexIndex):
        #Extract all 6 vertices of this hexTile
        hexTile = self.hexTileDict[hexIndex]
        vertexList = polygon_corners(self.flat, hexTile.hex)

        playersToRobDict = {}

        for vertex in vertexList:
            if(self.boardGraph[vertex].state['Player'] != None): #There is a settlement on this vertex
                playerToRob = self.boardGraph[vertex].state['Player']
                if(playerToRob not in playersToRobDict.keys()): #only add a player once with his/her first settlement/city
                    #playersToRobDict[playerToRob] = self.draw_possible_players_to_rob(vertex)
                    playersToRobDict[playerToRob] = vertex

        return playersToRobDict

    #Function to get a hexTile with a particular number
    def getHexResourceRolled(self, diceRollNum):
        hexesRolled = [] #Empty list to store the hex index rolled (min 1, max 2)
        for hexTile in self.hexTileDict.values():
            if hexTile.resource.num == diceRollNum:
                hexesRolled.append(hexTile.index)

        return hexesRolled

max_stars_repo_name: nazaruka/Catan-AI

max_issues_repo_name: nazaruka/Catan-AI

max_forks_repo_name: nazaruka/Catan-AI

reason: import numpy,import networkx


text:
import numpy as np
# from scipy.sparse import csc_matrix
# from scipy.linalg import qr
# from scipy.sparse.linalg import spsolve
# from numpy.dual import inv

def residual(func,*parameters,x,data):
    '''
    Evaluate the residual value of a fit compare to experimental value x and data.

    Script written by Timothée Chauviré (https://github.com/TChauvire/EPR_ESR_Suite/), 09/09/2020

    Parameters
    ----------
    func : user defined function/model to do the simulation
        TYPE : Python mathematical function definition.
    parameters : Best fit parameters of the user fitted function
        TYPE: Tuple of values
    x : experimental abscissa axis
        TYPE : Numpy data array
    data : experimental data
        TYPE : Numpy data array

    Returns
    -------
    residual : residual spectrum of a fitted spectrum compared to an experimental spectrume
        TYPE : Numpy data array

    '''
    y_model = func(x,*parameters)
    return data-y_model

def ComputeSER(func,*parameters,x,data):
    '''
    Compute the normalized residual sum of square of the residual of a function

    Script written by Timothée Chauviré (https://github.com/TChauvire/EPR_ESR_Suite/), 09/09/2020

    Parameters
    ----------
    func : user defined function/model to do the simulation
        TYPE : Python mathematical function definition.
    parameters : Best fit parameters of the user fitted function
        TYPE: Tuple of values
    x : experimental abscissa axis
        TYPE : Numpy data array
    data : experimental data
        TYPE : Numpy data array

    Returns
    -------
    SER : normalized residual sum of square
        TYPE : real float value

    '''
    y = func(x,*parameters)
    NumMeas=data.shape[0]
    NumParams=len(parameters)
    resnorm = np.sum((y-data)**2)
    SER = np.sqrt(resnorm/(NumMeas - NumParams))
    return SER

def ComputeJacobian(func, *params, x, eps = 1e-12):
    '''
    Compute the Jacobian of a function.

    Script written by Timothée Chauviré (https://github.com/TChauvire/EPR_ESR_Suite/), 09/09/2020

    Parameters
    ----------
    func : user defined function/model to do the simulation
        TYPE : python function definition with parameters *params
        DESCRIPTION.
    *params : Best fit parameters of the user fitted function
        TYPE: Tuple of values
    x : abscissa axis associated to the experimental data to be fitted.
        TYPE : Numpy data array
    eps : precision
        TYPE= real float, optional
        DESCRIPTION. The default is 1e-12.

    Returns
    -------
    J : jacobian matrix
        TYPE : Numpy data array

    '''

    # Not sure of this implementation, result are weird but with good dimension...
    nparam = len(params)
    npts = len(x)
    J = np.full((npts, nparam),np.nan)
    params2 = np.full((nparam,),np.nan)
    paraminit = np.full((nparam,),np.nan)
    for (i, value) in enumerate(params):
        paraminit[i] = value
    for i in range(nparam):
        params2 =  paraminit
        params2[i] = paraminit[i]+eps
        f1 = func(x,*params2)
        params2[i] = paraminit[i]-eps
        f2 = func(x,*params2)
        J[ : , i] = (f1 - f2) / (2 * eps)
    return J

# def ComputeJacobian(func,*parameters, x, dx=1e-8):
#     npts = len(x)
#     nparam = len(parameters)
#     jac = np.full((nparam, nparam),np.nan)
#     for i in range(nparam):  # through columns to allow for vector addition
#         for j in range(npts):
#             Dxj = (abs(x[j])*dx if x[j] != 0 else dx)
#             x_plus = np.array([(xi if k != j else xi + Dxj) for k, xi in enumerate(x)])
#             jac[:, j] = (func(x_plus,*parameters) - func(x,*parameters))/Dxj
#     return jac

def error_jacobian(SER=None,jac=None,*args,**kwargs):
    '''
    Compute the error on parameters associated to non linear fit based on Jacobian matrix.

    Script written by Timothée Chauviré (https://github.com/TChauvire/EPR_ESR_Suite/), 09/09/2020

    Parameters
    ----------
    SER : normalized residual sum of square
        TYPE : Numpy data array
        DESCRIPTION. The default is None.
    jac : jacobian matrix,
        TYPE : Numpy data array
        DESCRIPTION. The default is None.

    Returns
    -------
    error_parameters : error evaluation of the fitted parameters of a fit
        TYPE : Numpy data array

    covariance : covariance matrix
        TYPE : Numpy data array

    '''

    # Problem how to deal with singular matrix...
    covariance= np.dot(SER**2,np.linalg.inv(np.dot(jac.T,jac)))
    error_parameters=np.sqrt(np.diag(covariance))

    # # Other solution :
    # Q,R,P=qr(jac,mode='full',pivoting=True)
    # # R is upper triangular but not in size ==> reshaping is necessary :)
    # Rinv=inv(R[0:R.shape[-1],0:R.shape[-1]])
    # JTJInverse=np.matmul(Rinv,Rinv.T)
    # covariance=np.dot(SER ** 2,JTJInverse)
    return error_parameters,covariance

def cov_to_correlationmatrix(cov=None,*args,**kwargs):
    '''
    Convert covariance matrix to correlation matrix

    Script written by Timothée Chauviré (https://github.com/TChauvire/EPR_ESR_Suite/), 09/09/2020

    Parameters
    ----------
    cov = covariance matrix
        TYPE : Numpy data array
    Returns
    -------
    corr : correlation matrix
        TYPE : Numpy data array
    '''
    v = np.sqrt(np.diag(cov))
    outer_v = np.outer(v, v)
    corr = cov / outer_v
    corr[cov == 0] = 0
    return corr

def stddev_covnmatrix(cov=None,*args,**kwargs):
    '''
    Evaluate standard deviation of fitted parameters from covariance matrix from .

    Script written by Timothée Chauviré (https://github.com/TChauvire/EPR_ESR_Suite/), 09/09/2020

    Parameters
    ----------
    cov : covariance matrix,
        TYPE : Numpy data array
        DESCRIPTION. The default is None.

    Returns
    -------
    std : standard deviation associated to parameters obtained via a fit of experimental value.
        TYPE : Numpy data array
    '''
    std = np.sqrt(np.diag(cov))
    return std

def corr_to_covariancematrix(corr=None,std=None,*args,**kwargs):
    '''Convert correlation matrix to covariance matrix

    Script written by Timothée Chauviré (https://github.com/TChauvire/EPR_ESR_Suite/), 09/09/2020

    Parameters
    ----------
    corr = Correlation matrix
        TYPE : Numpy data array
    std = standard deviation of each variable
        TYPE : Numpy data array

    Returns
    -------
    cov : covariance matrix,
        TYPE : Numpy data array
    '''
    cov = corr*std*std.reshape(-1,1) # Calculate the covariance matrix
    return cov

def error_vandermonde(x,residuals=None,rank=None,*args,**kwargs):
    '''
    Function to generate 1) error estimation on parameters determined by the function
    coef, [residuals, rank, singular_values, rcond] = numpy.polynomial.polynomial.polyfit()
    2) covariance matrix associated.

    The Vandermonde matrix generated by vand = polyvander(x,rank)
    The Covariance matrix cov is obtained via the vandermonde matrix V
    via this numerical steps:
        1) compute np.dot(V.T,V).inv
        2) and multiply it by the residual/(nb of data points - nb of coefficients)
    The error parameters are then computed via : error_parameters = np.sqrt(np.diag(cov))

    Script written by Timothée Chauviré (https://github.com/TChauvire/EPR_ESR_Suite/), 09/09/2020

    Parameters
    ----------
    residuals : first value generated by polyfit in the list of the second output
        DESCRIPTION. float number. The default is None.
    vandermonde : Vandermonde matrix generated by vand = polyvander(x,rank)
        DESCRIPTION. The default is None.

    rank : necessary for multidimensional array
        DESCRIPTION. The default is None.

    Raises
    ------
    ValueError if rank is higher than the number of points
        "the number of data points must exceed order "
                                 "to scale the covariance matrix".

    Returns
    -------
    error_parameters : error uncertainties estimated on the parameters ordered
    from low to high polynomial order.
    By example for a linear model f(x)=a0+a1*x :
        error_parameters[0] = constant parameters a0,
        error_parameters[1] = slope a1.

    cov : covariance estimated for the Least-squares fit of a polynomial to data
        return a Matrix(Rank by rank)
    '''

    from numpy.polynomial.polynomial import polyvander

    if len(x) <= rank:
        raise ValueError("the number of data points must exceed order "
                                 "to scale the covariance matrix")
            # note, this used to be: fac = resids / (len(x) - order - 2.0)
            # it was deciced that the "- 2" (originally justified by "Bayesian
            # uncertainty analysis") is not was the user expects
            # (see gh-11196 and gh-11197)
    else:
        v = polyvander(x,rank) # generate the vandermonde matrix
        cov = residuals/(len(x) - rank)*np.linalg.inv(np.dot(v.T, v))
        error_parameters=np.sqrt(np.diag(cov))

    return error_parameters,cov

# Ideas To Do or test:
#     r,p = scipy.stats.pearsonr(x, y)  # and other optional statistical test
#     reference https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.pearsonr.html#scipy.stats.pearsonr
#     https://en.wikipedia.org/wiki/Pearson_correlation_coefficient

#     cov = np.cov(x,y)

# Code from scipy.optimize.curve_fit: (https://github.com/scipy/scipy/blob/v1.5.1/scipy/optimize/minpack.py#L532-L834)
# Do Moore-Penrose inverse discarding zero singular values.
#         _, s, VT = svd(res.jac, full_matrices=False)
#         threshold = np.finfo(float).eps * max(res.jac.shape) * s[0]
#         s = s[s > threshold]
#         VT = VT[:s.size]
#         pcov = np.dot(VT.T / s**2, VT)

# Check lmfit examples :
#    https://github.com/lmfit/lmfit-py/blob/master/examples/example_fit_with_derivfunc.py

#     bayesian analysis in python

reason: import numpy,from numpy,from scipy


text:
import numpy as np
from math import sin, cos, atan2, tan, sqrt, pi

import matplotlib.pyplot as plt
import time

from bdsim.components import TransferBlock, FunctionBlock
from bdsim.graphics import GraphicsBlock
# ------------------------------------------------------------------------ #

class MultiRotor(TransferBlock):
    """
    :blockname:`MULTIROTOR`

    .. table::
       :align: left

       +------------+---------+---------+
       | inputs     | outputs |  states |
       +------------+---------+---------+
       | 1          | 1       | 16      |
       +------------+---------+---------+
       | A(4,)      | dict    |         |
       +------------+---------+---------+
    """
    nin = 1
    nout = 1

	# Flyer2dynamics lovingly coded by Paul Pounds, first coded 12/4/04
	# A simulation of idealised X-4 Flyer II flight dynamics.
	# version 2.0 2005 modified to be compatible with latest version of Matlab
	# version 3.0 2006 fixed rotation matrix problem
	# version 4.0 4/2/10, fixed rotor flapping rotation matrix bug, mirroring
	# version 5.0 8/8/11, simplified and restructured
	# version 6.0 25/10/13, fixed rotation matrix/inverse wronskian definitions, flapping cross-product bug
	#
	# New in version 2:
	#   - Generalised rotor thrust model
	#   - Rotor flapping model
	#   - Frame aerodynamic drag model
	#   - Frame aerodynamic surfaces model
	#   - Internal motor model
	#   - Much coolage
	#
	# Version 1.3
	#   - Rigid body dynamic model
	#   - Rotor gyroscopic model
	#   - External motor model
	#
	# ARGUMENTS
	#   u       Reference inputs                1x4
	#   tele    Enable telemetry (1 or 0)       1x1
	#   crash   Enable crash detection (1 or 0) 1x1
	#   init    Initial conditions              1x12
	#
	# INPUTS
	#   u = [N S E W]
	#   NSEW motor commands                     1x4
	#
	# CONTINUOUS STATES
	#   z      Position                         3x1   (x,y,z)
	#   v      Velocity                         3x1   (xd,yd,zd)
	#   n      Attitude                         3x1   (Y,P,R)
	#   o      Angular velocity                 3x1   (wx,wy,wz)
	#   w      Rotor angular velocity           4x1
	#
	# Notes: z-axis downward so altitude is -z(3)
	#
	# CONTINUOUS STATE MATRIX MAPPING
	#   x = [z1 z2 z3 n1 n2 n3 z1 z2 z3 o1 o2 o3 w1 w2 w3 w4]
	#
	#
	# CONTINUOUS STATE EQUATIONS
	#   z` = v
	#   v` = g*e3 - (1/m)*T*R*e3
	#   I*o` = -o X I*o + G + torq
	#   R = f(n)
	#   n` = inv(W)*o
	#
    def __init__(self, model, *inputs, groundcheck=True, speedcheck=True, x0=None, **kwargs):
        r"""
        Create a a multi-rotor dynamic model block.

        :param model: Vehicle geometric and inertial parameters
        :type model: dict
        :param ``*inputs``: Optional incoming connections
        :type ``*inputs``: Block or Plug
        :param groundcheck: Prevent vehicle moving below ground, defaults to True
        :type groundcheck: bool
        :param speedcheck: Check for zero rotor speed, defaults to True
        :type speedcheck: bool
        :param x0: Initial state, defaults to None
        :type x0: float, optional
        :param ``**kwargs``: common Block options
        :return: a MULTIROTOR block
        :rtype: MultiRotor instance

        **Block ports**

            :input ω: a vector of input rotor speeds in (radians/sec).  These are,
                looking down, clockwise from the front rotor which lies on the x-axis.

            :output x: a dictionary signal with the following items:

                - ``x`` pose in the world frame as :math:`[x, y, z, \theta_Y, \theta_P, \theta_R]`
                - ``vb`` translational velocity in the world frame (metres/sec)
                - ``w`` angular rates in the world frame as yaw-pitch-roll rates (radians/second)
                - ``a1s`` longitudinal flapping angles (radians)
                - ``b1s`` lateral flapping angles (radians)

        Based on MATLAB code developed by Pauline Pounds 2004.
        """
        super().__init__(nin=1, nout=1, inputs=inputs, **kwargs)
        self.type = 'quadrotor'

        try:
            nrotors = model['nrotors']
        except KeyError:
            raise RuntimeError('vehicle model does not contain nrotors')
        assert nrotors % 2 == 0, 'Must have an even number of rotors'

        self.nstates = 12
        if x0 is not None:
            assert len(x0) == self.nstates, "x0 is the wrong length"
        else:
            x0 = np.zeros((self.nstates,))
        self._x0 = x0

        self.nrotors = nrotors
        self.model = model

        self.groundcheck = groundcheck
        self.speedcheck = speedcheck

        self.D = np.zeros((3,self.nrotors))
        self.theta = np.zeros((self.nrotors,))
        for i in range(0, self.nrotors):
            theta = i / self.nrotors * 2 * pi
            #  Di      Rotor hub displacements (1x3)
            # first rotor is on the x-axis, clockwise order looking down from above
            self.D[:,i] = np.r_[ model['d'] * cos(theta), model['d'] * sin(theta), model['h']]
            self.theta[i] = theta

        self.a1s = np.zeros((self.nrotors,))
        self.b1s = np.zeros((self.nrotors,))

    def output(self, t=None):

        model = self.model

        # compute output vector as a function of state vector
        #   z      Position                         3x1   (x,y,z)
        #   v      Velocity                         3x1   (xd,yd,zd)
        #   n      Attitude                         3x1   (Y,P,R)
        #   o      Angular velocity                 3x1   (Yd,Pd,Rd)

        n = self._x[3:6]   # RPY angles
        phi = n[0]         # yaw
        the = n[1]         # pitch
        psi = n[2]         # roll

        # rotz(phi)*roty(the)*rotx(psi)
        #  BBF > Inertial rotation matrix
        R = np.array([
                [cos(the) * cos(phi), sin(psi) * sin(the) * cos(phi) - cos(psi) * sin(phi), cos(psi) * sin(the) * cos(phi) + sin(psi) * sin(phi)],
                [cos(the) * sin(phi), sin(psi) * sin(the) * sin(phi) + cos(psi) * cos(phi), cos(psi) * sin(the) * sin(phi) - sin(psi) * cos(phi)],
                [-sin(the),           sin(psi) * cos(the),                                  cos(psi) * cos(the)]
            ])

        #inverted Wronskian
        iW = np.array([
                    [0,        sin(psi),             cos(psi)],
                    [0,        cos(psi) * cos(the), -sin(psi) * cos(the)],
                    [cos(the), sin(psi) * sin(the),  cos(psi) * sin(the)]
                ]) / cos(the)

        # return velocity in the body frame
        out = {}
        out['x'] = self._x[0:6]
        out['vb'] = np.linalg.inv(R) @ self._x[6:9]   # translational velocity mapped to body frame
        out['w'] = iW @ self._x[9:12]               # RPY rates mapped to body frame
        out['a1s'] = self.a1s
        out['b1s'] = self.b1s
        out['X'] = self._x

        return [out]

    def deriv(self):

        model = self.model

        # Body-fixed frame references
        #   ei      Body fixed frame references 3x1
        e3 = np.r_[0, 0, 1]

        # process inputs
        w = self.inputs[0]
        if len(w) != self.nrotors:
            raise RuntimeError('input vector wrong size')

        if self.speedcheck and np.any(w == 0):
            # might need to fix this, preculudes aerobatics :(
            # mu becomes NaN due to 0/0
            raise RuntimeError('quadrotor_dynamics: not defined for zero rotor speed');

        # EXTRACT STATES FROM X
        z = self._x[0:3]   # position in {W}
        n = self._x[3:6]   # RPY angles {W}
        v = self._x[6:9]   # velocity in {W}
        o = self._x[9:12]  # angular velocity in {W}

        # PREPROCESS ROTATION AND WRONSKIAN MATRICIES
        phi = n[0]    # yaw
        the = n[1]    # pitch
        psi = n[2]    # roll

        # rotz(phi)*roty(the)*rotx(psi)
        # BBF > Inertial rotation matrix
        R = np.array([
            [cos(the)*cos(phi), sin(psi)*sin(the)*cos(phi)-cos(psi)*sin(phi), cos(psi)*sin(the)*cos(phi)+sin(psi)*sin(phi)],
            [cos(the)*sin(phi), sin(psi)*sin(the)*sin(phi)+cos(psi)*cos(phi), cos(psi)*sin(the)*sin(phi)-sin(psi)*cos(phi)],
            [-sin(the),         sin(psi)*cos(the),                            cos(psi)*cos(the)]
            ])

        # Manual Construction
        #     Q3 = [cos(phi) -sin(phi) 0;sin(phi) cos(phi) 0;0 0 1];   % RZ %Rotation mappings
        #     Q2 = [cos(the) 0 sin(the);0 1 0;-sin(the) 0 cos(the)];   % RY
        #     Q1 = [1 0 0;0 cos(psi) -sin(psi);0 sin(psi) cos(psi)];   % RX
        #     R = Q3*Q2*Q1    %Rotation matrix
        #
        #    RZ * RY * RX
        # inverted Wronskian
        iW = np.array([
                    [0,        sin(psi),          cos(psi)],
                    [0,        cos(psi)*cos(the), -sin(psi)*cos(the)],
                    [cos(the), sin(psi)*sin(the), cos(psi)*sin(the)]
                ]) / cos(the)

        # ROTOR MODEL
        T = np.zeros((3,4))
        Q = np.zeros((3,4))
        tau = np.zeros((3,4))

        a1s = self.a1s
        b1s = self.b1s

        for i in range(0, self.nrotors):  # for each rotor

            # Relative motion
            Vr = np.cross(o, self.D[:,i]) + v
            mu = sqrt(np.sum(Vr[0:2]**2)) / (abs(w[i]) * model['r'])  # Magnitude of mu, planar components
            lc = Vr[2] / (abs(w[i]) * model['r'])                     # Non-dimensionalised normal inflow
            li = mu                                                  # Non-dimensionalised induced velocity approximation
            alphas = atan2(lc, mu)
            j = atan2(Vr[1], Vr[0])                                  # Sideslip azimuth relative to e1 (zero over nose)
            J = np.array([
                    [cos(j), -sin(j)],
                    [sin(j),  cos(j)]
                ])                                                   # BBF > mu sideslip rotation matrix

            # Flapping
            beta = np.array([
                    [((8/3*model['theta0'] + 2 * model['theta1']) * mu - 2 * lc * mu) / (1 - mu**2 / 2)], # Longitudinal flapping
                    [0]                                                              # Lattitudinal flapping (note sign)
                ])

                # sign(w) * (4/3)*((Ct/sigma)*(2*mu*gamma/3/a)/(1+3*e/2/r) + li)/(1+mu^2/2)];

            beta = J.T @ beta;                                    # Rotate the beta flapping angles to longitudinal and lateral coordinates.
            a1s[i] = beta[0] - 16 / model['gamma'] / abs(w[i]) * o[1]
            b1s[i] = beta[1] - 16 / model['gamma'] / abs(w[i]) * o[0]

            # Forces and torques

            # Rotor thrust, linearised angle approximations

            T[:,i] = model['Ct'] * model['rho'] * model['A'] * model['r']**2 * w[i]**2 * \
                np.r_[-cos(b1s[i]) * sin(a1s[i]), sin(b1s[i]), -cos(a1s[i])*cos(b1s[i])]

            # Rotor drag torque - note that this preserves w[i] direction sign

            Q[:,i] = -model['Cq'] * model['rho'] * model['A'] * model['r']**3 * w[i] * abs(w[i])* e3

            tau[:,i] = np.cross(T[:,i], self.D[:,i])    # Torque due to rotor thrust

        # RIGID BODY DYNAMIC MODEL
        dz = v
        dn = iW @ o

        dv = model['g'] * e3 + R @ np.sum(T, axis=1) / model['M']

        # vehicle can't fall below ground, remember z is down
        if self.groundcheck and z[2] > 0:
            z[0] = 0
            dz[0] = 0

        do = np.linalg.inv(model['J']) @ (np.cross(-o, model['J'] @ o) + np.sum(tau, axis=1) + np.sum(Q, axis=1)) # row sum of torques

        # # stash the flapping information for plotting
        # self.a1s = a1s
        # self.b1s = b1s

        return np.r_[dz, dn, dv, do]  # This is the state derivative vector

# ------------------------------------------------------------------------ #

class MultiRotorMixer(FunctionBlock):
    """
    :blockname:`MULTIROTORMIXER`

    .. table::
       :align: left

       +--------+---------+---------+
       | inputs | outputs |  states |
       +--------+---------+---------+
       | 4      | 1       | 0       |
       +--------+---------+---------+
       | float  |         |         |
       +--------+---------+---------+
    """

    nin = 4
    nout = 1
    inlabels = ('𝛕r', '𝛕p', '𝛕y', 'T')
    outlabels = ('ω',)

    def __init__(self, maxw=1000, minw=5, **kwargs):
        """
        Create a block that displays/animates a multi-rotor flying vehicle.

        :param maxw: maximum rotor speed in rad/s, defaults to 1000
        :type maxw: float
        :param minw: minimum rotor speed in rad/s, defaults to 5
        :type minw: float
        :param ``**kwargs``: common Block options
        :return: a MULTIROTORMIXER block
        :rtype: MultiRotorMixer instance

        **Block ports**

            :input 𝛕r: roll torque
            :input 𝛕p: pitch torque
            :input 𝛕y: yaw torque
            :input T: total thrust

            :output ω: 1D array of rotor speeds

        Derived from Simulink model by Pauline Pounds 2004
        """
        super().__init__(inputs=inputs, **kwargs)
        self.type = 'multirotormixer'
        self.minw = minw
        self.maxw = maxw

    def output(self, t):
        w = np.zeros((self.nrotors,))
        tau = self.inputs
        for i in self.nrotors:
            # roll and pitch coupling
            w[i] += -tau[0] * sin(self.theta[i]) + tau[1] * cos(self.theta[i])

            # yaw coupling
            sign = 1 if i % 1 == 0 else -1
            w[i] += sign * tau[2]

            # overall thrust
            w[i] += tau[3] / self.nrotors

        # clip the rotor speeds to the range [minw, maxw]
        w = np.clip(w, self.minw, self.maxw)

        # convert to thrust
        w = np.sqrt(w) / self.model['b']

        # negate alterate rotors to indicate counter-rotation
        for i in self.nrotors:
            if i % 1 == 0:
                w[i] = -w[i]

        return [w]

# ------------------------------------------------------------------------ #

class MultiRotorPlot(GraphicsBlock):
    """
    :blockname:`MULTIROTORPLOT`

    .. table::
       :align: left

       +--------+---------+---------+
       | inputs | outputs |  states |
       +--------+---------+---------+
       | 1      | 0       | 0       |
       +--------+---------+---------+
       | dict   |         |         |
       +--------+---------+---------+
    """

    nin = 1
    nout = 0
    inlabels = ('x',)

    # Based on code lovingly coded by Paul Pounds, first coded 17/4/02
    # version 2 2004 added scaling and ground display
    # version 3 2010 improved rotor rendering and fixed mirroring bug

    # Displays X-4 flyer position and attitude in a 3D plot.
    # GREEN ROTOR POINTS NORTH
    # BLUE ROTOR POINTS EAST

    # PARAMETERS
    # s defines the plot size in meters
    # swi controls flyer attitude plot; 1 = on, otherwise off.

    # INPUTS
    # 1 Center X position
    # 2 Center Y position
    # 3 Center Z position
    # 4 Yaw angle in rad
    # 5 Pitch angle in rad
    # 6 Roll angle in rad

    def __init__(self, model, *inputs, scale=[-2, 2, -2, 2, 10], flapscale=1, projection='ortho', **kwargs):
        """
        Create a block that displays/animates a multi-rotor flying vehicle.

        :param model: A dictionary of vehicle geometric and inertial properties
        :type model: dict
        :param ``*inputs``: Optional incoming connections
        :type ``*inputs``: Block or Plug
        :param scale: dimensions of workspace: xmin, xmax, ymin, ymax, zmin, zmax, defaults to [-2,2,-2,2,10]
        :type scale: array_like, optional
        :param flapscale: exagerate flapping angle by this factor, defaults to 1
        :type flapscale: float
        :param projection: 3D projection, one of: 'ortho' [default], 'perspective'
        :type projection: str
        :param ``**kwargs``: common Block options
        :return: a MULTIROTORPLOT block
        :rtype: MultiRotorPlot instance

        **Block ports**

            :input y: a dictionary signal that includes the item:

                - ``x`` pose in the world frame as :math:`[x, y, z, \theta_Y, \theta_P, \theta_R]`
                - ``X`` pose in the world frame as :math:`[x, y, z, \theta_Y, \theta_P, \theta_R]`
                - ``a1s``
                - ``b1s``

        .. figure:: ../../figs/multirotorplot.png
           :width: 500px
           :alt: example of generated graphic

           Example of quad-rotor display.

        Written by Pauline Pounds 2004
        """
        super().__init__(nin=1, inputs=inputs, **kwargs)
        self.type = 'quadrotorplot'
        self.model = model
        self.scale = scale
        self.nrotors = model['nrotors']
        self.projection = projection
        self.flapscale = flapscale

    def start(self, state):
        quad = self.model

        # vehicle dimensons
        d = quad['d'];  # Hub displacement from COG
        r = quad['r'];  # Rotor radius

        #C = np.zeros((3, self.nrotors))   ## WHERE USED?
        self.D = np.zeros((3,self.nrotors))

        for i in range(0, self.nrotors):
            theta = i / self.nrotors * 2 * pi
            #  Di      Rotor hub displacements (1x3)
            # first rotor is on the x-axis, clockwise order looking down from above
            self.D[:,i] = np.r_[ quad['d'] * cos(theta), quad['d'] * sin(theta), quad['h']]

        #draw ground
        self.fig = plt.figure()
        # no axes in the figure, create a 3D axes
        self.ax = self.fig.add_subplot(111, projection='3d', proj_type=self.projection)

        # ax.set_aspect('equal')
        self.ax.set_xlabel('X')
        self.ax.set_ylabel('Y')
        self.ax.set_zlabel('-Z (height above ground)')

        self.panel = self.ax.text2D(0.05, 0.95, '', transform=self.ax.transAxes,
            fontsize=10, family='monospace', verticalalignment='top',
            bbox=dict(boxstyle='round', facecolor='white', edgecolor='black'))

        # TODO allow user to set maximum height of plot volume
        self.ax.set_xlim(self.scale[0], self.scale[1])
        self.ax.set_ylim(self.scale[2], self.scale[3])
        self.ax.set_zlim(0, self.scale[4])

        # plot the ground boundaries and the big cross
        self.ax.plot([self.scale[0], self.scale[1]], [self.scale[2], self.scale[3]], [0, 0], 'b-')
        self.ax.plot([self.scale[0], self.scale[1]], [self.scale[3], self.scale[2]], [0, 0], 'b-')
        self.ax.grid(True)

        self.shadow, = self.ax.plot([0, 0], [0, 0], 'k--')
        self.groundmark, = self.ax.plot([0], [0], [0], 'kx')

        self.arm = []
        self.disk = []
        for i in range(0, self.nrotors):
            h, = self.ax.plot([0], [0], [0])
            self.arm.append(h)
            if i == 0:
                color = 'b-'
            else:
                color = 'g-'
            h, = self.ax.plot([0], [0], [0], color)
            self.disk.append(h)

        self.a1s = np.zeros((self.nrotors,))
        self.b1s = np.zeros((self.nrotors,))

    def step(self, state):

        def plot3(h, x, y, z):
            h.set_data_3d(x, y, z)
            # h.set_data(x, y)
            # h.set_3d_properties(np.r_[z])

        # READ STATE
        z = self.inputs[0]['x'][0:3]
        n = self.inputs[0]['x'][3:6]

        # TODO, check input dimensions, 12 or 12+2N, deal with flapping

        a1s = self.inputs[0]['a1s']
        b1s = self.inputs[0]['b1s']

        quad = self.model

        # vehicle dimensons
        d = quad['d']  # Hub displacement from COG
        r = quad['r']  # Rotor radius

        # PREPROCESS ROTATION MATRIX
        phi, the, psi = n    # Euler angles

        # BBF > Inertial rotation matrix
        R = np.array([
                [cos(the) * cos(phi), sin(psi) * sin(the) * cos(phi) - cos(psi) * sin(phi), cos(psi) * sin(the) * cos(phi) + sin(psi) * sin(phi)],
                [cos(the) * sin(phi), sin(psi) * sin(the) * sin(phi) + cos(psi) * cos(phi), cos(psi) * sin(the) * sin(phi) - sin(psi)*  cos(phi)],
                [-sin(the),           sin(psi)*cos(the),                                    cos(psi) * cos(the)]
            ])

        # Manual Construction
        #Q3 = [cos(psi) -sin(psi) 0;sin(psi) cos(psi) 0;0 0 1];   %Rotation mappings
        #Q2 = [cos(the) 0 sin(the);0 1 0;-sin(the) 0 cos(the)];
        #Q1 = [1 0 0;0 cos(phi) -sin(phi);0 sin(phi) cos(phi)];
        #R = Q3*Q2*Q1;    %Rotation matrix

        # CALCULATE FLYER TIP POSITONS USING COORDINATE FRAME ROTATION
        F = np.array([
                [1,  0,  0],
                [0, -1,  0],
                [0,  0, -1]
            ])

        # Draw flyer rotors
        theta = np.linspace(0, 2 * pi, 20)
        circle = np.zeros((3, 20))
        for j, t in enumerate(theta):
            circle[:,j] = np.r_[r * sin(t), r * cos(t), 0]

        hub = np.zeros((3, self.nrotors))
        tippath = np.zeros((3, 20, self.nrotors))
        for i in range(0, self.nrotors):
            hub[:,i] = F @ (z + R @ self.D[:,i])  # points in the inertial frame

            q = self.flapscale   # Flapping angle scaling for output display - makes it easier to see what flapping is occurring
            # Rotor -> Plot frame
            Rr = np.array([
                    [cos(q * a1s[i]),  sin(q * b1s[i]) * sin(q * a1s[i]),  cos(q * b1s[i]) * sin(q * a1s[i])],
                    [0,                cos(q * b1s[i]),                   -sin(q*b1s[i])],
                    [-sin(q * a1s[i]), sin(q * b1s[i]) * cos(q * a1s[i]),  cos(q * b1s[i]) * cos(q * a1s[i])]
                ])

            tippath[:,:,i] = F @ R @ Rr @ circle
            plot3(self.disk[i], hub[0,i] + tippath[0,:,i], hub[1,i] + tippath[1,:,i], hub[2,i] + tippath[2,:,i])

        # Draw flyer
        hub0 = F @ z  # centre of vehicle
        for i in range(0, self.nrotors):
            # line from hub to centre plot3([hub(1,N) hub(1,S)],[hub(2,N) hub(2,S)],[hub(3,N) hub(3,S)],'-b')
            plot3(self.arm[i], [hub[0,i], hub0[0]], [hub[1,i], hub0[1]], [hub[2,i], hub0[2]])

            # plot a circle at the hub itself
            #plot3([hub(1,i)],[hub(2,i)],[hub(3,i)],'o')

        # plot the vehicle's centroid on the ground plane
        plot3(self.shadow, [z[0], 0], [-z[1], 0], [0, 0])
        plot3(self.groundmark, z[0], -z[1], 0)

        textstr = f"t={state.t: .2f}\nh={z[2]: .2f}\nγ={n[0]: .2f}"
        self.panel.set_text(textstr)

        super().step(state=state)

max_stars_repo_name: tassos/robotics-toolbox-python

max_issues_repo_name: tassos/robotics-toolbox-python

max_forks_repo_name: tassos/robotics-toolbox-python


text:
#This function calculates abundance based chao similarity as a metric for numpy pdist function
#Created by Jeffrey R. Smith
#Last edited September 5, 2018

#Import necessary library
import numpy as np

def chao(j, k):

	#Find species shared between the two sites
	shared	=	np.argwhere((j>=1) & (k>=1)).flatten()

	#If there are no shared species return dissimilarity of 1
	if len(shared > 0):

		#Calculate the total number of individuals in site J and K of all species
		Nj	=	np.sum(j)
		Nk	=	np.sum(k)

		#Calculate the number of individuals in site J and K of species in both sites
		Cj	=	np.sum(j[shared])
		Ck	=	np.sum(k[shared])

		#Find speices that only occur with 1 individual in site in either J or K that are also present in the other site
		jones	=	np.argwhere((j>=1) & (k==1)).flatten()
		kones	=	np.argwhere((j==1) & (k>=1)).flatten()

		#Find speices that only occur with 2 individuals in site in either J or K that are also present in the other site
		jtwos	=	np.argwhere((j>=1) & (k==2)).flatten()
		ktwos	=	np.argwhere((j==2) & (k>=1)).flatten()

		#Sum up the number of singletons
		a1j	=	len(jones)
		a1k	=	len(kones)

		#Sum up doubletons (set to one if there are none as specified in the original paper)
		a2j	=	max(len(jtwos), 1)
		a2k	=	max(len(ktwos), 1)

		#Calculate the total number of individuals in site J of species that occur with 1 individual in site K and vice versa
		Sj	=	np.sum(j[jones])
		Sk	=	np.sum(k[kones])

		#Calculate indices
		U	=	min(Cj/Nj + (Nk -1)/Nk * a1j/(2*a2j) * Sj/Nj, 1)
		V	=	min(Ck/Nk + (Nj -1)/Nj * a1k/(2*a2k) * Sk/Nk, 1)

		z = U*V / (U + V - U*V)
		return 1-z

	else:
		return 1

#Sample command to run this script
#statDistances	=	scipy.spatial.distance.squareform(scipy.spatial.distance.pdist(communityMatrix, chao))

max_stars_repo_name: jeffreysmith-jrs/general

max_issues_repo_name: jeffreysmith-jrs/general

max_forks_repo_name: jeffreysmith-jrs/general


text:
#!/usr/bin/env python3
# manual

"""
This script allows you to manually control the simulator or Duckiebot
using the keyboard arrows.
"""

import sys
import argparse
import pyglet
from pyglet.window import key
import numpy as np
import gym
import gym_duckietown
from gym_duckietown.envs import DuckietownEnv
from gym_duckietown.wrappers import UndistortWrapper

####
from PIL import Image
import cv2
import math
from apriltag import Detector
import transformations as tf
####

# from experiments.utils import save_img

parser = argparse.ArgumentParser()
parser.add_argument('--env-name', default='Duckietown')
parser.add_argument('--map-name', default='udem1')
parser.add_argument('--distortion', default=False, action='store_true')
parser.add_argument('--draw-curve', action='store_true', help='draw the lane following curve')
parser.add_argument('--draw-bbox', action='store_true', help='draw collision detection bounding boxes')
parser.add_argument('--domain-rand', action='store_true', help='enable domain randomization')
parser.add_argument('--frame-skip', default=1, type=int, help='number of frames to skip')
parser.add_argument('--seed', default=1, type=int, help='seed')
args = parser.parse_args()

if args.env_name and args.env_name.find('Duckietown') != -1:
    env = DuckietownEnv(
        seed = args.seed,
        map_name = args.map_name,
        draw_curve = args.draw_curve,
        draw_bbox = args.draw_bbox,
        domain_rand = args.domain_rand,
        frame_skip = args.frame_skip,
        distortion = args.distortion,
    )
else:
    env = gym.make(args.env_name)

#env.reset()
env.render()

@env.unwrapped.window.event
def on_key_press(symbol, modifiers):
    """
    This handler processes keyboard commands that
    control the simulation
    """

    if symbol == key.BACKSPACE or symbol == key.SLASH:
        print('RESET')
        env.reset()
        env.render()
    elif symbol == key.PAGEUP:
        env.unwrapped.cam_angle[0] = 0
    elif symbol == key.ESCAPE:
        env.close()
        sys.exit(0)

    # Take a screenshot
    # UNCOMMENT IF NEEDED - Skimage dependency
    # elif symbol == key.RETURN:
    #     print('saving screenshot')
    #     img = env.render('rgb_array')
    #     save_img('screenshot.png', img)

# Register a keyboard handler
key_handler = key.KeyStateHandler()
env.unwrapped.window.push_handlers(key_handler)
def _draw_pose(overlay, camera_params, tag_size, pose, z_sign=1):

    opoints = np.array([
        -1, -1, 0,
         1, -1, 0,
         1,  1, 0,
        -1,  1, 0,
        -1, -1, -2*z_sign,
         1, -1, -2*z_sign,
         1,  1, -2*z_sign,
        -1,  1, -2*z_sign,
    ]).reshape(-1, 1, 3) * 0.5*tag_size

    edges = np.array([
        0, 1,
        1, 2,
        2, 3,
        3, 0,
        0, 4,
        1, 5,
        2, 6,
        3, 7,
        4, 5,
        5, 6,
        6, 7,
        7, 4
    ]).reshape(-1, 2)

    fx, fy, cx, cy = camera_params

    # matriz homogenea
    K = np.array([fx, 0, cx, 0, fy, cy, 0, 0, 1]).reshape(3, 3)
    print("matriz homogenea ", K)

    rvec, _ = cv2.Rodrigues(pose[:3,:3])
    tvec = pose[:3, 3]

    dcoeffs = np.zeros(5)

    ipoints, _ = cv2.projectPoints(opoints, rvec, tvec, K, dcoeffs)
    ipoints = np.round(ipoints).astype(int)
    ipoints = [tuple(pt) for pt in ipoints.reshape(-1, 2)]
    for i, j in edges:
        cv2.line(overlay, ipoints[i], ipoints[j], (0, 255, 0), 1, 16)

def global_pose(matrix,x_ob,y_ob,angle): # matrix es la pose del apriltag x_ob e y_ob es el x e y del apriltag
    tag_size = 0.18
    tile_size = 0.585

    T_a = tf.translation_matrix([
        -x_ob, -tag_size*3/4, y_ob]) # esto ya viene multiplicado por el tile_size
    R_a = tf.euler_matrix(0,angle,0)

    T_m_a = tf.concatenate_matrices(T_a, R_a)

    # pose tag con respecto al robot
    T_r_a = np.dot(matrix, tf.euler_matrix(0, np.pi, 0))
    # pose tag con respecto al mapa
    T_a_r = np.linalg.inv(T_r_a) # T_r_a-1
    T_m_r = np.dot(T_m_a, T_a_r)

    return T_m_r

def angle2(q,angle,euler):
    return q-(angle-yaw(euler))

def l1(x,y):
    return math.sqrt(x**2,y**2)

def yaw(euler_angles):
    return euler_angles[2]

def dist(matrix):
    return np.linalg.norm([matrix[0][3],matrix[1][3],matrix[2][3]])

def update(dt):
    """
    This function is called at every frame to handle
    movement/stepping and redrawing
    """
    action = np.array([0.0, 0.0])

    if key_handler[key.UP]:
        action = np.array([0.44, 0.0])
    if key_handler[key.DOWN]:
        action = np.array([-0.44, 0])
    if key_handler[key.LEFT]:
        action = np.array([0.35, +1])
    if key_handler[key.RIGHT]:
        action = np.array([0.35, -1])
    if key_handler[key.SPACE]:
        action = np.array([0, 0])

    # Speed boost
    if key_handler[key.LSHIFT]:
        action *= 1.5

    obs, reward, done, info = env.step(action)

    ### apriltags detector
    label = ""
    original = Image.fromarray(obs)
    cv_img = cv2.cvtColor(np.array(original), cv2.COLOR_RGB2BGR)

    detector = Detector()
    gray = cv2.cvtColor(np.array(original), cv2.COLOR_RGB2GRAY)
    detections, dimg = detector.detect(gray, return_image=True)
    camera = [305.57, 308.83, 303.07, 231.88]

    pose = [[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]]
    robot_pose = [[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]]
    for detection in detections:

        pose, e0, e1 = detector.detection_pose(detection, camera, 0.18 / 2)
        if not np.isnan(pose[0][0]):
            _draw_pose(cv_img,
                       camera,
                       0.18/ 2,
                       pose)
        # aqui si tiene que ir el 0.585
        robot_pose = global_pose(pose, 2.08*0.585, 4.05*0.585, math.pi)
        print(" ROBOT POSE ", robot_pose)
    label = 'detections = %d, dist = %.2f, pos = (%.2f, %.2f)' % (len(detections), pose[2][3], robot_pose[0][3], robot_pose[2][3])

    # cv2.imshow('win', cv_img)
    # cv2.waitKey(5)

    #if done:
    #   print('done!')
    #    env.reset()
    #    env.render()
    extra_label = pyglet.text.Label(
        font_name="Arial",
        font_size=14,
        x=5,
        y=600 - 19*2
    )
    env.render(mode="top_down")
    extra_label.text = label
    extra_label.draw()

pyglet.clock.schedule_interval(update, 1.0 / env.unwrapped.frame_rate)

# Enter main event loop
pyglet.app.run()

env.close()

max_stars_repo_name: amarufd/duckiego

max_issues_repo_name: amarufd/duckiego

max_forks_repo_name: amarufd/duckiego


text:
# -*- coding: utf-8 -*-
"""
Created on Thu Jun 14 12:15:28 2018

@author: peich
"""

import numpy as np
import matplotlib.pyplot as mpl
import csv
from matplotlib import colors

def abredatos(archivo='02_201606.csv'):
    with open(archivo, 'r') as f:
        reader = csv.reader(f)
        res = list(reader)
    return res

def abredatos1(archivo='02_201606.csv'):
    res = []
    with open(archivo, 'r') as f:
        reader = csv.reader(f)
        for row in reader:
            res.append(list(row))
    return resa

def columna(ncol,archivo='02_201606.csv'):
    with open(archivo, 'r') as f:
        reader = csv.reader(f)
        res = list(reader)
    col = [list(i) for i in zip(*res)][ncol]
    del(col[0:6])
    colflo = [float(i) for i in col]
    arraycol = np.array(colflo)
    return arraycol

def sacafila(nfil):
    fila = res[nfil]
    del(fila[0:7])
    filflo = [float(i) for i in fila]
    nombrefila = np.array(filflo)
    return nombrefila

def graf1(x,y,archivo='02_201606.csv'):
    with open(archivo, 'r') as f:
        reader = csv.reader(f)
        res = list(reader)
    columnax = [list(i) for i in zip(*res)][x]
    del(columnax[0:6])
    xflo = [float(i) for i in columnax]
    xar = np.array(xflo)
    columnay = [list(i) for i in zip(*res)][y]
    del(columnay[0:6])
    yflo = [float(i) for i in columnay]
    yar = np.array(yflo)
    mpl.scatter(xar,yar)
    mpl.xlabel(res[5][x])
    mpl.ylabel(res[5][y])
    mpl.xscale('log')
    mpl.yscale('log')
    mpl.yticks([0.1,1,10,100,1000,10000,100000])
    mpl.show()

def graf2(x,y1,y2,archivo='02_201606.csv'):
    with open(archivo, 'r') as f:
        reader = csv.reader(f)
        res = list(reader)
    columnax = [list(i) for i in zip(*res)][x]
    del(columnax[0:6])
    xflo = [float(i) for i in columnax]
    xar = np.array(xflo)
    columnay1 = [list(i) for i in zip(*res)][y1]
    del(columnay1[0:6])
    y1flo = [float(i) for i in columnay1]
    y1ar = np.array(y1flo)
    columnay2 = [list(i) for i in zip(*res)][y2]
    del(columnay2[0:6])
    y2flo = [float(i) for i in columnay2]
    y2ar = np.array(y2flo)
    mpl.scatter(xar,y1ar,c=(0,1,0))
    mpl.scatter(xar,y2ar,c=(1,0,0))
    mpl.xlabel(res[5][x])
    mpl.ylabel(res[5][y1]+' (verde) y '+res[5][y2]+' (rojo)')
    mpl.xscale('log')
    mpl.yscale('log')
    mpl.yticks([0.1,1,10,100,1000,10000,100000])
    mpl.show()

def semilog(x,y):
    mpl.scatter(x,y)
    mpl.xscale('log')
    mpl.show()

def loglog(x,y):
    mpl.scatter(x,y)
    mpl.xscale('log')
    mpl.yscale('log')
    mpl.show()

max_stars_repo_name: CurroPeich/electo-things

max_issues_repo_name: CurroPeich/electo-things

max_forks_repo_name: CurroPeich/electo-things


text:
# Copyright (c) OpenMMLab. All rights reserved.
import torch
import torch.nn as nn
import torch.nn.functional as F
import numpy as np

from mmseg.core import add_prefix
from mmseg.ops import resize
from .. import builder
from ..builder import SEGMENTORS
from .base import BaseSegmentor

@SEGMENTORS.register_module()
class EncoderDecoder(BaseSegmentor):
    """Encoder Decoder segmentors.

    EncoderDecoder typically consists of backbone, decode_head, auxiliary_head.
    Note that auxiliary_head is only used for deep supervision during training,
    which could be dumped during inference.
    """

    def __init__(self,
                 backbone,
                 decode_head,
                 neck=None,
                 auxiliary_head=None,
                 train_cfg=None,
                 test_cfg=None,
                 pretrained=None,
                 init_cfg=None):
        super(EncoderDecoder, self).__init__(init_cfg)
        if pretrained is not None:
            assert backbone.get('pretrained') is None, \
                'both backbone and segmentor set pretrained weight'
            backbone.pretrained = pretrained
        self.backbone = builder.build_backbone(backbone)
        if neck is not None:
            self.neck = builder.build_neck(neck)
        self._init_decode_head(decode_head)
        self._init_auxiliary_head(auxiliary_head)

        self.train_cfg = train_cfg
        self.test_cfg = test_cfg

        assert self.with_decode_head

    def _init_decode_head(self, decode_head):
        """Initialize ``decode_head``"""
        self.decode_head = builder.build_head(decode_head)
        self.align_corners = self.decode_head.align_corners
        self.num_classes = self.decode_head.num_classes

    def _init_auxiliary_head(self, auxiliary_head):
        """Initialize ``auxiliary_head``"""
        if auxiliary_head is not None:
            if isinstance(auxiliary_head, list):
                self.auxiliary_head = nn.ModuleList()
                for head_cfg in auxiliary_head:
                    self.auxiliary_head.append(builder.build_head(head_cfg))
            else:
                self.auxiliary_head = builder.build_head(auxiliary_head)

    def extract_feat(self, img):
        """Extract features from images."""
        x = self.backbone(img)
        if self.with_neck:
            x = self.neck(x)
        return x

    def encode_decode(self, img, img_metas):
        """Encode images with backbone and decode into a semantic segmentation
        map of the same size as input."""
        x = self.extract_feat(img)
        out = self._decode_head_forward_test(x, img_metas)
        out = resize(
            input=out,
            size=img.shape[2:],
            mode='bilinear',
            align_corners=self.align_corners)
        return out

    def _decode_head_forward_train(self, x, img_metas, gt_semantic_seg):
        """Run forward function and calculate loss for decode head in
        training."""
        losses = dict()
        loss_decode = self.decode_head.forward_train(x, img_metas,
                                                     gt_semantic_seg,
                                                     self.train_cfg)

        losses.update(add_prefix(loss_decode, 'decode'))
        return losses

    def _decode_head_forward_test(self, x, img_metas):
        """Run forward function and calculate loss for decode head in
        inference."""
        seg_logits = self.decode_head.forward_test(x, img_metas, self.test_cfg)
        return seg_logits

    def _auxiliary_head_forward_train(self, x, img_metas, gt_semantic_seg):
        """Run forward function and calculate loss for auxiliary head in
        training."""
        losses = dict()
        if isinstance(self.auxiliary_head, nn.ModuleList):
            for idx, aux_head in enumerate(self.auxiliary_head):
                loss_aux = aux_head.forward_train(x, img_metas,
                                                  gt_semantic_seg,
                                                  self.train_cfg)
                losses.update(add_prefix(loss_aux, f'aux_{idx}'))
        else:
            loss_aux = self.auxiliary_head.forward_train(
                x, img_metas, gt_semantic_seg, self.train_cfg)
            losses.update(add_prefix(loss_aux, 'aux'))

        return losses

    def forward_dummy(self, img):
        """Dummy forward function."""
        seg_logit = self.encode_decode(img, None)

        return seg_logit

    def forward_train(self, img, img_metas, gt_semantic_seg):
        """Forward function for training.

        Args:
            img (Tensor): Input images.
            img_metas (list[dict]): List of image info dict where each dict
                has: 'img_shape', 'scale_factor', 'flip', and may also contain
                'filename', 'ori_shape', 'pad_shape', and 'img_norm_cfg'.
                For details on the values of these keys see
                `mmseg/datasets/pipelines/formatting.py:Collect`.
            gt_semantic_seg (Tensor): Semantic segmentation masks
                used if the architecture supports semantic segmentation task.

        Returns:
            dict[str, Tensor]: a dictionary of loss components
        """

        x = self.extract_feat(img)

        losses = dict()

        loss_decode = self._decode_head_forward_train(x, img_metas,
                                                      gt_semantic_seg)
        losses.update(loss_decode)

        if self.with_auxiliary_head:
            loss_aux = self._auxiliary_head_forward_train(
                x, img_metas, gt_semantic_seg)
            losses.update(loss_aux)

        return losses

    # TODO refactor
    def slide_inference(self, img, img_meta, rescale):
        """Inference by sliding-window with overlap.

        If h_crop > h_img or w_crop > w_img, the small patch will be used to
        decode without padding.
        """

        h_stride, w_stride = self.test_cfg.stride
        h_crop, w_crop = self.test_cfg.crop_size
        batch_size, _, h_img, w_img = img.size()
        num_classes = self.num_classes
        h_grids = max(h_img - h_crop + h_stride - 1, 0) // h_stride + 1
        w_grids = max(w_img - w_crop + w_stride - 1, 0) // w_stride + 1
        preds = img.new_zeros((batch_size, num_classes, h_img, w_img))
        count_mat = img.new_zeros((batch_size, 1, h_img, w_img))
        for h_idx in range(h_grids):
            for w_idx in range(w_grids):
                y1 = h_idx * h_stride
                x1 = w_idx * w_stride
                y2 = min(y1 + h_crop, h_img)
                x2 = min(x1 + w_crop, w_img)
                y1 = max(y2 - h_crop, 0)
                x1 = max(x2 - w_crop, 0)
                crop_img = img[:, :, y1:y2, x1:x2]
                crop_seg_logit = self.encode_decode(crop_img, img_meta)
                preds += F.pad(crop_seg_logit,
                               (int(x1), int(preds.shape[3] - x2), int(y1),
                                int(preds.shape[2] - y2)))

                count_mat[:, :, y1:y2, x1:x2] += 1
        assert (count_mat == 0).sum() == 0
        if torch.onnx.is_in_onnx_export():
            # cast count_mat to constant while exporting to ONNX
            count_mat = torch.from_numpy(
                count_mat.cpu().detach().numpy()).to(device=img.device)
        preds = preds / count_mat
        if rescale:
            preds = resize(
                preds,
                size=img_meta[0]['ori_shape'][:2],
                mode='bilinear',
                align_corners=self.align_corners,
                warning=False)
        return preds

    def whole_inference(self, img, img_meta, rescale):
        """Inference with full image."""

        seg_logit = self.encode_decode(img, img_meta)
        if rescale:
            # support dynamic shape for onnx
            if torch.onnx.is_in_onnx_export():
                size = img.shape[2:]
            else:
                size = img_meta[0]['ori_shape'][:2]
            seg_logit = resize(
                seg_logit,
                size=size,
                mode='bilinear',
                align_corners=self.align_corners,
                warning=False)

        return seg_logit

    def inference(self, img, img_meta, rescale):
        """Inference with slide/whole style.

        Args:
            img (Tensor): The input image of shape (N, 3, H, W).
            img_meta (dict): Image info dict where each dict has: 'img_shape',
                'scale_factor', 'flip', and may also contain
                'filename', 'ori_shape', 'pad_shape', and 'img_norm_cfg'.
                For details on the values of these keys see
                `mmseg/datasets/pipelines/formatting.py:Collect`.
            rescale (bool): Whether rescale back to original shape.

        Returns:
            Tensor: The output segmentation map.
        """

        assert self.test_cfg.mode in ['slide', 'whole']
        ori_shape = img_meta[0]['ori_shape']
        assert all(_['ori_shape'] == ori_shape for _ in img_meta)
        if self.test_cfg.mode == 'slide':
            seg_logit = self.slide_inference(img, img_meta, rescale)
        else:
            seg_logit = self.whole_inference(img, img_meta, rescale)
        output = F.softmax(seg_logit, dim=1)
        flip = img_meta[0]['flip']
        if flip:
            flip_direction = img_meta[0]['flip_direction']
            assert flip_direction in ['horizontal', 'vertical']
            if flip_direction == 'horizontal':
                output = output.flip(dims=(3, ))
            elif flip_direction == 'vertical':
                output = output.flip(dims=(2, ))

        return output

    def simple_test(self, img, img_meta, rescale=True, return_probabilities=False):
        """Simple test with single image."""
        seg_logit = self.inference(img, img_meta, rescale)
        if not return_probabilities:
            seg_pred = seg_logit.argmax(dim=1)
            if torch.onnx.is_in_onnx_export():
                # our inference backend only support 4D output
                seg_pred = seg_pred.unsqueeze(0)
                return seg_pred
            seg_pred = seg_pred.cpu().numpy()
            # unravel batch dim
            seg_pred = list(seg_pred)
            return seg_pred
        else:
            seg_logit = seg_logit.cpu().numpy()
            seg_logit = np.transpose(seg_logit, axes=(0, 2 ,3, 1))
            return list(seg_logit)

    def aug_test(self, imgs, img_metas, rescale=True):
        """Test with augmentations.

        Only rescale=True is supported.
        """
        # aug_test rescale all imgs back to ori_shape for now
        assert rescale
        # to save memory, we get augmented seg logit inplace
        seg_logit = self.inference(imgs[0], img_metas[0], rescale)
        for i in range(1, len(imgs)):
            cur_seg_logit = self.inference(imgs[i], img_metas[i], rescale)
            seg_logit += cur_seg_logit
        seg_logit /= len(imgs)
        seg_pred = seg_logit.argmax(dim=1)
        seg_pred = seg_pred.cpu().numpy()
        # unravel batch dim
        seg_pred = list(seg_pred)
        return seg_pred

max_stars_repo_name: russelldj/mmsegmentation

max_issues_repo_name: russelldj/mmsegmentation

max_forks_repo_name: russelldj/mmsegmentation