text:
# -*- coding: utf-8 -*-
"""

Francisco Javier Caracuel Beltrán

VC - Visión por Computador

4º - GII - CSI - ETSIIT - UGR

Curso 2017/2018

Proyecto final - 7 - Implementar la creación de panoramas lineales con
proyección rectangular y mezcla de colores usando el algoritmo de Burt-Adelson

"""

import cv2
import numpy as np
import copy
import math
import random
from matplotlib import pyplot as plt
from matplotlib import image as mpimg

################################################################################
# Configuración general
#

# Ruta hacia las imágenes
path = "images/"

# Título para las imágenes
img_title = ""

# Número de columnas que tendrá la tabla con la salida de las imágenes
num_cols = 3

# Esquema de color que se utiliza por defecto para convertir las imágenes
cmap = cv2.COLOR_RGB2GRAY

# Se indica que se muestren la longitud completa de las matrices
np.set_printoptions(threshold=np.nan)

# Se inicializa SIFT
orb = cv2.ORB_create()
sift = cv2.xfeatures2d.SIFT_create()

#
################################################################################

################################################################################
# Utils
#

def set_c_map(imgs, cmap = cv2.COLOR_RGB2BGR):
    """
    Asigna un esquema de color a todas las imágenes que se reciben.

    Si se recibe una lista, las imágenes verán reflejadas su nuevo esquema de
    color directamente.
    Si se recibe una imagen, es necesario hacer una asignación con el resultado
    de aplicar esta función.

    ------------
    Para aplicar el color se puede utilizar:
    - cv2.COLOR_RGB2BGR (color)
    - cv2.COLOR_RGB2GRAY (escala de grises)
    - ...
    """

    # Se comprueba si el elemento que se recibe es una imagen o es una lista
    # con imágenes.
    if isinstance(imgs, list):

        # Se guarda la longitud de la lista de imágenes
        length = len(imgs)

        # Es una lista de imágenes, por lo que se recorre cada una para cambiar
        # el esquema de color
        for i in range(length):
            imgs[i] = cv2.cvtColor(imgs[i], cmap)

    else:

        # Si es una imagen se le aplica el esquema de color
        imgs = cv2.cvtColor(imgs, cmap)

    return imgs

def show_images(imgs, names = list(), cols = num_cols, title = "", gray = True):
    """
    Dada una lista de imágenes (imgs) y una lista de nombres (names), muestra en
    una tabla todas estas imágenes con su nombre correspondiente.
    Por defecto, el número de columnas que se van a mostrar es 3.
    Por defecto, el título que acompaña a cada imagen es "".
    """

    if gray:
        # Esquema de color que se utiliza por defecto en plt.imshow()
        plt.rcParams['image.cmap'] = 'gray'
    else:
        plt.rcParams['image.cmap'] = 'viridis'

    # Se guarda la cantidad de imágenes que se van a mostrar
    imgs_length = len(imgs)

    # Si la lista está vacía, se crea una con el tamaño de imágenes y rellena de
    # espacios en blanco para que no haya ningún error al mostrarlo
    if not names:
        names = [""]*imgs_length

    # Si hay menos imágenes que número de columnas se ha establecido, se
    # disminuye el número de columnas al número de imágenes y el número de filas
    # a 1.
    if imgs_length <= cols:

        cols = imgs_length
        rows = 1

    # Si hay más imágenes, el número de filas viene determinado por el número de
    # imágenes. Se debe redondear siempre al alza para poder colocar las
    # imágenes en la última fila
    else:

        rows = math.ceil(imgs_length/cols)

    # Se recorren todas las imágenes para colocarlas una a una en la posición
    # que les corresponde en la tabla
    for i, img in enumerate(imgs):

        # La imagen se recibe con flotantes, por lo que se hace el cambio a
        # enteros
        #img = copy.deepcopy(img).astype(int)

        # Se indica el número de filas y columnas y la posición que ocupa la
        # imagen actual en esa tabla
        plt.subplot(rows, cols, i+1)

        # Se cambia el esquema de color de la imagen
        #img = cv2.cvtColor(img, cmap)

        # Se indica el título que tendrá la imagen
        plt.title(title+names[i])

        # Se indica el valor que se mostrará en el eje X
        #plt.xlabel(i)

        # Se eliminan las marcas de los ejes X e Y
        plt.xticks([])
        plt.yticks([])

        # Se muestra la imagen en su lugar correspondiente
        plt.imshow(img)

    # Se visualiza la tabla con todas las imágenes
    plt.show()

def power_two(n):
    """
    Calcula el logaritmo en base 2 de un número, truncando los decimales
    """
    return int(math.log(n, 2))

def next_power_two(n):
    """
    Calcula el siguiente número potencia de 2 de otro número recibido por
    parámetro.
    """
    return pow(2, power_two(n)+1)

def convolution_c(img , kernel_x = None, kernel_y = None, sigma = 0, \
                    border = cv2.BORDER_DEFAULT, normalize = True, own = False):
    """

    Dada una imagen (img) y dos núcleos (kernel_x, kernel_y) realiza una
    convolución de la imagen utilizando dichos núcleos.
    Se aplicará sucesivamente el núcleo por todas las filas de la imagen y con
    esta transformación se vuelve a aplicar a todas las columnas.
    Si se recibe un kernel_x "None", se extrae el kernel a través de la
    función cv2.getGaussianKernel().
    Si se recibe un kernel_y "None", será el mismo kernel ya calculado en
    kernel_x.
    Se puede indicar un borde, que por defecto está deshabilitado. Nota: si se
    utiliza la propia función de convolución añadida a partir del Bonus 2, la
    opción de borde no se contempla al utilizarse por defecto un borde
    reflejado.
    Permite la opción de normalizar, que por defecto está activa. Se le puede
    enviar el parámetro normalize = False para que no la haga.
    En el Bonus 2 se desarrolla la propia función de convolución y se adapta
    en este apartado (ya que se utiliza para el resto de práctica) para poder
    utilizarla. El parámetro "own" indica si se utiliza o no la propia función
    de convolución. Por defecto no se utiliza.

    -------------------
    Opciones para el borde:
    - cv2.BORDER_REPLICATE
    - cv2.BORDER_REFLECT
    - cv2.BORDER_REFLECT_101
    - cv2.BORDER_WRAP
    - cv2.BORDER_CONSTANT
    -------------------

    Devuelve la imagen con la transformación realizada.

    """

    # Si se modifican los valores de img, se verán reflejados en el resto de
    # imágenes que hagan uso de la imagen base. Para evitar esto, se hace un
    # copiado en el que no es por referencia
    img = copy.deepcopy(img)

    # Se comprueba si se calcula el kernel
    if kernel_x is None:
        kernel_x = cv2.getGaussianKernel(6*sigma+1, sigma)

    # Se comprueba si se calcula el kernel_y
    if kernel_y is None:
        kernel_y = kernel_x

    # Se obtiene el número de filas y columnas que tiene la imagen. La función
    # shape devuelve también los canales siempre que la imagen no esté en
    # escala de grises, por lo que se comprueba cuantos valores devuelve la
    # función y en base a ello se guardan en dos o en tres variables (para
    # evitar errores de ejecución)
    shape = img.shape

    # Para no tener error en el ámbito de las filas y las columnas se declaran
    # en este punto
    # También se pueden obtener las filas como len(img) y las columnas como
    # len(img[0])
    rows = 0
    cols = 0

    # Si la imagen está en escala de grises devuelve dos valores
    if len(shape) == 2:
        rows, cols = shape

    # Si la imagen tiene un esquema de color devuelve tres valores
    elif len(shape) == 3:
        rows, cols, channels = shape

    # Es posible utilizar la función cv2.filter2D() sin tener en cuenta los
    # canales de la imagen, ya que esta función internamente se encarga de su
    # procesamiento.

    # Se recorren las filas de la imagen haciendo la convolución
    for i in range(rows):

        # Se modifica la fila "i" con la convolución, indicando la fila "i",
        # el núcleo que se quiere utilizar y la profundidad por defecto de la
        # imagen.
        # La fila de una imagen es una lista (vector) con "cols" enteros. El
        # resultado de aplicar filter2D devuelve una lista de "cols" listas,
        # por lo que al hacer la asignación da error. Se debe guardar el
        # resultado de hacer convolución y más tarde transformar la lista de
        # listas en una lista de enteros.

        ####
        # Ampliación Bonus 2
        #

        # A partir del Bonus 2, la función de convolución es propia, por lo que
        # en este punto se indica si se utiliza la propia o cv2.filter2D().
        if own:

            # La función de convolución propia devuelve una estructura
            # diferente a cv2.filter2D(). cv2.filter2D() devuelve una lista de
            # listas, la función de convolución propia devuelve una lista con
            # todos los elementos. El contenido de ambas es el mismo.
            img[i, :] = filter_2d(img[i, :], kernel_x)

        # Si no se utiliza la propia, se utiliza cv2.filter2D
        else:

            resConv = cv2.filter2D(img[i, :], -1, kernel_x, borderType=border)

            # Se guarda la lista de enteros creada en su fila correspondiente.
            # Como cv2.filter2D devuelve una lista de listas, de este modo se
            # convierte a lista de valores y se asigna a la fila correspondiente
            img[i, :] = [sublist[0] for sublist in resConv]

        #
        ####

    # Se recorren las columnas de la imagen haciendo la convolución
    for i in range(cols):

        # Se modifica la columna "i" con la convolución, indicando la columna
        # "i", el núcleo que se quiere utilizar y la profundidad por defecto de
        # la imagen
        # La columna de una imagen es una lista (vector) con "rows" enteros. El
        # resultado de aplicar filter2D devuelve una lista de "rows" listas,
        # por lo que al hacer la asignación da error. Se debe guardar el
        # resultado de hacer convolución y más tarde transformar la lista de
        # listas en una lista de enteros.

        ####
        # Ampliación Bonus 2
        #

        # A partir del Bonus 2, la función de convolución es propia, por lo que
        # en este punto se indica si se utiliza la propia o cv2.filter2D().
        if own:

            # La función de convolución propia devuelve una estructura
            # diferente a cv2.filter2D(). cv2.filter2D() devuelve una lista de
            # listas, la función de convolución propia devuelve una lista con
            # todos los elementos. El contenido de ambas es el mismo.
            img[:, i] = filter_2d(img[:, i], kernel_y)

        # Si no se utiliza la propia, se utiliza cv2.filter2D
        else:

            resConv = cv2.filter2D(img[:, i], -1, kernel_y, borderType=border)

            # Se guarda la lista de enteros creada en su fila correspondiente.
            # Como cv2.filter2D devuelve una lista de listas, de este modo se
            # convierte a lista de valores y se asigna a la fila correspondiente
            img[:, i] = [sublist[0] for sublist in resConv]

        #
        ####

    # Se normalizan los resultados calculados
    if normalize:
        cv2.normalize(img, img, 0, 255, cv2.NORM_MINMAX)

    # Se devuelve la imagen con la convolución realizada
    return img

