text:
from mod import basic_shapes as bs, easy_shaders as es, scene_graph as sg, transformations as tr
import numpy as np
from OpenGL.GL import *
from OpenGL.GL import glClearColor
import random

def grid(N, edge = 0.1):
    """
    Se calculan los puntos de la grilla (centro de cada cuadrado)
    """
    lado = (2 - edge)/N
    start = edge + lado/2 - 1
    end = 1 - edge - lado/2
    axis = np.linspace(start, end, N)
    return [lado, axis]

def enArco(v1, v2, punto):
    """
    Indica si punto está en el segmento que une v1 y v2
    """
    if v1[1] == v2[1]:
        if punto[1] != v2[1]:
            return False
        elif v2[0] > v1[0]:
            return v1[0] <= punto[0] <= v2[0]
        else:
            return v2[0] <= punto[0] <= v1[0]
    else:
        if punto[0] != v2[0]:
            return False
        elif v2[1] > v1[1]:
            return v1[1] <= punto[1] <= v2[1]
        else:
            return v2[1] <= punto[1] <= v1[1]

def enTrayecto(listaPares, punto):
    """
    Indica si punto se encuentra en el camino definido por los puntos de listaPares
    """
    n = len(listaPares)
    for i in range(n-1):
        if enArco(listaPares[i], listaPares[i + 1], punto):
            return True
    return False

def masCercano(number, arreglo):
    """
    Entrega el valor y posición del elemento de arreglo más cercano a number
    """
    array = np.abs(arreglo-number)
    valor = np.amin(array)
    index = np.where(array == valor)[0][0]
    valor = arreglo[index]
    return [valor, index]

def posInGrid(x,y,grilla):
    """
    Entrega la posición de la celda en la que se encuentra el punto x,y; según grilla
    """
    horizontal = masCercano(x,grilla)
    vertical = masCercano(y, grilla)
    return [horizontal[0],vertical[0]]

def largoDosPuntos(p1, p2):
    """
    Entrega la distancia entre 2 puntos, cuando estos tienen al menos una coordenada igual
    """
    return np.abs(np.sum(np.array(p2)-np.array(p1)))

def largo(posAd, posAt, nodos=None):
    """
    Entrega la distancia entre los puntos posAd y posAt, dada por el camino definido por nodos
    """
    if nodos is None:
        nodos = []
    n = len(nodos)
    if not n:
        return largoDosPuntos(posAd,posAt)
    else:
        large = largoDosPuntos(posAd,nodos[0])+largoDosPuntos(posAt,nodos[-1])
        for i in range(n-1):
            large += largoDosPuntos(nodos[i],nodos[i+1])
        return large

class Square:

    # Crea un cuadrado con textura

    def __init__(self,N):
        self.lado = grid(N)[0]

    def getElementoCuerpo(self, imagen):
        gpuSquare = es.toGPUShape(bs.createTextureQuad(imagen), GL_REPEAT, GL_NEAREST)

        body = sg.SceneGraphNode('body')
        body.transform = tr.matmul((tr.scale(self.lado, self.lado, 1), tr.rotationZ(-np.pi/2)))
        body.childs += [gpuSquare]
        return body

class Edge:

    # Se crea el borde del mapa

    def __init__(self,N, edge = 0.1):

        self.lado, self.grilla = grid(N)

        gpuBorde = es.toGPUShape(bs.createColorQuad(0, 0, 0))

        bord = sg.SceneGraphNode('borde')
        bord.transform = tr.scale(edge, 2, 1)
        bord.childs += [gpuBorde]

        bordeIzq = sg.SceneGraphNode('bordeIzq')
        bordeIzq.transform = tr.translate(-1+edge/2,0,0)
        bordeIzq.childs += [bord]

        bordeDer = sg.SceneGraphNode('bordeDer')
        bordeDer.transform = tr.translate(1 - edge / 2, 0, 0)
        bordeDer.childs += [bord]

        bordeSup = sg.SceneGraphNode('bordeSup')
        bordeSup.transform = tr.matmul([tr.translate(0,1-edge/2,0),tr.rotationZ(np.pi/2)])
        bordeSup.childs += [bord]

        bordeInf = sg.SceneGraphNode('bordeInf')
        bordeInf.transform = tr.matmul([tr.translate(0, -1 + edge / 2, 0), tr.rotationZ(np.pi/2)])
        bordeInf.childs += [bord]

        bordeMapa = sg.SceneGraphNode('bordeMapa')
        bordeMapa.childs = [bordeIzq, bordeDer, bordeSup, bordeInf]

        self.model = bordeMapa

    def draw(self, pipeline):
        glUseProgram(pipeline.shaderProgram)
        sg.drawSceneGraphNode(self.model, pipeline, 'transform')

class Cells:

    # Se crea el mapa cuadriculado

    def __init__(self, imagen, N, edge = 0.1):
        self.lado = grid(N)[0]

        gpuSquare = es.toGPUShape(bs.createTextureQuad(imagen, N/2, N/2), GL_REPEAT, GL_NEAREST)

        mapa = sg.SceneGraphNode('mapa')
        mapa.transform = tr.uniformScale(2-2*edge)
        mapa.childs += [gpuSquare]
        self.model = mapa

    def draw(self, pipeline):
        glUseProgram(pipeline.shaderProgram)
        glUniformMatrix4fv(glGetUniformLocation(pipeline.shaderProgram, "transform"), 1, GL_TRUE, tr.identity())
        sg.drawSceneGraphNode(self.model, pipeline, 'transform')

class GameOver:

    # Se crea la escena desplegada cuando el jugador pierde

    def __init__(self, goImage, blood, N):
        self.lado = grid(N)[0]
        gpuGameOver = es.toGPUShape(bs.createTextureQuad(goImage), GL_REPEAT, GL_NEAREST)
        gpuBlood = es.toGPUShape(bs.createTextureQuad(blood), GL_REPEAT, GL_NEAREST)

        go = sg.SceneGraphNode('go')
        go.childs += [gpuGameOver]
        self.pos = [0,0]

        blood = sg.SceneGraphNode('blood')
        blood.childs += [gpuBlood]

        escenaFinal = sg.SceneGraphNode('final')
        escenaFinal.childs += [go, blood]

        self.theta = 0
        self.size = 1
        self.muerto = False
        self.crece = True
        self.gira = True
        self.model = escenaFinal

    def update(self, dt):
        if self.muerto:
            if self.crece and self.gira:
                if self.size*0.1 >= 1.8 and np.abs(np.cos(self.theta)-1) <= 0.02 and np.abs(np.sin(self.theta)) <= 0.02:
                    self.crece = False
                    self.gira = False
                    self.theta = 0
                else:
                    self.theta += 3000*dt/(self.size**2.6)
                    self.size += 3*dt/(self.size**0.15)
        else:
            return

    def draw(self, pipeline):
        if self.muerto:
            glUseProgram(pipeline.shaderProgram)
            gameOverMove = sg.findNode(self.model, 'go')
            gameOverMove.transform = tr.matmul([tr.uniformScale(self.size*0.1),tr.rotationZ(self.theta), tr.translate(-0.015, 0, 0)])
            bloodPos = sg.findNode(self.model, 'blood')
            bloodPos.transform = tr.matmul((tr.translate(self.pos[0], self.pos[1], 0), tr.uniformScale(self.lado*3)))
            glUniformMatrix4fv(glGetUniformLocation(pipeline.shaderProgram, "transform"), 1, GL_TRUE, tr.identity())
            self.model.transform = tr.identity()
            sg.drawSceneGraphNode(self.model, pipeline, 'transform')

class Apple:

    # Se crea el modelo de manzana

    def __init__(self,N):
        self.lado = grid(N)[0]

        gpuCircle = es.toGPUShape(bs.createColorCircle(0.8, 0.05, 0.05))

        circulo = sg.SceneGraphNode('circulo')
        circulo.transform = tr.uniformScale(self.lado)
        circulo.childs += [gpuCircle]

        gpuLuz = es.toGPUShape(bs.createColorCircle(0.9, 0.7, 0.7))
        luz = sg.SceneGraphNode('luz')
        luz.transform = tr.matmul([tr.translate(-self.lado/5,self.lado/5,0),tr.uniformScale(self.lado/3)])
        luz.childs += [gpuLuz]

        gpuTallo = es.toGPUShape(bs.createColorQuad(120/255, 69/255, 19/255))
        tallo = sg.SceneGraphNode('tallo')
        tallo.transform = tr.matmul([tr.translate(0, self.lado/2 , 0), tr.scale(self.lado / 10, self.lado / 3, 1)])
        tallo.childs += [gpuTallo]

        gpuHoja = es.toGPUShape(bs.createColorCircle(50 / 255, 205 / 255, 50 / 255))
        hoja = sg.SceneGraphNode('hoja')
        sx = 5
        sy = 2
        hoja.transform = tr.matmul([tr.translate((self.lado / sx) * np.cos(np.pi / 3) + 0.006, self.lado / 2 +
                                                 self.lado / (2 * sy) * np.sin(np.pi / 3) + 0.005, 0),
                                    tr.rotationZ(-np.pi/3), tr.scale(self.lado / sx, self.lado / sy, 1)])
        hoja.childs += [gpuHoja]

        redManzana = sg.SceneGraphNode('redManzana')
        redManzana.childs += [circulo, luz, tallo, hoja]
        self.model = redManzana