def convolution_d(img , kernel_x = None, kernel_y = None, ksize = 3, sigma = 0,\
                    border = cv2.BORDER_DEFAULT, dx = 1, dy = 1, own = False):
    """
    Dada una imagen (img), realiza convolución con núcleo de primera derivada
    de tamaño ksize. Si se recibe un sigma, hace alisamiento con dicho sigma.
    Se puede indicar un borde, que por defecto está deshabilitado.
    Si se establece "own" a True utiliza la función de convolución propia.

    -------------------
    Opciones para el borde:
    - cv2.BORDER_REPLICATE
    - cv2.BORDER_REFLECT
    - cv2.BORDER_REFLECT_101
    - cv2.BORDER_WRAP
    - cv2.BORDER_CONSTANT
    -------------------

    Devuelve dos imágenes con la convolución con respecto a la x y con respecto
    a la y.
    """

    # Antes de hacer la convolución, se alisa para evitar ruido.
    if sigma > 0:
        img = convolution_c(img, sigma = sigma, own = own)

    # Si se modifican los valores de img, se verán reflejados en el resto de
    # imágenes que hagan uso de la imagen base. Para evitar esto, se hace un
    # copiado en el que no es por referencia
    img_x = copy.deepcopy(img)
    img_y = copy.deepcopy(img)

    # Se comprueba si se calcula el kernel_x
    if kernel_x is None:

        # Se calcula el kernel de la primera derivada con respecto a x.
        kernel_x = cv2.getDerivKernels(dx, 0, ksize)

        # Para aplicar convolución y no correlación se debe dar la vuelta al
        # kernel. Indicando -1 da la vuelta en el eje x e y, como solo tiene
        # uno, es correcto.
        #kernel_x = cv2.flip(kernel_x, -1)

    # Se comprueba si se calcula el kernel_y
    if kernel_y is None:

        # Se calcula el kernel de la primera derivada con respecto a y.
        kernel_y = cv2.getDerivKernels(0, dy, ksize)

        # Para aplicar convolución y no correlación se debe dar la vuelta al
        # kernel. Indicando -1 da la vuelta en el eje x e y, como solo tiene
        # uno, es correcto.
        #kernel_y = cv2.flip(kernel_y, -1)

    # Se hace la convolución de la 1º derivada con respecto a x.
    img_x = convolution_c(img_x, kernel_x[0], kernel_x[1], border = border, \
                                                                    own = own)

    # Se hace la convolución de la 1º derivada con respecto a y.
    img_y = convolution_c(img_y, kernel_y[0], kernel_y[1], border = border, \
                                                                    own = own)

    # Se devuelve la imagen con la convolución realizada
    return img_x, img_y

def get_block_size(sigma = 1.5):
    """
    Devuelve el tamaño de los vecinos (block_size) que se va a utilizar para
    obtener los puntos Harris. El valor se fija al valor correspondiente al uso
    de máscaras gaussianas de sigma 1.5. El tamaño de la máscara Gaussiana es
    6*1.5+1.
    """

    return 2*int(3*sigma+1)

def get_ksize(sigma = 1):
    """
    Devuelve el parámetro de apertura (ksize) que se va a utilizar para
    obtener los puntos Harris. El valor se fija al valor correspondiente al uso
    de máscaras gaussianas de sigma 1. El tamaño de la máscara Gaussiana es
    6*1+1.
    """

    return 2*int(3*sigma+1)

def selection_criteria_harris(l1, l2, k = 0.04):
    """
    Si se quiere obtener los puntos Harris de una imagen utilizando el operador
    Harris, es necesario seguir un criterio. En este caso se hará uso de los
    autovalores de M (l1, l2) y de la constante k.
    El criterio es: l1*l2 - k * (l1+l2)^2
    """

    # Se devuelve el valor dado siguiendo el criterio
    return ((l1*l2) - (k*((l1+l2)*(l1+l2))))

def is_center_local_max(data):
    """
    Recibe una matriz (data) y devuelve si el centro es el valor máximo de
    todos los elementos (True/False).
    """

    # Se guardan las dimensiones del objeto recibido para calcular su centro
    shape = data.shape

    # Se guarda el centro de los datos
    center = data[round(shape[0]/2)-1, round(shape[0]/2)-1]

    # Se devuelve si el centro es el máximo o no
    return center == np.max(data)

def not_max_suppression_harris(harris, threshold, env):
    """
    Suprime los valores no-máximos de una matriz (harris). Elimina los puntos
    que, aunque tengan un valor alto de criterio Harris, no son máximos locales
    de su entorno (env) para un tamaño de entorno prefijado.
    Permite la utilización de un umbral (threshold) para eliminar aquellos
    puntos Harris que se consideren elevados.
    Devuelve una lista con los keyPoints y su valor correspondiente.
    """

    # Se guardan las dimensiones de la imagen
    shape = harris.shape

    # Se consigue el mínimo valor para rellenar los bordes de la matriz con él
    min_harris = np.min(harris)

    # Se crea una nueva matriz donde se copiarán los puntos harris con un borde
    # de tamaño env relleno con el valor mínimo de los puntos Harris para que
    # no entre en conflicto al calcular el máximo local
    new_harris = np.ndarray(shape=(shape[0]+env+env, shape[1]+env+env))
    new_harris[:, :] = min_harris
    new_harris[env:shape[0]+env, env:shape[1]+env] = harris.copy()

    # Se crea una matriz rellena completamente a 255 con las mismas dimensiones
    # que la imagen actual
    res_harris = np.ndarray(shape=(shape[0]+env+env, shape[1]+env+env))
    res_harris[:, :] = 255

    # Los puntos Harris recibidos pueden tener multitud de valores. Si se
    # quieren obtener los más representativos se deben eliminar los que no
    # tienen un valor alto. En una zona completamente blanca pueden encontrarse
    # máximos locales, por lo que aparecerían en la imagen. Para evitar esto,
    # se indica un umbral, por el que por debajo de él no se tienen en cuenta
    # esos puntos
    threshold = 10**threshold

    # Lista con las coordenadas de los máximos locales y el valor que tiene
    # cada uno
    max_points = []

    # Se recorre cada posición de la matriz de puntos Harris (se ignora el
    # borde creado manualmente)
    for row in range(env, shape[0]+env):
        for col in range(env, shape[1]+env):

            # Si el punto actual tiene valor 255 o si su punto Harris está por
            # encima del umbral indicado se comprueba si es un máximo local
            if res_harris[row, col] == 255 and new_harris[row, col]>threshold:

                # Se obtiene el rango de datos que se deben conseguir para
                # comprobar que es máximo local, teniendo en cuenta el valor
                # de env
                row_insp_init = row - env
                row_insp_end = row + env
                col_insp_init = col - env
                col_insp_end = col + env

                # Se obtiene la matriz con los datos que rodean al punto actual
                data = new_harris[row_insp_init:row_insp_end+1, \
                                    col_insp_init:col_insp_end+1]

                # Se comprueba si el punto central es máximo local
                if is_center_local_max(data):

                    # En caso de ser máximo local, todo el rango de datos
                    # seleccionados se cambian a 0 para no volver a comprobarlos
                    res_harris[row_insp_init:row_insp_end+1, \
                                col_insp_init:col_insp_end+1] = 0

                    # Se guarda el punto actual real como máximo local
                    # utilizando la estructura cv2.KeyPoint
                    max_points.append([cv2.KeyPoint(row-env, col-env, _size=0, \
                                    _angle=0), harris[row-env, col-env]])

    # Se devuelven las coordenadas de los puntos máximos locales y su valor
    return max_points

def get_harris(img, sigma_block_size = 1.5, sigma_ksize = 1, k = 0.04, \
                threshold = -10000, env = 5, scale = -1):
    """
    Obtiene una lista potencial de los puntos Harris de una imagen (img).
    Los valores de los parámetros utilizados dependen del sigma que se
    recibe (sigma_block_size, sigma_ksize).
    Recibe (k), (threshol) y (env) utilizados en las funciones creadas para la
    obtención de los puntos.
    Se puede indicar la escala (scale) para devolver la escala a la que
    pertenecen los puntos generados junto a estos.
    """

    # Se obtiene block_size y ksize
    block_size = get_block_size(sigma_block_size)
    ksize = get_ksize(sigma_ksize)

    # Se calculan los autovalores y los autovectores de la imagen
    vals_and_vecs = cv2.cornerEigenValsAndVecs(img, block_size, ksize)

    # El resultado es una matriz de 3 dimensiones, donde se encuentra cada pixel
    # de la imagen y por cada pixel se tienen 6 valores:
    # - l1, l2: autovalores no ordenados de M
    # - x1, y1: autovalores correspondientes a l1
    # - x2, y2: autovalores correspondientes a l2
    # Se tiene que trabajar solo con los autovalores de M, por lo que se
    # separa en 6 canales el resultado y se guarda l1 y l2.
    vals_and_vecs = cv2.split(vals_and_vecs)
    l1 = vals_and_vecs[0]
    l2 = vals_and_vecs[1]

    # Se utiliza el criterio de Harris implementado para obtener la matriz
    # con el valor asociado a cada pixel
    harris = selection_criteria_harris(l1, l2, k)

    # Se suprimen los valores no máximos
    max_points = not_max_suppression_harris(harris, threshold, env)

    # Se guarda la longitud para recorrer todos los puntos y añadir la escala
    length = len(max_points)

    # Se añade la escala a cada punto
    for i in range(length):
        max_points[i][0].size = scale

    # Se devuelven los puntos
    return max_points

def get_best_harris(points, n = 1000):
    """
    Recibe una lista con los puntos Harris obtenidos formados por:
    (coordenadas (x, y), valor Harris, escala)
    Se puede especificar el número máximo de puntos Harris (n) que se quieren
    devolver, siendo elegidos los (n) de mayor valor.
    """

    # Se consigue una lista con los valores de los puntos
    values = [point[1] for point in points]

    # Se ordenan los puntos por el valor Harris de mayor a menor, guardando solo
    # los n (1000) puntos con valor más alto
    index_values_sorted = np.argsort(values)[::-1]
    index_values_sorted = index_values_sorted[0:n]

    # Se devuelven los puntos seleccionados
    return [points[i] for i in index_values_sorted]

def refine_harris_points(img, points, size, env = 5, zero_zone = -1, \
    criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 100,\
    0.001)):
    """
    Refina la posición de unos puntos (points) recibidos sobre una imagen (img)
    teniendo en cuenta el tamaño de la ventana (env) utilizada en la obtención
    de los puntos Harris.
    Zero_zone indica la mitad del tamaño de la zona muerta que se queda en el
    medio de la zona de búsqueda. Si se indica -1 no se le da tamaño.
    """

    # Se crea un array de numpy (requerido por la función cornerSubPix) con
    # las coordenadas de los puntos que se quieren refinar
    new_points = list()
    index = list()

    # Se recorren todos los puntos y solo se obtienen los que pertenezcan a la
    # al nivel de la pirámide correspondiente
    for i, point in enumerate(points):

        if point[0].size == size:
            new_points.append(point[0].pt)
            index.append(i)

    new_points = np.float32(new_points)

    # Se refinan los puntos (deben ser float)
    cv2.cornerSubPix(img, new_points, (env, env), \
                        (zero_zone, zero_zone), criteria)

    # Se modifican las coordenadas con los puntos refinados
    for i, point in enumerate(new_points):
        points[index[i]][0].pt = (point[0], point[1])

    return points

def get_orientations(img, points, sigma, own = True):
    """
    A partir de una imagen (img) y los puntos Harris calculados (points),
    obtiene la orientación que tendrá cada uno de los puntos.
    Se especifica un sigma (sigma) para aplicar antes de calcular las derivadas
    un alisamiento con ese sigma.
    Devuelve la orientación de cada punto.
    """

    # Se calculan las derivadas con respecto de x y de y de la imagen, aplicando
    # antes un alisamiento definido por sigma
    img_x, img_y = convolution_d(img, sigma = sigma, own = own)

    # Para aplicar el arcotangente es necesario modificar la estructura del
    # array de puntos de manera que existan dos filas (x, y) y tantas columnas
    # como puntos se tengan, por lo que se hace la transpuesta.
    new_points = np.array([(int(point[0].pt[0]), int(point[0].pt[1])) \
                                                        for point in points]).T

    # Se calcula la arcotangente de cada punto que devuelve la orientación que
    # tendra dicho punto. Esa operación se realiza para todos los puntos que se
    # tienen en el array
    orientations = np.arctan2(img_x, img_y)[new_points[0], new_points[1]]

    for point, orientation in zip(points, orientations):
        point[0].angle = orientation/np.pi*180

    # Se devuelven los puntos con sus orientaciones
    return points

def get_descriptors(img, keypoints):
    """
    Calcula los descriptores de los KeyPoints de la imagen y los devuelve junto
    a los KeyPoints en formato cv2.KeyPoint.
    """

    # Se convierte la imagen a enteros para que se puedan calcular los
    # descriptores
    img = np.uint8(img)

    # De la lista con los puntos que se tenía, se crea otra lista con el tipo
    # de datos cv2.KeyPoint para poder utilizar sift.compute()
    keypoints = [point[0] for point in keypoints]

    # Se calculan los descriptores
    return sift.compute(img, keypoints)

def get_keypoints_descriptors(img, own = False):
    """
    Dada una imagen, calcula sus keypoints y descriptores y los devuelve
    """

    if own:

        # Se obtienen los puntos Harris
        keypoints = get_harris(img, scale = 1)

        # Se guardan los mejores
        keypoints = get_best_harris(points)

        # Se actualizan los puntos con las orientaciones que tienen
        keypoints = get_orientations(img, points, sigma=5, own = False)

        # Se calculan los descriptores utilizando los keypoints generados
        # manualmente
        keypoints, descriptors = get_descriptors(img, keypoints)

    else:

        # Se extraen los keyPoints y los descriptores de las dos primeras
        # imágenes
        keypoints, descriptors = sift.detectAndCompute(img, None)

    # Se devuelven los keypoints y descriptores calculados
    return keypoints, descriptors

def get_matches_knn(img1, img2, k = 2, ratio = 0.75, n = 30, flag = 2, \
                        get_data = False, improve = False):
    """
    Dadas dos imágenes (img1, img2), calcula los keypoints y los descriptores
    para obtener los matches de ambas imágenes con la técnica
    "Lowe-Average-2NN", utilizando los k vecinos más cercanos (por defecto, 2).
    Para obtener los mejores matches, se calcula utilizando un ratio (ratio).
    Se puede indicar el número de matches que se quieren mostrar (n).
    El flag permite indicar si se quieren que se muestren los keypoints y los
    matches (0) o solo los matches (2).
    Si se indica get_data = True, se devolverán los keypoints, descriptores y
    los matches.
    Si se indica el flag "improve" como True, elegirá los mejores matches.
    Devuelve una imagen que se compone de las dos imágenes con sus matches.
    """

    # Se obtienen los keypoints y los descriptores de las dos imágenes
    keypoints1, descriptors1 = get_keypoints_descriptors(img1)
    keypoints2, descriptors2 = get_keypoints_descriptors(img2)

    # Se crea el objeto BFMatcher de OpenCV
    bf = cv2.BFMatcher()

    # Se consiguen los puntos con los que hace match indicando los vecinos más
    # cercanos con los que se hace la comprobación
    matches = bf.knnMatch(descriptors1, descriptors2, k)

    # Si se indica que se elijan los mejores matches
    if improve:

        # Solo se hace para hacer una mejora de los matches
        # Se guardan los puntos que cumplan con un radio en concreto
        good = []

        # Se recorren todos los matches
        for p1, p2 in matches:

            # Si la distancia del punto de la primera imagen es menor que la
            # distancia del segundo punto aplicándole un ratio, el punto se
            # considera bueno
            if p1.distance < ratio*p2.distance:
                good.append([p1])

        # Se ordenan los matches dependiendo de la distancia entre ambos puntos
        # guardando solo los n mejores
        matches = sorted(good, key = lambda x:x[0].distance)

    # Se guardan solo algunos puntos (n) aleatorios
    matches = random.sample(matches, n)

    # Se crea la imagen que se compone de ambas imágenes con los matches
    # generados. El sexto parámetro es la variable donde se quiere guardar la
    # salida
    matcher = cv2.drawMatchesKnn(img1, keypoints1, img2, keypoints2, matches, \
                                None, flags = flag)

    # El resultado normal será la imagen con los matches
    res = matcher

    # Si se ha indicado que se devuelvan los datos, el resultado serán todos Los
    # datos calculados
    if get_data:
        res = ((keypoints1, descriptors1), (keypoints2, descriptors2), matches)

    return res

def get_homography(img1, img2, improve = False):
    """
    Obtiene la homografía entre dos imágenes y la devuelve.
    Si se indica el flag "improve" como True, elegirá los mejores matches.
    """

    # Se calculan los keypoints, descriptores y matches de ambas imágenes
    data = get_matches_knn(img1, img2, get_data = True, improve = improve)

    # Se obtienen los keypoints, descriptores y matches entre ellos
    keypoints1 = data[0][0]
    descriptors1 = data[0][1]

    keypoints2 = data[1][0]
    descriptors2 = data[1][1]

    matches = data[2]

    # Se crean los puntos de la imagen fuente generando una matriz de filas
    # como matches hay y 2 columnas
    src_points = np.float32([keypoints1[point[0].queryIdx].pt \
                                    for point in matches]).reshape(-1, 1, 2)

    # Se crean los puntos de la imagen destino generando una matriz de filas
    # como matches hay y 2 columnas
    dst_points = np.float32([keypoints2[point[0].trainIdx].pt \
                                    for point in matches]).reshape(-1, 1, 2)

    # Para generar una homografía son necesarios al menos cuatro puntos, si no
    # se obtienen esos cuatro puntos, se devuelve una matriz 3x3 vacía
    if len(src_points) > 4:
        M, mask = cv2.findHomography(src_points, dst_points, cv2.RANSAC, 1)
    else:
        M = np.ndarray((3, 3), dtype=cv2.float32)

    # Se devuelve la homografía
    return M

def get_homography_to_center(img, width, height):
    """
    Calcula la homografía necesaria para llevar una imagen (img) al centro de
    un mosaico de tamaño (width x height)
    """

    # Se obtienen las dimensiones de la imagen central
    height_center, width_center, z_center = img.shape

    # Se calculan las coordenadas de inicio de la imagen central
    t_width = width/2 - width_center/2
    t_height = height/2 - height_center/2

    # Se crea una matriz 3x3 con los valores a 0 menos la diagonal a 1 y los
    # valores de la traslación
    m0 = np.array([[1, 0, t_width], [0, 1, t_height], [0, 0, 1]], \
                    dtype=np.float32)

    return m0

def get_mosaic(*imgs, crop = True, improve = False, blend = True, \
                                                            width_blend = 100):
    """
    Genera un mosaico a partir de las imágenes recibidas (imgs).
    Se puede especificar si se recortan los bordes de la imagen (crop = True).
    Si se indica el flag "improve" como True, elegirá los mejores matches.
    Con la variable "blend" se puede indicar si se quiere que los bordes de
    las imágenes se difuminen entre ellas.
    "width_blend" indica el ancho de la zona que se va a utilizar para la mezcla
    de las imágenes con el algoritmo de Burt-Adelson.
    Devuelve una imagen que contiene el mosaico.
    """

    # Se obtiene el siguiente tamaño potencia de 2 para poder hacer
    # correctamente las transformaciones con la pirámide laplaciana
    width_blend = next_power_two(width_blend)

    imgs_aux = copy.deepcopy(imgs)

    # Se calcula cuál es el número de la imagen que se encuentra en el centro
    index_img_center = int(len(imgs)/2)

    # Se obtiene la imagen que está en el centro
    img_center =  imgs[index_img_center]

    # Se calcula el ancho del mosaico como la suma de todos los anchos de las
    # imágenes
    width = sum([img_center.shape[1] for img in imgs])

    # Se calcula la altura del mosaico
    height = imgs[0].shape[0]*2

    # Para componer el mosaico hace falta una homografía principal que se
    # encargue de llevar las imágenes desde su posición hasta el centro.
    # Se debe calcular manualmente indicando el desplazamiento de la imagen
    # central hasta el centro del mosaico
    m0 = get_homography_to_center(img_center, width, height)

    # Una vez que se tiene la homografía principal, se coloca la imagen central
    # en el centro del mosaico
    #img = cv2.warpPerspective(img_center, m0, (width, height), \
                                #borderMode=cv2.BORDER_TRANSPARENT)

    # Se crea una lista con las homografías
    homographies = [None] * len(imgs)

    # Se añade la principal en la posición que ocupa la imagen central
    homographies[index_img_center] = m0

    # Se recorren todas las imágenes hacia la izquierda, haciendo las
    # homografías de unas con otras y colocándolas hacia el centro de la imagen
    # con la homografía principal
    for i in range(0, index_img_center)[::-1]:

        # Se calcula la homografía de la imagen actual con la de su derecha
        m = get_homography(imgs[i], imgs[i+1], improve = improve)

        # Se le aplica a la homografía la transformación que permite llevar
        # la imagen a la parte central
        m = np.dot(homographies[i+1], m)

        # Se guarda la homografía actual para que sea utilizada por el resto
        homographies[i] = m

        # Se aplica la transformación a la imagen para juntarla con las que ya
        # se encuentran en el mosaico
        #img = cv2.warpPerspective(imgs[i], m, (width, height), \
                                    #dst=img, borderMode=cv2.BORDER_TRANSPARENT)

    # Se recorren todas las imágenes hacia la derecha, haciendo las homografías
    # de unas con otras y colocándolas hacia el centro de la imagen con la
    # homografía principal
    for i in range(index_img_center+1, len(imgs)):

        # Se calcula la homografía de la imagen actual con la de su izquierda
        m = get_homography(imgs[i], imgs[i-1], improve = improve)

        # Se le aplica a la homografía la transformación que permite llevar
        # la imagen a la parte central
        m = np.dot(homographies[i-1], m)

        # Se guarda la homografía actual para que sea utilizada por el resto
        homographies[i] = m

        # Se aplica la transformación a la imagen para juntarla con las que ya
        # se encuentran en el mosaico
        #img = cv2.warpPerspective(imgs[i], m, (width, height), \
                                    #dst=img, borderMode=cv2.BORDER_TRANSPARENT)

    img = None

    img = cv2.warpPerspective(imgs[0], homographies[0], (width, height), \
                                dst=img, borderMode=cv2.BORDER_TRANSPARENT)

    for i in range(1, len(imgs)):

        # Si se quiere utilizar la mejora de transición entre imagen
        if blend:

            # Se calcula el punto hacia donde va la imagen de la derecha
            p1 = warp_point(0, 0, homographies[i])

            # Se obtiene la parte correspondiente de la imagen que ya se
            # encuentra en el canvas
            img_left = copy.deepcopy(img[:, \
                                int(p1[0]-width_blend):int(p1[0]+width_blend)])

        # Se lleva la siguiente imagen a su lugar correspondiente del canvas
        img = cv2.warpPerspective(imgs[i], homographies[i], (width, height), \
                                    dst=img, borderMode=cv2.BORDER_TRANSPARENT)

        # Si se quiere utilizar la mejora de transición entre imagen
        if blend:

            # Se recoge la otra parte de la imagen que se quiere mezcla, ésta
            # ya con el cambio fuerte de intensidad
            img_right = copy.deepcopy(img[:, \
                                int(p1[0]-width_blend):int(p1[0]+width_blend)])

            # Se mezclan las imágenes con el algoritmo de Burt-Adelson
            img_blend = blend_images(img_left, img_right, 0.5, 5, 1, \
                                                        return_blended = True)

            # Se sustituye el resultado en el canvas
            img[:, int(p1[0]-width_blend):int(p1[0]+width_blend)] = \
                                                    copy.deepcopy(img_blend)

    # Si se ha indicado, se recortan los bordes de la imagen
    if crop:
        img = crop_image(img)

    # Se devuelve la imagen
    return img