class AppleCreator:

    # Se crea el objeto que se encarga de dibujar las manzanas en lugares donde no esté la serpiente

    def __init__(self,N,n):
        self.lado, self.grilla = grid(N)
        self.n = n
        self.N = N
        paresSnake = []
        for i in range(self.n+1):
            par = [self.grilla[(N + N % 2) // 2+i], self.grilla[(N + N % 2) // 2]]
            paresSnake.append(par)
        while True:
            a = random.choice(self.grilla)
            b = random.choice(self.grilla)
            pos = [a, b]
            # Si la posición escogida al azar no está ocupada por la serpiente, se detiene el loop
            if pos not in paresSnake:
                self.pos = pos
                break
        self.model = Apple(N).model
        self.posSnake = None

    def verticesSerpiente(self,cola):
        vertices = [posInGrid(cola.pos_x, cola.pos_y, self.grilla)]
        parametro = cola.lead
        while not isinstance(parametro,Head):
            n = len(parametro.giro)
            if n > 0:
                for i in range(n//4):
                    par = [parametro.giro[i*4], parametro.giro[i*4+1]]
                    if par not in vertices:
                        vertices.append(par)
            parametro = parametro.lead
        posHead = posInGrid(parametro.pos_x, parametro.pos_y, self.grilla)
        vertices.append(posHead)
        self.posSnake = vertices

    def draw(self, pipeline):
        glUseProgram(pipeline.shaderProgram)
        self.model.transform = tr.translate(self.pos[0], self.pos[1], 0)
        sg.drawSceneGraphNode(self.model, pipeline, 'transform')

    def createNew(self):
        while True:
            a = random.choice(self.grilla)
            b = random.choice(self.grilla)
            pos = [a, b]
            if not enTrayecto(self.posSnake, pos):
                self.pos = pos
                break

class BodyPart:

    theta = 0
    direccion = 'left'
    giro = []
    estado = True

    def __init__(self, N, imagen, edge = 0.1):
        self.lado, self.grilla = grid(N)
        self.N = N
        partMoving = sg.SceneGraphNode('part')
        partMoving.childs += [Square(N).getElementoCuerpo(imagen)]
        self.model = partMoving
        self.pos_x = self.grilla[(N + N % 2) // 2]
        self.pos_y = self.grilla[(N + N % 2) // 2]
        self.edge = edge

    def update(self, dt):
        if self.estado:
            self.nextPos(dt)

    def nextPos(self, dt):
        if self.direccion == 'right':
            self.pos_x = self.pos_x + dt

        elif self.direccion == 'left':
            self.pos_x = self.pos_x - dt

        elif self.direccion == 'up':
            self.pos_y = self.pos_y + dt

        else:
            self.pos_y = self.pos_y - dt

    def draw(self, pipeline):
        glUseProgram(pipeline.shaderProgram)
        glUniformMatrix4fv(glGetUniformLocation(pipeline.shaderProgram, "transform"), 1, GL_TRUE, tr.identity())
        self.model.transform = tr.matmul([tr.translate(self.pos_x, self.pos_y,0),tr.rotationZ(self.theta)])
        sg.drawSceneGraphNode(self.model, pipeline, 'transform')

class Head(BodyPart):

    celdasSnake = None
    giroPropio = None

    def __init__(self, N, imagen, edge = 0.1):
        super().__init__(N, imagen, edge)

    def direction(self, r, theta1, vertice_x, vertice_y):
        if self.estado:
            self.direccion = r
            self.theta += theta1
            self.pos_x = vertice_x
            self.pos_y = vertice_y
            self.giro += [self.pos_x,self.pos_y, r, theta1]

    def updatePart(self, parts):
        vertices = []
        vertices += [posInGrid((parts[-1]).pos_x, (parts[-1]).pos_y, self.grilla)]
        nPartes = len(parts)
        for i in range(nPartes - 2):
            n = len(parts[-(i + 2)].giroPropio)
            for j in range(n // 4):
                par = [parts[-(i + 2)].giroPropio[j * 4], parts[-(i + 2)].giroPropio[j * 4 + 1]]
                if par not in vertices:
                    vertices.append(par)
        self.celdasSnake = vertices

    def comer(self, manzana, cuerpo, n, imagen):
        if self.estado:
            if posInGrid(self.pos_x, self.pos_y, self.grilla) == manzana.pos:
                manzana.verticesSerpiente(cuerpo[n-1])
                manzana.createNew()
                cuerpo[-1].tieneSeguidor()
                cuerpo.append(TailGrow(cuerpo[-1], self.N, imagen))
                n += 1
                return n, cuerpo
            else:
                return n, cuerpo
        else:
            return n, cuerpo

    def choque(self, seAcabo):
        if not self.estado:
            return
        maxim = 1 - self.edge
        minim = self.edge - 1
        posActual = posInGrid(self.pos_x,self.pos_y, self.grilla)
        if not minim <= self.pos_y <= maxim or not minim <= self.pos_x <= maxim:
            self.estado = False
            glClearColor(0.3, 0.3, 0.3, 1.0)
            seAcabo.pos = [self.pos_x, self.pos_y]
            seAcabo.muerto = True
        if enTrayecto(self.celdasSnake, posActual):
            self.estado = False
            glClearColor(0.3, 0.3, 0.3, 1.0)
            seAcabo.pos = [self.pos_x, self.pos_y]
            seAcabo.muerto = True

class MiddlePart(BodyPart):

    def __init__(self, sucesor, N, imagen, edge = 0.1):
        super().__init__(N, imagen, edge)
        self.lead = sucesor
        self.giro = self.lead.giroPropio
        self.giroPropio = []
        self.estado = self.lead.estado

        if self.lead.direccion == 'left':
            self.pos_x = self.lead.pos_x + self.lead.lado
            self.pos_y = self.lead.pos_y

        elif self.lead.direccion == 'right':
            self.pos_x = self.lead.pos_x - self.lead.lado
            self.pos_y = self.lead.pos_y

        elif self.lead.direccion == 'up':
            self.pos_x = self.lead.pos_x
            self.pos_y = self.lead.pos_y - self.lead.lado

        elif self.lead.direccion == 'down':
            self.pos_x = self.lead.pos_x
            self.pos_y = self.lead.pos_y + self.lead.lado

        self.theta = self.lead.theta
        self.direccion = self.lead.direccion

    def turnH(self, parteGiro):
        self.direccion = parteGiro[2]
        self.theta += parteGiro[3]
        self.pos_y = parteGiro[1]
        if parteGiro[2] == 'right':
            self.pos_x = self.lead.pos_x - self.lado
        else:
            self.pos_x = self.lead.pos_x + self.lado
        self.giroPropio += parteGiro
        del self.giro[0:4]

    def turnV(self, parteGiro):
        self.direccion = parteGiro[2]
        self.theta += parteGiro[3]
        self.giroPropio += parteGiro
        if parteGiro[2] == 'up':
            self.pos_y = self.lead.pos_y - self.lado
        else:
            self.pos_y = self.lead.pos_y + self.lado
        self.pos_x = parteGiro[0]
        del self.giro[0:4]

    def change(self):
        parteGiro = self.giro[0:4]
        if self.direccion == 'up':
            if self.pos_y >= parteGiro[1]:
                self.turnH(parteGiro)

        elif self.direccion == 'down':
            if self.pos_y <= parteGiro[1]:
                self.turnH(parteGiro)

        elif self.direccion == 'left':
            if self.pos_x <= parteGiro[0]:
                self.turnV(parteGiro)

        elif self.direccion == 'right':
            if self.pos_x >= parteGiro[0]:
                self.turnV(parteGiro)
        else:
            return

    def updateDirection(self):
        if not self.lead.estado:
            self.estado = False
            return
        if len(self.giro) > 0:
            self.change()

class Neck(MiddlePart):

    def __init__(self, sucesor, N, imagen, edge = 0.1):
        super().__init__(sucesor, N, imagen, edge)
        self.giro = self.lead.giro

class Tail(MiddlePart):

    def __init__(self, sucesor, N, imagen, edge = 0.1):
        super().__init__(sucesor, N, imagen, edge)
        self.seguidor = False

    def updateDirection(self):
        if not self.lead.estado:
            self.estado = False
            return
        if len(self.giro) > 0:
            self.change()

    def turnH(self, parteGiro):
        self.direccion = parteGiro[2]
        self.theta += parteGiro[3]

        if parteGiro[2] == 'right':
            self.pos_x = self.lead.pos_x - self.lado
        else:
            self.pos_x = self.lead.pos_x + self.lado

        self.pos_y = parteGiro[1]

        if self.seguidor:
            self.giroPropio += parteGiro
        del self.giro[0:4]

    def turnV(self, parteGiro):
        self.direccion = parteGiro[2]
        self.theta += parteGiro[3]

        if parteGiro[2] == 'up':
            self.pos_y = self.lead.pos_y - self.lado
        else:
            self.pos_y = self.lead.pos_y + self.lado

        if self.seguidor:
            self.giroPropio += parteGiro
        self.pos_x = parteGiro[0]
        del self.giro[0:4]

    def tieneSeguidor(self):
        self.seguidor = True

class TailGrow(Tail):

    def __init__(self, sucesor, N, imagen, edge = 0.1):
        super().__init__(sucesor, N, imagen, edge)
        self.pos_x = sucesor.pos_x
        self.pos_y = sucesor.pos_y
        self.direccion = self.lead.direccion
        self.moverse = False
        self.theta = sucesor.theta

    def still(self):
        nNodos = len(self.giro)
        if nNodos > 0:
            nodos = []
            for i in range(nNodos // 4):
                nodo = [self.giro[i * 4], self.giro[i * 4 + 1]]
                nodos.append(nodo)
            if largo([self.lead.pos_x, self.lead.pos_y], [self.pos_x, self.pos_y], nodos) >= self.lado:
                self.moverse = True
        else:
            if largo([self.lead.pos_x, self.lead.pos_y], [self.pos_x, self.pos_y]) >= self.lado:
                self.moverse = True

    def updateDirection(self):
        if not self.lead.estado:
            self.estado = False
            return
        if not self.moverse:
            self.still()
        else:
            if len(self.giro) > 0:
                self.change()

    def update(self, dt):
        if self.moverse and self.estado:
            self.nextPos(dt)


text:
"""
pytorch has a "functional" grad API [1,2] as of v1.5

    torch.autograd.functional

in addition to

    # jax.nn and jax.experimental.stax
    torch.nn.functional

However, unlike jax, torch.autograd.functional's functions don't return
functions. Also one has to supply the function to differentiate *and* the
input.

    torch.autograd.functional.hessian(func, x) -> Tensor

whereas

    jax.hessian(func) -> hess_func
    jax.hessian(func)(x) -> hess_func(x) -> DeviceArray

Note that this torch API pattern is similar to what we find with
torch.autograd.grad():

    x=torch.rand(3)
    x.requires_grad_()
    a=x.sin().sum()
    torch.autograd.grad(a, x) -> x.grad
    # same as
    #   a.backward()
    #   x.grad

only that here the first argument is a Tensor and not a function.

resources
    https://pytorch.org/tutorials/beginner/pytorch_with_examples.html#autograd
    https://pytorch.org/tutorials/beginner/blitz/autograd_tutorial.html#sphx-glr-beginner-blitz-autograd-tutorial-py
    https://pytorch.org/docs/stable/autograd.html
    https://pytorch.org/docs/stable/notes/autograd.html

[1] https://github.com/pytorch/pytorch/commit/1f4a4aaf643b70ebcb40f388ae5226a41ca57d9b
[2] https://pytorch.org/docs/stable/autograd.html#functional-higher-level-api
"""

import torch
import numpy as np

rand = torch.rand

### gradient of scalar result c w.r.t. x, evaluated at x
### step by step, see grad_fn
##x = torch.rand(3, requires_grad=True)
####c = x.sin().pow(2.0).sum()
##a = torch.sin(x)
##print(a)
##b = torch.pow(a, 2.0)
##print(b)
##c = torch.sum(b)
##print(c)
### same as torch.autograd.grad(c,x)
##c.backward()
##print(x.grad)
##
##
### VJP: extract one row of J
##x = torch.rand(3, requires_grad=True)
##v = torch.tensor([1.0,0,0])
##b = x.sin().pow(2.0)
##b.backward(v)
##print(x.grad)

#-----------------------------------------------------------------------------
# poor man's jax-like functional API
#-----------------------------------------------------------------------------

def _wrap_input(func):
    def wrapper(_x):
        if isinstance(_x, torch.Tensor):
            x = _x
        else:
            x = torch.Tensor(np.atleast_1d(_x))
        x.requires_grad = True
        if x.grad is not None:
            x.grad.zero_()
        return func(x)
    return wrapper

@_wrap_input
def cos(x):
    return torch.cos(x)

@_wrap_input
def func(x):
    return torch.sin(x).pow(2.0).sum()

def grad(func):
    @_wrap_input
    def _gradfunc(x):
        out = func(x)
        out.backward(torch.ones_like(out))
        # x.grad is a Tensor of x.shape which holds the derivatives of func
        # w.r.t each x[i,j,k,...] evaluated at x, got it?
        return x.grad
    return _gradfunc

elementwise_grad = grad

def test():
    assert torch.allclose(grad(torch.sin)(1.234), cos(1.234))
    x = rand(10)*5 - 5
    assert torch.allclose(elementwise_grad(torch.sin)(x), torch.cos(x))
    assert grad(func)(x).shape == x.shape

    # Different grad APIs
    x1 = rand(3, requires_grad=True)
    # This is how one copies an array in pytorch
    x2 = x1.clone().detach()
    x2.requires_grad = True
    c1 = func(x1)
    # same as torch.autograd.backward(c1)
    c1.backward()
    g1 = x1.grad
    c2 = func(x2)
    g2 = torch.autograd.grad(c2, x2)[0]
    assert (g1==g2).all()

if __name__ == '__main__':
    test()

ext: py

reason: import numpy


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Updated on Mon May  3 20:01:42 2021

@author: mpolonenko; rkmaddox
"""

import numpy as np
from expyfun.io import read_wav, write_wav, write_hdf5
from glob import glob
import tempfile
import shutil
import os
from scipy.fftpack import ifft, ifftshift
from scipy.interpolate import interp1d, RegularGridInterpolator
import scipy.signal as sig
import parselmouth
import matplotlib.pyplot as plt
plt.ion()

# %% User settings

band_type = 'original'  # original (4-band) or audiological (5-band)
n_ears = 1  # 1 (diotic) or 2 (dichotic)

only_bbn = False  # if you only want to make broadband, not multiband (faster)

save_unaltered = True
save_broadband = True
save_multiband = True

overwrite_file = True
save_hdf5 = True

story = 'alchemyst'  # alchemyst (male narrator) or wrinkle (female narrator)
narrator = 'male'
main_path = '/mnt/data/abr_pip_speech_bands/'
wav_in_path = main_path + 'audio_books/{}/'.format(story)
out_path = '/mnt/data/abr_peakyspeech/stimuli/{}/'.format(story)

start_file = 0
n_trials = np.inf

# %% Set up the bands and f_shifts

fs_filt = 44100  # used for main experiments
# fs_filt = 48000  # used for pilot and all next experiments
n_filt = int(fs_filt * 5e-3)
n_filt += ((n_filt + 1) % 2)  # must be odd
freq = np.arange(n_filt) / float(n_filt) * fs_filt

# also determines HP cutoff for mixing high-freq original speech (8 kHz)
fc_band_ori = 1e3 * 2 ** np.array([-0.5, 0.5, 1.5, 2.5, 3.5])
fc_band_aud = 1e3 * 2. ** np.array([-1, 0, 1, 2, 3, 4])

if band_type == 'original':
    fc_band = fc_band_ori
    name_band = ''
else:
    fc_band = fc_band_aud
    name_band = '_aud'

if n_ears == 1:
    name_ears = '_diotic'
    f0_band = fc_band
else:
    name_ears = '_dichotic'
    f0_band = np.repeat(fc_band, n_ears)

n_band = len(fc_band) - 1
if only_bbn:
    n_f0 = n_ears
    n_f0_fake = 0
else:
    n_f0 = n_ears * n_band  # would just be n_ears for only broadband
    n_f0_fake = 6  # number of fake pulse trains to calculate common component

# %% Make band filters

def amp_fun(f, top_width, trans_width):
    # f, top_width, trans_width all assume OCTAVES
    # top and trans width should add up to 1 for complementary filters
    return np.maximum(0, np.minimum(
            1, 1 + top_width / (2 * trans_width) - np.abs(f) / trans_width))

top_width = 0.5
trans_width = 0.5

amp_band = np.zeros((n_band + 1, n_filt))
for fi, f in enumerate(fc_band):
    amp_band[fi] = amp_fun(np.log2(freq / f), top_width, trans_width) ** 0.5
amp_band[0, freq <= fc_band[0]] = 1
amp_band[-1, freq >= fc_band[-1]] = 1
amp_band[:, :-n_filt // 2:-1] = amp_band[:, 1:n_filt // 2 + 1]

h_band = ifftshift(ifft(amp_band).real, axes=-1) * \
    np.atleast_2d(sig.nuttall(n_filt))

# Make single-band filter - based off the  multiband filters
amp_single = np.copy(amp_band[-2:])
amp_single[0, freq <= fc_band[-2]] = 1
amp_single[:, :-n_filt // 2:-1] = amp_single[:, 1:n_filt // 2 + 1]
h_single = ifftshift(ifft(amp_single).real, axes=-1) * \
    np.atleast_2d(sig.nuttall(n_filt))

# plot broadband filters
plt.figure(figsize=(3.5, 4))
plt.subplot(311)
plt.semilogx(freq, amp_single.T, '-')
plt.semilogx(freq, (amp_band ** 2).sum(0), 'k:')
plt.xlim([1, 12e3])
plt.ylim([0, 1.1])
plt.ylabel('Gain')
plt.title('Design frequency response')
plt.subplot(312)
plt.plot((np.arange(n_filt) - n_filt // 2) / float(fs_filt) * 1e3, h_single.T)
plt.xlabel('Time (ms)')
plt.title('Impulse response (windowed)')
plt.subplot(313)
plt.semilogx(freq, 20 * np.log10(np.abs(np.fft.fft(h_single))).T, '-')
plt.semilogx(freq, 10 * np.log10((np.abs(np.fft.fft(h_band)) ** 2).sum(0)),
             'k:')
plt.xlim([0.1, 12e3])
plt.ylim([-60, 6])
plt.ylabel('Gain (dB)')
plt.title('Actual frequency response')
plt.tight_layout(pad=0.2)
plt.savefig(out_path + 'filters_broadband' + name_band + '.png')
plt.close()

# plot multiband filters
plt.figure(figsize=(3.5, 4))
plt.subplot(311)
plt.semilogx(freq, amp_band.T)
plt.semilogx(freq, (amp_band ** 2).sum(0), 'k:')
plt.xlim([0, 12e3])
plt.ylim([0, 1.1])
plt.ylabel('Gain')
plt.title('Design frequency response')
plt.subplot(312)
plt.plot((np.arange(n_filt) - n_filt // 2) / float(fs_filt) * 1e3, h_band.T)
plt.xlabel('Time (ms)')
plt.title('Impulse response (windowed)')
plt.subplot(313)
plt.semilogx(freq, 20 * np.log10(np.abs(np.fft.fft(h_band))).T)
plt.semilogx(freq, 10 * np.log10((np.abs(np.fft.fft(h_band)) ** 2).sum(0)),
             'k:')
plt.xlim([0.1, 12e3])
plt.ylim([-60, 6])
plt.ylabel('Gain (dB)')
plt.title('Actual frequency response')
plt.tight_layout(pad=0.2)
plt.savefig(out_path + 'filters_multiband' + name_band + '.png')
plt.close()

# %% make the stimuli

files = sorted(glob(wav_in_path + story + '*.wav'))
if n_trials == np.inf:
    n_trials = len(files)

for fi, fn in enumerate(files[start_file:n_trials]):
    if narrator == 'male':
        f0_min, f0_max = 60, 350
    elif narrator == 'female':
        f0_min, f0_max = 90, 500
    slop = 1.6
    n_ep_max = 1
    print('\n===== %s =====\n' % fn.split('/')[-1])
    print('analyzing pulses')
    temp_dir = tempfile.mkdtemp(prefix='wav_to_pulses')
    try:
        shutil.copyfile(fn, temp_dir + '/stim.wav')
        # call the praat processing
        wav_path = temp_dir + '/stim'
        # pulls actual f0 min, max to send to next praat code
        parselmouth.praat.run_file('wav_to_pitch.praat', '%s' % wav_path,
                                   str(f0_min), str(f0_max))
        with open(wav_path + '.Pitch') as fid:
            pitch = list(np.copy(fid.readlines()))[0][:-1]
        f0_min, f0_max, f_mean, f_std = [float(d) for d in pitch.split('_')]
        parselmouth.praat.run_file('wav_to_pulses.praat', '%s' % wav_path,
                                   str(f0_min), str(f0_max))
        # read in the pulse times
        with open(wav_path + '.PointProcess') as fid:
            lines = np.copy(fid.readlines())
        assert len(lines) > 7
        pulse_times = np.array([float(l.replace(' ', '').strip().split('=')[1])
                                for l in lines[7:]])
    except Exception:
        shutil.rmtree(temp_dir)
        raise Exception
    shutil.rmtree(temp_dir)
    x, fs = read_wav(fn)
    assert(fs == fs_filt)
    x = x[0]
    x *= 0.01 / x.std()  # normalize to rms = 0.01

    # %% pad stimulus to make sure long enough for f shifts
    # Pad/Concatenate zeros to improve frequency resolution
    # Desired duration is the based on the max_f of 0.05
    # Need twice the desired duration of 20s minus the stimulus duration
    stim_dur = len(x)
    desired_dur = 40 * fs

    if stim_dur < desired_dur:
        pts_rqd = desired_dur - stim_dur
        x = np.pad(x, (0, pts_rqd), 'constant', constant_values=0)

    # %% smooth the pulses
    pulse_times_fix = np.copy(pulse_times)
    n_smooth_iter = 10
    for _ in range(n_smooth_iter):
        pulse_times_fix_start = np.copy(pulse_times_fix)
        for pii in np.arange(1, len(pulse_times) - 1):
            if np.abs(np.log2(np.diff(pulse_times_fix_start[[pii - 1, pii]])) -
                      np.log2(np.diff(pulse_times_fix_start[[pii, pii + 1]])))\
                    < np.log2(slop):
                pulse_times_fix[pii] = np.mean(
                    pulse_times_fix_start[pii - 1:pii + 1 + 1])
    # unique - due to error when 2 pulse times round to the same integer
    pulse_inds = np.unique(np.round(pulse_times_fix * fs).astype(int))
    pulse_inds_unfix = np.round(pulse_times * fs).astype(int)
    # Need to ensure pulse_inds and pulse_times_fix are the same length
    pulse_times_fix = np.unique(pulse_times_fix)

    # %% find the gaps for sign reversals (helps mitigate artifact)
    b_env, a_env = sig.butter(1, 6 / (fs / 2.))
    env = sig.filtfilt(b_env, a_env, np.abs(x))
    flip_regions = np.where(env[:stim_dur] < .01 * np.median(
        env[:stim_dur]))[0]
    flip_inds = flip_regions[np.where(np.diff(flip_regions) > 1)[0] - 1]
    flip_spikes = np.zeros(env.shape)
    flip_spikes[flip_inds] = 1
    b_smooth, a_smooth = sig.butter(1, 10e3 / (fs / 2.))
    flip_sign = sig.filtfilt(b_smooth, a_smooth,
                             (-1) ** (np.cumsum(flip_spikes) +
                                      files.index(fn) % 2))

    # %% analyze the spectrogram and create mixer (when voiced/voiceless)
    print('analyzing spectrogram')
    nperseg = int(fs / 50.)  # denominator is freq in Hz
    noverlap = nperseg - int(fs / f0_max)
    f_sg, t_sg, sg = sig.spectrogram(x, fs, nperseg=nperseg, noverlap=noverlap,
                                     scaling='spectrum', window='hamming')
    phase = np.nan * np.ones((n_f0 + n_f0_fake, x.shape[-1]))
    breaks = np.where(np.diff(pulse_times_fix) >= 1. / f0_min)[0]

    if 0 not in breaks:
        # add first and last element in pulse_times_fix
        breaks = np.insert(breaks, [0], [0])
        print('Adding first element of pulse_times to breaks.')
    if pulse_times_fix.size - 1 not in breaks:
        # add first and last element in pulse_times_fix
        breaks = np.insert(breaks, [breaks.size], [pulse_times_fix.size - 1])
        print('Adding last element of pulse_times to breaks.')

    mixer = np.zeros(x.shape)
    for bii in range(len(breaks) - 1):
        inds = np.arange(pulse_inds[breaks[bii] + 1],
                         pulse_inds[breaks[bii + 1] - 1])
        if len(inds):
            try:
                ifun = interp1d(pulse_inds[breaks[bii] + 1:breaks[bii + 1]],
                                2 * np.pi * np.arange(breaks[bii + 1] -
                                                      breaks[bii] - 1),
                                kind='cubic')
                print('cubic', len(inds))
            except:
                try:
                    ifun = interp1d(pulse_inds[breaks[bii] +
                                               1:breaks[bii + 1]],
                                    2 * np.pi * np.arange(breaks[bii + 1] -
                                                          breaks[bii] - 1),
                                    kind='quadratic')
                    print('quadratic', len(inds))
                except:
                    ifun = interp1d(pulse_inds[breaks[bii] +
                                               1:breaks[bii + 1]],
                                    2 * np.pi * np.arange(breaks[bii + 1] -
                                                          breaks[bii] - 1),
                                    kind='linear')
                    print('linear', len(inds))
            phase[0, inds] = ifun(inds)
        if len(inds) > 1:
            start_inds = np.arange(pulse_inds[breaks[bii] + 1],
                                   pulse_inds[breaks[bii] + 2])
            stop_inds = np.arange(pulse_inds[breaks[bii + 1] - 2],
                                  pulse_inds[breaks[bii + 1] - 1])
            mixer[start_inds] = np.sin(np.arange(len(start_inds)) * 0.5 *
                                       np.pi / len(start_inds)) ** 2
            mixer[stop_inds] = np.sin(np.arange(len(stop_inds)) * 0.5 *
                                      np.pi / len(stop_inds))[::-1] ** 2
            mixer[start_inds[-1]:stop_inds[0]] = 1
    mixer = np.atleast_2d(1 - mixer)

    xsg, ysg = np.meshgrid(t_sg, f_sg)
    interp = RegularGridInterpolator([f_sg, t_sg], sg, method='linear',
                                     bounds_error=False, fill_value=0)
    phase_orig = np.copy(phase[0])

    # %% make the new band phases and pulse_inds
    # This part of the coded was edited to improve computation time and
    # allow shorter stimuli to be processed.
    print('make new band phases and pulse_inds')
    f_shift_max = 1
    f_shift_f_min = 0.0  # if 0 will use next highest frequency bin
    f_shift_f_max = 0.05
    f_shift_ind_min = int(np.maximum(1, np.round(f_shift_f_min * len(x) / fs)))
    f_shift_ind_max = int(np.round(f_shift_f_max * len(x) / fs))
    n_comp = f_shift_ind_max - f_shift_ind_min

    f_shift_fft = np.zeros((n_f0 + n_f0_fake, len(x)), dtype=np.complex)
    f_shift_fft[:, f_shift_ind_min:f_shift_ind_max] = np.exp(
        1j * 2 * np.pi * np.random.rand(n_f0 + n_f0_fake, n_comp))
    f_shift_fft[:, :-f_shift_ind_max:-1] = f_shift_fft[
        :, 1:f_shift_ind_max].conj()
    f_shift = ifft(f_shift_fft).real
    f_shift *= f_shift_max / np.abs(f_shift).max(axis=-1, keepdims=True)

    for f0_ind in np.arange(n_f0 + n_f0_fake):
        phase[f0_ind] = phase_orig + (
                2 * np.pi * np.cumsum(f_shift[f0_ind]) / float(fs)
                + np.random.rand() * 2 * np.pi)

    # split into true (for stimuli) & fake f0 (used to derive common component)
    phase_fake_pulses = np.copy(phase)[n_f0:]
    phase = phase[:n_f0]

    # extract f0 and t0 from the pulse inds
    f0 = np.diff(phase_orig, axis=-1) * fs / 2 / np.pi
    t0 = np.arange(phase.shape[-1] - 1) / float(fs)
    f0 = f0[pulse_inds]
    t0 = t0[pulse_inds]
    t0 = t0[np.invert(np.isnan(f0))]
    f0 = f0[np.invert(np.isnan(f0))]

    print('Computing harmonic amplitude envelopes')
    f_harm_max = np.minimum(fc_band[-1] * 2, fs / 2. - 1)
    n_harm = int(np.floor(f_harm_max / f0_min))
    amp0 = np.zeros((n_harm, len(x)))
    for hi in range(n_harm):
        f = f0 * (hi + 1)
        amp0_temp = interp(np.transpose([f, t0])) ** 0.5
        amp0_interp = interp1d(t0, amp0_temp, kind='cubic', bounds_error=False,
                               fill_value=0)
        amp0_new = np.arange(stim_dur) / float(fs)
        amp0[hi, :stim_dur] = amp0_interp(amp0_new)
    # No longer need these lines - interp already smooths the amp0 waveform
    # amp0[:, np.isnan(amp0[0])] = 0
    # bplp, aplp = sig.butter(1, 70 / (fs / 2.))
    # amp0 = np.array([sig.filtfilt(bplp, aplp, amp) for amp in amp0])

    print('Generating the harmonics')
    x_harm = np.zeros(phase.shape)
    for f0_ind in range(n_f0):
        for hi in range(n_harm):
            x_harm[f0_ind] += np.nan_to_num(np.cos(phase[f0_ind] * (hi + 1))
                                            * amp0[hi])
        """
        MODIFY THE PHASE BY A NONLINEAR FUNCTION THAT SHIFTS when
        f0 * hi + 1 is higher than a certain cutoff, e.g. 2k. this where there
        should be no beating, and the frequency transition can be broad enough
        to avoid super nasty phase skip artifacts--need to be careful about
        pulse times for reconstruction though...
        """
        x_harm[f0_ind, np.isnan(x_harm[f0_ind])] = 0
    x_harm *= (1 - mixer)
    x_harm *= ((1 - mixer) * x).std() / x_harm.std()
    x_harm *= -1
    x *= -1  # this is to make it the same sign (cond) as x_harm (inspection)

    # filter the different f0_bands
    x_harm_band = np.zeros((n_band, n_ears, len(x)))
    if not only_bbn:
        f0_ind = 0
        for band_ind in range(n_band):
            for ear_ind in range(n_ears):
                x_harm_band[band_ind, ear_ind] = sig.fftconvolve(
                    x_harm[f0_ind], h_band[band_ind], 'same')
                f0_ind += 1
    if n_ears == 1:
        x_harm_band = x_harm_band.mean(1)
    # %% combine the different bands, as well as voiced/voiceless segments

    x_harm_mix = x_harm_band.sum(0)

    x_play_band = (mixer[0] * x + x_harm_mix + (1 - mixer[0]) *
                   sig.fftconvolve(x, h_band[-1], 'same'))
    if n_ears == 1:
        x_play_single = (mixer[0] * x + sig.fftconvolve(x_harm[0],
                                                        h_single[0], 'same') +
                         (1 - mixer[0]) * sig.fftconvolve(x, h_single[1],
                                                          'same'))
    else:
        x_play_single = (mixer[0] * x + sig.fftconvolve(
            x_harm[:2], np.atleast_2d(h_single[0]), 'same') + (1 - mixer[0]) *
            sig.fftconvolve(x, h_single[1], 'same'))

    x_play_band *= flip_sign
    x_play_single *= flip_sign
    x *= flip_sign

    # %% Revert back to original length
    if stim_dur < desired_dur:
        x = x[:stim_dur]
        x_play_single = x_play_single[..., :stim_dur]
        x_play_band = x_play_band[..., :stim_dur]
        phase_orig = phase_orig[:stim_dur]
        phase = phase[..., :stim_dur]
        phase_fake_pulses = phase_fake_pulses[..., :stim_dur]
        mixer = mixer[..., :stim_dur]
        x_harm = x_harm[..., :stim_dur]
        x_harm_band = x_harm_band[..., :stim_dur]
        x_harm_mix = x_harm_mix[..., :stim_dur]
        flip_sign = flip_sign[..., :stim_dur]

    # %% Create pulse trains for deriving peaky speech responses
    pulse_inds = [np.where(np.diff(np.mod(phase[bi], 2 * np.pi)) < 0)[0]
                  for bi in range(n_f0)]
    fake_pulse_inds = [np.where(np.diff(np.mod(
        phase_fake_pulses[bi], 2 * np.pi)) < 0)[0] for bi in range(n_f0_fake)]

    # %% if selected only_bbn
    if only_bbn:
        x_harm_band = []
        x_harm_mix = []
        x_play_band = []
        phase_fake_pulses = []

    # %% save files
    print('saving')
    save_paths = [out_path + '{}_wav/multiband/'.format(story),
                  out_path + '{}_wav/broadband/'.format(story),
                  out_path + '{}_wav/unaltered/'.format(story),
                  out_path + '{}_hdf5/hdf5_full/'.format(story),
                  out_path + '{}_hdf5/hdf5_reduced/'.format(story)]
    for sp in save_paths:
        if not os.path.exists(sp):
            os.makedirs(sp)

    if save_multiband:
        write_wav(out_path + '{}_wav/multiband/'.format(story) +
                  fn.split('/')[-1][:-4] + '_multiband' + name_band +
                  name_ears + fn.split('/')[-1][-4:], x_play_band, fs,
                  overwrite=overwrite_file)
    if save_broadband:
        write_wav(out_path + '{}_wav/broadband/'.format(story) +
                  fn.split('/')[-1][:-4] + '_broadband' + name_band +
                  name_ears + fn.split('/')[-1][-4:], x_play_single, fs,
                  overwrite=overwrite_file)
    if save_unaltered:
        write_wav(out_path + '{}_wav/unaltered/'.format(story) +
                  fn.split('/')[-1][:-4] + '_unaltered' + name_ears +
                  fn.split('/')[-1][-4:], x, fs, overwrite=overwrite_file)
    if save_hdf5:
        write_hdf5(out_path + '{}_hdf5/hdf5_full/'.format(story) +
                   fn.split('/')[-1][:-4] + '_bands' + name_band + '.hdf5',
                   dict(
            pulse_inds=pulse_inds,
            fake_pulse_inds=fake_pulse_inds,
            pulse_inds_unfix=pulse_inds_unfix,
            mixer=mixer,
            x_harm_band=x_harm_band,
            x_harm_mix=x_harm_mix,
            x_harm=x_harm,
            x_play_band=x_play_band,
            x_play_single=x_play_single,
            x=x,
            fs=fs,
            flip_inds=flip_inds,
            flip_sign=flip_sign,
            f0_min=f0_min,
            f0_max=f0_max,
            n_smooth_iter=n_smooth_iter,
            fc_band=fc_band,
            amp_band=amp_band,
            amp_single=amp_single,
            h_band=h_band,
            h_single=h_single,
            f_harm_max=f_harm_max,
            n_harm=n_harm,
            f_shift=f_shift,
            f_shift_f_min=f_shift_f_min,
            f_shift_f_max=f_shift_f_max,
            f_shift_ind_min=f_shift_ind_min,
            f_shift_ind_max=f_shift_ind_max,
            phase_orig=phase_orig,
            phase=phase,
            phase_fake_pulses=phase_fake_pulses,
            n_f0=n_f0,
            n_f0_fake=n_f0_fake,
            slop=slop,
            top_width=top_width,
            trans_width=trans_width,
            n_filt=n_filt), overwrite=overwrite_file)

        write_hdf5(out_path + '{}_hdf5/hdf5_reduced/'.format(story) +
                   fn.split('/')[-1][:-4] + '_bands' + name_band +
                   '_reduced.hdf5', dict(
            pulse_inds=pulse_inds,
            fake_pulse_inds=fake_pulse_inds,
            mixer=mixer,
            x_play_band=x_play_band,
            x_play_single=x_play_single,
            x=x,
            fs=fs,
            fc_band=fc_band), overwrite=overwrite_file)

max_stars_repo_name: maddoxlab/peaky-speech

max_issues_repo_name: maddoxlab/peaky-speech

max_forks_repo_name: maddoxlab/peaky-speech

reason: import numpy,import scipy,from scipy


text:
import csv
import os
import pickle
import numpy, scipy.io
from collections import defaultdict
__author__ = 'diegogaleano'
__email__  = 'Diego.Galeano.2014@rhul.live.ac.uk'
__date__  = '19-10-2016'

class EasyParsers(object):
	"""
	Broad class that parser different types of files.
	"""
	def __init__(self, directory_file=''):

		self.directory_file = directory_file

	def get_data_directory(self):
		'''
		This function set the directory_file were the data folder is.
		:param filename: name of the file in data.
		'''
		my_dir = os.getcwd()
		#code_dir = os.path.abspath(os.path.join(my_dir, os.pardir))
		#main_dir = os.path.abspath(os.path.join(code_dir, os.pardir))
		#print data_dir
		data_dir = my_dir + '/data/databases/'

		return data_dir

	def set_file_directory(self, filename):
		'''
		This function set the data directory as files directory
		:param filename:
		:return:
		'''
		self.directory_file = filename

	def get_results_directory(self):
		'''
		This function set the directory_file were the results folder is.
		:param filename: name of the file in results.
		'''
		my_dir = os.getcwd()
		code_dir = os.path.abspath(os.path.join(my_dir, os.pardir))
		main_dir = os.path.abspath(os.path.join(code_dir, os.pardir))
		result_dir = main_dir + '/data/results/'

		return result_dir
	def get_images_directory(self):
		'''
		This function set the directory_file were the results folder is.
		:param filename: name of the file in results.
		'''
		my_dir = os.getcwd()
		code_dir = os.path.abspath(os.path.join(my_dir, os.pardir))
		main_dir = os.path.abspath(os.path.join(code_dir, os.pardir))
		result_dir = main_dir + '/data/images/'

		return result_dir

	def parse_csv(self, filename, quote = '|'):
		'''
		Parse the csv file in plain way.
		:return: list
		'''
		data = list()
		with open(filename, 'rb') as csvfile:
			spamreader = csv.reader(csvfile, delimiter=',', quotechar=quote)
			for row in spamreader:
				data.append(row)

		return data

	def parse_tsv(self, filename):
		'''
		Parse the tsv file in plain way.
		:return: list
		'''
		data = list()
		with open(filename, 'rt') as tsvfile:
			spamreader = csv.reader(tsvfile, delimiter='\t')
			for row in spamreader:
				data.append(row)
		return data

	def parse_csv_header(self, filename, column_id = 0):
		'''
		Parse the csv file.
		:param columnID: which column should be use as first key.
		:return: dictionary, keys = headers from first line.
		'''
		my_dict = defaultdict()
		header = list()
		# Plain parsing
		data = self.parse_csv(filename)

		# We get the others headers.
		for h in data[0]:
			header.append(h)

		for d in data[1::]:
			if d[column_id] not in my_dict:
				my_dict[d[column_id]] = defaultdict(list)

			for idx,h in enumerate(header):
				if idx != column_id:
					my_dict[d[column_id]][h].append(d[idx].replace('"',''))

		return my_dict

	def parse_tsv_header(self, filename, column_id = 0):
		'''
		Parse the csv file.
		:param columnID: which column should be use as first key.
		:return: dictionary, keys = headers from first line.
		'''
		my_dict = defaultdict()
		header = list()
		# Plain parsing
		data = self.parse_tsv(filename)

		# We get the others headers.
		for h in data[0]:
			header.append(h)

		for d in data[1::]:
			if d[column_id] not in my_dict:
				my_dict[d[column_id]] = defaultdict(list)

			for idx,h in enumerate(header):
				if idx != column_id:
					my_dict[d[column_id]][h].append(d[idx].replace('"',''))

		return my_dict

	def save_pickle(self, directory,variable_name, variable):
		'''
		Save the data in pickle format.
		:param variable_name:
		:param variable:
		:return:
		'''
		output = open(directory + variable_name, 'wb')
		pickle.dump(variable, output)
		output.close()

	def save_matlab(self, directory,variable_name, variable):
		'''
		Save the data in .mat format.
		:param variable_name:
		:param variable:
		:return:
		'''
		scipy.io.savemat(directory + variable_name,mdict ={variable_name: variable})

	def read_pickle(self,directory, variable_name):
		'''
		Read pickle files
		:param variable_name:
		:return: data.
		'''
		pkl_file = open(directory + variable_name, 'rb')
		data = pickle.load(pkl_file)
		return data

	def features_dictionary_to_npmatrix(self, my_dict, list_order, list_features, list_fcfp):
		'''
		Convert a dictionary into a np matrix according to the list that states the order of the indexes.
		Only will take the features indicated in list_features.
		:param my_dict:
		:param list_order:
		:param list_features
		:return:
		'''
		Nfeatures = len(list_features)
		Ndata = len(list_order)
		Nfcfp = len(list_fcfp)
		MatrixF1 = numpy.zeros(shape=(Ndata, Nfeatures))
		MatrixF2 = numpy.zeros(shape=(Ndata, Nfcfp))

		for idx, k in enumerate(list_order):
			for idy, f in enumerate(list_features):
				MatrixF1[idx][idy] = my_dict[k][f]

			for idy, f in enumerate(list_fcfp):
				if f in my_dict[k]['fcfp']:
					MatrixF2[idx][idy] = my_dict[k]['fcfp'][f]

		return MatrixF1, MatrixF2

	#if __name__ == '__main__':
	# Example
	# csv_parser = EasyParsers(' ')
	# csv_parser.default_data_directory('BroadScreeningMoleculesFeatures.csv')
	# data = csv_parser.parse_csv_header()

	#for k,v in data.iteritems():
	#print k,v

	#print len(data)


text:
#!/usr/bin/env python3 -u
# Copyright (c) 2017-present, Facebook, Inc.
# All rights reserved.
#
# This source code is licensed under the license found in the LICENSE file in
# the root directory of this source tree.
from __future__ import print_function

import argparse
import csv
import os

import numpy as np
import torch
from torch.autograd import Variable
import torch.backends.cudnn as cudnn
import torch.nn as nn
import torch.optim as optim
import torchvision.transforms as transforms
import torchvision.datasets as datasets
from torch.optim import lr_scheduler

PACKAGE_PARENT = '..'
SCRIPT_DIR = os.path.dirname(os.path.realpath(os.path.join(os.getcwd(), os.path.expanduser(__file__))))
os.sys.path.append(os.path.normpath(os.path.join(SCRIPT_DIR, PACKAGE_PARENT)))
from glico_model.cifar100 import get_cifar100
from models import vgg19_bn
from create_dataset_cifar100 import CompatibleDataset
from glico_model.utils import _load_lowshot_cifar, AverageMeter, get_classifier

# from utils import progress_bar

parser = argparse.ArgumentParser(description='PyTorch CIFAR10 Training')
parser.add_argument('--lr', default=0.1, type=float, help='learning rate')
parser.add_argument('--resume', '-r', action='store_true',
                    help='resume from checkpoint')
parser.add_argument('--model', default="wideresnet", type=str,
                    help='model type (default: ResNet18)')
parser.add_argument('--name', default='0', type=str, help='name of run')
parser.add_argument('--seed', default=0, type=int, help='random seed')
parser.add_argument('--batch-size', default=128, type=int, help='batch size')
parser.add_argument('--epoch', default=200, type=int,
                    help='total epochs to run')
parser.add_argument('--no-augment', dest='augment', action='store_false',
                    help='use standard augmentation (default: True)')
parser.add_argument('--decay', default=1e-4, type=float, help='weight decay')
parser.add_argument('--alpha', default=1., type=float,
                    help='mixup interpolation coefficient (default: 1)')
parser.add_argument('--shot', default=0, type=int,
                    help='shot smaple')
parser.add_argument('--pretrained', '-pre', action='store_true',
                    help='pretrained')
args = parser.parse_args()

use_cuda = torch.cuda.is_available()

best_acc = 0  # best test accuracy
start_epoch = 0  # start from epoch 0 or last checkpoint epoch

if args.seed != 0:
	torch.manual_seed(args.seed)
# Data
print('==> Preparing data..')
if args.augment:
	transform_train = transforms.Compose([
		transforms.RandomCrop(32, padding=4),
		transforms.RandomHorizontalFlip(),
		transforms.ToTensor(),
		transforms.Normalize((0.507, 0.487, 0.441), (0.267, 0.256, 0.276))
	])
else:
	transform_train = transforms.Compose([
		transforms.ToTensor(),
		transforms.Normalize((0.507, 0.487, 0.441), (0.267, 0.256, 0.276)),

	])

transform_test = transforms.Compose([
	transforms.ToTensor(),
	transforms.Normalize((0.507, 0.487, 0.441), (0.2023, 0.1994, 0.2010)),
])

# trainset = datasets.CIFAR10(root='~/data', train=True, download=False,
#                             transform=transform_train)
# trainloader = torch.utils.data.DataLoader(trainset,
#                                           batch_size=args.batch_size,
#                                           shuffle=True, num_workers=8)
#

# train_data_imgs, train_lables = _load_lowshot_cifar(data_dir="../../data/cifar-100", split="train", num_shot=args.shot)
# train_data = CompatibleDataset(train_data_imgs, train_lables, transform=transform_train)
# train_data = datasets.CIFAR100(root='../../data/', train=True, download=False,
#                             transform=transform_train)
#
# testset = datasets.CIFAR100(root='../../data/', train=False, download=False,
#                             transform=transform_test)
# testloader = torch.utils.data.DataLoader(testset, batch_size=100,
#                                          shuffle=False, num_workers=8)
# testset = _load_lowshot_cifar(data_dir="../../data/cifar-100", split="test",
#                                                           num_shot=args.shot)
cifar_dir_cs = '/cs/dataset/CIFAR/'
train_labeled_dataset, train_unlabeled_dataset, _, test_data = get_cifar100(cifar_dir_cs, n_labeled=args.shot,
                                                                            n_unlabled=10)
trainloader = torch.utils.data.DataLoader(train_labeled_dataset,
                                          batch_size=args.batch_size,
                                          shuffle=True, num_workers=8)
testloader = torch.utils.data.DataLoader(test_data, batch_size=100,
                                         shuffle=False, num_workers=8)
# Model
if args.resume:
	# Load checkpoint.
	print('==> Resuming from checkpoint..')
	assert os.path.isdir('checkpoint'), 'Error: no checkpoint directory found!'
	checkpoint = torch.load('./checkpoint/ckpt.t7' + args.name + '_'
	                        + str(args.seed))
	net = checkpoint['net']
	best_acc = checkpoint['acc']
	start_epoch = checkpoint['epoch'] + 1
	rng_state = checkpoint['rng_state']
	torch.set_rng_state(rng_state)
else:
	print('==> Building model..')
	# net = vgg19_bn(num_classes=100)
	net = get_classifier(100, args.model, args.pretrained)

if not os.path.isdir('results'):
	os.mkdir('results')
logname = ('results/log_' + net.__class__.__name__ + '_' + args.name + '_'
           + str(args.seed) + '.csv')

if use_cuda:
	net.cuda()
	net = torch.nn.DataParallel(net)
	print(torch.cuda.device_count())
	cudnn.benchmark = True
	print('Using CUDA..')

criterion = nn.CrossEntropyLoss()
optimizer = optim.SGD(net.parameters(), lr=args.lr, momentum=0.9, weight_decay=args.decay, nesterov=True)
milestones = [60, 120, 160]
scheduler = lr_scheduler.MultiStepLR(optimizer, milestones=milestones, gamma=0.2)

def mixup_data(x, y, alpha=1.0, use_cuda=True):
	'''Returns mixed inputs, pairs of targets, and lambda'''
	if alpha > 0:
		lam = np.random.beta(alpha, alpha)
	else:
		lam = 1

	batch_size = x.size()[0]
	if use_cuda:
		index = torch.randperm(batch_size).cuda()
	else:
		index = torch.randperm(batch_size)

	mixed_x = lam * x + (1 - lam) * x[index]
	y_a, y_b = y, y[index]
	return mixed_x, y_a, y_b, lam

def mixup_criterion(criterion, pred, y_a, y_b, lam):
	return lam * criterion(pred, y_a) + (1 - lam) * criterion(pred, y_b)

def train(epoch):
	print('\nEpoch: %d' % epoch)
	net.train()
	train_loss = 0
	reg_loss = 0
	correct = 0
	total = 0
	for batch_idx, (inputs, targets) in enumerate(trainloader):
		if use_cuda:
			inputs, targets = inputs.cuda(), targets.cuda()

		inputs, targets_a, targets_b, lam = mixup_data(inputs, targets,
		                                               args.alpha, use_cuda)
		inputs, targets_a, targets_b = map(Variable, (inputs,
		                                              targets_a, targets_b))
		outputs = net(inputs)
		loss = mixup_criterion(criterion, outputs, targets_a, targets_b, lam)
		train_loss += loss.item()
		_, predicted = torch.max(outputs.data, 1)
		total += targets.size(0)
		correct += (lam * predicted.eq(targets_a.data).cpu().sum().float()
		            + (1 - lam) * predicted.eq(targets_b.data).cpu().sum().float())

		optimizer.zero_grad()
		loss.backward()
		optimizer.step()

		# progress_bar(batch_idx, len(trainloader),
		#              'Loss: %.3f | Reg: %.5f | Acc: %.3f%% (%d/%d)'
		#              % (train_loss/(batch_idx+1), reg_loss/(batch_idx+1),
		#                 100.*correct/total, correct, total))
		print('Loss: %.3f | Reg: %.5f | Acc: %.3f%% (%d/%d)'
		      % (train_loss / (batch_idx + 1), reg_loss / (batch_idx + 1),
		         100. * correct / total, correct, total))
	return (train_loss / batch_idx, reg_loss / batch_idx, 100. * correct / total)

def test(epoch):
	global best_acc
	net.eval()
	test_loss = 0
	correct = 0
	total = 0
	top1 = AverageMeter('Acc@1', ':6.2f')
	top5 = AverageMeter('Acc@5', ':6.2f')
	topk = (1, 5)
	for batch_idx, (inputs, targets) in enumerate(testloader):
		if use_cuda:
			inputs, targets = inputs.cuda(), targets.cuda()
		outputs = net(inputs)

		maxk = max(topk)
		batch_size = targets.size(0)
		# target = validate_loader_consistency(batch_size, idx, target, test_data)
		_, pred = outputs.topk(maxk, 1, True, True)
		pred = pred.t()
		correct = pred.eq(targets.view(1, -1).expand_as(pred))
		res = []
		for k in topk:
			correct_k = correct[:k].view(-1).float().sum(0, keepdim=True)
			res.append(correct_k.mul_(100.0 / batch_size))
		acc = res[0].item()
		top1.update(acc)
		top5.update(res[1].item())

		# progress_bar(batch_idx, len(testloader),
		#              'Loss: %.3f | Acc: %.3f%% (%d/%d)'
		#              % (test_loss/(batch_idx+1), 100.*correct/total,
		#                 correct, total))
		print('Loss: %.3f | Acc: %.3f' % (test_loss / (batch_idx + 1), acc))
	if epoch == start_epoch + args.epoch - 1 or acc > best_acc:
		checkpoint(acc, epoch)
	if acc > best_acc:
		best_acc = acc
	# return (test_loss/batch_idx, 100. * correct / total)
	return (test_loss / batch_idx, top1.avg, top5.avg)

def checkpoint(acc, epoch):
	# Save checkpoint.
	print('Saving..')
	state = {
		'net': net,
		'acc': acc,
		'epoch': epoch,
		'rng_state': torch.get_rng_state()
	}
	if not os.path.isdir('checkpoint'):
		os.mkdir('checkpoint')
	torch.save(state, './checkpoint/ckpt.t7' + args.name + '_'
	           + str(args.seed))

def adjust_learning_rate(optimizer, epoch):
	"""decrease the learning rate at 100 and 150 epoch"""
	lr = args.lr
	if epoch >= 100:
		lr /= 10
	if epoch >= 150:
		lr /= 10
	for param_group in optimizer.param_groups:
		param_group['lr'] = lr

if not os.path.exists(logname):
	with open(logname, 'w') as logfile:
		logwriter = csv.writer(logfile, delimiter=',')
		logwriter.writerow(['epoch', 'train loss', 'reg loss', 'train acc',
		                    'test loss', 'test acc'])

scores_test_acc1_fewshot = dict()
scores_test_acc5_fewshot = dict()
for seed in [71, 132, 54]:
	torch.manual_seed(seed)
	for epoch in range(start_epoch, args.epoch):
		train_loss, reg_loss, train_acc = train(epoch)
		if epoch % 10 == 0:
			test_loss, test_acc1, test_acc5 = test(epoch)
		scheduler.step()
		# adjust_learning_rate(optimizer, epoch)
		with open(logname, 'a') as logfile:
			logwriter = csv.writer(logfile, delimiter=',')
			logwriter.writerow([epoch, train_loss, reg_loss, train_acc, test_loss,
			                    test_acc1])

	print('accuracy@1', test_acc1)
	print('accuracy@5', test_acc5)

	scores_test_acc1_fewshot[seed] = test_acc1
	scores_test_acc5_fewshot[seed] = test_acc5

	w = csv.writer(open(f"mixup_shot_{args.shot}_acc1.csv", "w+"))
	for key, val in scores_test_acc1_fewshot.items():
		w.writerow([key, val])
	w = csv.writer(open(f"mixup_shot_{args.shot}_acc5.csv", "w+"))
	for key, val in scores_test_acc5_fewshot.items():
		w.writerow([key, val])

max_stars_repo_name: iazuri/glico-learning-small-sample

max_issues_repo_name: iazuri/glico-learning-small-sample

max_forks_repo_name: iazuri/glico-learning-small-sample


text:
# PROJETO Govbot

# PARTE 1
# Operações Básicas

import pandas as pd
govbot = pd.read_csv ('govbot_v_01.csv', sep=',', encoding='utf-8')

print(govbot)

govbot.info()

govbot.head()

govbot.sort_values(by='INÍCIO', ascending=True)

govbot['PAÍS'].value_counts()

govbot['PAÍS'].value_counts(normalize=True)

# PARTE 2
# Operações com Gráficos

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

govbot['PAÍS'].value_counts().plot.bar(title='Percentual de chatbots por governo / país')

govbot['PAÍS'].value_counts().plot.bar()
plt.title('Percentual de chatbots por governo / país',fontsize=14,color='black')
plt.xlabel('Países avaliados',fontsize=10,color='green')
plt.ylabel('Níveis de adesão à tecnologia chatbot',fontsize=10,color='green')
plt.rcParams ['xtick.labelsize'] = 14
plt.rcParams ['ytick.labelsize'] = 14

govbot.boxplot(column='INÍCIO', color='blue')
plt.rcParams ['xtick.labelsize'] = 10
plt.rcParams ['ytick.labelsize'] = 10
plt.show()

sns.catplot(x="INÍCIO",kind='box',data=govbot,color='blue')
plt.title("Concentração de chatbots por ano de início de atividade")
plt.xlabel('INÍCIO')
plt.ylabel('ATIVIDADE')
plt.show()

sns.catplot(x="NÍVEL",data=govbot)
plt.title('Concentração de chatbots por nível de governo',fontsize=14,color='black')
plt.show()

sns.displot(govbot['PÚBLICO'],bins= 30,color='blue')
plt.title('Atendimento por tipo de público',fontsize=14,color='black')
plt.show()

sns.stripplot(data=govbot, x='SETOR',y='PAÍS')
plt.title('Países e concentração de chatbots governamentais por setor',fontsize=14,color='black')
plt.show()

# Conclusão

max_stars_repo_name: pedroidsg/dataprojects

max_issues_repo_name: pedroidsg/dataprojects

max_forks_repo_name: pedroidsg/dataprojects


text:
#!/usr/bin/python
from __future__ import division

import os
import sys
import tempfile
from scipy import stats
from math import isnan, isinf
from numpy import isclose
from subprocess import check_call, call

root = os.path.abspath( sys.argv[ 0 ] )
for i in range( 3 ):
    root = os.path.dirname( root )
sys.path.append( os.path.join( root, "lib" ) )
from testutil import Multitmp, ParallelTest
from util import mktemp

root = os.path.abspath( sys.argv[ 0 ] )
for i in range( 3 ):
    root = os.path.dirname( root )
sys.path.append( os.path.join( root, "lib" ) )
from testutil import ParallelTest

# Run test.py -h to get usage information

class FerretTest( ParallelTest ):
    def __init__( self, *args ):
        self.__backup = None
        ParallelTest.__init__( self, *args )

    def __del__( self ):
        if self.__backup is not None:
            call( [ "rm", "-rf", self.__backup ] )

    def _getInputDir( self ):
        return {
            "test":   "inputs/input_test",
            "tiny":   "inputs/input_dev",
            "small":  "inputs/input_small",
            "medium": "inputs/input_medium",
            "large":  "inputs/input_large",
            "huge":   "inputs/input_native",
        }[ self.size ]

    def checkArgs( self, *args ):
        result = ParallelTest.checkArgs( self, *args )
        self.__backup = tempfile.mkdtemp()
        check_call( [
            "rsync", "-a", self._getInputDir() + "/", self.__backup
        ] )
        return result

    def getCommand( self, outfile ):
        inputdir = self._getInputDir()
        cmd = [ self.exe, "%s/corel" % inputdir, "lsh", "%s/queries" % inputdir ]
        cmd += {
            "test":   [  "5",  "5", "1" ],
            "tiny":   [  "5",  "5", "1" ],
            "small":  [ "10", "20", "1" ],
            "medium": [ "10", "20", "1" ],
            "large":  [ "10", "20", "1" ],
            "huge":   [ "50", "20", "1" ],
        }[ self.size ]
        cmd += [ outfile ]
        return cmd, dict()

    def readFile( self, infile ):
        queries = list()
        with open( infile ) as fh:
            for i, line in enumerate( fh ):
                cur_rank = list()
                terms = line.split()
                if len( terms ) == 0:
                    queries.append( ( "", i, [] ) )
                    continue
                query = terms[0]
                ranks = terms[1:]
                for j, rank in enumerate( ranks ):
                    terms = rank.split(":")
                    if len( terms ) != 2:
                        cur_rank.append( ( "", j, 0.0 ) )
                        continue
                    image, rank = terms
                    try:
                        rank = float( rank )
                    except ValueError:
                        cur_rank.append( ( "", j, 0.0 ) )
                        continue
                    if not isnan( rank ):
                        cur_rank.append( ( image, j, rank ) )
                cur_rank.sort()
                queries.append( ( query, i, cur_rank ) )
        queries.sort()
        return queries

    def comm( self, l1, l2 ):
        i, j = 0, 0
        a, b = list(), list()
        a_extra, b_extra = list(), list()
        while i < len( l1 ) and j < len( l2 ):
            if l1[ i ][ 0 ] == l2[ j ][ 0 ]:
                a.append( l1[ i ][ 1: ] )
                b.append( l2[ j ][ 1: ] )
                i += 1
                j += 1
            elif l1[ i ][ 0 ] < l2[ j ][ 0 ]:
                a_extra.append( l1[ i ][ 1: ] )
                i += 1
            else:
                b_extra.append( l2[ j ][ 1: ] )
                j += 1
        while i < len( l1 ):
            a_extra.append( l1[ i ][ 1: ] )
            i += 1
        while j < len( l2 ):
            b_extra.append( l2[ j ][ 1: ] )
            j += 1
        return ( a_extra, a, b, b_extra )

    def validateCorrectness( self, outfile ):
        correctness = ParallelTest.validateCorrectness( self, outfile )
        if not correctness:
            return False
        if self.options.error:
            golden = self.getGolden()
            gold_queries = self.readFile( golden )
            errors = list()
            # Error function for missing/extra things. Max good is 0, max bad is 1
            def errorFun( missing, extra ):
                return 1 - ( 1 / ( 2 + ( 2 * missing ) ) ) - ( 1 / ( 2 + ( 2 * extra ) ) )
            for fname in outfile:
                # No output file is max error
                if not os.path.isfile(fname):
                    return False
                test_queries = dict()
                test_queries = self.readFile( fname )
                t1 = 0          # Kendall tau penalty for weights
                t2 = 0          # Kendall tau penalty for rank output order
                t3 = 0          # Kendall tau penalty for query output order
                w = 0           # Weighting error penalty
                r = 0           # Penalty for missing or extra ranks
                extra_queries, test_queries_int, gold_queries_int, missing_queries = \
                    self.comm( test_queries, gold_queries )
                # Penalty for missing or extra queries
                q = errorFun( len( missing_queries ), len( extra_queries ) )
                if len( gold_queries_int ) < 2:
                    t3 = 1
                else:
                    t3, _ = stats.kendalltau( [ x for x, _ in test_queries_int ],
                                              [ x for x, _ in gold_queries_int ] )
                    t3 = (.5) - (.5 * t3) # Change tau scale to [0-1] where 0 is good
                if isclose( t3, 0 ):
                    t3 = 0.0
                for test_query, gold_query in zip( test_queries_int, gold_queries_int ):
                    extra_ranks, test_ranks_int, gold_ranks_int, missing_ranks = \
                        self.comm( test_query[1], gold_query[1] )
                    r += errorFun( len( missing_ranks ), len( extra_ranks ) )
                    if len( gold_ranks_int ) < 2:
                        tau1 = 1
                    else:
                        tau1, _ = stats.kendalltau( [x for _, x in gold_ranks_int],
                                                    [x for _, x in test_ranks_int] )
                        tau1 = (.5) - (.5 * tau1) # Change tau scale to [0-1] where 0 is good
                    if isclose( tau1, 0 ):
                        tau1 = 0.0
                        absolute_error = sum( [ abs( a[1] - b[1] ) for a, b in
                                                zip( gold_ranks_int, test_ranks_int ) ] )
                        if isnan( absolute_error ) or isinf( absolute_error ):
                            w += 1
                        else:
                            w += absolute_error / ( absolute_error + 1 )
                    t1 += tau1
                    if len( gold_ranks_int ) < 2:
                        tau2 = 1
                    else:
                        tau2, _ = stats.kendalltau( [x for x, _ in gold_ranks_int],
                                                    [x for x, _ in test_ranks_int] )
                        tau2 = (.5) - (.5 * tau2) # Change tau scale to [0-1] where 0 is good
                    if isclose( tau2, 0 ):
                        tau2 = 0.0
                    t2 += tau2
                t1 = t1 / len( gold_queries_int )
                t2 = t2 / len( gold_queries_int )
                w = w / len( gold_queries_int )
                r = r / len( gold_queries_int )
                error = 1000 * q + 100 * r + 10 * t1 + 10 * t2 + 5 * t3 + w
                errors.append( 1 / ( error + 1 ) )
            self.error = errors
            return True
        else:
            return correctness

    def diff( self, golden, actual):
        self.error = list()
        if self.options.error:
            return True
        else:
            return ParallelTest.diff( self, golden, actual )

    def getParallelFitness( self, root, metrics ):
        results = ParallelTest.getParallelFitness( self, root, metrics )
        check_call( [
            "rsync", "-a", self.__backup + "/", self._getInputDir()
        ])
        if self.options.error:
            if results == [ [ 0 ] ]:
                return [ [ 0 ], [ 0 ] ]
            results.append( self.error )
            return results
        else:
            return results

FerretTest().run( root )

# exit non-zero so GenProg won't consider this a "repair"
exit( 1 )

max_issues_repo_name: dornja/powergauge

max_forks_repo_name: dornja/powergauge

reason: from numpy,from scipy


text:
from typing import Optional
import numpy as np

from ..state import LandmarkStateVar
from . import EvalTreeNode, Evaluator, TransformEvaluator, ComposeLandmarkEvaluator

class CameraIntrinsics:

  def __init__(self, fu: float, fv: float, cu: float, cv: float, b: float) -> None:
    """Simple class to hold the stereo camera intrinsics.
    Args:
      fu (float): Focal length in the u-coordinate (horizontal)
      fv (float): Focal length in the v-coordinate (vertical)
      cu (float): Focal center offset in the u-coordinate (horizontal)
      cv (float): Focal center offset in the v-coordinate (vertical)
      b (float): Stereo baseline
    """
    self.fu = fu
    self.fv = fv
    self.cu = cu
    self.cv = cv
    self.b = b

def camera_model(ints: CameraIntrinsics, p: np.ndarray) -> np.ndarray:
  g = np.zeros(p.shape[:-2] + (4, 1))
  g[..., 0, 0] = ints.fu * p[..., 0, 0] / p[..., 2, 0] + ints.cu
  g[..., 1, 0] = ints.fv * p[..., 1, 0] / p[..., 2, 0] + ints.cv
  g[..., 2, 0] = ints.fu * (p[..., 0, 0] - p[..., 3, 0] * ints.b) / p[..., 2, 0] + ints.cu
  g[..., 3, 0] = ints.fv * p[..., 1, 0] / p[..., 2, 0] + ints.cv
  return g

def camera_model_jac(ints: CameraIntrinsics, p: np.ndarray) -> np.ndarray:
  dgdp = np.zeros(p.shape[:-2] + (4, 4))
  dgdp[..., 0, 0] = ints.fu / p[..., 2, 0]
  dgdp[..., 0, 2] = -ints.fu * p[..., 0, 0] / (p[..., 2, 0]**2)
  dgdp[..., 1, 1] = ints.fv / p[..., 2, 0]
  dgdp[..., 1, 2] = -ints.fv * p[..., 1, 0] / (p[..., 2, 0]**2)
  dgdp[..., 2, 0] = ints.fu / p[..., 2, 0]
  dgdp[..., 2, 2] = -ints.fu * (p[..., 0, 0] - p[..., 3, 0] * ints.b) / (p[..., 2, 0]**2)
  dgdp[..., 2, 3] = -ints.fu * ints.b / p[..., 2, 0]
  dgdp[..., 3, 1] = ints.fv / p[..., 2, 0]
  dgdp[..., 3, 2] = -ints.fv * p[..., 1, 0] / (p[..., 2, 0]**2)
  return dgdp

class StereoCameraErrorEval(Evaluator):
  """Stereo camera error function evaluator."""

  def __init__(self, meas: np.ndarray, intrinsics: CameraIntrinsics, T_cam_landmark: TransformEvaluator,
               landmark: LandmarkStateVar):
    super().__init__()

    self._meas: np.ndarray = meas
    self._intrinsics: CameraIntrinsics = intrinsics
    self._error_eval: TransformEvaluator = ComposeLandmarkEvaluator(T_cam_landmark, landmark)

  def is_active(self):
    return self._error_eval.is_active()

  def evaluate(self, lhs: Optional[np.ndarray] = None):

    tree: EvalTreeNode = self._error_eval.get_eval_tree()
    point_in_cam_frame: np.ndarray = tree.value

    error = self._meas - camera_model(self._intrinsics, point_in_cam_frame)

    if lhs is None:
      return error

    lhs = -lhs @ camera_model_jac(self._intrinsics, point_in_cam_frame)
    jacs = self._error_eval.compute_jacs(lhs, tree)

    return error, jacs

max_stars_repo_name: utiasASRL/pysteam

max_issues_repo_name: utiasASRL/pysteam

max_forks_repo_name: utiasASRL/pysteam


text:
import argparse
import numpy as np
import cv2

def mse(imageA, imageB):
    prediction = cv2.imread(imageA,0)
    actual = cv2.imread(imageB,0)

    err = np.sum((prediction.astype("float") - actual.astype("float")) ** 2)
    err /= float(prediction.shape[0] * prediction.shape[1])
    print("RMSE = ", err)

    cv2.imshow('prediction',prediction)
    cv2.imshow('actual',actual)

    cv2.waitKey(10000)
    cv2.destroyAllWindows()

    return err

if __name__ == '__main__':
    parser = argparse.ArgumentParser(
        description='Simple testing funtion for Monodepthv2 models.')

    parser.add_argument('--imageA', type=str,
                        help='path to a test image or folder of images', required=True)
    parser.add_argument('--imageB', type=str,
                        help='path to a test image or folder of images', required=True)

    args = parser.parse_args()
    val = mse(args.imageA,args.imageB)
    # print(val)


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-

import sys, os
sys.path.append(os.pardir)  # 親ディレクトリのファイルをインポートするための設定
import numpy as np
from abc import abstractmethod
from alifebook_lib.simulators import VehicleSimulator

class SubsumptionModule(object):
    def __init__(self):
        super(SubsumptionModule, self).__init__()
        self.__inputs = {}
        self.__outputs = {}
        self.child_modules = {}
        self.set_active_module_name("")
        self.on_init()

    def add_child_module(self, name, module):
        self.child_modules[name] = module

    def set_inputs(self, inputs):
        self.__inputs = inputs
        for m in self.child_modules.values():
            m.set_inputs(inputs)

    def get_input(self, name):
        return self.__inputs.get(name, None)

    def set_output(self, name, val):
        self.__outputs[name] = val

    def get_output(self, name):
        return self.__outputs.get(name, None)

    def update(self):
        for m in self.child_modules.values():
            m.update()
        self.on_update()

    def set_active_module_name(self, name):
        self.__active_module_name = name

    def get_active_module_name(self):
        return self.__active_module_name

    @abstractmethod
    def on_init(self):
        pass

    @abstractmethod
    def on_update(self):
        pass

class AvoidModule(SubsumptionModule):
    def on_init(self):
        pass

    def on_update(self):
        self.set_output("left_wheel_speed",  10 + 30 * self.get_input("left_distance"))
        self.set_output("right_wheel_speed", 10 + 30 * self.get_input("right_distance"))
        self.set_active_module_name(self.__class__.__name__)

class WanderModule(SubsumptionModule):
    TURN_START_STEP = 100
    TURN_END_STEP = 180
    def on_init(self):
        self.counter = 0
        self.add_child_module('avoid', AvoidModule())

    def on_update(self):
        if self.get_input("right_distance") < 0.001 and self.get_input("left_distance") < 0.001:
            self.counter = (self.counter + 1) % self.TURN_END_STEP
        else:
            self.counter = 0

        if self.counter < self.TURN_START_STEP:
            # counterがTURN_START_STEPに達するまでは下位のモジュールを抑制しない
            self.set_output("left_wheel_speed",  self.child_modules['avoid'].get_output("left_wheel_speed"))
            self.set_output("right_wheel_speed", self.child_modules['avoid'].get_output("right_wheel_speed"))
            self.set_active_module_name(self.child_modules['avoid'].get_active_module_name())
        elif self.counter == self.TURN_START_STEP:
            # ランダムに左回りか右回りを決定して車輪の速度をセットする
            if np.random.rand() < 0.5:
                self.set_output("left_wheel_speed",  15)
                self.set_output("right_wheel_speed", 10)
            else:
                self.set_output("left_wheel_speed",  10)
                self.set_output("right_wheel_speed", 15)
            self.set_active_module_name(self.__class__.__name__)
        else:
            # counterがリセットされるまでは車輪の速度はそのまま
            pass

from t3 import T3

class ChaosWanderModule(SubsumptionModule):
    def on_init(self):
        self.add_child_module('avoid', AvoidModule())
        self.t3 = T3(omega0 = 0.9, omega1 = 0.3, epsilon = 0.1)
        self.t3.set_parameters(omega0 = np.random.rand())
        self.t3.set_parameters(omega1 = np.random.rand())

    def on_update(self):
        x, y = self.t3.next()  # update chaos dynamics
        if self.get_input("right_distance") < 0.001 and self.get_input("left_distance") < 0.001:
            # 距離センサーが触れない間はカオスを利用して動き回る
            left_wheel_speed = 50 * x
            right_wheel_speed = 50 * y
            self.set_output("left_wheel_speed",  left_wheel_speed)
            self.set_output("right_wheel_speed", right_wheel_speed)
            self.set_active_module_name(self.__class__.__name__)
        else:
            # 距離センサーが検知したら、avoid層で回避して、カオスのパラメターも変更して別の振る舞いを獲得する
            self.set_output("left_wheel_speed",  self.child_modules['avoid'].get_output("left_wheel_speed"))
            self.set_output("right_wheel_speed", self.child_modules['avoid'].get_output("right_wheel_speed"))
            self.t3.set_parameters(omega0 = np.random.rand())
            self.t3.set_parameters(omega1 = np.random.rand())
            self.set_active_module_name(self.child_modules['avoid'].get_active_module_name())

class ExploreModule(SubsumptionModule):
    def on_init(self):
        # カオスバージョンのwander moduleに切替可能
        self.add_child_module('wander', WanderModule())
        #self.add_child_module('wander', ChaosWanderModule())

    def on_update(self):
        if self.get_input('feed_touching'):
            # エサを検知したので下位のモジュールは抑制してスピードダウン
            self.set_output("left_wheel_speed",  0)
            self.set_output("right_wheel_speed", 0)
            self.set_active_module_name(self.__class__.__name__)
        else:
            # エサがない時は下位のモジュールは抑制せずにそのままアウトプットとする
            self.set_output("left_wheel_speed",  self.child_modules['wander'].get_output("left_wheel_speed"))
            self.set_output("right_wheel_speed", self.child_modules['wander'].get_output("right_wheel_speed"))
            self.set_active_module_name(self.child_modules['wander'].get_active_module_name())

######################
# change architecture
######################
controller = AvoidModule()
#controller = WanderModule()
#controller = ChaosWanderModule()  # ワンダーモジュールの内部にカオスを入れる
#controller = ExploreModule()

# simulatorの初期化 (Appendix参照)
simulator = VehicleSimulator(obstacle_num=5, feed_num=40)

while simulator:
    # 現在のセンサー情報を取得
    sensor_data = simulator.get_sensor_data()
    # サブサンプションのコントローラをアップデート
    controller.set_inputs(sensor_data)
    controller.update()

    # アクションを生成してアップデート
    left_wheel_speed  = controller.get_output('left_wheel_speed')
    right_wheel_speed = controller.get_output('right_wheel_speed')
    action = [left_wheel_speed, right_wheel_speed]
    color = (0, 0, 255)
    active_module = controller.get_active_module_name()
    if active_module is "AvoidModule":
        simulator.set_bodycolor((255, 0, 0))
    elif active_module in ("WanderModule", "ChaosWanderModule"):
        simulator.set_bodycolor((0, 255, 0))
    elif active_module in ("ExploreModule", "ChaosExploreModule"):
        simulator.set_bodycolor((0, 0, 255))
    simulator.update(action)


text:
import datetime
import os
import pandas as pd
import numpy as np
from pathlib import Path
from sklearn.model_selection import train_test_split

ROOT_DIR = Path(os.path.dirname(os.path.abspath(__file__))).parent.parent
DATA_DIR = os.path.join(ROOT_DIR, 'data')

class Dataset:

    @classmethod
    def sample_dataframe(cls, df, pct):
        if pct in np.linspace(0.1, 100, 1000):
            num_samples = int(len(df) * (pct / 100))
            return df.sample(num_samples, random_state=2145)
        else:
            raise ValueError('Percentage parameter should be between 0.1% and 100%')

    @classmethod
    def save_dataframe(cls, df, file_name, file_path, add_date=True):
        try:
            file_name_only, file_ext = file_name.split('.')
        except ValueError:
            raise ValueError('Expected file extension in file name.')
        if file_ext in ['pkl', 'csv']:
            if add_date:
                current_time = datetime.datetime.utcnow().isoformat().split('.')[0].replace('-','').replace('T','_').replace(':','')
                file_name = current_time + ' - ' + file_name

            file_name = os.path.join(DATA_DIR, file_path, file_name)
            if file_ext == 'csv':
                df.to_csv(file_name)
            else:
                df.to_pickle(file_name)
        else:
            raise ValueError('Supported file extensions: csv and pkl.')

    @classmethod
    def load_dataframe(cls, file_name, file_path):
        try:
            file_name_only, file_ext = file_name.split('.')
        except ValueError:
            raise ValueError('Expected file extension in file name.')
        if file_ext in ['pkl', 'csv']:
            file_name = os.path.join(DATA_DIR, file_path, file_name)
            if file_ext == 'csv':
                df = pd.read_csv(file_name)
            else:
                df = pd.read_pickle(file_name)
            return df
        else:
            raise ValueError('Supported file extensions: csv and pkl.')

    @classmethod
    def save_vect_corpus(cls, X, y, file_name, file_path, add_date=True):
        try:
            file_name_only, file_ext = file_name.split('.')
        except ValueError:
            raise ValueError('Expected file extension in file name.')
        if file_ext == 'npy':
            if add_date:
                current_time = datetime.datetime.utcnow().isoformat().split('.')[0].replace('-','').replace('T','_').replace(':','')
                file_name = current_time + ' - ' + file_name_only

            X_file_name = os.path.join(DATA_DIR, file_path, file_name + 'X.' + file_ext)
            y_file_name = os.path.join(DATA_DIR, file_path, file_name + 'y.' + file_ext)
            np.save(file=X_file_name, arr=X)
            np.save(file=y_file_name, arr=y)
        else:
            raise ValueError('Supported file extension: pkl.')

    @classmethod
    def load_vect_corpus(cls, X_file_name, y_file_name, file_path):
        try:
            X_file_name_only, X_file_ext = X_file_name.split('.')
            y_file_name_only, y_file_ext = y_file_name.split('.')
        except ValueError:
            raise ValueError('Expected file extension in file name.')
        if X_file_ext == 'npy' and y_file_ext == 'pkl':
            X_file_name = os.path.join(DATA_DIR, file_path, X_file_name)
            y_file_name = os.path.join(DATA_DIR, file_path, y_file_name)
            return np.load(X_file_name), np.load(y_file_name)
        else:
            raise ValueError('Supported file extension: pkl.')

    @classmethod
    def split_dataset(cls, X, y, test_size):
        return train_test_split(X, y, test_size=test_size, random_state=57, stratify=y)


text:
# -*- coding: utf-8 -*-
"""
This code was designed to plot outputs of the NetworkSedimentTransporter
landlab component.

This code plots:
    - the network, with option to color each link according to a link attribute.
    - the parcels, with option to color and size each parcel according to
    parcel attributes.

Authors: Katy Barnhart, Jon Czuba, Allison Pfeiffer
"""
import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.collections import LineCollection
from matplotlib.colors import Normalize

from landlab.plot.network_sediment_transporter.locate_parcel_xy import locate_parcel_xy
from landlab.utils.return_array import return_array_at_link

def plot_network_and_parcels(
    grid,
    parcels,
    parcel_time_index=None,
    map_buffer=0.1,
    parcel_filter=None,
    network_color=None,
    link_attribute=None,
    link_attribute_title=None,
    network_cmap="cividis",
    network_norm=None,
    network_linewidth=None,
    parcel_color=None,
    parcel_color_attribute=None,
    parcel_color_attribute_title=None,
    parcel_color_cmap="plasma",
    parcel_color_norm=None,
    parcel_size=None,
    parcel_size_attribute=None,
    parcel_size_attribute_title=None,
    parcel_size_norm=None,
    parcel_size_min=5,
    parcel_size_max=40,
    parcel_alpha=0.5,
    fig=None,
    **kwargs
):
    """Plot a river network and parcels on the river network. Intended to
    display the results of the NetworkSedimentTransporter component.

    The river network (an instance of NetworkModelGrid) is plotted either as
    straight links between grid nodes, or (if the network was created using a
    shapefile to set network topology) as sinuous lines representing the actual
    link geometry.

    The parcels (an instance of DataRecord) are represented as dot markers
    along the links, with the marker location set by parcel attribute
    `location_at_link`. The default is to plot the parcel locations at the
    last timestep in DataRecord, though any time index may be specified.

    Use of this plotting tool is described in detail in a landlab tutorial.

    Parameters
    ----------
    grid : NetworkModelGrid
        Instance of NetworkModelGrid.
    parcels : DataRecord
        Instance of Landlab DataRecord, with the same attribute requirements as
        NetworkSedimentTransporter.
    parcel_time_index : int, time index of parcels DataRecord
        Parcel time index to plot. Default is last timestep in parcels
        DataRecord.
    map_buffer : 0.1
        Increase the plot extent by at least this much (default 0.1). Note,
        b/c of axis equal, may be more.
    parcel_filter : boolean array of shape (number_of_parcels, )
        Filter to plot only a selection of the parcels.

    ########################
    ## Part 1: Network. To set the link colors provide either:

    network_color="k"
        Uniform color for network links.

    # or

    link_attribute : array or field name at link
        Value used to set link color. Categorical options not supported. Must
        be continuous.
    link_attribute_title : str
        String to use as the title, if link_attribute is a string, it is
        used as the default.
    network_cmap : "cividis"
        Name of colormap for network.
    network_norm : matplotlib color normalizer.
        https://matplotlib.org/3.1.1/tutorials/colors/colormapnorms.html
        Default is linear between min and max of link_attribute.

    # linewidth will be recognized by either link coloring option.

    network_linewidth : float
        Width of network lines (default 0.5).

    ########################
    ## Part 2: Parcels. To set the parcel color, provide either:

    parcel_color : color str
        Constant color used for parcel markers (default "k").

    # or

    parcel_color_attribute : parcel attribute name.
        Categorical options not supported. Must be continuous.
    parcel_color_attribute_title : str
        String to use as the legend title. If parcel_color_attribute is a
        string, it is used as the default.
    parcel_color_cmap : cmap str
        Name of colormap for variable parcel color (default "plasma").
    parcel_color_norm : matplotlib color normalizer.
        https://matplotlib.org/3.1.1/tutorials/colors/colormapnorms.html
        Default is linear between min and max of parcel_color_attribute.

    # for parcel size use either:

    parcel_size : float
        Marker size, in points.

    # or

    parcel_size_attribute: parcel atribute name.
        Categorical options not supported. Must be continuous.
    parcel_size_attribute_title : str
        string to use as the title, if parcel_size_attribute is a string, it is
        used as the default.
    parcel_size_norm : par
        matplotlib color normalizer.
        https://matplotlib.org/3.1.1/tutorials/colors/colormapnorms.html
        Default is linear between min and max of parcel_size_attribute.,
    parcel_size_min : float
        Specify the smallest size of the dot markers plotted, in
        units of points (default 5). Use with parcel_size_max. They will be
        aligned with the limits of parcel_size_norm.
    parcel_size_max : float
        Specify the largest size of the dot markers plotted, in
        units of points (default 40). Use with parcel_size_min. They will be
        aligned with the limits of parcel_size_norm.

    # with constant or attribute, can set parcel transparency

    parcel_alpha : float, between 0 and 1
        Specify parcel marker transparency (default 0.5).

    ################### Miscellaneous

    fig :  figure object
        Default is to create a new figure object.

    **kwargs Anything else to pass to figure creation.

    Returns
    -------
    fig :
        Figure object.

    """
    # part 0 checking and default setting.

    # only network color/linewidth provided OR link attribute.
    if (link_attribute is not None) and (network_color is not None):
        msg = "Only one of link_attribute and network_color can be provided."
        raise ValueError(msg)

    if link_attribute is None:
        network_color = network_color or "c"
        network_linewidth = network_linewidth or 0.5
        legend_link = False
    else:
        legend_link = True
        if link_attribute_title is None:
            if isinstance(link_attribute, str):
                link_attribute_title = link_attribute
            else:
                link_attribute_title = ""

    # only parcel color OR parcel_color_attribute.
    if (parcel_color_attribute is not None) and (parcel_color is not None):
        msg = "Only one of parcel_color_attribute and parcel_color can be provided."
        raise ValueError(msg)

    if parcel_color_attribute is None:
        parcel_color = parcel_color or "k"
        legend_parcel_color = False
    else:
        legend_parcel_color = True
        if parcel_color_attribute_title is None:
            parcel_color_attribute_title = parcel_color_attribute

    # only parcel size or parcel_size_attribute
    if (parcel_size_attribute is not None) and (parcel_size is not None):
        msg = "Only one of parcel_size_attribute and parcel_size can be provided."
        raise ValueError(msg)

    if parcel_size_attribute is None:
        parcel_size = parcel_size or 1.0
        legend_parcel_size = False
    else:
        legend_parcel_size = True

        if parcel_size_attribute_title is None:
            parcel_size_attribute_title = parcel_size_attribute

    # parcel time:
    # cant use standard value or default because a value of 0 is valid.
    if parcel_time_index is None:
        parcel_time_index = -1

    # Figure out whether the legend will have one, two, or three
    # parts (linewidth, parcel size, parcel color)
    n_legends = legend_link + legend_parcel_size + legend_parcel_color

    # set up figure, label and legend gridspecs.
    if fig is None:
        fig = plt.figure(**kwargs)

    spec = gridspec.GridSpec(
        ncols=1,
        nrows=3,
        left=0,
        right=1,
        top=1,
        bottom=0,
        figure=fig,
        height_ratios=[1, 0.1, 0.2],
    )
    ax = fig.add_subplot(spec[0, 0])
    if n_legends > 0:
        label_spec = spec[1, 0].subgridspec(
            ncols=2 * n_legends - 1,
            nrows=1,
            wspace=0.1,
            hspace=1,
            width_ratios=[1] + (n_legends - 1) * [0.2, 1],
        )
        legend_spec = spec[2, 0].subgridspec(
            ncols=2 * n_legends - 1,
            nrows=1,
            wspace=0.1,
            hspace=1,
            width_ratios=[1] + (n_legends - 1) * [0.2, 1],
        )
        legend_idx = 0

    # SET UP LINESEGMENTS FOR NETWORK. If polylines exist use, otherwise use
    # endpoints. Also get the ranges so a buffer can be placed around the
    # network.

    segments = []
    if "x_of_polyline" in grid.at_link:

        x_of_polylines = grid["link"]["x_of_polyline"]
        y_of_polylines = grid["link"]["y_of_polyline"]
        for x, y in zip(x_of_polylines, y_of_polylines):
            segment = np.array((np.array(x), np.array(y))).T
            segments.append(segment)

        xmin, ymin = np.concatenate(segments).min(axis=0)
        xmax, ymax = np.concatenate(segments).max(axis=0)

        xbuffer = map_buffer * (xmax - xmin)
        ybuffer = map_buffer * (ymax - ymin)

        xlim = (xmin - xbuffer, xmax + xbuffer)
        ylim = (ymin - ybuffer, ymax + ybuffer)

    else:

        for i in range(grid.size("link")):
            nal = grid.nodes_at_link[i]
            segment = np.array(
                (np.array(grid.x_of_node[nal]), np.array(grid.y_of_node[nal]))
            ).T
            segments.append(segment)

        xbuffer = map_buffer * (grid.x_of_node.max() - grid.x_of_node.min())
        ybuffer = map_buffer * (grid.y_of_node.max() - grid.y_of_node.min())
        xlim = (grid.x_of_node.min() - xbuffer, grid.x_of_node.max() + xbuffer)
        ylim = (grid.y_of_node.min() - ybuffer, grid.y_of_node.max() + ybuffer)

    # Add Linesegments and Configure.

    # if there is a link attribute.
    if link_attribute is not None:
        line_segments = LineCollection(
            segments,
            cmap=network_cmap,
            norm=network_norm,
            linewidth=network_linewidth,
            zorder=1,
        )
        line_segments.set_array(return_array_at_link(grid, link_attribute))
        ax.add_collection(line_segments)

        # create label
        lax = fig.add_subplot(label_spec[0, legend_idx])
        lax.text(
            0.5,
            0.0,
            "Line Color",
            transform=lax.transAxes,
            color="k",
            ha="center",
            va="top",
            size=plt.rcParams["font.size"] + 2,
        )
        lax.axis("off")

        # add legend.
        lax = fig.add_subplot(legend_spec[0, legend_idx])
        legend_idx += 2
        fig.colorbar(
            line_segments, cax=lax, orientation="horizontal", label=link_attribute_title
        )

    # if link values are constant.
    else:
        line_segments = LineCollection(
            segments, colors=network_color, linewidth=network_linewidth, zorder=1
        )
        ax.add_collection(line_segments)

    # Part 2: Add Parcels.
    X = np.empty(len(parcels.dataset.element_id))
    Y = np.empty(len(parcels.dataset.element_id))

    # Locate parcel XY for each parcel at a particular time
    # some aspects of this may be possible to speed up, but at minimum
    # locate_parcel_xy must be called for each link (since calculating location)
    # along link requires interpoloation.
    # if this occurs we must also ensure the parcel order is maintained b/c of
    # color and shape formatting.

    for parcel_idx in range(parcels.dataset.item_id.size):
        XY = locate_parcel_xy(grid, parcels, parcel_time_index, parcel_idx)
        X[parcel_idx] = XY[0]
        Y[parcel_idx] = XY[1]

    # plot X,Y point on delineated network and color/size point according to a
    # certain attribute of the parcel or the link in which the parcel resides

    # if a parcel color attribute is provided.
    if parcel_color_attribute is not None:
        # if this is true, then instead of none we need to get the right
        # values from the parcels and scale/normalize them correctly. At present
        # plan to support only continuous values. Can be extended to strs as
        # categorical.
        if parcel_color_attribute in parcels.dataset:
            if "time" in parcels.dataset[parcel_color_attribute].dims:
                parcel_color = parcels.dataset[parcel_color_attribute].values[
                    :, parcel_time_index
                ]
            else:
                parcel_color = parcels.dataset[parcel_color_attribute].values
        else:
            msg = "Parcel color attribute {attribute} not present in parcels.".format(
                attribute=parcel_color_attribute
            )
            raise ValueError(msg)

        if parcel_filter is not None:
            parcel_color = parcel_color[parcel_filter]

    if parcel_size_attribute is not None:
        # if this is true, then instead of none we need to get the right
        # values from the parcels and scale/normalize them correctly. At present
        # plan to support only continuous values. Can be extended to strs as
        # categorical.
        if parcel_size_attribute in parcels.dataset:
            if "time" in parcels.dataset[parcel_size_attribute].dims:
                parcel_size_values = parcels.dataset[parcel_size_attribute].values[
                    :, parcel_time_index
                ]
            else:
                parcel_size_values = parcels.dataset[parcel_size_attribute].values

            if parcel_size_norm is None:
                parcel_size_norm = Normalize(
                    vmin=parcel_size_values.min(), vmax=parcel_size_values.max()
                )

            parcel_size = parcel_size_min + (
                parcel_size_max - parcel_size_min
            ) * parcel_size_norm(parcel_size_values)
        else:
            msg = "Parcel size attribute {attribute} not present in parcels.".format(
                attribute=parcel_size_attribute
            )
            raise ValueError(msg)

        if parcel_filter is not None:
            parcel_size = parcel_size[parcel_filter]

    # add scatter, filter x and y if necessary.
    if parcel_filter is not None:
        X = X[parcel_filter]
        Y = Y[parcel_filter]

    scatter = ax.scatter(
        X,
        Y,
        s=parcel_size,
        c=parcel_color,
        alpha=parcel_alpha,
        cmap=parcel_color_cmap,
        norm=parcel_color_norm,
        zorder=2,
    )

    # create legends.
    if legend_parcel_color:
        lax = fig.add_subplot(label_spec[0, legend_idx])
        lax.text(
            0.5,
            0.0,
            "Parcel Color",
            transform=lax.transAxes,
            color="k",
            ha="center",
            va="top",
            size=plt.rcParams["font.size"] + 2,
        )
        lax.axis("off")

        lax = fig.add_subplot(legend_spec[0, legend_idx])
        legend_idx += 2
        fig.colorbar(
            scatter,
            cax=lax,
            orientation="horizontal",
            label=parcel_color_attribute_title,
        )

    if legend_parcel_size:
        lax = fig.add_subplot(label_spec[0, legend_idx])
        lax.text(
            0.5,
            0.0,
            "Parcel Size",
            transform=lax.transAxes,
            color="k",
            ha="center",
            va="top",
            size=plt.rcParams["font.size"] + 2,
        )
        lax.axis("off")

        lax = fig.add_subplot(legend_spec[0, legend_idx])
        handles, _ = scatter.legend_elements(prop="sizes", alpha=0.6)

        if len(handles) - 1 == 0:
            han = handles
            lab = [parcel_size_values.min()]
        else:
            han = handles[:: len(handles) - 1]
            lab = [parcel_size_values.min(), parcel_size_values.max()]

        lax.legend(
            han, lab, title=parcel_size_attribute_title, loc="center", frameon=False,
        )

        plt.axis("off")

    # Set the plot limits
    ax.set_xlim(xlim)
    ax.set_ylim(ylim)

    # make axes equal
    ax.axis("equal")

    return fig

max_stars_repo_name: elbeejay/landlab

max_issues_repo_name: elbeejay/landlab

max_forks_repo_name: elbeejay/landlab