def crop_image(img):
    """
    Recorta todos los bordes negros a una imagen y la devuelve.
    Código completamente obtenido de:
    - https://stackoverflow.com/questions/13538748/crop-black-edges-with-opencv
    """

    # Se convierte la imagen a escala de grises
    gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)

    # Interesa guardar el segundo valor devuelto que contiene una matriz del
    # tamaño de la imagen con valores binarios, diferenciando el negro del
    # resto de colores
    _, thresh = cv2.threshold(gray, 1, 255, cv2.THRESH_BINARY)

    # Se llama a la función de OpenCV que devuelve las coordenadas de todos Los
    # puntos que forman el contorno
    contours = cv2.findContours(thresh, cv2.RETR_EXTERNAL, \
                                    cv2.CHAIN_APPROX_SIMPLE)[0]

    # Se obtienen las esquinas de los puntos obtenidos anteriormente
    x, y, w, h = cv2.boundingRect(contours)

    # Se devuelve la imagen delimitada por las esquinas
    return img[y:y+h,x:x+w]

def warp_point(x, y, m):
    """
    Función que calcula el lugar al que se mueve un punto (x, y) de una imagen
    al aplicar una homografía (m)
    """

    # Se calcula el desplazamiento con la tercera fila de la matriz que contiene
    # la homografía
    div = m[2,0]*x + m[2,1]*y + m[2,2];

    # Se calculan los nuevos puntos x, y con la transformación correspondiente
    # de la homografía
    x = (m[0,0]*x + m[0,1]*y + m[0,2])/div;
    y = (m[1,0]*x + m[1,1]*y + m[1,2])/div;

    return x, y

def blend_images(img1, img2, half = 0.5, num_levels = 5, \
                                        threshold = 1, return_blended = False):
    """
    Función que recibe dos imágenes (img1, img2) y, mediante el algoritmo de
    Burt-Adelson, las mezcla de manera que existe una suave transición entre
    ambas. Las imágenes deben tener el mismo tamaño para un correcto
    funcionamiento.
    Se puede especificar la zona de la imagen donde se quiere hacer la mezcla.
    Si half = 0.5, la mezcla se realiza a la mitad. Si half > 0.5, la mezcla se
    realiza gradualmente hacia la derecha. Si half < 0.5, la mezcla se realiza
    a la izquierda.
    num_levels indica el número de niveles que se quiere que tengan las
    pirámides gaussiana y laplaciana utilizadas para la realización de la
    mezcla.
    threshold indica el umbral de mezcla que se va utilizar. 1 es el umbral por
    defecto. Si se quiere que el área de mezcla sea más extensa, se puede
    aumentar dicho umbral.
    Se puede devolver la imagen mezclada (return_blended = True) o separada
    en dos mitades (return_blended = False).
    """

    # Se crea una máscara del tamaño de las imágenes donde se colocará la
    # mezcla de ambas
    mask = np.zeros_like(img1, dtype='float32')

    # Dependiendo de la zona donde se quiere hacer la mezcla, se rellena de 1's
    # una parte u otra de la máscara. Por defecto, la mezcla se realiza en la
    # mitad de la máscara (0.5)
    mask[:, int(img1.shape[1] * half):] = 1

    # Se guardan las dimensiones para saber la parte que se debe devolver
    img_shape = img1.shape

    # Se copian las imágenes y las máscaras para no modificar las originales
    img1 = img1.copy()
    img2 = img2.copy()
    mask = mask.copy()

    # Se genera la pirámide gaussiana para la imagen 1, 2 y la máscara que
    # contiene el resultado. La primera imagen que compone la pirámide gaussiana
    # es la original
    gaussian_pyramid_img1 = [img1]
    gaussian_pyramid_img2 = [img2]
    gaussian_pyramid_mask = [mask]

    # Se recorre tantas veces como niveles se quiera que tenga la pirámide y se
    # van creando las imágenes más pequeñas
    for i in range(num_levels):

        # Se reduce la imagen a la mitad
        img1 = cv2.pyrDown(img1)
        img2 = cv2.pyrDown(img2)
        mask = cv2.pyrDown(mask)

        # Se añade la imagen reducida a la pirámide
        gaussian_pyramid_img1.append(np.float32(img1))
        gaussian_pyramid_img2.append(np.float32(img2))
        gaussian_pyramid_mask.append(np.float32(mask))

    # Se genera la pirámide laplaciana para la imagen 1, 2 y máscara. La imagen
    # del último nivel de la pirámide gaussiana coincide con el último nivel de
    # la pirámide laplaciana
    laplacian_pyramid_img1 = [gaussian_pyramid_img1[num_levels-1]]
    laplacian_pyramid_img2 = [gaussian_pyramid_img2[num_levels-1]]
    laplacian_pyramid_mask = [gaussian_pyramid_mask[num_levels-1]]

    # Se recorren los niveles al contrario para componer con las altas
    # frecuencias la pirámide laplaciana
    for i in range(num_levels-1, 0, -1):

        # Para obtener las altas frecuencias se resta la imagen del nivel
        # superior al actual con el resultado de redimensionar el tamaño de la
        # imagen en el nivel actual
        laplacian_img1 = np.subtract(gaussian_pyramid_img1[i-1], \
                                            cv2.pyrUp(gaussian_pyramid_img1[i]))
        laplacian_img2 = np.subtract(gaussian_pyramid_img2[i-1], \
                                            cv2.pyrUp(gaussian_pyramid_img2[i]))

        # Se añaden las imágenes calculadas a la pirámide laplaciana
        laplacian_pyramid_img1.append(laplacian_img1)
        laplacian_pyramid_img2.append(laplacian_img2)

        # Se invierta la imagen correspondiente de la máscara
        laplacian_pyramid_mask.append(gaussian_pyramid_mask[i-1])

    # Se mezclan las imágenes correspondientes en cada nivel
    blended_imgs = []

    for l_img1, l_img2, g_mask in zip(laplacian_pyramid_img1, \
                                laplacian_pyramid_img2, laplacian_pyramid_mask):

        # Se mezclan ambas imágenes
        blended_imgs.append(l_img1 * g_mask + l_img2 * (threshold - g_mask))

    # Se reconstruye la imagen mezclada en base a las pirámides
    blended = blended_imgs[0]

    # Se recorren los niveles
    for i in range(1, num_levels):

        # Se redimensiona la imagen al doble
        blended = cv2.pyrUp(blended)

        # Se le suma las altas frecuencias para terminar de componer la imagen
        # sin perder calidad
        blended = cv2.add(blended, blended_imgs[i])

    # Si se debe devolver la imagen mezclada no se hace nada. Si se quiere
    # obtener las dos mitades, se separan
    if not return_blended:

        # Se crea una lista
        blended_imgs =  []

        # Se separan las imágenes en dos mitades
        blended_imgs.append(blended[:, 0:int(img_shape[1]/2)])
        blended_imgs.append(blended[:, int(img_shape[1]/2):])

        # Se sustituye la variable que se devuelve con las dos mitades
        blended = blended_imgs

    return blended

#
################################################################################

################################################################################
# Ejecuciones
#

# En caso de utilizar este fichero de manera externa para el uso de sus
# funciones no se ejecutará ningún código
if __name__ == "__main__":

    # Se insertan los ejercicios que se quieren ejecutar
    ex = list()

    # Ejecución del proyecto:
    ex.append(1)

    # Ejecución de la prueba del algoritmo Burt-Adelson
    #ex.append(2)

    # Listas para la visualización de las imágenes
    imgs = []
    imgs_title = []

    # Se cargan las imágenes
    img1 = mpimg.imread(path+'img_1.jpg')
    img2 = mpimg.imread(path+'img_2.jpg')
    img3 = mpimg.imread(path+'img_3.jpg')
    img4 = mpimg.imread(path+'img_4.jpg')
    img5 = mpimg.imread(path+'img_5.jpg')
    img6 = mpimg.imread(path+'img_6.jpg')
    img7 = mpimg.imread(path+'img_7.jpg')
    img8 = mpimg.imread(path+'img_8.jpg')
    img9 = mpimg.imread(path+'img_9.jpg')
    img10 = mpimg.imread(path+'img_10.jpg')
    img11 = mpimg.imread(path+'img_11.jpg')
    img12 = mpimg.imread(path+'img_12.jpg')
    img13 = mpimg.imread(path+'img_13.jpg')

    ############################################################################
    # Proyecto final
    #

    """
    Ejecución del proyecto
    """

    if 1 in ex:

        # Se añade el mosaico de las imágenes que se quieran
        imgs.append(get_mosaic(img1, img2, img3, img4, img5, img6, img7, img8, \
        img9, img10, img11, img12, img13,\
        improve = True, crop = True, blend = False, width_blend = 200))

        # Se añade un título a la imagen
        imgs_title.append("")

        # Se muestran las imágenes que se leen inicialmente
        show_images(imgs, imgs_title, cols = 1, gray = False)

        # Se limpian las listas para la próxima ejecución
        imgs.clear()
        imgs_title.clear()

    if 2 in ex:

        apple = cv2.imread(path+"apple.jpg")
        orange = cv2.imread(path+"orange.jpg")

        #imgs.append(apple)
        #imgs.append(orange)
        #imgs_title.append("")
        #imgs_title.append("")

        #show_images(imgs, imgs_title, cols = 1, gray = False)

        blended = blend_images(apple, orange, 0.5, 5, 1, return_blended = True)

        cv2.imwrite(path+"out_blended.png", blended)

    #
    ############################################################################

#
################################################################################

max_stars_repo_path: Proyecto final/proyecto_final.py

max_issues_repo_path: Proyecto final/proyecto_final.py

max_forks_repo_path: Proyecto final/proyecto_final.py


text:
#!/usr/bimn/env python

import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from scipy import interpolate
from scipy.interpolate import UnivariateSpline
import glob
from openpyxl import Workbook
from openpyxl import load_workbook
from mpl_toolkits.mplot3d import Axes3D
from sklearn import preprocessing, cross_validation, svm
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import StratifiedShuffleSplit
from sklearn.preprocessing import LabelBinarizer
from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import RandomForestRegressor
import pickle
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split

# # Star CCM Data Frame
# # df = pd.read_csv('adiabatic_sides_rear_temp.csv', index_col = 'Time')
# # df['Full rad'] = pd.read_csv('rear_temp_full_radiation.csv', index_col = 'Time')
# # df['Adiabatic'] = pd.read_csv('adiabatic_rear_temp.csv', index_col = 'Time')

# Analytical Solution Data Frame
# # df1 = pd.read_csv('analytical_updated.csv', index_col = 'Time')
# # df1['Parker'] = df1['Parker']
# # df1['Cowen'] = df1['Cowen']

# Post Processing Data Frame
# # df2 = pd.DataFrame()
# # df2['Adiabatic/parker'] = (df['Adiabatic']/df1['Parker'])
# # df2['Adiabatic_sides/cowen'] = (df['Adiabatic sides']/df1['Cowen'])

def listify(gen):
    # Convert a generator to a list
    return [k for k in gen]

def load_excel_data(filename):
	"""
	Passed a full excel file path, load a list of thermophysical properties.
	Functionality depends on proper ordering of values in the correct
	column and row - see example sheet.
	Args:
		-filename: A valid .xlsx or .xls file containing data
	"""
	wb = load_workbook(filename, data_only = True)
	ws = get_ws("Sheet1", wb)
	package = [float(k.value) for k in listify(ws.columns)[3][1:]]
	return package

def get_ws(ws, wb):
	if ws not in wb:
		print("Error: Could not find worksheet %s in workbook. Aborting")
		sys.exit()
	return wb.get_sheet_by_name(ws)

#half rise time function
def calc_halfrise(dataframe):
	#max temp

	# try:
		# dataframe = dataframe['Cowen']
	# except:
		# print('did not work')

	max_temp = dataframe.max()

	# starting temp and time
	start_temp = dataframe.iloc[0]

	#start_temp = np.argmax(dataframe > (initial_temp + 0.0001))
	#start_temp =  [i for i in dataframe if i > (initial_temp + 0.0001)][0]

	start_time = dataframe[dataframe == start_temp].index

	half_temp = (max_temp - start_temp) / 2 + start_temp

	x = dataframe.index
	y = dataframe

	#cubic spline fit
	s = UnivariateSpline(x, y, s=0, k=5)

	#creating x and ys
	xs = np.arange(0, x[-1] ,.000001)
	ys = s(xs)

	#finding half rise time by interpolating on the spline
	half_x = np.interp(half_temp, ys, xs)

	root = half_x #- start_time

	return root

def calc_actual_alpha():
	#for squares
	#pathMeta =r'C:\Users\grech\Documents\Python Scripts\senior\Square_meta'

	#for cylinder
	pathMeta =r'C:\Users\grech\Documents\Python Scripts\senior\Cylinder_meta'

	allMeta = glob.glob(pathMeta + "/*.xlsx")

	names = []
	alphas = []
	Qs = []
	aspects = []
	temps = []

	for file in allMeta:
		data = load_excel_data(file)

		#appends file name for square
		data.append(str(file)[-30:-10])

		name = data[-1]

		#defines the values in the data
		temp = data[0] #[K]
		rho = data[5] #[kg/m^3]
		cp = data[6] # [j/kg K]
		k = data[4] #[w/m K]
		l = data[2] #[m]
		r = data[3] #[m]
		Q = data[1] # [W/m^2]

		aspect = (2 * r) / l
		alpha = k / (rho * cp) #[m^2/s]

		names.append(name)
		alphas.append(alpha)
		Qs.append(Q)
		aspects.append(aspect)
		temps.append(temp)

	main = pd.DataFrame({'Actual Alpha' : alphas,
						'Name' : names,
						'Heat Deposited W/m^2' : Qs,
						'Aspect Ratio' : aspects,
						'Temeprature' : temps})
	#for squares
	#path =r'C:\Users\grech\Documents\Python Scripts\senior\Square'

	#for analytical squares
	#path = r'C:\Users\grech\Documents\Python Scripts\senior\Analytical_square'

	#for cylinders
	path =r'C:\Users\grech\Documents\Python Scripts\senior\Cylinder'

	#for analytical cylinders
	#path = r'C:\Users\grech\Documents\Python Scripts\senior\Analytical_cylinder'

	allFiles = glob.glob(path + "/*.csv")

	_names = []
	_alphas = []

	for file in allFiles:
		try:
			df_all = pd.read_csv(file, index_col = 'Physical Time: Physical Time (s)', nrows = 600)
		except:
			df_all = pd.read_csv(file, index_col = 'Time', nrows = 600)

		#use this name for square
		name = str(file)[-24:-4]

		#call half rise function
		half_time = calc_halfrise(df_all)

		l = .001

		alpha = (.1388 * l ** 2 )/ half_time

		_names.append(name)
		_alphas.append(alpha)

	_main = pd.DataFrame({'Calculated Alpha' : _alphas,
						 'Name' : _names})

	result = _main.merge(main, on = 'Name', how = 'left')

	result['calculated/actual'] = (result['Calculated Alpha'] / result['Actual Alpha']).astype(float)

	#result.to_csv('Square_results.csv')
	result.to_csv('Cylinder_results.csv')
	#result.to_csv('Analytical_square.csv')
	#result.to_csv('Analytical_cylinder.csv')
	return result

df = calc_actual_alpha()

# df = pd.concat([df1,df2])
# df.to_csv('Analytical_results.csv')

def plot_learning_curves(model, X, y):
    X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2, random_state=10)
    train_errors, val_errors = [], []
    for m in range(1, len(X_train)):
        model.fit(X_train[:m], y_train[:m])
        y_train_predict = model.predict(X_train[:m])
        y_val_predict = model.predict(X_val)
        train_errors.append(mean_squared_error(y_train_predict, y_train[:m]))
        val_errors.append(mean_squared_error(y_val_predict, y_val))

    plt.plot(np.sqrt(train_errors), "r-+", linewidth=2, label="Training")
    plt.plot(np.sqrt(val_errors), "b-", linewidth=3, label="Validating")
    plt.legend(loc="upper right", fontsize=14)   # not shown in the book
    plt.xlabel("Training set size (as a % of total data)", fontsize=14) # not shown
    plt.title('Learning Curve')
    plt.ylabel("RMSE", fontsize=14)              # not shown

df = pd.read_csv('star_results.csv')
data_cat = df['shape']
encoder = LabelBinarizer()
data_cat_1hot = encoder.fit_transform(data_cat)
df['shape'] = data_cat_1hot

split = StratifiedShuffleSplit(n_splits=1, test_size=0.25, random_state=42)
for train_index, test_index in split.split(df, df["Aspect Ratio"]):
    strat_train_set = df.loc[train_index]
    strat_test_set = df.loc[test_index]

X = np.array(strat_train_set.drop(['calculated/actual'], 1))
#X_train = preprocessing.scale(X)

y_train = strat_train_set['calculated/actual']

X_test = np.array(strat_test_set.drop(['calculated/actual'], 1))
#X_test = preprocessing.scale(X_test)

y_test = strat_test_set['calculated/actual']

clf = RandomForestRegressor()

clf.fit(X, y_train)
#clf.fit(X_t, t_y)

with open('r_forest_regressor.pickle','wb') as f:
    pickle.dump(clf, f)

# pickle_in = open('r_forest_regressor.pickle','rb')
# clf = pickle.load(pickle_in)

print(clf.score(X_test, y_test))

#This will do the mean squared error from 10 cross validations
# from sklearn.model_selection import cross_val_score
# forest_scores = cross_val_score(clf, X_train, y_train,
                                # scoring="neg_mean_squared_error", cv=10)
# forest_rmse_scores = np.sqrt(-forest_scores)
# print(forest_rmse_scores.mean())

#This will plot the learning curve error
#plot_learning_curves(clf, X_train, y_train)

# This code makes the 3D plots

# df1 = pd.read_csv('Analytical_cylinder.csv')
# df1 = df1[df1.Temeprature != 1573]
# df2 = pd.read_csv('Analytical_square.csv')

# df1 = pd.read_csv('Cylinder_results.csv')
# df1 = df1[df1.Temeprature != 1573]
# df2 = pd.read_csv('Square_results.csv')

# x1 = df1['Aspect Ratio'].values
# y1 = df1['Temeprature'].values
# z1 = df1['calculated/actual'].values

# x2 = df2['Aspect Ratio'].values
# y2 = df2['Temeprature'].values
# z2 = df2['calculated/actual'].values

# fig = plt.figure()
# ax = fig.gca(projection='3d')
# ax.scatter(x1, y1, z1, s=50, c = 'r', alpha = 1.0, marker = 'o')
# ax.scatter(x2, y2, z2, s=50, c = 'b', alpha = 0.5, marker = 's')
# ax.set_xlabel(r'Aspect Ratio$(\frac{diameter}{thckness})$')
# ax.set_ylabel('Temperature [K]')
# ax.set_zlabel(r' $(\frac{Calculated}{Actual}) \alpha$')
# ax.set_title(r'Deviation of $\alpha$ vs. Temeprature & Aspect Ratio')
# ax.tick_params(axis='y', width=10, labelsize=10, pad=0)
# plt.savefig('Analytical.jpg', dpi=1000)
# plt.show()

ext: py

reason: import numpy,from scipy


text:
# -*- coding: utf-8 -*-
"""
Created on Thu Feb 16 16:29:55 2017

@author: Administrator
"""
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from pandas import Series
from pandas import DataFrame

import os
import plotly.offline as pyof
import plotly.graph_objs as go
from plotly import figure_factory as ff

class Plotly_PyQt5():
    def __init__(self):
        plotly_dir = 'plotly_html'
        if not os.path.isdir(plotly_dir):
            os.mkdir(plotly_dir)

        self.path_dir_plotly_html = os.getcwd() + os.sep + plotly_dir

        # path_plotly_four_line = os.getcwd() + os.sep + plotly_dir + os.sep + 'mpl_two_line.html'
        # pyof.plot_mpl(a, resize=True, auto_open=False, filename=path_plotly_four_line)

    def get_plotly_path_if_hs300_bais(self, file_name='if_hs300_bais.html'):
        path_plotly = self.path_dir_plotly_html + os.sep + file_name
        df = pd.read_excel(r'if_index_bais.xlsx')

        '''绘制散点图'''
        line_main_price = go.Scatter(
            x=df.index,
            y=df['main_price'],
            name='main_price',
            connectgaps=True,  # 这个参数表示允许连接数据缺口
        )

        line_hs300_close = go.Scatter(
            x=df.index,
            y=df['hs300_close'],
            name='hs300_close',
            connectgaps=True,
        )
        data = [line_hs300_close, line_main_price]

        layout = dict(title='if_hs300_bais',
                      xaxis=dict(title='Date'),
                      yaxis=dict(title='Price'),
                      )

        fig = go.Figure(data=data, layout=layout)
        pyof.plot(fig, filename=path_plotly, auto_open=False)
        return path_plotly

    def get_plot_path_matplotlib_plotly(self, file_name='matplotlib_plotly.html'):
        path_plotly = self.path_dir_plotly_html + os.sep + file_name

        N = 50
        x = np.random.rand(N)
        y = np.random.rand(N)
        colors = np.random.rand(N)
        area = np.pi * (15 * np.random.rand(N)) ** 2  # 0 to 15 point radii
        scatter_mpl_fig = plt.figure()
        plt.scatter(x, y, s=area, c=colors, alpha=0.5)

        pyof.plot_mpl(scatter_mpl_fig, filename=path_plotly, resize=True, auto_open=False)
        return path_plotly

    def get_plotly_path_monte_markovitz(self, file_name='monte_markovitz.html', monte_count=400, risk_free=0.03):
        """
        """
        path_plotly = self.path_dir_plotly_html + os.sep + file_name
        df = pd.read_excel(r'data\组合.xlsx', index_col=[0])
        returns = df.pct_change()
        returns.dropna(inplace=True)
        noa = 3

        # 蒙特卡洛随机模拟结果
        port_returns = []
        port_variance = []

        for p in range(monte_count):
            weights = np.random.random(noa)
            weights /= np.sum(weights)
            port_returns.append(np.sum(returns.mean() * 50 * weights))  # 加入模拟的均值
            port_variance.append(np.sqrt(np.dot(weights.T, np.dot(returns.cov() * 50, weights))))  # 加入模拟的标准差

        port_returns = np.array(port_returns)
        port_variance = np.array(port_variance)
        color_array = (port_returns - risk_free) / port_variance  # 夏普比，不同的夏普比对应的颜色是不同的。

        trace1 = go.Scatter(
            x=port_variance,
            y=port_returns,
            mode='markers',
            marker=dict(
                size='6',
                color=color_array,  # 通过一个可变的变量表示颜色，结果是绘图颜色可变。
                colorscale='Viridis',
                # 设置colorbar
                colorbar=dict(
                    tickmode='linear',
                    tick0=color_array.min(),
                    dtick=(color_array.max() - color_array.min()) / 5,
                ),
                showscale=True,
            )
        )
        data = [trace1]

        pyof.plot(data, filename=path_plotly, auto_open=False)
        return path_plotly

    def get_plotly_path_combination_table(self, file_name='产品组合信息表.html', df=None, w=None):
        '''画产品组合信息表
        df为dataframe，末尾一列为组合。
        w为权重
        '''
        path_plotly = self.path_dir_plotly_html + os.sep + file_name
        risk_free = 0.03
        std_year = df.std() * 50
        mean_year = df.pct_change().mean() * 50 - risk_free
        sharp_year = mean_year / std_year
        max_drawback = (df / df.cummax()).min()

        _temp = {i[0]: i[1] for i in zip(['风险', '均值', '夏普比', '最大回撤'], [std_year, mean_year, sharp_year, max_drawback])}
        df_info = pd.concat(_temp, axis=1).T
        df_info.index.name = '指标'
        df_info = df_info.round(decimals=3)
        table = ff.create_table(df_info, index=True, index_title='指标')

        pyof.plot(table, filename=path_plotly, auto_open=False)
        return path_plotly

    def get_plotly_path_combination_pie(self, file_name='产品组合成分饼图.html', df=None, w=None):
        '''
        画产品组合成分饼图
        df为dataframe，末尾一列为组合。
        w为权重
        :return:
        '''
        path_plotly = self.path_dir_plotly_html + os.sep + file_name

        labels = df.columns[:-1]
        values = w
        trace = go.Pie(labels=labels, values=values, text='哈哈')
        layout = dict(title='产品组合成分图')
        fig = dict(data=[trace], layout=layout)

        pyof.plot(fig, filename=path_plotly, auto_open=False)
        return path_plotly

    def get_plotly_path_combination_versus(self, file_name='产品组合VS沪深300.html', df=None, df_base=None, w=None):
        '''
        画产品组合VS沪深300
        df为dataframe，末尾一列为组合。
        w为权重
        :return:
        '''
        path_plotly = self.path_dir_plotly_html + os.sep + file_name

        df_contra = pd.concat([df_base.close, df['组合']], axis=1, join='inner')
        df_contra = df_contra / df_contra.iloc[0, :]
        df_contra.rename(columns={'close': '沪深300'}, inplace=True)

        trace1 = go.Scatter(x=df_contra.index, y=df_contra.iloc[:, 0], mode='lines+markers',
                            name='沪深300', marker=dict(color='blue'))
        trace2 = go.Scatter(x=df_contra.index, y=df_contra.iloc[:, 1], mode='lines+markers',
                            name='产品组合', marker=dict(color='red'))
        data = [trace1, trace2]
        layout = {'title': '产品组合VS沪深300'}
        fig = dict(data=data, layout=layout)

        pyof.plot(fig, filename=path_plotly, auto_open=False)
        return path_plotly

max_stars_repo_name: allen-zqh/plotly

max_issues_repo_name: allen-zqh/plotly

max_forks_repo_name: allen-zqh/plotly

reason: import numpy


text:
import os
from numpy.core.multiarray import empty
import pandas as pd
import numpy as np
from keras.models import Sequential
from keras.layers import Dense
from keras.layers import LSTM
from keras.layers import Dropout
from keras.layers import *
from sklearn.preprocessing import MinMaxScaler

num_time_series = 100

csv_path = os.path.join(os.path.abspath(__file__), "../../../dir/nasdaq2007_17.csv")

model_path = os.path.join(os.path.abspath(__file__), "../../../models/Forecast_all_model.h5")

# Read the input file
df = pd.read_csv(csv_path, header=None, delimiter='\t')
file_ids = df.iloc[:, [0]].values
df = df.drop(df.columns[0], axis=1)
df = df.transpose()

# Training data: 80%, Test data: 20%
train_size = int(len(df) * 0.80)
test_size = len(df) - train_size

sc = MinMaxScaler(feature_range = (0, 1))

# Creating a data structure with 60 time-steps and 1 output
X_train = []
y_train = []
for step in range (0, num_time_series):

    # Scale and concatenate all time series
    training_set = df.iloc[:train_size, [step]].values
    training_set_scaled = sc.fit_transform(training_set)

    for i in range(60, train_size):
        X_train.append(training_set_scaled[i-60:i, 0])
        y_train.append(training_set_scaled[i, 0])

X_train, y_train = np.array(X_train), np.array(y_train)
X_train = np.reshape(X_train, (X_train.shape[0], X_train.shape[1], 1))
print(X_train.shape)

model = Sequential()
# Adding the first LSTM layer and some Dropout regularisation
model.add(LSTM(units = 50, return_sequences = True, input_shape = (X_train.shape[1], X_train.shape[2])))
model.add(Dropout(0.2))

# Adding a second LSTM layer and some Dropout regularisation
model.add(LSTM(units = 50, return_sequences = True))
model.add(Dropout(0.2))

# Adding a third LSTM layer and some Dropout regularisation
model.add(LSTM(units = 50, return_sequences = True))
model.add(Dropout(0.2))

# Adding a fourth LSTM layer and some Dropout regularisation
model.add(LSTM(units = 50))
model.add(Dropout(0.2))

# Adding the output layer
model.add(Dense(units = 1))

# Compiling the RNN
model.compile(optimizer = 'adam', loss = 'mean_squared_error')

# Fitting the RNN to the training set
model.fit(X_train, y_train, epochs = 60, batch_size = 2048, validation_split=0.1)

# Save model
if os.path.isfile(model_path) is False:
    model.save(model_path)

reason: import numpy,from numpy


text:
#!/usr/bin/env python
#Modified: April 12, 2021: VVK
# Extract images from rosbag files into .jpeg for training. Bag file names and
#   export paths are hardcoded.
import rospy
import std_msgs
from cv_bridge import CvBridge, CvBridgeError
import cv2
from sensor_msgs.msg import Image
import numpy as np

def main():
    rospy.init_node('image_extractor', anonymous=True)
    extractor=image_extractor()
    try:
        rospy.spin()
    except KeyboardInterrupt:
        print("Shutting down")

class image_extractor:

    def __init__(self):
        self.bridge=CvBridge()
        rospy.Subscriber('/Thermal_Panorama', Image, self.extract)
        self.fileCounter=1

    ## FOR EACH IMAGE:
    # def extract(self,data):
    #     imageFromBag=self.bridge.imgmsg_to_cv2(data, "mono8")
    #     # imageRGB = cv2.cvtColor(imageFromBag,cv2.COLOR_GRAY2RGB)
    #     # The export path is hardcoded:
    #     print(self.fileCounter)
    #     fileNumStr=str(self.fileCounter)
    #     fileNumStr=fileNumStr.zfill(5)
    #     cv2.imwrite('/home/vamsi/RainDataTrain/bag13/'+fileNumStr+'.jpeg',imageFromBag)
    #     self.fileCounter=self.fileCounter+1

    ## FOR ALTERNATE IMAGES:
    def extract(self,data):
        imageFromBag=self.bridge.imgmsg_to_cv2(data, "mono8")
        # imageRGB = cv2.cvtColor(imageFromBag,cv2.COLOR_GRAY2RGB)
        # The export path is hardcoded:
        if  (self.fileCounter % 2) == 0:
            print(self.fileCounter)
            fileNumStr=str(self.fileCounter)
            fileNumStr=fileNumStr.zfill(5)
            cv2.imwrite('/home/vamsi/RainDataTrain/Jul12Bag8/'+fileNumStr+'.jpeg',imageFromBag)
            self.fileCounter=self.fileCounter+1
        else:
            self.fileCounter=self.fileCounter+1

if __name__=='__main__':
    main()
#######################


text:
from . line import Line, SyntheticLine
from . circuit_element import CircuitElement
import numpy as np
from . utils import parse_dss_bus_name, parse_dss_phases, parse_phase_matrix
from typing import Tuple

class VoltageRegulator(CircuitElement):
    dss_module_name = 'RegControls'

    def __init__(self, name: str, dss):
        super().__init__(name, dss)
        self._set_lines()
        # set this regcontrol and its transformer as active,
        # in order to get the tap number
        dss.RegControls.Name(self.__name__)
        dss.Transformers.Name(self.transformer_name)
        self.tap = dss.RegControls.TapNumber()
        self.set_gamma(self.tap)
        self.Itx = np.zeros((3,), dtype=complex)
        # complex current entering/leaving the regControl node
        self.Ireg = np.zeros((3,), dtype=complex)

    def set_gamma(self, tapNumber: int) -> None:
        """
        Takes a Tap Number from opendss and maps it to a voltage ratio.
        sets self.gamma to the result.
        There are 33 default integer Tap Numbers, from -16 to +16,
        corresponding to the defaults of 0.90 to 1.10 voltage ratio
        """
        VRMIN = .90
        VRMAX = 1.10
        # move range [-16,16] to [0,1]
        result = (tapNumber + 16) / 32
        # compress to VRMAX - VRMIN
        result *= VRMAX - VRMIN
        # move to [VRMIN, VRMAX]
        result += VRMIN
        # assign to self
        self.gamma = result

    def _set_related_bus(self, dss):
        """
        Overrides superclass method to get both tx and regControl Bus.
        Set base vals, phases, and self.bus_name based on this bus
        """
        dss.RegControls.Name(self.__name__)
        self.transformer_name = dss.RegControls.Transformer()
        dss.Transformers.Name(self.transformer_name)
        tx, regControl_bus = dss.CktElement.BusNames()  # get upstream, regcontrol buses
        self.tx, self.rx = parse_dss_bus_name(
            tx),  parse_dss_bus_name(regControl_bus)
        self.regControl_bus = self.rx
        self.related_bus = self.regControl_bus  # alias for related bus
        self.key = (self.tx, self.rx)

    def _set_phases(self, dss):
        """
        override super() to set name from CktElement.BusNames
        """
        dss.RegControls.Name(self.__name__)
        self.phases = parse_dss_phases(dss.CktElement.BusNames()[0])
        self.phase_matrix = parse_phase_matrix(self.phases)

    def _set_lines(self):
        """
        Each voltage regulator is modeled with two lines:
        1. txBus --> regControlBus
        2. regControlBus--> txBus
        Finds or creates this pair of SyntheticLines
        from the line_group,  and assigns voltage regulators to both lines
        Add the upstream line to the topology for the main line_group
        and to the key_to_element_dict
        """
        self.line_tx_to_reg = SyntheticLine(unique_key=False,
                                            name=self.__name__ + '_to_reg',
                                            key=self.key)
        self.line_reg_to_tx = SyntheticLine(unique_key=False,
                                            name=self.__name__ + '_to_tx',
                                            key=self.key[-1::-1])
        # self.downstream_line = self._find_downstream_line(line_group)
        self.line_tx_to_reg.add_voltage_regulator(self)
        self.line_reg_to_tx.add_voltage_regulator(self)
        # self.downstream_line.add_voltage_regulator(self)

    # def _find_downstream_line(self, line_group):
    #     for line in line_group.get_elements():
    #         if line.tx == self.regControl_bus:
    #             return line
    #     return None

max_stars_repo_name: elaguerta/gigpower

max_issues_repo_name: elaguerta/gigpower

max_forks_repo_name: msankur/gigpower


text:
import os
import sys
import pandas as pd
from datetime import datetime
from numpy import isnan
import streamlit as st

p = os.path.abspath('../..')
if p not in sys.path:
    sys.path.insert(0,p)

# Functions
def str_time_to_sec(x):
    if x == 0:
        return 0
    pt = datetime.strptime(x,'%H:%M:%S')
    return pt.second + pt.minute*60 + pt.hour*3600

# Load class
@st.cache
def load_data():
    df = pd.read_excel(f'./bjj_videos.xlsx')
    df['tags'] = df['From'] + ', ' + df['To'].fillna('') + ', ' + df['Other'].fillna('')
    df['End_time'] = df['End_time'].fillna(0).apply(str_time_to_sec)
    df['Start_time'] = df['Start_time'].fillna(0).apply(str_time_to_sec)
    return df

st.set_page_config(layout="wide")
df = load_data()
st.title("Brazilian Jiu-Jitsu Instructional Videos")
c1, c2 = st.beta_columns((1, 2))
st.sidebar.header("Search Filters:")

radio = st.sidebar.radio('Language', ['English', 'Portuguese', 'All'])

# Side Pannel---------------------------
if radio == 'All':
    df_subset = df.copy()
else:
    df_subset = df[df['Language'] == radio]
type_options = list(df_subset['Type'].unique())
pos_type = st.sidebar.selectbox("Select Type", ['All'] + type_options)
if pos_type != 'All':
    df_subset = df_subset[df_subset['Type'] == pos_type]
teacher_options = list(df_subset['Teacher'].unique())
teacher = st.sidebar.selectbox("Select Teacher", ['All'] + teacher_options)
if teacher != 'All':
    df_subset = df_subset[df_subset['Teacher'] == teacher]

title = st.sidebar.text_input('Search using tags', 'Type any tag word')
if title != 'Type any tag word':
    text_matches = df_subset['tags'].str.contains(title).shape[0]
    df_subset = df_subset[df_subset['tags'].str.contains(title)]
else:
    text_matches = 0
df_subset = df_subset.reset_index(drop=True)

# Main page-------------------------------

# c1, c2, c3, c4 = st.beta_columns((2, 1, 1, 1))
with c1:
    st.subheader(f'{df_subset.shape[0]} videos found')
    with st.beta_expander('See all results'):
        st.table(df_subset[['Teacher', 'Type', 'From', 'To']])
    selected_indices = st.multiselect('Select only 1 result:', df_subset.index)
    selected_rows = df_subset.loc[selected_indices]
    st.write('### Selected video:', selected_rows[['Teacher', 'Type', 'From', 'To']].head(1))
with c2:
    if selected_rows.shape[0] > 0:
        url = selected_rows['Link'].iloc[0]
        start_time = selected_rows['Start_time'].iloc[0]
        st.video(url, start_time=start_time)
    else:
        st.write('No video selected')

reason: from numpy


text:
# -*- coding: utf-8 -*-
# @Author: harshit
# @Date:   2019-05-14 07:03:11
# @Last Modified by:   harshit
# @Last Modified time: 2019-05-14 08:02:59

from keras.layers import Dense, Input
from keras.models import Model
from keras.optimizers import Adam
import numpy as np
import matplotlib.pyplot as plt
import os

batch_size = 1024
epochs = 10
sessions = 5

def sin():
    inputs = Input((1,))
    l0 = Dense(
        512, activation='relu', kernel_initializer="random_uniform")(inputs)
    l1 = Dense(
        1024, activation='relu', kernel_initializer="random_uniform")(l0)
    l2 = Dense(
        1024, activation='relu', kernel_initializer="random_uniform")(l1)
    output = Dense(1, activation='linear')(l2)

    model = Model(inputs=[inputs], outputs=[output])
    adm = Adam(lr=0.001)
    model.compile(optimizer=adm, loss='mean_absolute_error', metrics=['mae'])
    model.summary()
    return model

def main():

    model = sin()
    # setting model callbacks

    for sess in range(sessions):
        X = np.random.random_sample((batch_size,)) * 10 * np.pi - (5 * np.pi)
        Y = np.sin(X)
        for epoch in range(epochs):
            print("Epoch", epoch)
            model.fit(X, Y, validation_split=0.1)
            y_ = model.predict(X)

            fig, axs = plt.subplots(1, 2)
            axs[0].scatter(X, Y, color="blue")
            axs[1].scatter(X, y_, color="red")
            name = str(sess) + "_" + str(epoch) + ".png"
            fig.savefig(os.path.join("sinx_graphs", name))
            plt.close()

if __name__ == '__main__':
    main()


text:
# -*- coding: utf-8 -*-
"""
Created on Sat Feb 27 20:16:39 2021

@author: jkcle
"""
from numpy.linalg import det
import numpy as np
from math import gamma

def coupled_normal_entropy(Sigma, kappa):
    """
    This function calculates the coupled entropy of a coupled Gaussian
    distribution using its Sigma matrix and kappa value.

    Parameters
    ----------
    Sigma : numpy ndarray
        The equivalent of a covariance matrix for a coupled Gaussian
        distribution.
    kappa : float
        A positive coupling value.

    Returns
    -------
    entropy : float
        The coupled entropy of the coupled Gaussian distribution with the
        covariance matrix equivalent of Sigma and coupling value kappa.

    """

    assert (((type(kappa) == float)
            | (type(kappa) == int))
            & (kappa > 0.)), "kappa must be positive!"

    assert ((type(Sigma) == np.ndarray)
            & (Sigma.shape[0] == Sigma.shape[1])), "Sigma is a square matrix!"

    # Find the number of dimensions using the square matrix Sigma.
    dim = Sigma.shape[0]

    # If the distribution is 1-D, the determinant is just the single value in
    # Sigma.
    if dim == 1:
        determinant = Sigma[0, 0]
    # Otherwise, calculate the determinant of the Sigma matrix.
    else:
        determinant = det(Sigma)

    # The coupled entropy calculation is broken up over several lines.
    entropy = (((np.pi/kappa)**dim) * determinant)**(kappa/(1+dim*kappa))
    entropy *= (1+dim*kappa)
    entropy *= (gamma(1/(2*kappa))/gamma(0.5*(dim + 1/kappa)))**(2*kappa
                                                                /(1+dim*kappa))
    entropy += -1
    entropy /= (2*kappa)

    # Return the coupled entropy.
    return entropy

scale = 1.

cauchy_entropy = np.log(4 * np.pi * scale)
coupled_entropy = coupled_normal_entropy(np.array([[scale]]), 1)

print(cauchy_entropy)
print(coupled_entropy)


text:
import numpy as np

from crf_dual_layer import DualCRFLayer,get_crf_training_loss
import lasagne,time,theano,logging,sys
import theano.tensor as T
from bionlp.utils.utils import theano_logsumexp
logging.basicConfig()
logger= logging.getLogger(__name__)

def setup_NN(worker,x_in,u_in,mask_in,y_in,params,numTags,emb_w):
    logger.info('Loading CRF-RNN network format with pairwise modeling. Forward Backward will be used for calculating the marginals while training. The output sequence will be calculated using posterior decoding.')
    batch_size = x_in.shape[0]
    print('X train batch size {0}'.format(x_in.shape))
    print('U train batch size {0}'.format(u_in.shape))
    print('Y train batch size {0}'.format(y_in.shape))
    print('mask train batch size {0}'.format(mask_in.shape))
    xt=theano.shared(x_in)
    mt=theano.shared(mask_in)
    ut=theano.shared(u_in)
    yt=theano.shared(y_in)
    premodel=time.time()

    l_in=lasagne.layers.InputLayer(shape=(None,params['maxlen'],1),input_var=xt.astype('int32'))
    l_u_in=lasagne.layers.InputLayer(shape=(None,params['maxlen'],u_in.shape[2]))
    l_mask=lasagne.layers.InputLayer(shape=(None,params['maxlen']))
    #---------- Replaced -------------

    l_u_reshape=lasagne.layers.ReshapeLayer(l_u_in,(-1,u_in.shape[2]))
    l_e_out=lasagne.layers.get_output(l_u_reshape)
    l_o_f = theano.function([l_u_in.input_var],l_e_out)
    logger.info("output shape for reshaped feature net = {0}".format(l_o_f(u_in.astype('float32')).shape))
    logger.info('Using {0} dimensional layer for extra features'.format(params['feature1']))
    l_u_dense=lasagne.layers.DenseLayer(l_u_reshape,params['feature1'])
    l_u_proc=lasagne.layers.ReshapeLayer(l_u_dense,(batch_size,params['maxlen'],1,params['feature1']))
    if params['word2vec']==1:
        l_emb=lasagne.layers.EmbeddingLayer(l_in,emb_w.shape[0],emb_w.shape[1],W=emb_w.astype('float32'))
        l_emb.add_param(l_emb.W, l_emb.W.get_value().shape, trainable=params['emb1'])
    else:
        l_emb=lasagne.layers.EmbeddingLayer(l_in,emb_w.shape[0],emb_w.shape[1])

    if params['emb2']>0:
        l_emb1=lasagne.layers.EmbeddingLayer(l_in,emb_w.shape[0],params['emb2'])
        l_emb=lasagne.layers.ConcatLayer([l_emb,l_emb1],axis=3)

    l_e_out=lasagne.layers.get_output(l_emb)
    l_p_out=lasagne.layers.get_output(l_u_proc)
    l_o_f = theano.function([l_in.input_var],l_e_out)
    l_p_f = theano.function([l_u_in.input_var],l_p_out)
    logger.info("output shape for emb net ={0}".format(l_o_f(x_in.astype('int32')).shape))
    logger.info("output shape for feature input={0}".format(l_p_f(u_in.astype('float32')).shape))

    l_emb=lasagne.layers.ConcatLayer([l_emb,l_u_proc],axis=3)
    l_e_out=lasagne.layers.get_output(l_emb)
    l_o_f = theano.function([l_in.input_var,l_u_in.input_var],l_e_out)
    logger.info("output shape for emb+feature net={0}".format(l_o_f(x_in.astype('int32'),u_in.astype('float32')).shape))

    #-----------Replaced end ---------

    dropout_backward=lasagne.layers.DropoutLayer(l_emb,params['noise1'])
    dropout_forward=lasagne.layers.DropoutLayer(l_emb,params['noise1'])

    rnn_nonlinearity = lasagne.nonlinearities.tanh
    #rnn_nonlinearity = lasagne.nonlinearities.elu
    logger.info('Using RNN nonlinearity {0}'.format(rnn_nonlinearity))

    backward1= lasagne.layers.LSTMLayer(dropout_backward,params['hidden1'],mask_input=l_mask,peepholes=False,forgetgate=lasagne.layers.Gate(b=lasagne.init.Constant(1.)),nonlinearity=rnn_nonlinearity,backwards=True,precompute_input=True)
    forward1= lasagne.layers.LSTMLayer(dropout_forward,params['hidden1'],mask_input=l_mask,peepholes=False,forgetgate=lasagne.layers.Gate(b=lasagne.init.Constant(1.)),nonlinearity=rnn_nonlinearity,precompute_input=True)

    crf_layer=lasagne.layers.ConcatLayer([forward1,backward1],axis=2)

    mid_out=lasagne.layers.get_output(crf_layer,deterministic=True)
    mid_output = theano.function([l_in.input_var,l_u_in.input_var,l_mask.input_var],mid_out)
    logger.info("output shape for for unary mid layer {0}".format(mid_output(x_in.astype('int32'),u_in.astype('float32'),mask_in.astype('float32')).shape))
    logger.info("output sum for for unary mid layer {0}".format(np.sum(mid_output(x_in.astype('int32'),u_in.astype('float32'),mask_in.astype('float32'))[0,:,:],axis=1)))

    dropout1=lasagne.layers.DropoutLayer(crf_layer,p=params['noise1'])
    dimshuffle_layer=lasagne.layers.DimshuffleLayer(dropout1,(0,2,1))
    pairwise = lasagne.layers.Conv1DLayer(dimshuffle_layer,numTags*numTags,2)
    pairwise = lasagne.layers.DimshuffleLayer(pairwise,(0,2,1))
    pairwise = lasagne.layers.batch_norm(pairwise)

    reshape_for_dense=lasagne.layers.ReshapeLayer(dropout1,(-1,params['hidden1']*2))
    output_nonlinearity=lasagne.nonlinearities.tanh
    dense1=lasagne.layers.DenseLayer(reshape_for_dense,numTags,nonlinearity=output_nonlinearity)
    dense1 = lasagne.layers.batch_norm(dense1)
    crf_layer=lasagne.layers.ReshapeLayer(dense1,(batch_size,params['maxlen'],numTags))

    mid_out=lasagne.layers.get_output(crf_layer,deterministic=True)
    mid_output = theano.function([l_in.input_var,l_u_in.input_var,l_mask.input_var],mid_out)
    logger.info("output shape for for unary mid layer {0}".format(mid_output(x_in.astype('int32'),u_in.astype('float32'),mask_in.astype('float32')).shape))
    logger.info("output sum for for unary mid layer {0}".format(np.sum(mid_output(x_in.astype('int32'),u_in.astype('float32'),mask_in.astype('float32'))[0,:,:],axis=1)))

    sum_layer=DualCRFLayer([crf_layer,pairwise,l_mask],mask_input=True)

    unaryPotential = lasagne.layers.get_output(crf_layer,deterministic = False)
    outp=lasagne.layers.get_output(sum_layer,deterministic = False)
    eval_out=lasagne.layers.get_output(sum_layer,deterministic=True)
    crf_output = theano.function([l_in.input_var,l_u_in.input_var,l_mask.input_var],eval_out)

    lstm_output = crf_output # Included for future functionality
    print("output shape for theano net",lstm_output(x_in.astype('int32'),u_in.astype('float32'),mask_in.astype('float32')).shape)
    t_out=T.tensor3()
    eval_cost=T.mean((eval_out-t_out)**2)

    all_params = lasagne.layers.get_all_params(sum_layer,trainable=True)
    print('Params :{0}'.format(all_params))
    logger.info('\'l1\' and \'l2\' Regularization is applied to all trainable parameters in the network')
    l2_cost= params['l2'] *lasagne.regularization.apply_penalty(all_params,lasagne.regularization.l2)
    l1_cost= params['l1'] *lasagne.regularization.apply_penalty(all_params,lasagne.regularization.l1)
    regularization_losses=l2_cost+l1_cost

    num_params = lasagne.layers.count_params(sum_layer,trainable = True)
    print('Number of parameters: {0}'.format(num_params))

    crf_train_loss=get_crf_training_loss(sum_layer,pairwise,t_out,numTags,params,x_in,u_in,y_in,mask_in,l_in,l_u_in,l_mask)
    cost=T.mean(crf_train_loss)+ regularization_losses
    cost_loss=T.mean(crf_train_loss)
    cost_regularization= regularization_losses

    updates = lasagne.updates.adagrad(cost, all_params,learning_rate=params['learning-rate'])
    #updates = lasagne.updates.adam(cost, all_params,learning_rate=params['learning-rate'])
    if params['momentum']==2:
        updates = lasagne.updates.apply_nesterov_momentum(updates,all_params,momentum=0.9)
    if params['momentum']==1:
        updates = lasagne.updates.apply_momentum(updates,all_params,momentum=0.9)

    train = theano.function([l_in.input_var,l_u_in.input_var,t_out,l_mask.input_var],cost,updates=updates)

    compute_cost = theano.function([l_in.input_var,l_u_in.input_var,t_out,l_mask.input_var],cost)
    compute_cost_loss = theano.function([l_in.input_var,l_u_in.input_var,t_out,l_mask.input_var],cost_loss)
    compute_cost_regularization = theano.function([],cost_regularization)
    acc_=T.sum(T.eq(T.argmax(eval_out,axis=2),T.argmax(t_out,axis=2))*l_mask.input_var)/T.sum(l_mask.input_var)
    compute_acc = theano.function([l_in.input_var,l_u_in.input_var,t_out,l_mask.input_var],acc_)
    print('Time to build and compile model {0}'.format(time.time()-premodel))

    return {'crf_output':crf_output,'lstm_output':lstm_output,'train':train,'compute_cost':compute_cost,'compute_acc':compute_acc,'compute_cost_loss':compute_cost_loss,'compute_cost_regularization':compute_cost_regularization,'final_layers':sum_layer}

max_stars_repo_name: abhyudaynj/LSTM-CRF-models

max_issues_repo_name: gabrielStanovsky/LSTM-CRF-models

max_forks_repo_name: gabrielStanovsky/LSTM-CRF-models

reason: import numpy,import theano


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

xm = np.array([78, 82, 72, 76, 74, 69])
xp = np.array([65, 70, 62, 82])

ym = np.array([100, 108, 117, 90, 87, 83])
yp = np.array([74, 79, 74, 111])

plt.scatter(xp, yp,c = 'c', marker = 'o', label = "Puiši")
plt.scatter(xm, ym,c = 'm', marker = 'o', label = "Meitenes")
plt.title("Stress/miera pulss")
plt.xlabel("Miera stāvokļa pulss")
plt.ylabel("Vidējais pulss stresa laikā")
x = np.array([78, 82, 76, 74, 69, 65, 70, 62, 82])
y = np.array([100, 108, 90, 87, 83, 74, 79, 74, 111])
m, b = np.polyfit(x, y, 1)
plt.plot(x, m*x + b)
plt.legend(loc = 2)

plt.legend(loc = "center left", bbox_to_anchor = (1, 0.5), numpoints = 1)

plt.show()

max_stars_repo_name: MartaRudzate/pulss

max_issues_repo_name: MartaRudzate/pulss

max_forks_repo_name: MartaRudzate/pulss


text:
from typing import Union

import numpy as np
from casadi import MX, SX, DM

class Mapping:
    """
    Mapping of index set to a different index set

    Example of use:
        - to_map = Mapping([0, 1, 1, 3, -1, 1], [3])
        - obj = np.array([0.1, 0.2, 0.3, 0.4, 0.5, 0.6])
        - mapped_obj = to_map.map(obj)
    Expected result:
        - mapped_obj == np.array([0.1, 0.2, 0.2, -0.4, 0, 0.1])

    Attributes
    ----------
    map_idx: list[int]
        The actual index list that links to the other set, an negative value links to a numerical 0

    Methods
    -------
    map(self, obj: list) -> list
        Apply the mapping to an obj
    len(self) -> int
        Get the len of the mapping
    """

    def __init__(self, map_idx: Union[list, tuple, range, np.ndarray]):
        """
        Parameters
        ----------
        map_idx: Union[list, tuple, range]
            The actual index list that links to the other set
        """
        self.map_idx = map_idx

    def map(self, obj: Union[tuple, list, np.ndarray, MX, SX, DM]) -> Union[np.ndarray, MX, SX, DM]:
        """
        Apply the mapping to an matrix object. The rows are mapped while the columns are preserved as is

        Parameters
        ----------
        obj: Union[np.ndarray, MX, SX, DM]
            The matrix to map

        Returns
        -------
        The list mapped
        """
        # Declare a zero filled object
        if isinstance(obj, (tuple, list)):
            obj = np.array(obj)

        if isinstance(obj, np.ndarray):
            if len(obj.shape) == 1:
                obj = obj[:, np.newaxis]
            mapped_obj = np.zeros((len(self.map_idx), obj.shape[1]))
        elif isinstance(obj, (MX, SX, DM)):
            mapped_obj = type(obj).zeros(len(self.map_idx), obj.shape[1])
        else:
            raise RuntimeError("map must be applied on np.ndarray, MX or SX")

        # Fill the positive values
        index_plus_in_origin = [abs(v) for v in self.map_idx if v is not None and v >= 0]
        index_plus_in_new = [i for i, v in enumerate(self.map_idx) if v is not None and v >= 0]
        mapped_obj[index_plus_in_new, :] = obj[index_plus_in_origin, :]  # Fill the non zeros values

        # Fill the negative values
        index_minus_in_origin = [abs(v) for v in self.map_idx if v is not None and v < 0]
        index_minus_in_new = [i for i, v in enumerate(self.map_idx) if v is not None and v < 0]
        mapped_obj[index_minus_in_new, :] = -obj[index_minus_in_origin, :]  # Fill the non zeros values

        return mapped_obj

    @property
    def len(self) -> int:
        """
        Get the len of the mapping

        Returns
        -------
        The len of the mapping
        """

        return len(self.map_idx)

class BiMapping:
    """
    Mapping of two index sets between each other

    Attributes
    ----------
    to_second: Mapping
        The mapping that links the first variable to the second
    to_first: Mapping
        The mapping that links the second variable to the first
    """

    def __init__(
        self, to_second: Union[Mapping, int, list, tuple, range], to_first: Union[Mapping, int, list, tuple, range]
    ):
        """
        Parameters
        ----------
        to_second: Union[Mapping, list[int], tuple[int], range]
            The mapping that links the first index set to the second
        to_first: Union[Mapping, list[int], tuple[int], range]
            The mapping that links the second index set to the first
        """
        if isinstance(to_second, (list, tuple, range)):
            to_second = Mapping(to_second)
        if isinstance(to_first, (list, tuple, range)):
            to_first = Mapping(to_first)

        if not isinstance(to_second, Mapping):
            raise RuntimeError("to_second must be a Mapping class")
        if not isinstance(to_first, Mapping):
            raise RuntimeError("to_first must be a Mapping class")

        self.to_second = to_second
        self.to_first = to_first

max_stars_repo_name: Sina-Mehdiz/Bioptim

max_issues_repo_name: Sina-Mehdiz/Bioptim

max_forks_repo_name: Sina-Mehdiz/Bioptim


text:
__author__ = 'max'

import re
import numpy as np

def is_uni_punctuation(word):
    match = re.match("^[^\w\s]+$]", word, flags=re.UNICODE)
    return match is not None

def is_punctuation(word, pos, punct_set=None):
    if punct_set is None:
        return is_uni_punctuation(word)
    else:
        return pos in punct_set or pos == 'PU' # for chinese

def eval(batch_size, words, postags, heads_pred, types_pred, heads, types,lengths,
         punct_set=None, symbolic_root=False, symbolic_end=False):
    ucorr = 0.
    lcorr = 0.
    total = 0.
    ucomplete_match = 0.
    lcomplete_match = 0.

    ucorr_nopunc = 0.
    lcorr_nopunc = 0.
    total_nopunc = 0.
    ucomplete_match_nopunc = 0.
    lcomplete_match_nopunc = 0.

    corr_root = 0.
    total_root = 0.
    start = 1 if symbolic_root else 0
    end = 1 if symbolic_end else 0
    for i in range(batch_size):
        ucm = 1.
        lcm = 1.
        ucm_nopunc = 1.
        lcm_nopunc = 1.
        for j in range(start, lengths[i] - end):
            word = words[i][j]

            pos = postags[i][j]

            total += 1
            if heads[i][j] == heads_pred[i][j]:
                ucorr += 1
                if types[i][j] == types_pred[i][j]:
                    lcorr += 1
                else:
                    lcm = 0
            else:
                ucm = 0
                lcm = 0

            if not is_punctuation(word, pos, punct_set):
                total_nopunc += 1
                if heads[i][j] == heads_pred[i][j]:
                    ucorr_nopunc += 1
                    if types[i][j] == types_pred[i][j]:
                        lcorr_nopunc += 1
                    else:
                        lcm_nopunc = 0
                else:
                    ucm_nopunc = 0
                    lcm_nopunc = 0

            if heads_pred[i][j] == 0:
                total_root += 1
                corr_root += 1 if heads_pred[i][j] == 0 else 0

        ucomplete_match += ucm
        lcomplete_match += lcm
        ucomplete_match_nopunc += ucm_nopunc
        lcomplete_match_nopunc += lcm_nopunc

    return (ucorr, lcorr, total, ucomplete_match, lcomplete_match), \
           (ucorr_nopunc, lcorr_nopunc, total_nopunc, ucomplete_match_nopunc, lcomplete_match_nopunc), \
           (corr_root, total_root), batch_size


text:
#!/usr/bin/env python3

import fnmatch
import csv
from datetime import datetime, timedelta
import sys
import os
import numpy as np

DATETIME_FORMAT = "%Y-%m-%d %H:%M:%S" # Strip milliseconds ".%f"

def GetMSE(input_path):
    dates = []
    actual = []
    predict = []
    with open(input_path) as input_file:
        reader = csv.reader(input_file)
        for line in reader:
            #dates.append(datetime.strptime(line[0], DATETIME_FORMAT))
            actual.append(float(line[1]))
            predict.append(float(line[2]))

    y = np.array(actual)
    y_hat = np.array(predict)

    data_min = 2 - np.min([np.min(y), np.min(y_hat)])
    se = (np.log(y + data_min) - np.log(y_hat + data_min)) ** 2
    print("MSE of %s: %s" % (input_path, np.mean(se)))

    return se

# ==============================================
# main
# ==============================================
if __name__ == '__main__':
    losses = np.array([])
    for root, dirnames, filenames in os.walk(sys.argv[1]):
        for filename in sorted(fnmatch.filter(filenames, '*.csv')):
            print(filename)
            file_path = os.path.join(root, filename)
            losses = np.append(losses, GetMSE(file_path))

    print(np.mean(losses))